TSP Examples#
This section includes source code for all of the Gurobi TSP examples.
The same source code can be found in the examples
directory of the
Gurobi distribution.
/* 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 node. The lazy constraint callback
adds new constraints to cut them off.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "gurobi_c.h"
/* Define structure to pass data to the callback function */
struct callback_data {
int n;
};
/* Given an integer-feasible solution 'sol', find the smallest
sub-tour. Result is returned in 'tour', and length is
returned in 'tourlenP'. */
static void
findsubtour(int n,
double *sol,
int *tourlenP,
int *tour)
{
int i, node, len, start;
int bestind, bestlen;
int *seen = NULL;
seen = (int *) malloc(n*sizeof(int));
if (seen == NULL) {
fprintf(stderr, "Out of memory\n");
exit(1);
}
for (i = 0; i < n; i++)
seen[i] = 0;
start = 0;
bestlen = n+1;
bestind = -1;
while (start < n) {
for (node = 0; node < n; node++)
if (seen[node] == 0)
break;
if (node == n)
break;
for (len = 0; len < n; len++) {
tour[start+len] = node;
seen[node] = 1;
for (i = 0; i < n; i++) {
if (sol[node*n+i] > 0.5 && !seen[i]) {
node = i;
break;
}
}
if (i == n) {
len++;
if (len < bestlen) {
bestlen = len;
bestind = start;
}
start += len;
break;
}
}
}
for (i = 0; i < bestlen; i++)
tour[i] = tour[bestind+i];
*tourlenP = bestlen;
free(seen);
}
/* Subtour elimination callback. Whenever a feasible solution is found,
find the shortest subtour, and add a subtour elimination constraint
if that tour doesn't visit every node. */
int __stdcall
subtourelim(GRBmodel *model,
void *cbdata,
int where,
void *usrdata)
{
struct callback_data *mydata = (struct callback_data *) usrdata;
int n = mydata->n;
int *tour = NULL;
double *sol = NULL;
int i, j, len, nz;
int error = 0;
if (where == GRB_CB_MIPSOL) {
sol = (double *) malloc(n*n*sizeof(double));
tour = (int *) malloc(n*sizeof(int));
if (sol == NULL || tour == NULL) {
fprintf(stderr, "Out of memory\n");
exit(1);
}
GRBcbget(cbdata, where, GRB_CB_MIPSOL_SOL, sol);
findsubtour(n, sol, &len, tour);
if (len < n) {
int *ind = NULL;
double *val = NULL;
ind = (int *) malloc(len*(len-1)/2*sizeof(int));
val = (double *) malloc(len*(len-1)/2*sizeof(double));
if (ind == NULL || val == NULL) {
fprintf(stderr, "Out of memory\n");
exit(1);
}
/* Add subtour elimination constraint */
nz = 0;
for (i = 0; i < len; i++)
for (j = i+1; j < len; j++)
ind[nz++] = tour[i]*n+tour[j];
for (i = 0; i < nz; i++)
val[i] = 1.0;
error = GRBcblazy(cbdata, nz, ind, val, GRB_LESS_EQUAL, len-1);
free(ind);
free(val);
}
free(sol);
free(tour);
}
return error;
}
/* Euclidean distance between points 'i' and 'j'. */
static double
distance(double *x,
double *y,
int i,
int j)
{
double dx = x[i] - x[j];
double dy = y[i] - y[j];
return sqrt(dx*dx + dy*dy);
}
int
main(int argc,
char *argv[])
{
GRBenv *env = NULL;
GRBmodel *model = NULL;
int i, j, len, n, solcount;
int error = 0;
char name[100];
double *x = NULL;
double *y = NULL;
int *ind = NULL;
double *val = NULL;
struct callback_data mydata;
if (argc < 2) {
fprintf(stderr, "Usage: tsp_c size\n");
exit(1);
}
n = atoi(argv[1]);
if (n == 0) {
fprintf(stderr, "Argument must be a positive integer.\n");
} else if (n > 100) {
printf("It will be a challenge to solve a TSP this large.\n");
}
x = (double *) malloc(n*sizeof(double));
y = (double *) malloc(n*sizeof(double));
ind = (int *) malloc(n*sizeof(int));
val = (double *) malloc(n*sizeof(double));
if (x == NULL || y == NULL || ind == NULL || val == NULL) {
fprintf(stderr, "Out of memory\n");
exit(1);
}
/* Create random points */
for (i = 0; i < n; i++) {
x[i] = ((double) rand())/RAND_MAX;
y[i] = ((double) rand())/RAND_MAX;
}
/* Create environment */
error = GRBloadenv(&env, "tsp.log");
if (error) goto QUIT;
/* Create an empty model */
error = GRBnewmodel(env, &model, "tsp", 0, NULL, NULL, NULL, NULL, NULL);
if (error) goto QUIT;
/* Add variables - one for every pair of nodes */
/* Note: If edge from i to j is chosen, then x[i*n+j] = x[j*n+i] = 1. */
/* The cost is split between the two variables. */
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
sprintf(name, "x_%d_%d", i, j);
error = GRBaddvar(model, 0, NULL, NULL, distance(x, y, i, j)/2,
0.0, 1.0, GRB_BINARY, name);
if (error) goto QUIT;
}
}
/* Degree-2 constraints */
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
ind[j] = i*n+j;
val[j] = 1.0;
}
sprintf(name, "deg2_%d", i);
error = GRBaddconstr(model, n, ind, val, GRB_EQUAL, 2, name);
if (error) goto QUIT;
}
/* Forbid edge from node back to itself */
for (i = 0; i < n; i++) {
error = GRBsetdblattrelement(model, GRB_DBL_ATTR_UB, i*n+i, 0);
if (error) goto QUIT;
}
/* Symmetric TSP */
for (i = 0; i < n; i++) {
for (j = 0; j < i; j++) {
ind[0] = i*n+j;
ind[1] = i+j*n;
val[0] = 1;
val[1] = -1;
error = GRBaddconstr(model, 2, ind, val, GRB_EQUAL, 0, NULL);
if (error) goto QUIT;
}
}
/* Set callback function */
mydata.n = n;
error = GRBsetcallbackfunc(model, subtourelim, (void *) &mydata);
if (error) goto QUIT;
/* Must set LazyConstraints parameter when using lazy constraints */
error = GRBsetintparam(GRBgetenv(model), GRB_INT_PAR_LAZYCONSTRAINTS, 1);
if (error) goto QUIT;
/* Optimize model */
error = GRBoptimize(model);
if (error) goto QUIT;
/* Extract solution */
error = GRBgetintattr(model, GRB_INT_ATTR_SOLCOUNT, &solcount);
if (error) goto QUIT;
if (solcount > 0) {
int *tour = NULL;
double *sol = NULL;
sol = (double *) malloc(n*n*sizeof(double));
tour = (int *) malloc(n*sizeof(int));
if (sol == NULL || tour == NULL) {
fprintf(stderr, "Out of memory\n");
exit(1);
}
error = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n*n, sol);
if (error) goto QUIT;
/* Print tour */
findsubtour(n, sol, &len, tour);
printf("Tour: ");
for (i = 0; i < len; i++)
printf("%d ", tour[i]);
printf("\n");
free(tour);
free(sol);
}
QUIT:
/* Free data */
free(x);
free(y);
free(ind);
free(val);
/* Error reporting */
if (error) {
printf("ERROR: %s\n", GRBgeterrormsg(env));
exit(1);
}
/* Free model */
GRBfreemodel(model);
/* Free environment */
GRBfreeenv(env);
return 0;
}
/* 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 node. The lazy constraint callback
adds new constraints to cut them off. */
#include "gurobi_c++.h"
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <sstream>
using namespace std;
string itos(int i) {stringstream s; s << i; return s.str(); }
double distance(double* x, double* y, int i, int j);
void findsubtour(int n, double** sol, int* tourlenP, int* tour);
// Subtour elimination callback. Whenever a feasible solution is found,
// find the smallest subtour, and add a subtour elimination constraint
// if the tour doesn't visit every node.
class subtourelim: public GRBCallback
{
public:
GRBVar** vars;
int n;
subtourelim(GRBVar** xvars, int xn) {
vars = xvars;
n = xn;
}
protected:
void callback() {
try {
if (where == GRB_CB_MIPSOL) {
// Found an integer feasible solution - does it visit every node?
double **x = new double*[n];
int *tour = new int[n];
int i, j, len;
for (i = 0; i < n; i++)
x[i] = getSolution(vars[i], n);
findsubtour(n, x, &len, tour);
if (len < n) {
// Add subtour elimination constraint
GRBLinExpr expr = 0;
for (i = 0; i < len; i++)
for (j = i+1; j < len; j++)
expr += vars[tour[i]][tour[j]];
addLazy(expr <= len-1);
}
for (i = 0; i < n; i++)
delete[] x[i];
delete[] x;
delete[] tour;
}
} catch (GRBException e) {
cout << "Error number: " << e.getErrorCode() << endl;
cout << e.getMessage() << endl;
} catch (...) {
cout << "Error during callback" << endl;
}
}
};
// Given an integer-feasible solution 'sol', find the smallest
// sub-tour. Result is returned in 'tour', and length is
// returned in 'tourlenP'.
void
findsubtour(int n,
double** sol,
int* tourlenP,
int* tour)
{
bool* seen = new bool[n];
int bestind, bestlen;
int i, node, len, start;
for (i = 0; i < n; i++)
seen[i] = false;
start = 0;
bestlen = n+1;
bestind = -1;
node = 0;
while (start < n) {
for (node = 0; node < n; node++)
if (!seen[node])
break;
if (node == n)
break;
for (len = 0; len < n; len++) {
tour[start+len] = node;
seen[node] = true;
for (i = 0; i < n; i++) {
if (sol[node][i] > 0.5 && !seen[i]) {
node = i;
break;
}
}
if (i == n) {
len++;
if (len < bestlen) {
bestlen = len;
bestind = start;
}
start += len;
break;
}
}
}
for (i = 0; i < bestlen; i++)
tour[i] = tour[bestind+i];
*tourlenP = bestlen;
delete[] seen;
}
// Euclidean distance between points 'i' and 'j'.
double
distance(double* x,
double* y,
int i,
int j)
{
double dx = x[i]-x[j];
double dy = y[i]-y[j];
return sqrt(dx*dx+dy*dy);
}
int
main(int argc,
char *argv[])
{
if (argc < 2) {
cout << "Usage: tsp_c++ size" << endl;
return 1;
}
int n = atoi(argv[1]);
double* x = new double[n];
double* y = new double[n];
int i;
for (i = 0; i < n; i++) {
x[i] = ((double) rand())/RAND_MAX;
y[i] = ((double) rand())/RAND_MAX;
}
GRBEnv *env = NULL;
GRBVar **vars = NULL;
vars = new GRBVar*[n];
for (i = 0; i < n; i++)
vars[i] = new GRBVar[n];
try {
int j;
env = new GRBEnv();
GRBModel model = GRBModel(*env);
// Must set LazyConstraints parameter when using lazy constraints
model.set(GRB_IntParam_LazyConstraints, 1);
// Create binary decision variables
for (i = 0; i < n; i++) {
for (j = 0; j <= i; j++) {
vars[i][j] = model.addVar(0.0, 1.0, distance(x, y, i, j),
GRB_BINARY, "x_"+itos(i)+"_"+itos(j));
vars[j][i] = vars[i][j];
}
}
// Degree-2 constraints
for (i = 0; i < n; i++) {
GRBLinExpr expr = 0;
for (j = 0; j < n; j++)
expr += vars[i][j];
model.addConstr(expr == 2, "deg2_"+itos(i));
}
// Forbid edge from node back to itself
for (i = 0; i < n; i++)
vars[i][i].set(GRB_DoubleAttr_UB, 0);
// Set callback function
subtourelim cb = subtourelim(vars, n);
model.setCallback(&cb);
// Optimize model
model.optimize();
// Extract solution
if (model.get(GRB_IntAttr_SolCount) > 0) {
double **sol = new double*[n];
for (i = 0; i < n; i++)
sol[i] = model.get(GRB_DoubleAttr_X, vars[i], n);
int* tour = new int[n];
int len;
findsubtour(n, sol, &len, tour);
assert(len == n);
cout << "Tour: ";
for (i = 0; i < len; i++)
cout << tour[i] << " ";
cout << endl;
for (i = 0; i < n; i++)
delete[] sol[i];
delete[] sol;
delete[] tour;
}
} catch (GRBException e) {
cout << "Error number: " << e.getErrorCode() << endl;
cout << e.getMessage() << endl;
} catch (...) {
cout << "Error during optimization" << endl;
}
for (i = 0; i < n; i++)
delete[] vars[i];
delete[] vars;
delete[] x;
delete[] y;
delete env;
return 0;
}
/* 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 node. The lazy constraint callback
// adds new constraints to cut them off.
using System;
using Gurobi;
class tsp_cs : GRBCallback {
private GRBVar[,] vars;
public tsp_cs(GRBVar[,] xvars) {
vars = xvars;
}
// Subtour elimination callback. Whenever a feasible solution is found,
// find the smallest subtour, and add a subtour elimination
// constraint if the tour doesn't visit every node.
protected override void Callback() {
try {
if (where == GRB.Callback.MIPSOL) {
// Found an integer feasible solution - does it visit every node?
int n = vars.GetLength(0);
int[] tour = findsubtour(GetSolution(vars));
if (tour.Length < n) {
// Add subtour elimination constraint
GRBLinExpr expr = 0;
for (int i = 0; i < tour.Length; i++)
for (int j = i+1; j < tour.Length; j++)
expr.AddTerm(1.0, vars[tour[i], tour[j]]);
AddLazy(expr <= tour.Length-1);
}
}
} catch (GRBException e) {
Console.WriteLine("Error code: " + e.ErrorCode + ". " + e.Message);
Console.WriteLine(e.StackTrace);
}
}
// Given an integer-feasible solution 'sol', return the smallest
// sub-tour (as a list of node indices).
protected static int[] findsubtour(double[,] sol)
{
int n = sol.GetLength(0);
bool[] seen = new bool[n];
int[] tour = new int[n];
int bestind, bestlen;
int i, node, len, start;
for (i = 0; i < n; i++)
seen[i] = false;
start = 0;
bestlen = n+1;
bestind = -1;
node = 0;
while (start < n) {
for (node = 0; node < n; node++)
if (!seen[node])
break;
if (node == n)
break;
for (len = 0; len < n; len++) {
tour[start+len] = node;
seen[node] = true;
for (i = 0; i < n; i++) {
if (sol[node, i] > 0.5 && !seen[i]) {
node = i;
break;
}
}
if (i == n) {
len++;
if (len < bestlen) {
bestlen = len;
bestind = start;
}
start += len;
break;
}
}
}
for (i = 0; i < bestlen; i++)
tour[i] = tour[bestind+i];
System.Array.Resize(ref tour, bestlen);
return tour;
}
// Euclidean distance between points 'i' and 'j'
protected static double distance(double[] x,
double[] y,
int i,
int j) {
double dx = x[i]-x[j];
double dy = y[i]-y[j];
return Math.Sqrt(dx*dx+dy*dy);
}
public static void Main(String[] args) {
if (args.Length < 1) {
Console.WriteLine("Usage: tsp_cs nnodes");
return;
}
int n = Convert.ToInt32(args[0]);
try {
GRBEnv env = new GRBEnv();
GRBModel model = new GRBModel(env);
// Must set LazyConstraints parameter when using lazy constraints
model.Parameters.LazyConstraints = 1;
double[] x = new double[n];
double[] y = new double[n];
Random r = new Random();
for (int i = 0; i < n; i++) {
x[i] = r.NextDouble();
y[i] = r.NextDouble();
}
// Create variables
GRBVar[,] vars = new GRBVar[n, n];
for (int i = 0; i < n; i++) {
for (int j = 0; j <= i; j++) {
vars[i, j] = model.AddVar(0.0, 1.0, distance(x, y, i, j),
GRB.BINARY, "x"+i+"_"+j);
vars[j, i] = vars[i, j];
}
}
// Degree-2 constraints
for (int i = 0; i < n; i++) {
GRBLinExpr expr = 0;
for (int j = 0; j < n; j++)
expr.AddTerm(1.0, vars[i, j]);
model.AddConstr(expr == 2.0, "deg2_"+i);
}
// Forbid edge from node back to itself
for (int i = 0; i < n; i++)
vars[i, i].UB = 0.0;
model.SetCallback(new tsp_cs(vars));
model.Optimize();
if (model.SolCount > 0) {
int[] tour = findsubtour(model.Get(GRB.DoubleAttr.X, vars));
Console.Write("Tour: ");
for (int i = 0; i < tour.Length; i++)
Console.Write(tour[i] + " ");
Console.WriteLine();
}
// Dispose of model and environment
model.Dispose();
env.Dispose();
} catch (GRBException e) {
Console.WriteLine("Error code: " + e.ErrorCode + ". " + e.Message);
Console.WriteLine(e.StackTrace);
}
}
}
/* 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 node. The lazy constraint callback
// adds new constraints to cut them off.
import com.gurobi.gurobi.*;
public class Tsp extends GRBCallback {
private GRBVar[][] vars;
public Tsp(GRBVar[][] xvars) {
vars = xvars;
}
// Subtour elimination callback. Whenever a feasible solution is found,
// find the subtour that contains node 0, and add a subtour elimination
// constraint if the tour doesn't visit every node.
protected void callback() {
try {
if (where == GRB.CB_MIPSOL) {
// Found an integer feasible solution - does it visit every node?
int n = vars.length;
int[] tour = findsubtour(getSolution(vars));
if (tour.length < n) {
// Add subtour elimination constraint
GRBLinExpr expr = new GRBLinExpr();
for (int i = 0; i < tour.length; i++)
for (int j = i+1; j < tour.length; j++)
expr.addTerm(1.0, vars[tour[i]][tour[j]]);
addLazy(expr, GRB.LESS_EQUAL, tour.length-1);
}
}
} catch (GRBException e) {
System.out.println("Error code: " + e.getErrorCode() + ". " +
e.getMessage());
e.printStackTrace();
}
}
// Given an integer-feasible solution 'sol', return the smallest
// sub-tour (as a list of node indices).
protected static int[] findsubtour(double[][] sol)
{
int n = sol.length;
boolean[] seen = new boolean[n];
int[] tour = new int[n];
int bestind, bestlen;
int i, node, len, start;
for (i = 0; i < n; i++)
seen[i] = false;
start = 0;
bestlen = n+1;
bestind = -1;
node = 0;
while (start < n) {
for (node = 0; node < n; node++)
if (!seen[node])
break;
if (node == n)
break;
for (len = 0; len < n; len++) {
tour[start+len] = node;
seen[node] = true;
for (i = 0; i < n; i++) {
if (sol[node][i] > 0.5 && !seen[i]) {
node = i;
break;
}
}
if (i == n) {
len++;
if (len < bestlen) {
bestlen = len;
bestind = start;
}
start += len;
break;
}
}
}
int result[] = new int[bestlen];
for (i = 0; i < bestlen; i++)
result[i] = tour[bestind+i];
return result;
}
// Euclidean distance between points 'i' and 'j'
protected static double distance(double[] x,
double[] y,
int i,
int j) {
double dx = x[i]-x[j];
double dy = y[i]-y[j];
return Math.sqrt(dx*dx+dy*dy);
}
public static void main(String[] args) {
if (args.length < 1) {
System.out.println("Usage: java Tsp ncities");
System.exit(1);
}
int n = Integer.parseInt(args[0]);
try {
GRBEnv env = new GRBEnv();
GRBModel model = new GRBModel(env);
// Must set LazyConstraints parameter when using lazy constraints
model.set(GRB.IntParam.LazyConstraints, 1);
double[] x = new double[n];
double[] y = new double[n];
for (int i = 0; i < n; i++) {
x[i] = Math.random();
y[i] = Math.random();
}
// Create variables
GRBVar[][] vars = new GRBVar[n][n];
for (int i = 0; i < n; i++)
for (int j = 0; j <= i; j++) {
vars[i][j] = model.addVar(0.0, 1.0, distance(x, y, i, j),
GRB.BINARY,
"x"+String.valueOf(i)+"_"+String.valueOf(j));
vars[j][i] = vars[i][j];
}
// Degree-2 constraints
for (int i = 0; i < n; i++) {
GRBLinExpr expr = new GRBLinExpr();
for (int j = 0; j < n; j++)
expr.addTerm(1.0, vars[i][j]);
model.addConstr(expr, GRB.EQUAL, 2.0, "deg2_"+String.valueOf(i));
}
// Forbid edge from node back to itself
for (int i = 0; i < n; i++)
vars[i][i].set(GRB.DoubleAttr.UB, 0.0);
model.setCallback(new Tsp(vars));
model.optimize();
if (model.get(GRB.IntAttr.SolCount) > 0) {
int[] tour = findsubtour(model.get(GRB.DoubleAttr.X, vars));
assert tour.length == n;
System.out.print("Tour: ");
for (int i = 0; i < tour.length; i++)
System.out.print(String.valueOf(tour[i]) + " ");
System.out.println();
}
// Dispose of model and environment
model.dispose();
env.dispose();
} catch (GRBException e) {
System.out.println("Error code: " + e.getErrorCode() + ". " +
e.getMessage());
e.printStackTrace();
}
}
}
#!/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
constraints are added if needed."""
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)
return tour, m.ObjVal
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("")
' 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 node. The lazy constraint callback
' adds new constraints to cut them off.
Imports Gurobi
Class tsp_vb
Inherits GRBCallback
Private vars As GRBVar(,)
Public Sub New(xvars As GRBVar(,))
vars = xvars
End Sub
' Subtour elimination callback. Whenever a feasible solution is found,
' find the smallest subtour, and add a subtour elimination constraint
' if the tour doesn't visit every node.
Protected Overrides Sub Callback()
Try
If where = GRB.Callback.MIPSOL Then
' Found an integer feasible solution - does it visit every node?
Dim n As Integer = vars.GetLength(0)
Dim tour As Integer() = findsubtour(GetSolution(vars))
If tour.Length < n Then
' Add subtour elimination constraint
Dim expr As GRBLinExpr = 0
For i As Integer = 0 To tour.Length - 1
For j As Integer = i + 1 To tour.Length - 1
expr.AddTerm(1.0, vars(tour(i), tour(j)))
Next
Next
AddLazy(expr <= tour.Length - 1)
End If
End If
Catch e As GRBException
Console.WriteLine("Error code: " & e.ErrorCode & ". " & e.Message)
Console.WriteLine(e.StackTrace)
End Try
End Sub
' Given an integer-feasible solution 'sol', returns the smallest
' sub-tour (as a list of node indices).
Protected Shared Function findsubtour(sol As Double(,)) As Integer()
Dim n As Integer = sol.GetLength(0)
Dim seen As Boolean() = New Boolean(n - 1) {}
Dim tour As Integer() = New Integer(n - 1) {}
Dim bestind As Integer, bestlen As Integer
Dim i As Integer, node As Integer, len As Integer, start As Integer
For i = 0 To n - 1
seen(i) = False
Next
start = 0
bestlen = n+1
bestind = -1
node = 0
While start < n
For node = 0 To n - 1
if Not seen(node)
Exit For
End if
Next
if node = n
Exit While
End if
For len = 0 To n - 1
tour(start+len) = node
seen(node) = true
For i = 0 To n - 1
if sol(node, i) > 0.5 AndAlso Not seen(i)
node = i
Exit For
End If
Next
If i = n
len = len + 1
If len < bestlen
bestlen = len
bestind = start
End If
start = start + len
Exit For
End If
Next
End While
For i = 0 To bestlen - 1
tour(i) = tour(bestind+i)
Next
System.Array.Resize(tour, bestlen)
Return tour
End Function
' Euclidean distance between points 'i' and 'j'
Protected Shared Function distance(x As Double(), y As Double(), _
i As Integer, j As Integer) As Double
Dim dx As Double = x(i) - x(j)
Dim dy As Double = y(i) - y(j)
Return Math.Sqrt(dx * dx + dy * dy)
End Function
Public Shared Sub Main(args As String())
If args.Length < 1 Then
Console.WriteLine("Usage: tsp_vb nnodes")
Return
End If
Dim n As Integer = Convert.ToInt32(args(0))
Try
Dim env As New GRBEnv()
Dim model As New GRBModel(env)
' Must set LazyConstraints parameter when using lazy constraints
model.Parameters.LazyConstraints = 1
Dim x As Double() = New Double(n - 1) {}
Dim y As Double() = New Double(n - 1) {}
Dim r As New Random()
For i As Integer = 0 To n - 1
x(i) = r.NextDouble()
y(i) = r.NextDouble()
Next
' Create variables
Dim vars As GRBVar(,) = New GRBVar(n - 1, n - 1) {}
For i As Integer = 0 To n - 1
For j As Integer = 0 To i
vars(i, j) = model.AddVar(0.0, 1.0, distance(x, y, i, j), _
GRB.BINARY, "x" & i & "_" & j)
vars(j, i) = vars(i, j)
Next
Next
' Degree-2 constraints
For i As Integer = 0 To n - 1
Dim expr As GRBLinExpr = 0
For j As Integer = 0 To n - 1
expr.AddTerm(1.0, vars(i, j))
Next
model.AddConstr(expr = 2.0, "deg2_" & i)
Next
' Forbid edge from node back to itself
For i As Integer = 0 To n - 1
vars(i, i).UB = 0.0
Next
model.SetCallback(New tsp_vb(vars))
model.Optimize()
If model.SolCount > 0 Then
Dim tour As Integer() = findsubtour(model.Get(GRB.DoubleAttr.X, vars))
Console.Write("Tour: ")
For i As Integer = 0 To tour.Length - 1
Console.Write(tour(i) & " ")
Next
Console.WriteLine()
End If
' Dispose of model and environment
model.Dispose()
env.Dispose()
Catch e As GRBException
Console.WriteLine("Error code: " & e.ErrorCode & ". " & e.Message)
Console.WriteLine(e.StackTrace)
End Try
End Sub
End Class