#!/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.importsysimportloggingimportmathimportrandomfromcollectionsimportdefaultdictfromitertoolsimportcombinationsimportgurobipyasgpfromgurobipyimportGRBdefshortest_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 neighboursnode_neighbors=defaultdict(list)fori,jinedges:node_neighbors[i].append(j)assertall(len(neighbors)==2forneighborsinnode_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=Nonewhileunvisited:cycle=[]neighbors=list(unvisited)whileneighbors:current=neighbors.pop()cycle.append(current)unvisited.remove(current)neighbors=[jforjinnode_neighbors[current]ifjinunvisited]ifshortestisNoneorlen(cycle)<len(shortest):shortest=cycleassertshortestisnotNonereturnshortestclassTSPCallback:"""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=nodesself.x=xdef__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."""ifwhere==GRB.Callback.MIPSOL:try:self.eliminate_subtours(model)exceptException:logging.exception("Exception occurred in MIPSOL callback")model.terminate()defeliminate_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),vinvalues.items()ifv>0.5]tour=shortest_subtour(edges)iflen(tour)<len(self.nodes):# add subtour elimination constraint for every pair of cities in tourmodel.cbLazy(gp.quicksum(self.x[i,j]fori,jincombinations(tour,2))<=len(tour)-1)defsolve_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.    """withgp.Env()asenv,gp.Model(env=env)asm:# 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):vfor(i,j),vinx.items()})# Create degree 2 constraintsforiinnodes:m.addConstr(gp.quicksum(x[i,j]forjinnodesifi!=j)==2)# Optimize model using lazy constraints to eliminate subtoursm.Params.LazyConstraints=1cb=TSPCallback(nodes,x)m.optimize(cb)# Extract the solution as a touredges=[(i,j)for(i,j),vinx.items()ifv.X>0.5]tour=shortest_subtour(edges)assertset(tour)==set(nodes)returntour,m.ObjValif__name__=="__main__":# Parse argumentsiflen(sys.argv)<2:print("Usage: tsp.py npoints <seed>")sys.exit(0)npoints=int(sys.argv[1])seed=int(sys.argv[2])iflen(sys.argv)>2else1# Create n random points in 2Drandom.seed(seed)nodes=list(range(npoints))points=[(random.randint(0,100),random.randint(0,100))foriinnodes]# Dictionary of Euclidean distance between each pair of pointsdistances={(i,j):math.sqrt(sum((points[i][k]-points[j][k])**2forkinrange(2)))fori,jincombinations(nodes,2)}tour,cost=solve_tsp(nodes,distances)print("")print(f"Optimal tour:{tour}")print(f"Optimal cost:{cost:g}")print("")