Transport Analytics Training Series - Last Revision: October 2022

The Sweep Method¶

   

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.

In [4]:
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)
In [5]:
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)
In [6]:
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)
In [7]:
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))
In [8]:
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)
In [9]:
depot = facil_coord
seed = cust_coord[0]
vertex = cust_coord[2]

We combine all nodes into a single list.

In [10]:
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
Out[10]:
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
In [11]:
df = df.sort_values(by=["angle"])
df
Out[11]:
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
In [12]:
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
Out[12]:
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
In [13]:
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)
In [14]:
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
Out[14]:
484.89689520868654