Notice: Trying to get property of non-object in /var/www/examples/traveling-salesman-problem/index.php on line 51

Notice: Trying to get property of non-object in /var/www/examples/traveling-salesman-problem/index.php on line 80

Notice: Trying to get property of non-object in /var/www/examples/traveling-salesman-problem/index.php on line 94

The Traveling Salesman Problem

with integer programming and Gurobi

In this example we'll solve the Traveling Salesman Problem.

We'll construct a mathematical model of the problem, implement this model in Gurobi's Python interface, and compute and visualize an optimal solution.

Although your own business may not involve traveling salesmen, the same basic techniques used in this example can be used for many other applications like vehicle rounting, circuit design and DNA sequencing.

Problem Description

The Traveling Salesman Problem (TSP) is a classic problem in combinatorial optimization. It was first formulated as an integer program by Dantzig, Fulkerson and Johnson in 1954.

In this example, we consider a salesman traveling in the US. The salesman starts in New York and has to visit a set of cities on a business trip before returning home. The problem then consists of finding the shortest tour which visits every city on the itinerary.

We will formulate this problem as an integer program and implement it in Gurobi. The implementation will also demonstrate the use of lazy constraints in Gurobi.

Mathematical Model

We let the $n$ selected cities in the salesman's tour be the set of vertices $V$ of a graph. The set of edges $E$ of the graph corresponds to the different connections between each city. Since we can travel from any city to another, the graph is complete. That is, there is an edge between every pair of nodes. For each edge in the the graph we associate a binary variable $x_{ij} = \left\{\begin{array}{ll} 1 & \text{if edge (i,j) \in E is in tour }\\ 0 & \mathrm{otherwise.} \end{array}\right.$ Since the edges are undirected, we have that $x_{ij} = x_{ji}$, and it suffices to only include edges with $i < j$ in the model.

We want to minimize the total distance travelled during the tour. Therefore, we calculate the distance $d_{ij}$ between each pair of nodes $i$ and $j$. The total distance travelled is then the sum of the distances of the edges included in the tour $\text{distance travelled} = \sum_{(i,j) \in E} d_{ij} x_{ij}.$

The tour should only pass through each city once. Therefore, each node in the graph should have exactly one incoming edge and one outgoing edge. In other words, for every node $i$ exactly two of the $x_{ij}$ binary variables should be equal to 1. We write this constraint as $\sum_{j \in V} x_{ij} = 2 \quad \forall i \in V.$ This constraint means that the saleman should enter and leave each city exactly once.

With just this constraint on the number of edges in the tour entering and exit each node, we may produce solutions that are not connected tours. For example, the figure below shows a graph with 6 nodes and two disconnected subtours. The first subtour goes through nodes 1, 2, and 3, while the second subtour goes through nodes 4, 5, 6, and 7. Note that in this solution each node has exactly two edges in the tour incident to it. But there is no path for the salesman to travel between subtour. So we must add extra constraints to our model to eliminate these solutions.

To eliminate the subtours we add the following constraints: $\sum_{i,j \in S, \, i \neq j} x_{ij} \leq \left\vert{S}\right\vert - 1, \quad \forall S \subset V, S \ne \emptyset$ These constraints require that for each proper (nonempty) subset $S$ of the set of cities $V$, the number of edges between the nodes of $S$ must be at most $\left\vert{S}\right\vert - 1$.

Indeed if the number of edges were equal to $S$, it would be possible to form a subtour. For example, in the figure above the subset of nodes $S=\{1,2,3\}$ has 3 edges in the tour: $x_{13} = 1, x_{12} = 1, x_{23} = 1$. So $\sum_{i,j \in \{1,2,3\}, i \ne j} x_{ij} = 3 > 2 = |\{1, 2, 3\}| - 1$ Thus, the subtour elimination constraint above is violated.

