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.

  1. Brief Introduction
  2. Subtour Elimination
  3. Implementation

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:

  1. enter each city exactly once.
  2. leave each city excatly once.
  3. 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]