The Traveling Salesman Problem

Exact Formulations for Symmetric and Asymmetric TSP

Connor Owen
cowen33@gatech.edu

In the previous post, we introduced the traveling salesman problem, and then explored some heuristics that give us a solution for the problem. While these heuristics are generally computationally easy and fast, we have no guarantee on how good the solution will be, and no guarantee on optimality.

In this notebook, we will look at an exact formulation of the TSP using integer programming. This will allow us to obtain the optimal solution (guaranteed!). There are two versions of the TSP we will solve: the symmetric and the asymmetric TSP. The symmetric version is one in which the travel costs across routes are the same for both directions, i.e., the cost from A to B is the same as the cost from B to A. The asymmetric version is a TSP in which the cost from one node to the next need not be the same value in reverse.

                             *                *                 *

Let's start by looking at a small example of how to solve each version of the TSP given a list of cities and a distance matrix. Here, our cities are named for letters, and the 'distances' matrix shows the distance taking the (row) → (column) edge. You will notice that this matrix is symmetric. For instance, A → E is 25, and E → A is also 25.

In [134]:
distDict = {}
cityList = ['A', 'B', 'C', 'D', 'E', 'F']

##symmetric distance matrix
distances = [
    [0,  35, 20, 21, 25, 16],
    [35, 0,  14, 22, 9,  20],
    [20, 14, 0,  16, 6,  7],
    [21, 22, 16, 0,  21, 21],
    [25, 9,  6,  21, 0,  10],
    [16, 20, 7,  21, 10, 0]
    ]

For ease of use, we will combine these two pieces of information into a dictionary structure. This will be useful to use later when we create the optimization model.

In [135]:
##Create a dictionary to map cities to distances
for i, city in enumerate(cityList):
    distDict[city] = {x:distances[i][j] for j,x in enumerate(cityList) if x != city}

For some basic visualizations in this notebook, we will use the Networkx library from Python. This library has great graph functions, and a basic drawing function. Some of the data structures need to be changed for using this library; however, I will not cover in depth this process, but it can be figured out from the code and the Networkx documentation.

In [87]:
%matplotlib inline
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt

#create empty graph object
G = nx.Graph()

##change distDict structure for networkx graph drawing feature
graphList = [
    (city1, city2, {'weight': distDict[city1][city2]}) 
    for city1 in cityList for city2 in cityList if city1 != city2]
G.add_edges_from(graphList)

#draw graph in a circular layout
nx.draw_circular(G, with_labels=True, node_color='lightblue', node_size=1750,width=2.0)
labels = nx.get_edge_attributes(G,'weight')
#draw on edge labels; create variable so return doesn't print to notebook
draw = nx.draw_networkx_edge_labels(G, pos=nx.circular_layout(G), edge_labels=labels)

Now, we will begin the formulation of the symmetric TSP. In this notebook, we will be using the Python library "PuLP". PuLP is an optimization modeler, and can be used to create LPs and MIPs. While PuLP can connect to several different solvers (including commercial solvers), we will use the default Cbc solver.

We begin by importing all from PuLP, and then creating a problem object variable. Later, we will add constraints and an objective function to this variable.

In [60]:
from pulp import *

#create problem object 
model = pulp.LpProblem("Symmetric TSP", pulp.LpMinimize)

Next we will create all of our decision variables. In our model, we will only have binary variables, which correspond to every possible arc in our graph. A 1 for this variable will mean that we use the arc, while a 0 means that the arc is not used. There exist other formulations for this problem, but in my opinion, this one is easy to conceptualize, and thus is suitable for exploring a MIP formulation for the TSP.

In [61]:
##create arcs list
##arcs correspond to every edge without being mirrored
arcs = list(itertools.combinations(cityList, 2))

##create binary variables corresponding to every edge in the graph
arc_vars = LpVariable.dicts("Arcs", [edge for edge in arcs], 0, 1, LpInteger)

Now we will add the objective function to our problem variable. The objective of the TSP is to minimize the total distance traveled, so we will multiply the weights of the edges times the arc variables.

