# Solving CVRP TW HF in Python

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')

#### 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.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]
# 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