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

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

Leave a Comment