• Gurobi Staff

Hi,

What exactly do you have difficulties with? Could you post a minimal working example of what you already tried?

I would recommend having a look at our Introduction to Modeling with Python webinar. Additionally, having a glance at our Python examples should help.

To model a constraint of the form $$\sum_{j\in I \\\{i\}} x_{jit} = v_i^t$$, you can use the following snippet

import gurobipy as gpm = gp.Model()I = [0,1,2,3]J = [0,1,2,3]T = [0,1,2,3]x = m.addVars(I,J,T,name="x")v = m.addVars(I,T,name="v")m.addConstrs((gp.quicksum(x[j,i,t] for j in J if i != j) == v[i,t] for i in I for t in T), name="constr")

Best regards,
Jaromił

Dear Jaromił,

In my constraint, there is a decision variable (x[ijt]) which I want to model its input and output flow. And I don’t know how to code it.

v_{i}^{t}=f_{i}^{t}+\sum_{j \in R \setminus \{i\}}x_{jit}+y_{i}^{t}-\sum_{j \in R \setminus \{i\}} x_{ijt}\quad \forall i \in R, t \in T
\tag{3}

, in which, (i, j) are indices of the locations and t is time period.

Here is some parts of my code:

input data________________

# number of regions
n = 5
regions = [*range(1,n)]

# planning periods

t = 3
days = [*range(t)]

# Create a dictionary to capture the coordinates of regions and their demand
coordinates = {
1:(1, 2),  2:(2.5, 1),  3:(5, 1), 4:(6.5, 3.5)}

# Compute pairwise distance matrix
# numpy linalg norm = euclidean n=2

def distance(city1, city2):
c1 = coordinates[city1]
c2 = coordinates[city2]
diff = (c1[0]-c2[0], c1[1]-c2[1])
return math.sqrt(diff[0]*diff[0]+diff[1]*diff[1])

dist = {(c1, c2): distance(c1, c2) for c1, c2 in combinations(regions, 2)}

model_________________________

# Transshipment amount from region i to region j in period t
x = m.addVars(dist, days, vtype=GRB.CONTINUOUS, name='transshipment')

# Inventory level at region i at the end of period t
v = m.addVars(regions, days, vtype=GRB.CONTINUOUS, name='inventory_end')

# Relief supply receive from Centeral Agency to region i in period t

# Inventory level at region in at the beginning of period t
f = m.addVars(regions, days, vtype=GRB.CONTINUOUS, name='inventory_beginning')

# 3. Transshipment and Inventory  (my difficulty)

m.addConstrs((v[region, day] == f[region, day] + x.sum('*', day) + y[region, day] - x.sum('*', day)
for region in regions for day in days), name='balance')

• Gurobi Staff

Your sums have additional conditions so I would prefer using the quicksum function here.

m.addConstrs((v[region, day] == f[region, day] + gp.quicksum(x[j,region,day] for j in regions if j != region)               + y[region, day] - gp.quicksum(x[region,j,day] for j in regions if j != region)               for region in regions for day in days), name='balance')

Please note that this code results in a KeyError

KeyError: (2, 1, 0)

That is because $$\texttt{x[2,1,0]}$$ is not present. However, according to your formulation, it should be there. Your $$\texttt{dist}$$ dictionary looks like

{(1, 2): 1.8027756377319946, (1, 3): 4.123105625617661, (1, 4): 5.70087712549569, (2, 3): 2.5, (2, 4): 4.716990566028302, (3, 4): 2.9154759474226504}

so it looks like there is indeed something missing. You don't have backward arcs like $$\texttt{(2,1),(3,1),...}$$ and you don't have any arcs with $$\texttt{i=j}$$ which makes the condition $$j \in R \\\{i\}$$ redundant in your formulation.

Yes, that’s right. Could you please guide me how I can fix it?

I didn’t define my arcs correctly, resulting in key errors in other constraints too.

• Gurobi Staff

You are probably looking for itertools.product.

dist = {(c1, c2): distance(c1, c2) for c1, c2 in product(regions,regions)}print(dist)

results in

{(1, 1): 0.0, (1, 2): 1.8027756377319946, (1, 3): 4.123105625617661, (1, 4): 5.70087712549569, (2, 1): 1.8027756377319946, (2, 2): 0.0, (2, 3): 2.5, (2, 4): 4.716990566028302, (3, 1): 4.123105625617661, (3, 2): 2.5, (3, 3): 0.0, (3, 4): 2.9154759474226504, (4, 1): 5.70087712549569, (4, 2): 4.716990566028302, (4, 3): 2.9154759474226504, (4, 4): 0.0}

How I can output x[ijt] values in Gurobi?

I also tried to output other decision variables like v[it], but the output is not correct. Here is my code.

rows = days.copy()
columns = regions.copy()
v_plan = pd.DataFrame(columns=columns, index=rows, data=0.0)

for region, day in v.keys():
if (abs(v[region, day].x) > 1e-6):
v_plan.loc[region, day] = np.round(v[region, day].x, 1)
v_plan

• Gurobi Staff

What is the error you get?

You can only access solution values after a successful optimization run, i.e., after calling model.optimize and finding at least 1 feasible solution.

I call model.optimize() and get obj value, but when using the code for output display, it isn't right.

The output for one of my variables is here. I don't know why there is column [0]

• Gurobi Staff

Your days go from 0 to 2

# planning periodst = 3days = [*range(t)]print(days)>[0, 1, 2]

Columns are regions and rows are planning periods. So why there is a column[0] and a row [3]?

• Gurobi Staff

That's because $$\texttt{v.keys()}$$ looks like this

print(v.keys()) ( 1 , 0 ) ( 1 , 1 ) ( 1 , 2 ) ( 2 , 0 ) ( 2 , 1 ) ( 2 , 2 ) ( 3 , 0 ) ( 3 , 1 ) ( 3 , 2 ) ( 4 , 0 ) ( 4 , 1 ) ( 4 , 2 )

so in the $$\texttt{for}$$-loop you add the 0 column and 3 row. You have to swap columns and rows

rows = days.copy()columns = regions.copy()v_plan = pd.DataFrame(columns=rows, index=columns, data=0.0)

Thanks, could you help me how I can output x[ijt] values?

• Gurobi Staff

You can access the solution values of the x variables via

for i in regions:    for j in regions:        for d in days:            print(x[i,j,d])

I am not aware of a good way to show 3D data in a table via pandas. You could try asking in a forum dedicated for pandas, e.g., stackoverflow, or someone in this forum with more knowledge might help.

Thank you so much. I genuinely appreciate all your help.