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

Implementation

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]

Table Of Contents

Previous topic

Working with Solver Options

This Page