Hide navigation sidebar
Hide table of contents sidebar
Back to top
Toggle table of contents sidebar

# tsp.py

    #!/usr/bin/env python3

    # Copyright 2025, Gurobi Optimization, LLC

    # Solve a traveling salesman problem on a randomly generated set of points
    # using lazy constraints.  The base MIP model only includes 'degree-2'
    # constraints, requiring each node to have exactly two incident edges.
    # Solutions to this model may contain subtours - tours that don't visit every
    # city.  The lazy constraint callback adds new constraints to cut them off.

    import sys
    import logging
    import math
    import random
    from collections import defaultdict
    from itertools import combinations

    import gurobipy as gp
    from gurobipy import GRB

    def shortest_subtour(edges):
        """Given a list of edges, return the shortest subtour (as a list of nodes)
        found by following those edges. It is assumed there is exactly one 'in'
        edge and one 'out' edge for every node represented in the edge list."""

        # Create a mapping from each node to its neighbours
        node_neighbors = defaultdict(list)
        for i, j in edges:
            node_neighbors[i].append(j)
        assert all(len(neighbors) == 2 for neighbors in node_neighbors.values())

        # Follow edges to find cycles. Each time a new cycle is found, keep track
        # of the shortest cycle found so far and restart from an unvisited node.
        unvisited = set(node_neighbors)
        shortest = None
        while unvisited:
            cycle = []
            neighbors = list(unvisited)
            while neighbors:
                current = neighbors.pop()
                cycle.append(current)
                unvisited.remove(current)
                neighbors = [j for j in node_neighbors[current] if j in unvisited]
            if shortest is None or len(cycle) < len(shortest):
                shortest = cycle

        assert shortest is not None
        return shortest

    class TSPCallback:
        """Callback class implementing lazy constraints for the TSP.  At MIPSOL
        callbacks, solutions are checked for subtours and subtour elimination
        constraints are added if needed."""

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

        def __call__(self, model, where):
            """Callback entry point: call lazy constraints routine when new
            solutions are found. Stop the optimization if there is an exception in
            user code."""
            if where == GRB.Callback.MIPSOL:
                try:
                    self.eliminate_subtours(model)
                except Exception:
                    logging.exception("Exception occurred in MIPSOL callback")
                    model.terminate()

        def eliminate_subtours(self, model):
            """Extract the current solution, check for subtours, and formulate lazy
            constraints to cut off the current solution if subtours are found.
            Assumes we are at MIPSOL."""
            values = model.cbGetSolution(self.x)
            edges = [(i, j) for (i, j), v in values.items() if v > 0.5]
            tour = shortest_subtour(edges)
            if len(tour) < len(self.nodes):
                # add subtour elimination constraint for every pair of cities in tour
                model.cbLazy(
                    gp.quicksum(self.x[i, j] for i, j in combinations(tour, 2))
                    <= len(tour) - 1
                )

    def solve_tsp(nodes, distances):
        """
        Solve a dense symmetric TSP using the following base formulation:

        min  sum_ij d_ij x_ij
        s.t. sum_j x_ij == 2   forall i in V
             x_ij binary       forall (i,j) in E

        and subtours eliminated using lazy constraints.
        """

        with gp.Env() as env, gp.Model(env=env) as m:
            # Create variables, and add symmetric keys to the resulting dictionary
            # 'x', such that (i, j) and (j, i) refer to the same variable.
            x = m.addVars(distances.keys(), obj=distances, vtype=GRB.BINARY, name="e")
            x.update({(j, i): v for (i, j), v in x.items()})

            # Create degree 2 constraints
            for i in nodes:
                m.addConstr(gp.quicksum(x[i, j] for j in nodes if i != j) == 2)

            # Optimize model using lazy constraints to eliminate subtours
            m.Params.LazyConstraints = 1
            cb = TSPCallback(nodes, x)
            m.optimize(cb)

            # Extract the solution as a tour
            edges = [(i, j) for (i, j), v in x.items() if v.X > 0.5]
            tour = shortest_subtour(edges)
            assert set(tour) == set(nodes)

            return tour, m.ObjVal

    if __name__ == "__main__":
        # Parse arguments
        if len(sys.argv) < 2:
            print("Usage: tsp.py npoints <seed>")
            sys.exit(0)
        npoints = int(sys.argv[1])
        seed = int(sys.argv[2]) if len(sys.argv) > 2 else 1

        # Create n random points in 2D
        random.seed(seed)
        nodes = list(range(npoints))
        points = [(random.randint(0, 100), random.randint(0, 100)) for i in nodes]

        # Dictionary of Euclidean distance between each pair of points
        distances = {
            (i, j): math.sqrt(sum((points[i][k] - points[j][k]) ** 2 for k in range(2)))
            for i, j in combinations(nodes, 2)
        }

        tour, cost = solve_tsp(nodes, distances)

        print("")
        print(f"Optimal tour: {tour}")
        print(f"Optimal cost: {cost:g}")
        print("")

Help and Feedback

Press `/` or click to search this documentation | Gurobi website | Gurobi community forum
Go to other Gurobi documentation | Gurobi download | Contact support