Solving TSPs to Optimality with Integer Linear Programming in Python

Posted on 28. August 202428. August 2024 by Fabian Leuthold

In this post we will learn how to optimally solve the TSP problem using integer linear programming using Google OR-Tools for mathmatical modelling in Python. In previous posts I have already presented two ways of solving the TSP using heuristic approaches: Construct solutions using the Nearest Neighbor construction heuristic and improve solutions using the 2-opt based local search. In this post we will see a solution approach based on mathematical programming to always find the optimal solution.

The Problem

The Traveling Salesman Problem (TSP) is defined as follows: Given a set of locations specified by their geographical coordinates, the goal is to determine the shortest possible route that visits each location exactly once and returns to the starting point. The order of visiting the locations is flexible, and the objective is to minimize the total travel distance.

Example for a TSP trip (by Xypron, public domain)

In practice, we usually start with a set of addresses for all locations. To solve the problem, we need a distance matrix that provides the distance in kilometers between every pair of locations $(i,j)$ . One approach is to use an OSRM server with OpenStreetMap data for the desired country. First, geocode the addresses to obtain latitude and longitude coordinates, and then generate the distance matrix $d_{ij}$ .

Since generating distance matrices is not our primary focus, we will work with a TSP problem provided in the standardized TSP file format. We will use a simple command to convert this file into a spatial distance matrix.

The Model

Let’s start to model the problem as an Mixed Integer Linear Program (MILP). When designing an MILP, we usually follow the following recipe:

Choose the decision variables
Choose the auxiliary variables (if required)
Define the indices and sets (if required)
Define the parameters
Define the constraints
Define the objective function
Validate the model

1) Choose the Decision Variables

To model the sequence of visits among the points, we typically use a decision variable with two indices, $i$ and $j$ . This approach helps in representing the relationship between any pair of locations, which is useful for tracking the sequence of the tour. Remember, it’s important to model the connection between each location and its previous or next location in the TSP tour.

Therefore we define the decision variable as follows: $x_{ij}$ = 1, if location $j$ is visited in our TSP trip just after location $i$ , 0 otherwise.

2) Choose the Auxiliary Variables

At the moment, we don’t yet have a need for auxiliary variables.

3) Define the Indices and Sets

We have two indices, $i$ and $j$ , both are coming from the same set $I={1, 2, 3, ..., n}$ with $n$ being the number of locations to be visited in our TSP.

4) Define the Parameters

The only parameter we need here is the distance matrix $d_{ij}$ from above, providing the distance in km between location $i$ and $j$ . Since our problem is symmetrical, the vaues in the distance matrix are symmetrical with respect to the diagonal (lower / upper triangular matrix are identical).

5) Define the Constraints

Leave every point for exactly one successor:
$\forall i \in I: \sum_j x_{ij} = 1$
Reach every point from exactly one predecessor:
$\forall j \in I: \sum_i x_{ij} = 1$

At first sight, these constraints might seem to be complete. For now, let’s leave it that way. We’ll comeback to this later, when following the 7-step recipe presented above.

6) Objective Function

Since we want to minimize the total distance, the objective function is straight forward:

$min \sum_{(i,j) \in (I \times I)} x_{ij} d_{ij}$

We just have to minimize the sum of all sections in the trip.

7) Validate the Model

To validate the model, we should first check it at a logical level. This involves verifying that the model handles special cases and boundary conditions correctly.

In particular, the model must ensure that each location is visited exactly once and returned to, resulting in a closed tour. However, the current model may allow multiple closed tours, which is not permitted. This issue was discussed in more detail in a previous post. We need to address this to complete the model and eliminate any impermissible configurations.

Our model has failed the logical test, so we need to enhance it to address this issue. To correct this, we will revise the model by introducing the missing subtour elimination constraints (For further details, please refer to the following link).

2) Choose the auxiliary Variables

The Miller-Tucker-Zemlin subtour elimination constraint requires an additional decision variable $u_i$ to enumerate the stops in sequential order.

5) Define the Constraints

Now based on the new auxiliary variable $u_i$ we can formulate the subtour elimination constraints. Since it doesn’t matter, where to start, let’s start at the first location:

$u_1 = 1$

$\forall i \ne 1: \;\; 2\le u_i \le n$

$\forall i \ne 1, \forall j \ne 1: \;\; u_i - u_j +1 \le (n-1)(1-x_{ij})$

Since we have n locations, all subsequent $u_i$ must have assigned a vaule between 2 and $n$ . The last constraint makes sure, that the index of every succeeding stop gets assigned an index which is at least 1 larger than the predecessor’s index. Together with the above condition, each successor can only be assigned an index that is exactly one greater than that of its predecessor.

7) Validate the Model

Now let’s validate the model once more. From the logical standpoint of view, the model seems to be correct and complete now. Now it’s time to check wheter the model solves the problem instances as expected in practice.

