The CVRP model that we introduced in the previous notebook can take a lot of time to solve, especially in larger problem instances.
As was also the case with the TSP problem, in real-world applications we have no option but to use a heauristic approach. In this notebook we are going to implement the Sweep heuristic - a simple, geometric approach that can produce valid VRP solutions very quickly.
As it is a quite crude heuristic, these solutions are far from perfect, but can serve as a very good starting point.
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from sklearn.datasets import make_blobs
from scipy.spatial import distance
%matplotlib inline
plot_size = 20
plot_width = 10
plot_height = 10
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}
plt.rcParams.update(params)
num_customer = 15
rnd_state = 5
seed_index = 1
cust_demand_low = 100
cust_demand_high = 200
vehicle_capacity = 1000
np.random.seed(seed=rnd_state)
center_box = (125, 150)
cust_coord, _ = make_blobs(n_samples=num_customer, centers=2, cluster_std=20,
center_box=center_box, random_state = rnd_state)
cust_coord = np.round(cust_coord,0)
facil_coord = (150, 150)
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');
for i in range(num_customer): plt.annotate(i+1, (cust_coord[i, 0]+0.3, cust_coord[i, 1]+0.3))
import math
def get_angle(a, b, c): #a: seed, b: depot, c: vertex
angle = math.degrees(math.atan2(c[1]-b[1], c[0]-b[0]) - math.atan2(a[1]-b[1], a[0]-b[0]))
return round(angle + 360 if angle < 0 else angle, 1)
depot = facil_coord
seed = cust_coord[0]
vertex = cust_coord[2]
We combine all nodes into a single list.
cust_name = [('C'+str(i+1)) for i in range(num_customer)]
cust_demand = np.random.randint(low=cust_demand_low, high=cust_demand_high, size = num_customer)
cust_angles = [get_angle(seed, depot, cust_coord[i]) for i in range(num_customer)]
cust = list(zip(cust_name, cust_coord[:, 0], cust_coord[:, 1], cust_angles, cust_demand))
df = pd.DataFrame(cust, columns = ['cust', 'x', 'y', 'angle', 'demand'])
df
| cust | x | y | angle | demand | |
|---|---|---|---|---|---|
| 0 | C1 | 107.0 | 143.0 | 0.0 | 199 |
| 1 | C2 | 134.0 | 140.0 | 22.8 | 178 |
| 2 | C3 | 100.0 | 161.0 | 338.3 | 161 |
| 3 | C4 | 133.0 | 178.0 | 292.0 | 116 |
| 4 | C5 | 112.0 | 135.0 | 12.3 | 173 |
| 5 | C6 | 130.0 | 146.0 | 2.1 | 108 |
| 6 | C7 | 150.0 | 162.0 | 260.8 | 162 |
| 7 | C8 | 131.0 | 141.0 | 16.1 | 127 |
| 8 | C9 | 97.0 | 133.0 | 8.5 | 130 |
| 9 | C10 | 146.0 | 135.0 | 65.8 | 180 |
| 10 | C11 | 113.0 | 140.0 | 5.9 | 107 |
| 11 | C12 | 179.0 | 142.0 | 155.3 | 176 |
| 12 | C13 | 154.0 | 184.0 | 254.0 | 115 |
| 13 | C14 | 123.0 | 159.0 | 332.3 | 153 |
| 14 | C15 | 111.0 | 131.0 | 16.7 | 180 |
df = df.sort_values(by=["angle"])
df
| cust | x | y | angle | demand | |
|---|---|---|---|---|---|
| 0 | C1 | 107.0 | 143.0 | 0.0 | 199 |
| 5 | C6 | 130.0 | 146.0 | 2.1 | 108 |
| 10 | C11 | 113.0 | 140.0 | 5.9 | 107 |
| 8 | C9 | 97.0 | 133.0 | 8.5 | 130 |
| 4 | C5 | 112.0 | 135.0 | 12.3 | 173 |
| 7 | C8 | 131.0 | 141.0 | 16.1 | 127 |
| 14 | C15 | 111.0 | 131.0 | 16.7 | 180 |
| 1 | C2 | 134.0 | 140.0 | 22.8 | 178 |
| 9 | C10 | 146.0 | 135.0 | 65.8 | 180 |
| 11 | C12 | 179.0 | 142.0 | 155.3 | 176 |
| 12 | C13 | 154.0 | 184.0 | 254.0 | 115 |
| 6 | C7 | 150.0 | 162.0 | 260.8 | 162 |
| 3 | C4 | 133.0 | 178.0 | 292.0 | 116 |
| 13 | C14 | 123.0 | 159.0 | 332.3 | 153 |
| 2 | C3 | 100.0 | 161.0 | 338.3 | 161 |
not_visited = df.index.tolist() # initialize list of unvisited nodes
cur_tour = 1
cur_visit = 1
cap_left = vehicle_capacity # keeping track of remaining vehicle capacity
while not_visited: # while unvisited customer exist
cust = not_visited[0] # go to the next unvisited customer
cust_demand = df.loc[cust, "demand"] # obtain their demand
if cust_demand > cap_left: # if we do not have enough capacity left
cur_tour += 1 # start a new tour (increment counter)
cur_visit = 1 # reset our vist counter to 1
cap_left = vehicle_capacity # and reset our vehicle capacity
df.loc[cust, "tour"] = cur_tour # assign customer to the current vehicle
df.loc[cust, "visit"]= cur_visit # keep track of the visit number
cur_visit += 1 # increment our visit counter
cap_left -= cust_demand # reduce remaining capacity
df.loc[cust, "cap left"]= cap_left # keep track of remaining capacity
not_visited.remove(cust) # remove customer from unvisited nodes
total_tours = cur_tour
df
| cust | x | y | angle | demand | tour | visit | cap left | |
|---|---|---|---|---|---|---|---|---|
| 0 | C1 | 107.0 | 143.0 | 0.0 | 199 | 1.0 | 1.0 | 801.0 |
| 5 | C6 | 130.0 | 146.0 | 2.1 | 108 | 1.0 | 2.0 | 693.0 |
| 10 | C11 | 113.0 | 140.0 | 5.9 | 107 | 1.0 | 3.0 | 586.0 |
| 8 | C9 | 97.0 | 133.0 | 8.5 | 130 | 1.0 | 4.0 | 456.0 |
| 4 | C5 | 112.0 | 135.0 | 12.3 | 173 | 1.0 | 5.0 | 283.0 |
| 7 | C8 | 131.0 | 141.0 | 16.1 | 127 | 1.0 | 6.0 | 156.0 |
| 14 | C15 | 111.0 | 131.0 | 16.7 | 180 | 2.0 | 1.0 | 820.0 |
| 1 | C2 | 134.0 | 140.0 | 22.8 | 178 | 2.0 | 2.0 | 642.0 |
| 9 | C10 | 146.0 | 135.0 | 65.8 | 180 | 2.0 | 3.0 | 462.0 |
| 11 | C12 | 179.0 | 142.0 | 155.3 | 176 | 2.0 | 4.0 | 286.0 |
| 12 | C13 | 154.0 | 184.0 | 254.0 | 115 | 2.0 | 5.0 | 171.0 |
| 6 | C7 | 150.0 | 162.0 | 260.8 | 162 | 2.0 | 6.0 | 9.0 |
| 3 | C4 | 133.0 | 178.0 | 292.0 | 116 | 3.0 | 1.0 | 884.0 |
| 13 | C14 | 123.0 | 159.0 | 332.3 | 153 | 3.0 | 2.0 | 731.0 |
| 2 | C3 | 100.0 | 161.0 | 338.3 | 161 | 3.0 | 3.0 | 570.0 |
vect_orig = []
vect_dest = []
arrw_width = []
dist_tot = 0
for t in range(total_tours):
d_tour = df[df['tour']==t+1]
arrw_width.append(vehicle_capacity)
v_from = depot
for index, row in d_tour.iterrows():
v_to = (row['x'],row['y'])
arrw_width.append(row['cap left'])
vect_orig.append(v_from)
vect_dest.append(v_to)
dist_tot += distance.euclidean(v_from, v_to)
v_from = v_to
vect_orig.append(v_from)
vect_dest.append(depot)
dist_tot += distance.euclidean(v_from, depot)
vect_orig = np.array(vect_orig)
vect_dest = np.array(vect_dest)
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([]);
dist_tot
484.89689520868654
