Transport Analytics Training Series - Last Revision: October 2022

The Travelling Salesman Problem - MTZ Model¶

image-4.png

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.

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 = 15
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');
In [5]:
from pulp import *

prob = LpProblem('prob', LpMinimize)

We are back again in using PuLP.

In [6]:
from scipy.spatial import distance

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

dist_matrix.shape
Out[6]:
(15, 15)
In [7]:
N = [('N'+str(i+1)) for i in range(num_cities)]
In [8]:
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)}
In [9]:
x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)
In [10]:
prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) 
                                      for i in N])
In [11]:
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
In [12]:
%time status = prob.solve()

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

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/pa01/opt/anaconda3/lib/python3.9/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/1ecfb06cefb24b2cb6c52e0284a6b62b-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/1ecfb06cefb24b2cb6c52e0284a6b62b-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 35 COLUMNS
At line 1086 RHS
At line 1117 BOUNDS
At line 1328 ENDATA
Problem MODEL has 30 rows, 210 columns and 420 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 283.65 - 0.00 seconds
Cgl0004I processed model has 30 rows, 210 columns (210 integer (210 of which binary)) and 420 elements
Cutoff increment increased from 1e-05 to 0.00999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 283.65
Cbc0038I Before mini branch and bound, 210 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 283.65 - took 0.00 seconds
Cbc0012I Integer solution of 283.65 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0001I Search completed - best objective 283.65, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 283.65 to 283.65
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                283.65000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.02

CPU times: user 6.36 ms, sys: 6.36 ms, total: 12.7 ms
Wall time: 177 ms

STATUS
Optimal
In [13]:
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)
In [14]:
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
In [15]:
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
In [16]:
for i in N:
    for j in N:
        if i!=j:
            prob += x[i][j] + x[j][i] <= 1
In [17]:
%time status = prob.solve()

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

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/pa01/opt/anaconda3/lib/python3.9/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/ecf92a2c7a7e4120a7f065d87136d1e5-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/ecf92a2c7a7e4120a7f065d87136d1e5-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 245 COLUMNS
At line 1716 RHS
At line 1957 BOUNDS
At line 2168 ENDATA
Problem MODEL has 240 rows, 210 columns and 840 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 294.801 - 0.00 seconds
Cgl0004I processed model has 135 rows, 210 columns (210 integer (210 of which binary)) and 630 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 294.801
Cbc0038I Before mini branch and bound, 210 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 294.801 - took 0.00 seconds
Cbc0012I Integer solution of 294.80063 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective 294.800627027689, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 294.801 to 294.801
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                294.80062703
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

CPU times: user 14 ms, sys: 31.1 ms, total: 45.1 ms
Wall time: 27.3 ms

