import copy

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import numpy.random as rnd
import tsplib95
import tsplib95.distances as distances

from alns import ALNS
from alns.accept import HillClimbing
from alns.select import RouletteWheel
from alns.stop import MaxRuntime

%matplotlib inline

SEED = 7654

The travelling salesman problem

The travelling salesman problem (TSP) is a classic problem in operations research. It asks how to construct the minimum distance tour between a number of nodes, such that each node is visited once and the tour concludes at the starting city (that is, it forms a cycle). It is perhaps the best-known problem in the class of NP-hard problems.

Data

There are a considerable number of test data sets available for the TSP, varying in size from a hundred or so locations to many hundreds of thousands. For the sake of exposition, we shall use one of the smaller data sets: the data from the XQF131 VLSI instance, made available here. It consists of ‘only’ 131 nodes, with an optimal tour length of 564.

DATA = tsplib95.load('data/xqf131.tsp')
CITIES = list(DATA.node_coords.keys())

# Precompute the distance matrix - this saves a bunch of time evaluating moves.
# + 1 since the cities start from one (not zero).
COORDS = DATA.node_coords.values()
DIST = np.empty((len(COORDS) + 1, len(COORDS) + 1))

for row, coord1 in enumerate(COORDS, 1):
    for col, coord2 in enumerate(COORDS, 1):
        DIST[row, col] = distances.euclidean(coord1, coord2)

SOLUTION = tsplib95.load('data/xqf131.opt.tour')
OPTIMAL = DATA.trace_tours(SOLUTION.tours)[0]

print(f"Total optimal tour length is {OPTIMAL}.")

Total optimal tour length is 564.

def draw_graph(graph, only_nodes=False):
    """
    Helper method for drawing TSP (tour) graphs.
    """
    fig, ax = plt.subplots(figsize=(12, 6))

    if only_nodes:
        nx.draw_networkx_nodes(graph, DATA.node_coords, node_size=25, ax=ax)
    else:
        nx.draw_networkx(graph, DATA.node_coords, node_size=25, with_labels=False, ax=ax)

draw_graph(DATA.get_graph(), only_nodes=True)

To use the ALNS meta-heuristic, we need to have destroy and repair operators that work on a proposed solution, and a way to describe such a solution in the first place. Let’s start with the solution state.

Solution state

class TspState:
    """
    Solution class for the TSP problem. It has two data members, nodes, and edges.
    nodes is a list of IDs. The edges data member, then, is a mapping from each node
    to their only outgoing node.
    """

    def __init__(self, nodes, edges):
        self.nodes = nodes
        self.edges = edges

    def objective(self):
        """
        The objective function is simply the sum of all individual edge lengths,
        using the rounded Euclidean norm.
        """
        return sum(DIST[node, self.edges[node]] for node in self.nodes)

    def to_graph(self):
        """
        NetworkX helper method.
        """
        graph = nx.Graph()

        for node in self.nodes:
            graph.add_node(node, pos=DATA.node_coords[node])

        for node_from, node_to in self.edges.items():
            graph.add_edge(node_from, node_to)

        return graph

Destroy operators

Destroy operators break parts of a solution down, leaving an incomplete state. This is the first part of each iteration of the ALNS meta-heuristic; the incomplete solution is subsequently repaired by any one repair operator. We will consider three destroy operators: worst removal, path removal and random removal. We will also use a separate parameter, the degree of destruction, to control the extent of the damage done to a solution in each step.

DEGREE_OF_DESTRUCTION = 0.1

def edges_to_remove(state):
    return int(len(state.edges) * DEGREE_OF_DESTRUCTION)

def worst_removal(current, rnd_state):
    """
    Worst removal iteratively removes the 'worst' edges, that is,
    those edges that have the largest distance.
    """
    destroyed = copy.deepcopy(current)

    worst_edges = sorted(destroyed.nodes,
                         key=lambda node: DIST[node, destroyed.edges[node]])

    for idx in range(edges_to_remove(current)):
        del destroyed.edges[worst_edges[-(idx + 1)]]

    return destroyed

