Untitled
unknown
c_cpp
a year ago
14 kB
4
Indexable
#include "bendershelper.h"
BendersHelper::BendersHelper(SVRPInstance &instance, SVRPSolution &solution,
SCIP *scip, ArcSCIPVarMap &x, DNodeSCIPVarMap &y,
Params ¶ms)
: instance(instance), x(x), y(y), params(params), networkCost(network),
networkRecourseCost(network), networkIdMap(network),
networkDemandMap(network), model(env), flowVar(network),
xConsGEQ(instance.d_g), xConsLEQ(instance.d_g), nodeCons(network) {
EpsForIntegrality = 1e-5;
nCuts = 0;
nRounds = 0;
createPoissonNetwork(scip);
createGurobiModel();
}
// Create basic gurobi model.
void BendersHelper::createGurobiModel() {
// Initialize Gurobi formulation
model.set(GRB_StringAttr_ModelName, "Benders-Formulation");
model.set(GRB_IntAttr_ModelSense, GRB_MINIMIZE);
// model.getEnv().set(GRB_IntParam_LazyConstraints, 1);
// model.set(GRB_IntParam_InfUnbdInfo, 1);
model.set(GRB_IntParam_LogToConsole, 0);
model.set(GRB_DoubleParam_MIPGapAbs, 0.01);
model.set(GRB_IntParam_Method, 1);
// model.set(GRB_IntParam_DualReductions, 0);
// Initialize gurobi variables.
z = model.addVar(0.0, GRB_INFINITY, 1.0, GRB_CONTINUOUS, "z");
for (ArcIt a(network); a != INVALID; ++a) {
flowVar[a] = model.addVar(
0.0, GRB_INFINITY, 0.0, GRB_CONTINUOUS,
"f_(" + to_string(networkIdMap[network.source(a)]) + "," +
to_string(networkDemandMap[network.source(a)]) + ")->(" +
to_string(networkIdMap[network.target(a)]) + "," +
to_string(networkDemandMap[network.target(a)]) + ")");
}
// Flow conservation constraints.
for (DNodeIt v(network); v != INVALID; ++v) {
double RHS = 0.0;
if (v == source) {
RHS = -instance.k;
} else if (v == target) {
RHS = instance.k;
}
GRBLinExpr flowExpr = 0;
for (InArcIt a(network, v); a != INVALID; ++a) {
flowExpr += flowVar[a];
}
for (OutArcIt a(network, v); a != INVALID; ++a) {
flowExpr -= flowVar[a];
}
nodeCons[v] = model.addConstr(flowExpr == RHS);
}
// x-map constraints.
ArcGRBExprMap xExpr(instance.d_g, 0);
for (ArcIt a(network); a != INVALID; ++a) {
Arc originalArc =
findArc(instance.d_g,
instance.d_g.nodeFromId(networkIdMap[network.source(a)]),
instance.d_g.nodeFromId(networkIdMap[network.target(a)]));
xExpr[originalArc] += flowVar[a];
}
for (ArcIt a(instance.d_g); a != INVALID; ++a) {
xConsGEQ[a] = model.addConstr(xExpr[a] + z >= 0.0);
xConsLEQ[a] = model.addConstr(xExpr[a] - z <= 0.0);
}
// Recourse constraint.
GRBLinExpr recourseExpr = 0;
for (ArcIt a(network); a != INVALID; ++a) {
recourseExpr += networkRecourseCost[a] * flowVar[a];
}
recourseCons = model.addConstr(recourseExpr - z <= 0.0);
model.update();
}
// Create network for computing normal recourse costs.
void BendersHelper::createPoissonNetwork(SCIP *scip) {
source = network.addNode();
target = network.addNode();
networkIdMap[source] = 0;
networkIdMap[target] = 0;
networkDemandMap[source] = 0;
networkDemandMap[target] = instance.capacity;
map<tuple<int, int>, DNode> stateMap;
stateMap[{0, 0}] = source;
int addedArcs = 0;
// Do a dfs from the state (0,0)
std::stack<tuple<int, int>> queue;
queue.push({0, 0});
while (queue.size() > 0) {
tuple<int, int> currState = queue.top();
DNode u = instance.d_g.nodeFromId(get<0>(currState));
int accDemand = get<1>(currState);
DNode currNode = stateMap[currState];
queue.pop();
for (OutArcIt a(instance.d_g, u); a != INVALID; ++a) {
DNode v = instance.d_g.target(a);
int nextDemand = int(instance.d_demand[v]);
if (accDemand + nextDemand > instance.capacity) {
continue;
}
Arc depotArc = findArc(instance.d_g,
instance.d_g.nodeFromId(instance.depot), v);
double recourseCost = 2.0 * instance.d_weight[depotArc];
// Compute the failure probability
double prob = 0.0;
double aux;
int t = 1;
while (true) {
aux = instance.poissonFailureProbability(t, accDemand,
nextDemand);
if (aux >= 1e-5) {
prob += aux;
t++;
} else {
break;
}
}
tuple<int, int> nextState = {instance.d_g.id(v),
accDemand + nextDemand};
DNode nextNode = INVALID;
if (instance.d_g.id(v) == instance.depot) {
nextNode = target;
} else if (stateMap.find(nextState) != stateMap.end()) {
nextNode = stateMap[nextState];
} else {
nextNode = network.addNode();
networkIdMap[nextNode] = instance.d_g.id(v);
networkDemandMap[nextNode] = accDemand + nextDemand;
stateMap[nextState] = nextNode;
queue.push(nextState);
}
Arc networkArc = network.addArc(currNode, nextNode);
networkCost[networkArc] = 0.0;
networkRecourseCost[networkArc] = prob * recourseCost;
addedArcs++;
}
}
std::cout << "number of nodes in the network: " << stateMap.size()
<< std::endl;
std::cout << "number of arcs in the network: " << addedArcs << std::endl;
}
SCIP_RETCODE BendersHelper::addBendersInequalities(SCIP *scip,
SCIP_CONSHDLR *conshdlr,
SCIP_SOL *sol,
SCIP_RESULT *result) {
auto started = chrono::high_resolution_clock::now();
ArcValueMap alphaValues(instance.d_g);
// Vector to store found rows.
std::vector<SCIP_ROW *> addedRows;
// Start probing mode.
unsigned int lperror, cutoff;
SCIP_CALL(SCIPstartProbing(scip));
// Create a new node in probing mode and solve it.
SCIP_CALL(SCIPnewProbingNode(scip));
SCIP_NODE *probingNode = SCIPgetCurrentNode(scip);
SCIP_CALL(SCIPsolveProbingLP(scip, -1, &lperror, &cutoff));
try {
while (addedRows.size() < 200) {
// Get current recourse cost.
double currRecourse = 0.0;
for (DNodeIt v(instance.d_g); v != INVALID; ++v) {
currRecourse += SCIPgetSolVal(scip, sol, y[v]);
}
// Change variable bounds.
// for (ArcIt a(network); a != INVALID; ++a) {
// Arc originalArc = findArc(
// instance.d_g,
// instance.d_g.nodeFromId(networkIdMap[network.source(a)]),
// instance.d_g.nodeFromId(networkIdMap[network.target(a)]));
// if (SCIPgetSolVal(scip, sol, x[originalArc]) <=
// EpsForIntegrality) {
// flowVar[a].set(GRB_DoubleAttr_UB, 0.0);
// } else {
// flowVar[a].set(GRB_DoubleAttr_UB, 1.0);
// }
// }
// Change constraints bounds.
for (ArcIt a(instance.d_g); a != INVALID; ++a) {
xConsGEQ[a].set(GRB_DoubleAttr_RHS,
SCIPgetSolVal(scip, sol, x[a]));
xConsLEQ[a].set(GRB_DoubleAttr_RHS,
SCIPgetSolVal(scip, sol, x[a]));
}
recourseCons.set(GRB_DoubleAttr_RHS, currRecourse);
// Update model and solve.
model.update();
model.write("prlp.lp");
model.optimize();
if (model.get(GRB_IntAttr_Status) != GRB_OPTIMAL) {
std::cout
<< "Benders error: problem was not solved to optimality"
<< std::endl;
exit(1);
}
for (ArcIt a(network); a != INVALID; ++a) {
DNode u = network.source(a);
DNode v = network.target(a);
double pi_u = nodeCons[u].get(GRB_DoubleAttr_Pi);
double pi_v = nodeCons[v].get(GRB_DoubleAttr_Pi);
Arc originalArc = findArc(
instance.d_g, instance.d_g.nodeFromId(networkIdMap[u]),
instance.d_g.nodeFromId(networkIdMap[v]));
double LHS = xConsGEQ[originalArc].get(GRB_DoubleAttr_Pi) +
xConsLEQ[originalArc].get(GRB_DoubleAttr_Pi) +
recourseCons.get(GRB_DoubleAttr_Pi) *
networkRecourseCost[a] -
pi_u + pi_v;
if (LHS >= EpsForIntegrality) {
// Check dual feasibility.
std::cout << "------------------------" << std::endl;
std::cout << "constraint for "
<< flowVar[a].get(GRB_StringAttr_VarName)
<< std::endl;
std::cout << "pi_u: " << pi_u << std::endl;
std::cout << "pi_v: " << pi_v << std::endl;
std::cout << "alpha1: "
<< xConsGEQ[originalArc].get(GRB_DoubleAttr_Pi)
<< std::endl;
std::cout << "alpha2: "
<< xConsLEQ[originalArc].get(GRB_DoubleAttr_Pi)
<< std::endl;
std::cout
<< "gamma: " << recourseCons.get(GRB_DoubleAttr_Pi)
<< std::endl;
std::cout << "recourse: " << networkRecourseCost[a]
<< std::endl;
std::cout << "LHS: " << LHS << std::endl;
std::cout
<< "reduced cost: " << flowVar[a].get(GRB_DoubleAttr_RC)
<< std::endl;
}
}
// Check objective.
if (model.get(GRB_DoubleAttr_ObjVal) <= EpsForIntegrality) {
break;
}
// Build alpha values.
for (ArcIt a(instance.d_g); a != INVALID; ++a) {
alphaValues[a] = xConsGEQ[a].get(GRB_DoubleAttr_Pi) +
xConsLEQ[a].get(GRB_DoubleAttr_Pi);
}
double RHS = -instance.k * nodeCons[source].get(GRB_DoubleAttr_Pi) +
instance.k * nodeCons[target].get(GRB_DoubleAttr_Pi);
double gammaValue = recourseCons.get(GRB_DoubleAttr_Pi);
// Create row to add to master LP.
SCIP_ROW *row;
SCIP_CALL(SCIPcreateEmptyRowConshdlr(
scip, &row, conshdlr, "Benders cut", -SCIPinfinity(scip), -RHS,
FALSE, FALSE, TRUE));
SCIP_CALL(SCIPcacheRowExtensions(scip, row));
// Create probing constraint.
SCIP_CONS *cons;
SCIP_CALL(SCIPcreateConsLinear(scip, &cons, "probingConstraint", 0,
NULL, NULL, -SCIPinfinity(scip),
-RHS, TRUE, FALSE, TRUE, TRUE, FALSE,
TRUE, FALSE, FALSE, FALSE, TRUE));
double LHS = 0.0;
for (DNodeIt v(instance.d_g); v != INVALID; ++v) {
SCIPaddVarToRow(scip, row, y[v], gammaValue);
SCIPaddCoefLinear(scip, cons, y[v], gammaValue);
LHS += gammaValue * SCIPgetSolVal(scip, sol, y[v]);
}
for (ArcIt a(instance.d_g); a != INVALID; ++a) {
SCIPaddVarToRow(scip, row, x[a], alphaValues[a]);
SCIPaddCoefLinear(scip, cons, x[a], alphaValues[a]);
LHS += alphaValues[a] * SCIPgetSolVal(scip, sol, x[a]);
}
// Add row to the vector.
addedRows.push_back(row);
// Add constraint to probing node.
SCIP_CALL(SCIPaddConsNode(scip, probingNode, cons, NULL));
// Solve probing LP.
SCIP_CALL(SCIPsolveProbingLP(scip, -1, &lperror, &cutoff));
if (cutoff) {
std::cout << "cutoff node" << std::endl;
break;
}
}
} catch (GRBException e) {
std::cout << "Gurobi error! Error number: " << e.getErrorCode()
<< std::endl;
std::cout << e.getMessage() << std::endl;
} catch (...) {
std::cout << "Unknown error!" << std::endl;
}
// Exit probing mode.
SCIP_CALL(SCIPbacktrackProbing(scip, 0));
SCIP_CALL(SCIPendProbing(scip));
// Add rows to the master LP.
for (SCIP_ROW *row : addedRows) {
// Add the cut to the LP.
SCIP_CALL(SCIPflushRowExtensions(scip, row));
SCIP_Bool infeasible = false;
SCIP_CALL(SCIPaddPoolCut(scip, row));
SCIP_CALL(SCIPaddRow(scip, row, TRUE, &infeasible));
// SCIP_CALL(SCIPprintRow(scip, row, NULL));
if (*result != SCIP_CUTOFF) {
*result = SCIP_SEPARATED;
}
if (infeasible) {
*result = SCIP_CUTOFF;
}
SCIP_CALL(SCIPreleaseRow(scip, &row));
nCuts++;
}
auto done = chrono::high_resolution_clock::now();
std::cout << "Benders execution time: "
<< double(
chrono::duration_cast<chrono::microseconds>(done - started)
.count()) /
1e6
<< "s" << std::endl;
std::cout << "Added benders cuts: " << addedRows.size() << std::endl;
nRounds++;
return SCIP_OKAY;
}Editor is loading...
Leave a Comment