Hide navigation sidebar
Hide table of contents sidebar
Back to top
Toggle table of contents sidebar

# 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 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 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 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 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 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 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 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 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

    # 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.

    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 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 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

Help and Feedback

Press `/` or click to search this documentation | Gurobi website | Gurobi community forum
Go to other Gurobi documentation | Gurobi download | Contact support