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

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.

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 binding and needed. The idea is to
add those potentially binding 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

This is how I have it implemented:

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 | ```
dm = (
( 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)
E = [(i,j) for i in V for j in V if i!=j]
print "there are %d cities"%n
from pymprog import *
beginModel('TSP by row generation')
x = var(E, 'x', bool)
minimize( sum(dm[i][j]*x[i,j] for i,j in E), 'TravelDist' )
st([sum( x[k,j] for j in V if j!=k ) == 1 for k in V], 'leave')
st([sum( x[i,k] for i in V if i!=k ) == 1 for k in V], 'enter')
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:", vobj()
print "Optimal tour:\n", subt
break
print "New subtour: ", 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]
st( sum(x[i,j] for i in subt for j in nots) >= 1 )
solve() #solve the IP problem again
``` |

And here is the output:

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