def path_removal(current, rnd_state):
    """
    Removes an entire consecutive sub-path, that is, a series of
    contiguous edges.
    """
    destroyed = copy.deepcopy(current)

    node_idx = rnd_state.choice(len(destroyed.nodes))
    node = destroyed.nodes[node_idx]

    for _ in range(edges_to_remove(current)):
        node = destroyed.edges.pop(node)

    return destroyed

def random_removal(current, rnd_state):
    """
    Random removal iteratively removes random edges.
    """
    destroyed = copy.deepcopy(current)

    for idx in rnd_state.choice(len(destroyed.nodes),
                                edges_to_remove(current),
                                replace=False):
        del destroyed.edges[destroyed.nodes[idx]]

    return destroyed

Repair operators

We implement a simple, greedy repair strategy. It determines a set of nodes that are currently not visited, and then links these up to the tour such that it forms one cycle.

def would_form_subcycle(from_node, to_node, state):
    """
    Ensures the proposed solution would not result in a cycle smaller
    than the entire set of nodes. Notice the offsets: we do not count
    the current node under consideration, as it cannot yet be part of
    a cycle.
    """
    for step in range(1, len(state.nodes)):
        if to_node not in state.edges:
            return False

        to_node = state.edges[to_node]

        if from_node == to_node and step != len(state.nodes) - 1:
            return True

    return False

def greedy_repair(current, rnd_state):
    """
    Greedily repairs a tour, stitching up nodes that are not departed
    with those not visited.
    """
    visited = set(current.edges.values())

    # This kind of randomness ensures we do not cycle between the same
    # destroy and repair steps every time.
    shuffled_idcs = rnd_state.permutation(len(current.nodes))
    nodes = [current.nodes[idx] for idx in shuffled_idcs]

    while len(current.edges) != len(current.nodes):
        node = next(node for node in nodes
                    if node not in current.edges)

        # Computes all nodes that have not currently been visited,
        # that is, those that this node might visit. This should
        # not result in a subcycle, as that would violate the TSP
        # constraints.
        unvisited = {other for other in current.nodes
                     if other != node
                     if other not in visited
                     if not would_form_subcycle(node, other, current)}

        # Closest visitable node.
        nearest = min(unvisited,
                      key=lambda other: DIST[node, other])

        current.edges[node] = nearest
        visited.add(nearest)

    return current

Initial solution

We start from a very simple, greedily constructed initial solution. This solution is not good (as can clearly be seen below), but it is feasible.

random_state = rnd.RandomState(SEED)
state = TspState(CITIES, {})

init_sol = greedy_repair(state, random_state)
print(f"Initial solution objective is {init_sol.objective()}.")

Initial solution objective is 870.0.

draw_graph(init_sol.to_graph())

Heuristic solution

Here we perform the ALNS procedure. The heuristic is given 60 seconds of runtime.

alns = ALNS(random_state)

alns.add_destroy_operator(random_removal)
alns.add_destroy_operator(path_removal)
alns.add_destroy_operator(worst_removal)

alns.add_repair_operator(greedy_repair)

select = RouletteWheel([3, 2, 1, 0.5], 0.8, 3, 1)
accept = HillClimbing()
stop = MaxRuntime(60)

result = alns.iterate(init_sol, select, accept, stop)

solution = result.best_state
objective = solution.objective()
pct_diff = 100 * (objective - OPTIMAL) / OPTIMAL

print(f"Best heuristic objective is {objective}.")
print(f"This is {pct_diff:.1f}% worse than the optimal solution, which is {OPTIMAL}.")

Best heuristic objective is 576.0.
This is 2.1% worse than the optimal solution, which is 564.

_, ax = plt.subplots(figsize=(12, 6))
result.plot_objectives(ax=ax, lw=2)

Let’s have a look at the solution:

draw_graph(solution.to_graph())

Conclusions

In the code above we implemented a very simple heuristic for the TSP, using the ALNS meta-heuristic framework. We did not tinker too much with the various hyperparameters available on the ALNS implementation, but even for these relatively basic heuristic methods and workflow we find a very good result - just 2% worse than the optimal tour.

This notebook showcases how the ALNS library may be put to use to construct powerful, efficient heuristic pipelines from simple, locally greedy operators.