Constraint indexes are not read by the model
AnsweredHi,
I have a model for periodic vehicle routing problem which is supposed to determine visiting schedules for each node i \in I. The code is currently working, but in the wrong way. I guess the model does not read the
for t in T2:
for i in I2:
mdl.addConstr(v[i,t+20] >= v[i,t])
part properly, because when I run it, i get results like that:
<gurobi.Var v[1,1] (value 1.0)>
<gurobi.Var v[2,39] (value 1.0)>
<gurobi.Var v[2,40] (value 1.0)>
<gurobi.Var v[3,1] (value 1.0)>
<gurobi.Var v[4,1] (value 1.0)>
<gurobi.Var v[5,1] (value 1.0)>
<gurobi.Var v[6,1] (value 1.0)>
<gurobi.Var v[7,1] (value 1.0)>
<gurobi.Var v[8,1] (value 1.0)>
<gurobi.Var v[9,39] (value 1.0)>
<gurobi.Var v[9,40] (value 1.0)>
<gurobi.Var v[10,1] (value 1.0)>
I try to forbid visiting the same node consecutively before 20 days are passed if the node belongs to I2 set. How can I fix my model to achieve that?
Thank you in advance for your help.
Best,
import random
import gurobipy as gp
from gurobipy import*
mdl = Model("490_deneme")
'''Sets and parameters'''
N = []
for i in range(11):
N.append(i)
I = []
for i in range(1,11):
I.append(i)
#I2 = [3,9,12,15,16,21,23,27,28]
I2 = [2,9]
#I4 = [11,14]
I4 = []
K = [1, 2, 3]
T = []
for i in range(1,41):
T.append(i)
T2 = []
for i in range(1,21):
T2.append(i)
T2_ = []
for i in range(21,41):
T2_.append(i)
T4 = []
for i in range(1,11):
T4.append(i)
C=5
#q = {1:0.72,2:0.58,3:0.93,4:0.68,5:0.52,6:0.38,7:0.55,8:0.28,9:1.25,10:0.2,11:1.77,12:0.92,13:0.43,14:1.73,15:1.2,16:1.1,17:0.45,18:0.53,19:0.43,20:0.5,21:0.92,22:0.57,23:0.77,24:0.35,25:0.71,26:0.63,27:0.77,28:0.87,29:0.38,30:0.61}
q = {1:0.72,2:0.58,3:0.93,4:0.68,5:0.52,6:0.38,7:0.55,8:0.28,9:1.25,10:0.2}
FREQ = {1:1,2:2,3:1,4:1,5:1,6:1,7:1,8:1,9:2,10:1}
DIST = {(0,0):0,(0,1):53,(0,2):112,(0,3):65,(0,4):64,(0,5):98,(0,6):51,(0,7):65,(0,8):59,(0,9):29,(0,10):21,
(1,0):53,(1,1):0,(1,2):108,(1,3):73,(1,4):105,(1,5):74,(1,6):16,(1,7):19,(1,8):80,(1,9):34,(1,10):44,
(2,0):112,(2,1):108,(2,2):0,(2,3):48,(2,4):96,(2,5):47,(2,6):93,(2,7):96,(2,8):60,(2,9):125,(2,10):93,
(3,0):65,(3,1):73,(3,2):48,(3,3):0,(3,4):59,(3,5):52,(3,6):58,(3,7):67,(3,8):18,(3,9):81,(3,10):47,
(4,0):64,(4,1):105,(4,2):96,(4,3):59,(4,4):0,(4,5):110,(4,6):96,(4,7):109,(4,8):41,(4,9):92,(4,10):62,
(5,0):98,(5,1):74,(5,2):47,(5,3):52,(5,4):110,(5,5):0,(5,6):60,(5,7):58,(5,8):70,(5,9):100,(5,10):78,
(6,0):51,(6,1):16,(6,2):93,(6,3):58,(6,4):96,(6,5):60,(6,6):0,(6,7):15,(6,8):67,(6,9):42,(6,10):37,
(7,0):65,(7,1):19,(7,2):96,(7,3):67,(7,4):109,(7,5):58,(7,6):15,(7,7):0,(7,8):78,(7,9):51,(7,10):51,
(8,0):59,(8,1):80,(8,2):60,(8,3):18,(8,4):41,(8,5):70,(8,6):67,(8,7):78,(8,8):0,(8,9):80,(8,10):45,
(9,0):29,(9,1):34,(9,2):125,(9,3):81,(9,4):92,(9,5):100,(9,6):42,(9,7):51,(9,8):80,(9,9):0,(9,10):36,
(10,0):21,(10,1):44,(10,2):93,(10,3):47,(10,4):62,(10,5):78,(10,6):37,(10,7):51,(10,8):45,(10,9):36,(10,10):0}
'''Decision Variables'''
pairs= [(n,m,k,t) for n in N for m in N for k in K for t in T]
x = mdl.addVars(pairs, vtype=GRB.BINARY, name ="x")
bins=[(i,k,t) for i in I for k in K for t in T]
y = mdl.addVars(bins, vtype=GRB.BINARY, name ="y")
visits=[(i,t) for i in I for t in T]
v = mdl.addVars(visits, vtype=GRB.BINARY, name ="v")
mtzs=[(i,k,t) for i in I for k in K for t in T]
u = mdl.addVars(mtzs, name ="u")
'''Constraints'''
for i in I:
for k in K:
for t in T:
mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == y[i,k,t])
for i in I:
for k in K:
for t in T:
if t !=40:
mdl.addConstr(quicksum(x[i,n,k,t+1] for n in N if n !=i) == y[i,k,t])
mdl.addConstrs(quicksum(y[i,k,t] for k in K) == v[i,t] for i in I for t in T)
for n in I:
for i in I:
for t in T:
if n !=i:
mdl.addConstr(quicksum(x[n,i,k,t] for k in K) <= (v[n,t]+v[i,t])/2)
for t in T:
for j in I:
mdl.addConstr(quicksum(x[n,i,k,t] for n in N for k in K) == v[i,t])
mdl.addConstrs(quicksum(y[i,k,t] for k in K for t in T) == FREQ[i] for i in I)
mdl.addConstrs(x[i,i,k,t] == 0 for i in I for k in K for t in T)
mdl.addConstrs(quicksum(x[0,i,k,40] for i in I) == 0 for k in K)
mdl.addConstrs(u[j,k,t]  u[i,k,t] >= q[j]/FREQ[j]C*(1x[i,j,k,t]) for i in I for j in I for k in K for t in T if i!=j)
mdl.addConstrs(q[i] <= u[i,k,t] for i in I for k in K for t in T)
mdl.addConstrs(u[i,k,t] <= C for i in I for k in K for t in T)
for t in T2:
for i in I2:
mdl.addConstr(v[i,t+20] >= v[i,t])
'''Objective Function'''
obj=quicksum(x[n,m,k,t]*DIST[(n,m)] for n in N for m in N for k in K for t in T)
mdl.setObjective(obj, GRB.MINIMIZE)
'''Solve'''
mdl.optimize()
print('Obj: %g' % mdl.objVal)
x_cor = []
y_cor = []
for t in T:
for i in I:
for n in N:
for k in K:
if x[i,n,k,t].X == 1:
print(i,n,k,t)
for v in mdl.getVars():
if v.x==1:
print(v)

