Cut Generation for an economic lot-sizing (ELS) problem

The two programs demonstrate two methods for adding cutting planes using a cut round callback. They both formulate and solve an economic lot-sizing problem and adds valid (l-S)-inequalities as cutting planes within a cut round callback.

els_addcuts.c uses XPRSaddcuts to directly add the new cutting planes to a node problem of the branch-and-bound search. The cutting planes must be presolved first before they can be added.

els_managedcuts.c uses XPRSmanagedcuts to provide cutting planes to the Optimizer, which will managed these cuts automatically. The cutting planes can be provided in the original space and do not need to be presolved.

/***********************************************************************
  Xpress Optimizer Examples
  =========================

  file els.c
  ``````````
  Demonstrates how to implement cutting planes as part of a MIP branch-and=
  bound search using the cutround callback.

  Cuts are added as Optimizer managed cuts using XPRSaddmanagedcuts.

  Economic lot sizing problem. Solved by adding (l,S)-inequalities
  in a branch-and-cut heuristic (using the cutround callback).

  Problem description:
  ````````````````````

  ELS considers production planning over a horizon of T periods. In period t,
  t=1,...,T, there is a given demand DEMAND[p,t] that must be satisfied by
  production produce[p,t] in period t and by inventory carried over from previous
  periods. There is a set-up cost SETUPCOST[t] associated with production in
  period t and the total production capacity per period is limited. The unit
  production cost in period t is PRODCOST[p,t]. There is no
  inventory or stock-holding cost.

  Cutting planes:
  ```````````````

  A well-known class of valid inequalities for ELS are the
  (l,S)-inequalities.  Let D(p,q) denote the demand in periods p
  to q and y(t) be a binary variable indicating whether there is any
  production in period t.  For each period l and each subset of periods S
  of 1 to l, the (l,S)-inequality is

      sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(t,l) * y(t)
          >= D(1,l)

  It says that actual production x(t) in periods included S plus maximum
  potential production D(t,l)*y(t) in the remaining periods (those not in
  S) must at least equal total demand in periods 1 to l.  Note that in
  period t at most D(t,l) production is required to meet demand up to
  period l.

  Based on the observation that

      sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(t,l) * y(t)
          >= sum (t=1:l) min(x(t), D(t,l) * y(t))
          >= D(1,l)

  it is easy to develop a separation algorithm and thus a cutting plane
  algorithm based on these (l,S)-inequalities.

  (c) 2025 Fair Isaac Corporation

*****************************************************************/

#include <stdlib.h>
#include <stdio.h>

#include "xprs.h"

/*
 * Calls an Xpress optimizer function and checks the return code.
 * If the call fails then the function
 * - prints a short error message to stderr,
 * - sets variable 'returnCode' to the error,
 * - and branches to label 'cleanup'.
 */
#define CHECK_RETURN(call) do {                         \
    int result_ = call;                                 \
    if ( result_ != 0 ) {                               \
      fprintf(stderr, "Line %d: %s failed with %d\n",   \
              __LINE__, #call, result_);                \
      returnCode = result_;                             \
      goto cleanup;                                     \
    }                                                   \
  } while (0)

static int callbackError = 0;       /* Indicates an error in a callback. */

/* XPRS optimizer message callback */
static void XPRS_CC cbMessage(XPRSprob cbprob, void *cbdata, const char *msg, int len, int msgtype)
{
  (void)cbprob; /* unused */
  (void)cbdata; /* unused */
  if (msgtype > 0) printf("%*s\n", len, msg);
}

/* Tolerance for satisfiability */
#define EPS 1e-6

/* Number of time periods. Reduce to create an easier problem. */
#define DIM_TIMES 15
/* Number of products to produce. */
#define DIM_PRODUCTS 4

static int DEMAND[DIM_PRODUCTS][DIM_TIMES] = {
  {2, 3, 5, 3, 4, 2, 5, 4, 1, 3, 4, 2, 3, 5, 2},
  {3, 1, 2, 3, 5, 3, 1, 2, 3, 3, 4, 5, 1, 4, 1},
  {3, 5, 2, 1, 2, 1, 3, 3, 5, 2, 2, 1, 3, 2, 3},
  {2, 2, 1, 3, 2, 1, 2, 2, 3, 3, 2, 2, 3, 1, 2}
};

static int SETUPCOST[DIM_TIMES] = {
  17, 14, 11, 6, 9, 6, 15, 10, 8, 7, 12, 9, 10, 8, 12
};