STATUS
Optimal
In [18]:
plot_solution(N, x, cities_coord)
In [19]:
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
In [20]:
#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)
In [21]:
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
In [22]:
%time status = prob.solve()

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

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/pa01/opt/anaconda3/lib/python3.9/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/faf5162448a94ff79c45e2e0ab94ecee-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/kk/14m2mg8s0w9c_tln749srnrw0000gp/T/faf5162448a94ff79c45e2e0ab94ecee-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 217 COLUMNS
At line 1842 RHS
At line 2055 BOUNDS
At line 2280 ENDATA
Problem MODEL has 212 rows, 224 columns and 966 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 285.134 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 182 strengthened rows, 0 substitutions
Cgl0004I processed model has 212 rows, 224 columns (224 integer (210 of which binary)) and 2604 elements
Cutoff increment increased from 1e-05 to 0.00999
Cbc0038I Initial state - 18 integers unsatisfied sum - 1.71429
Cbc0038I Pass   1: suminf.    2.57011 (22) obj. 303.739 iterations 46
Cbc0038I Pass   2: suminf.    1.71429 (6) obj. 536.944 iterations 85
Cbc0038I Pass   3: suminf.    1.60000 (8) obj. 558.544 iterations 17
Cbc0038I Pass   4: suminf.    3.06893 (20) obj. 585.323 iterations 58
Cbc0038I Pass   5: suminf.    3.46589 (21) obj. 573.154 iterations 65
Cbc0038I Pass   6: suminf.    1.94286 (9) obj. 531.984 iterations 70
Cbc0038I Pass   7: suminf.    2.03218 (15) obj. 531.516 iterations 35
Cbc0038I Pass   8: suminf.    4.21212 (17) obj. 561.054 iterations 43
Cbc0038I Pass   9: suminf.    4.31233 (20) obj. 574.244 iterations 29
Cbc0038I Pass  10: suminf.    2.79592 (9) obj. 560.37 iterations 49
Cbc0038I Pass  11: suminf.    2.86667 (18) obj. 564.76 iterations 41
Cbc0038I Pass  12: suminf.    2.34286 (14) obj. 587.182 iterations 41
Cbc0038I Pass  13: suminf.    1.60000 (11) obj. 586.884 iterations 37
Cbc0038I Pass  14: suminf.    5.17333 (19) obj. 581.205 iterations 70
Cbc0038I Pass  15: suminf.    2.16364 (20) obj. 579.758 iterations 31
Cbc0038I Pass  16: suminf.    2.67816 (26) obj. 543.585 iterations 67
Cbc0038I Pass  17: suminf.    3.00000 (16) obj. 565.187 iterations 86
Cbc0038I Pass  18: suminf.    6.40000 (16) obj. 559.422 iterations 17
Cbc0038I Pass  19: suminf.    4.93333 (19) obj. 547.673 iterations 53
Cbc0038I Pass  20: suminf.    1.71429 (4) obj. 569.347 iterations 93
Cbc0038I Pass  21: suminf.    1.60000 (8) obj. 560.336 iterations 35
Cbc0038I Pass  22: suminf.    2.02727 (19) obj. 545.188 iterations 70
Cbc0038I Pass  23: suminf.    4.41395 (19) obj. 539.073 iterations 45
Cbc0038I Pass  24: suminf.    3.80000 (17) obj. 543.116 iterations 40
Cbc0038I Pass  25: suminf.    2.64762 (12) obj. 548.049 iterations 50
Cbc0038I Pass  26: suminf.    3.15107 (31) obj. 540.719 iterations 64
Cbc0038I Pass  27: suminf.    1.71429 (4) obj. 525.841 iterations 86
Cbc0038I Pass  28: suminf.    5.13325 (23) obj. 538.144 iterations 43
Cbc0038I Pass  29: suminf.    3.09333 (17) obj. 568.171 iterations 69
Cbc0038I Pass  30: suminf.    2.16889 (11) obj. 552.75 iterations 47
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 141 integers at bound fixed and 0 continuous
Cbc0038I Full problem 212 rows 224 columns, reduced to 177 rows 83 columns - too large
Cbc0038I Mini branch and bound did not improve solution (0.08 seconds)
Cbc0038I After 0.08 seconds - Feasibility pump exiting - took 0.03 seconds
Cbc0031I 12 added rows had average density of 83.833333
Cbc0013I At root node, 12 cuts changed objective from 285.24 to 318.09266 in 100 passes
Cbc0014I Cut generator 0 (Probing) - 488 row cuts average 5.3 elements, 0 column cuts (0 active)  in 0.085 seconds - new frequency is 1
Cbc0014I Cut generator 1 (Gomory) - 1445 row cuts average 194.3 elements, 0 column cuts (0 active)  in 0.050 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.019 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.002 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 16 row cuts average 173.9 elements, 0 column cuts (0 active)  in 0.031 seconds - new frequency is -100
Cbc0014I Cut generator 5 (FlowCover) - 67 row cuts average 2.0 elements, 0 column cuts (0 active)  in 0.027 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 240 row cuts average 81.8 elements, 0 column cuts (0 active)  in 0.016 seconds - new frequency is 1
Cbc0014I Cut generator 7 (ZeroHalf) - 53 row cuts average 35.9 elements, 0 column cuts (0 active)  in 0.260 seconds - new frequency is -100
Cbc0010I After 0 nodes, 1 on tree, 1e+50 best solution, best possible 318.09266 (0.80 seconds)
Cbc0016I Integer solution of 404.7 found by strong branching after 3847 iterations and 28 nodes (0.93 seconds)
Cbc0038I Full problem 212 rows 224 columns, reduced to 31 rows 14 columns
Cbc0038I Full problem 212 rows 224 columns, reduced to 163 rows 22 columns - 4 fixed gives 145, 15 - ok now
Cbc0038I Full problem 212 rows 224 columns, reduced to 24 rows 15 columns
Cbc0016I Integer solution of 396.94 found by strong branching after 4301 iterations and 46 nodes (0.98 seconds)
Cbc0038I Full problem 212 rows 224 columns, reduced to 72 rows 23 columns
Cbc0016I Integer solution of 393.68 found by strong branching after 5511 iterations and 104 nodes (1.10 seconds)
Cbc0016I Integer solution of 365.19 found by strong branching after 6173 iterations and 138 nodes (1.16 seconds)
Cbc0038I Full problem 212 rows 224 columns, reduced to 186 rows 20 columns - 3 fixed gives 174, 15 - still too large
Cbc0038I Full problem 212 rows 224 columns, reduced to 44 rows 15 columns
Cbc0016I Integer solution of 364.18 found by strong branching after 9168 iterations and 208 nodes (1.35 seconds)
Cbc0016I Integer solution of 363.31 found by strong branching after 9170 iterations and 208 nodes (1.35 seconds)
Cbc0004I Integer solution of 327.08 found after 9424 iterations and 213 nodes (1.36 seconds)
Cbc0038I Full problem 212 rows 224 columns, reduced to 181 rows 25 columns - 1 fixed gives 175, 15 - still too large
Cbc0001I Search completed - best objective 327.0799999999999, took 17280 iterations and 328 nodes (1.64 seconds)
Cbc0032I Strong branching done 4686 times (67799 iterations), fathomed 54 nodes and fixed 56 variables
Cbc0035I Maximum depth 22, 3969 variables fixed on reduced cost
Cuts at root node changed objective from 285.24 to 318.093
Probing was tried 558 times and created 856 cuts of which 0 were active after adding rounds of cuts (0.096 seconds)
Gomory was tried 557 times and created 1476 cuts of which 0 were active after adding rounds of cuts (0.090 seconds)
Knapsack was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.019 seconds)
Clique was tried 100 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.002 seconds)
MixedIntegerRounding2 was tried 100 times and created 16 cuts of which 0 were active after adding rounds of cuts (0.031 seconds)
FlowCover was tried 100 times and created 67 cuts of which 0 were active after adding rounds of cuts (0.027 seconds)
TwoMirCuts was tried 557 times and created 623 cuts of which 0 were active after adding rounds of cuts (0.065 seconds)
ZeroHalf was tried 100 times and created 53 cuts of which 0 were active after adding rounds of cuts (0.260 seconds)

Result - Optimal solution found

Objective value:                327.08000000
Enumerated nodes:               328
Total iterations:               17280
Time (CPU seconds):             1.60
Time (Wallclock seconds):       1.65

Option for printingOptions changed from normal to all
Total time (CPU seconds):       1.60   (Wallclock seconds):       1.65

CPU times: user 7.63 ms, sys: 29.2 ms, total: 36.8 ms
Wall time: 1.67 s

STATUS
Optimal
In [23]:
plot_solution(N, x, cities_coord)
In [24]:
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