Transport Analytics Training Series - Last Revision: October 2022

The Capacitated Vehicle Routing Problem¶

image.png

It appears that we are on to something!

The mTSP model gave us some promising results, so we decided to go ahead with the multiple delivery campaigns that will run in parallel. But it turns out that the sleighs have limited capacity.

Can you help?

We knew that the mTSP model is not the ideal solution. Not only are the tours uneven, but we have no way to take into account the capacity of the vehicle, or the possibility that the nodes have different levels of demand.

If we would like to incorporate these aspects in our analysis, then we would need to follow a VRP-based approach.

In [1]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from sklearn.datasets import make_blobs

%matplotlib inline

plot_size   = 30
plot_width  = 38
plot_height = 20

params = {'legend.fontsize': 'large',
          'figure.figsize': (plot_width,plot_height),
          'axes.labelsize': plot_size,
          'axes.titlesize': plot_size,
          'xtick.labelsize': plot_size*0.75,
          'ytick.labelsize': plot_size*0.75,
          'axes.titlepad': 25}
plt.rcParams.update(params)
In [2]:
num_customer = 10
num_facility = 1
num_vehicles = 3
In [3]:
center_box = (100, 200)

cust_coord, _  = make_blobs(n_samples=num_customer, 
                            centers=2, 
                            cluster_std=20, 
                            center_box=center_box,  
                            random_state = 2)

facil_coord, _ = make_blobs(n_samples=num_facility,  
                            centers=1, 
                            cluster_std=1, 
                            center_box=center_box, 
                            random_state = 50)
In [4]:
plt.scatter(cust_coord[:, 0], 
            cust_coord[:, 1], 
            s=plot_size*2, 
            cmap='viridis');

plt.scatter(facil_coord[:, 0], 
            facil_coord[:, 1], 
            s=plot_size*5, 
            cmap='viridis');

We combine all nodes into a single list.

In [5]:
all_coord = list(facil_coord) + list(cust_coord)
all_coord
Out[5]:
[array([148.83923613, 121.34372996]),
 array([165.7474142 , 131.60904527]),
 array([154.62857111, 148.43678344]),
 array([137.40408992, 140.40355585]),
 array([154.58363786, 167.03226365]),
 array([140.0088288 , 143.71274428]),
 array([107.73077851,  85.75767587]),
 array([153.65711856,  77.68686145]),
 array([160.09765683, 123.75665829]),
 array([144.43027807,  80.23411428]),
 array([122.44044584,  84.41247088])]
In [6]:
from scipy.spatial import distance

dist_matrix = distance.cdist(all_coord, all_coord, 'euclidean')
In [7]:
payload = 4
cost = 100
budget = 100000
In [8]:
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('D'+str(i+1)) for i in range(num_customer)]
N = S + D

We are using the above parameter variable to set the number of facilities to be opened.

In [9]:
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(N)} for i,customer in enumerate(N)}
In [11]:
from pulp import *

prob = LpProblem('prob', LpMinimize)
In [12]:
x = LpVariable.dicts('x', (N,N), lowBound = 0, upBound = 1, cat = LpInteger)
f = LpVariable.dicts('f', (N,N), lowBound = 0, cat = LpInteger)
In [13]:
prob += lpSum([x[i][j]*d[i][j] for j in N for i in N]) + lpSum([cost * x[i][j] for i in S for j in D])
In [14]:
for i in D:
    prob += lpSum([x[j][i] for j in N]) - lpSum([x[i][j] for j in N]) == 0

for i in D:
    prob += lpSum([f[j][i] for j in N]) - lpSum([f[i][j] for j in N]) == 1
    
for i in N:
    for j in N:
        prob += f[i][j] <= payload * x[i][j]

for i in N:
    for j in S:
        prob += f[i][j] == 0

for i in N:
    prob += x[i][i] == 0

for i in S:
    prob += lpSum([x[i][j] for j in D]) >= 1

for i in D:
    prob += lpSum([x[i][j] for j in N]) == 1
    
for i in D:
    for j in D:
        prob += x[i][j] + x[j][i] <= 1 

prob += lpSum([cost * x[i][j] for i in S for j in D]) <= budget

Note that in the 2nd constraint we use n to define the number of facilities to be opened.

In [15]:
status = prob.solve()

print(f'STATUS\n{LpStatus[status]}\n')

STATUS
Optimal
In [16]:
route_num = 1
for i,facility in enumerate(S):
    for j,customer in enumerate(D):
        if x[facility][customer].varValue >= 1:
            dest = customer
            print(f'Route #{route_num}:')
            print(f'{facility}, {customer}: {x[facility][customer].varValue}, {f[facility][customer].varValue}')
            route_num += 1
            while dest != facility:
                orig = dest
                for k in N:
                    if x[orig][k].varValue >= 1:
                        dest = k
                        print(f'{orig}, {dest}: {x[orig][dest].varValue}, {f[orig][dest].varValue}')

