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:
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]