I am trying to implement this gurobi code for CVRP TW HF in Pulp but I had no luck, the problem is called Capacitated Vehicle Routing Problem with Time Windows and Heterogeneous Fleet
I implemented in gurobi, but I had no luck in Pul, any help?
Here is the gurobi code:
from gurobipy import *
import numpy as np
import matplotlib.pyplot as plt
############## Data #############
n= 10
clients= [i for i in range(n) if i != 0]
nodes= [0] + clients
arcs= [(i,j) for i in nodes for j in nodes if i!=j]
## Demand
np.random.seed(0)
q= {n: np.random.randint(10,15) for n in clients}
q[0]=0
## Time Windows
e= {0:0,1:10,2:10,3:10,4:20,5:20,6:20,7:40,8:40,9:40,10:40} ## start time
l= {0:200,1:100,2:100,3:100,4:150,5:150,6:150,7:180,8:180,9:180,10:180} ## Max time
## Service Time
s= {n:np.random.randint(3,5) for n in clients} ##service time
s[0]=0
## Vehicles
vehicles= [1,2,3,4]
Q= {1:50,2:50,3:25,4:25}
## Coordinates
X= np.random.rand(len(nodes))*100
Y= np.random.rand(len(nodes))*100
### Distance & Time to Travel
distance= {(i,j): np.hypot(X[i]-X[j],Y[i]-Y[j]) for i in nodes for j in nodes if i !=j}
time = {(i,j): np.hypot(X[i]-X[j],Y[i]-Y[j]) for i in nodes for j in nodes if i !=j}
############### Plot ###################
plt.figure(figsize=(10,5))
plt.scatter(X,Y,color="blue")
plt.scatter(X[0],Y[0], color="red", marker="D")
plt.annotate("DC|$t_{%d}$=(%d$,%d$)" %(0,e[0],l[0]),(X[0]-1,Y[0]-5.5))
for i in clients:
plt.annotate('$q_{%d}=%d$|$t_{%d}$=(%d$,%d$)' %(i,q[i],i,e[i],l[i]),(X[i]+1,Y[i]))
plt.title("Capacitated Vehicle Routing Problem with Time Windows and Heterogeneous Fleet")
####################### Model ###############################
arc_var = [(i,j,k) for i in nodes for j in nodes for k in vehicles if i!=j]
arc_time= [(i,k) for i in nodes for k in vehicles ]
model= Model('VRPTWHF')
x= model.addVars(arc_var, vtype=GRB.BINARY, name="x")
t= model.addVars(arc_time, vtype=GRB.CONTINUOUS, name="t")
#### Objective
model.setObjective(quicksum(distance[i,j]*x[i,j,k] for i,j,k in arc_var), GRB.MINIMIZE)
#### Constraints
model.addConstrs(quicksum(x[0,j,k] for j in clients)<=1 for k in vehicles)
model.addConstrs(quicksum(x[i,0,k] for i in clients)<=1 for k in vehicles)
model.addConstrs(quicksum(x[i,j,k] for j in nodes for k in vehicles if i != j)==1 for i in clients)
model.addConstrs(quicksum(x[i,j,k] for j in nodes if i!=j) - quicksum(x[j,i,k] for j in nodes if i!=j) == 0 for i in nodes for k in vehicles)
model.addConstrs(quicksum(q[i]*quicksum(x[i,j,k] for j in nodes if i!=j) for i in clients )<= Q[k] for k in vehicles)
model.addConstrs((x[i,j,k]==1)>> (t[i,k]+s[i]+time[i,j] == t[j,k]) for i in clients for j in clients for k in vehicles if i!=j)
model.addConstrs(t[i,k] >= e[i] for i,k in arc_time)
model.addConstrs(t[i,k] <= l[i] for i,k in arc_time)
# model.Params.timeLimit= 60
# model.Params.MIPGAP= 0.1
model.optimize()
########### Print the output ################
print("Objective Function:" ,str(round(model.ObjVal,2)))
for v in model.getVars():
if v.x >0.9:
print(str(v.VarName)+"="+str(v.x))
#### Routes and Trucks
routes=[]
trucks=[]
K= vehicles
N= nodes
for k in vehicles:
for i in nodes:
if i != 0 and x[0,i,k].x>0.9:
aux=[0,i]
while i != 0:
j=i
for h in nodes:
if j != h and x[j,h,k].x>0.9:
aux.append(h)
i=h
routes.append(aux)
trucks.append(k)
print(routes)
print(vehicles)
#### Time Calculation
time_acum= list()
for n in range(len(routes)):
for k in range(len(routes[n])-1):
if k == 0:
aux=[0]
else:
i=routes[n][k]
j=routes[n][k+1]
t=time[i,j]+s[i]+aux[-1]
aux.append
time_acum.append(aux)
##### Plot
import matplotlib.patches as mpatches
plt.figure(figsize=(10,5))
plt.scatter(X,Y,color="blue")
plt.scatter(X[0],Y[0], color="red", marker="D")
plt.annotate("DC", (X[0]-1,Y[0]-5.5))
for r in range(len(routes)):
for n in range(len(routes[r])-1):
i = routes[r][n]
j = routes[r][n+1]
plt.plot([X[i],X[j]], [Y[i],Y[j]], color="r", alpha=.4)
for r in range(len(time_acum)):
for n in range(len(time_acum[r])):
i = routes[r][n]
plt.annotate('$q_{%d}=%d$ |$t_{%d}=%d$' %(i,q[i],i,time_acum[r][n]),(X[i]+1,Y[i]))
patch= [mpatches.Patch(label="Vehicle "
+str(trucks[n])+"|cap="+str(Q[trucks[n]])) for n in range(len(trucks))]
plt.legend(handles=patch,loc="best")
plt.title('Capacitated Vehicle Routing Problem with Time Windows and Heterogeneous Fleet')
I tried to formulate the problem using Pulp, here is my formulation
import pandas as pd
import numpy as np
from pulp import *
import matplotlib.pyplot as plt
#------------------------ Data --------------------------------
df= pd.read_excel('Data Egypt.xlsx')
df.columns
locations= df['Location']
n= df.shape[0]-1
clients= [i for i in range(n) if i != 0]
nodes= [0] + clients
arcs= [(i,j) for i in nodes for j in nodes if i!=j]
## Demand
Demand= list(df['Demand'][1:].values)
q= dict(zip(clients,Demand))
## Time Windows
start= list(df['From'].values)
end= list(df['To'].values)
e= dict(zip(nodes,start))
l= dict(zip(nodes,end))
speed= 35
## Service Time
service= (df['Service'].values)
s= dict(zip(clients,service))
s[0]=0
## Vehicles
vehicles= list(df['Vehicles'].values[:5])
capacity= list(df['Capacity'].values[:5])
Q= dict(zip(vehicles,capacity))
## Coordinates
X= df['X'].values
Y= df['Y'].values
## Circuty Factor
cf= 1.66
### Distance & Time to Travel
distance= {(i,j): (np.sqrt((X[i]-X[j])**2 + (Y[i]-Y[j])**2)*100*cf) for i in nodes for j in nodes if i !=j}
time = {i: distance[i]/speed for i in distance}
############### Plot ###################
# plt.figure(figsize=(10,7))
# plt.scatter(X,Y,color="blue")
# plt.scatter(X[0],Y[0], color="red", marker="D")
# plt.annotate("DC|$t_{%d}$=(%d$,%d$)" %(0,e[0],l[0]),(X[0],Y[0]))
# texts= [plt.annotate(f'cust{i}|q={q[i]}|t={e[i]}:{l[i]}', (X[i], Y[i])) for i in clients]
# adjust_text(texts,arrowprops={'arrowstyle':'fancy','alpha':.5 ,'color':'crimson'},expand_points=(1.5,1.5))
# plt.title("Capacitated Vehicle Routing Problem with Time Windows and Heterogeneous Fleet")
#------------------------- Model -----------------------------------------
arc_var= [(i,j,k) for i in nodes for j in nodes for k in vehicles if i != j]
arc_time= [(i,k) for i in nodes for k in vehicles]
model= LpProblem('CVRPTWHF', LpMinimize)
x= LpVariable.dicts('x',arc_var, cat="Binary")
t= LpVariable.dicts('t', arc_time, cat="Continuous")
### Objective Function
model+= lpSum(distance[i,j]*x[i,j,k] for i,j,k in arc_var) # Minimize Distance
# model+= lpSum(x[i,j,k] for k in vehicles for i in [0] for j in clients if j != i) # Minimize Number of Trucks
### Constraints
for j in clients:
for k in vehicles:
model+= lpSum(x[0,j,k]) <=1
for i in clients:
for k in vehicles:
model+= lpSum(x[i,0,k]) <=1
for j in nodes:
for k in vehicles:
for i in clients:
if i != j:
model += lpSum(x[i,j,k]) == 1
for j in nodes:
for i in nodes:
for k in vehicles:
if i != j:
model+= lpSum(x[i,j,k]) - lpSum(x[j,i,k]) ==0
for j in nodes:
for i in clients:
for k in vehicles:
if i != j:
model+= lpSum(q[i]) * lpSum(x[i,j,k]) <= Q[k]
for i in clients:
for j in clients:
for k in vehicles:
if i != j:
model+= (x[i,j,k]==1) >= (t[i,k]+s[i]+time[i,j] == t[j,k]) #We Didin't Follow Exactly the Code
for i,k in arc_time:
model+= t[i,k] >= e[i]
for i,k in arc_time:
model+= t[i,k] <= l[i]
solver= model.solve()
print(solver)
Output= -1
I didn’t have any luck, kindly support