Tsp.java#

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