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!
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)
num_cities = 22
center_box = (100, 200)
cities_coord, _ = make_blobs(n_samples=num_cities,
centers=2,
cluster_std=20,
center_box=center_box,
random_state = 2)
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))
from scipy.spatial import distance
dist_matrix = distance.cdist(cities_coord, cities_coord, 'euclidean')
from pulp import *
prob = LpProblem('prob', LpMinimize)
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)
d = 'N03'
m = 3
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
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
u = LpVariable.dicts('u', N, 0, num_cities-1, LpInteger)
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
%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
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)
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