One of the difficulties we encounter in solving the Traveling Salesman Problem (TSP) is that, for even a small numer of cities, a complete description of the problem requires a factorial number of constraints. This is apparent in the standard formulation used to teach the TSP to OR students. Consider a set of $n$ cities with the distance from city $i$ to city $j$ denoted $d_{ij}$. We attempt to minimize the total distance of a tour entering and leaving each city exactly once. $x_{ij} = 1$ if the edge from city $i$ to city $j$ is included in the tour, $0$ otherwise:
$$ \small \begin{align*} \min\quad & z = \sum_i \sum_{j\ne i} d_{ij} x_{ij}\\ \text{s.t.}\quad& \sum_{j\ne i} x_{ij} = 1 &\quad\forall&\ i & \text{leave each city once}\\ & \sum_{i\ne j} x_{ij} = 1 &\quad\forall&\ j & \text{enter each city once}\\ & x_{ij} \in \{0,1\} &\quad\forall&\ i,j \end{align*} $$
This appears like a reasonable formulation until we solve it and see that our solution contains disconnected subtours. Suppose we have four cities, labeled $A$ through $D$. Connecting $A$ to $B$, $B$ to $A$, $C$ to $D$ and $D$ to $C$ provides a feasible solution to our formulation, but does not constitute a cycle. Here is a more concrete example of two disconnected subtours $\{(1,5),(5,1)\}$ and $\{(2,3),(3,4),(4,2)\}$ over five cities:
ampl: display x;
x [*,*]
: 1 2 3 4 5 :=
1 0 0 0 0 1
2 0 0 1 0 0
3 0 0 0 1 0
4 0 1 0 0 0
5 1 0 0 0 0
;
Realizing we just solved the Assignment Problem, we now add subtour elimination constraints. These require that any proper, non-null subset of our $n$ cities is connected by at most $n-1$ active edges:
$$ \sum_{i \in S} \sum_{j \in S} x_{ij} \leq |S|-1 \quad\forall\ S \subset {1, …, n}, S \ne O $$
Indexing subtour elimination constraints over a power set of the cities completes the formulation. However, this requires an additional $\sum_{k=2}^{n-1} \begin{pmatrix} n \\ k \end{pmatrix}$ rows tacked on the end of our matrix and is clearly infeasible for large $n$. The most current computers can handle using this approach is around 19 cities. It remains an instructive tool for understanding the combinatorial explosion that occurs in problems like TSP and is worth translating into a modeling language. So how does one model it on a computer?
Unfortunately, AMPL, the gold standard in mathematical modeling languages, is unable to index over sets. Creating a power set in AMPL requires going through a few contortions. The following code demonstrates power and index sets over four cities:
set cities := 1 .. 4 ordered;
param n := card(cities);
set indices := 0 .. (2^n - 1);
set power {i in indices} := {c in cities: (i div 2^(ord(c) - 1)) mod 2 = 1};
display cities;
display n;
display indices;
display power;
This yields the following output:
set cities := 1 2 3 4;
n = 4
set indices := 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15;
set power[0] := ; # empty
set power[1] := 1;
set power[2] := 2;
set power[3] := 1 2;
set power[4] := 3;
set power[5] := 1 3;
set power[6] := 2 3;
set power[7] := 1 2 3;
set power[8] := 4;
set power[9] := 1 4;
set power[10] := 2 4;
set power[11] := 1 2 4;
set power[12] := 3 4;
set power[13] := 1 3 4;
set power[14] := 2 3 4;
set power[15] := 1 2 3 4;
Note how the index set contains an index for each row in our power set. We can now generate the subtour elimination constraints:
var x {cities cross cities} binary;
s.t. subtours {i in indices: card(power[i]) > 1 and card(power[i]) < card(cities)}:
sum {(c,k) in power[i] cross power[i]: k != c} x[c,k] <= card(power[i]) - 1;
expand subtours;
subject to subtours[3]: x[1,2] + x[2,1] <= 1;
subject to subtours[5]: x[1,3] + x[3,1] <= 1;
subject to subtours[6]: x[2,3] + x[3,2] <= 1;
subject to subtours[7]: x[1,2] + x[1,3] + x[2,1] + x[2,3] + x[3,1] + x[3,2] <= 2;
subject to subtours[9]: x[1,4] + x[4,1] <= 1;
subject to subtours[10]: x[2,4] + x[4,2] <= 1;
subject to subtours[11]: x[1,2] + x[1,4] + x[2,1] + x[2,4] + x[4,1] + x[4,2] <= 2;
subject to subtours[12]: x[3,4] + x[4,3] <= 1;
subject to subtours[13]: x[1,3] + x[1,4] + x[3,1] + x[3,4] + x[4,1] + x[4,3] <= 2;
subject to subtours[14]: x[2,3] + x[2,4] + x[3,2] + x[3,4] + x[4,2] + x[4,3] <= 2;
While this does work, the code for generating the power set looks like voodoo. Understanding it required piece-by-piece decomposition, an exercise I suggest you go through yourself if you have a copy of AMPL and 15 minutes to spare:
set foo {c in cities} := {ord(c)};
set bar {c in cities} := {2^(ord(c) - 1)};
set baz {i in indices} := {c in cities: i div 2^(ord(c) - 1)};
set qux {i in indices} := {c in cities: (i div 2^(ord(c) - 1)) mod 2 = 1};
display foo;
display bar;
display baz;
display qux;
This may be an instance where open source leads commercial software. The good folks who produce the SCIP Optimization Suite provide an AMPL-like language called ZIMPL with a few additional useful features. One of these is power sets. Compared to the code above, doesn’t this look refreshing?
set cities := {1 to 4};
set power[] := powerset(cities);
set indices := indexset(power);