Route #1:
S1, D2: 1.0, 4.0
D2, D4: 1.0, 3.0
D4, D5: 1.0, 2.0
D5, D3: 1.0, 1.0
D3, S1: 1.0, 0.0
Route #2:
S1, D6: 1.0, 4.0
D6, D10: 1.0, 3.0
D10, D9: 1.0, 2.0
D9, D7: 1.0, 1.0
D7, S1: 1.0, 0.0
Route #3:
S1, D8: 1.0, 2.0
D8, D1: 1.0, 1.0
D1, S1: 1.0, 0.0
In [17]:
vect_orig = []
vect_dest = []
arrw_width = []
arrw_color = []

u = v = np.zeros((len(vect_orig),len(vect_orig)))

for i,orig in enumerate(N):
    for j,dest in enumerate(N):
        if x[orig][dest].varValue >= 1:
            vect_orig.append(all_coord[j])
            vect_dest.append(all_coord[i])
            arrw_width.append(f[orig][dest].varValue)

vect_orig = np.array(vect_orig)
vect_dest = np.array(vect_dest)

In the above code we are iterating through our sets of facilities and customers.

Wherever we find a pair that is connected, we will append the origin and destination coordinates of a link between them to vect_orig and vect_dest, respectively.

In [18]:
plt.scatter(cust_coord[:, 0], cust_coord[:, 1], s=plot_size*10, color='red', zorder = 10000);
plt.scatter(facil_coord[:, 0], facil_coord[:, 1], s=plot_size*20, color='yellow', zorder = 10000);

for i in range(len(vect_orig)):
    plt.arrow(vect_orig[i][0],vect_orig[i][1],vect_dest[i][0]-vect_orig[i][0],vect_dest[i][1]-vect_orig[i][1],
              width=0.2,#arrw_width[i] / 2000,
              head_width=0,color='green',zorder = 1)

frame1 = plt.gca(); frame1.axes.xaxis.set_ticklabels([]); frame1.axes.yaxis.set_ticklabels([]);

Summary¶

from pulp import *
from scipy.spatial import distance

# Calculate distance matrix 
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')

# Initialise problem
prob = LpProblem('prob', LpMinimize)

# Define sets of customers and facilities
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('D'+str(i+1)) for i in range(num_customer)]
H = [('H'+str(i+1)) for i in range(num_vehicles)]
N = S + D

# Prepare distance matrix parameters
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(N)} for i,customer in enumerate(N)}
p = {vehicle:payload[i] for i,vehicle in enumerate(H)}
c = {vehicle:cost_h[i] for i,vehicle in enumerate(H)}

# Define decision variables
x = LpVariable.dicts('x', (N,N), lowBound = 0, upBound = 1, cat = LpInteger)
f = LpVariable.dicts('f', (N,N), lowBound = 0, cat = LpInteger)

# Objective function
prob += lpSum([lpSum([x[i][j]*d[i][j]]) for j in N for i in N]) + lpSum([cost * x[i][j] for i in S for j in D])

# Constraint 1: Vehicle Equilibrium
for i in D:
    prob += lpSum([x[j][i] for j in N]) - lpSum([x[i][j] for j in N]) == 0

# Constraint 2: Flow Equilibrium (assuming demand of one unit of cargo at each demand point)
for i in D:
    prob += lpSum([f[j][i] for j in N]) - lpSum([f[i][j] for j in N]) == 1

# Constraint 3: Cargo flow must not exceed vehicle payload
for i in N:
    for j in N:
        prob += f[i][j] <= payload * x[i][j]

# Constraint 4: Return to the supply depot with no cargo
for i in N:
    for j in S:
        prob += f[i][j] == 0

# Constraint 5: Ensure no vehicle waits at any node
for i in N:
    prob += x[i][i] == 0

# Constraint 6: Ensure that at least one vehicle leaves the depot
for i in S:
    prob += lpSum([x[i][j] for j in D]) >= 1

# Constraint 7: All nodes must be visited
for i in D:
    prob += lpSum([x[i][j] for j in N]) == 1

# Constraint 8: Sub-tour elimination
for i in D:
    for j in D:
        prob += x[i][j] + x[j][i] <= 1 

# Constraint 9: Budget constraint
prob += lpSum([cost * x[i][j] for i in S for j in D]) <= budget

# Solve the model
status = prob.solve()