The Traveling Salesman Problem
Updated:
The Traveling Salesman Problem (TSP), a captivating conundrum in mathematical optimization, has seamlessly integrated itself into a multitude of industries, reshaping the way we approach efficiency and problem-solving.
The problem consists of finding the path that minizes the overall traveled distance between locations such that all locations are visited only once.
In the realm of logistics, the TSP takes center stage, offering a beacon of optimization for supply chains and delivery routes. Giants like Amazon leverage TSP-inspired algorithms to orchestrate last-mile deliveries, aligning packages with real-time traffic data and delivery priorities. This dynamic approach shaves miles off routes, minimizes fuel consumption, and ensures prompt deliveries.
Manufacturing, too, bows to the TSP’s prowess. Within bustling factories, where precision and organization reign supreme, the TSP guides the sequencing of production steps. Whether in the assembly of intricate automobiles or the manufacturing of various goods, the TSP minimizes idle time, reduces bottlenecks, and ultimately enhances productivity.
Surprisingly, the TSP extends its influence beyond the realms of logistics and manufacturing. The world of DNA sequencing benefits from its route-optimizing abilities, accelerating genetic research by reducing sequencing time and costs. Additionally, the TSP finds its place in the intricate landscape of circuit design, optimizing signal propagation in integrated circuits and exemplifying the synergy between mathematics and engineering.
Modules
1
2
3
4
5
6
7
8
9
10
import pulp
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import math
import networkx as nx
import itertools
%load_ext watermark
%watermark -n -u -v -iv -w
1
2
3
4
5
6
7
8
9
10
11
12
Last updated: Mon Aug 28 2023
Python implementation: CPython
Python version : 3.11.1
IPython version : 8.12.0
matplotlib: 3.7.2
networkx : 3.1
numpy : 1.25.2
pulp : 2.7.0
Watermark: 2.3.1
Data
The data used for this exercise is a dictionary containing the name of the cities and distances between them.
1
2
3
4
5
6
7
8
9
10
11
12
13
cities = {
'name': []
,'distance': []
}
# We will generate synthetic coordinates for n_cities
n_cities = 4
M = np.random.rand(n_cities,n_cities)
cities['distance'] = (M + M.T)/2
for i in range(n_cities):
cities['name'].append(f'city{i}')
cities['distance'][i, i] = 0
cities
1
2
3
4
5
{'name': ['city0', 'city1', 'city2', 'city3'],
'distance': array([[0. , 0.80122558, 0.50666924, 0.13568748],
[0.80122558, 0. , 0.64279003, 0.31393952],
[0.50666924, 0.64279003, 0. , 0.51592309],
[0.13568748, 0.31393952, 0.51592309, 0. ]])}
We use networkx to plot the graph generated by the distance matrix.
1
2
3
4
5
6
fig, ax = plt.subplots(figsize=(20, 10))
G = nx.DiGraph(cities['distance'])
pos = nx.spring_layout(G)
nx.draw_networkx_nodes(G, pos, ax=ax);
nx.draw_networkx_labels(G, pos, ax=ax);
nx.draw_networkx_edges(G, pos, ax=ax, width=2);
Formulation
For each $i,j=1,\ldots,n$, let $x_{ij}$ be a binary variable defined as 1 if there exists a path between city $i$ to $j$ and 0 otherwise.
1
x = pulp.LpVariable.dicts("x", [(i, j) for i in range(len(cities['name'])) for j in range(len(cities['name'])) if i != j], cat='Binary')
Minimize the distance $d_{ij}$ between cities $i$ and $j$
\[\min\sum_{i=1}^n\sum_{j\neq i,j=1}^nd_{ij}x_{ij}\]1
2
3
4
5
# Define the TSP problem
prob = pulp.LpProblem("TSP", pulp.LpMinimize)
# Define the objective function
prob += pulp.lpSum([cities['distance'][i, j] * x[(i, j)] for i in range(len(cities['name'])) for j in range(len(cities['name'])) if i != j])
For each city $j$, the salesman must arrive exactly one time
\[\sum_{i=1,i\neq j}^n x_{ij}=1\quad j=1,\ldots,n\]For each city $i$, the salesman must leave exactly one time:
\[\sum_{j=1,j\neq i}^n x_{ij}=1\quad i=1,\ldots,n\]1
2
3
for i in range(len(cities['name'])):
prob += pulp.lpSum([x[(i, j)] for j in range(len(cities['name'])) if i != j]) == 1
prob += pulp.lpSum([x[(j, i)] for j in range(len(cities['name'])) if i != j]) == 1
Subtour elimination constraint—ensures no proper subset $Q$ can form a sub-tour, so the solution returned is a single tour and not the union of smaller tours
\[\sum_{i\in Q}\sum_{j\neq i,j\in Q}^n x_{ij}\leq|Q|-1\quad \forall Q\subsetneq\{1,\ldots,n\},|Q|\geq2\]1
2
3
4
for k in range(len(cities['name'])):
for S in range(2, len(cities['name'])):
for subset in itertools.combinations([i for i in range(len(cities['name'])) if i != k], S):
prob += pulp.lpSum([x[(i, j)] for i in subset for j in subset if i != j]) <= len(subset) - 1
Solver
Here, we used the PuLP solver to obtain a solution to the formulated problem.
1
2
3
4
5
6
7
8
# Solve the problem using the CBC solver
prob.solve(pulp.PULP_CBC_CMD())
# Print the status of the solution
print("Status:", pulp.LpStatus[prob.status])
# Print the optimal objective value
print("Total distance traveled:", pulp.value(prob.objective))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
command line - /home/docker.datascience/.pyenv/versions/3.11.1/lib/python3.11/site-packages/pulp/solverdir/cbc/linux/64/cbc /tmp/111feaee2b8746bebc1aed7c856d7248-pulp.mps timeMode elapsed branch printingOptions all solution /tmp/111feaee2b8746bebc1aed7c856d7248-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 29 COLUMNS
At line 138 RHS
At line 163 BOUNDS
At line 176 ENDATA
Problem MODEL has 24 rows, 12 columns and 72 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 1.59909 - 0.02 seconds
Cgl0004I processed model has 18 rows, 12 columns (12 integer (12 of which binary)) and 60 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 1.59909
Cbc0038I Before mini branch and bound, 12 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.06 seconds)
Cbc0038I After 0.06 seconds - Feasibility pump exiting with objective of 1.59909 - took 0.00 seconds
Cbc0012I Integer solution of 1.5990863 found by feasibility pump after 0 iterations and 0 nodes (0.06 seconds)
Cbc0001I Search completed - best objective 1.5990862673485, took 0 iterations and 0 nodes (0.06 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 1.59909 to 1.59909
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: 1.59908627
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.07
Time (Wallclock seconds): 0.07
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.08 (Wallclock seconds): 0.08
Status: Optimal
Total distance traveled: 1.5990862673485606
Solution Analysis
Let’s understand what the solution is.
The minimal tour is the minimal distance necessary to visit all the nodes once.
1
print("minimal tour: ", prob.objective.value())
1
minimal tour: 1.5990862673485606
The followin code extracts the route of the optimal solution.
1
2
3
4
5
6
7
8
9
10
11
12
# Extract the solution
solution = []
start_city = 0
next_city = start_city
while True:
for j in range(len(cities['name'])):
if j != next_city and x[(next_city, j)].value() == 1:
solution.append((next_city, j))
next_city = j
break
if next_city == start_city:
break
1
2
3
4
# Print the solution
print("Route:")
for i in range(len(solution)):
print(str(solution[i][0]) + " -> " + str(solution[i][1]))
1
2
3
4
5
Route:
0 -> 2
2 -> 1
1 -> 3
3 -> 0
Plot Solution
Here, we can plot the solution as a graph. The optimal solution gives the distance from city $i$ (row) to city $j$ (column).
1
2
3
4
sol_matrix = np.zeros((len(solution), len(solution)))
for i in range(len(solution)):
sol_matrix[solution[i]] = cities['distance'][solution[i]]
sol_matrix
1
2
3
4
array([[0. , 0. , 0.50666924, 0. ],
[0. , 0. , 0. , 0.31393952],
[0. , 0.64279003, 0. , 0. ],
[0.13568748, 0. , 0. , 0. ]])
The following function allow us to draw this graph using the distances as weights on the vertices.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
class Draw_Graph():
def __init__(self, G, pos=None, arc: float = 0.25, edge_attribute: str='weight'):
self.G = G
self.arc = arc
self.edge_attribute = edge_attribute
self.pos = pos
def draw_edges(self):
self.curved_edges = [edge for edge in self.G.edges() if reversed(edge) in self.G.edges()]
self.straight_edges = list(set(self.G.edges()) - set(self.curved_edges))
nx.draw_networkx_edges(self.G, self.pos, ax=ax, edgelist=self.straight_edges, width=2)
nx.draw_networkx_edges(self.G, self.pos, ax=ax, edgelist=self.curved_edges, connectionstyle=f'arc3, rad = {self.arc}', width=2)
def draw_labels(self):
self.edge_weights = nx.get_edge_attributes(self.G, self.edge_attribute)
self.curved_edge_labels = {edge: self.edge_weights[edge] for edge in self.curved_edges}
self.straight_edge_labels = {edge: self.edge_weights[edge] for edge in self.straight_edges}
self.draw_networkx_edge_labels(self.G, self.pos, ax=ax, edge_labels=self.curved_edge_labels, rotate=False,rad = self.arc)
nx.draw_networkx_edge_labels(self.G, self.pos, ax=ax, edge_labels=self.straight_edge_labels, rotate=False)
@staticmethod
def draw_networkx_edge_labels(
G,
pos,
edge_labels=None,
label_pos=0.5,
font_size=10,
font_color="k",
font_family="sans-serif",
font_weight="normal",
alpha=None,
bbox=None,
horizontalalignment="center",
verticalalignment="center",
ax=None,
rotate=True,
clip_on=True,
rad=0
):
"""Draw edge labels.
Parameters
----------
G : graph
A networkx graph
pos : dictionary
A dictionary with nodes as keys and positions as values.
Positions should be sequences of length 2.
edge_labels : dictionary (default={})
Edge labels in a dictionary of labels keyed by edge two-tuple.
Only labels for the keys in the dictionary are drawn.
label_pos : float (default=0.5)
Position of edge label along edge (0=head, 0.5=center, 1=tail)
font_size : int (default=10)
Font size for text labels
font_color : string (default='k' black)
Font color string
font_weight : string (default='normal')
Font weight
font_family : string (default='sans-serif')
Font family
alpha : float or None (default=None)
The text transparency
bbox : Matplotlib bbox, optional
Specify text box properties (e.g. shape, color etc.) for edge labels.
Default is {boxstyle='round', ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0)}.
horizontalalignment : string (default='center')
Horizontal alignment {'center', 'right', 'left'}
verticalalignment : string (default='center')
Vertical alignment {'center', 'top', 'bottom', 'baseline', 'center_baseline'}
ax : Matplotlib Axes object, optional
Draw the graph in the specified Matplotlib axes.
rotate : bool (deafult=True)
Rotate edge labels to lie parallel to edges
clip_on : bool (default=True)
Turn on clipping of edge labels at axis boundaries
Returns
-------
dict
`dict` of labels keyed by edge
Examples
--------
>>> G = nx.dodecahedral_graph()
>>> edge_labels = nx.draw_networkx_edge_labels(G, pos=nx.spring_layout(G))
Also see the NetworkX drawing examples at
https://networkx.org/documentation/latest/auto_examples/index.html
See Also
--------
draw
draw_networkx
draw_networkx_nodes
draw_networkx_edges
draw_networkx_labels
"""
if ax is None:
ax = plt.gca()
if edge_labels is None:
labels = {(u, v): d for u, v, d in G.edges(data=True)}
else:
labels = edge_labels
text_items = {}
for (n1, n2), label in labels.items():
(x1, y1) = pos[n1]
(x2, y2) = pos[n2]
(x, y) = (
x1 * label_pos + x2 * (1.0 - label_pos),
y1 * label_pos + y2 * (1.0 - label_pos),
)
pos_1 = ax.transData.transform(np.array(pos[n1]))
pos_2 = ax.transData.transform(np.array(pos[n2]))
linear_mid = 0.5*pos_1 + 0.5*pos_2
d_pos = pos_2 - pos_1
rotation_matrix = np.array([(0,1), (-1,0)])
ctrl_1 = linear_mid + rad*rotation_matrix@d_pos
ctrl_mid_1 = 0.5*pos_1 + 0.5*ctrl_1
ctrl_mid_2 = 0.5*pos_2 + 0.5*ctrl_1
bezier_mid = 0.5*ctrl_mid_1 + 0.5*ctrl_mid_2
(x, y) = ax.transData.inverted().transform(bezier_mid)
if rotate:
# in degrees
angle = np.arctan2(y2 - y1, x2 - x1) / (2.0 * np.pi) * 360
# make label orientation "right-side-up"
if angle > 90:
angle -= 180
if angle < -90:
angle += 180
# transform data coordinate angle to screen coordinate angle
xy = np.array((x, y))
trans_angle = ax.transData.transform_angles(
np.array((angle,)), xy.reshape((1, 2))
)[0]
else:
trans_angle = 0.0
# use default box of white with white border
if bbox is None:
bbox = dict(boxstyle="round", ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0))
if not isinstance(label, str):
label = f"{label:0.2f}" # this makes "1" and 1 labeled the same
t = ax.text(
x,
y,
label,
size=font_size,
color=font_color,
family=font_family,
weight=font_weight,
alpha=alpha,
horizontalalignment=horizontalalignment,
verticalalignment=verticalalignment,
rotation=trans_angle,
transform=ax.transData,
bbox=bbox,
zorder=1,
clip_on=clip_on,
)
text_items[(n1, n2)] = t
ax.tick_params(
axis="both",
which="both",
bottom=False,
left=False,
labelbottom=False,
labelleft=False,
)
return text_items
def plot(self):
self.__call__()
def __call__(self):
self.draw_edges()
self.draw_labels()
The graph can be shown below
1
2
3
4
5
6
7
G = nx.DiGraph(sol_matrix)
pos = nx.spring_layout(G)
fig, ax = plt.subplots(figsize=(20, 10))
p = Draw_Graph(G, pos)
p.plot()
nx.draw_networkx_nodes(G, pos, ax=ax);
nx.draw_networkx_labels(G, pos, ax=ax);
References
[1] https://soumenatta.medium.com/solving-the-traveling-salesman-problem-using-pulp-in-python-edd23a6aee4d
[2] https://towardsdatascience.com/solving-geographic-travelling-salesman-problems-using-python-e57284b14cd7
[3] https://soumenatta.medium.com/solving-the-p-median-problem-using-pulp-in-python-31d9bc13cc2d
Comments