Transport Analytics Training Series - Last Revision: October 2022

The Multiple Travelling Salesman Problem¶

image.png

Santa realised that he could never be able to deliver all the toys by himself, and he therefore decided that the only solution would be to clone himself!

But the question still remains... In what sequence would the Santas visit the chimneys?

If we assume that the capacity of the sleigh is not a problem, we can always try to apply the mTSP model.

Our mTSP implementation builds upon the classical TSP model that we presented in the first notebook. We had to make some minor modifications to the MTZ SECs to ensure that they can handle the multiple tours that begin from our depot.

If we had left things as before, the additional tours would have been treated as subtours!

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   = 15
plot_width  = 16
plot_height = 8

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_cities = 22
In [3]:
center_box = (100, 200) 

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

for i in range(num_cities):
    plt.annotate(i+1, (cities_coord[i, 0]+0.3, cities_coord[i, 1]+0.3))
In [5]:
from scipy.spatial import distance
dist_matrix = distance.cdist(cities_coord, cities_coord, 'euclidean')
In [6]:
from pulp import *

prob = LpProblem('prob', LpMinimize)
In [7]:
N = [(f'N{(i+1):02}') for i in range(num_cities)]
c = {n_f:{n_t: round(dist_matrix[i][j],2) for j,n_t in enumerate(N)} for i,n_f in enumerate(N)}
x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)
In [8]:
d = 'N03'

m = 3
In [9]:
prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) 
                                      for i in N])

for i in N:
    if i!=d:
        prob += lpSum([x[i][j] for j in N  if i!=j]) == 1

for j in N:
    if j!=d:
        prob += lpSum([x[i][j] for i in N  if i!=j]) == 1
In [10]:
prob += lpSum([x[d][j] for j in N if d!=j]) == m

prob += lpSum([x[i][d] for i in N if i!=d]) == m
In [11]:
u = LpVariable.dicts('u', N, 0, num_cities-1, LpInteger)
In [12]:
H = 8  # Lower  number to be visited by salesmen
L = 2  # Higher number to be visited by salesmen
 
for i in N:
    if i != d:
        prob += u[i] + (H-2)*x[d][i] - x[i][d] <= H-1

for i in N:
    if i != d:
        prob += u[i] + x[d][i] + (2-L)*x[i][d] >= 2

for i in N:
    if i != d:
        prob += x[d][i] + x[i][d] <=1

for i in N:
    for j in N:
        if i != j and (i != d and j!= d):
            prob += u[i] - u[j] +H*x[i][j] +(H-2)*x[j][i] <= H-1
In [13]:
%time status = prob.solve()

print(f'\nSTATUS\n{LpStatus[status]}\n')
CPU times: user 30.6 ms, sys: 16.3 ms, total: 46.9 ms
Wall time: 14min 49s

STATUS
Optimal

In [14]:
def plot_solution(_N, _x, _cities_coord):
    vect_orig = []
    vect_dest = []
    u = v = np.zeros((len(vect_orig),len(vect_orig)))

    for j,n_t in enumerate(_N):
        for i,n_f in enumerate(_N):
            if i!=j:
                if _x[n_f][n_t].varValue > 0:
                    vect_orig.append(_cities_coord[i])
                    vect_dest.append(_cities_coord[j])

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

    plt.scatter(_cities_coord[:, 0], 
                _cities_coord[:, 1], 
                s=plot_size*2, 
                cmap='viridis',
                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],
                  color='black')
    
    for i in range(num_cities):
        plt.annotate(i+1, (cities_coord[i, 0]+0.3, cities_coord[i, 1]+0.3))

plot_solution(N, x, cities_coord)
In [15]:
for j,n_t in enumerate(N):
    for i,n_f in enumerate(N):
        if i!=j and x[n_f][n_t].varValue > 0:
                print(f'{n_f} -> {n_t}')
N22 -> N01
N17 -> N02
N01 -> N03
N04 -> N03
N08 -> N03
N12 -> N04
N03 -> N05
N03 -> N06
N21 -> N07
N07 -> N08
N06 -> N09
N02 -> N10
N20 -> N11
N19 -> N12
N18 -> N13
N16 -> N14
N11 -> N15
N05 -> N16
N03 -> N17
N10 -> N18
N13 -> N19
N09 -> N20
N15 -> N21
N14 -> N22