Hi Benay,The real issue is that the constraint below does not guarantee that the consecutive visits to node \(i\) are at least 20 days apart. The constraint below just ensures that if node \(i\) is visited at time \(t\), it will be visited at time \(t+20\) as well. It does not guarantee anything for the intermediate times between time \(t\) and time \(t+20\).
mdl.addConstr(v[i, t+20] >= v[i, t])
Given your explanation, if node \(i\) is visited at time \(t\), i.e., \(v_{it} =1\), it should not be visited in the next 20 days, i.e., \(\sum_{j = t+1}^{t+20} v_{ij} = 0\).
You can use the Gurobi indicator constraint API to implement:
\[\mbox{if } v_{it} =1 \Rightarrow \sum_{j = t+1}^{t+20} v_{ij} = 0\]
Specifically, you can replace the above constraint with the constraint below in your code snippet:
start_time, end_time = T[0], T[1] + 1
mdl.addConstrs(
(v[i, t] == 1)
>> (
quicksum(
v[i, j] for j in range(max(t + 1, start_time), min(t + 20, end_time))
) == 0
)
for i in I2
for t in T
)I tested this with your code snippet and I got:
<gurobi.Var v[1,1] (value 1.0)>
<gurobi.Var v[2,1] (value 1.0)>
<gurobi.Var v[2,39] (value 1.0)>
<gurobi.Var v[3,1] (value 1.0)>
<gurobi.Var v[4,1] (value 1.0)>
<gurobi.Var v[5,1] (value 1.0)>
<gurobi.Var v[6,1] (value 1.0)>
<gurobi.Var v[7,1] (value 1.0)>
<gurobi.Var v[8,1] (value 1.0)>
<gurobi.Var v[9,1] (value 1.0)>
<gurobi.Var v[9,39] (value 1.0)>
<gurobi.Var v[10,1] (value 1.0)>Best regards,
Maliheh
0 
Thanks a lot, Maliheh, you really saved my day.
Now that I handled the visiting period constraint, I got errors from the basic VRP constraints of my model. The main logic of my model is like this I have 10 nodes and a depot, 3 vehicles, and 40 days in my planning horizon. I try to construct a periodic vehicle routing schedule in which the visiting times of each node (regarding their frequency of visit) is satisfied, and a vehicle completes its tour (i.e, depart from depot (node 0), go to its assigned nodes for that period, and return to the depot at the end of the same period.). But each vehicle does not suppose to depart from depot at a given period, neither they should do that in every period. Everything should supposed to happen just when it is needed (and the visiting frequencies should determine this condition), and in the remaining times vehicles can remain idle in the depot, as this is what we essentially want when we try to minimize distance. similarly, not all the vehicles are supposed to be utilized if I can handle the tours with the same vehicle. I plan to solve this problem in a very large network consisting of 175 nodes, and I know when I do that I'm gonna need 3 vehicles or more. But since I cannot run my code in a simple pilot network, I cannot move forward. Here is my current code:
# * coding: utf8 *
"""
Created on Thu Jan 6 21:58:28 2022
@author: User
"""
import random
import pandas as pd
import gurobipy as gp
from gurobipy import*
import matplotlib.pyplot as plt
mdl = Model("periodic vrp")
'''Sets and parameters'''
N = []
for i in range(11):
N.append(i)
I = []
for i in range(1,11):
I.append(i)
I2 = [3,9]
K = [1, 2, 3]
T = []
for i in range(1,41):
T.append(i)
T2 = []
for i in range(1,21):
T2.append(i)
C=5
random.seed(1)
def distance(p1,p2):
x_dist=(p1[0]p2[0])
y_dist=(p1[1]p2[1])
euclidian=(x_dist**2+y_dist**2)**0.5
return euclidian
D = [] # set of main distribution centers
D_loc = {}
for i in range(11):
D.append(i)
D_loc[i] = ((random.randint(0,100)),random.randint(0,100))
DIST = {}
for i in D:
for j in D:
DIST[(i,j)] = distance(D_loc[i],D_loc[j])
q = {1:0.72,2:0.58,3:0.93,4:0.68,5:0.52,6:0.38,7:0.55,8:0.28,9:1.25,10:0.2}
FREQ = {1:1,2:1,3:2,4:1,5:1,6:1,7:1,8:1,9:2,10:1}
'''Decision Variables'''
pairs= [(n,m,k,t) for n in N for m in N for k in K for t in T]
x = mdl.addVars(pairs, vtype=GRB.BINARY, name ="x")
bins=[(i,k,t) for i in I for k in K for t in T]
y = mdl.addVars(bins, vtype=GRB.BINARY, name ="y")
visits=[(i,t) for i in I for t in T]
v = mdl.addVars(visits, vtype=GRB.BINARY, name ="v")
mtzs=[(i,k,t) for i in I for k in K for t in T]
u = mdl.addVars(mtzs, name ="u")
'''Constraints'''
#flow balance constraints and connection of decision variables to each other
# for i in I:
# for k in K:
# for t in T:
# mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == quicksum(x[i,n,k,t] for n in N if n !=i))
for i in I:
for k in K:
for t in T:
mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == y[i,k,t])
for i in I:
for k in K:
for t in T:
if t !=40:
mdl.addConstr(quicksum(x[i,n,k,t+1] for n in N if n !=i) == y[i,k,t])
#connection of decision variables to each other
mdl.addConstrs(quicksum(y[i,k,t] for k in K) == v[i,t] for i in I for t in T)
for n in I:
for i in I:
for t in T:
if n !=i:
mdl.addConstr(quicksum(x[n,i,k,t] for k in K) <= (v[n,t]+v[i,t])/2)
for t in T:
for j in I:
mdl.addConstr(quicksum(x[n,i,k,t] for n in N for k in K) == v[i,t])
#
mdl.addConstrs(quicksum(y[i,k,t] for k in K for t in T) == FREQ[i] for i in I) #every node should be visited as much as its required frequency
mdl.addConstrs(quicksum(x[0,i,k,40] for i in I) == 0 for k in K)
mdl.addConstrs(quicksum(x[0,i,k,t] for i in I) <= 1 for t in T for k in K)
#mdl.addConstrs(quicksum(x[i,0,k,t] for i in I) == 1 for t in T for k in K)
mdl.addConstrs(u[j,k,t]  u[i,k,t] >= q[j]/FREQ[j]C*(1x[i,j,k,t]) for i in I for j in I for k in K for t in T if i!=j)
mdl.addConstrs(q[i] <= u[i,k,t] for i in I for k in K for t in T)
mdl.addConstrs(u[i,k,t] <= C for i in I for k in K for t in T)
start_time, end_time = T[0], T[1] + 1
mdl.addConstrs(
(v[i, t] == 1)
>> (
quicksum(
v[i, j] for j in range(max(t + 1, start_time), min(t + 20, end_time))
) == 0
)
for i in I2
for t in T
)
# mdl.addConstrs(v[i,t+20]+v[i,t+30] >= v[i,t]+v[i,t+10] for i in I4 for t in T4)
'''Objective Function'''
obj=quicksum(x[n,m,k,t]*DIST[(n,m)] for n in N for m in N for k in K for t in T)
mdl.setObjective(obj, GRB.MINIMIZE)
'''Solve'''
mdl.optimize()
print('Obj: %g' % mdl.objVal)
for v in mdl.getVars():
if v.x==1:
print(v)and when I enforce the currently commented constraint of returning to the depot at the end of each period, which is:
mdl.addConstrs(quicksum(x[i,0,k,t] for i in I) == 1 for t in T for k in K)
i get infeasibility I guess because gurobi says that
AttributeError: Unable to retrieve attribute 'objVal'
I really tried my best to work it out but could not find any way for a month. If you have any suggestions, I would be more than happy to hear that.
PS: I already tried using gurobi's cbLazy constraints but it did not change anything for my case, or I could not implement it successfully, I am not sure to be honest. I directly retrieve the code
def subtourelim(model, where):
if where == GRB.Callback.MIPSOL:
# make a list of edges selected in the solution
vals = model.cbGetSolution(model._vars)
selected = gp.tuplelist((i, j) for i, j in model._vars.keys()
if vals[i, j] > 0.5)
# find the shortest cycle in the selected edge list
tour = subtour(selected)
if len(tour) < 40:
# add subtour elimination constr. for every pair of cities in tour
model.cbLazy(gp.quicksum(model._vars[i, j]
for i, j in combinations(tour, 2))
<= len(tour)1)from the tsp example in here. It gives the same results with my subtour elimination constraints written with u_ijk variables.
Thank you in advance for your help and cooperation.
Best,
0 
Hi Benay,To debug your code and narrow down the infeasibility in your model, you can use the Gurobi Model.computeIIS() method to find an irreducible inconsistent subsystem (IIS) and write it into a file as below:
mdl.computeIIS()
mdl.write("model.ilp")To make things easier, it would be best to use the smallest problem size that results in infeasibility.Since your problem is an extension of the classical VRP problem, you might want to consider first implementing a classical VRP problem using one of the formulations available in the literature. You can then build on top of the base VRP, the different constraints that you have in your model.Best regards,Maliheh0
Please sign in to leave a comment.
Comments
3 comments