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