#!/bin/env python # # Solve an instance of the TSP with Xpress using callbacks # # 2016 (C) FICO # # # Retrieve an example from # # http://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. # import networkx as nx import xpress as xp import re, math, sys from matplotlib import pyplot as plt if sys.version_info >= (3,): # Import with Python 3 import urllib.request as ul else: # Use Python 2 import urllib as ul # # 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 = 'dj38.tsp' ul.urlretrieve ('http://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 check_tour (prob, G, isheuristic, 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 = [] try: prob.getlpsol (s, None, None, None) except: print ("Can't get LP solution at this node, bailing out") return 1 # bail out orignode = 1 nextnode = 1 card = 0 while nextnode != orignode or card == 0: FS = [j for j in V if j != nextnode and s[prob.getIndex (x[nextnode,j])] == 1] # forward star card += 1 if len (FS) < 1: return (True, None) # reject solution if we can't close the loop nextnode = FS [0] # If there are n arcs in the loop, the solution is feasible return (card < n, None) # to accept the cutoff, return second element of tuple as None # # Callback for adding subtour elimination constraints # # Return nonzero if the node is infeasible, 0 otherwise # def eliminate_subtour (prob, G): """ Function to insert subtour elimination constraints """ s = [] # initialize s to an empty list to provide it as an output parameter try: prob.getlpsol (s, None, None, None) except: print ("Can't get LP solution at this node, bailing out") return 0 # bail out # Starting from node 1, gather all connected nodes of a loop in # set M. if M == V, then the solution is valid if integer, # otherwise add a subtour elimination constraint orignode = 1 nextnode = 1 connset = [] while nextnode != orignode or len (connset) == 0: connset.append (nextnode) FS = [j for j in V if j != nextnode and s [prob.getIndex (x [nextnode, j])] == 1] # forward star if len (FS) < 1: return 0 nextnode = FS [0] if len (connset) < n: # Add a subtour elimination using the nodes in connset (or, if # card (connset) > n/2, its complement) if len (connset) <= n/2: columns = [x[i,j] for i in connset for j in connset if i != j] nArcs = len (connset) else: columns = [x[i,j] for i in V for j in V if not i in connset and not j in connset and i != j] nArcs = n - len (connset) nTerms = len (columns) prob.addcuts ([1], ['L'], [nArcs - 1], [0, nTerms], columns, [1] * nTerms) return 0 # return nonzero for infeasible # # Formulate problem, set callback function and solve # n = len (points) # number of nodes V = range (1, n+1) # set of nodes A = [(i,j) for i in V for j in V if i != j] # set of arcs (i.e. all pairs since it is a complete graph) x = {(i,j): xp.var (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 = xp.problem () p.addVariable (x) p.addConstraint (conservation_in, conservation_out) xind = {(i,j): p.getIndex (x[i,j]) for (i,j) in x.keys ()} # 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)) p.controls.maxtime = -200 # negative for "stop even if no solution is found" p.addcboptnode (eliminate_subtour, G, 1) p.addcbpreintsol (check_tour, G, 1) p.solve () try: sol = p.getSolution () # Read solution and store it in the graph for (i,j) in A: if sol [p.getIndex (x[i,j])] > 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 get a solution in the time limit or could not draw solution, bailing out')