Obtaining wrong solutions for VRPTW with SDLs
AnsweredDear Gurobi community,
currently I am trying to implement a VRP with (hard) time windows extended by the opportunity to deliver packages to Shared Delivery Locations (SDL) instead of delivering customers who are served by these SDLs.
My implementation does provide feasible solutions, but these solutions are obviously totally wrong. I tried a lot, but I couldn't find the reason for these wrong results.
This is the Python code of the implemented model:
import numpy as np
import gurobipy as gp
import random
import time
random.seed(10)
def VRPTW_SDL(N,N_D,I,F,V,e,Tmax,a,l,delta,gamma,t,c,C):
start_time = time.time()
mod = gp.Model()
x = mod.addVars(N_D,N_D,vtype = 'B', name = 'x') #binary variable indicating whether j is visited directly after node i (1) or not (0)
Y = mod.addVars(N,F,vtype = 'B', name = 'Y') #binary variable indicating whether i is delivered to SDL (Shared Delivery Location) f (1) or not (0)
z = mod.addVars(F,vtype ='B', name = 'z') #binary variable indicating whether SDL f is visited (1) or not
s = mod.addVars(N_D,vtype = 'C',name = 's') #non-negative variable indicating the visit time at node i
mod.setObjective(gp.quicksum(c[i,j]*x[i,j] for i in N_D for j in N_D if j!=i)
+ gp.quicksum(delta*Y[i,f] for i in I for f in F)
+ gp.quicksum(gamma*x[0,j] for j in N), gp.GRB.MINIMIZE)
mod.addConstrs((gp.quicksum(x[i,j] for i in N_D if i!=j)
+ gp.quicksum(Y[i,f] for f in V[j-1])
== 1 for j in I),name = 'NB_1')
mod.addConstrs((gp.quicksum(x[j,i] for j in N_D ) ==
gp.quicksum(x[i,j] for j in N_D ) for i in N_D), name = 'NB_2')
mod.addConstrs((z[f] >= (1/len(I))*gp.quicksum(Y[i,f] for i in I)
for f in F), name = 'NB_3')
mod.addConstrs((gp.quicksum(x[i,f] for i in N_D) == z[f] for f in F), name = 'NB_4')
mod.addConstrs((s[i] + t[i,j] + e[i] - (2*Tmax)*(1-x[i,j]) <= s[j] for i in N_D for j in N if i!=j),name = 'NB_5')
mod.addConstrs(((-Tmax*(gp.quicksum(Y[i,f] for f in F) + a[i])) <= s[i] for i in I), name= 'NB_6a')
mod.addConstrs(((Tmax*(gp.quicksum(Y[i,f] for f in F) + l[i])) >= s[i] for i in I), name= 'NB_6b')
mod.addConstrs((s[j] + e[j] + t[j,0] <= Tmax for j in N), name = 'NB_7')
mod.addConstrs((gp.quicksum(Y[i,f] for i in I) <= B[f-C-1] for f in F), name = 'NB_8')
mod.addConstrs((Y[i,f] == 0 for i in I for f in F if f not in V[i-1]), name = 'NB_9')
mod.optimize()
mod.write('VRPTW_SDL.mps')
#If an optimal is found, print values of decision variables and objective function
if mod.status == gp.GRB.OPTIMAL:
print(mod.printAttr('X'))
print(mod.ObjVal)
#if model is infeasible, print those constraints that are not satisfied
elif mod.status == gp.GRB.Status.INFEASIBLE:
print('Optimization was stopped with status %d' % mod.status)
mod.computeIIS()
for c in mod.getConstrs():
if c.IISConstr:
print('%s' % c.constrName)
end_time = time.time()
elapsed_time = end_time - start_time
print("Elapsed time: ", elapsed_time)
C = 6 #Number of customers/requests
Num_F = 2 #Number of SDLs
F = [i for i in range(C+1,C+1+Num_F)] #Set of SDL-Nodes
I = [i for i in range(1,C+1)] #Set of AHD-customer nodes
N = I + F #Set of AHD-customer nodes and SDL nodes
N_D = N.copy()# Set N with depot node
N_D.insert(0,0)
V = [] #Set of SDL-nodes for each request i
for v in I:
F_cus_no = int(random.uniform(1,3))
F_cus = random.sample(F,F_cus_no)
F_cus.sort()
V.append(F_cus)
delta = 1 #Compensation for delivery to SDL
gamma = 3 #Fixed Cost for each Vehicle used
B = [] #Capacity of SDL-Boxes
for b in F:
B.append(int(random.uniform(1,3)))
#Distance matrix cij
c = np.zeros((len(N_D),len(N_D)))
for i in N_D:
for j in N_D:
if j!=i:
c[i,j] = int(random.uniform(5,8))
else:
c[i,j] = float("inf")
#travel time matrix tij
t = c #equals distance from i to j
e = [0]
e_c = [3 for i in I] #Service time of each customer i
e_f = [5 for f in F] #Service time of each SDL
e = e + e_c + e_f #Service time of customers and SDLs
Tmax = 35 #Maximum travel time for vehicles used
a = [0 for i in I] #lower bound of time windows of customers i
a.insert(0,0) #earliest time where vehicles are allowed to leave the depot
l = [int(random.uniform(20,30)) for i in I] #upper bound of time windows of customers i
l.insert(0,Tmax) #latest time to get back to the depot
Rundreise = VRPTW_SDL(N,N_D,I,F,V,e,Tmax,a,l,delta,gamma,t,c,C)
The result:
Variable X
-------------------------
x[0,0] 1
x[0,5] 1
x[1,1] 1
x[2,2] 1
x[3,3] 1
x[4,4] 1
x[5,0] 1
x[5,5] 1
x[6,6] 1
Y[7,7] 1
Y[7,8] 1
Y[8,8] 1
s[5] 5
14.0
I would be very happy if you could help me to find the reason for these wrong results.
Kind regards
Alex
-
Official comment
-
Gurobi Staff
This post is more than two years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum, or try Gurobot, our chatbot interface offering instant, expert-level support. -
-
-
Hi Alex,
The first issue is quite obvious: You have created variables x(i,i) that do not exist in the problem. We recommend only creating those variables that have a representation in your optimization problem. So, first create all feasible tuples (i,j) and use this list as parameter in addVars().
To debug model issues, it is recommended to construct a very small problem instance, only with a few nodes in this case, and write out the model as LP file with mod.write("model.lp"). This file format can be easily read with a text viewer.
Best regards,
Mario0 -
-
Hi Mario,
thanks for your help, now it works.
Kind regards
Alex
0 -
Post is closed for comments.
Comments
3 comments