Transport Analytics Training Series - Last Revision: October 2022
Transport Analytics Training Series - Last Revision: October 2022
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')
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
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')
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
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')
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
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