KeiruaProd

Traveling salesman problem with 2-opt

I’ve come across the Traveling Salesman Problem (TSP) again by chance recently. Here are some notes of an evening down this rabbit hole − elegant or simple algorithms to seemingly simple problems have a beauty to which I cannot resist :

[The TSP] asks the following question: “Given a list of cities and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?”

The fun fact here is that, despite this very simple description and the fact that us, human, can visually find a good approximation easily, it’s very hard for a computer to find the optimal solution. NP-hard, actually: there is no way to find a solution in polynomial time. And no matter how much you dig, there is no way out of this.

Despite this sad acknowledgment, it turns that the 2-opt algorithm is both a simple and fast heuristic to approximate the best tour. The idea of 2-opt is simple: starting from a random tour, exchange 2 edges so that, when they are swapped, they produce a shorter tour. Repeat until you’ve reached the necessary precision.

(image from https://arthur.maheo.net/python-local-tsp-heuristics/)

There are some variants, like 3-opt (exchange 3 edges instead of 2) all the way to k-opt, and a very good heuristic is Lin-Kerninghan algorithm, which is an adaptative algorithm that finds which k and which k-opt to use in order to best improve the current tour.

The TSP is well researched; here are some easy to read articles on the topic (research papers are linked inside, but they often are much less accessible) :

I’ve played a bit with a sample implementation of 2-opt:

# tsp.py
import numpy as np

# Calculate the euclidian distance in n-space of the route r traversing cities c, ending at the path start.
path_distance = lambda r,c: np.sum([np.linalg.norm(c[r[p]]-c[r[p-1]]) for p in range(len(r))])
# Reverse the order of all elements from element i to element k in array r.
two_opt_swap = lambda r,i,k: np.concatenate((r[0:i],r[k:-len(r)+i-1:-1],r[k+1:len(r)]))

def two_opt(cities,improvement_threshold): # 2-opt Algorithm adapted from https://en.wikipedia.org/wiki/2-opt
    route = np.arange(cities.shape[0]) # Make an array of row numbers corresponding to cities.
    improvement_factor = 1 # Initialize the improvement factor.
    best_distance = path_distance(route,cities) # Calculate the distance of the initial path.
    while improvement_factor > improvement_threshold: # If the route is still improving, keep going!
        distance_to_beat = best_distance # Record the distance at the beginning of the loop.
        for swap_first in range(1,len(route)-2): # From each city except the first and last,
            for swap_last in range(swap_first+1,len(route)): # to each of the cities following,
                new_route = two_opt_swap(route,swap_first,swap_last) # try reversing the order of these cities
                new_distance = path_distance(new_route,cities) # and check the total distance with this modification.
                if new_distance < best_distance: # If the path distance is an improvement,
                    route = new_route # make this the accepted best route
                    best_distance = new_distance # and update the distance corresponding to this route.
        improvement_factor = 1 - best_distance/distance_to_beat # Calculate how much the route has improved.
    return route # When the route is no longer improving substantially, stop searching and return the route.

From this, we can plot a route:

import matplotlib.pyplot as plt
import numpy as np
from tsp import two_opt, path_distance

# 2-opt algorithm implementation
# https://stackoverflow.com/questions/25585401/travelling-salesman-in-scipy

if __name__ == '__main__':
    # Create a matrix of cities, with each row being a location in 2-space (function works in n-dimensions).
    cities = np.random.RandomState(42).rand(70,2)
    # Find a good route with 2-opt ("route" gives the order in which to travel to each city by row number.)
    route = two_opt(cities,0.001)

    # Reorder the cities matrix by route order in a new matrix for plotting.
    new_cities_order = np.concatenate((np.array([cities[route[i]] for i in range(len(route))]),np.array([cities[0]])))
    # Plot the cities.
    plt.scatter(cities[:,0],cities[:,1])
    # Plot the path.
    plt.plot(new_cities_order[:,0],new_cities_order[:,1])
    plt.show()
    # Print the route as row numbers and the total distance travelled by the path.
    print("Route: " + str(route) + "\n\nDistance: " + str(path_distance(route,cities)))

or do some low-tech profiling :

import matplotlib.pyplot as plt
import numpy as np
from tsp import two_opt, path_distance

import time

nb_iterations = 10

if __name__ == '__main__':
    for nb_cities in range(70, 100):
        avg = 0
        for i in range(nb_iterations):
            cities = np.random.rand(nb_cities,2)
            start_time = time.$()
            route = two_opt(cities,0.001)
            elapsed = time.perf_counter_ns() - start_time
            avg += elapsed

        avg = avg / nb_iterations
        # average time in milliseconds
        print(f"{nb_cities},{avg/10**6}")

[time.perf_counter](https://docs.python.org/3/library/time.html#time.perf_counter) and perf_counter_ns are well.

We can see that the execution time follows an exponential (actually, it’s O(n^2))

See a typo ? You can suggest a modification on Github.