The Travelling Salesman Problem - MTZ Model

Santa Claus needs to visit 15 chimneys as quickly as possible to deliver presents to the children who have been good in 2020 and stuck to social distancing rules.

Can you help him find the shortest route crossing each home exactly once and return to the first home to then go back to the North Pole from there?

Santa has asked you to obtain the best possible solution to this problem.

Since Santa has asked for the best solution, we need to use an exact solution technique. We can therefore proceed to implement a mathematical model that implements the Miller, Tucker and Zemlin subtour elimination constraints (SECs).

We proceed to implement the model using PuLP. For the purposes of this notebook, we obtain the locations of the chimneys using our trusty blob maker.

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 = 15

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');

from pulp import *

prob = LpProblem('prob', LpMinimize)

We are back again in using PuLP.

from scipy.spatial import distance

dist_matrix = distance.cdist(cities_coord, cities_coord, 'euclidean')

dist_matrix.shape

(15, 15)

N = [('N'+str(i+1)) 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)

prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) 
                                      for i in N])

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

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

%time status = prob.solve()

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

CPU times: user 3.88 ms, sys: 8.05 ms, total: 11.9 ms
Wall time: 23.3 ms

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')

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}')

N13 -> N1
N12 -> N2
N7 -> N3
N1 -> N4
N8 -> N5
N5 -> N6
N3 -> N7
N6 -> N8
N10 -> N9
N9 -> N10
N2 -> N11
N11 -> N12
N4 -> N13
N15 -> N14
N14 -> N15

from pulp import *

prob = LpProblem('prob', LpMinimize)

N = [('N'+str(i+1)) for i in range(num_cities)]

x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)

c = {n_f:{n_t: dist_matrix[i][j] for j,n_t in enumerate(N)} 
                                 for i,n_f in enumerate(N)}

prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) 
                                      for i in N])

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

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

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

%time status = prob.solve()

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

CPU times: user 6.75 ms, sys: 9.54 ms, total: 16.3 ms
Wall time: 23.1 ms

STATUS
Optimal

plot_solution(N, x, cities_coord)

prob = LpProblem('prob', LpMinimize)

N = [('N'+str(i+1)) 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)

prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) for i in N])

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

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

#we need to keep track of the order in the tour to eliminate the possibility of subtours
u = LpVariable.dicts('u', N, 0, num_cities-1, LpInteger)

U=num_cities

for i in N:
    for j in N:
        if i != j and (i != 'N1' and j!= 'N1'):
            prob += u[i] - u[j] <= U*(1-x[i][j]) - 1

%time status = prob.solve()

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

CPU times: user 5.61 ms, sys: 8.38 ms, total: 14 ms
Wall time: 2.23 s

STATUS
Optimal

plot_solution(N, x, cities_coord)

N_unvisited = N.copy()
current = 'N1'
tour=[]
tour.append(N_unvisited.pop(N_unvisited.index(current)))

while len(N_unvisited) > 0:  
    for k in N_unvisited:
        if x[current][k].varValue ==1:
            tour.append( N_unvisited.pop(N_unvisited.index(k)))
            current=k
            break
            
tour.append('N1')

print('Optimal tour!')
print(' -> '.join(tour))

Optimal tour!
N1 -> N13 -> N12 -> N14 -> N9 -> N10 -> N7 -> N3 -> N15 -> N2 -> N11 -> N6 -> N8 -> N5 -> N4 -> N1