Try our new documentation site (beta).

Filter Content By
Version

### tsp.py

#!/usr/bin/env python3.11

# Copyright 2024, 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

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)

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("")


Choose the evaluation license that fits you best, and start working with our Expert Team for technical guidance and support.

Get a free, full-featured license of the Gurobi Optimizer to experience the performance, support, benchmarking and tuning services we provide as part of our product offering.
Gurobi supports the teaching and use of optimization within academic institutions. We offer free, full-featured copies of Gurobi for use in class, and for research.
##### Cloud Trial

Request free trial hours, so you can see how quickly and easily a model can be solved on the cloud.