# 10.9 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. 9.

Its corresponding adjacency and cost matrices \(A\) and \(c\) are:

Typically, the problem is modelled 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 (1) can be easily implemented in *Fusion*:

```
Model::t M = new Model();
auto M_ = finally([&]() { M->dispose(); });
auto x = M->variable(new_array_ptr<int, 1>({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.sum`

function; - we use
`Expr.dot`

to compute the objective function.

The solution to problem (1) is not necessarily a closed tour. In fact (1) 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. 9, i.e. there are two loops, namely 0->3->0 and 1->2->1.

A solution to (1) 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 (1) 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 (3),
**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 (2):

```
for (auto c : cycles)
{
int csize = c.size();
auto tmp = std::shared_ptr<monty::ndarray<int, 2> >(new ndarray<int, 2>( shape(csize, 2)) );
for (auto i = 0; i < csize; ++i)
{
(*tmp)(i, 0) = std::get<0>(c[i]);
(*tmp)(i, 1) = std::get<1>(c[i]);
}
M->constraint(Expr::sum(x->pick(tmp)), Domain::lessThan( 1.0 * csize - 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. 9.

Formulation (3) 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 (4) are easily implemented as follows:

```
M->constraint(x->diag(), Domain::equalsTo(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 (5) are easily implemented as follows:

```
M->constraint(Expr::add(x, x->transpose()), Domain::lessThan(1.0));
```

The complete working example

```
template<typename MT>
void tsp(int n, MT A, MT C, bool remove_1_hop_loops, bool remove_2_hop_loops)
{
Model::t M = new Model();
auto M_ = finally([&]() { M->dispose(); });
auto x = M->variable(new_array_ptr<int, 1>({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_1_hop_loops)
M->constraint(x->diag(), Domain::equalsTo(0.));
if (remove_2_hop_loops)
M->constraint(Expr::add(x, x->transpose()), Domain::lessThan(1.0));
int it = 0;
while (true)
{
M->solve();
it++;
typedef std::vector< std::tuple<int, int> > cycle_t;
std::list< cycle_t > cycles;
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
{
if ( (*(x->level()))[i * n + j] <= 0.5 )
continue;
bool found = false;
for (auto && c : cycles)
{
for (auto && cc : c)
{
if ( i == std::get<0>(cc) || i == std::get<1>(cc) ||
j == std::get<0>(cc) || j == std::get<1>(cc) )
{
c.push_back( std::make_tuple(i, j) );
found = true;
break;
}
}
if (found) break;
}
if (!found)
cycles.push_back(cycle_t(1, std::make_tuple(i, j)));
}
std::cout << "Iteration " << it << "\n";
for (auto c : cycles) {
for (auto cc : c)
std::cout << "(" << std::get<0>(cc) << "," << std::get<1>(cc) << ") ";
std::cout << "\n";
}
if (cycles.size() == 1) break;
for (auto c : cycles)
{
int csize = c.size();
auto tmp = std::shared_ptr<monty::ndarray<int, 2> >(new ndarray<int, 2>( shape(csize, 2)) );
for (auto i = 0; i < csize; ++i)
{
(*tmp)(i, 0) = std::get<0>(c[i]);
(*tmp)(i, 1) = std::get<1>(c[i]);
}
M->constraint(Expr::sum(x->pick(tmp)), Domain::lessThan( 1.0 * csize - 1 ));
}
}
try {
std::cout << "Solution\n";
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
std::cout << (int) (*(x->level()))(i * n + j);
std::cout << "\n";
}
} catch (...) {}
}
int main()
{
auto A_i = new_array_ptr<int, 1>({0, 1, 2, 3, 1, 0, 2, 0});
auto A_j = new_array_ptr<int, 1>({1, 2, 3, 0, 0, 2, 1, 3});
auto C_v = new_array_ptr<double, 1>({1., 1., 1., 1., 0.1, 0.1, 0.1, 0.1});
int n = 4;
tsp(n, Matrix::sparse(n, n, A_i, A_j, 1.), Matrix::sparse(n, n, A_i, A_j, C_v), true, false);
tsp(n, Matrix::sparse(n, n, A_i, A_j, 1.), Matrix::sparse(n, n, A_i, A_j, C_v), true, true);
return 0;
}
```