static int PRODCOST[DIM_PRODUCTS][DIM_TIMES] = {
  {5, 3, 2, 1, 3, 1, 4, 3, 2, 2, 3, 1, 2, 3, 2},
  {1, 4, 2, 3, 1, 3, 1, 2, 3, 3, 3, 4, 4, 2, 2},
  {3, 3, 3, 4, 4, 3, 3, 3, 2, 2, 1, 1, 3, 3, 3},
  {2, 2, 2, 3, 3, 3, 4, 4, 4, 3, 3, 2, 2, 2, 3}
};

static int CAP[DIM_TIMES] = {
  12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12
};

struct prob_data_s {
  int produce[DIM_PRODUCTS][DIM_TIMES];
  int setup[DIM_PRODUCTS][DIM_TIMES];
  double D[DIM_PRODUCTS][DIM_TIMES][DIM_TIMES];
};

/* Creates an economic lot sizing problem. */
static int LoadELSProblem(XPRSprob prob, struct prob_data_s *probData)
{
  int returnCode = 0;

  int p, s, t;
  int col_next;

  int xprs_cols;

  int *row_inds = NULL;
  int row_start = 0;
  double *row_coefs = NULL;
  double row_rhs;
  char row_type;
  int row_length;

  /* Calculate demand D(p, s, t) as the demand for product p from time s to time t(inclusive). */
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (s = 0; s < DIM_TIMES; s++) {
      double sum_demand = 0.0;
      for (t = s; t < DIM_TIMES; t++) {
        sum_demand += DEMAND[p][t];
        probData->D[p][s][t] = sum_demand;
      }
    }
  }

  /* Start with an empty problem. */
  CHECK_RETURN(XPRSloadlp(prob, "ELS", 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));

  /* Define the variables. */
  CHECK_RETURN(XPRSgetintattrib(prob, XPRS_COLS, &xprs_cols));
  col_next = xprs_cols;
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (t = 0; t < DIM_TIMES; t++) {
      probData->produce[p][t] = col_next++;
    }
  }
  CHECK_RETURN(XPRSaddcols(prob, col_next - xprs_cols, 0, NULL, NULL, NULL, NULL, NULL, NULL));
  CHECK_RETURN(XPRSgetintattrib(prob, XPRS_COLS, &xprs_cols));
  col_next = xprs_cols;
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (t = 0; t < DIM_TIMES; t++) {
      probData->setup[p][t] = col_next++;
    }
  }
  CHECK_RETURN(XPRSaddcols(prob, col_next - xprs_cols, 0, NULL, NULL, NULL, NULL, NULL, NULL));

  /* Allocate space to define constraints. */
  CHECK_RETURN(XPRSgetintattrib(prob, XPRS_COLS, &xprs_cols));
  row_coefs = malloc(xprs_cols * sizeof(double));
  row_inds = malloc(xprs_cols * sizeof(int));
  if ((row_coefs == NULL) || (row_inds == NULL)) {
    fprintf(stderr, "malloc failure");
    returnCode = -1;
    goto cleanup;
  }

  /* Add the objective function :
   *  MinCost:= sum(t in TIMES) (SETUPCOST(t) * sum(p in PRODUCTS) setup(p,t) +
   *                            sum(p in PRODUCTS) PRODCOST(p, t) *produce(p, t) )
   */
  row_length = 0;
  for (t = 0; t < DIM_TIMES; t++) {
    for (p = 0; p < DIM_PRODUCTS; p++) {
      row_inds[row_length] = probData->setup[p][t]; row_coefs[row_length] = SETUPCOST[t]; row_length++;
      row_inds[row_length] = probData->produce[p][t]; row_coefs[row_length] = PRODCOST[p][t]; row_length++;
    }
  }
  CHECK_RETURN(XPRSchgobj(prob, row_length, row_inds, row_coefs));

  /* Add constraints. */

  /* Production in periods 0 to t must satisfy the total demand
   * during this period of time, for all t in [0,T[
   *    forall(p in PRODUCTS, t in TIMES)
   *      Dem(t) : = sum(s in 1..t) produce(p,s) >= sum(s in 1..t) DEMAND(s)
   */
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (t = 0; t < DIM_TIMES; t++) {
      row_length = 0;
      for (s = 0; s <= t; s++) {
        row_inds[row_length] = probData->produce[p][s]; row_coefs[row_length] = 1.0; row_length++;
      }
      row_rhs = 0.0;
      for (s = 0; s <= t; s++) row_rhs += DEMAND[p][s];
      row_type = 'G';
      CHECK_RETURN(XPRSaddrows(prob, 1, row_length, &row_type, &row_rhs, NULL, &row_start, row_inds, row_coefs));
    }
  }

  /* If there is production during t then there is a setup in t :
   *    forall(p in PRODUCTS, t in TIMES)
   *       ProdSetup(t) : = produce(t) <= D(t,TIMES) * setup(t)
   */
  for (p = 0; p < DIM_PRODUCTS; p++) {
    double sum_demand = 0.0;
    for (t = DIM_TIMES - 1; t >= 0; t--) {
      sum_demand += DEMAND[p][t];
      row_length = 0;
      row_inds[row_length] = probData->produce[p][t]; row_coefs[row_length] = 1.0; row_length++;
      row_inds[row_length] = probData->setup[p][t]; row_coefs[row_length] = -sum_demand; row_length++;
      row_rhs = 0.0;
      row_type = 'L';
      CHECK_RETURN(XPRSaddrows(prob, 1, row_length, &row_type, &row_rhs, NULL, &row_start, row_inds, row_coefs));
    }
  }

  /* Capacity limits :
   *    forall(t in TIMES) Capacity(t) : = sum(p in PRODUCTS) produce(p, t) <= CAP(t)
   */
  for (t = 0; t < DIM_TIMES; t++) {
    row_length = 0;
    for (p = 0; p < DIM_PRODUCTS; p++) {
      row_inds[row_length] = probData->produce[p][t]; row_coefs[row_length] = 1.0; row_length++;
    }
    row_rhs = CAP[t];
    row_type = 'L';
    CHECK_RETURN(XPRSaddrows(prob, 1, row_length, &row_type, &row_rhs, NULL, &row_start, row_inds, row_coefs));
  }

  /* setup variables are binary:
   *   forall(p in PRODUCTS, t in TIMES) setup(p,t) is_binary
   */
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (t = 0; t < DIM_TIMES; t++) {
      char col_type = 'B';
      CHECK_RETURN(XPRSchgcoltype(prob, 1, &probData->setup[p][t], &col_type));
    }
  }

  /* Add names for the variables. */
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (t = 0; t < DIM_TIMES; t++) {
      char name[32];
      snprintf(name, sizeof(name), "setup(%i,%i)", p, t);
      CHECK_RETURN(XPRSaddnames(prob, XPRS_NAMES_COLUMN, name, probData->setup[p][t], probData->setup[p][t]));
      snprintf(name, sizeof(name), "produce(%i,%i)", p, t);
      CHECK_RETURN(XPRSaddnames(prob, XPRS_NAMES_COLUMN, name, probData->produce[p][t], probData->produce[p][t]));
    }
  }

