TSP / ATSP algorithm.

This post will be dedicated to the Travelling salesman problem (TSP), one of the most known combinatorial optimization problem. The goal of this problem is simple: A salesman has to visit a number of places (nodes). Which route would visits them all in the shortest total distance (time, cost, etc...)
The simplest way to solve the problem is list and try all possible routes and them select the shortest one. The problem here is the computational time which will grow faster than exponential (Stirling's formula: $n! \approx \sqrt{2\pi n}(\frac{n}{e})^n$)

There are other solutions which can reduce this computational time, for example dynamic programming which solves the problem in exponential time E. The first implemented method using that more efficient method was the Held-Karp algorithm which finds the shortest route in $2^n n^2$
Other methods are: k-opt, Greedy algorithm, Greedy k-opt, Genetic algorithm, Simulating annealing, neuronal network approach among others.

On this post I will show an implementation for the simplest and naive method, an implementation for dealing with asymmetric input matrices and finally we will compare these methods with the genetic algorithm(this implementation has too many lines of code for one post)

import itertools
import numpy as np

INF = 1000000


class TSPExact:
    def __init__(self, matrix, redundant=False):
        """
        input: Adjacency matrix N x N
        len(nodes) = shape(input)
        @param matrix: N x N matrix
        @param redundant: false by default (n-1)!
        """
        # 1: Check if the adjacency matrix is symmetric
        self.cost_matrix = np.array(matrix)
        self.symmetry = None
        if self.symmetric() is False:
            self.symmetry = False
            self.as_symmetric()

        self.size = len(self.cost_matrix)
        self.nodes = range(self.size)
        self.all_tours = itertools.permutations

        if redundant is False:
            self.tours = self.non_redundant()
        else:
            self.tours = list(self.all_tours(self.nodes))

    def as_symmetric(self):
        """
        Reformulate an asymmetric TSP as a symmetric TSP: 
        "Jonker and Volgenant 1983"
        This is possible by doubling the number of nodes. For each city a dummy 
        node is added: (a, b, c) => (a, a', b, b', c, c')

        distance = "value"
        distance (for each pair of dummy nodes and pair of nodes is INF)
        distance (for each pair node and its dummy node is -INF)
        ------------------------------------------------------------------------
          |A'   |B'   |C'   |A    |B    |C    |
        A'|0    |INF  |INF  |-INF |dBA  |dCA  |
        B'|INF  |0    |INF  |dAB  |-INF |dCB  |
        C'|INF  |INF  |0    |dAC  |dBC  |-INF |
        A |-INF |dAB  |dAC  |0    |INF  |INF  |
        B |dBA  |-INF |dBC  |INF  |0    |INF  |
        C |dCA  |dCB  |-INF |INF  |INF  |0    |

        For large matrix an exact solution is infeasible
        if N > 5 (N = N*2 > 10) then use other techniques: 
        Heuristics and relaxation methods
        @return: new symmetric matrix

        [INF][A.T]
        [A  ][INF]
        """

        shape = len(self.cost_matrix)
        mat = np.identity(shape) * - INF + self.cost_matrix

        new_shape = shape * 2
        new_matrix = np.ones((new_shape, new_shape)) * INF
        np.fill_diagonal(new_matrix, 0)

        # insert new matrices
        new_matrix[shape:new_shape, :shape] = mat
        new_matrix[:shape, shape:new_shape] = mat.T
        # new cost matrix after transformation
        self.cost_matrix = new_matrix

    def total_tour(self, tour):
        total = sum(self.cost_matrix[tour[node], tour[(node + 1) % self.size]] for node in self.nodes)
        return total, tour

    def shortest(self):
        min_tour = min(self.total_tour(tour) for tour in self.tours)
        if self.symmetry is False:
            min_tour_real = min_tour[0] + self.size / 2 * INF
            min_cycle_real = min_tour[1][0::2]
            return min_tour_real, min_cycle_real
        else:
            return min_tour

    def non_redundant(self):
        start_node = self.nodes[0]
        new_set = list(self.nodes)[1:]
        return [list(np.append(start_node, tour)) for tour in self.all_tours(new_set)]

    def symmetric(self):
        #@return: boolean symmetric | asymmetric
        mat = [tuple(row) for row in self.cost_matrix]
        return mat == zip(*mat)

    @property
    def tour_iterations(self):
        return list(self.tours)

    @property
    def tour_size(self):
        return self.size

The code above tried to be self explanatory, but it might be some points unclear. The trickiest part is the transformation from asymmetric matrix to symmetric matrix. As shown in the code, we are using the implementation suggested by Jonker R, Volgenant T (1983) (“Transforming asymmetric into symmetric traveling salesman problems.” Operations Research Letters, 2, 161–163). The process behind is simple, doubling the number of nodes and transposing the original matrix it is possible get a symmetric matrix which can be solvable applying the current known methods. Finally this algorithm requires the removal of all the dummy nodes and a subtle operation in order to get the real shortest tour value.

$shortestRoute (real) = shortestRoute + \frac{Size(InputMatrix)}{2} \times INF.Value $
We can use slicing in Numpy for removing dummy nodes.

min_tour_real = min_tour[0] + self.size / 2 * INF
min_cycle_real = min_tour[1][0::2]

Unit testing

Now it's time to compare some algorithms. I have hard coded three different matrices for testing purpose. We must keep in mind that the matrix size will be multiplied by two if input matrix is asymmetric. In addition, one must point out that the basic algorithm becomes infeasible when the number of nodes is higher than 12 (6 - asymmetric matrix)

from Algorithms.scripts.TSP_Algorithm import *
import unittest
import time

MATRIX = [[0, 750, 43, 402], [750, 0, 736, 671], [43, 736, 0, 117], [402, 671, 117, 0]]
MATRIX5 = [[0, 750, 43, 402, 133], [750, 0, 736, 671, 85], [43, 736, 0, 117, 934], 
     [402, 671, 117, 0, 982], [133, 85, 934, 982, 0]]

MATRIX8 = [[0, 750, 431, 402, 133, 489, 310, 566], [750, 0, 736, 671, 85, 174, 870, 927],
           [431, 736, 0, 117, 934, 65, 939, 305], [402, 671, 117, 0, 982, 727, 911, 870],
           [133, 85, 934, 982, 0, 956, 834, 118], [489, 174, 65, 727, 956, 0, 322, 397],
           [310, 870, 939, 911, 834, 322, 0, 127], [566, 927, 305, 870, 118, 397, 127, 0]]


class TestAlgorithms(unittest.TestCase):
    @staticmethod
    def test_exact():
        start_time = time.time()
        tsp = TSPExact(MATRIX5)
        total_value, tour = tsp.shortest()
        print "__exact__"
        print total_value
        print tour
        print "Time (s)", time.time() - start_time
  
    @staticmethod
    def test_genetic():
        start_time = time.time()
        tsp = TSPGenetic(MATRIX5)
        total, route = tsp.calculate()
        print "__genetic__"
        print total
        print route
        print "Time (s)", time.time() - start_time

Some results

MATRIX(Asymmetric)
__exact__
1580.0
[0, 2, 3, 1]
Time (s) 0.130000114441
__genetic__
1580
[0, 2, 3, 1]
Time (s) 0.00100016593933

MATRIX5(Asymmetric)
__exact__
1049.0
[0, 2, 3, 1, 4]
Time (s) 10.1160001755
__genetic__
1049
[0, 2, 3, 1, 4]
Time (s) 0.000999927520752

MATRIX8(Asymmetric)
__exact__
16 nodes - time...
__genetic__
1398
[0, 3, 2, 5, 1, 4, 7, 6]
Time (s) 0.00399994850159


Comments

As shown on the previous tests, the problem becomes intractable for the Exact algorithm as the number of nodes increases. Using the Genetic algorithm we got the same results, however this algorithm might not obtain the exact solution in large graphs.