// (c) 2024-2024 Fair Isaac Corporation
/**
* Economic lot sizing, ELS, problem. Solved by adding (l,S)-inequalities in
* several rounds looping over the root node.
*
* ELS considers production planning over a horizon of T periods. In period t,
* t=1,...,T, there is a given demand DEMAND[t] that must be satisfied by
* production prod[t] in period t and by inventory carried over from previous
* periods. There is a set-up up cost SETUPCOST[t] associated with production in
* period t. The unit production cost in period t is PRODCOST[t]. There is no
* inventory or stock-holding cost.
*/
#include <chrono>
#include <iostream>
#include <xpress.hpp>
using namespace xpress;
using namespace xpress::objects;
using xpress::objects::utils::scalarProduct;
using xpress::objects::utils::sum;
using std::chrono::duration_cast;
using std::chrono::steady_clock;
double const EPS = 1e-6;
int const T = 6; /* Number of time periods */
/* Data */
std::vector<double> DEMAND{1, 3, 5, 3, 4, 2}; /* Demand per period */
std::vector<double> SETUPCOST{17, 16, 11, 6, 9, 6}; /* Setup cost / period */
std::vector<double> PRODCOST{5, 3, 2, 1, 3, 1}; /* Prod. cost / period */
std::vector<std::vector<double>>
D(T, std::vector<double>(T)); /* Total demand in periods t1 - t2 */
/** Create the model in P.
* @param p The Xpress Solver instance in which the model is created.
* @param prod Array to store the production variables.
* @param setup Array to store the setup variables.
*/
void modEls(XpressProblem &p, std::vector<Variable> &prod,
std::vector<Variable> &setup) {
for (int s = 0; s < T; s++)
for (int t = 0; t < T; t++)
for (int k = s; k <= t; k++)
D[s][t] += DEMAND[k];
// Variables
prod = p.addVariables(T)
.withType(ColumnType::Continuous)
.withName([](auto t) { return "prod" + std::to_string(t + 1); })
.toArray();
setup = p.addVariables(T)
.withType(ColumnType::Binary)
.withName([](auto t) { return "setup" + std::to_string(t + 1); })
.toArray();
// Objective: Minimize total cost
p.setObjective(
sum(scalarProduct(setup, SETUPCOST), scalarProduct(prod, PRODCOST)),
ObjSense::Minimize);
// Constraints
// Production in period t must not exceed the total demand for the
// remaining periods; if there is production during t then there
// is a setup in t
// for all t in [0,T[
// prod[t] <= setup[t] * D[t][T-1]
p.addConstraints(T,
[&](auto t) { return prod[t] <= setup[t] * D[t][T - 1]; });
// Production in periods 0 to t must satisfy the total demand
// during this period of time
// for all t in [0,T[
// sum(s in [0,t+1[) prod[s] >= D[0][t]
p.addConstraints(T, [&](auto t) {
return sum(t + 1, [&](auto s) { return prod[s]; }) >= D[0][t];
});
p.writeProb("ELS.lp", "l");
}
/**
* Cut generation loop at the top node:
* solve the LP and save the basis
* get the solution values
* identify and set up violated constraints
* load the modified problem and load the saved basis
* @param p The Xpress Solver instance that is populated with the model.
* @param prod Production variables.
* @param setup Setup variables.
*/
void solveEls(XpressProblem &p, std::vector<Variable> &prod,
std::vector<Variable> &setup) {
p.callbacks.addMessageCallback(XpressProblem::console);
/* Disable automatic cuts - we use our own */
p.controls.setCutStrategy(XPRS_CUTSTRATEGY_NONE);
/* Switch presolve off */
p.controls.setPresolve(XPRS_PRESOLVE_NONE);
int ncut = 0, npass = 0, npcut = 0;
steady_clock::time_point starttime = steady_clock::now();
std::vector<double> sol;
do {
p.writeProb("model" + std::to_string(npass) + ".lp");
npass++;
npcut = 0;
// Solve the LP-problem
p.lpOptimize();
if (p.attributes.getSolStatus() != SolStatus::Optimal)
throw std::runtime_error("failed to optimize with status " +
to_string(p.attributes.getSolStatus()));
// Get the solution values:
sol = p.getSolution();
// Search for violated constraints:
for (int l = 0; l < T; l++) {
double ds = 0.0;
for (int t = 0; t <= l; t++) {
if (prod[t].getValue(sol) < D[t][l] * setup[t].getValue(sol) + EPS) {
ds += prod[t].getValue(sol);
} else {
ds += D[t][l] * setup[t].getValue(sol);
}
}
// Add the violated inequality: the minimum of the actual production
// prod[t] and the maximum potential production D[t][l]*setup[t] in
// periods 0 to l must at least equal the total demand in periods 0 to l.
// sum(t=1:l) min(prod[t], D[t][l]*setup[t]) >= D[0][l]
if (ds < D[0][l] - EPS) {
Expression cut = sum(l + 1, [&](auto t) {
if (prod[t].getValue(sol) < D[t][l] * setup[t].getValue(sol) + EPS)
return 1.0 * prod[t];
else
return D[t][l] * setup[t];
});
p.addConstraint(cut >= D[0][l]);
ncut++;
npcut++;
}
}
steady_clock::time_point endtime = steady_clock::now();
std::cout << "Iteration " << npass << ", "
<< duration_cast<std::chrono::milliseconds>(endtime - starttime)
.count() *
1e-3
<< " sec"
<< ", objective value: " << p.attributes.getObjVal()
<< ", cuts added: " << npcut << " (total " << ncut << ")"
<< std::endl;
if (npcut == 0)
std::cout << "Optimal integer solution found:" << std::endl;
} while (npcut > 0);
// Print out the solution:
for (int t = 0; t < T; t++) {
std::cout << "Period " << (t + 1) << ": prod" << prod[t].getValue(sol)
<< " (demand: " << DEMAND[t] << ", cost: " << PRODCOST[t] << ")"
<< ", setup " << setup[t].getValue(sol) << " (cost "
<< SETUPCOST[t] << ")" << std::endl;
}
}
int main() {
XpressProblem prob;
std::vector<Variable> prod;
std::vector<Variable> setup;
modEls(prob, prod, setup); // Model the problem
solveEls(prob, prod, setup); // Solve the problem
return 0;
}
|
// (c) 2024-2024 Fair Isaac Corporation
/** Economic lot sizing, ELS, problem.
* Solved by adding (l,S)-inequalities in a branch-and-cut heuristic
* (using the cut manager).
*
* ELS considers production planning over a horizon of T periods. In period t,
* t=1,...,T, there is a given demand DEMAND[t] that must be satisfied by
* production prod[t] in period t and by inventory carried over from previous
* periods. There is a set-up up cost SETUPCOST[t] associated with production in
* period t. The unit production cost in period t is PRODCOST[t]. There is no
* inventory or stock-holding cost.
*
* *** This model cannot be run with a Community Licence ***
*/
#include <iostream>
#include <xpress.hpp>
using namespace xpress;
using namespace xpress::objects;
using xpress::objects::utils::scalarProduct;
using xpress::objects::utils::sum;
double const EPS = 1e-6;
int const T = 6; /* Number of time periods */
/* Data */
std::vector<double> DEMAND{1, 3, 5, 3, 4, 2}; /* Demand per period */
std::vector<double> SETUPCOST{17, 16, 11, 6, 9, 6}; /* Setup cost / period */
std::vector<double> PRODCOST{5, 3, 2, 1, 3, 1}; /* Prod. cost / period */
std::vector<std::vector<double>>
D(T, std::vector<double>(T)); /* Total demand in periods t1 - t2 */
void printProblemStatus(XpressProblem const &prob) {
std::cout << "Problem status:" << std::endl
<< "\tSolve status: " << prob.attributes.getSolveStatus()
<< std::endl
<< "\tLP status: " << prob.attributes.getLpStatus() << std::endl
<< "\tMIP status: " << prob.attributes.getMipStatus() << std::endl
<< "\tSol status: " << prob.attributes.getSolStatus() << std::endl;
}
/**
* Cut generation algorithm:
* - get the solution values
* - identify and set up violated constraints
* - add cuts to the problem
* @param p Xpress solver instance as passed to the callback.
* @param prod Production variables.
* @param setup Setup variables.
*/
void separate(XpressProblem &p, std::vector<Variable> &prod,
std::vector<Variable> &setup) {
int ncut = 0;
// Add cut only to optimal relaxations
if (p.attributes.getLpStatus() != LPStatus::Optimal)
return;
auto sol = p.getCallbackSolution();
auto slack = p.getCallbackSlacks();
auto duals = p.getCallbackDuals();
auto djs = p.getCallbackRedCosts();
// Search for violated constraints:
for (int l = 0; l < T; l++) {
double ds = 0.0;
for (int t = 0; t <= l; t++) {
if (prod[t].getValue(sol) < D[t][l] * setup[t].getValue(sol) + EPS) {
ds += prod[t].getValue(sol);
} else {
ds += D[t][l] * setup[t].getValue(sol);
}
}
// Add the violated inequality: the minimum of the actual production
// prod[t] and the maximum potential production D[t][l]*setup[t]
// in periods 0 to l must at least equal the total demand in periods
// 0 to l.
// sum(t=1:l) min(prod[t], D[t][l]*setup[t]) >= D[0][l] */
if (ds < D[0][l] - EPS) {
Expression cut = sum(l + 1, [&](auto t) {
if (prod[t].getValue(sol) < D[t][l] * setup[t].getValue(sol) + EPS)
return 1.0 * prod[t];
else
return D[t][l] * setup[t];
});
p.addCut(0, cut >= D[0][1]);
ncut++;
}
}
if (ncut > 0) {
std::cout << "Cuts added: " << ncut << "(depth "
<< p.attributes.getNodeDepth() << ", node "
<< p.attributes.getNodes() << ")" << std::endl;
}
}
/** Create the model in P.
* @param p The Xpress Solver instance in which the model is created.
* @param prod Array to store the production variables.
* @param setup Array to store the setup variables.
*/
void modEls(XpressProblem &p, std::vector<Variable> &prod,
std::vector<Variable> &setup) {
for (int s = 0; s < T; s++)
for (int t = 0; t < T; t++)
for (int k = s; k <= t; k++)
D[s][t] += DEMAND[k];
// Variables
prod = p.addVariables(T)
.withType(ColumnType::Continuous)
.withName([](auto t) { return "prod" + std::to_string(t + 1); })
.toArray();
setup = p.addVariables(T)
.withType(ColumnType::Binary)
.withName([](auto t) { return "setup" + std::to_string(t + 1); })
.toArray();
// Objective: Minimize total cost
p.setObjective(
sum(scalarProduct(setup, SETUPCOST), scalarProduct(prod, PRODCOST)),
ObjSense::Minimize);
// Constraints
// Production in period t must not exceed the total demand for the
// remaining periods; if there is production during t then there
// is a setup in t
// for all t in [0,T[
// prod[t] <= setup[t] * D[t][T-1]
p.addConstraints(T,
[&](auto t) { return prod[t] <= setup[t] * D[t][T - 1]; });
// Production in periods 0 to t must satisfy the total demand
// during this period of time
// for all t in [0,T[
// sum(s in [0,t+1[) prod[s] >= D[0][t]
p.addConstraints(T, [&](auto t) {
return sum(t + 1, [&](auto s) { return prod[s]; }) >= D[0][t];
});
p.writeProb("ELSCut.lp", "l");
}
/** Solve the model.
* Cuts are added dynamically with the optnode callback.
* @param p The Xpress Solver instance that is populated with the model.
* @param prod Production variables.
* @param setup Setup variables.
*/
void solveEls(XpressProblem &p, std::vector<Variable> &prod,
std::vector<Variable> &setup) {
p.callbacks.addMessageCallback(XpressProblem::console);
p.controls.setLpLog(0);
p.controls.setMipLog(3);
// Disable automatic cuts - we use our own
p.controls.setCutStrategy(XPRS_CUTSTRATEGY_NONE);
// Switch presolve off
p.controls.setPresolve(XPRS_PRESOLVE_NONE);
p.controls.setMipPresolve(0);
p.callbacks.addOptnodeCallback(
[&](auto &prob, auto) { separate(prob, prod, setup); });
// Solve the MIP
p.optimize();
if (p.attributes.getSolStatus() != SolStatus::Optimal)
throw std::runtime_error("optimization failed with status " +
to_string(p.attributes.getSolStatus()));
// Get the solution values:
auto sol = p.getSolution();
// Print out the solution:
for (int t = 0; t < T; t++) {
std::cout << "Period " << (t + 1) << ": prod" << prod[t].getValue(sol)
<< " (demand: " << DEMAND[t] << ", cost: " << PRODCOST[t] << ")"
<< ", setup " << setup[t].getValue(sol) << " (cost "
<< SETUPCOST[t] << ")" << std::endl;
}
printProblemStatus(p);
}
int main() {
XpressProblem prob;
std::vector<Variable> prod;
std::vector<Variable> setup;
modEls(prob, prod, setup); // Model the problem
solveEls(prob, prod, setup); // Solve the problem
return 0;
}
|
// (c) 2023-2025 Fair Isaac Corporation
/** Demonstrates how to implement cutting planes as part of a MIP
* branch-and-bound search using the cutround callback.
*
* Cuts are added as user cuts using XpressProblem.addManagedCut().
*
* Economic lot sizing problem. Solved by adding (l,S)-inequalities
* in a branch-and-cut heuristic (using the cutround callback).
*
* 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.
*
* 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
* <pre>
* 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)
* </pre>
*
* 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
* <pre>
* 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)
* </pre>
* it is easy to develop a separation algorithm and thus a cutting plane
* algorithm based on these (l,S)-inequalities.
*/
#include <xpress.hpp>
#include <xpress_objects.hpp>
#include <iostream>
using namespace xpress;
using namespace xpress::objects;
using namespace xpress::objects::utils;
/** Tolerance for satisfiability. */
#define EPS 1e-6
/** Number of time periods. */
#define DIM_TIMES 15
/** Number of products to produce. */
#define DIM_PRODUCTS 4
/** Demand per product (first dim) and time period (second dim). */
static int const 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}
};
/** Setup cost. */
static double const SETUPCOST[DIM_TIMES] = {
17, 14, 11, 6, 9, 6, 15, 10, 8, 7, 12, 9, 10, 8, 12
};
/** Production cost per product (first dim) and time period (second dim). */
static double const 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}
};
/** Capacity. */
static const int CAP[DIM_TIMES] = {
12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12
};
/** Cut round callback for separating our cutting planes. */
class Callback {
std::vector<std::vector<Variable>> const &produce;
std::vector<std::vector<Variable>> const &setup;
double const (&sumDemand)[DIM_PRODUCTS][DIM_TIMES][DIM_TIMES];
public:
Callback(std::vector<std::vector<Variable>> const &p,
std::vector<std::vector<Variable>> const &s,
double const (&d)[DIM_PRODUCTS][DIM_TIMES][DIM_TIMES])
: produce(p)
, setup(s)
, sumDemand(d)
{
}
void operator()(XpressProblem &prob, int, int *) const {
// 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.
if (prob.attributes.getCutRounds() >= 1)
return;
// Get the solution vector.
// Xpress will only fire the CutRound callback when a solution is
// available, so there is no need to check whether a solution is
// available.
auto sol = prob.getCallbackSolution();
// 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]
int nCutsAdded = 0;
for (int p = 0; p < DIM_PRODUCTS; p++) {
for (int l = 0; l < DIM_TIMES; l++) {
double sum = 0.0;
for (int t = 0; t <= l; t++) {
if (produce[p][t].getValue(sol) < sumDemand[p][t][l] * setup[p][t].getValue(sol) + EPS) {
sum += produce[p][t].getValue(sol);
}
else {
sum += sumDemand[p][t][l] * setup[p][t].getValue(sol);
}
}
if (sum < sumDemand[p][0][l] - EPS) {
// Create the violated cut.
LinExpression cut = LinExpression::create();
for (int t = 0; t <= l; t++) {
if (produce[p][t].getValue(sol) < sumDemand[p][t][l] * setup[p][t].getValue(sol)) {
cut.addTerm(1.0, produce[p][t]);
}
else {
cut.addTerm(sumDemand[p][t][l], setup[p][t]);
}
}
// Give the cut to Xpress to manage.
// It will automatically be presolved.
prob.addManagedCut(true, cut >= sumDemand[p][0][l]);
++nCutsAdded;
// If we modified the problem in the callback, Xpress
// will automatically trigger another roound of cuts,
// so there is no need to set the action return
// argument.
}
}
}
if (nCutsAdded > 0) {
std::cout << "Cuts added: " << nCutsAdded << std::endl;
}
}
};
int main(void) {
XpressProblem prob;
prob.callbacks.addMessageCallback(XpressProblem::console);
// Create an economic lot sizing problem.
// Calculate demand D(p, s, t) as the demand for product p from
// time s to time t(inclusive).
double sumDemand[DIM_PRODUCTS][DIM_TIMES][DIM_TIMES];
for (int p = 0; p < DIM_PRODUCTS; p++) {
for (int s = 0; s < DIM_TIMES; s++) {
double thisSumDemand = 0.0;
for (int t = s; t < DIM_TIMES; t++) {
thisSumDemand += DEMAND[p][t];
sumDemand[p][s][t] = thisSumDemand;
}
}
}
auto produce = prob.addVariables(DIM_PRODUCTS, DIM_TIMES)
.withName([](auto p, auto t){ return xpress::format("produce(%d,%d)", p, t); })
.toArray();
auto setup = prob.addVariables(DIM_PRODUCTS, DIM_TIMES)
.withType(ColumnType::Binary)
.withName([](auto p,auto t){ return xpress::format("setup(%d,%d)", p, t); })
.toArray();
// 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) )
prob.setObjective(sum(DIM_TIMES, DIM_PRODUCTS,
[&](auto t, auto p){
return SETUPCOST[t] * setup[p][t] + PRODCOST[p][t] * produce[p][t];
}));
// 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)
prob.addConstraints(DIM_PRODUCTS, DIM_TIMES,
[&](auto p, auto t){
return sum(t + 1,
[&](auto s){
return produce[p][s];
})
>= sumDemand[p][0][t];
});
// 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 (int p = 0; p < DIM_PRODUCTS; ++p) {
for (int t = DIM_TIMES - 1; t >= 0; --t) {
prob.addConstraint(produce[p][t] <= sumDemand[p][t][DIM_TIMES - 1] * setup[p][t]);
}
}
// Capacity limits :
// forall(t in TIMES) Capacity(t) : = sum(p in PRODUCTS) produce(p, t) <= CAP(t)
prob.addConstraints(DIM_TIMES,
[&](auto t){
return sum(DIM_PRODUCTS,
[&](auto p) {
return produce[p][t];
})
<= CAP[t];
});
// Add a CutRound callback for separating our cuts.
Callback cb(produce, setup, sumDemand);
prob.callbacks.addCutRoundCallback(cb);
prob.writeProb("elsmanagedcuts.lp");
prob.optimize();
if (prob.attributes.getSolveStatus() == SolveStatus::Completed &&
prob.attributes.getSolStatus() == SolStatus::Optimal) {
std::cout << "Solved problem to optimality." << std::endl;
}
else {
std::cout << "Failed to solve problem with solvestatus "
<< to_string(prob.attributes.getSolveStatus())
<< " and solstatus "
<< to_string(prob.attributes.getSolStatus())
<< std::endl;
}
return 0;
}
|
