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

# tsp_c.c

    /* 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;
    }

Help and Feedback

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