# Solve an instance of the TSP with Xpress using callbacks.
# Usage: Retrieve an example from
# https://www.math.uwaterloo.ca/tsp/world/countries.html
# and load the TSP instance, then solve it using the Xpress Optimizer
# library with the appropriate callback. Once the optimization is over
# (i.e. the time limit is reached, or we find an optimal solution) the
# optimal tour is displayed using matplotlib.
#
# (C) 1983-2025 Fair Isaac Corporation
import networkx as nx
import xpress as xp
import re
import math
from matplotlib import pyplot as plt
from urllib.request import urlretrieve
# Download instance from TSPLib
#
# Replace with any of the following for a different instance:
#
# ar9152.tsp (9125 nodes)
# bm33708.tsp (33708 nodes)
# ch71009.tsp (71009 nodes)
# dj38.tsp (38 nodes)
# eg7146.tsp (7146 nodes)
# fi10639.tsp (10639 nodes)
# gr9882.tsp (9882 nodes)
# ho14473.tsp (14473 nodes)
# ei8246.tsp (8246 nodes)
# ja9847.tsp (9847 nodes)
# kz9976.tsp (9976 nodes)
# lu980.tsp (980 nodes)
# mo14185.tsp (14185 nodes)
# nu3496.tsp (3496 nodes)
# mu1979.tsp (1979 nodes)
# pm8079.tsp (8079 nodes)
# qa194.tsp (194 nodes)
# rw1621.tsp (1621 nodes)
# sw24978.tsp (24978 nodes)
# tz6117.tsp (6117 nodes)
# uy734.tsp (734 nodes)
# vm22775.tsp (22775 nodes)
# wi29.tsp (29 nodes)
# ym7663.tsp (7663 nodes)
# zi929.tsp (929 nodes)
# ca4663.tsp (4663 nodes)
# it16862.tsp (16862 nodes)
filename = 'wi29.tsp'
urlretrieve('https://www.math.uwaterloo.ca/tsp/world/' + filename, filename)
# Read file consisting of lines of the form "k: x y" where k is the
# point's index while x and y are the coordinates of the point. The
# distances are assumed to be Euclidean.
instance = open(filename, 'r')
coord_section = False
points = {}
G = nx.Graph()
#
# Coordinates of the points in the graph.
#
for line in instance.readlines():
if re.match('NODE_COORD_SECTION.*', line):
coord_section = True
continue
elif re.match('EOF.*', line):
break
if coord_section:
coord = line.split(' ')
index = int(coord[0])
cx = float(coord[1])
cy = float(coord[2])
points[index] = (cx, cy)
G.add_node(index, pos=(cx, cy))
instance.close()
print("Downloaded instance, created graph.")
# Callback for checking if the solution forms a tour.
#
# Returns a tuple (a,b) with:
# a: True if the solution is to be rejected, False otherwise
# b: real cutoff value
def cbpreintsol(prob, G, soltype, cutoff):
"""
Use this function to refuse a solution unless it forms a tour.
"""
# Obtain solution, then start at node 1 to see if the solutions at
# one form a tour. The vector s is binary as this is a preintsol()
# callback.
s = prob.getCallbackSolution()
reject = False
nextnode = 1
tour = []
while nextnode != 1 or len(tour) == 0:
# Find the edge leaving nextnode.
edge = None
for j in V:
if j != nextnode and s[x[nextnode, j].index] > 0.5:
edge = x[nextnode, j]
nextnode = j
break
if edge is None:
break
tour.append(edge)
# If there are n arcs in the loop, the solution is feasible.
if len(tour) < n:
# The tour given by the current solution does not pass through
# all the nodes and is thus infeasible.
# If soltype is non-zero then we reject by setting reject=True.
# If instead soltype is zero then the solution came from an
# integral node. In this case we can reject by adding a cut
# that cuts off that solution. Note that we must NOT set
# reject=True in that case because that would result in just
# dropping the node, no matter whether we add cuts or not.
if soltype != 0:
reject = True
else:
# The solution is infeasible and it was obtained from an integral
# node. In this case we can generate and inject a cut that cuts
# off this solution so that we don't find it again.
# Presolve cut in order to add it to the presolved problem.
colind, rowcoef, drhsp, status = prob.presolveRow(rowtype='L',
origcolind=tour,
origrowcoef=[1] * len(tour),
origrhs=len(tour) - 1)
# Since mipdualreductions=0, presolving the cut must succeed, and
# the cut should never be relaxed as this would imply that it did
# not cut off a subtour.
assert status == 0
prob.addCuts(cuttype=[1],
rowtype=['L'],
rhs=[drhsp],
start=[0, len(colind)],
colind=colind,
cutcoef=rowcoef)
# To accept the cutoff, return second element of tuple as None.
return (reject, None)
# Formulate problem, set callback function and solve.
n = len(points) # Number of nodes.
V = range(1, n+1) # Set of nodes.
# Set of arcs (i.e. all pairs since it is a complete graph).
A = [(i, j) for i in V for j in V if i != j]
p = xp.problem()
x = {(i, j): p.addVariable(name='x_{0}_{1}'.format(i, j),
vartype=xp.binary) for (i, j) in A}
conservation_in = [xp.Sum(x[i, j] for j in V if j != i) == 1 for i in V]
conservation_out = [xp.Sum(x[j, i] for j in V if j != i) == 1 for i in V]
p.addConstraint(conservation_in, conservation_out)
# Objective function: total distance travelled.
p.setObjective(xp.Sum(math.sqrt((points[i][0] - points[j][0])**2 +
(points[i][1] - points[j][1])**2) * x[i, j]
for (i, j) in A))
# Should find a reasonable solution within 20 seconds.
p.controls.timelimit = 20
p.addPreIntsolCallback(cbpreintsol, G, 1)
# Disable dual reductions (in order not to cut optimal solutions)
# and nonlinear reductions, in order to be able to presolve the
# cuts.
p.controls.mipdualreductions = 0
p.optimize()
if p.attributes.solstatus not in [xp.SolStatus.OPTIMAL, xp.SolStatus.FEASIBLE]:
print("Solve status:", p.attributes.solvestatus.name)
print("Solution status:", p.attributes.solstatus.name)
else:
# Read solution and store it in the graph.
sol = p.getSolution()
try:
for (i, j) in A:
if sol[x[i, j].index] > 0.5:
G.add_edge(i, j)
# Display best tour found.
pos = nx.get_node_attributes(G, 'pos')
nx.draw(G, points) # Create a graph with the tour.
plt.show() # Display it interactively.
except:
print('Could not draw solution')
|
# TSP example using numpy functions (for efficiency).
#
# (C) 1983-2025 Fair Isaac Corporation
import xpress as xp
import numpy as np
def cb_preintsol(prob, data, soltype, cutoff):
'''Callback for checking if solution is acceptable.'''
n = data
xsol = prob.getCallbackSolution()
xsolf = np.array(xsol) # flattened
xsol = xsolf.reshape(n, n) # matrix-shaped
nextc = np.argmax(xsol, axis=1)
i = 0
ncities = 1
# Scan cities in order until we get back to 0 or the solution is
# wrong and we're diverging.
while nextc[i] != 0 and ncities < n:
ncities += 1
i = nextc[i]
reject = False
if ncities < n:
# The tour given by the current solution does not pass through
# all the nodes and is thus infeasible.
# If soltype is non-zero then we reject by setting reject=True.
# If instead soltype is zero then the solution came from an
# integral node. In this case we can reject by adding a cut
# that cuts off that solution. Note that we must NOT set
# reject=True in that case because that would result in just
# dropping the node, no matter whether we add cuts or not.
if soltype != 0:
reject = True
else:
# Obtain an order by checking the maximum of the variable matrix
# for each row.
unchecked = np.zeros(n)
ngroup = 0
# Initialize the vectors to be passed to addcuts.
cut_mstart = [0]
cut_ind = []
cut_coe = []
cut_rhs = []
nnz = 0
ncuts = 0
while np.min(unchecked) == 0 and ngroup <= n:
'''Seek a tour
'''
ngroup += 1
firstcity = np.argmin(unchecked)
assert (unchecked[firstcity] == 0)
i = firstcity
ncities = 0
# Scan cities in order.
while True:
unchecked[i] = ngroup # Mark city i with its new group, to be used in addcut.
ncities += 1
i = nextc[i]
if i == firstcity or ncities > n + 1:
break
assert ncities < n # We know solution is infeasible.
# Unchecked[unchecked == ngroup] marks nodes to be made part of
# subtour elimination inequality.
# Find indices of current subtour. S is the set of nodes
# traversed by the subtour, compS is its complement.
S = np.where(unchecked == ngroup)[0].tolist()
compS = np.where(unchecked != ngroup)[0].tolist()
indices = [i*n+j for i in S for j in compS]
# Check if solution violates the cut, and if so add the cut to
# the list.
if sum(xsolf[i] for i in indices) < 1 - 1e-3:
# Presolve cut in order to add it to the presolved problem.
mcolsp, dvalp, drhsp, status = prob.presolveRow(rowtype='G',
origcolind=indices,
origrowcoef=np.ones(len(indices)),
origrhs=1)
# Since mipdualreductions=0, presolving the cut must succeed, and the cut should
# never be relaxed as this would imply that it did not cut off a subtour.
assert status == 0
nnz += len(mcolsp)
ncuts += 1
cut_ind.extend(mcolsp)
cut_coe.extend(dvalp)
cut_rhs.append(drhsp)
cut_mstart.append(nnz)
if ncuts > 0:
assert (len(cut_mstart) == ncuts + 1)
assert (len(cut_ind) == nnz)
prob.addCuts(cuttype=[0] * ncuts,
rowtype=['G'] * ncuts,
rhs=cut_rhs,
start=cut_mstart,
colind=cut_ind,
cutcoef=cut_coe)
return (reject, None)
def print_sol(p, n):
'''Print the solution: order of nodes and cost.'''
xsol = np.array(p.getSolution()).reshape(n, n)
nextc = np.argmax(xsol, axis=1)
i = 0
# Scan cities in order.
tour = []
while i != 0 or len(tour) == 0:
tour.append(str(i))
i = nextc[i]
print('->'.join(tour), '->0; cost: ', p.attributes.objval, sep='')
def create_initial_tour(n):
'''Returns a permuted trivial solution 0->1->2->...->(n-1)->0.'''
sol = np.zeros((n, n))
p = np.random.permutation(n)
for i in range(n):
sol[p[i], p[(i + 1) % n]] = 1
return sol.flatten()
def solve_opttour():
'''Create a random TSP problem.'''
n = 50
CITIES = range(n) # Set of cities: 0..n-1.
np.random.seed(3)
X = 100 * np.random.rand(n)
Y = 100 * np.random.rand(n)
# Compute distance matrix.
dist = np.ceil(np.sqrt ((X.reshape(n,1) - X.reshape(1,n))**2 +
(Y.reshape(n,1) - Y.reshape(1,n))**2))
p = xp.problem()
# Create variables as a square matrix of binary variables.
fly = xp.array([p.addVariable(vartype=xp.binary, name='x_{0}_{1}'.format(i,j))
for i in CITIES for j in CITIES]).reshape(n,n)
# Degree constraints.
p.addConstraint(xp.Sum(fly[i,:]) - fly[i,i] == 1 for i in CITIES)
p.addConstraint(xp.Sum(fly[:,i]) - fly[i,i] == 1 for i in CITIES)
# Fix diagonals (i.e. city X -> city X) to zero.
p.addConstraint(fly[i,i] == 0 for i in CITIES)
# Objective function.
p.setObjective (xp.Sum((dist * fly).flatten()))
# Add callbacks.
p.addPreIntsolCallback(cb_preintsol, n)
# Disable dual reductions (in order not to cut optimal solutions)
# and nonlinear reductions, in order to be able to presolve the
# cuts.
p.controls.mipdualreductions = 0
# Create 10 trivial solutions: simple tour 0->1->2...->n->0
# randomly permuted.
for k in range(10):
InitTour = create_initial_tour(n)
p.addMipSol(solval=InitTour, name="InitTour_{}".format(k))
p.optimize()
if p.attributes.solstatus not in [xp.SolStatus.OPTIMAL, xp.SolStatus.FEASIBLE]:
print("Solve status:", p.attributes.solvestatus.name)
print("Solution status:", p.attributes.solstatus.name)
else:
print_sol(p,n) # Print solution and cost.
if __name__ == '__main__':
solve_opttour()
|