Hint: The TSP is one of the hardest combinatorial problems and the model presented here allows to find optimal solutions. However, it’s by far not the best model, since the integer relaxation of the MTZ-constraint is very weak. Even one of the smallest national problem instances provided in the university of waterloo, the Djibouti – instance, makes gurobi – the world’s best MIP solver – sweat!

The Code

Now let’s implement the solution in python. The cool thing is, that using google or-tools we have a powerful tool at hand to translate the model into code. Let’s first see, how to do this:

    # constraint 1: leave every point exactly once
    log.info('Creating ' + str(num_nodes) + ' Constraint 1... ')
    for i in all_nodes:
        model.Add(sum(x[(i, j)] for j in all_nodes) == 1)

As we can see, that the constraint can easily be translated. The for-all clause ends up in the loop, the sigma-formula is translated directly with binding the other left indices using for … in … clauses and the term can be take exactly from the formal constraint definition.

Let’s translate the other constraints:

    # constraint 2: reach every point from exactly one other point
    log.info('Creating ' + str(num_nodes) + ' Constraint 2... ')
    for j in all_nodes:
        model.Add(sum(x[(i, j)] for i in all_nodes) == 1)

    # constraint 3.1: subtour elimination constraints (Miller-Tucker-Zemlin) part 1
    log.info('Creating 1 Constraint 3.1... ')
    model.Add(u[0] == 1)

    # constraint 3.2: subtour elimination constraints (Miller-Tucker-Zemlin) part 2
    log.info('Creating ' + str(len(all_but_first_nodes)) + ' Constraint 3.2... ')
    for i in all_but_first_nodes:
        model.Add(2 <= u[i])
        model.Add(u[i] <= num_nodes)

    # constraint 3.3: subtour elimination constraints (Miller-Tucker-Zemlin) part 3
    log.info('Creating ' + str(len(all_but_first_nodes)) + ' Constraint 3.2... ')
    for i in all_but_first_nodes:
        for j in all_but_first_nodes:
            model.Add(u[i] - u[j] + 1 <= (num_nodes - 1) * (1 - x[(i, j)]))

Finally we have to translate the objective function:

    # Minimize the total distanceg der Anzahl geänderte Dienste und Anzahl Dienstfahrten
    model.Minimize(sum(x[(i, j)] * dima[(i, j)] for i in all_nodes for j in all_nodes))

That’s it! The rest of the code is about loading the TSP data, generating the distance matrix, calling the solver – i.e. gurobi – to solve the model and then output the solution:

import math
from builtins import min

import scipy.spatial as sp
import pandas as pd
import numpy as np
import logging as log
import matplotlib.pyplot as plt
import json
import sys
# also required packages: openpyxl
from ortools.linear_solver import pywraplp

def solve_OrTools(dima):
    '''
    generate mip model using google or-tools and solve it

    :param dima: the distance matrix
    :return:  solution X, model, status
    '''

    if dima.ndim != 2 or dima.shape[0] != dima.shape[1]:
        raise ValueError("Invalid dima dimensions detected. Square matrix expected.")

    # determine number of nodes
    num_nodes = dima.shape[0]
    all_nodes = range(0, num_nodes)
    all_but_first_nodes = range(1, num_nodes)

    # Create the model.
    solver_name = 'CP-SAT'
    log.info('Instantiating solver ' + solver_name)
    model = pywraplp.Solver.CreateSolver(solver_name)
    model.EnableOutput()

    log.info('Defining MIP model... ')
    # generating decision variables X_ij
    log.info('Creating ' + str(num_nodes * num_nodes) + ' boolean x_ij variables... ')
    x = {}
    for i in all_nodes:
        for j in all_nodes:
            x[(i, j)] = model.BoolVar('x_i%ij%i' % (i, j))

    log.info('Creating ' + str(num_nodes) + ' boolean u_i variables... ')
    u = {}
    for i in all_nodes:
        u[i] = model.IntVar(0, num_nodes, 'u_i%i' % i)

    # constraint 1: leave every point exactly once
    log.info('Creating ' + str(num_nodes) + ' Constraint 1... ')
    for i in all_nodes:
        model.Add(sum(x[(i, j)] for j in all_nodes) == 1)

    # constraint 2: reach every point from exactly one other point
    log.info('Creating ' + str(num_nodes) + ' Constraint 2... ')
    for j in all_nodes:
        model.Add(sum(x[(i, j)] for i in all_nodes) == 1)

    # constraint 3.1: subtour elimination constraints (Miller-Tucker-Zemlin) part 1
    log.info('Creating 1 Constraint 3.1... ')
    model.Add(u[0] == 1)

    # constraint 3.2: subtour elimination constraints (Miller-Tucker-Zemlin) part 2
    log.info('Creating ' + str(len(all_but_first_nodes)) + ' Constraint 3.2... ')
    for i in all_but_first_nodes:
        model.Add(2 <= u[i])
        model.Add(u[i] <= num_nodes)

    # constraint 3.3: subtour elimination constraints (Miller-Tucker-Zemlin) part 3
    log.info('Creating ' + str(len(all_but_first_nodes)) + ' Constraint 3.2... ')
    for i in all_but_first_nodes:
        for j in all_but_first_nodes:
            model.Add(u[i] - u[j] + 1 <= (num_nodes - 1) * (1 - x[(i, j)]))

    # Minimize the total distance
    model.Minimize(sum(x[(i, j)] * dima[(i, j)] for i in all_nodes for j in all_nodes))

    log.info('Solving MIP model... ')
    status = model.Solve()

    return x, model, status

