•  Gurobi Staff
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, T[-1] + 1mdl.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

•   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: utf-8 -*-"""Created on Thu Jan 6 21:58:28 2022@author: User"""import randomimport pandas as pdimport gurobipy as gpfrom gurobipy import*import matplotlib.pyplot as pltmdl = 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=5random.seed(1)def distance(p1,p2):x_dist=(p1-p2)y_dist=(p1-p2)euclidian=(x_dist**2+y_dist**2)**0.5return euclidianD = [] # set of main distribution centersD_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 othermdl.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 frequencymdl.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*(1-x[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, T[-1] + 1mdl.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 I2for 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 solutionvals = 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 listtour = subtour(selected)if len(tour) < 40:# add subtour elimination constr. for every pair of cities in tourmodel.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.

Best,

•  Gurobi Staff
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,
Maliheh