/*********************************************************************** Xpress Optimizer Examples ========================= file els_managedcuts.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; 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; 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; 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; 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, "els_managedcuts", "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; }