Using affine decision rules to solve a two-stage facility location problem with facility disruption

Problem Definition

We consider a facility location problem with facility disruption. Given a set of potential facility locations \(V_1\) and a set of customers \(V_2\), the goal is to select a subset of facility locations to activate in order to serve all customers’ demand, while minimizing the total cost. This version considers random disruptions of the opened facilities making them unusable.

Note that there is also an example for the deterministic version of the FLP using Column Generation.

Each facility \(i\in V_1\) has an opening cost \(f_i\) and a maximum capacity \(q_i\). Each customer \(j\in V_2\) has a demand \(d_j\). The unitary cost for serving customer \(j\in V_2\) from facility \(i\in V_1\) is \(t_{ij}\). The maximum number of facilities that can be disrupted is noted \(\Gamma\).

The two-stage robust facility location problem can be formulated as

\[\begin{align} \min_{ x\in \{0,1\}^{|V_1|} } \quad & \left\{ \sum_{i\in V_1} f_i x_i + \max_{\xi\in \Xi} \min_{y\in Y(x,\xi)} \quad \sum_{i\in V_1} \sum_{j\in V_2} t_{ij} y_{ij} \right\}, \end{align}\]

in which \(\Xi := \left\{ \xi\in\{0,1\}^{|V_1|} : e^\top\xi \le \Gamma \right\}\) and \(Y(x,\xi)\) is the feasible set of the second-stage problem defined as

\[\begin{split}\begin{align} & \sum_{j\in V_2} d_j y_{ij} \le q_i, && i\in V_1,\\ & \sum_{i\in V_1} y_{ij} \ge 1, && j\in V_2, \\ & y_{ij} \le x_i && i\in V_1, j\in V_2, \\ & y_{ij} \le 1 - \xi_i && i\in V_1, j\in V_2, \\ & x_i \in \{0,1\}, && i\in V_1, \\ & y_{ij} \in \{0,1\}, && i\in V_1, j\in V_2. \end{align}\end{split}\]

In this example, we will approximate the robust problem using affine decision rules and solve it using Gurobi. To this end, we will restrict the adjustable decisions to be affine functions of the uncertain parameters, i.e., \(y(\xi) = \bar y + Y\xi\).

Implementation

//
// Created by henri on 28.11.24.
//
#include <iostream>
#include "idol/modeling.h"
#include "idol/mixed-integer/problems/facility-location-problem/FLP_Instance.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/BranchAndBound.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/branching-rules/factories/PseudoCost.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/node-selection-rules/factories/BestEstimate.h"
#include "idol/mixed-integer/optimizers/wrappers/HiGHS/HiGHS.h"
#include "idol/mixed-integer/optimizers/wrappers/GLPK/GLPK.h"
#include "idol/mixed-integer/optimizers/wrappers/Gurobi/Gurobi.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/node-selection-rules/factories/BestBound.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/branching-rules/factories/MostInfeasible.h"
#include "idol/mixed-integer/optimizers/callbacks/ReducedCostFixing.h"
#include "idol/robust/modeling/Description.h"
#include "idol/robust/optimizers/deterministic/Deterministic.h"
#include "idol/robust/optimizers/affine-decision-rule/AffineDecisionRule.h"

using namespace idol;

int main(int t_argc, const char** t_argv) {

    Env env;

    // Read instance
    const auto instance = Problems::FLP::read_instance_1991_Cornuejols_et_al("adr-facility-location.data.txt");
    const unsigned int n_customers = instance.n_customers();
    const unsigned int n_facilities = instance.n_facilities();

    // Uncertainty set
    Model uncertainty_set(env);
    const double Gamma = .2;
    const auto xi = uncertainty_set.add_vars(Dim<2>(n_facilities, n_customers), 0., 1., Continuous, 0., "xi");
    uncertainty_set.add_ctr(idol_Sum(i, Range(n_facilities), idol_Sum(j, Range(n_customers), xi[i][j])) <= Gamma);

    // Make model
    Model model(env);
    Robust::Description description(uncertainty_set);

    auto x = model.add_vars(Dim<1>(n_facilities), 0., 1., Binary, 0., "x");
    auto y = model.add_vars(Dim<2>(n_facilities, n_customers), 0., 1., Continuous, 0., "y");

    for (unsigned int i = 0 ; i < n_facilities ; ++i) {
        model.add_ctr(idol_Sum(j, Range(n_customers), instance.demand(j) * y[i][j]) <= instance.capacity(i));
    }

    for (unsigned int j = 0 ; j < n_customers ; ++j) {
        model.add_ctr(idol_Sum(i, Range(n_facilities), y[i][j]) >= 1);
    }

    for (unsigned int i = 0 ; i < n_facilities ; ++i) {
        for (unsigned int j = 0 ; j < n_customers ; ++j) {
            model.add_ctr(y[i][j] <= x[i]);
        }
    }

    model.set_obj_expr(idol_Sum(i, Range(n_facilities),
                                instance.fixed_cost(i) * x[i]
                                + idol_Sum(j, Range(n_customers),
                                           instance.per_unit_transportation_cost(i, j) *
                                           instance.demand(j) *
                                           y[i][j]
                                )
                       )
    );

    for (unsigned int i = 0 ; i < n_facilities ; ++i) {
        for (unsigned int j = 0; j < n_customers; ++j) {
            const auto c = model.add_ctr(y[i][j] <= 1);
            description.set_uncertain_rhs(c, -xi[i][j]); // models y_ij <= 1 - xi_j
            description.set_stage(y[i][j], 1); // models y_ij as a second-stage variable
        }
    }

    model.use(Gurobi());
    model.optimize();
    std::cout << "Deterministic Problem has value: " << model.get_best_obj() << std::endl;

    /*
     * Alternatively, you could also do
     * const auto adr_result = Robust::AffineDecisionRule::make_model(model, description);
     * std::cout << Robust::Description::View(adr_result.model, description) << std::endl;
     */

    model.use(
            Robust::AffineDecisionRule(description)
                    .with_deterministic_optimizer(Gurobi().with_logs(true))
    );
    model.optimize();
    std::cout << "Affine Decision Rule Problem has value: " << model.get_best_obj() << std::endl;

    model.use(
            Robust::Deterministic(description)
                    .with_deterministic_optimizer(Gurobi())
    );
    model.optimize();
    std::cout << "Robust Problem has value: " << model.get_best_obj() << std::endl;

    return 0;
}