// (c) 2024-2024 Fair Isaac Corporation
#include <xpress.hpp>
#include <cstdlib>
#include <numeric>
using namespace xpress;
using namespace xpress::objects;
using xpress::objects::utils::sum;
/**
* Example for solving a MIP using lazily separated cuts/constraints.
*
* We solve a random instance of the symmetric TSP using lazily separated
* cuts/constraints.
*
* <p>
* The model is based on a directed graph G = (V,E).
* We have one binary variable x[e] for each edge e in E. That variable
* is set to 1 if edge e is selected in the optimal tour and 0 otherwise.
* </p>
* <p>
* The model contains only two explicit constraints:
* <pre>
for each v in V: sum(u in V : u != v) x[uv] == 1
for each v in V: sum(u in V : u != v) x[vu] == 1
</pre>
* These state that node u should have exactly one outgoing and exactly
* one incoming edge in a tour.
* </p>
* <p>
* The above constraints ensures that the selected edges form tours. However,
* it allows multiple tours, also known as subtours. So we need a constraint
* that requires that there is only one tour (which then necessarily hits
* all the nodes). This constraint is known as a subtour elimination constraint
* and is
* <pre>
sum(e in S) x[e] <= |S|-1 for each subtour S
* </pre>
*
* Since there are exponentially many subtours in a graph, this constraint is
* not stated explicitly. Instead we check for any solution that the optimizer
* finds, whether it satisfies the subtour elimination constraints. If it does
* then we accept the solution. Otherwise we reject the solution and augment the
* model by the violated subtours eliminiation constraint.
* </p>
* <p>
* This lazy addition of constraints is implemented using a preintsol callback
* that rejects any solution that violates a subtour elimination constraint and
* injects a violated subtour elimination constraint in case the solution
* candidate came from an integral node.
* </p>
* <p>
* An important thing to note about this strategy is that dual reductions have
* to be disabled. Since the optimizer does not see the whole model (subtour
* elimination constraints are only generated on the fly), dual reductions may
* cut off the optimal solution.
* </p>
*/
class TravelingSalesPerson {
/** Number of nodes in the instance. */
unsigned const nodes;
/** X coordinate of nodes. */
std::vector<double> nodeX;
/** Y coordinate of nodes. */
std::vector<double> nodeY;
/** Variables the edges. */
std::vector<std::vector<Variable>> x;
public:
/**
* Construct a new random instance.
*
* @param nodes The number of nodes in the instance.
* @param seed Random number seed.
*/
TravelingSalesPerson(unsigned nodes, unsigned seed = 0)
: nodes(nodes)
, nodeX(nodes)
, nodeY(nodes)
{
// Fill coordinates with random values.
std::srand(seed);
auto randCoord = [](){ return 4.0*static_cast<double>(rand())/RAND_MAX; };
std::generate(nodeX.begin(), nodeX.end(), randCoord);
std::generate(nodeY.begin(), nodeY.end(), randCoord);
}
/**
* Get the distance between two nodes.
*
* @param u First node.
* @param v Second node.
* @return The distance between <code>u</code> and <code>v</code>. The distance
* is symmetric.
*/
double distance(int u, int v) const {
return sqrt((nodeX[u] - nodeX[v]) * (nodeX[u] - nodeX[v]) + (nodeY[u] - nodeY[v]) * (nodeY[u] - nodeY[v]));
}
private:
/**
* Find the tour rooted at 0 in a solution. As a side effect, the tour is
* printed to the console.
*
* @param sol The current solution.
* @param from Stores the tour. <code>from[u]</code> yields the predecessor of
* <code>u</code> in the tour. If <code>from[u]</code> is negative
* then <code>u</code> is not in the tour. This parameter can be
* <code>null</code>.
* @return The length of the tour.
*/
unsigned findTour(std::vector<double> const &sol, std::vector<int> * from) {
std::unique_ptr<std::vector<int>> tmpFrom(nullptr);
if (!from) {
tmpFrom = std::make_unique<std::vector<int>>();
from = tmpFrom.get();
}
from->assign(nodes, -1);
unsigned node = 0;
unsigned used = 0;
std::cout << "0";
while (node != 0 || used == 0) {
// Find the edge leaving node
Variable edge = XpressProblem::NULL_VARIABLE;
for (unsigned i = 0; i < nodes; ++i) {
if (i != node && x[node][i].getValue(sol) > 0.5) {
std::cout << " -> " << i;
edge = x[node][i];
(*from)[i] = node;
node = i;
++used;
break;
}
}
if (!edge)
break;
}
std::cout << std::endl;
return used;
}
/**
* Integer solution check callback.
*/
void preIntsol(XpressProblem &prob, int soltype, int *p_reject, double * /*p_cutoff*/) {
std::cout << "Checking candidate solution ..." << std::endl;
// Get current solution and check whether it is feasible
std::vector<double> sol(nodes * nodes);
prob.getLpSol(sol, nullptr, nullptr, nullptr);
std::vector<int> from(nodes);
unsigned used = findTour(sol, &from);
std::cout << "Solution is ";
if (used < nodes) {
// The tour given by the current solution does not pass through
// all the nodes and is thus infeasible.
// If soltype is non-zero then we reject by setting
// p_reject.value=1.
// If instead soltype is zero then the solution came from an
// integral node. In this case we can reject by adding a cut
// that cuts off that solution. Note that we must NOT reject
// in that case because that would result in just dropping
// the node, no matter whether we add cuts or not.
std::cout << "infeasible (" << used << " edges)" << std::endl;
if (soltype != 0) {
*p_reject = 1;
} else {
// The tour is too short. Get the edges on the tour and
// add a subtour elimination constraint
LinExpression subtour = LinExpression::create();
for (unsigned u = 0; u < nodes; ++u) {
if (from[u] >= 0)
subtour.addTerm(x[from[u]][u]);
}
// We add the constraint. The solver must translate the
// constraint from the original space into the presolved
// space. This may fail (in theory). In that case the
// addCut() function will return non-zero.
if (prob.addCut(1, subtour.leq(used - 1)) != 0)
throw std::runtime_error("failed to presolve subtour elimination constraint");
}
} else {
std::cout << "feasible" << std::endl;
}
}
/** Create a feasible tour and add this as initial MIP solution. */
void createInitialTour(XpressProblem &prob) {
std::vector<Variable> variable(nodes);
std::vector<double> value(nodes);
// Create a tour that just visits each node in order.
for (unsigned i = 0; i < nodes; ++i) {
variable[i] = x[i][(i + 1) % nodes];
value[i] = 1.0;
}
prob.addMipSol(value, variable, "init");
}
public:
/**
* Solve the TSP represented by this instance.
*/
void solve() {
XpressProblem prob;
// Create variables. We create one variable for each edge in
// the (complete) graph. That is, we create variables from each
// node u to all other nodes v. We even create a variable for
// the self-loop from u to u, but that variable is fixed to 0.
// x[u][v] gives the variable that represents edge uv.
// All variables are binary.
x = prob.addVariables(nodes, nodes)
.withType(ColumnType::Binary)
.withName([](auto i, auto j){ return format("x_%d_%d", i, j); })
.withUB([](auto i, auto j){ return (i == j) ? 0.0 : 1.0; })
.toArray();
// Objective. All variables are in the objective and their
// respective coefficient is the distance between the two nodes.
prob.setObjective(sum(nodes,
[&](auto u){
return sum(nodes,
[&](auto v){ return x[u][v] * distance(u, v); }); }),
ObjSense::Minimize);
// Create a vector with indices 0, ..., nodes-1
std::vector<int> indices(nodes);
std::iota(indices.begin(), indices.end(), 0);
// Constraint: In the graph that is induced by the selected
// edges, each node should have exactly one outgoing
// and exactly one incoming edge.
// These are the only constraints we add explicitly.
// Subtour elimination constraints are added
// dynamically via a callback.
prob.addConstraints(nodes,
[&](auto u){ return sum(indices,
[&](auto v){ return (u == v ? 0.0 : 1.0) * x[u][v]; })
== 1.0; });
prob.addConstraints(nodes,
[&](auto u){ return sum(indices,
[&](auto v){ return (u == v ? 0.0 : 1.0) * x[v][u]; })
== 1.0; });
// Create a starting solution.
// This is optional but having a feasible solution available right
// from the beginning can improve optimizer performance.
createInitialTour(prob);
// Write out the model in case we want to look at it.
prob.writeProb("travelingsalesperson.lp", "l");
// We don't have all constraints explicitly in the matrix, hence
// we must disable dual reductions. Otherwise MIP presolve may
// cut off the optimal solution.
prob.controls.setMipDualReductions(0);
// Add a callback that rejects solutions that do not satisfy
// the subtour constraints.
prob.callbacks.addPreIntsolCallback(std::bind(&TravelingSalesPerson::preIntsol, this,
std::placeholders::_1,
std::placeholders::_2,
std::placeholders::_3,
std::placeholders::_4
));
// Add a message listener to display log information.
prob.callbacks.addMessageCallback(XpressProblem::console);
prob.optimize();
if (prob.attributes.getSolStatus() != SolStatus::Optimal)
throw std::runtime_error("failed to solve");
auto sol = prob.getSolution();
// Print the optimal tour.
std::cout << "Tour with length " << prob.attributes.getMipBestObjVal() << std::endl;
findTour(sol, nullptr);
x.clear(); // We are done with the variables
}
};
int main() {
try {
TravelingSalesPerson(10).solve();
return 0;
}
catch (std::exception& e) {
std::cout << "Exception: " << e.what() << std::endl;
return -1;
}
}
|
/***********************************************************************
Xpress Optimizer Examples
=========================
file txp.cc
```````````
Solve a MIP using cuts/constraints that are lazily separated.
We take a random instance of the symmetric TSP and solve that using
lazily separated constraints.
(c) 2024-2024 Fair Isaac Corporation
***********************************************************************/
#include <xpress.hpp>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <iostream>
using xpress::XPRSProblem;
/** Example for solving a MIP using lazily separated cuts/constraints.
*
* We solve a random instance of the symmetric TSP using lazily separated
* cuts/constraints.
*
* <p>
* The model is based on a graph G = (V,E).
* We have one binary variable x[e] for each edge e in E. That variable
* is set to 1 if edge e is selected in the optimal tour and 0 otherwise.
* </p>
* <p>
* The model contains only one explicit set of constraints:
* <pre>
for each v in V: sum(u in V : u != v) x[uv] == 2
</pre>
* This states that from all the edges incident to a node u, exactly two
* must be selected in the optimal tour.
* </p>
* <p>
* The above constraints ensures that the selected edges form tours. However,
* it allows multiple tours, also known as subtours. So we need a constraint
* that requires that there is only one tour (which then necessarily hits
* all the nodes). This constraint is known as subtour elimination constraint
* and is
* <pre>
sum(e in S) x[e] <= |S|-1 for each subtour S
</pre>
* Since there are exponentially many subtours in a graph, this constraint
* is not stated explicitly. Instead we check for any solution that the
* optimizer finds, whether it satisfies the subtour elimination constraints.
* If it does then we accept the solution. Otherwise we reject the solution
* and augment the model by the violated subtour eliminiation constraints.
* </p>
* <p>
* This lazy addition of constraints is implemented using
* a preintsol callback that rejects any solution that violates a
* subtour elimination constraint and injects a violated subtour elimination
* constraint in case the solution candidate came from an integral node.
* </p>
* <p>
* An important thing to note about this strategy is that dual reductions
* have to be disabled. Since the optimizer does not see the whole model
* (subtour elimination constraints are only generated on the fly), dual
* reductions may cut off the optimal solution.
* </p>
*/
class TSP {
public:
/** Number of nodes in the instance. */
unsigned const nodes;
/** Number of edges in the instance. */
unsigned const edges;
/** X coordinate of nodes. */
std::vector<double> nodeX;
/** Y coordinate of nodes. */
std::vector<double> nodeY;
/** Variable indices for the edges. */
std::vector<std::vector<int>> x;
public:
/** Construct a new random instance.
* @param nodes The number of nodes in the instance.
* @param seed Random number seed.
*/
TSP(unsigned n, int seed = 0)
: nodes(n)
, edges((n * (n - 1)) / 2)
, nodeX(n)
, nodeY(n)
, x(n, std::vector<int>(n))
{
// Fill coordinates with random values.
std::srand(seed);
auto randCoord = [](){ return 4.0*static_cast<double>(rand())/RAND_MAX; };
std::generate(nodeX.begin(), nodeX.end(), randCoord);
std::generate(nodeY.begin(), nodeY.end(), randCoord);
}
/** Get the distance between two nodes.
* @param u First node.
* @param v Second node.
* @return The distance between <code>u</code> and <code>v</code>.
* The distance is symmetric.
*/
double distance(unsigned u, unsigned v) const {
return sqrt((nodeX[u] - nodeX[v]) * (nodeX[u] - nodeX[v]) +
(nodeY[u] - nodeY[v]) * (nodeY[u] - nodeY[v]));
}
/** Find the tour rooted at 0 in a solution.
* As a side effect, the tour is printed to the console.
* @param sol The current solution.
* @param from Stores the tour. <code>from[u]</code> yields the
* predecessor of <code>u</code> in the tour. If
* <code>from[u]</code> is negative then <code>u</code>
* is not in the tour.
* This parameter can be <code>null</code>.
* @return The length of the tour.
*/
unsigned findTour(std::vector<double> const &sol, std::vector<int> *from) const {
std::unique_ptr<std::vector<int>> allocatedFrom;
if (!from) {
allocatedFrom.reset(new std::vector<int>(nodes));
from = allocatedFrom.get();
}
from->assign(nodes, -1);
std::vector<bool> inTour(edges); // Marks edges on the subtour
unsigned u = 0;
unsigned used = 0;
std::cout << "0";
do {
for (unsigned v = 0; v < nodes; ++v) {
if (u == v) // no self-loops
continue;
else if ((*from)[v] != -1) // node already on tour
continue;
else if (sol[x[u][v]] < 0.5) // edge not selected in solution
continue;
else if (inTour[x[u][v]]) // edge already on tour
continue;
else {
std::cout << " -> " << v;
inTour[x[u][v]] = true;
(*from)[v] = u;
used += 1;
u = v;
break;
}
}
} while (u != 0);
std::cout << std::endl;
return used;
}
/** Integer solution check callback.
*/
class PreIntsolCallback {
TSP const &data;
public:
PreIntsolCallback(TSP const &d) : data(d) {}
void operator()(XPRSProblem &prob, int soltype, int *p_reject, double *) {
std::cout << "Checking candidate solution ..." << std::endl;
// Get current solution and check whether it is feasible
std::vector<double> sol(data.edges);
prob.getLpSol(sol, nullptr, nullptr, nullptr);
std::vector<int> from(data.nodes);
unsigned used = data.findTour(sol, &from);
std::cout << "Solution is " << std::flush;
if (used < data.nodes) {
// The tour given by the current solution does not pass through
// all the nodes and is thus infeasible.
// If soltype is non-zero then we reject.
// If instead soltype is zero then the solution came from an
// integral node. In this case we can reject by adding a cut
// that cuts off that solution. Note that we must NOT return
// null in that case because that would result in just dropping
// the node, no matter whether we add cuts or not.
std::cout <<"infeasible (" << used << " edges)" << std::endl;
if (soltype != 0) {
*p_reject = 1;
}
else {
// The tour is too short. Get the edges on the tour and
// add a subtour elimination constraint
std::vector<int> ind(used);
for (unsigned u = 0, next = 0; u < data.nodes; ++u) {
if (from[u] >= 0)
ind[next++] = data.x[u][from[u]];
}
std::vector<double> val(used, 1.0);
// Since we created the constraint in the original space,
// we must crush it to the presolved space before we can
// add it.
int status = -1;
double prerhs = 0.0;
int maxcoefs = prob.attributes.getCols();
int coefs = -1;
std::vector<int> preind(maxcoefs);
std::vector<double> preval(maxcoefs);
prob.presolveRow('L', used, ind, val, used - 1,
maxcoefs, &coefs, preind, preval, &prerhs,
&status);
if (status == 0)
prob.addCuts(1, std::vector<int>{0}, "L", &prerhs,
std::vector<int>{0, coefs},
preind, preval);
else
throw std::runtime_error("failed to presolve subtour elimination constraint");
}
}
else {
std::cout << "feasible" << std::endl;
// We accept the solution, so nothing to do
}
}
};
/** Create a feasible tour and add this as initial MIP solution. */
void createInitialTour(XPRSProblem &prob) {
std::vector<int> ind(nodes);
std::vector<double> val(nodes);
// Create a tour that just visits each node in order.
for (unsigned i = 0; i < nodes; ++i) {
ind[i] = x[i][(i + 1) % nodes];
val[i] = 1.0;
}
prob.addMipSol(nodes, val, ind, "initial");
}
/** Solve the TSP represented by this instance.
*/
void solve() {
XPRSProblem prob;
// Create variables. We create one variable for each edge in
// the (complete) graph. x[u][v] gives the index of the variable
// that represents edge uv. Since edges are undirected, x[v][u]
// gives the same variable.
// All variables are binary.
for (unsigned i = 0; i < nodes; ++i) {
for (unsigned j = i + 1; j < nodes; ++j) {
x[j][i] = x[i][j] = prob.attributes.getCols();
prob.addCols(1, 0, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr);
prob.chgColType(1, &x[i][j], "B");
}
}
// Objective. All variables are in the objective and their
// respective coefficient is the distance between the two nodes.
std::vector<int> objInd(edges);
std::vector<double> objVal(edges);
for (unsigned i = 0, nz = 0; i < nodes; ++i) {
for (unsigned j = i + 1; j < nodes; ++j) {
objInd[nz] = x[i][j];
objVal[nz] = distance(i, j);
++nz;
}
}
prob.chgObj(edges, objInd, objVal);
// Constraint: In the graph that is induced by the selected
// edges, each node should have degree 2.
// This is the only constraint we add explicitly.
// Subtour elimination constraints are added
// dynamically via a callback.
for (unsigned u = 0; u < nodes; ++u) {
std::vector<int> ind(nodes - 1);
for (unsigned v = 0, next = 0; v < nodes; ++v) {
if (u != v)
ind[next++] = x[u][v];
}
std::vector<double> val(nodes - 1, 1.0);
prob.addRows(1, nodes - 1, "E", std::vector<double>({2}), nullptr, std::vector<int>({0}), ind, val);
}
// Create a starting solution.
// This is optional but having a feasible solution available right
// from the beginning can improve optimizer performance.
createInitialTour(prob);
// Write out the model in case we want to look at it.
prob.writeProb("tsp.lp", "l");
// We don't have all constraints explicitly in the matrix, hence
// we must disable dual reductions. Otherwise MIP presolve may
// cut off the optimal solution.
prob.controls.setMipDualReductions(0);
// Add a callback that rejects solutions that do not satisfy
// the subtour constraints.
PreIntsolCallback cb(*this);
prob.addPreIntsolCallback(cb);
// Add a message listener to display log information.
prob.addMessageCallback(XPRSProblem::console);
prob.optimize("", nullptr, nullptr);
auto sol = prob.getSolution();
// Print the optimal tour.
std::cout << "Tour with length " << prob.attributes.getObjVal() << std::endl;
findTour(sol, nullptr);
}
};
int main(void) {
try {
TSP(10).solve();
return 0;
}
catch (std::exception& e) {
std::cout << "Exception: " << e.what() << std::endl;
return -1;
}
}
|
