# Case: Subtour Elimination¶

This example is adapted from the presentation
titled *Mathematical Programming: Modeling and
Applications* by Giacomo Nannicini. Thanks are
hereby extended to him.

## Brief Introduction¶

The Traveling Salesman Problem (TSP) is a very well known problem in the literature. Applications of TSP include: logistics, crane control, placing circuits on a board minimizing the required time, and many more. Unfortunately, it is a very difficult problem. For not too large instances, it can be done on a desktop machine.

Here is a definition of a TSP problem: A salesman must visit all cities to see his customers, and return to the starting point. He wants to minimize the total travel distance. Here are are going to play with a small example of TSP, assuming that the distance between any two cities is symmetric.

## Subtour Elimination¶

A subtour is also a round tour that returns back to where you start, but does not visit all the cities. A formulation of TSP is this:

- enter each city exactly once.
- leave each city excatly once.
- make sure there is no subtour.

To make sure there is no subtour, we must consider
*all* subset of cities, and make sure that there
is an arc leaving a city in the subset and entering
a city NOT in the subset. So there are exponential
number of subtour elimination constraints.
Obviously, only a small number of them will be
actually needed to eliminate subtours.
The idea is to start out without them and then
add those violated ones gradually,
until the solution contains no subtour.
For a more detailed discussion on TSP, please see
http://www.tsp.gatech.edu/methods/opt/subtour.htm

## Implementation¶

This is how I have it implemented using PyMathProg:

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 | ```
dm = ( # distance matrix
( 0,86,49,57,31,69,50),
(86, 0,68,79,93,24, 5),
(49,68, 0,16, 7,72,67),
(57,79,16, 0,90,69, 1),
(31,93, 7,90, 0,86,59),
(69,24,72,69,86, 0,81),
(50, 5,67, 1,59,81, 0))
n = len(dm) #how many cities
V = range(n) # set of cities
E = [(i,j) for i in V for j in V if i!=j]
print("there are %d cities"%n)
from pymprog import *
begin('subtour elimination')
x = var('x', E, bool)
minimize(sum(dm[i][j]*x[i,j] for i,j in E), 'dist')
for k in V:
sum( x[k,j] for j in V if j!=k ) == 1
sum( x[i,k] for i in V if i!=k ) == 1
solver(float, msg_lev=glpk.GLP_MSG_OFF)
solver(int, msg_lev=glpk.GLP_MSG_OFF)
solve() #solve the IP problem
def subtour(x):
"find a subtour in current solution"
succ = 0
subt = [succ] #start from node 0
while True:
succ=sum(x[succ,j].primal*j
for j in V if j!=succ)
if succ == 0: break #tour found
subt.append(int(succ+0.5))
return subt
while True:
subt = subtour(x)
if len(subt) == n:
print("Optimal tour length: %g"%vobj())
print("Optimal tour:"); print(subt)
break
print("New subtour: %r"% subt)
if len(subt) == 1: break #something wrong
#now add a subtour elimination constraint:
nots = [j for j in V if j not in subt]
sum(x[i,j] for i in subt for j in nots) >= 1
solve() #solve the IP problem again
end()
``` |

And here is the output:

```
there are 7 cities
New subtour: [0, 4, 2]
New subtour: [0, 6, 3, 2, 4]
Optimal tour length: 153.0
Optimal tour:
[0, 5, 1, 6, 3, 2, 4]
```