cleanup:

  if (row_inds) free(row_inds);
  if (row_coefs) free(row_coefs);

  return returnCode;
}

/* Cut round callback for separating our cutting planes. */
void XPRS_CC cbCutRound(XPRSprob prob, void *cbdata, int ifxpresscuts, int *p_action)
{
  (void)ifxpresscuts; /* unused */
  (void)p_action; /* unused */

  int returnCode = 0;

  struct prob_data_s *probData = (struct prob_data_s *)cbdata;

  int p, t, l;

  int xprs_inputcols;
  int xprs_cutrounds;

  double *sol = NULL;
  int isAvailable;

  char    cutSense;
  int     cutLen;
  int    *cutInds = NULL;
  double *cutCoefs = NULL;
  double  cutRhs;
  int cutsStart[2];
  int nCutsAdded = 0;

  /* Apply only one round of cutting on each node. */
  /* Because the CutRound callback is fired *before* a round of cutting, */
  /* the CUTROUNDS attribute will start from 0 on the first invocation of */
  /* the callback. */
  CHECK_RETURN(XPRSgetintattrib(prob, XPRS_CUTROUNDS, &xprs_cutrounds));
  if (xprs_cutrounds >= 1) {
    goto cleanup;
  }

  /* We need the dimension of the solution vector. */
  CHECK_RETURN(XPRSgetintattrib(prob, XPRS_INPUTCOLS, &xprs_inputcols));

  /* Get the solution vector. */
  sol = malloc(xprs_inputcols * sizeof(*sol));
  if (!sol) {
    fprintf(stderr, "malloc failure");
    returnCode = -1;
    goto cleanup;
  }
  CHECK_RETURN(XPRSgetcallbacksolution(prob, &isAvailable, sol, 0, xprs_inputcols-1));
  /* Xpress will only fire the CutRound callback when a solution is available, */
  /* so there is no need to check the isAvailable return argument. */

  /* Search for violated constraints : the minimum of the actual production
   * produce[p][t] and the maximum potential production D[p][t][l]*setup[p][t]
   * in periods 0 to l must at least equal the total demand in periods
   * 0 to l.
   *    sum(t=1:l) min(prod[p][t], D[p][t][l]*setup[p][t]) >= D[p][0][l]
   */
  for (p = 0; p < DIM_PRODUCTS; p++) {
    for (l = 0; l < DIM_TIMES; l++) {
      double sum = 0.0;
      for (t = 0; t <= l; t++) {
        if (sol[probData->produce[p][t]] < probData->D[p][t][l] * sol[probData->setup[p][t]] + EPS) {
          sum += sol[probData->produce[p][t]];
        } else {
          sum += probData->D[p][t][l] * sol[probData->setup[p][t]];
        }
      }
      if (sum < probData->D[p][0][l] - EPS) {
        // Create the violated cut.
        if (!cutInds || !cutCoefs) {
          cutInds = malloc(xprs_inputcols * sizeof(*cutInds));
          cutCoefs = malloc(xprs_inputcols * sizeof(*cutCoefs));
          if (!cutInds || !cutCoefs) {
            fprintf(stderr, "malloc failure");
            returnCode = -1;
            goto cleanup;
          }
        }
        cutLen = 0;
        cutSense = 'G';
        for (t = 0; t <= l; t++) {
          if (sol[probData->produce[p][t]] < probData->D[p][t][l] * sol[probData->setup[p][t]] + EPS) {
            cutCoefs[cutLen] = 1.0;
            cutInds[cutLen] = probData->produce[p][t];
            cutLen += 1;
          } else {
            cutCoefs[cutLen] = probData->D[p][t][l];
            cutInds[cutLen] = probData->setup[p][t];
            cutLen += 1;
          }
        }
        cutRhs = probData->D[p][0][l];

        /* Give the cut to Xpress to manage. It will automatically be presolved. */
        cutsStart[0] = 0;
        cutsStart[1] = cutLen;
        CHECK_RETURN(XPRSaddmanagedcuts(prob, 1 /*globalvalid*/, 1, &cutSense, &cutRhs, cutsStart, cutInds, cutCoefs));
        nCutsAdded += 1;
        /* If we add any cuts in the callback, Xpress will automatically trigger another roound of cuts, */
        /* so there is no need to set the p_action return argument. */
      }
    }
  }

  if (nCutsAdded > 0) {
    printf("Cuts added: %i\n", nCutsAdded);
  }

cleanup:

  if (sol) free(sol);
  if (cutInds) free(cutInds);
  if (cutCoefs) free(cutCoefs);

  if (returnCode) {
    /* There was an error in the callback.Stop the Optimizer and remember the error. */
    XPRSinterrupt(prob, XPRS_STOP_USER);
    callbackError = 1;
  }
}

