Solving a robust facility location problem by its deterministic equivalent

Problem Definition

We consider a facility location problem with uncertain demand. 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 introduces uncertainty in the customers’ demands.

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 uncertainty in customer demands is controlled by a parameter \(\Gamma\).

In this robust variant, we consider that the demands are uncertain and can be expressed as \(d_j(\xi) = d_j(1 + 0.2 \xi_j)\) with \(\xi\) being an unknown vector taken in the uncertainty set

\[\Xi := \left\{ \xi\in[ 0, 1 ]^{|V_2|} : \sum_{j\in V_2} \xi_j \le \Gamma \right\}.\]

The goal is to minimize the total cost of opening facilities and serving customers, considering the worst-case demand scenario.

The robust version can be formulated as

\[\begin{split}\begin{align*} \min_{x,y} \quad & \sum_{i\in V_1} f_i x_i + \sum_{i\in V_1} \sum_{j\in V_2} t_{ij} y_{ij} \\ \text{s.t.} \quad & \sum_{j\in V_2} d_j(\xi) y_{ij} \le q_i, && i\in V_1, \quad \text{for all }\xi\in\Xi, \\ & \sum_{i\in V_1} y_{ij} = 1, && j\in V_2, \\ & y_{ij} \le x_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 reformulate this robust facility location problem as a deterministic one. Then, we will solve it using Gurobi.

Hint

Here, the deterministic reformulation reads as follows.

\[\begin{split}\begin{align*} \min_{x,y,u.v} \quad & \sum_{i\in V_1} f_i x_i + \sum_{i\in V_1} \sum_{j\in V_2} t_{ij} y_{ij} \\ \text{s.t.} \quad & \sum_{j\in V_2} d_j y_{ij} + \Gamma u + \sum_{j\in V_2} v_j \le q_i, && i\in V_1, \\ & u + v_j \ge 0.2 d_j y_j, && j\in V_2, \\ & \sum_{i\in V_1} y_{ij} = 1, && j\in V_2, \\ & u, v_j \ge 0, && j\in V_2, \\ & y_{ij} \le x_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}\]

Implementation

//
// Created by henri on 06/04/23.
//
#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"

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("robust-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 = 1;
    const auto xi = uncertainty_set.add_vars(Dim<1>(n_customers), 0., 1., Continuous, 0., "xi");
    uncertainty_set.add_ctr(idol_Sum(i, Range(n_customers), xi[i]) <= 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) {

        const auto c = 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) {
            description.set_uncertain_mat_coeff(c, y[i][j], 0.2 * instance.demand(j) * xi[j]);
        }

    }

    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) {
            description.set_uncertain_obj(y[i][j], 0.2 * instance.per_unit_transportation_cost(i, j) * instance.demand(j) * xi[j]);
        }
    }

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

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

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

    return 0;
}