So finally the integer program formulation becomes $\begin{array}{ll} \text{minimize} & \sum_{(i,j) \in E} d_{ij} x_{ij}} \\ \text{subject to} & \sum_{j \in V} x_{ij}} = 2 \quad \forall i \in V \\ & \sum_{i,j \in S, \, i \neq j} x_{ij}} \leq \left\vert{S}\right\vert - 1 \quad \forall S \subset V, S \ne \emptyset \\ & x_{ij} \in \{ 0, 1 \} \end{array}$

If the set of cities $V$ is of size $n$, there are $2^n - 2$ subsets $S$ of $V$, excluding $S = V$ and $S = \emptyset$. Instead of explicitally including a constraint $\sum_{i,j \in S, \, i \neq j} x_{ij} \leq \left\vert{S}\right\vert - 1$ for each $S$. We include above constraints implicitly as lazy constraints. That is, we generate and add these constraints to our model in a lazy fashion.

Initially we will have no subtour elimination constraints in our model. When Gurobi finds a feasible solution that satisifies the other constraints, we compute the shortest cycle in edges included in the tour so far. If this cycle is of length $n$, then the model is solved. Otherwise, a cycle of length $m < n$ defines a subtour, and a set $S$ with $|S| = m$. We then add the corresponding constraint to eliminate this subtour as a lazy constraint, and Gurobi continues solving this new modified model. This process continues until the shortest cycle is of length $n$, implying that all subtour elimination constraints have been satisfied.

Implementation

Below is the full implementation of the model (and the associated data) in Gurobi's Python interface:

import math import random from gurobipy import * # Callback - use lazy constraints to eliminate sub-tours def subtourelim(model, where): if where == GRB.callback.MIPSOL: selected = [] # make a list of edges selected in the solution for i in range(n): sol = model.cbGetSolution([model._vars[i,j] for j in range(n)]) selected += [(i,j) for j in range(n) if sol[j] > 0.5] # find the shortest cycle in the selected edge list tour = subtour(selected) if len(tour) < n: # add a subtour elimination constraint expr = 0 for i in range(len(tour)): for j in range(i+1, len(tour)): expr += model._vars[tour[i], tour[j]] model.cbLazy(expr <= len(tour)-1) # Euclidean distance between two points def distance(points, i, j): dx = points[i] - points[j] dy = points[i] - points[j] return math.sqrt(dx*dx + dy*dy) # Given a list of edges, finds the shortest subtour def subtour(edges): visited = [False]*n cycles = [] lengths = [] selected = [[] for i in range(n)] for x,y in edges: selected[x].append(y) while True: current = visited.index(False) thiscycle = [current] while True: visited[current] = True neighbors = [x for x in selected[current] if not visited[x]] if len(neighbors) == 0: break current = neighbors thiscycle.append(current) cycles.append(thiscycle) lengths.append(len(thiscycle)) if sum(lengths) == n: break return cycles[lengths.index(min(lengths))] n = 50 # Create n random points random.seed(1) points = [] for i in range(n): points.append((random.randint(0,100),random.randint(0,100))) m = Model() # Create variables vars = {} for i in range(n): for j in range(i+1): vars[i,j] = m.addVar(obj=distance(points, i, j), vtype=GRB.BINARY, name='e'+str(i)+'_'+str(j)) vars[j,i] = vars[i,j] m.update() # Add degree-2 constraint, and forbid loops for i in range(n): m.addConstr(quicksum(vars[i,j] for j in range(n)) == 2) vars[i,i].ub = 0 m.update() # Optimize model m._vars = vars m.params.LazyConstraints = 1 m.optimize(subtourelim) solution = m.getAttr('x', vars) selected = [(i,j) for i in range(n) for j in range(n) if solution[i,j] > 0.5] assert len(subtour(selected)) == n

Live Demo

The cities which are part of the tour have been highlighted: New York, Houston, San Francisco, Seattle, Minneapolis and Denver. Click "Compute Tour" to find the optimal tour using Gurobi.

You can hover over cities to show their names and click to add and remove them from the tour.