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