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

Its corresponding adjacency and cost matrices $$A$$ and $$c$$ are:

$\begin{split}\begin{array}{cc} A= \left[ \begin{array}{cccc} 0 & 1 & 1 & 1\\ 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1\\ 1 & 0 & 0 & 0 \end{array} \right], & c= \left[ \begin{array}{cccc} 0 & 1 & 0.1 & 0.1\\ 0.1 & 0 & 1 & 0\\ 0 & 0.1 & 0 & 1\\ 1 & 0 & 0 & 0 \end{array} \right]. \end{array}\end{split}$

Typically, the problem is modeled introducing a set of binary variables $$x_{ij}$$ such that

$\begin{split}x_{ij} = \left\lbrace \begin{array}{cl} 0 & \text{if arc } (i,j) \text{ is in the tour}, \\ 1 & \text{otherwise}. \end{array}\right.\end{split}$

Now we can introduce the following simple model:

(11.31)$\begin{split}\begin{array}{lll} \min & \sum_{i,j} c_{ij}x_{ij} &\\ \st & \sum_i x_{ij} = 1 & \forall j=1,\ldots,n, \\ & \sum_j x_{ij} = 1 & \forall i=1,\ldots,n, \\ & x_{ij}\leq A_{ij} & \forall i,j,\\ & x_{ij}\in \{0,1\} & \forall i,j. \end{array}\end{split}$

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::t M = new Model();
auto M_ = finally([&]() { M->dispose(); });

auto x = M->variable("x", Set::make(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:

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:

(11.32)$\sum_{(i,j)\in c} x_{ij} \leq |c| -1 \quad \forall c\in C.$

Thus the problem we want to solve at each step is

(11.33)$\begin{split}\begin{array}{lll} \min & \sum_{i,j} c_{ij}x_{ij} &\\ \st & \sum_i x_{ij} = 1 & \forall j=1,\ldots,n, \\ & \sum_j x_{ij} = 1 & \forall i=1,\ldots,n, \\ & x_{ij}\leq A_{ij} & \forall i,j,\\ & x_{ij}\in \{0,1\} & \forall i,j,\\ & \sum_{(i,j)\in c} x_{ij} \leq |c| -1 & \forall c\in C, \end{array}\end{split}$

where $$C$$ is the set of cycles in all the cycle covers we have seen so far. The overall solution scheme is the following:

1. set $$C$$ as the empty set,

2. solve problem (11.33),

3. if $$x$$ has only one cycle stop,

4. 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):

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

(11.34)$x_{ii} = 0\quad \forall i=1,\ldots,n.$

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.));


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:

(11.35)$x_{ij} + x_{ji} \leq 1 \quad \forall i,j=1,\ldots,n.$

Constraints (11.35) are easily implemented as follows:

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


The complete working example

Listing 11.16 The complete code for the TSP examples. Click here to download.
void tsp(int n, Matrix::t A, Matrix::t 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("x", Set::make(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)

int it = 0;
while (true)
{
M->solve();
it++;

typedef std::vector< std::tuple<int, int> > cycle_t;
std::list< cycle_t > cycles;

auto xlevel = x->level();

for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
{

if ( (*xlevel)[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 {
auto xlevel = x->level();
std::cout << "Solution\n";
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
std::cout << (int) (*xlevel)(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;
}