In [62]:
model += sum(arc_vars[edge]*distDict[edge[0]][edge[1]] for edge in arcs)

Next up we can begin adding constraints. Apart from our binary constraints on the arc variables, we have 2 conditions that must be met in the graph in order to make a valid TSP:

  • Each node has exactly two adjacent edges: one edge coming in, and one edge leaving.

  • There can be no subtours. This means, that there is no closed loop while traveling (except for the final tour).

In [63]:
import itertools

##every city has exactly two connected edges
for city in cityList:  
    ##get all of the edges that contain the city
    edge_combinations = []
    for edge in arcs:
        if city in edge:
            edge_combinations.append(edge)
    model += sum(arc_vars[edge] for edge in edge_combinations) == 2

##create all possible subtours in a list
subtours = []
for L in range(1, len(cityList)):
    for subset in itertools.combinations(cityList, L):
        subtours.append(subset)

##we can ensure no subtour exists by making sure that the number of edges
##leaving from a set of nodes is greater than the number of nodes
##subtour elimination
for subtour in subtours: 
    subtour_edges = list(itertools.combinations(subtour,2))
    model += sum(arc_vars[edge] for edge in subtour_edges) <= len(subtour) - 1

And that's all there is to the formulation! Now, we tell the model to solve itself, and then we can access the information from the optimal solution.

In [76]:
model.solve()

print(model.objective.value())
solList = []
for edge in arcs:
    if arc_vars[edge].value() == 1:
        solList.append(edge)
        print(edge)

81.0
('A', 'D')
('A', 'F')
('B', 'D')
('B', 'E')
('C', 'E')
('C', 'F')

Now that we have the optimal solution, we can once again use Networkx to visualize our solution. Again, we will not cover here the details behind implenting the visualization, but all is available in the Networkx documentation. Note that in a properly scaled TSP, arcs do not cross; however, for our network visualization, it is not to scale, so some arcs do cross.

In [109]:
colorList = []
# for edge in arcs:
#     if arc_vars[edge].value() == 1:
#         colorList.append(1)
#     else:
#         colorList.append(0)
for pair in G.edges():
    if (pair[0], pair[1]) in solList:
        #append -1 for 'connected'
        #better for color mapping
        colorList.append(-1)
    else:
        colorList.append(0)
        
nx.draw_circular(G, with_labels=True, node_color='lightblue',
            node_size=1500, edge_color=colorList, 
            width=10.0, edge_cmap=plt.cm.Pastel1)
draw = nx.draw_networkx_edge_labels(G, pos=nx.circular_layout(G), edge_labels=labels)

Now we will cover the asymmetric TSP. The asymmetric version differs from the symmetric TSP in that the distance from node A to node B is not necessarily the same as the distance from node B to node A. Take a look at the asymmetric distance matrix (corresponding to the same city list as before) to understand:

In [144]:
asymmetricDistances = [
    [0, 35, 20, 21, 25, 16],
    [27, 0, 14, 22, 9, 20],
    [15, 17, 0, 16, 6, 7],
    [28, 22, 15, 0, 21, 21],
    [22, 14, 8, 21, 0, 12],
    [11, 23, 7, 22, 10, 0]
]

As it turns out, we can actually make a transformation to this distance matrix, and preserve the symmetric problem formualtion. The conversion duplicates each node, calling it a "ghost node", and then creates a no-cost (or negative cost) weight to its twin original node. The edges from a ghost node to an original non-twin node should not exist, or should be a very high weight (so as not to be selected). This formulation can be seen here. The function below takes an asymmetric distance matrix and then converts it to a symmetric distance matrix, which can then be used with the previous formulation.

In [145]:
import copy

def asymmetricConversion(cities, distanceMatrix):
    
    matrix = copy.deepcopy(distanceMatrix)
    cityList = cities[:]
    ##create ghost nodes
    cityMirror = ['{}1'.format(city) for city in cityList]
    
    distances = []
    for i, city in enumerate(cityList):
        row = [1000 for i  in range(len(cityList))]
        row.extend([matrix[j][i] for j in range(len(matrix))])
        distances.append(row)
    for i, city in enumerate(cityList): 
        row = matrix[i]
        row.extend([1000 for i in range(len(matrix))])
        distances.append(row)
        
    cityList.extend(cityMirror)
    return cityList, distances