int main(void) {

  int returnCode = 0;

  XPRSprob prob = NULL;

  int solvestatus, solstatus;

  struct prob_data_s probData = { 0 };

  /* Initialize the Optimizer. */
  if (XPRSinit("") != 0) {
    char message[512];
    XPRSgetlicerrmsg(message, sizeof(message));
    fprintf(stderr, "Licensing error: %s\n", message);
    return 1;
  }

  /* Create a problem with a message handler so we can see what is going on. */
  CHECK_RETURN(XPRScreateprob(&prob));
  CHECK_RETURN(XPRSaddcbmessage(prob, cbMessage, NULL, 0));

  /* Create our Economic Lot Sizing problem. */
  CHECK_RETURN(LoadELSProblem(prob, &probData));
  CHECK_RETURN(XPRSwriteprob(prob, "problem", "l"));

  /* Add a CutRound callback for separating our cuts. */
  CHECK_RETURN(XPRSaddcbcutround(prob, cbCutRound, &probData, 0));

  /* Solve the ELS problem. */
  CHECK_RETURN(XPRSoptimize(prob, "", &solvestatus, &solstatus));
  if (callbackError) {
    returnCode = -3;
    goto cleanup;
  }

  if ((solvestatus == XPRS_SOLVESTATUS_COMPLETED) && (solstatus == XPRS_SOLSTATUS_OPTIMAL)) {
    printf("Solved problem to optimality.\n");
  } else {
    printf("Failed to solve problem with solvestatus %i and solstatus %i\n", solvestatus, solstatus);
  }

cleanup:

  XPRSdestroyprob(prob);
  XPRSfree();

  return returnCode;
}