def print_solution(x):
    num_nodes = len(x)
    all_nodes = range(0, int(math.sqrt(num_nodes)))
    for i in all_nodes:
        for j in all_nodes:
            if int(x[(i,j)].solution_value()) == 1:
                log.info(f'x({i}/{j})=' + str(int(x[i,j].solution_value())))


def extract_tsp_sequence(x, num_nodes):
    """Extrahiert die Reihenfolge der Knoten aus den Entscheidungsvariablen x[i,j].

    Args:
        x: Ein Dictionary mit den 2D Entscheidungsvariablen x[i,j].
        num_nodes: Die Anzahl der Knoten im TSP.

    Returns:
        Eine Liste, die die Reihenfolge der besuchten Knoten enthält.
    """
    # Umgekehrte Zuordnung für den nächsten Knoten: von Knoten zu Knoten
    next_node = {}

    # Fülle die next_node Map basierend auf den Entscheidungsvariablen
    for i in range(num_nodes):
        for j in range(num_nodes):
            if x[i, j].solution_value() == 1:
                next_node[i] = j

    # Finde die Reihenfolge der besuchten Knoten
    sequence = []
    current_node = 0  # Wir starten mit einem Anfangsknoten, z.B. 0
    for _ in range(num_nodes):
        sequence.append(current_node)
        current_node = next_node[current_node]

    return sequence

def plot_tsp_solution(tsp, solution, title):
    # Extract coordinates
    x_coords = [tsp.lat[node] for node in solution]
    y_coords = [tsp.lng[node] for node in solution]

    # Add the start point to the end to complete the tour
    x_coords.append(x_coords[0])
    y_coords.append(y_coords[0])

    # Plot the nodes
    plt.scatter(x_coords, y_coords, c='gray')

    # Plot the path
    plt.plot(x_coords, y_coords, c='green')

    # Annotate the nodes
    if len(tsp) < 200:
        for idx, node in enumerate(solution):
            plt.annotate(str(node), (x_coords[idx], y_coords[idx]))

    plt.xlabel('X Coordinate')
    plt.ylabel('Y Coordinate')
    plt.title('TSP Visualization (' + title + ')')
    plt.show()


def main():
    # configure logger for info level
    log.basicConfig(
        format='%(asctime)s %(levelname)-8s %(message)s',
        level=log.INFO,
        datefmt='%Y-%m-%d %H:%M:%S',
        stream=sys.stdout)

    # load tsp instance
    tsp_problem = 'berlin52.tsp'

    log.info("Reading TSP problem Instance " + tsp_problem)
    tsp = pd.read_csv('./data/' + tsp_problem, sep=' ', skiprows=6, dtype = float,
                      names=['nodeId', 'lat', 'lng'], skipfooter=1, engine='python')
    tsp = tsp.sort_values(by='nodeId', inplace=False)

    A = tsp[['lat', 'lng']].to_numpy()
    dima = sp.distance_matrix(A, A)

    # now solve problem
    x, model, status = solve_OrTools(dima)

    # check problem response
    if status == pywraplp.Solver.OPTIMAL:
        log.info('Solution:')
        log.info('optimal solution found.')
        log.info('Objective value =' + str(model.Objective().Value()))
        tour = extract_tsp_sequence(x, len(A))
        print(f"Optimal Tour: {tour}")
        plot_tsp_solution(tsp, tour, "berlin52")
    elif status == pywraplp.Solver.INFEASIBLE:
        log.info('The problem is infeasible.')
    else:
        log.info('The problem could not be solved. Return state was: ' + status)

# Press the green button in the gutter to run the script.
if __name__ == '__main__':
    main()

The complete source can be downloaded from my github account here. To solve the berlin52 TSP problem instance, I used the free and excellent Google constraint programming solver, called CP-SAT. You could as well use ‘SCIP’ or ‘GUROBI’ – for larger instances – instead. There are also plenty of other solvers you can use here, some of which are provided by google’s or-tools, others are provided from other open source projects and published under vairous licenses.

The following plot shows the optimal tour resulting from the MILP presented above for the berlin52 TSP instance:

The optimal tour distance turns out to be 7544.365901904088 units.