Now we use the asymmetricConversion function to create the "symmetric" matrix, and then create a distance dictionary again.

In [156]:
asymCityList = copy.deepcopy(cityList)
(asymCityList, distances) = asymmetricConversion(asymCityList, asymmetricDistances)

distDict2 = {}
for i, city in enumerate(asymCityList):
    distDict2[city] = {x:distances[i][j] for j,x in enumerate(asymCityList) if x != city}

Here is a look inside the final dictionary for those curious:

In [158]:
print(distDict2)

{'C': {'B': 1000, 'E1': 8, 'A': 1000, 'A1': 20, 'D': 1000, 'B1': 14, 'F': 1000, 'F1': 7, 'D1': 15, 'E': 1000, 'C1': 0}, 'D': {'B': 1000, 'E1': 21, 'C': 1000, 'A': 1000, 'A1': 21, 'B1': 22, 'F': 1000, 'F1': 22, 'D1': 0, 'E': 1000, 'C1': 16}, 'E1': {'B': 14, 'C': 8, 'A': 22, 'A1': 1000, 'D': 21, 'B1': 1000, 'F': 12, 'F1': 1000, 'D1': 1000, 'E': 0, 'C1': 1000}, 'F1': {'B': 23, 'E1': 1000, 'C': 7, 'A': 11, 'A1': 1000, 'D': 22, 'B1': 1000, 'F': 0, 'D1': 1000, 'E': 10, 'C1': 1000}, 'E': {'B': 1000, 'E1': 0, 'C': 1000, 'A': 1000, 'A1': 25, 'D': 1000, 'B1': 9, 'F1': 10, 'D1': 21, 'F': 1000, 'C1': 6}, 'C1': {'B': 17, 'E1': 1000, 'C': 0, 'A': 15, 'A1': 1000, 'D': 16, 'B1': 1000, 'F': 7, 'F1': 1000, 'D1': 1000, 'E': 6}, 'B': {'E1': 14, 'C': 1000, 'A': 1000, 'A1': 35, 'D': 1000, 'B1': 0, 'E': 1000, 'F1': 23, 'D1': 22, 'F': 1000, 'C1': 17}, 'B1': {'B': 0, 'C': 14, 'A': 27, 'A1': 1000, 'D': 22, 'E1': 1000, 'F': 20, 'F1': 1000, 'D1': 1000, 'E': 9, 'C1': 1000}, 'A': {'B': 1000, 'E1': 22, 'C': 1000, 'A1': 0, 'D': 1000, 'B1': 27, 'E': 1000, 'F1': 11, 'D1': 28, 'F': 1000, 'C1': 15}, 'A1': {'B': 35, 'E1': 1000, 'C': 20, 'A': 0, 'D': 21, 'B1': 1000, 'F': 16, 'F1': 1000, 'D1': 1000, 'E': 25, 'C1': 1000}, 'D1': {'B': 22, 'E1': 1000, 'C': 15, 'A': 28, 'A1': 1000, 'D': 0, 'B1': 1000, 'F': 21, 'F1': 1000, 'E': 21, 'C1': 1000}, 'F': {'B': 1000, 'E1': 12, 'C': 1000, 'A': 1000, 'A1': 16, 'D': 1000, 'B1': 20, 'F1': 0, 'D1': 21, 'E': 1000, 'C1': 7}}

Below is a function form of the symmetric TSP that was broken down earlier. We will take the asymmetric TSP matrix and use it with this function.

In [224]:
from pulp import *
import itertools

