| |||||||||||
| |||||||||||
Python notebooks Description Python notebooks available in the GitHub repository python-notebooks.
Source Files By clicking on a file name, a preview is opened at the bottom of this page. tsp_callbacks.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# **Solving a TSP problem using callbacks**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***tsp_callbacks.ipynb***\n",
"\n",
"Solve an instance of the [Traveling Salesperson Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem) with Xpress using callbacks and *NumPy* arrays. \n",
"\n",
"Generate random TSP data, then solve the problem using the FICO® 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*.\n",
"\n",
"*This example requires a full license of the FICO® Xpress Optimizer. Click on [this link](https://www.fico.com/en/fico-xpress-trial-and-licensing-options) for more information about trial and licensing options.*\n",
"\n",
"© Copyright 2025 Fair Isaac Corporation\n",
"\n",
"Licensed under the Apache License, Version 2.0 (the \"License\"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.\n",
" \n",
"Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.\n",
"\n",
"This example uses FICO® Xpress software. By running it, you agree to the Community License terms of the [Xpress Shrinkwrap License Agreement](https://community.fico.com/s/contentdocument/06980000002h0i5AAA) with respect to the FICO® Xpress software. See the [licensing options](https://www.fico.com/en/fico-xpress-trial-and-licensing-options) overview for additional details and information about obtaining a paid license."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Install the xpress package\n",
"%pip install -q xpress"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Standard formulation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's start with the standard formulation, where binary variables $use_{ij} \\in \\{0,1\\}, \\forall i,j \\in \\mathcal{N}$ represent the decision of whether the tour uses the arc $(i,j)$ (i.e. if we go from city $i$ to $j$) or not. An optimal tour can be found by solving:\n",
"$$\n",
"\\min \\sum_{i,j \\in \\mathcal{N}} dist_{ij} \\cdot use_{ij}\n",
"$$\n",
"\n",
"subject to:\n",
"\n",
"* We have to enter and leave every city, and source node cannot be the same as the destination node: \n",
"$$\n",
"\\sum_{j \\in \\mathcal{N}} use_{ij} = 1, \\quad \\forall i \\in \\mathcal{N} \\\\\n",
"\\sum_{j \\in \\mathcal{N}} use_{ji} = 1, \\quad \\forall i \\in \\mathcal{N} \\\\\n",
"use_{ii} = 0, \\quad \\forall i \\in \\mathcal{N}\n",
"$$\n",
"\n",
"Where $dist_{ij}, \\forall i,j \\in \\mathcal{N}$ represents the distance (or cost) associated with traveling on an arc $(i,j)$. We overlook subtour elimination constraints for now."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xpress as xp\n",
"import numpy as np\n",
"import networkx as nx \n",
"import matplotlib.pyplot as plt\n",
"import itertools\n",
"\n",
"rndseed = 10\n",
"\n",
"np.random.seed(rndseed)\n",
"\n",
"'''Create random TSP problem data'''\n",
"n = 17\n",
"CITIES = range(n) # set of cities: 0..n-1\n",
"\n",
"np.random.seed(0)\n",
"\n",
"X = 100 * np.random.rand(n)\n",
"Y = 100 * np.random.rand(n)\n",
"\n",
"XY = (X, Y)\n",
"\n",
"# Compute distance matrix\n",
"dist = np.ceil(np.sqrt ((X.reshape(n,1) - X.reshape(1,n))**2 +\n",
" (Y.reshape(n,1) - Y.reshape(1,n))**2))\n",
"\n",
"# Create problem\n",
"p = xp.problem()\n",
"\n",
"# Create variables as a square matrix of binary variables. Note\n",
"# the use of p.addVariables to ensure NumPy uses the Xpress operations for handling these vectors.\n",
"use = p.addVariables(n,n, vartype=xp.binary, name='x')\n",
"\n",
"# Degree constraints\n",
"p.addConstraint(xp.Sum(use[i,:]) == 1 for i in CITIES)\n",
"p.addConstraint(xp.Sum(use[:,i]) == 1 for i in CITIES)\n",
"\n",
"# Fix diagonals (i.e. city X -> city X) to zero\n",
"p.addConstraint(use[i,i] == 0 for i in CITIES)\n",
"\n",
"# Objective function\n",
"p.setObjective(xp.Sum((dist * use).flatten()))\n",
"\n",
"p.optimize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's use the *matplotlib* and *networkx* Python packages to visualize the solution obtained by Xpress and check if there are any subtours:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Visualize solution using networkx and matplotlib\n",
"def plot_sol(p):\n",
"\n",
" # collect the edges: if the value of x[i,j] is 1, then the edge (i,j) is in the solution\n",
" edges = []\n",
" sol = p.getSolution(use) # get solution values\n",
" for (i,j) in itertools.permutations(CITIES, 2):\n",
" if sol[i,j] > 0.5: # variable is binary so > 0.5 --> is 1\n",
" edges.append( (i,j) )\n",
"\n",
" # create dictionary with coordinates for each node to plot the graph\n",
" xy = {}\n",
" for i in range(n):\n",
" xy[i] = (X[i], Y[i])\n",
"\n",
" # make figure look nicer\n",
" plt.figure(figsize=(10,6), dpi=100)\n",
"\n",
" # create empty graph\n",
" optgraph = nx.Graph()\n",
"\n",
" # add edges\n",
" optgraph.add_edges_from(edges)\n",
"\n",
" # draw the cities, with labels in the position xy (see when we read the instance)\n",
" nx.draw(optgraph, node_size=300, pos=xy, with_labels=True, node_color='lightblue')\n",
"\n",
" # show drawing\n",
" plt.show()\n",
"\n",
"# Plot solution\n",
"if p.attributes.solstatus not in [xp.SolStatus.OPTIMAL, xp.SolStatus.FEASIBLE]:\n",
" print(\"Solve status:\", p.attributes.solvestatus.name)\n",
" print(\"Solution status:\", p.attributes.solstatus.name)\n",
"else:\n",
" plot_sol(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, the solution contains subtours, and therefore it is not yet a valid solution for the problem yet. We need to add constraints that eliminate any such subtour.\n",
"\n",
"A straightforward approach for formulating subtour elimination constraints is to add a constraint for each subtour contained in the set of all subtours $\\mathcal{S}$ as follows:\n",
"\n",
"$$\n",
"\\sum_{i,j \\in S} use_{ij} \\leq |\\mathcal{S}| - 1 \\quad \\forall \\mathcal{S} \\subsetneq \\mathcal{N},\\ |\\mathcal{S}| \\geq 2\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Add subtour elimination constraints, solve and plot solution\n",
"p.addConstraint(xp.Sum(use[i,j] for i in subset for j in subset) <= len(subset) - 1 \n",
" for L in range(2,len(CITIES)) \n",
" for subset in itertools.combinations(CITIES, L)) \n",
"\n",
"p.optimize()\n",
"\n",
"# Plot solution\n",
"plot_sol(p) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is now valid and complete, since an Hamiltonian tour (no subtours) has been formed. However, you could experience that the problem building and solving times are rather high for a TSP with only 17 cities.\n",
"\n",
"This is due to the complete enumeration of all possible subtours, whose number grows exponentially with the instance size making the problem creation and solving times prohibitively long when problems grow just somewhat larger."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Formulation using the Miller, Tucker, Zemlin subtour elimination constraints"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we will model the same problem but using the Miller, Tucker, Zemlin subtour elimination constraints, which requires a new set of continuous variables:\n",
"\n",
"$step_i$ = step at which node $i$ is visited, $\\forall i = {2,...,|\\mathcal{N}|}$\n",
"\n",
"which are used as auxiliary variables for the formulation of subtour elimination constraints. Although a new set of real variables of size $|\\mathcal{N}|$ is introduced, we are able to prevent all subtours by introducing a set of constraints of size $(|\\mathcal{N}|-1)^2$, considerably reducing the problem size and complexity when comparing with the previous formulation:\n",
"\n",
"$$\n",
"\\min \\sum_{i,j \\in \\mathcal{N}} cost_{ij} \\cdot use_{ij}\n",
"$$\n",
"\n",
"subject to:\n",
"\n",
"* We have to enter and leave every city, and source node cannot be the destination node for any move: \n",
"$$\n",
"\\sum_{j \\in \\mathcal{N}} use_{ij} = 1, \\quad \\forall i \\in \\mathcal{N} \\\\\n",
"\\sum_{j \\in \\mathcal{N}} use_{ji} = 1, \\quad \\forall i \\in \\mathcal{N} \\\\\n",
"use_{ii} = 0, \\quad \\forall i \\in \\mathcal{N}\n",
"$$\n",
"* Miller, Tucker, Zemlin subtour elimination constraints:\n",
"$$\n",
"step_j \\geq step_i + 1 - (n-1) \\cdot (1 - use_{ij}), \\forall i,j = {2,..,n}\n",
"$$\n",
"\n",
"with $n = |\\mathcal{N}|$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create problem\n",
"p = xp.problem()\n",
"\n",
"# Create variables as a square matrix of binary variables. Note\n",
"# the use of p.addVariables to ensure NumPy uses the Xpress operations for handling these vectors.\n",
"use = p.addVariables(n,n, vartype=xp.binary, name='x')\n",
"step = p.addVariables(n, name='t')\n",
"\n",
"# Degree constraints\n",
"p.addConstraint(xp.Sum(use[i,:]) == 1 for i in CITIES)\n",
"p.addConstraint(xp.Sum(use[:,i]) == 1 for i in CITIES)\n",
"\n",
"# Fix diagonals (i.e. city X -> city X) to zero\n",
"p.addConstraint(use[i,i] == 0 for i in CITIES)\n",
"\n",
"# Miller, Tucker, Zemlin subtour elimination constraints\n",
"p.addConstraint([step[j] >= step[i] + 1 - (n-1)*(1 - use[i,j]) for i in range(1,n) for j in range(1,n)])\n",
"\n",
"# Objective function\n",
"p.setObjective (xp.Sum((dist * use).flatten()))\n",
"\n",
"p.optimize()\n",
"\n",
"# Plot solution\n",
"if p.attributes.solstatus not in [xp.SolStatus.OPTIMAL, xp.SolStatus.FEASIBLE]:\n",
" print(\"Solve status:\", p.attributes.solvestatus.name)\n",
" print(\"Solution status:\", p.attributes.solstatus.name)\n",
"else:\n",
" plot_sol(p) # print solution and cost"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can confirm by looking at the graph and solver logs, this formulation yields the same solution as the previous one, but solves the problem much faster. Let's try to beat that in the next section!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using callbacks"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's try to speed up the problem creation and solving times even more by using **solver callbacks** to add only those subtour elimination constraints that are really needed during the solving process using the standard formulation.\n",
"\n",
"Callback functions are user–defined routines to be specified to the Optimizer which are called at specific stages during the optimization process, prompting the Optimizer to return to the user's program before continuing with the solution algorithm. They are a useful tool particularly for Mixed Integer Programming (MIP) problems by allowing the user to interact with the solver during the solution search process (Branch-and-Bound). There are specific callbacks for other problem types, too.\n",
"\n",
"[Using callback functions](https://www.fico.com/fico-xpress-optimization/docs/latest/solver/optimizer/python/HTML/chCallbacks.html) is simple: \n",
" * the user first defines a function (say *myfunction*) with a specific signature that is to be run every time the branch-and-bound reaches a well-specified point\n",
" * secondly, the user specifies that the function is to be invoked every time the algorithm reaches a well-specified point calls a function (such as [problem.addcbpreintsol](https://www.fico.com/fico-xpress-optimization/docs/latest/solver/optimizer/python/HTML/problem.addcbpreintsol.html)) with *myfunction* as its argument\n",
" * finally, the user runs the *optimize* command that launches the branch-and-bound, the simplex solver, or the barrier solver; it is while these are run that myfunction is called at the intended points (depending on the selected callback type, check the [list of callbacks accepted by the Xpress Optimizer](https://www.fico.com/fico-xpress-optimization/docs/latest/solver/optimizer/HTML/chapter5.html?scroll=section5002))\n",
"\n",
"In this TSP example, we will relax (leave out) the subtour elimination constraints, and instead add a callback funtion that is triggered every time an integer (feasible) solution to the relaxed problem is found during the B&B process. By using the [problem.addcbpreintsol](https://www.fico.com/fico-xpress-optimization/docs/latest/solver/optimizer/python/HTML/problem.addcbpreintsol.html) solver callback, we will instruct the solver whether to accept (if no subtours) or discard (subtours exist) the solution found. In the latter case we add suitable subtour elimination constraints to the problem by calling the function [problem.addcuts](https://www.fico.com/fico-xpress-optimization/docs/dms2024-03/solver/optimizer/python/HTML/problem.addcuts.html) with the appropriate arguments. Therefore, the user-defined callback function needs to perform two main tasks:\n",
" * checking whether a given solution forms a Hamiltonian tour\n",
" * applying subtour elimination constraints to exclude the current node solution subtours from being found again"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Define callback function\n",
"def cb_preintsol(prob, data, soltype, cutoff):\n",
" '''Callback for checking if solution is acceptable\n",
" prob: Xpress problem object\n",
" data: data object = number of cities\n",
" soltype: type of MIP solution that has been found. 0 - LP relaxation is integer feasible, 1 - MIP solution found by a heuristic, 2 - MIP solution provided by the user\n",
" '''\n",
"\n",
" n = data # number of cities\n",
" xsol=prob.getCallbackSolution(use) # get array with values of x variables in the current context-specific solution\n",
" xsol = np.array(xsol).reshape(n,n) # convert into matrix-shaped np array\n",
" tour = np.argmax(xsol, axis=1) # get indices of the max (non-zero) element in each column = resulting tour\n",
"\n",
" i = 0 # to store next city\n",
" ncities = 1 # to count number of cities visited\n",
"\n",
" # Scan cities in order until we get back to 0 or the solution is wrong.\n",
" while tour[i] != 0 and ncities < n:\n",
" ncities += 1\n",
" i = tour[i]\n",
"\n",
" reject = False # at this point do not reject the solution\n",
" if ncities < n:\n",
" # The tour given by the current solution does not pass through\n",
" # all the cities and thus contains subtours.\n",
" # If soltype is non-zero, the solution was found by a heuristic\n",
" # or provided by the user, so we reject it by setting reject=True without further branching.\n",
" # If instead soltype is zero, the solution came from the LP relaxation\n",
" # of the current B&B node which was found to be integer feasible. In this case\n",
" # we will reject it by adding a cut that cuts off that solution.\n",
" # Note that we must NOT set reject=True in that case because that would result in just\n",
" # dropping the node, no matter whether we add cuts or not.\n",
" if soltype != 0:\n",
" reject = True\n",
" else:\n",
" # Obtain an order by checking the maximum of the variable matrix\n",
" # for each row\n",
" unchecked = np.zeros(n) # np array of zeros of size n to assign each city to a subtour\n",
" nsubtour = 0 # to number and identify subtours\n",
"\n",
" # Initialize the vectors to be passed to problem.addcuts\n",
" cut_mstart = [0]\n",
" cut_ind = []\n",
" cut_coe = []\n",
" cut_rhs = []\n",
"\n",
" nnz = 0 # to mark the start of the variables indices for each cut\n",
" ncuts = 0 # to count number of cuts (subtours)\n",
"\n",
" # while there are cities that are not assigned to a subtour\n",
" while np.min(unchecked) == 0:\n",
" nsubtour += 1 # current subtour id\n",
" firstcity = np.argmin(unchecked) # first city (index) in unchecked that is not assigned a subtour yet (i.e., =0)\n",
" i = firstcity \n",
"\n",
" # Scan cities in order for each subtour\n",
" while True:\n",
" unchecked[i] = nsubtour # mark city i as belonging to current subtour\n",
" i = tour[i] # proceed to next city in the current subtour\n",
"\n",
" if i == firstcity: # if next city is the first city, subtour is finished, stop\n",
" break\n",
"\n",
" # Find indices of all variables with origin in this subtour and destination\n",
" # outside this subtour. We will force the sum of those variables to be >= 1\n",
" # (i.e. subtour will be excluded after the cut is added). S is the set of cities\n",
" # traversed by the subtour, compS is its complement (list of cities not contained in subtour S)\n",
" S = np.where(unchecked == nsubtour)[0].tolist()\n",
" compS = np.where(unchecked != nsubtour)[0].tolist()\n",
" indices = [i*n+j for i in S for j in compS] # convert into decision variables array indices\n",
"\n",
" # We need to presolve new row to adapt its data (rhs, column indices and row coefficients) for the presolved problem (instead of the original problem)\n",
" mcolsp, dvalp = [], [] # lists to be populated by the 'presolverow' method\n",
"\n",
" # Presolve new row (cut) in order to add it to the presolved problem via addcuts\n",
" drhsp, status = prob.presolverow(rowtype='G',\n",
" origcolind=indices,\n",
" origrowcoef=np.ones(len(indices)),\n",
" origrhs=1,\n",
" maxcoefs=prob.attributes.cols,\n",
" colind=mcolsp,\n",
" rowcoef=dvalp)\n",
"\n",
" nnz += len(mcolsp)\n",
" ncuts += 1\n",
"\n",
" cut_ind.extend(mcolsp) # array which will be filled with the columns of the presolved row\n",
" cut_coe.extend(dvalp) # array which will be filled with the coefficients of the presolved row\n",
" cut_rhs.append(drhsp) # array containing the right hand side elements for the cuts\n",
" cut_mstart.append(nnz) # array containing offset into the colind and cutcoef arrays indicating the start of each cut\n",
"\n",
" # add cuts\n",
" if ncuts > 0:\n",
" # call prob.addcuts and pass the necessary arguments to add the subtour elimination constraints (list of cuts)\n",
" prob.addcuts(cuttype=[0] * ncuts,\n",
" rowtype=['G'] * ncuts,\n",
" rhs=cut_rhs,\n",
" start=cut_mstart,\n",
" colind=cut_ind,\n",
" cutcoef=cut_coe) \n",
"\n",
" return (reject, None)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we will solve the problem with the standard formulation (without subtour elimination constraints), but **adding the callback function** and plotting the solution to validate the approach. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create problem\n",
"p = xp.problem()\n",
"\n",
"# Create variables as a square matrix of binary variables. Note\n",
"# the use of p.addVariables to ensure NumPy uses the Xpress operations for handling these vectors.\n",
"use = p.addVariables(n,n, vartype=xp.binary, name='x')\n",
"\n",
"# Degree constraints\n",
"p.addConstraint(xp.Sum(use[i,:]) == 1 for i in CITIES)\n",
"p.addConstraint(xp.Sum(use[:,i]) == 1 for i in CITIES)\n",
"\n",
"# Fix diagonals (i.e. city X -> city X) to zero\n",
"p.addConstraint(use[i,i] == 0 for i in CITIES)\n",
"\n",
"# Objective function\n",
"p.setObjective(xp.Sum((dist * use).flatten()))\n",
"\n",
"# Add callback function\n",
"p.addcbpreintsol(cb_preintsol, n)\n",
"\n",
"# Solve the problem and print solution\n",
"p.optimize()\n",
"\n",
"print(p.attributes.solstatus.name)\n",
"\n",
"if p.attributes.solstatus not in [xp.SolStatus.OPTIMAL, xp.SolStatus.FEASIBLE]:\n",
" print(\"Solve status:\", p.attributes.solvestatus.name)\n",
" print(\"Solution status:\", p.attributes.solstatus.name)\n",
"else:\n",
" plot_sol(p) # plot solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can observe, the problem has been solved even faster than using the MTZ formulation, yielding the same solution and objective value in considerably less time, making the problem even more scalable. \n",
"\n",
"As an optional exercise, you can experiment with increasing the number the number of cities (*n* below) in the cell below and run the previous code cell again to see how high the number of cities <tt>n</tt> can go before the problem becomes computationally too expensive (let's say, takes more than 1 minute to solve to proven optimality)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# Create random TSP problem data\n",
"n = 150\n",
"\n",
"CITIES = range(n) # set of cities: 0..n-1\n",
"\n",
"np.random.seed(0)\n",
"\n",
"X = 100 * np.random.rand(n)\n",
"Y = 100 * np.random.rand(n)\n",
"\n",
"XY = (X, Y)\n",
"\n",
"# Compute distance matrix\n",
"dist = np.ceil(np.sqrt ((X.reshape(n,1) - X.reshape(1,n))**2 +\n",
" (Y.reshape(n,1) - Y.reshape(1,n))**2))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
| |||||||||||
| © Copyright 2025 Fair Isaac Corporation. |
