Untitled

mail@pastecode.io avatar
unknown
c_cpp
2 months ago
14 kB
2
Indexable
Never
#include "bendershelper.h"

BendersHelper::BendersHelper(SVRPInstance &instance, SVRPSolution &solution,
                             SCIP *scip, ArcSCIPVarMap &x, DNodeSCIPVarMap &y,
                             Params &params)
    : 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;
}
Leave a Comment