def symmetric_TSP_Pulp(cityList, distDict):
    model = pulp.LpProblem("Symmetric TSP", pulp.LpMinimize)
    
    ##SET CONSTRUCTION##
    
    ##create all possible subtours
    subtours = []
    for L in range(1, len(cityList)):
        for subset in itertools.combinations(cityList, L):
            subtours.append(subset)
    
    ##arcs correspond to every edge without being mirrored
    arcs = list(itertools.combinations(cityList, 2))
    
    ##DECISION VARIABLES##
    
    ##create binary variables corresponding to every edge in the graph
    arc_vars = LpVariable.dicts("Arcs", [edge for edge in arcs], 0, 1, LpInteger)
    
    ##OBJECTIVE FUNCTION##
    model += sum(arc_vars[edge]*distDict[edge[0]][edge[1]] for edge in arcs)
    
    ##CONSTRAINTS##
    
    ##every city has exactly two connected edges
    for city in cityList:  
        ##get all of the edges that contain the city
        edge_combinations = []
        for edge in arcs:
            if city in edge:
                edge_combinations.append(edge)
        model += sum(arc_vars[edge] for edge in edge_combinations) == 2
    
    ##subtour elimination
    for subtour in subtours: 
        subtour_edges = list(itertools.combinations(subtour,2))
        model += sum(arc_vars[edge] for edge in subtour_edges) <= len(subtour) - 1
        
    model.solve()
    
    solList = []
    print(model.objective.value())
    for edge in arcs:
        if arc_vars[edge].value() == 1:
            if edge[0] == edge[1][0:1]:
                pass
            elif edge[0] != edge[1][0:1]:
                solList.append((edge[0], edge[1][0:1]))
            print(edge)  
            
    return solList

Now we run the function with the asymmetric city List and the corresponding distance dictionary.

In [212]:
solList = symmetric_TSP_Pulp(asymCityList, distDict2)

78.0
('A', 'A1')
('A', 'F1')
('B', 'B1')
('B', 'D1')
('C', 'C1')
('C', 'E1')
('D', 'A1')
('D', 'D1')
('E', 'B1')
('E', 'E1')
('F', 'C1')
('F', 'F1')

From here, we can create the visualization again using Networkx, this time without edge weights. Since Networkx automatically takes out repeated edges, we have to reverse the solution list to ensure that all selected edges are colored.

In [213]:
colorList = []
solListReverse = [(i,j) for (j, i) in solList]
for pair in G.edges():
    if (pair[0], pair[1]) in solList or (pair[0], pair[1]) in solListReverse:
        #append -1 for 'connected'
        #better for color mapping
        colorList.append(-1)
    else:
        colorList.append(0)
        
nx.draw_circular(G, with_labels=True, node_color='lightblue',
            node_size=1500, edge_color=colorList, 
            width=10.0, edge_cmap=plt.cm.Pastel1)
# draw = nx.draw_networkx_edge_labels

A better option for a visualization would be a directed graph, meaning that the edges can only be traversed in the indicated direction. Networkx has the ability to create and draw a directed graph. The below visualization shows how you can implement a directed graph from our solution information.

In [223]:
G_Dir = nx.DiGraph()
G_Dir.add_edges_from(solList)

colorList = []
solListReverse = [(i,j) for (j, i) in solList]
for pair in G_Dir.edges():
    if (pair[0], pair[1]) in solList or (pair[0], pair[1]) in solListReverse:
        #append -1 for 'connected'
        #better for color mapping than 1
        colorList.append(-1)
    else:
        colorList.append(0)

nx.draw_circular(G_Dir, with_labels=True, node_color='lightblue',
            node_size=1500, edge_color=colorList, 
            width=5.0, edge_cmap=plt.cm.Pastel1)

This is a basic formulation for the Traveling Salesman Problem. While it is guaranteed to return an optimal solution, it is still not enough with current technology to be able to solve large instances of the problem in a short amount of time. One simpler improvement to this formulation is to modify the subtour elimination constraint. With the current formulation, we have to consider every possible subtour, which means our set of subtours grows exponentially with the number of nodes in the problem. A smarter formulation generates these constraints as needed, after a subproblem is solved. This subproblem can either be a relaxation of all integer constraints, or just a relaxation of the subtour constraint.