11.8 Travelling Salesman Problem (TSP)¶
The Travelling Salesman Problem is one of the most famous and studied problems in combinatorics and integer optimization. In this case study we shall:
show how to compactly define a model with Fusion;
implement an iterative algorithm that solves a sequence of optimization problems;
modify an optimization problem by adding more constraints;
show how to access the solution of an optimization problem.
The material presented in this section draws inspiration from [Pat03].
In a TSP instance we are given a directed graph \(G=(N,A)\), where \(N\) is the set of nodes and \(A\) is the set of arcs. To each arc \((i,j)\in A\) corresponds a nonnegative cost \(c_{ij}\). The goal is to find a minimum cost Hamilton cycle in \(G\), that is a closed tour passing through each node exactly once. For example, consider the small directed graph in Fig. 11.9.
Fig. 11.9 (Left) a directed graph with costs. (Middle) The minimum cycle cover found in the first iteration. (Right) The minimum cost travelling salesman tour.¶
Its corresponding adjacency and cost matrices \(A\) and \(c\) are:
Typically, the problem is modeled introducing a set of binary variables \(x_{ij}\) such that
Now we can introduce the following simple model:
It describes the constraint that every vertex has exactly one incoming and one outgoing arc in the tour, and that only arcs present in the graph can be chosen. Problem (11.31) can be easily implemented in Fusion:
Model M = new Model();
Variable x = M.Variable(new int[] {n, n}, Domain.Binary());
M.Constraint(Expr.Sum(x, 0), Domain.EqualsTo(1.0));
M.Constraint(Expr.Sum(x, 1), Domain.EqualsTo(1.0));
M.Constraint(x, Domain.LessThan( A ));
M.Objective(ObjectiveSense.Minimize, Expr.Dot(C, x));
Note in particular how:
we can sum over rows and/or columns using the
Expr.Sumfunction;we use
Expr.Dotto compute the objective function.
The solution to problem (11.31) is not necessarily a closed tour. In fact (11.31) models another problem known as minimum cost cycle cover, whose solution may consist of more than one cycle. In our example we get the solution depicted in Fig. 11.9, i.e. there are two loops, namely 0->3->0 and 1->2->1.
A solution to (11.31) solves the TSP problem if and only if it consists of a single cycle. One classical approach ensuring this is the so-called subtour elimination: once we found a solution of (11.31) composed of at least two cycles, we add constraints that explicitly avoid that particular solution:
Thus the problem we want to solve at each step is
where \(C\) is the set of cycles in all the cycle covers we have seen so far. The overall solution scheme is the following:
set \(C\) as the empty set,
solve problem (11.33),
if \(x\) has only one cycle stop,
else add the cycles of \(x\) to \(C\) and goto 2.
Cycle detection is a fairly easy task and we omit the procedure here for the sake of simplicity. Now we show how to add a constraint for each cycle. Since we have the list of arcs, and each one corresponds to a variable \(x_{ij}\), we can use the function Variable.Pick to compactly define constraints of the form (11.32):
foreach (ArrayList c in cycles) {
int [,] cc = new int[c.Count, 2];
int i = 0;
foreach (int[] a in c) {
cc[i, 0] = a[0];
cc[i, 1] = a[1];
i++;
}
M.Constraint(Expr.Sum(x.Pick(cc)), Domain.LessThan( 1.0 * c.Count - 1 ));
}
Executing our procedure will yield the following output:
it #1 - solution cost: 2.200000
cycles:
[0,3] - [3,0] -
[1,2] - [2,1] -
it #2 - solution cost: 4.000000
cycles:
[0,1] - [1,2] - [2,3] - [3,0] -
solution:
0 1 0 0
0 0 1 0
0 0 0 1
1 0 0 0
Thus we first discover the two-cycle solution; then the second iteration is forced not to include those cycles, and a new solution is located. This time it consists of one loop, and as expected the cost is higher. The solution is depicted in Fig. 11.9.
Formulation (11.33) can be improved in some cases by exploiting the graph structure. Some simple tricks follow.
Self-loops
Self-loops are never part of a TSP tour. Typically self-loops are removed by penalizing them with a huge cost \(c_{ii}\). Although this works in practice, it is more advisable to just fix the corresponding variables to zero, i.e.
This removes redundant variables, and avoids unnecessarily large coefficients that can negatively affect the solver.
Constraints (11.34) are easily implemented as follows:
M.Constraint(x.Diag(), Domain.EqualsTo(0.0));
Two-arc loops removal
In networks with more than two nodes two-loop arcs can also be ignored. They are simple to detect and their number is of the same order as the size of the graph. The constraints we need to add are:
Constraints (11.35) are easily implemented as follows:
M.Constraint(Expr.Add(x, x.Transpose()), Domain.LessThan(1.0));
The complete working example
public class tsp {
public static void Main(string [] args) {
int []A_i = new int[] {0, 1, 2, 3, 1, 0, 2, 0};
int []A_j = new int[] {1, 2, 3, 0, 0, 2, 1, 3};
double []C_v = new double[] {1.0, 1.0, 1.0, 1.0, 0.1, 0.1, 0.1, 0.1};
int n = 4;
Matrix adj = Matrix.Sparse(n, n, A_i, A_j, 1.0);
Matrix costs = Matrix.Sparse(n, n, A_i, A_j, C_v);
tsp_fusion(n, adj, costs , true, true);
tsp_fusion(n, adj, costs , true, false);
}
public static void tsp_fusion(int n, Matrix A, Matrix C,
bool remove_loops,
bool remove_2_hop_loops) {
Model M = new Model();
Variable x = M.Variable(new int[] {n, n}, Domain.Binary());
M.Constraint(Expr.Sum(x, 0), Domain.EqualsTo(1.0));
M.Constraint(Expr.Sum(x, 1), Domain.EqualsTo(1.0));
M.Constraint(x, Domain.LessThan( A ));
M.Objective(ObjectiveSense.Minimize, Expr.Dot(C, x));
if (remove_2_hop_loops)
M.Constraint(Expr.Add(x, x.Transpose()), Domain.LessThan(1.0));
if (remove_loops)
M.Constraint(x.Diag(), Domain.EqualsTo(0.0));
for (int it = 0; true; it++)
{
M.Solve();
try
{
Console.WriteLine("\nit #{0} - solution cost: {1}", it, M.PrimalObjValue());
} catch (SolutionError e) {}
ArrayList cycles = new ArrayList();
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
try {
if ( x.Level()[i * n + j] <= 0.5) continue;
} catch (SolutionError e) {}
bool found = false;
foreach (ArrayList c in cycles) {
foreach (int[] a in c)
if ( i == a[0] || j == a[0] ||
i == a[1] || j == a[1] )
{
c.Add( new int[] {i, j} );
found = true;
break;
}
if (found == true) break;
}
if (found != true) {
ArrayList ll = new ArrayList();
ll.Add( new int[] {i, j} );
cycles.Add(ll);
}
}
}
Console.WriteLine("\ncycles:");
foreach (ArrayList c in cycles) {
foreach (int[] a in c)
Console.Write("[{0},{1}] - ", a[0], a[1]);
Console.Write("\n");
}
if (cycles.Count == 1) break;
foreach (ArrayList c in cycles) {
int [,] cc = new int[c.Count, 2];
int i = 0;
foreach (int[] a in c) {
cc[i, 0] = a[0];
cc[i, 1] = a[1];
i++;
}
M.Constraint(Expr.Sum(x.Pick(cc)), Domain.LessThan( 1.0 * c.Count - 1 ));
}
}
try {
Console.WriteLine("\nsolution:");
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
Console.Write(" {0} ", ((int) x.Level()[i * n + j]));
Console.Write("\n");
}
} catch (SolutionError e) {}
}
}
