# Blood supply chain design in disaster management: a comparison of risk-averse measures in robust optimization approaches

Komlanvi Parfait Ametana, Mehdi Amiri-Aref, Olga Battaïa, François Clautiaux

Boris Detienne, Laurent Facq

► To cite this version:

Komlanvi Parfait Ametana, Mehdi Amiri-Aref, Olga Battaïa, François Clautiaux, Boris Detienne, et al.. Blood supply chain design in disaster management: a comparison of risk-averse measures in robust optimization approaches. 2025. hal-04634185v2

HAL Id: hal-04634185

https://hal.science/hal-04634185v2

Preprint submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Blood supply chain design in disaster management: a comparison of risk-averse measures in robust optimization approaches

Komlanvi Parfait Ametana¹

François Clautiaux¹

Mehdi Amiri-Aref²

Boris Detienne¹

Olga Battaïa²

Laurent Facq¹

January 16, 2025

¹IMB UMR CNRS 5251, Inria Bordeaux Sud-Ouest, Université de Bordeaux,
200 Avenue de la Vieille Tour, Talence, 33405, France

² Kedge Business School, Bordeaux, France

## Abstract

Disasters, both natural and man-made, pose significant challenges to healthcare systems, particularly in managing blood supply chains effectively. Ensuring an adequate and timely supply of blood products becomes crucial during these crises to save lives and mitigate the impact of catastrophic events. This paper develops two-stage and multi-period risk aversion models for the blood supply chain in disaster management, focusing on mean-CVaR and worst-case criterion measures. We leverage advanced algorithms based on Benders decomposition and column-and-row generation techniques to produce practical solution methods which enable solving large-scale problem instances effectively. Numerical results demonstrate that these methods not only reduce the computational burden but also significantly enhance the solver's ability to explore feasible solutions efficiently. Our computational experiments show that a coarse approximation of uncertainty is more effective than a fine approximation and that good preparation enables a more effective response to an emergency. Our findings provide insights for policymakers, healthcare practitioners, and logistics professionals to improve blood supply chain management strategies in disaster response and preparedness efforts.

**Keywords**— Blood supply chain; Disaster management; Risk measurement; Two-stage robust optimization.

## 1 Introduction

Since the 1950s, both the number and magnitude of disasters triggered by natural hazards have continuously increased, with the number of affected people rising proportionally (Boonmee et al., 2017). In 2023, 339 such disasters were registered by the Emergency Event Database EM-DAT with a massive impact on the populations: 81,576 people were killed, 277,382 injured, and 76.52 million people affected. Economic damages from these disasters were estimated at US$295 billion (Beinsure, 2023).

Email addresses: komlanvi-parfait.ametana@u-bordeaux.fr, boris-detienne@u-bordeaux.fr, francois.clautiaux@u-bordeaux.fr, laurent.facq@u-bordeaux.fr, olga.battaia@kedgebs.com, mehdi.amiri-aref@kedgebs.com
Disasters present unique challenges to blood supply chains that differ significantly from normal operations, including sudden demand surges, infrastructure damage, donor availability issues, and the perishability of blood products (Van Denakker et al., 2023). In the context of disasters like earthquakes, the blood supply chain plays a crucial role, and studying it is essential for several reasons. Primarily, a resilient and well-coordinated blood supply chain can save lives and mitigate the impact of disasters, as the demand for blood often surpasses regular stock due to injuries and trauma. To meet this critical need, optimizing resource allocation is key—such as determining optimal locations for blood collection centers and medical facilities, developing models to balance supply and demand while minimizing wastage, and creating strategies for efficient transportation and distribution in damaged areas. Additionally, blood products have specific characteristics that require specialized management. Analyzing past disasters and developing blood supply chain models also enhances overall disaster preparedness (see Asadpour et al. (2022) for more details).

In this study, we focus on establishing optimal locations for blood collection facilities to minimize unmet demand. We employ a two-stage, multi-period model to address the problem, considering hospital blood demand as the uncertain parameter. From a mathematical perspective, demand uncertainty for rare events like earthquakes is often modeled through a discrete set of blood demand scenarios (Jabbarzadeh et al., 2014). However, questions remain about the suitability of optimization model under uncertainty in life-or-death situations. Traditional facility location problems typically prioritize economic performance indicators, but such criteria may not be ethically appropriate in humanitarian contexts. On the other hand, their proximity and accessibility to the demand points can ensure more efficient disaster management operations (UNOCHA, 2018).

To mitigate the potential catastrophic loss of life, we propose an optimization model that integrates blood collection operations into strategic facility location decisions, aiming to minimize the number of people left unsaved in affected areas. This study seeks to provide blood collection facility solutions that aid decision-makers in navigating the uncertainties of post-disaster scenarios. Additionally, this paper contributes to the broader discussion on ethically appropriate optimization in disaster management by comparing various risk-averse objective functions and analyzing their impact on the estimated number of lives saved. Key contributions include:

* A comprehensive evaluation of risk-averse models for the blood supply chain, utilizing mean-CVAR and worst-case risk measures to address uncertainties.

* Development of a novel mathematical model prioritizing humanitarian goals while balancing supply chain costs.

* Implementation of advanced decomposition techniques (Benders Decomposition and Column-and-Constraint Generation) to solve large-scale instances efficiently.

* Comparison of the efficiency of the results obtained with coarse and fine uncertainty approximations in the studied context.

* Managerial insights on the budget management and on the importance of the preparation phase for practitioners.

These contributions aim to enhance emergency blood supply chain management, integrating humanitarian priorities while maintaining economic and operational resilience.

We applied our proposed mathematical model to the city of Paris, France, as a hypothetical case study for a potential future catastrophic event, highlighting the critical need for an immediate and effective humanitarian response plan in the region. The analysis identifies population centers and their respective populations, along with currently available hospitals, as the foundation of the numerical study. Population data is utilized to estimate potential blood demand and supply, while hospital locations serve as clustering reference points. The results provide valuable insights for policymakers, healthcare professionals, and
logistics experts, offering guidance to enhance blood supply chain management strategies in disaster
preparedness and response efforts.

Our study presentation is organized as follows: We analyze the relevant literature in section 2. In section 3, we formally introduce the considered optimization problem. In section 4, we present different mathematical models for risk aversion evaluated in this study. We propose advanced methods for solving these large MILP models in section 5.

We describe the methodology used to compare the models on an out-of-sample set in section 6 and the procedure used to generate the data in section 6.2. The results are analyzed in section 7. A conclusion of the study is presented in section 8.

## 2 Literature review

The specificity of the blood supply chain lies in its complexity and sensitivity, which result from the critical and vital nature of blood products. Each stage of this chain, from the blood collection from the donor to its transfusion to the recipient, requires careful management to respect strict regulations designed to ensure blood products' safety, quality, and traceability. Its sensitivity also stems from its dependence on voluntary donors and the need to maintain sufficient stocks to meet demand, especially in the case of disasters triggered by natural hazards.

### 2.1 Blood supply chain management

In addressing the dynamic challenges of the blood supply chain during crises such as disasters triggered by natural hazards and pandemics, it becomes imperative to ensure its flexibility and responsiveness, as highlighted by Tirkolaee et al. (2023). In such circumstances, the demand for blood products surges, requiring swift mobilization of resources and efficient coordination with humanitarian organizations, hospitals, and emergency services. In practice, the management of the blood supply chain is a complex process including several interdependent planning tasks (Meneses et al., 2023) which have to be addressed at strategic, tactical, and operational decision-making levels (Pirabán et al., 2019).

At the strategic level, the most studied problem concerns strategic network design for the blood supply chain (Chaiwuttisak et al., 2016; Bruno et al., 2019). At the tactical level, inventory management and planning are the most frequently considered in the literature (Hemmelmayr et al., 2010; Zahiri et al., 2018; Rajendran and Ravi Ravindran, 2019; Samani and Hosseini-Motlagh, 2019; Ahmadimanesh et al., 2020). Shih and Rajendran (2020) developed inventory models for the blood supply chain including hospitals and blood centers taking into account the platelet demand and supply uncertainty while considering two types of demand at the blood center: regular demand from the hospitals and the emergency demand when hospitals experience a shortage. Zhao et al. (2021) compared centralized and decentralized strategies for inventory management in the blood supply chain amid fluctuations in supply and demand. Finally, at the operational level, the response to the disaster is managed. A more detailed description of different planning levels and corresponding modeling approaches can be found in two recent comprehensive reviews of optimization models for blood supply chain management published by Pirabán et al. (2019) and Meneses et al. (2023). In the following, we focus on the studies considering the uncertainty in the design of the blood supply chain.

### 2.2 Blood supply chain under uncertainty

Jabbarzadeh et al. (2014) pioneered work on the emergency blood supply chain, studying the case of an earthquake. In their work, they used a *solution robustness* and *model robustness* approach, which is a way of looking at robustness introduced by Mulvey and Vanderbei (1995). A solution is said to be *solution robust* if the value of the objective function remains close to the optimal value associated with
each scenario (as obtained when solving the deterministic model with one scenario). A solution is said
to be *model robust* if feasible for any scenario. The authors proposed a two-stage robust model where
the decisions were related to determining the number and location of permanent and temporary blood
collection facilities, the allocation of these facilities to blood donors and the blood inventory levels at
each period. The planning horizon was divided into several periods characterized by its own decisions.
The objective was to minimize the expectation of the network cost plus a factor times the network cost
variance corresponding to the *solution robustness* described below. Another objective was to minimize
a factor (penalty) multiplied by the expectation of the unsatisfied demand corresponding to the *model
robustness*. The authors applied their model to an earthquake in Iran.

Fahimnia et al. (2017) presented a stochastic mathematical model for designing a supply chain to
minimize total costs and maximize suppliers’ social and environmental scores in the face of disruptions.
The decisions concerned the determination of the number of permanent and temporary blood collection
facilities, the quantity of blood to be collected and transported, the blood inventory levels, and the
quantity of blood transported between facilities. Uncertain parameters in the study included the demand
for blood, the supply of donors, and the fixed and operational costs associated with the blood supply
chain.

Zahiri and Pishvaee (2017) developed a bi-objective mathematical programming model for optimizing
the blood supply chain network to minimize its total cost and maximize demand satisfaction according
to blood group compatibility. The decisions concerned the location of temporary and permanent blood
collection sites, selecting blood testing laboratories, and assigning donors and blood products to demand
zones with perish time considerations. The mathematical model was based on a robust possibilistic
programming approach to address uncertain input parameters such as moving facilities, locating centers,
establishing labs, transportation costs, blood demand, perishability, and blood group compatibility.

Ghatreh Samani et al. (2018) proposed a tri-objective model using a possibilistic programming ap-
proach while considering uncertainty on all the parameters and the perishability of blood products. They
took into account in their model the establishment of blood collection facilities and regional blood centers,
and the assignment between the different facilities of the network. They manage also the flow of blood
in the network. The objective function is composed of the economic cost, the maximum unmet demand,
and the time span between blood production in regional blood centers and demand nodes.

Salehi et al. (2019) studied a two-stage multi-period stochastic model, considering the uncertainty
about blood demand and the potential transfusion of different blood groups based on medical needs. The
objective was to minimize the total costs associated with the blood supply chain network, considering the
average and the cost variability. The decision-making process integrated strategic and tactical decisions
such as the location of temporary and permanent blood collection facilities, the capacity of temporary
blood collection facilities, assignment of donors to facilities, inventory, and backlogging.

Fazli-Khalaf et al. (2019) considered a stochastic tri-objective model based on possibilistic linear
programming tackles parameters such as blood demand in hospitals, transportation time, and laboratory
testing reliability. The first objective was to minimize the total supply chain costs, the second targeted
the minimization of total transportation time between facilities, and the third aimed to maximize the
total reliability of the tested blood at the laboratories. The decisions in the proposed mathematical model
included the location of temporary and permanent blood collection facilities and testing laboratories, the
distribution of blood units to hospitals, and the selection of transport modes for blood delivery.

Shih and Rajendran (2020) proposed a stochastic mixed-integer linear programming model for the
blood supply chain to manage blood supply under uncertainty in demand and supply. The decisions
included the ordering policy for platelet units by hospitals, the number of platelet units procured by the
blood center, the number of units shipped to hospitals, and the inventory and shortage of platelet units
at hospitals. The objective was to minimize the total cost incurred over the blood supply chain.

Khalilpourazari and Hashemi Doulabi (2023) introduced a two-stage multi-period stochastic model
based on fuzzy programming for designing an emergency blood supply chain network, emphasizing the potential impact of natural disasters like earthquakes. The model considered uncertain parameters related to blood demand, derivatives, and medical transfusions. The decisions included strategic preparedness and operational response, as in the papers cited above, addressing issues like the location of temporary and permanent blood collection facilities, the capacity of temporary blood collection facilities, assignment of donors to facilities, selection of transport modes for blood delivery, procurement, inventory, and backlogging. The objective was to minimize the supply chain's total costs and reduce transportation time.

In the literature, uncertainty has been modeled in two ways. The first is scenario modeling, where several possible scenarios are developed to form a discrete set of uncertainties. The second approach is possibility distribution, a concept used in fuzzy logic and possibilistic modeling to represent uncertainty differently from classical probability distributions. Whereas probability distributions quantify the likelihood of various events occurring, possibility distributions measure the possibility of an event occurring. Models based on a discrete set of uncertainties are two-stage adaptive models with variables that can be adapted to the scenarios. In contrast, models based on possibility distributions are static, deterministic models based on fuzzy programming. To guarantee the feasibility of solutions, they use the random constraint approach coupled with a possibility programming measure that evaluates the certainty degree (Inuiguchi and Ramík, 2000).

Risk measures have been used in objective functions to guide the direction of optimization and to evaluate models. A risk measure is a quantitative indicator that assesses the uncertainty or variability associated with a decision. It enables the estimation of potential losses and assists in management and planning decisions. A measure is risk-neutral if it does not give any particular preference to the level of risk associated with different options. In other words, it considers all sources of risk equally and frequently considering the average value of the objective function over all scenarios. A risk-averse measure reflects an aversion to risk, indicating a preference for less risky situations. A less risky situation refers to a context or circumstances where uncertainty and the probability of loss are minimized. Thus, in the previous studies, primarily risk-neutral measures (mean value of the objective function) were used in models based on a discrete uncertainty set or fuzzy programming, except Zahiri and Pishvaee (2017) and Salehi et al. (2019), which used a risk-averse measure to evaluate the unsatisfied demand and the supply chain cost, respectively. The first minimizes the unsatisfied demand for the highest possible value of uncertainty. The second minimizes the maximum regret, where the regret of a scenario is the difference between the cost of the solution found while considering all scenarios and the cost of the solution of the deterministic model with the considered scenario.

For more on this subject, refer to the work of Asadpour et al. (2022), which presents a comprehensive literature review on the blood supply chain in a disaster context. Their study thoroughly examines the various layers involved in the blood supply chain, the different objectives studied, the solution methods used, and the types of uncertainty considered. They also present the different features that have been analyzed and suggest future research perspectives.

The comparison of our study with the previous research is summarized in Table 1.
<table><thead><tr><th rowspan="2">Articles</th><th colspan="2">Objective</th><th rowspan="2">Risk</th><th rowspan="2">Uncertainty</th><th rowspan="2">Solution Method</th></tr><tr><th>Cost</th><th>Demand</th></tr></thead><tbody><tr><td>Jabbarzadeh et al. (2014)</td><td>Yes</td><td>No</td><td>neutral</td><td>scenarios</td><td>B&B</td></tr><tr><td>Fahimnia et al. (2017)</td><td>Yes</td><td>No</td><td>neutral</td><td>scenarios</td><td>LR</td></tr><tr><td>Zahiri and Pishvaee (2017)</td><td>Yes</td><td>Yes</td><td>averse</td><td>fuzzy number</td><td>B&B</td></tr><tr><td>Ghatreh Samani et al. (2018)</td><td>Yes</td><td>Yes</td><td>averse</td><td>fuzzy number</td><td>B&B</td></tr><tr><td>Salehi et al. (2019)</td><td>Yes</td><td>No</td><td>averse</td><td>scenarios</td><td>B&C</td></tr><tr><td>Fazli-Khalaf et al. (2019)</td><td>Yes</td><td>No</td><td>neutral</td><td>fuzzy numbers</td><td>B&B</td></tr><tr><td>Shih and Rajendran (2020)</td><td>Yes</td><td>No</td><td>neutral</td><td>scenarios</td><td>B&B</td></tr><tr><td>Khalilpourazari and Hashemi Doulabi (2023)</td><td>Yes</td><td>No</td><td>neutral</td><td>fuzzy numbers</td><td>B&B</td></tr><tr><td>This research</td><td>No</td><td>Yes</td><td>averse</td><td>scenarios</td><td>B&B, CCG, Benders</td></tr></tbody></table>

B&B: Branch and Bound; B&C: Branch and Cut; CCG: Column and Constraint Generation; LR: Lagrangian relaxation; He:Heuristic/Meta-heuristic

<caption>Table 1: Comparison of the present and previous studies in Blood Supply Chain</caption>

As can be seen from the literature review, such well-known risk-averse measures as the min-max criterion (worst case optimization) and the combination of Conditional Value at Risk (CVAR) with the mean have not been considered yet in the blood supply chain design. The min-max robust optimization involves optimizing the worst-case objective over an uncertainty set. However, min-max models are often criticized for their over-conservative nature, favoring an excessively cautious approach. In contrast, CVAR, a risk measure widely used in finance, represents the average loss above a certain risk level. In our study, we conduct an original evaluation of these risk-averse measures to advise decision-makers on which one best satisfies the demand for blood in disaster management. Given the original contribution of this analysis within the context of blood supply chains, we subsequently examine the application of risk aversion measures as utilized in other contexts within the literature.

## 2.3 Risk-averse measures

Two risk-averse measures are commonly employed in optimization under uncertainty: the robust worst case (min-max) and mean-CVAR combination of CVAR with the expected value. As with the *solution robustness* and *model robustness* approaches, these risk-averse measures provide solutions that are feasible in any scenario (model robustness) but do not use the same selection criteria in the set of feasible solutions. In particular, the min-max criterion is oriented to minimize the losses in the worst case over an uncertainty set.

Ni et al. (2018) proposed a min-max robust model for optimizing facility location, emergency inventory pre-positioning, and relief delivery operations in a disaster relief network. The model considered uncertainties in demands, usable proportions of pre-positioned inventories, and road link capacities. The objective was to minimize the first- and worst-case second-stage costs among all possible realizations of the uncertain parameters falling into the uncertainty set. A case study of the Yushu earthquake was used to demonstrate the application and advantages of the proposed model.

Chen et al. (2022) proposed a min-max robust optimization framework for energy management in combined cooling, heating, and power systems. The authors established an integrated demand response (IDR) model for day-ahead scheduling and intra-day real-time regulation, focusing on uncertainties in renewable energy output and electric, cooling, and heating load. The two-stage robust scheduling optimization model incorporated the cost and constraints of IDR resources, aiming to enhance the system's resilience to uncertain risks and improve the economy and self-sufficiency of micro-energy grids.

Najafi et al. (2022) presented a novel linear max-min-max robust optimization tool for operators in multi-energy systems. This tool accounted for electricity market prices and wind generation uncertainties, incorporating a power-to-gas storage system. The study introduced an uncertainty budget model, enhancing system robustness against various forecasting uncertainties. The study's primary objective was to minimize the total operational cost of procuring energy carriers while meeting short-term demands through a robust optimization model. This approach was designed to ensure feasibility under various uncertainty scenarios and provide optimal solutions for worst-case parameter realizations.
Zhang et al. (2023) addressed optimizing relief kit assembly and distribution in post-disaster scenarios using a min-max robust model. The key objectives were to minimize total costs and maximize demand satisfaction, considering uncertainties in both demand and travel time. The proposed model encompassed facility location and relief kit assembly in the first stage, followed by relief kit distribution in the second stage. The study validated the model through computational experiments and a case study based on earthquakes in Yunnan Province, demonstrating its effectiveness in achieving cost efficiency and meeting demand in post-disaster relief operations.

Nevertheless, min-max models frequently draw criticism for their excessively conservative nature, often adopting an overly cautious approach. CVAR exhibits less conservatism than the worst-case approach, accounting for risks in extreme scenarios. Widely applied across various studies, CVAR proves to be more pertinent than conventional risk measures like expectation in certain contexts.

Noyan et al. (2022) studied risk-averse two-stage stochastic programming with CVAR as the risk measure for disaster management and developed decomposition algorithms for solving such problems. The model focused on determining response facility locations and inventory levels of relief supplies amidst demand and damage level uncertainties.

Yu et al. (2017) addressed an uncapacitated facility location problem incorporating random facility disruptions with independent and correlated disruption scenarios. The objective was to minimize the expected costs, including setup, day-to-day transportation, and penalties. CVAR and absolute-semi-deviation risk measures were used to express the associated risks. As a result, the model was designed to control transportation cost risks under facility disruptions while optimizing facility location and customer assignments. Numerical results showcased the superiority of risk-averse models in enhancing reliability over classic risk-neutral counterparts.

The CVAR risk measure can be found in many other works on various issues (see, for example Pisciella et al. (2016); Bushaj et al. (2022); Afsari et al. (2024)).

In this study, we opt for scenario-based uncertainty modeling and implement the CVAR and worst-case risk measures for the supply chain design problem described in the next section. As shown in Table 2, our contribution involves integrating the routing problem with the specificities of the blood supply chain (donor supply) while focusing on the humanitarian objective. We also used decomposition algorithms to improve the handling of large-scale instances.

<table><caption>Table 2: Comparison of the present and previous studies that used risk-averse measure</caption><thead><tr><th>Articles</th><th>Objective Eco.</th><th>Objective Hum.</th><th>Worst Case</th><th>CVAR</th><th>Disaster</th><th>BSC</th><th>Facility</th><th>Routing</th><th>Solution Method</th></tr></thead><tbody><tr><td>Ni et al. (2018)</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>Benders</td></tr><tr><td>Chen et al. (2022)</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>No</td><td>No</td><td>No</td><td>No</td><td>CCG</td></tr><tr><td>Najafi et al. (2022)</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>No</td><td>No</td><td>No</td><td>No</td><td>B&B</td></tr><tr><td>Zhang et al. (2023)</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>B&B</td></tr><tr><td>Noyan et al. (2022)</td><td>Yes</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>B&C</td></tr><tr><td>Yu et al. (2017)</td><td>Yes</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td><td>B&C, LD</td></tr><tr><td>This research</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>B&B, CCG, Benders</td></tr></tbody><tfoot><tr><td colspan="10">B&B : Branch and Bound; B&C : Branch and Cut; CCG : Column and Constraint Generation; LD : Lagrangian decomposition</td></tr></tfoot></table>

## 3 Problem description

This work focuses on developing a strategic plan for managing blood supply-chain logistics during natural disasters, particularly earthquakes. The primary objective of this supply chain is to transport blood from donors to hospitals where patients require blood transfusions. Blood is collected from two types of facilities: fixed collection centers and mobile collection facilities. Fixed collection centers operate within a specific location and collect blood from nearby donors, while mobile collection facilities can be moved to
different places during the response phase. Since the number of patients requiring treatment is unknown until the disaster occurs, we employ a two-stage approach. In the preparedness phase, we set up permanent blood collection centers in specific locations and created initial blood stocks in hospitals. In the response phase, hospitals receive patients who require treatment on an ongoing basis. We decide where to dispatch mobile collection facilities and the amount of blood each facility delivers to each hospital for each period.

In the remainder, we make the following assumptions:

* We consider a two-stage setting where the uncertain parameters related to all periods are revealed simultaneously. Lead time for the order processing is assumed to be negligible.

* We do not consider individual donors but donor groups, each associated with a possible location and a deterministic blood supply. A single blood type is considered.

* We assume that hospitals and blood collection centers have limited storage capacity. The perishable aspect of blood is not taken into account.

* The Health Insurance Portability and Accountability Act (HIPAA) regulations imply that each hospital can receive blood only from a designated blood center and cannot share or procure blood from other hospitals.

To define the problem, we model the time horizon $\mathcal{T}$ as a set of $T+1$ discrete periods $\{0, \dots, T\}$, where the preparedness phase occurs at period 0, and the response phase takes place during $\mathcal{T}^* := \mathcal{T}\setminus\{0\}$. We denote by $\mathcal{H}$ the set of hospitals where patients are treated and by $\mathcal{I}$ and $\mathcal{J}$ the set of donor groups and the set of potential locations for fixed blood collection facilities, respectively. We also define $\mathcal{E} \subseteq \mathcal{I} \times \mathcal{J}$ as the set of pairs $(i, j)$ such that blood of donor group $i$ can be collected at location $j$ (i.e., the two locations are close enough).

During the preparedness phase, two types of costs are considered: $f_j$, the cost of opening a fixed blood collection site at $j \in \mathcal{J}$, and $\bar{c}_h$, the cost of maintaining each unit of blood in stock at hospital $h \in \mathcal{H}$ during the preparedness.

During the response phase, the cost of collecting one unit of blood from donor group $i \in \mathcal{I}$ at a (fixed or mobile) facility located at $j \in \mathcal{J}$ is $o_{i,j}$. Additionally, $a_{j,h}$ is the transport cost from blood collection site $j \in \mathcal{J}$ to hospital $h \in \mathcal{H}$ and $c_h$ is the cost of storing one unit of blood at hospital $h \in \mathcal{H}$ during one time period. Finally, moving one mobile collection facility from location $j \in \mathcal{J}$ to another location $j' \in \mathcal{J}$, at the beginning of period $t \in \mathcal{T}^*$, costs $b_{j,j',t}$. The deterministic parameters of the problem are summarized in Table 3.
<table><caption>Table 3: Deterministic parameters</caption><thead><tr><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td><i>H</i></td><td>set of hospitals</td></tr><tr><td><i>I</i></td><td>set of donor groups</td></tr><tr><td><i>J</i></td><td>set of potential locations for blood collection facilities</td></tr><tr><td><i>D</i><sub>i,j</sub>, <i>D</i><sub>i,h</sub></td><td>distance between locations <i>i</i> &isin; <i>I</i> and <i>j</i> &isin; <i>J</i> or <i>h</i> &isin; <i>H</i></td></tr><tr><td><i>d</i></td><td>maximum acceptable distance for assigning a donor group to a location</td></tr><tr><td><i>E</i></td><td>set of pairs (<i>i</i>, <i>j</i>) such that <i>j</i> &isin; <i>J</i> can collect blood from <i>i</i> &isin; <i>I</i>, <i>E</i> = {(<i>i</i>, <i>j</i>) &isin; <i>I</i> &times; <i>J</i> : <i>D</i><sub>i,j</sub> &le; <i>d</i>}</td></tr><tr><td><i>T</i></td><td>the latest response period</td></tr><tr><td><i>T</i></td><td>time horizon (<i>T</i> = {0, ..., <i>T</i>}</td></tr><tr><td><i>&rho;</i><sub><i>i</i></sub></td><td>supply of donor group <i>i</i> &isin; <i>I</i></td></tr><tr><td><i>&gamma;</i><sub><i>h</i></sub></td><td>capacity for blood inventory of hospital <i>h</i> &isin; <i>H</i></td></tr><tr><td><i>&gamma;&#773;</i><sub><i>h</i></sub></td><td>maximum blood inventory before the response phase of hospital <i>h</i> &isin; <i>H</i></td></tr><tr><td><i>&beta;</i></td><td>capacity for blood inventory of each fixed collection facility</td></tr><tr><td><i>&delta;</i></td><td>capacity for blood inventory of each mobile collection facility</td></tr><tr><td><i>c</i><sub><i>h</i></sub></td><td>unit blood storage cost at <i>h</i> &isin; <i>H</i> during the response phase</td></tr><tr><td><i>c&#773;</i><sub><i>h</i></sub></td><td>unit blood storage cost at <i>h</i> &isin; <i>H</i> before the response phase</td></tr><tr><td><i>b</i><sub><i>j</i>,<i>j</i>'<sub>,<i>t</i></sub></sub></td><td>cost of moving a mobile collection facility from <i>j</i> &isin; <i>J</i> to <i>j</i>' &isin; <i>J</i> at the beginning of <i>t</i> &isin; <i>T</i><sup>*</sup></td></tr><tr><td><i>f</i><sub><i>j</i></sub></td><td>cost of opening a fixed blood collection facility at <i>j</i> &isin; <i>J</i></td></tr><tr><td><i>a</i><sub><i>j</i>,<i>h</i></sub></td><td>blood unit transport cost from <i>j</i> &isin; <i>J</i> to <i>h</i> &isin; <i>H</i></td></tr><tr><td><i>o</i><sub><i>i</i>,<i>j</i></sub></td><td>blood unit collection cost from <i>i</i> &isin; <i>I</i> at <i>j</i> &isin; <i>J</i></td></tr><tr><td><i>B</i></td><td>available budget</td></tr></tbody></table>

We consider the quantity of blood needed in hospitals to make treatments as uncertain input data of our problem. A realization of the uncertain parameters is a vector $\xi$ in $\mathbb{R}_+^{|H| \times |T^*|}$, where $\xi_{h,t}$ is the quantity of blood needed during period $t \in T^*$ to treat patients at hospital $h \in H$ under realisation $\xi$. Let $\Xi$ be the set of all possible realizations for the uncertain parameters. The uncertain parameters of the problem are summarized in Table 4.

<table><caption>Table 4: Uncertain data</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td>$\Xi$</td><td>Uncertainty set</td></tr><tr><td>$\xi \in \Xi$</td><td>realisation of the uncertainty</td></tr><tr><td>$\xi_{h,t}$</td><td>blood demand for hospital $h$ at period $t$ in uncertainty realization $\xi$</td></tr></tbody></table>

We will now introduce the problem we are considering using the abovementioned notations. In the preparedness phase, we decide on the location of the fixed collection facilities and the amount of initial bloodstock in each hospital before the demand for blood is known. The initial stock at each hospital $h$ is limited to $\bar{\gamma}_h$. The total cost of preparedness comprises the cost of the opened fixed collection facilities and the cost of maintaining the initial stock at hospitals.

In the response phase, which begins once the blood demand $\xi$ is revealed, we determine the location of each mobile collection facility at each period, the assignment of donor groups to collection facilities, the quantity of blood to collect from donors, and the quantity of blood to send to each hospital. The capacity $\gamma_h$ of each hospital $h \in H$ to store blood must not be exceeded at any time. Additionally, each
fixed (resp. mobile) facility cannot collect more blood than its capacity $\beta$ (resp. $\delta$). Donor group $i \in \mathcal{I}$
can give blood at any facility located at site $j \in \mathcal{J}$ if $(i,j) \in \mathcal{E}$. The total amount of blood collected
from donor group $i$ during the entire response phase must not exceed its supply $\rho_i$. Finally, at most one
collection facility can be assigned to location $j \in \mathcal{J}$ during a period.

The cost incurred by the response phase is related to the movements of mobile collection facilities,
blood collection and transportation costs, and the maintenance of the stocks. The total cost of the network
and its operation, including preparation and response costs, must not exceed the available budget of $B$.

The objective of this optimization problem is to minimize the unfulfilled blood demand during the response period $\mathcal{T}^*$. This is a stochastic problem: the evaluation of a design solution depends on the realized uncertain parameter $\xi \in \Xi$. The next section discusses different risk measures and the associated mathematical programming models.

# 4 Mathematical programming models

In this section, we discuss several mathematical programming models that are designed to solve the blood supply chain design problem. Specifically, we present a model that minimizes a risk measure, a robust model, and a model from the literature. We begin by presenting in section 4.1 modeling blocks for the preparedness and response phases common to all our models. To account for the stochastic nature of the problem, we then provide in section 4.2 three two-stage models for specific variants of our problem under uncertainty.

## 4.1 Basic modeling blocks

### 4.1.1 The preparedness phase

Decisions in the preparation phase are modeled by variables $\boldsymbol{x}$ and $\bar{\boldsymbol{s}}$:

* for each location $j \in \mathcal{J}$, the binary variable $x_j$ is equal to 1 if a fixed collection facility is built at $j$, and 0 otherwise.

* for each hospital $h \in \mathcal{H}$, $\bar{s}_h \in \mathbb{R}_+$ represents the initial blood supply in $h$. It is limited to $\bar{s}_h \le \bar{\gamma}_h$, $h \in \mathcal{H}$.

In the following, we denote by $\mathcal{X}$ the set of valid decisions from the preparation phase, i.e., $\mathcal{X} = \{(\boldsymbol{x}, \bar{\boldsymbol{s}}) \in \{0, 1\}^{|\mathcal{J}|} \times \mathbb{R}_+^{|\mathcal{H}|} : \bar{s}_h \le \bar{\gamma}_h, \forall h \in \mathcal{H}\}$.

### 4.1.2 The response phase

The response phase problem is parameterized by the decisions $(\boldsymbol{x}, \bar{\boldsymbol{s}})$ made in the preparedness phase and by the random realization $\boldsymbol{\xi} \in \Xi$. We denote the value of the associated optimal response as $Q(\boldsymbol{x}, \bar{\boldsymbol{s}}, \boldsymbol{\xi})$. To manage the mobile collection facilities and blood flow in the response phase, we use four types of variables:

* The binary variable $y_{j,j',t}$ is equal to 1 if and only if a mobile blood collection facility located at $j \in \mathcal{J}$ at the end of period $t-1$ moves to location $j' \in \mathcal{J}$ and is ready to serve at period $t > 1$. For the special case of $t=1$, we set $y_{j,j',1} = 0$ for $j' \ne j$, while variable $y_{j,j,1}$ is equal to 1 if and only if a mobile facility is present in $j$ at the first period. Only those facilities present at some location during the first time can be used in the rest of the response phase.

* The real variable $q_{i,j,h,t} \in \mathbb{R}_+$ defines the amount of blood to be collected at time $t$ from donor group $i \in \mathcal{I}$ at location $j \in \mathcal{J}$, for $(i,j) \in \mathcal{E}$, to supply hospital $h \in \mathcal{H}$.

* The real variable $s_{h,t} \in \mathbb{R}_+$ defines the amount of blood stored in hospital $h \in \mathcal{H}$ in period $t$.

* The real variable $\theta_{h,t} \in \mathbb{R}_+$ defines the unsatisfied blood demand for hospital $h \in \mathcal{H}$ in period $t$.
Given a first-stage solution $(\boldsymbol{x}, \bar{\boldsymbol{s}}) \in \mathcal{X}$ and a specific $\boldsymbol{\xi} \in \Xi$, the response phase problem is represented by equations (1)-(13).

$$Q(\boldsymbol{x}, \bar{\boldsymbol{s}}, \boldsymbol{\xi}) := \min \sum_{t \in \mathcal{T}^*} \sum_{h \in \mathcal{H}} \theta_{h,t} \quad (1)$$

The objective function (1) minimizes the unsatisfied demand during the response phase, where $\theta_{h,t}$ represents the blood shortage at hospital $h \in \mathcal{H}$ in period $t \in \mathcal{T}^*$.

$$\sum_{j \in \mathcal{J}} f_j x_j + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h + \sum_{t \in \mathcal{T}^*} \left( \sum_{j \in \mathcal{J}} \sum_{j' \in \mathcal{J}} b_{j,j',t} y_{j,j',t} + \sum_{h \in \mathcal{H}} (c_h s_{h,t} + \sum_{(i,j) \in \mathcal{E}} (a_{j,h} + o_{i,j}) q_{i,j,h,t}) \right) \le B \quad (2)$$

Constraint (2) ensures that the total network costs (facility installation, initial storage, mobile facility movement, hospital storage, and blood transportation) do not exceed the available budget $B$.

$$x_j + \sum_{j' \in \mathcal{J}} y_{j',j,t} \le 1 \quad \forall j \in \mathcal{J}, t \in \mathcal{T}^* \quad (3)$$

Constraints (3) specify that only one collection facility (fixed or mobile) can be present at a given location at a given time.

$$\sum_{j' \in \mathcal{J}} y_{j',j,t-1} = \sum_{j' \in \mathcal{J}} y_{j,j',t} \quad \forall j \in \mathcal{J}, t \in \mathcal{T}^* \setminus \{1\} \quad (4)$$

Constraints (4) ensure the continuity of mobile collection facility movement: a mobile facility can leave location $j$ at the end of period $t$ only if it was present there during that period.

$$\sum_{(i,j) \in \mathcal{E}} \sum_{h \in \mathcal{H}} q_{i,j,h,t} \le \beta x_j + \sum_{j' \in \mathcal{J}} \delta y_{j',j,t} \quad \forall j \in \mathcal{J}, t \in \mathcal{T}^* \quad (5)$$

Constraints (5) limit the amount of blood collected at a facility in each period to its capacity, which depends on the type of facility (fixed or mobile).

$$\sum_{\substack{j \in \mathcal{J} \\ (i,j) \in \mathcal{E}}} \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} q_{i,j,h,t} \le \rho_i \quad \forall i \in \mathcal{I} \quad (6)$$

Constraints (6) ensure that the maximum blood supply for each donor group is not exceeded.

$$\sum_{(i,j) \in \mathcal{E}} q_{i,j,h,1} + \bar{s}_h + \theta_{h,1} = s_{h,1} + \xi_{h,1} \quad \forall h \in \mathcal{H} \quad (7)$$

Constraints (7) guarantee the conservation of blood flows in hospitals for the first period.

$$\sum_{(i,j) \in \mathcal{E}} q_{i,j,h,t} + s_{h,t-1} + \theta_{h,t} = s_{h,t} + \xi_{h,t} \quad \forall h \in \mathcal{H}, t \in \mathcal{T}^* \setminus \{1\} \quad (8)$$

Constraints (8) guarantee the conservation of blood flows in hospitals for subsequent periods.

$$s_{h,t} \le \gamma_h \quad \forall h \in \mathcal{H}, t \in \mathcal{T} \quad (9)$$
Constraints (9) ensure that the blood storage capacity in hospitals is not exceeded.

$$y_{j,j',t} \in \{0,1\} \qquad \forall (j,j') \in \mathcal{J} \times \mathcal{J}, t \in \mathcal{T}^* \qquad (10)$$

$$q_{i,j,h,t} \in \mathbb{R}_+ \qquad \forall (i,j) \in \mathcal{E}, h \in \mathcal{H}, t \in \mathcal{T}^* \qquad (11)$$

$$s_{h,t} \in \mathbb{R}_+ \qquad \forall h \in \mathcal{H}, t \in \mathcal{T}^* \qquad (12)$$

$$\theta_{h,t} \in \mathbb{R}_+ \qquad \forall h \in \mathcal{H}, t \in \mathcal{T}^* \qquad (13)$$

The domains of variables in the response phase are defined by constraints (10)-(13).

## 4.2 Discrete scenario-based two-stage models

We first present the general problem considered.

### 4.2.1 General mathematical model

Our focus is on minimizing the impact of the disaster on the population's chances of surviving in good health, given a financial constraint dictated beforehand by exogenous considerations. Since human lives are at stake, we employ a risk-averse strategy. As mentioned in Section 2, the model developed by Jabbarzadeh et al. (2014) optimizes a linear combination of the average unserved demand (through model robustness) and the average and variance of the economic cost (via solution robustness). Hence, this approach is risk-neutral regarding humanitarian considerations: it does not differentiate between a solution that performs well across all scenarios and one that performs exceptionally in some cases but could result in numerous casualties in others. Also, the risk-averse component of the objective function employed for the economic criterion does not constitute a coherent risk measure in the sense of Artzner et al. (1999). Indeed, it does not satisfy the monotonicity property: a solution with small variability in the cost might be preferred to another solution whose cost is smaller in all scenarios and possesses larger variability.

We elected the conditional value-at-risk as our preferred risk measure to model our aversion to humanitarian risk. Given a probability threshold $\alpha \in (0, 1)$, the conditional value at risk $\alpha$ of a random variable $X$ can be intuitively interpreted as the expected value of $X$ in the $100(1-\alpha)\%$ largest outcomes<fnref n="1" />. Formally, the conditional value-at-risk is defined as:

$$\mathrm{CVAR}_{\alpha}(X) = \inf_{\eta \in \mathbb{R}} \left\{ \eta + \frac{1}{1-\alpha} \mathbb{E}_{\omega}[X - \eta]_+ \right\}.$$

Besides, CVAR possesses nice theoretical properties, as its objective satisfies the axioms of a coherent risk measure. Additionally, it has favorable computational properties, as a linear formulation of the model is readily available. However, optimizing CVAR only would completely discard the most favorable outcomes from the evaluation of preparedness solutions, leading to potentially severely suboptimal strategies in many scenarios. To mitigate this undesirable effect, we combine CVAR with the expected value of the unserved demand over all the scenarios with the help of a real weight $\lambda \in [0, 1]$.

Hence, the general mathematical model can be cast as:

$$(\mathrm{MCVAR}^*) : \min (1 - \lambda) \mathbb{E}_{\xi \in \Xi} [Q((\boldsymbol{x}, \bar{\boldsymbol{s}}), \boldsymbol{\xi})] + \lambda \mathrm{CVAR}_{\alpha, \xi \in \Xi} [Q((\boldsymbol{x}, \bar{\boldsymbol{s}}), \boldsymbol{\xi})] \qquad (14)$$

$$\text{s.t.} \ (\boldsymbol{x}, \bar{\boldsymbol{s}}) \in \mathcal{X} \qquad (15)$$

However, this general model remains intractable in practice. In the following, we propose several approximate models with different balances between the model's accuracy and the computational effort required to solve it. The computationally tractable models are based on sampling a finite family $(\boldsymbol{\xi}^\omega)$:

<footnote n="1">¹This interpretation holds for continuous random variables but can be inaccurate in the case of discrete random variables. For more details, we refer to Rockafellar and Uryasev (2002); Schultz and Tiedemann (2006).</footnote>
$\omega \in \Omega$) of scenarios from $\Xi$, where $\Omega$ is the set of scenario indices. The probability of scenario $\omega \in \Omega$ is denoted by $p_\omega$. Section 4.2.2 presents models based on the straightforward application of this idea. Section 4.2.3 introduces a model with a different objective function, which accommodates a larger number of scenarios and allows the computation of near-optimal solutions to this proxy problem with a larger sample size. Finally, Section 4.2.4 recalls the model proposed in Jabbarzadeh et al. (2014), which we include for comparison purposes.

### 4.2.2 Mean-CVAR models

Model (MCVAR) is the sampling-based approximation of (MCVAR*). We derive a mixed-integer linear programming (MILP) formulation using the classical linear formulation of CVAR, proven valid for discrete random variables by Rockafellar and Uryasev (2002) and independently proven to be valid in the context of two-stage stochastic optimization with mixed-integer recourse by Schultz and Tiedemann (2006) and Noyan (2012). To this end, we introduce new variables $\eta \in \mathbb{R}$ and $v^\omega \in \mathbb{R}_+$ for each $\omega \in \Omega$. Given a risk threshold $\alpha \in (0, 1)$ and a weight $\lambda \in [0, 1]$, we obtain:

$$ (\text{MCVAR}) : \min (1 - \lambda) \sum_{\omega \in \Omega} p_{\omega} \mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \xi^{\omega}) + \lambda \left( \eta + \frac{1}{1 - \alpha} \sum_{\omega \in \Omega} p_{\omega} v^{\omega} \right) \quad (16) $$

$$ \text{s.t.} \ (\mathbf{x}, \bar{\mathbf{s}}) \in \mathcal{X} \qquad (17) $$

$$ \eta \in \mathbb{R} \qquad (18) $$

$$ v^{\omega} \ge \mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \xi^{\omega}) - \eta \qquad \forall \omega \in \Omega \quad (19) $$

$$ v^{\omega} \ge 0 \qquad \forall \omega \in \Omega \quad (20) $$

The objective function (16) minimizes a convex combination of the expected unmet demand and the CVAR$_\alpha$ of the unmet demand, and Constraints (19) and (20) allow the linearization of CVAR. The MILP formulation is obtained by writing $\mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \xi^\omega)$ using a block of second-stage variables and constraints of the form (1)-(13) for each $\omega \in \Omega$.

The size of (MCVAR) grows linearly with $|\Omega|$, and solving (MCVAR) is usually computationally challenging when a large number of scenarios is considered. That is why we introduce another approximation of (MCVAR*) obtained by the linear relaxation of the second-stage problem. Specifically, model (MCVAR) is defined by relaxing the integrality in Constraints (10), in each scenario block. This considerably reduces the number of binary variables in the model and makes it amenable to linear programming-based decomposition methods (see Section 5).

### 4.2.3 Robust worst-case model

As the number of scenarios included in $\Omega$ grows, the (MCVAR) model becomes more challenging to solve. In this section, we propose another model whose size also grows linearly with $|\Omega|$, but tends to be less difficult to solve in practice.

Robust optimization finds solutions that maintain their effectiveness even in adverse situations. In other words, this approach aims to solve models that perform well despite unfavorable scenarios, thereby enhancing the reliability of the solutions. The primary objective is to minimize the sensitivity to disturbances and ensure acceptable performance even in critical circumstances. The robust two-stage model we propose employs a worst-case approach consistent with this paradigm. The objective function follows a min-max-min logic, where the outer min operates on the first-stage variables, the inner min minimizes the value of the recourse objective function, and the max identifies the worst-case scenario for the selected first-stage solution.

$$ \min_{(\mathbf{x}, \bar{\mathbf{s}}) \in \mathcal{X}} \max_{\omega \in \Omega} \mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \xi^{\omega}) \qquad (21) $$

Our worst-case solution approaches are based on the following equivalent model, which yields an MILP
model when $Q(x, \bar{s}, \xi^\omega)$ is replaced with its MILP expression (1)-(13):

$$ (\text{WORST-CASE}) : \min z \qquad (22) $$

$$ \text{s.t. } (x, \bar{s}) \in \mathcal{X} \qquad (23) $$

$$ z \ge Q(x, \bar{s}, \xi^\omega) \qquad \forall \omega \in \Omega \qquad (24) $$

$$ z \in \mathbb{R} \qquad (25) $$

The epigraph variable $z \in \mathbb{R}$ equals the worst-case recourse cost at optimality thanks to constraints (24). The objective (22) is to minimize the maximum value of the possible recourse to the first-stage solution (23).

### 4.2.4 The model of Jabbarzadeh et al. (2014)

We rewrite the model proposed in Jabbarzadeh et al. (2014) using the same variables as defined in Section 4.1, duplicating the recourse variables $s, y, q$ and $\theta$ for each $\omega \in \Omega$. We also use some additional notation and variables: the binary variable $k_{i,j,t}^\omega$ is equal to 1 if the donors of group $i \in \mathcal{I}$ are assigned to the collection facility present in $j \in \mathcal{J}$ in period $t \in \mathcal{T}^*$, under realization $\omega \in \Omega$. As proposed by Yu and Li (2000), the Average Absolute Deviation of the economic cost is computed to assess its variability. They suggest using non-negative variables $\Theta^\omega \in \mathbb{R}^+$, $\omega \in \Omega$ to make this term linear. For the sake of conciseness, we denote by $\psi(x, \bar{s}, y^\omega, k^\omega, s^\omega, q^\omega)$ the cost of the policy $(x, \bar{s}, y^\omega, k^\omega, s^\omega, q^\omega)$ under realization $\xi^\omega$, $\omega \in \Omega$:

$$ \psi(x, \bar{s}, y^\omega, k^\omega, s^\omega, q^\omega) := $$

$$ \sum_{j \in \mathcal{J}} f_j x_j + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h + \sum_{t \in \mathcal{T}^*} \left[ \sum_{j \in \mathcal{J}} \sum_{j' \in \mathcal{J}} b_{j,j',t} y_{j,j',t}^\omega + \sum_{h \in \mathcal{H}} c_h s_{h,t}^\omega + \sum_{h \in \mathcal{H}} \sum_{(i,j) \in \mathcal{E}} (a_{j,h} + o_{i,j}) q_{i,j,h,t}^\omega \right] $$

In constraints (40), $M_{i,j,h,t}^\omega := \min\{\rho_i, \max\{\beta, \delta\}\}$ is a constant parameter that is sufficiently large to ensure that $k_{i,j,t}^\omega$ is equal to one whenever $q_{i,j,h,t}^\omega > 0$ for some $h \in \mathcal{H}$, while not discarding any feasible value of $q_{i,j,h,t}^\omega$. In objective function (26), the first and second terms minimize the mean and variability of the cost, respectively; this corresponds to the concept of *solution robustness*. The third term minimizes the average unsatisfied demand; this corresponds to the idea of *model robustness*. The penalty weight $\Gamma \in \mathbb{R}_+$ for not meeting the demand for blood is used to explore the trade-off between the *solution robustness* and the *model robustness*.

$$ (J) : \min \sum_{\omega \in \Omega} p_{\omega} \left[ \psi(x, \bar{s}, y^{\omega}, k^{\omega}, s^{\omega}, q^{\omega}) - \sum_{\omega' \in \Omega} p_{\omega'} \psi(x, \bar{s}, y^{\omega'}, k^{\omega'}, s^{\omega'}, q^{\omega'}) + 2\Theta^{\omega} \right] \\ + \sum_{\omega \in \Omega} p_{\omega} \psi(x, \bar{s}, y^{\omega}, k^{\omega}, s^{\omega}, q^{\omega}) + \Gamma \sum_{\omega \in \Omega} \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} p_{\omega} \theta_{h,t}^{\omega} \qquad (26) $$

$$ \text{s.t. } x_j + \sum_{j' \in \mathcal{J}} y_{j',j,t}^\omega \le 1 \qquad \forall j \in \mathcal{J}, t \in \mathcal{T}^*, \omega \in \Omega \qquad (27) $$

$$ \sum_{i \in \mathcal{I}} \sum_{h \in \mathcal{H}} q_{i,j,h,t}^\omega \le \beta x_j + \sum_{j' \in \mathcal{J}} \delta y_{j',j,t}^\omega \qquad \forall j \in \mathcal{J}, t \in \mathcal{T}^*, \omega \in \Omega \qquad (28) $$

$$ \sum_{j \in \mathcal{J}} \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} q_{i,j,h,t}^\omega \le \rho_i \qquad \forall i \in \mathcal{I}, \omega \in \Omega \qquad (29) $$

$$ \sum_{i \in \mathcal{I}} \sum_{j \in \mathcal{J}} q_{i,j,h,1}^\omega + \bar{s}_h + \theta_{h,1}^\omega = s_{h,1}^\omega + \xi_{h,1}^\omega \qquad \forall h \in \mathcal{H}, \omega \in \Omega \qquad (30) $$

$$ \sum_{i \in \mathcal{I}} \sum_{j \in \mathcal{J}} q_{i,j,h,t}^\omega + s_{h,t-1}^\omega + \theta_{h,t}^\omega = s_{h,t}^\omega + \xi_{h,t}^\omega \qquad \forall h \in \mathcal{H}, t \in \mathcal{T}^* \setminus \{1\}, \omega \in \Omega \qquad (31) $$

$$ s_{h,t}^\omega \le \gamma_h \qquad \forall h \in \mathcal{H}, t \in \mathcal{T}, \omega \in \Omega \qquad (32) $$

$$ y_{j,j',t}^\omega \in \{0, 1\} \qquad \forall (j, j') \in \mathcal{J} \times \mathcal{J}, t \in \mathcal{T}^*, \omega \in \Omega \qquad (33) $$

$$ q_{i,j,h,t}^\omega \in \mathbb{R}_+ \qquad \forall i \in \mathcal{I}, j \in \mathcal{J}, h \in \mathcal{H}, t \in \mathcal{T}^*, \omega \in \Omega \qquad (34) $$
$$s_{h,t}^{\omega} \in \mathbb{R}_{+} \qquad \forall h \in \mathcal{H}, t \in \mathcal{T}^{*}, \omega \in \Omega \qquad (35)$$

$$\theta_{h,t}^{\omega} \in \mathbb{R}_{+} \qquad \forall h \in \mathcal{H}, t \in \mathcal{T}^{*}, \omega \in \Omega \qquad (36)$$

$$\sum_{j' \in \mathcal{J}} y_{j',j,t-1}^{\omega} \geq \sum_{j' \in \mathcal{J}} y_{j,j',t}^{\omega} \qquad \forall j \in \mathcal{J}, t \in \mathcal{T}^{*} \setminus \{1\}, \omega \in \Omega \qquad (37)$$

$$k_{i,j,t}^{\omega} \leq x_j + \sum_{j' \in \mathcal{J}} y_{j',j,t}^{\omega} \qquad \forall i \in \mathcal{I}, j \in \mathcal{J}, t \in \mathcal{T}^{*}, \omega \in \Omega \qquad (38)$$

$$D_{i,j} k_{i,j,t}^{\omega} \leq d \qquad \forall (i,j) \in \mathcal{E}, t \in \mathcal{T}^{*}, \omega \in \Omega \qquad (39)$$

$$q_{i,j,h,t}^{\omega} \leq M_{i,j,h,t}^{\omega} k_{i,j,t}^{\omega} \qquad \forall i \in \mathcal{I}, j \in \mathcal{J}, h \in \mathcal{H}, t \in \mathcal{T}^{*}, \omega \in \Omega \qquad (40)$$

$$\psi(\mathbf{x}, \bar{\mathbf{s}}, \mathbf{y}^{\omega}, \mathbf{k}^{\omega}, \mathbf{s}^{\omega}, \mathbf{q}^{\omega}) - \sum_{\omega' \in \Omega} p_{\omega'} \psi(\mathbf{x}, \bar{\mathbf{s}}, \mathbf{y}^{\omega'}, \mathbf{k}^{\omega'}, \mathbf{s}^{\omega'}, \mathbf{q}^{\omega'}) + \Theta^{\omega} \geq 0 \qquad \forall \omega \in \Omega \qquad (41)$$

$$k_{i,j,t}^{\omega} \in \{0, 1\} \qquad \forall i \in \mathcal{I}, j \in \mathcal{J}, t \in \mathcal{T}^{*}, \omega \in \Omega \qquad (42)$$

$$\Theta^{\omega} \in \mathbb{R}^{+} \qquad \forall \omega \in \Omega \qquad (43)$$

Apart from the objective function, there are a few noticeable differences between the model of Jabbarzadeh et al. (2014) and the models we developed for this study. First, Constraint (37), which models the path of mobile collection facilities, ensures that a mobile facility leaves one location for another only if it was there in the previous period. Using an inequality allows for solutions where the number of mobile collection facilities decreases over time. Considering that the cost of keeping a mobile facility in the same location for two consecutive periods is zero, in our models, we have adopted a flow conservation constraint expressed by Constraint (4) for the number of mobile collection facilities.

Second, Constraints (38) ensure the assignment of donors only to existing facilities, while Constraints (39) ensure that each donor group can only be assigned to facilities within an acceptable distance. Constraints (40) ensure that no blood is collected at a location with no donor group assigned. We have not included these constraints in our models since Constraints (5) and (6) linearly imply them.

Third, no variable related to donor groups and locations that are out of the acceptable range is present in our models thanks to the definition of $\mathcal{E}$. Note that modern MILP solvers automatically remove such variables from model (J) during the presolve procedure.

Fourth, we do not need to use Constraint (41) resulting from the linearization of the Average Absolute Deviation in the objective function, nor Constraints (42) and (43) defining the domain of the auxiliary variables. On the opposite, the model of Jabbarzadeh et al. (2014) does not use the budget constraint (2).

# 5 Solution Methods

One of the main challenges in studying blood supply chains, as noted by Ghatreh Samani et al. (2018), is solving large-scale problems where branch-and-bound algorithms often reach their limitations. Decomposition methods are techniques that iteratively solve relaxations of the original problem while ensuring optimality. They involve handling a master problem (a relaxation of the main problem) and subproblems that help improve this relaxation. For our problem, we propose two decomposition methods: the Benders algorithm to solve model (MCVAR) and the CCG algorithm for model (WORST-CASE).

## 5.1 Benders Decomposition

In this section, we present a Benders decomposition to solve model (MCVAR). This solution algorithm has been independently elaborated in Benders (1962) and Slyke and Wets (1969). The interested reader can refer to Rahmaniani et al. (2017) for a review of the method and its possible improvements. The general principle is to iteratively build an outer approximation of the second-stage value function as
a function of the first-stage variables. The second-stage value function is convex piecewise linear and separable according to the scenario indices. The approach takes advantage of these properties to design a linear programming representation of the function as the sum of the value functions of each scenario, each represented through its subgradients.

Given a first-stage solution $(\boldsymbol{x}^*, \bar{\boldsymbol{s}}^*, \eta^*)$ of $(\overline{\text{MCVA}\bar{\text{R}}})$, the second-stage value function is equal to $\sum_{\omega \in \Omega} p_{\omega} z_{\omega}$, where $z_{\omega}$ is the specific contribution of scenario $\omega \in \Omega$ and is the optimal value of the Benders separation subproblem:

$$ (SP_{\omega}^{B}(\boldsymbol{x}^*, \bar{\boldsymbol{s}}^*, \eta^*)) : \min (1 - \lambda) \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \theta_{h,t} + \frac{\lambda}{1 - \alpha} v \quad (44) $$

s.t. (2) – (9) with $\boldsymbol{x}^*$, $\bar{\boldsymbol{s}}^*$ and $\boldsymbol{\xi}^*$

$$ v - \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \theta_{h,t} \geq -\eta^* \quad (45) $$

$$ v \geq 0 \quad (46) $$

$$ y_{j,j',t} \in [0, 1] \quad \forall (j, j') \in \mathcal{J} \times \mathcal{J}, t \in \mathcal{T}^* \quad (47) $$

$$ q_{i,j,h,t} \in \mathbb{R}_+ \quad \forall (i, j) \in \mathcal{E}, h \in \mathcal{H}, t \in \mathcal{T}^* \quad (48) $$

$$ s_{h,t} \in \mathbb{R}_+ \quad \forall h \in \mathcal{H}, t \in \mathcal{T}^* \quad (49) $$

$$ \theta_{h,t} \in \mathbb{R}_+ \quad \forall h \in \mathcal{H}, t \in \mathcal{T}^* \quad (50) $$

To enforce the relatively complete recourse property, we add the following induced feasibility constraint to the first-stage set of constraints: $\sum_{j \in \mathcal{J}} f_j x_j + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h \leq B$. This redundant constraint ensures that any solution $(\boldsymbol{x}^*, \bar{\boldsymbol{s}}^*, \eta^*)$ that is feasible for the first stage is also feasible for the second stage. It follows that, for a first-stage feasible $(\boldsymbol{x}^*, \bar{\boldsymbol{s}}^*, \eta^*)$ and scenario $\omega \in \Omega$, the optimal value of $(SP_{\omega}^{B}(\boldsymbol{x}^*, \bar{\boldsymbol{s}}^*, \eta^*))$ equals the optimal value of its dual linear program, leading to the Benders reformulation of $(\overline{\text{MCVA}\bar{\text{R}}})$:

$$ (MP^B) : \min \lambda\eta + \sum_{\omega \in \Omega} p_{\omega} z_{\omega} \quad (51) $$

$$ \text{s.t.} \quad \sum_{j \in \mathcal{J}} f_j x_j + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h \leq B \quad (52) $$

$$ \begin{aligned} z_{\omega} &\geq \pi^{*\omega} \left( \sum_{j \in \mathcal{J}} f_j x_j + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h - B \right) - \nu^{*\omega} \eta + \sum_{j \in \mathcal{J}} \sum_{t \in \mathcal{T}^*} \zeta_{j,t}^{*\omega} (x_j - 1) \\ &\quad - \sum_{j \in \mathcal{J}} \sum_{t \in \mathcal{T}^*} \beta \chi_{j,t}^{*\omega} x_j + \sum_{h \in \mathcal{H}} \Lambda_{h,1}^{*\omega} (\xi_{h,1}^{\omega} - \bar{s}_h) + \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^* - \{1\}} \Lambda_{h,t}^{*\omega} \xi_{h,t}^{\omega} \\ &\quad - \sum_{i \in \mathcal{I}} \psi_i^{*\omega} \rho_i - \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \Upsilon_{h,t}^{*\omega} \gamma_h \\ &\qquad \forall \omega \in \Omega, (\pi^{*\omega}, \zeta^{*\omega}, \chi^{*\omega}, \Lambda^{*\omega}, \nu^{*\omega}, \psi^{*\omega}, \Upsilon^{*\omega}) \in Q_D^{\omega} \quad (53) \end{aligned} $$

$$ (\boldsymbol{x}, \bar{\boldsymbol{s}}) \in \mathcal{X} \quad (54) $$

$$ \eta \in \mathbb{R}_+ \quad (55) $$

$$ z_{\omega} \geq 0 \quad \forall \omega \in \Omega \quad (56) $$

Constraints (52) are the induced feasibility constraints, and Constraints (53) are the so-called Benders optimality cuts. The set $Q_D^{\omega}$ is the projection of the set of extreme points of the dual linear program of $(SP_{\omega}^{B}(\boldsymbol{x}^*, \bar{\boldsymbol{s}}^*, \eta^*))$ on the components associated with non-zero right-hand-side constraints: $\pi^{*\omega}$ and $\zeta^{*\omega}$ are associated with (2) and (3), $\chi^{*\omega}$, $\Lambda^{*\omega}$, $\psi^{*\omega}$, and $\Upsilon^{*\omega}$ are dual values associated with the constraints (6) to (9), respectively, while $\nu^{*\omega}$ is associated with (45). Note that $\Lambda^{*\omega}$, indexed by $t$, is used for both constraints (7) and (8), and that the upper bound constraint on $y$ in Constraint (47) is redundant because of Constraints (3).

The Benders decomposition algorithm solves $(MP^B)$ iteratively, starting from a relaxation where only a subset of the optimality cuts (53) are considered. This relaxation, called the relaxed master program
(RMP$^B$), is progressively refined until we can prove that its optimal solution is also optimal for (MP$^B$).
When a candidate solution $(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*, \mathbf{z}^*)$ of (RMP$^B$) is obtained, its feasibility is checked by solving
the separation problem $(SP_{\omega}^B(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*))$ for each $\omega \in \Omega$. If $z_{\omega}^* \le \text{Opt}(SP_{\omega}^B(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*))$ for all $\omega \in \Omega$,
then $(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*)$ is optimal and the algorithm stops. Otherwise, the relaxed master program is refined
by adding a new optimality cut defined with optimal dual variables of $(SP_{\omega}^B(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*))$. In practice,
vanilla implementations of Benders decomposition often suffer from convergence issues. To improve the
convergence of the solution process and speed up each iteration, we implemented the following acceleration
techniques.

**Initialization of the relaxed master program** Our preliminary experiments show that starting with a good subset of optimality cuts is crucial for the convergence of the Benders algorithm: even small-size instances of (MCVAR) cannot be solved to optimality when starting from an empty set of optimality cuts. We use the initialization strategy proposed in Fischetti et al. (2017). It consists of solving the linear relaxation of (MCVAR), i.e. where the integrality of variables $\mathbf{x}$ are relaxed, using Benders decomposition. Some optimality cuts generated during the process are kept to define the initial relaxed master program (RMP$^B$). A simple version of the in-out stabilization technique for cutting planes Ben-Ameur and Neto (2007) is employed. The detailed algorithm is provided in Appendix A.

**Bunching** The bunching technique (see e.g. Birge and Louveaux (2011)) takes advantage of similarities between scenarios to quickly derive optimality cuts for some of them without solving the associated separation subproblems. In problem (MCVAR), uncertainty only appears in the right-hand side of the constraints. Formally, it follows that the separation problem for scenario $\omega \in \Omega$ can be written in the synthetic form $(SP_{\omega}^B(\mathbf{x}, \mathbf{s}, \eta)) : \min\{g^{\top}\mathbf{r} : W\mathbf{r} \ge \mathbf{d}_{\omega} - \mathbf{T}[\mathbf{x}, \mathbf{s}, \eta], \mathbf{r} \ge \mathbf{0}\}$, and the subproblems of all scenarios share the same elements $g$, $W$ and $\mathbf{T}$. Let us consider a given first-stage solution $(\hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\eta})$, a scenario $\omega \in \Omega$, an optimal basis $\mathbf{B}^*$ of the linear program $(SP_{\omega}^B(\hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\eta}))$ and the associated primal-dual solution pair $(\mathbf{r}^*, \varphi^*)$. Since $\varphi^*$ is dual feasible for $(SP_{\omega}^B(\hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\eta}))$, we have that $W^{\top}\varphi^* \le g$. Matrix $W$ and vector $g$ being independent on the scenario, $\varphi^*$ and $\mathbf{B}^*$ are also dual feasible for the separation subproblems of all scenarios $\omega' \in \Omega$. Hence, if $\mathbf{B}^*$ is also primal feasible for a scenario $\omega' \in \Omega$, then $\varphi^*$ is an optimal solution of the dual of $(SP_{\omega}^B(\hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\eta}))$ which can be used to define an optimality cut associated with scenario $\omega'$. To check the primal feasibility of $\mathbf{B}^*$, one can verify the non-negativity of $(\mathbf{B}^*)^{-1}(\mathbf{d}_{\omega'} - \mathbf{T}[\mathbf{x}, \mathbf{s}, \eta])$ as suggested in Birge and Louveaux (2011). However, modern commercial LP solvers do not provide a straightforward way of retrieving $(\mathbf{B}^*)^{-1}$ or checking the feasibility of $\mathbf{B}^*$ by other simple means. That is why we opted for a simple and very fast heuristic procedure. From the primal optimal solution $\mathbf{r}^*$ of $\omega$, we try to build a primal feasible solution $\mathbf{r}'$ for $\omega'$ using a simple greedy repairing heuristic (described in Appendix A). If $g^{\top}\mathbf{r}' = \varphi^{*\top}(\mathbf{d}_{\omega'} - \mathbf{T}[\mathbf{x}, \mathbf{s}, \eta])$, then $\varphi^*$ is proven optimal for $\omega'$ by the weak linear programming duality theorem. In this case, we do not solve the separation problem $(SP_{\omega}^B(\hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\eta}))$ and define the optimality cut associated with $\omega'$ and first-stage solution $(\hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\eta})$. This simple heuristic reduces the computing time of the separation phase of the Benders algorithm by around 10% on average and up to 30% for some instances.

**Early termination of the separation procedure** At a given iteration of the Benders algorithm, let $(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*, \mathbf{z}^*)$ be an optimal solution of the relaxed master program (RMP$^B$). Then $\lambda\eta^* + \sum_{\omega \in \Omega} p_{\omega} z_{\omega}^*$ is a lower bound on the expected value of the first-stage solution $(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*)$. This lower bound can trivially be improved when some separation subproblems in $\Omega' \subseteq \Omega$ have already been solved to obtain $\underline{v} = \lambda\eta^* + \sum_{\omega \in \Omega \setminus \Omega'} p_{\omega} z_{\omega}^* + \sum_{\omega \in \Omega'} \text{Opt}(SP_{\omega}^B(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*))$. If $\underline{v}$ exceeds a known upper bound $\bar{v}$ on the optimal value of (MP$^B$), then $(\mathbf{x}^*, \bar{\mathbf{s}}^*, \eta^*)$ is proven to be non-optimal, and there is no need to separate the optimality cuts associated with scenarios in $\omega \in \Omega \setminus \Omega'$. Note that this case can occur when implementing the Benders decomposition as a branch-and-cut algorithm (see Appendix A).
## 5.2 Column-and-Constraint Generation algorithm for the robust worst-case

In this section, we present a Column-and-Constraint Generation (CCG) algorithm for solving model (WORST-CASE). The CCG algorithm has been introduced by Zeng and Zhao (2013) to solve robust optimization problems with recourse. It relies on an iterative framework, where a master problem and a subproblem are solved successively. The master problem is a relaxation of the original problem, which includes only the variables and constraints associated with a subset of the scenarios. The subproblem is designed to identify the missing variables and constraints that must be incorporated into the master problem to obtain a feasible (and optimal) solution.

Let $(WC_R(\bar{\Omega}))$ denote the relaxed master program obtained from (WORST-CASE) by replacing the scenario set $\Omega$ with $\bar{\Omega} \subseteq \Omega$. Algorithm 1 starts by initializing $\bar{\Omega}$ with an arbitrary scenario. The main loop at iteration $k$ solves the relaxed master program, yielding a first-stage solution $(x^{k*}, \bar{s}^{k*})$ together with its approximate worst-case cost $z^{k*}$. Then, the separation phase checks whether this approximate cost is exact. To this end, it suffices to solve the second-stage problem (1)-(13) for all scenarios in $\omega \in \Omega \setminus \bar{\Omega}$ and check their optimal values $Q(x^{k*}, \bar{s}^{k*}, \xi^\omega)$, against $z^{k*}$. However, these problems are mixed-integer linear programs whose solutions can be computationally challenging. That is why we first solve their linear relaxation, whose value is denoted by $Q_R(x^{k*}, \bar{s}^{k*}, \xi^\omega)$ in the pseudocode. If scenarios $\omega \in \Omega \setminus \bar{\Omega}$ such that $Q_R(x^{k*}, \bar{s}^{k*}, \xi^\omega) > z^{k*}$ are identified, they are added to $\bar{\Omega}$ to refine the relaxation, and the algorithm proceeds to the next iteration. Otherwise, the exact second-stage values of the scenarios not included in $\bar{\Omega}$ yet are computed. As soon as one scenario $\omega \in \Omega \setminus \bar{\Omega}$ such that $Q(x^{k*}, \bar{s}^{k*}, \xi^\omega) > z^{k*}$ is identified, it is added to $\bar{\Omega}$ and the algorithm proceeds to the next iteration without solving other subproblems. If Constraints (24) are satisfied for all $\omega \in \Omega \setminus \bar{\Omega}$, then $(x^{k*}, \bar{s}^{k*}, z^{k*})$ is both a feasible solution of (WORST-CASE) and an optimal solution of a relaxation, hence it is optimal for (WORST-CASE).

Since all problems solved by Algorithm 1 are (mixed-integer) linear programs, the finite convergence of the algorithm is implied by the finite number of scenarios.

# 6 Computational experiments

To compare the different models, we use the Sample Average Approximation (SAA) method (Ahmed and Shapiro, 2002; Kim et al., 2015), commonly used in stochastic programming. It approximates the expectation of a random function by averaging over a finite number of randomly generated samples. SAA can be leveraged to solve a stochastic optimization problem by reformulating it as a deterministic optimization problem, which captures the uncertainty through sampled scenarios. The decision variables are then optimized based on the sample average objective function, which approximates the actual objective function. SAA is typically used when it is computationally infeasible to solve the original stochastic optimization problem exactly or when the true distribution of the uncertain variables is unknown and can only be approximated through sampling. However, SAA introduces approximation errors due to the finite size of the sample used, and the quality of the approximation depends on the number of scenarios and their representativeness of the actual distribution (Ahmed and Shapiro, 2002; Kim et al., 2015).

We use SAA to approximate the mean-CVAR measure of the unmet demand during the response phase to evaluate the performance of a first-stage decision. To obtain a fair comparison, these scenarios differ from those used to optimize the models. Let $\hat{\Omega}$ be the set of scenarios used to compute this measure (called out-of-sample set hereafter). For a first-stage solution $(x, \bar{s})$ provided by the resolution of a model, we compute the value of the recourse for all scenarios $\hat{\omega} \in \hat{\Omega}$ by solving model (1)-(12). Then, for two parameters $\alpha$ and $\lambda$ such that $0 \le \alpha, \lambda \le 1$, we obtain an evaluation of this solution for the out-of-sample set $\hat{\Omega}$ as follows.

$$ (1 - \lambda) \mathbb{E}_{\hat{\omega}}[Q(x, \bar{s}, \xi^\omega)] + \lambda \text{ CVaR}_\alpha(Q(x, \bar{s}, \xi^\omega)). \quad (57) $$
<caption>**Algorithm 1:** CCG Algorithm</caption>

**Data:** $\Omega$

**Result:** Optimal solution $(\mathbf{x}, \bar{\mathbf{s}})$ of (WORST-CASE)

1. $\bar{\Omega} \leftarrow \{\text{arbitrary } \omega \in \Omega\}$;
2. converged $\leftarrow$ false;
3. $k \leftarrow 0$;
4. **while** $\neg$converged **do**
5. &nbsp;&nbsp;&nbsp;&nbsp;$(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, z^{k*}) \leftarrow \text{Solve } (WC_R(\bar{\Omega}))$;
6. &nbsp;&nbsp;&nbsp;&nbsp;L $\leftarrow$ $\emptyset$;
7. &nbsp;&nbsp;&nbsp;&nbsp;$\hat{\Omega} \leftarrow \Omega \setminus \bar{\Omega}$;
8. &nbsp;&nbsp;&nbsp;&nbsp;**for** $\omega \in \hat{\Omega}$ **do**
9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $Q_R(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \boldsymbol{\xi}^{\omega})$;
10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $Q_R(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \boldsymbol{\xi}^{\omega}) > z^{k*}$ **then**
11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;L $\leftarrow$ L $\cup$ $\{\omega\}$;
12. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
13. &nbsp;&nbsp;&nbsp;&nbsp;**end**
14. &nbsp;&nbsp;&nbsp;&nbsp;**while** L = $\emptyset$ and $\hat{\Omega} \neq \emptyset$ **do**
15. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Choose $\omega \in \hat{\Omega}$;
16. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\hat{\Omega} \leftarrow \hat{\Omega} - \{\omega\}$;
17. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $Q(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \boldsymbol{\xi}^{\omega})$;
18. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $Q(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \boldsymbol{\xi}^{\omega}) > z^{k*}$ **then**
19. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;L $\leftarrow$ L $\cup$ $\{\omega\}$;
20. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
21. &nbsp;&nbsp;&nbsp;&nbsp;**end**
22. &nbsp;&nbsp;&nbsp;&nbsp;**if** L = $\emptyset$ **then**
23. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;converged $\leftarrow$ true;
24. &nbsp;&nbsp;&nbsp;&nbsp;**end**
25. &nbsp;&nbsp;&nbsp;&nbsp;**else**
26. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\bar{\Omega} \leftarrow \bar{\Omega} \cup L$;
27. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k + 1$;
28. &nbsp;&nbsp;&nbsp;&nbsp;**end**
29. **end**

Similarly to Shapiro (2003), we derive a confidence interval of the objective function based on the
out-of-sample evaluation. In this purpose, we define random variable $\mathcal{Y}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi})$ as follows.

$$
\mathcal{Y}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}) = (1 - \lambda)\mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}) + \lambda \left( \text{VaR}_{\alpha}(\mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}})) + \frac{1}{1-\alpha} \left[ \mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}) - \text{VaR}_{\alpha}(\mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}})) \right]_+ \right)
$$

We then compute the 95%-confidence interval for $\mathbb{E}_{\hat{\omega}}[\mathcal{Y}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi})]$:

$$
CI_{95\%} = \mathbb{E}_{\hat{\omega}}[\mathcal{Y}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}})] \pm 1.96 \frac{\sigma_{\hat{\omega}}(\mathcal{Y}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}}))}{\sqrt{|\hat{\Omega}|}}
$$

where $\sigma_{\hat{\omega}}^2(\mathcal{Y}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}}))$ is the observed variance of $\mathcal{Y}$ over $\hat{\Omega}$ and 1.96 is the value of the 2.5 percentile of the normal distribution.
<table><caption>Table 5: Data from Jabbarzadeh et al. (2014) and Fazli-Khalaf et al. (2019)</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>number of hospitals |H|</td><td>2</td></tr><tr><td>number of donor groups |I|</td><td>22</td></tr><tr><td>number of potential locations for collection facilities |J|</td><td>22</td></tr><tr><td>number of time periods T</td><td>2</td></tr><tr><td>supply of each donor group ρ<sub>i</sub></td><td>{166, 399, 198, 543, 500, 145, 195, 238, 99, 191, 181, 151, 174, 305, 402, 181, 156, 246, 154, 214, 102, 81}</td></tr><tr><td>blood storage capacity</td><td></td></tr><tr><td>at each hospital γ<sub>h</sub></td><td>{1000, 2200} blood units</td></tr><tr><td>at each hospital before the response phase γ&#773;<sub>h</sub></td><td>0 blood unit</td></tr><tr><td>at each fixed facility β</td><td>500 blood units</td></tr><tr><td>at each mobile facility δ</td><td>200 blood units</td></tr><tr><td>blood storage cost during the response phase c<sub>h</sub></td><td>1 per blood unit</td></tr><tr><td>blood storage cost before the response phase c&#773;<sub>h</sub></td><td>0</td></tr><tr><td>cost of moving a mobile collection facility b<sub>j,j',t</sub></td><td>322.98 for t = 1 and 1 per kilometer for t > 1</td></tr><tr><td>cost of opening a fixed blood collection facility f<sub>j</sub></td><td>1518.23</td></tr><tr><td>blood unit transport cost a<sub>j,h</sub></td><td>0.02 × D(j,h) per blood unit</td></tr><tr><td>collecting cost for one unit of blood o<sub>i,j</sub></td><td>0.0690567 per blood unit</td></tr><tr><td>probability of each scenario p<sub>ω</sub></td><td>1/18</td></tr></tbody></table>

To use the same comparison baseline, we compare the first-stage solutions found by the models for the same given budget. In our models, this budget value is provided as a constraint. This is not the case for the model of Jabbarzadeh et al. (2014), which does not include the budget constraint. To overcome this difference, we first solve the model of Jabbarzadeh et al. (2014) for a set of values of Γ. Second, we evaluate the first-stage solutions that are feasible within the given budget. Finally, we keep the best-performing solution, i.e., the one with the lowest score of (57). If no solution found satisfies the budget constraint, it is considered impossible to design the network under that budget.

## 6.1 Benchmark data

The studies of Jabbarzadeh et al. (2014) and Fazli-Khalaf et al. (2019) provide data with some possible earthquake scenarios over Tehran, Iran. We consolidate this data in Table 5.

The discrete uncertainty set provided by Jabbarzadeh et al. (2014) contains 18 hospital blood demand scenarios for each period. They use the approach proposed by Tabatabaie et al. (2010) to estimate blood demand based on historical data from the Iranian Blood Transfusion Organization (IBTO). The limited scope of data is not enough for a comprehensive evaluation of the models to assess the impact of the chosen objectives on the estimated number of saved lives to guide the choices of decision-makers. Therefore, we developed a novel, more diversified dataset, presented here below.

## 6.2 Data generation

An instance of the considered problem is characterized by geographic data (locations of population, hospitals, distances) and costs. An instance is parameterized by a tuple $(\mathcal{H}, \mathcal{J}, T)$ where $\mathcal{H}$ is the set of hospitals, $\mathcal{J}$ is the set of potential locations for blood facilities, and $T$ is the number of periods in the response phase. Each instance is associated with a family of five sets of scenarios, which differ by their
sizes, representing uncertainty set $\Omega$ and an out-of-sample set of scenarios $\hat{\Omega}$ for evaluating solutions. We consider equiprobable scenarios in both $\Omega$ and $\hat{\Omega}$. We assume a donor group at each potential location; therefore, $\mathcal{I} = \mathcal{J}$.

We divide the considered geographical area into $|\mathcal{I}|$ groups of population $P_i$, $i \in \mathcal{I}$, and simulate the impact of different scenarios of earthquakes on each group. An initial shake and aftershocks define an earthquake scenario according to the time horizon. Let $\mathcal{M}$ be the set of magnitudes, $\mathcal{G}$, the geographical area, and $\mathcal{T}^*$, the response time horizon. A scenario is, therefore, a sequence of events $(e_t)_{t \in \mathcal{T}^*}$, where an event $e_t$ is defined by a pair $(m_t, g_t) \in \mathcal{M} \times \mathcal{G}$ such as $e_1$ is the initial shake event and aftershocks happen in a radius $d \in \mathbb{R}^+$ around $g_1$. Their geographic coordinates identify elements of $\mathcal{G}$.

The total potential blood supply $\rho_i$ of group $i$ is calculated on the basis suggested (Jabbarzadeh et al., 2014), i.e., 22 units of blood per population of 1000 individuals. We also assume that the impact of an earthquake scenario on a group cannot create a blood demand exceeding 22 units for a population of 1000. The demand of each group is assigned to a specific hospital, meaning that after the earthquake, victims of the same group will be brought to the same hospital. Therefore, after selecting subset $\mathcal{I}$ of locations and subset $\mathcal{H}$ of hospitals, we run a clustering procedure to assign victims of each group $i \in \mathcal{I}$ to a hospital $h \in \mathcal{H}$. We ensure that each cluster contains at least $\mu\%$ of groups.

To evaluate the impact of an earthquake scenario $\omega \in \Omega$ on the population, we denote by $P_{i,t}^{\omega}$ the population of group $i \in \mathcal{I}$ during period $t \in \mathcal{T}$, potentially impacted by the event of scenario $\omega$. According to the assumption made before, $P_{i,0}^{\omega} = \rho_i = 22 \times P_i/1000$, $i \in \mathcal{I}$, and $\omega \in \Omega$, where $P_{i,0}^{\omega}$ is the initial population ($P_{i,0}^{\omega} = P_{i,0}^{\omega'}, \forall \omega, \omega' \in \Omega$). This population will decrease with the occurrence of events in each scenario. For a given event $e_t$, the impact on a population $i \in \mathcal{I}$ is given by function $\phi(e_t, i)$ described in Equation (58) with $d_{max} \ge D(g, i)$, $\forall (g, i) \in \mathcal{G} \times \mathcal{I}$, and $D(g, i)$ is the function providing the distance between the geographical points $g$ and $i$.

$$ \phi(e_t, i) = m_t \left[ 1 - \frac{D^2(g_t, i)}{d_{max}^2} \right] \qquad (58) $$

Therefore, we compute the potentially impacted population of the group $P_{i,t}^{\omega} = P_{i,t-1}^{\omega}(1 - \phi(e_t, i))$. The demand for blood in hospital $\xi_{h,t}^{\omega}$ is then given by $\xi_{h,t}^{\omega} = \sum_{i \in \mathcal{I}_h} \phi(e_t, i) P_{i,t-1}^{\omega}$, where $\mathcal{I}_h$ is the set of groups assigned to hospital $h \in \mathcal{H}$ by clustering.

We use Algorithm 2 to generate the discrete set of scenarios $\Omega$. It is also used to create the out-of-sample set. We assume that the territory can be hit by one of $L$ principal shakes that differ in magnitude and location, which we call initial events. We generate $|\Omega|$ scenarios starting with these initial events by extending each initial event with $|\Omega|/L$ sequences of $|\mathcal{T}| - 1$ events.

In our numerical experiment, we use the city of Paris as a basis for the geographic data. We did a web scraping to obtain the geographical coordinates and population of the districts of Paris and the locations of the hospitals. We generated different instances from these data by selecting several values for the different parameters, resulting in 60 different instances. Then, we separately generated two sets of 20,000 scenarios for each instance, each with $L = 2000$. We constructed the uncertainty sets $\Omega$ of different sizes from the first set of scenarios by randomly selecting scenarios. The second set of 20,000 scenarios is the out-of-sample scenarios $\hat{\Omega}$. For the creation of these new instances, we used the following parameter values: $|\mathcal{H}|$ in $\{1, 2, 3, 4, 5\}$, $|\mathcal{J}|$ in $\{20, 40, 80\}$, $T$ in $\{2, 3, 4, 5\}$, and $|\Omega|$ in $\{20, 50, 100, 200, 300\}$.

Two types of hospitals are considered: small and large. The latter are assumed to be the most important hospitals in the territory, with large blood storage capacities. The former are less important and possess only half of the capacity of large hospitals. Let $n_b$ be the number of large hospitals and $n_s$ be the number of small ones. In our implementation, if $|\mathcal{H}|$ is even then $n_b = n_s$, otherwise $n_b = n_s + 1$. We assume that the total blood storage capacity of the hospitals is equal to the total supply of donors group: $\sum_{h \in \mathcal{H}} \gamma_h = \sum_{i \in \mathcal{I}} \rho_i$. The blood storage capacity of a small hospital $h \in \mathcal{H}$ is calculated as follows.

$$ \gamma_h = \sum_{i \in \mathcal{I}} \frac{\rho_i}{2 * n_b + n_s} $$
<caption>**Algorithm 2:** Generation of scenarios</caption>

**Data:** territory $\mathcal{G}$, set of magnitudes $\mathcal{M}$, time horizon $\mathcal{T}$, the number of locations $\mathcal{I}$, number of hospitals $\mathcal{H}$, the number of initiating events $L$, the set of indices of scenarios $\Omega$, population of each group $(P_i)_{i \in \mathcal{I}}$, $(\mathcal{I}_h)_{h \in \mathcal{H}}$, impact estimation function $\phi$ defined in Equation (58), maximum distance between a principal shake point and its aftershocks points $d$

**Result:** $(\xi_{h,t}^{\omega})_{t \in \mathcal{T}^*, \omega \in \Omega}$, blood demand scenarios

1 /* initialisation */;
2 for $(\omega, i) \in \Omega \times \mathcal{I}$ do
3 $P_{i,0}^{\omega} = 22P_i/1000$
4 end
5 /* generation of scenarios*/;
6 for $s \in \{1, \dots, L\}$ do
7 $e_1 \leftarrow$ random choice of an initial event;
8 for $\omega \in \{\frac{|\Omega|}{L}(s-1) + 1, \dots, \frac{|\Omega|}{L}s\}$ do
9 for $t \in \mathcal{T}^*$ do
10 for $h \in \mathcal{H}$ do
11 $\xi_{h,t}^{\omega} \leftarrow 0$;
12 for $i \in \mathcal{I}_h$ do
13 $\xi_{h,t}^{\omega} \leftarrow \xi_{h,t}^{\omega} + \phi(e_t, i)P_{i,t-1}^{\omega}$;
14 $P_{i,t}^{\omega} \leftarrow (1 - \phi(e_t, i))P_{i,t-1}^{\omega}$;
15 end
16 end
17 /* choosing of the next replicas event of period $t+1$*/;
18 $e_{t+1} \leftarrow$ random choice of an event such that $D(g_1, g_t) \le d$;
19 end
20 end
21 end

For the permanent blood collection facilities, we assume that their capacity depends on the population's size in the area. Moreover, we note that in Fazli-Khalaf et al. (2019), the capacity of a large hospital was 4.4 times larger than the capacity of a permanent blood collection facility. As the capacity of the hospital depends already on the population size, we use the same ratio to compute the capacity of the permanent blood collection facility $\beta$. Let us compare the ratio cost/capacity of permanent and mobile blood collection facilities. We note that for one unit of blood, using a permanent facility is 1.88 times more expensive than using a mobile facility. Taking this ratio into account, we compute the establishment cost of a permanent facility as $1.88 * \beta/200$. Other cost parameters are kept the same as in (Jabbarzadeh et al., 2014; Fazli-Khalaf et al., 2019). The values of parameters are summarized in Table 6.
<table><caption>Table 6: Deterministic data</caption><thead><tr><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>|H|</td><td>{1, 2, 3, 4, 5}</td></tr><tr><td>|I|</td><td>{20, 40, 80}</td></tr><tr><td>|J|</td><td>{20, 40, 80}</td></tr><tr><td>T</td><td>{2, 3, 4, 5}</td></tr><tr><td>ρ<sub>i</sub></td><td>22 * P<sub>i</sub>/1000</td></tr><tr><td>γ<sub>h</sub></td><td>{∑<sub>i∈I</sub> ρ<sub>i</sub>/(2 * n<sub>b</sub> + n<sub>s</sub>), 2 * ∑<sub>i∈I</sub> ρ<sub>i</sub>/(2 * n<sub>b</sub> + n<sub>s</sub>)} blood units</td></tr><tr><td>β</td><td>∑<sub>i∈I</sub> (2 * ρ<sub>i</sub>)/(4 * (2 * n<sub>b</sub> + n<sub>s</sub>)) blood units</td></tr><tr><td>δ</td><td>200 blood units</td></tr><tr><td>c<sub>h</sub></td><td>1 per blood unit</td></tr><tr><td>b<sub>j,j',t</sub></td><td>322.98 for t = 1 and for other t, 1 per kilometer</td></tr><tr><td>f<sub>j</sub></td><td>1.88 * β/200</td></tr><tr><td>a<sub>j,h</sub></td><td>0.02 per kilometer and per blood unit</td></tr><tr><td>o<sub>i,j</sub></td><td>0.0690567 per blood unit</td></tr><tr><td>B</td><td>{10%, 20%, 30%, 40%, 50%} of the estimated upper limit of the budget to be used</td></tr><tr><td>Γ</td><td>{1, 5, 10, 15, 25, 35, 45} times the estimated lower bound</td></tr></tbody></table>

The problem to solve can be represented by a triplet (|H|, |I|, |J|, |T|, Ω, B) where |H|, |I|, |J|, |T| is an instance, Ω the uncertainty set and B ∈ ℝ_+ the budget for the preparedness and the response to the earthquake.

$$B^{sup} = \sum_{j \in \mathcal{J}} f_j + \beta(o_{i,j} + a_{j,h*}) |\mathcal{T}^*| \quad (59)$$

where h* ∈ H with D(j, h*) ≥ D(j, h) ∀h ∈ H and o<sub>i,j</sub> = o<sub>i',j</sub> ∀(i, i') ∈ I<sup>2</sup>

To obtain values for the budget parameter, we first calculate an upper bound $B^{sup}$ (59) on the budget that can be required to build the supply chain under the condition that there is a permanent facility at each location that collects its total capacity at each period and delivers it to the most distant hospital. We then take 10%, 20%, 30%, 40%, and 50% of this upper bound to form the set of values for budget B associated with each instance |H|, |I|, |J|, |T|.

The percentages we use have been identified experimentally on instances 1_20_2 and 2_20_2 for which we resolved the model of Jabbarzadeh et al. (2014) with Ω of size 100 and Γ taking values in {1, 10, 50, 70, 100, 160, 200, 500, 1000, 3200}. We analyzed the average returned cost and set the possible budget values to 20,000, 50,000, 100,000, 150,000, 200,000, and 250,000. In a second step, we solved the mean-CVAR model (λ = 0.5 and α = 0.8) on the two instances considered with an Ω of 100 and successively with the different values of budget found above. By calculating the ratio of the budget used in the found solutions to the upper bound of the budget calculated for each instance (430,174 for the first instance and 371,304 for the second), we obtain the following percentages:

* 1_20_2: 4.65% - 11.62% - 23.25% - 34.87% - 46.49% - 58.12%

* 2_20_2: 5.39% - 13.47% - 26.93% - 40.40% - 53.86% - 67.33%

From these results, we have decided to use 10%, 20%, 30%, 40%, and 50% of the upper limit of each instance. It should be noted that for these two cases, with a budget of 250,000, the optimal solution found by the mean-CVAR model satisfies all the demand for blood so that it will be the same for any higher budget.
Concerning the parameter Γ used in the objective function in Jabbarzadeh et al. (2014), we have analyzed the results provided and have derived a rule for generating Γ values to solve the model of Jabbarzadeh et al. (2014) on new instances as shown in (60):

$$ \Gamma_{inf} = \frac{B^{sup}}{\max_{\omega \in \Omega} \left\{ \sum_{t \in T^*} \sum_{h \in H} \xi_{h,t}^{\omega} \right\}} \qquad (60) $$

Taking into account Formula (60), the following seven values of Γ were used for each new instance in numerical tests: $1\Gamma_{inf}$, $5\Gamma_{inf}$, $10\Gamma_{inf}$, $15\Gamma_{inf}$, $25\Gamma_{inf}$, $35\Gamma_{inf}$ and $45\Gamma_{inf}$.

## 6.3 Experimental setting

Models and algorithms were written in C++, and IBM ILOG Cplex 12.10 was employed as a solver for the MILP models. The evaluation process includes training and testing phases. In the training phase, each instance is solved with the tested method with the scenario set Ω and the specified budget. The solution obtained for the preparedness phase (decisions about fixed facilities) is stored and used as input for the testing phase. The testing phase involves solving the recourse problem for the first-stage solution (obtained during training) for each scenario from set $\hat{\Omega}$ and the given budget. For the training phase, we set a time limit of 1 hour, used two cores, and had a memory limit of 16 GB. For the testing phase on out-of-sample scenarios, we set a time limit of 15 minutes and used eight cores.

For each instance, we used five different training sets with the size of the training set Ω taking the following values {20, 50, 100, 200, 300}. The out-of-sample scenario set $\hat{\Omega}$ of each instance contained 20,000 scenarios. The following methods were tested: Jabbarzadeh et al. (2014) model (described in Section 4.2.4) labeled as jab, the robust worst-case model (WORST-CASE)(described in section 4.2.3) labeled as wc, the mean-CVAR model (MCVAR)(described in Section 4.2.2) labeled as mc8 for $\alpha = 0.8$, the model (MCVAR) defined by relaxing the integrality in Constraints (10), labeled as mcr8 for $\alpha = 0.8$, the Benders decomposition algorithm (described in section 5.1) labeled as ben and the CCG algorithm (described in section 5.2) labeled as ccg. Table 7 summarizes the number of runs in the computational experiments.

<table><caption>Table 7: Total number of runs for our computational experiments</caption><thead><tr><th>Model</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th><th>ben</th><th>ccg</th><th>Total</th></tr></thead><tbody><tr><td>Training</td><td>2,100</td><td>1,500</td><td>1,500</td><td>1,500</td><td>1,500</td><td>1,500</td><td>9,600</td></tr><tr><td>Testing</td><td colspan="7">32,700,000</td></tr></tbody></table>

# 7 Results and Analysis

## 7.1 Results for the training phase

First, we report the results obtained during the training phase. We compare the capability of different models to produce feasible solutions and to achieve convergence across various instance sizes. Additionally, we analyze the impact of different instance parameters on the difficulty of solving the tested models. Aggregated results are shown in Table 8. For each model and training set, we report the percentage of instances solved optimally (columns Optimal), the percentage of instances for which a feasible non-optimal solution was found (columns Non-optimal), and the percentage of instances that could not be built due to exceeding the memory limit (columns Memory error). We aggregate all results related to the same number of scenarios (20, 50, 100, 200, and 300). The average results are reported in the last line.
<table><thead><tr><th rowspan="2">|Ω|</th><th colspan="4">Optimal (%)</th><th colspan="4">Non-optimal (%)</th><th colspan="4">Memory error (%)</th></tr><tr><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th></tr></thead><tbody><tr><td>20</td><td>59.28</td><td>44</td><td>20</td><td>52.67</td><td>10.23</td><td>47.33</td><td>68</td><td>47.33</td><td>30.47</td><td>8.66</td><td>12</td><td>0</td></tr><tr><td>50</td><td>36.19</td><td>33</td><td>12.66</td><td>45</td><td>3.33</td><td>42.33</td><td>60.33</td><td>48</td><td>60.47</td><td>24.66</td><td>27</td><td>7</td></tr><tr><td>100</td><td>11.90</td><td>21</td><td>6.33</td><td>34.33</td><td>2.14</td><td>38</td><td>51.33</td><td>41</td><td>85.95</td><td>41</td><td>42.33</td><td>24.67</td></tr><tr><td>200</td><td>0</td><td>11.66</td><td>3</td><td>24.67</td><td>0</td><td>33.33</td><td>38.66</td><td>32.33</td><td>100</td><td>55</td><td>58.33</td><td>43</td></tr><tr><td>300</td><td>0</td><td>6.66</td><td>1</td><td>19</td><td>0</td><td>28.33</td><td>34</td><td>29.67</td><td>100</td><td>65</td><td>65</td><td>51.33</td></tr><tr><td>avg.</td><td>21.47</td><td>23.26</td><td>8.60</td><td>35.13</td><td>3.14</td><td>37.86</td><td>50.46</td><td>39.67</td><td>75.38</td><td>38.86</td><td>40.93</td><td>25.2</td></tr></tbody></table>

<table><caption>Table 8: Percentage of instances solved optimally, having reached a time limit, and having raised a memory error by model and training set.</caption><thead><tr><th rowspan="2">|Ω|</th><th colspan="2">Optimal (%)</th><th colspan="2">Non-optimal (%)</th><th colspan="2">Memory error (%)</th></tr><tr><th>ben</th><th>ccg</th><th>ben</th><th>ccg</th><th>ben</th><th>ccg</th></tr></thead><tbody><tr><td>20</td><td>91</td><td>50.67</td><td>9</td><td>49.33</td><td>0</td><td>0</td></tr><tr><td>50</td><td>88</td><td>46.67</td><td>12</td><td>53.33</td><td>0</td><td>0</td></tr><tr><td>100</td><td>83.33</td><td>46.33</td><td>16.67</td><td>52.33</td><td>0</td><td>1.33</td></tr><tr><td>200</td><td>80.33</td><td>40</td><td>17.67</td><td>51.33</td><td>2</td><td>8.67</td></tr><tr><td>300</td><td>79.33</td><td>34.67</td><td>19.33</td><td>44.33</td><td>1.34</td><td>21</td></tr><tr><td>avg.</td><td>84.4</td><td>43.67</td><td>14.93</td><td>50.13</td><td>0.67</td><td>6.2</td></tr></tbody></table>

From the obtained results, we can conclude that the number of scenarios significantly impacts the solver's ability to tackle each model. The percentages of Optimal and Non-optimal solutions decrease as the number of scenarios increases, while the percentage of memory errors rises. This is unsurprising since the number of scenarios affects the number of variables and constraints in all models, leading to larger constraint matrices. Additionally, even when the models can be loaded into memory, they contain binary variables indexed by the scenarios (except for mcr8 model), resulting in a larger branch-and-bound tree. For the jab model, the solver could not load any instance with 200 or 300 scenarios. For other models, feasible solutions were found for instances of all sizes, but the percentage of solved instances decreases significantly as the number of scenarios increases. We can notice that when the solver can load model jab, it can generally converge (only 3% of the instances are loaded but not solved optimally).

Table 9 is a disaggregated version of Table 8, where results are split by the number of locations. This table illustrates the impact of the number of locations on the problem's difficulty. The columns have the same meaning as in Table 8.

It can be observed from Table 9 that the difficulty also increases with the number of locations. For 80 locations, the jab model was able to solve only 10% of the 20 scenarios cases, and could not build the models for 50, 100, 200, and 300 scenario cases. On the other hand, the other models solved at least 64% of the 80-location and 20-scenarios instances.

The number of variables is of the order of $O(|\Omega||T||\mathcal{J}|^2|\mathcal{H}|)$ for all models, but the number of constraints in jab model is larger in comparison to other models: it is of the order of $O(|\Omega||T||\mathcal{J}|^2|\mathcal{H}|)$. In contrast, the number of constraints in the other models is of the order of $O(|\Omega||T|(3|\mathcal{J}| + |\mathcal{H}|))$. This explains the fact why jab is more affected by the number of locations.

The analysis of results in Table 9 shows the superiority of decomposition-based methods over mono-lithic approaches in solving these large-scale optimization problems. This comparison, involving four different methods mcr8, ben, wc, and ccg demonstrates the notable advantages in both optimality and robustness offered by decomposition strategies.

First, we observe that ben outperforms mcr8. The former achieves an average optimality rate of 84.4%
<table><thead><tr><th rowspan="2">|J|</th><th rowspan="2">|Ω|</th><th colspan="4">Optimal (%)</th><th colspan="4">Non-optimal (%)</th><th colspan="4">Memory error (%)</th></tr><tr><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th></tr></thead><tbody><tr><td rowspan="5">20</td><td>20</td><td>97.86</td><td>66</td><td>16</td><td>72</td><td>2.14</td><td>34</td><td>84</td><td>28</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>50</td><td>93.57</td><td>63</td><td>16</td><td>67</td><td>6.42</td><td>37</td><td>83</td><td>33</td><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td>100</td><td>35.71</td><td>52</td><td>12</td><td>58</td><td>6.42</td><td>48</td><td>88</td><td>42</td><td>57.85</td><td>0</td><td>0</td><td>0</td></tr><tr><td>200</td><td>0</td><td>33</td><td>7</td><td>53</td><td>0</td><td>67</td><td>93</td><td>47</td><td>100</td><td>0</td><td>0</td><td>0</td></tr><tr><td>300</td><td>0</td><td>20</td><td>2</td><td>41</td><td>0</td><td>70</td><td>88</td><td>59</td><td>100</td><td>10</td><td>10</td><td>0</td></tr><tr><td rowspan="5">40</td><td>20</td><td>80</td><td>38</td><td>17</td><td>36</td><td>18.57</td><td>62</td><td>83</td><td>64</td><td>1.42</td><td>0</td><td>0</td><td>0</td></tr><tr><td>50</td><td>15</td><td>29</td><td>15</td><td>32</td><td>3.57</td><td>71</td><td>83</td><td>68</td><td>81.42</td><td>0</td><td>2</td><td>0</td></tr><tr><td>100</td><td>0</td><td>10</td><td>6</td><td>27</td><td>0</td><td>62</td><td>64</td><td>73</td><td>100</td><td>28</td><td>30</td><td>0</td></tr><tr><td>200</td><td>0</td><td>2</td><td>2</td><td>18</td><td>0</td><td>33</td><td>23</td><td>49</td><td>100</td><td>65</td><td>75</td><td>33</td></tr><tr><td>300</td><td>0</td><td>0</td><td>1</td><td>16</td><td>0</td><td>15</td><td>14</td><td>30</td><td>100</td><td>85</td><td>85</td><td>54</td></tr><tr><td rowspan="5">80</td><td>20</td><td>0</td><td>28</td><td>27</td><td>50</td><td>10</td><td>46</td><td>37</td><td>50</td><td>90</td><td>26</td><td>36</td><td>0</td></tr><tr><td>50</td><td>0</td><td>7</td><td>7</td><td>36</td><td>0</td><td>19</td><td>15</td><td>43</td><td>100</td><td>74</td><td>78</td><td>21</td></tr><tr><td>100</td><td>0</td><td>1</td><td>1</td><td>18</td><td>0</td><td>4</td><td>2</td><td>8</td><td>100</td><td>95</td><td>97</td><td>74</td></tr><tr><td>200</td><td>0</td><td>0</td><td>0</td><td>3</td><td>0</td><td>0</td><td>0</td><td>1</td><td>100</td><td>100</td><td>100</td><td>96</td></tr><tr><td>300</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>100</td><td>100</td><td>100</td><td>100</td></tr></tbody><tfoot><tr><td>avg.</td><td></td><td>21.48</td><td>23.26</td><td>8.60</td><td>35.13</td><td>3.14</td><td>37.86</td><td>50.46</td><td>39.67</td><td>75.38</td><td>38.86</td><td>40.93</td><td>25.2</td></tr></tfoot></table>

<table><caption>Table 9: Percentage of instances solved optimally, having reached a time limit, and having raised a memory error, organized by model, training set, and location.</caption><thead><tr><th rowspan="2">|J|</th><th rowspan="2">|Ω|</th><th colspan="2">Optimal (%)</th><th colspan="2">Non-optimal (%)</th><th colspan="2">Memory error (%)</th></tr><tr><th>ben</th><th>ccg</th><th>ben</th><th>ccg</th><th>ben</th><th>ccg</th></tr></thead><tbody><tr><td rowspan="5">20</td><td>20</td><td>100</td><td>76</td><td>0</td><td>24</td><td>0</td><td>0</td></tr><tr><td>50</td><td>100</td><td>68</td><td>0</td><td>32</td><td>0</td><td>0</td></tr><tr><td>100</td><td>100</td><td>64</td><td>0</td><td>36</td><td>0</td><td>0</td></tr><tr><td>200</td><td>100</td><td>59</td><td>0</td><td>41</td><td>0</td><td>0</td></tr><tr><td>300</td><td>100</td><td>52</td><td>0</td><td>48</td><td>0</td><td>0</td></tr><tr><td rowspan="5">40</td><td>20</td><td>83</td><td>35</td><td>17</td><td>65</td><td>0</td><td>0</td></tr><tr><td>50</td><td>81</td><td>33</td><td>19</td><td>67</td><td>0</td><td>0</td></tr><tr><td>100</td><td>74</td><td>35</td><td>26</td><td>65</td><td>0</td><td>0</td></tr><tr><td>200</td><td>74</td><td>26</td><td>26</td><td>74</td><td>0</td><td>0</td></tr><tr><td>300</td><td>74</td><td>24</td><td>26</td><td>72</td><td>0</td><td>4</td></tr><tr><td rowspan="5">80</td><td>20</td><td>90</td><td>41</td><td>10</td><td>59</td><td>0</td><td>0</td></tr><tr><td>50</td><td>83</td><td>39</td><td>17</td><td>61</td><td>0</td><td>0</td></tr><tr><td>100</td><td>76</td><td>40</td><td>24</td><td>56</td><td>0</td><td>4</td></tr><tr><td>200</td><td>67</td><td>35</td><td>27</td><td>39</td><td>6</td><td>26</td></tr><tr><td>300</td><td>64</td><td>28</td><td>32</td><td>13</td><td>4</td><td>59</td></tr></tbody><tfoot><tr><td>avg.</td><td></td><td>84.4</td><td>43.67</td><td>14.93</td><td>50.13</td><td>0.67</td><td>6.2</td></tr></tfoot></table>
across all tested instances, even reaching a perfect 100% for smaller instances (|J| = 20 and |Ω| ≤ 300). Conversely, mcr8 achieves a significantly lower average optimality rate of 35.13%, with performance deteriorating as the problem size increases. For large instances (|J| = 80 and |Ω| = 300), mcr8 struggles, yielding 0% optimality. Furthermore, mcr8 suffers from memory errors in 25.2% of instances, whereas ben handles memory efficiently, with only 0.67% recorded faults.

Similarly, the second comparison between `wc` and `ccg` highlights the advantages of decomposition. While `wc` achieves a low average optimality rate of 23.26%, `ccg` nearly doubles this, reaching 43.67%. Additionally, `ccg` exhibits more stable performance across varying problem sizes, while `wc` faces significant challenges, with its optimality rate dropping to 1% for larger instances (|J| = 80 and |Ω| = 100). Regarding memory errors, `ccg` is more robust, with lower fault rates than `wc`, which suffers from 100% memory errors for large instances (|J| = 80 and |Ω| ≥ 200).

The analysis also reveals the limitations of the decomposition methods. While `ccg` performs well overall, it occasionally exhibits higher non-optimality rates than `wc` for larger instances. For example, when |J| = 40 and |Ω| = 200 or 300, `ccg` shows non-optimality rates of 74% and 72%, respectively, compared to 33% and 15% for `wc`. This indicates that `ccg` excels for smaller instances but may lose its advantage in finding optimal solutions as the problem size grows. Moreover, while `ccg` encounters fewer memory errors than `wc`, it still experiences significant fault rates for large instances (|J| = 80 and |Ω| = 200 or 300), with rates of 26% and 59%, respectively.

In Table 10, we report the computing time required to solve small to medium instances for the tested models. For each method, we report the average optimality gap indicated by the MILP solver after one hour of computing time (columns Average optimality gap) and the average solution time (columns Average resolution time). If convergence is not reached after one hour, we report 100% and 3600 seconds, respectively, for the two indicators.

We observe from Table 10 that among the monolithic methods, `mcr8` outperforms the others in terms of both average solution time and average optimality gap. This highlights the efficiency of the reformulated constraints in `mcr8`, which help to reduce the computational burden on the solver significantly. Additionally, we note a substantial reduction in the optimality gap by a factor of at least 40 when switching to the Benders decomposition method. Moreover, both decomposition methods, `ben` and `ccg`, exhibit lower average solution times compared to other methods, with a minimum reduction by a factor of 2 in the case of `ben`. When combining this observation with the data reported in Table 9, we can conclude that, for this phase, a coarse approximation of uncertainty proves more efficient than a fine approximation.

In conclusion, for the training phase of our experiments, while the `jab` model is easier to handle for small instances, it suffers from poor scalability and yields the lowest number of feasible solutions overall. Consequently, if a medium number of scenarios must be considered, the `jab` method becomes impractical. Furthermore, none of the monolithic methods managed to provide solutions for instances involving 300 scenarios and 80 locations.

Overall, these comparisons highlight the advantages of decomposition-based methods, as evidenced by the superior performance of `ben` and `ccg` regarding computational efficiency and robustness. Decomposition techniques excel by breaking down complex problems into smaller sub-problems, thereby better managing complexity and improving scalability. However, they still face challenges, particularly for large-scale instances, where further optimization efforts are necessary to improve performance and prevent memory issues. This analysis underscores the potential of decomposition approaches for large-scale optimization while identifying areas for future enhancement.

## 7.2 Results for the testing phase

This section reports the results of the computational experiments for the testing phase to assess the quality of the first-stage solutions produced by each method in the training phase. Since the size of the instances significantly impacts the computational results, we split the analysis into three parts, related
<table><thead><tr><th rowspan="2">|J|</th><th rowspan="2">|Ω|</th><th colspan="4">Average optimality gap (%)</th><th colspan="4">Average resolution time (second)</th></tr><tr><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th></tr></thead><tbody><tr><td rowspan="4">20</td><td>20</td><td>0</td><td>0.07</td><td>0.75</td><td>0.01</td><td>148.57</td><td>1553.29</td><td>3080.28</td><td>1244.86</td></tr><tr><td>50</td><td>0.03</td><td>0.16</td><td>1.86</td><td>0</td><td>498.44</td><td>1657.38</td><td>3124.99</td><td>1455.16</td></tr><tr><td>100</td><td>57.9</td><td>0.16</td><td>1.42</td><td>0</td><td>2528.99</td><td>2127.8</td><td>3242.5</td><td>1749.76</td></tr><tr><td>200</td><td>100</td><td>22.72</td><td>30.58</td><td>0.01</td><td>3600</td><td>2954.74</td><td>3480.07</td><td>2173.1</td></tr><tr><td></td><td>300</td><td>100</td><td>60.92</td><td>71.41</td><td>0.01</td><td>3600</td><td>3278.67</td><td>3564.36</td><td>2667.79</td></tr><tr><td rowspan="4">40</td><td>20</td><td>1.63</td><td>1.33</td><td>2.45</td><td>0.02</td><td>1026.67</td><td>2408.38</td><td>3009.04</td><td>2408.3</td></tr><tr><td>50</td><td>81.44</td><td>18.09</td><td>11.65</td><td>0.03</td><td>3141.61</td><td>2920.86</td><td>3220.32</td><td>2565.72</td></tr><tr><td>100</td><td>100</td><td>76.9</td><td>75.46</td><td>14.02</td><td>3600</td><td>3362.23</td><td>3470.9</td><td>2788.29</td></tr><tr><td>200</td><td>100</td><td>97.02</td><td>98</td><td>37.08</td><td>3600</td><td>3553.39</td><td>3581.46</td><td>3091.63</td></tr><tr><td></td><td>300</td><td>100</td><td>100</td><td>99</td><td>69.03</td><td>3600</td><td>3600</td><td>3598.5</td><td>3241.97</td></tr><tr><td rowspan="4">80</td><td>20</td><td>90.5</td><td>54.4</td><td>52.75</td><td>0.04</td><td>3600</td><td>2877.67</td><td>2892.76</td><td>2041.87</td></tr><tr><td>50</td><td>100</td><td>92.03</td><td>93.0</td><td>44.04</td><td>3600</td><td>3444.93</td><td>3415.99</td><td>2552.35</td></tr><tr><td>100</td><td>100</td><td>99.0</td><td>99.0</td><td>78.02</td><td>3600</td><td>3593.26</td><td>3582.03</td><td>3109.59</td></tr><tr><td>200</td><td>100</td><td>100</td><td>100</td><td>96</td><td>3600</td><td>3600</td><td>3600</td><td>3524.74</td></tr><tr><td></td><td>300</td><td>100</td><td>100</td><td>100</td><td>100</td><td>3600</td><td>3600</td><td>3600</td><td>3600</td></tr></tbody><tfoot><tr><td colspan="2">avg.</td><td>75.43</td><td>54.85</td><td>55.82</td><td>29.22</td><td>2889.63</td><td>2968.86</td><td>3364.21</td><td>2547.68</td></tr></tfoot></table>

<table><caption>Table 10: Average optimality gap and average resolution time per number of locations and number of scenarios</caption><thead><tr><th rowspan="2">|J|</th><th rowspan="2">|Ω|</th><th colspan="2">Average optimality gap (%)</th><th colspan="2">Average resolution time (second)</th></tr><tr><th>ben</th><th>ccg</th><th>ben</th><th>ccg</th></tr></thead><tbody><tr><td rowspan="4">20</td><td>20</td><td>0.0</td><td>-</td><td>43.06</td><td>1188.84</td></tr><tr><td>50</td><td>0.0</td><td>-</td><td>104.19</td><td>1350.32</td></tr><tr><td>100</td><td>0.0</td><td>-</td><td>208.29</td><td>1497.71</td></tr><tr><td>200</td><td>0.0</td><td>-</td><td>411.39</td><td>1787.98</td></tr><tr><td></td><td>300</td><td>0.0</td><td>-</td><td>637.44</td><td>2057.49</td></tr><tr><td rowspan="4">40</td><td>20</td><td>0.01</td><td>-</td><td>692.31</td><td>2408.63</td></tr><tr><td>50</td><td>0.01</td><td>-</td><td>892.16</td><td>2579.01</td></tr><tr><td>100</td><td>0.01</td><td>-</td><td>1308.25</td><td>2468.11</td></tr><tr><td>200</td><td>0.02</td><td>-</td><td>1626.79</td><td>2343.73</td></tr><tr><td></td><td>300</td><td>0.01</td><td>-</td><td>1964.1</td><td>2307.13</td></tr><tr><td rowspan="4">80</td><td>20</td><td>0.0</td><td>-</td><td>825.11</td><td>2147.69</td></tr><tr><td>50</td><td>0.0</td><td>-</td><td>1604.88</td><td>2108.32</td></tr><tr><td>100</td><td>0.0</td><td>-</td><td>2136.94</td><td>1905.98</td></tr><tr><td>200</td><td>6.0</td><td>-</td><td>2296.64</td><td>1976.45</td></tr><tr><td></td><td>300</td><td>4.02</td><td>-</td><td>2606.48</td><td>2545.23</td></tr></tbody><tfoot><tr><td colspan="2">avg.</td><td>0.67</td><td>-</td><td>1157.2</td><td>2044.84</td></tr></tfoot></table>
to instances with 20, 40, and 80 locations, respectively. For each instance size and each method, we first determine the best number of scenarios to use. The results of these preliminary computational experiments are presented in Appendix B. Based on the obtained results, we concluded that increasing the number of scenarios generally improves the solution quality, provided the solver can converge to an optimal solution. However, for a large number of scenarios, it becomes challenging to achieve convergence. Therefore, we prioritized convergence over a larger number of scenarios. The best configuration for each model was used for comparison, utilizing the mean-CVAR objective.

To compare the different models, we use the notion of performance profile introduced by Dolan and Moré (2002). In their work, they proposed a framework to evaluate and compare the performance of the set of solvers $S$ on a test set of problems $P$. They denoted by $w_{p,u}$ the computing time required to solve problem $p \in P$ by solver $u \in S$ and defined a baseline for comparisons name performance ratio $r_{p,u}$ computed as follows:

$$r_{p,u} = \frac{w_{p,u}}{\min\{w_{p,u'} : u' \in S\}}$$

This quantity compares the performance of a solver $u$ on a problem $p$, with the best performance of any solver on this problem. For the solver which does not solve problem $p \in P$, the ratio equals to $\hat{r}$, with $\hat{r} \ge r_{p,u}, \forall(p, u) \in P \times S$. They proved that the choice of $\hat{r}$ does not affect the performance evaluation. Therefore, the performance profile $F_u$ of a solver $u \in S$ is the distribution function of a performance metric of this solver. It is a non-decreasing, piece-wise constant function, continuous from the right at each break-point defined by the following relation:

$$F_u: \mathbb{R} \to [0, 1] \\ \tau \mapsto \frac{1}{|P|} |\{p \in P : r_{p,u} \le \tau\}|$$

$F_u(\tau)$ is the percentage of problems on which the solver $u \in S$ performs at most $\tau$ times worst than the best solver.

The performance profiles of jab, wc, and mc8 models are plotted in Figures 1a, 2a and as well as in figures given in Appendix B for the preliminary results. In these figures, the x-axis represents the performance of the method compared to the best one. Value 1 means that the best solution was found while value 5 is a special value indicating that the ratio was at least 5 or no solution was found. The y-axis represents the cumulative frequency of the corresponding ratio.

Figure 1a shows the comparative performance profile of jab, wc, mc8, and mcr8 models trained with respectively 50, 20, 100, and 200 scenarios. Figure 1b provides the distribution of the confidence intervals for these models.

<figure><img src="image_2.png" /><figcaption>Figure 1: Results for instances with 20 locations</figcaption></figure>
For these instances, mcr8 model has the best performance (abscissa equal to 1 for 99% of instances), followed by mc8 and wc models (respectively 98% and 81%), and these three approaches outperform jab model (6%). For failures (abscissa equal to 5), we can see 1% of instances for mcr8, mc8, and wc, and 32% for jab. The confidence intervals shown in Figure 1b confirm the clear dominance of wc, mc8, and mcr8 over jab while the results obtained for wc, mc8, and mcr8 are relatively close.

For 40 locations, the best results for all methods were obtained with 20 scenarios only except for mcr8 which performs for 50 scenarios. Figure 2 shows the comparative performance profile of jab, wc, and mc8 models trained with 20 scenarios. Figure 2a provides the distribution of the confidence intervals for these models. For 40 locations, we can observe similar results as for 20 locations: mcr8 provides the best results in 95% of the instances, compared to 92% for mc8, 77% for wc and 7% for jab. In terms of failures (abscissa equal to 5), wc and mc8 account for 3% of instances, while jab accounts for 42%, and mcr8 accounts for only 1%.

<figure><img src="image_3.png" /><figcaption>Figure 2: Results for instances with 40 locations</figcaption></figure>

For 80 locations, the best results for all methods were obtained with 20 scenarios only, they are reported in Figure 3. For these instances, the mcr8 model still exhibits the best performance (with an abscissa equal to 1 for 97% of instances), showing a considerable gap compared to the wc and mc8 models (respectively 43% and 49%). All three models significantly outperform the jab model, which achieves only 1%. Regarding failures, the mcr8 model registers failures for 3% of instances, the wc and mc8 records respectively 45% and 48%, and jab accounts for 99%.

<figure><img src="image_4.png" /><figcaption>Figure 3: Results for instances with 80 locations</figcaption></figure>
The obtained results demonstrate that employing risk-averse objective functions leads to more robust solutions during the preparedness phase, thereby achieving better overall performance in meeting demand when mitigating earthquake damages. This ultimately enhances the effectiveness of emergency response.

From an optimization standpoint, one important conclusion is, for approximating mean-CVAR in this problem, relaxing the integrality constraints in the second stage proves to be more effective than adopting a worst-case objective approximation, both in terms of computational time (when using a suitable decomposition algorithm, see Table 10) and solution quality.

## 7.3 Managerial insights on the budget utilisation

### Performance Analysis of the Models by Budget

Figure 4 shows the average mean-CVAR of unsatisfied demand by model and budget value used in the experiments.

<figure><img src="image_5.png" /><figcaption>Figure 4: Percentage difference from the best percentage (calculated as mean-CVAR) of unmet demand relative to demand, by model and budget</figcaption></figure>

$$ \mathrm{Perc}(|\mathcal{H}|, |\mathcal{J}|, |\mathcal{T}^*|, B, (\mathbf{x}^*, \bar{\mathbf{s}})_u) = \frac{(1-\lambda) \mathbb{E}_{\hat{\omega}}[\mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}})] + \lambda \mathrm{CVaR}_{\alpha}(\mathcal{Q}(\mathbf{x}, \bar{\mathbf{s}}, \boldsymbol{\xi}^{\hat{\omega}}))}{(1-\lambda) \mathbb{E}_{\hat{\omega}}[\sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \boldsymbol{\xi}_{h,t}^{\hat{\omega}}] + \lambda \mathrm{CVaR}_{\alpha}(\sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \boldsymbol{\xi}_{h,t}^{\hat{\omega}})} \quad (61) $$

$$ \text{gap}(|\mathcal{H}|, |\mathcal{J}|, |\mathcal{T}^*|, B, u) = \text{Perc}(|\mathcal{H}|, |\mathcal{J}|, |\mathcal{T}^*|, B, (\mathbf{x}^*, \bar{\mathbf{s}})_u) - \min_{u' \in \mathcal{S}} \{\text{Perc}(|\mathcal{H}|, |\mathcal{J}|, |\mathcal{T}^*|, B, (\mathbf{x}^*, \bar{\mathbf{s}})_{u'})\} \quad (62) $$
To plot Figure 4, we calculated the rate of unsatisfied demand as a percentage of total demand for each model $u \in S$, instance and budget using formula (61) where $(x^*, \bar{s})_u$ is the solution of model $u$, found during training phase and evaluated during testing phase. We calculated the gap (62) between each model's rate and the best rate using these values. Finally, we calculated and plotted the average gap for each model calculated over all instances and budgets.

We can notice that on average, **mc8** and **wc** models better satisfy the demand than **jab** model for the same amount of the available budget. We can also observe that this difference is even more significant for lower amounts of the available budget. This observation confirms our hypothesis that the risk-averse objective functions have a positive impact on the estimated number of saved lives and should be preferred in the design of a blood supply chain for ethically acceptable optimization in the context of mitigation of the damages caused by an earthquake.

## Analysis of Budget Allocation Between Preparation and Response

As illustrated in Figure 5, the **jab** model consistently allocates a low proportion of the budget to the preparation phase, particularly at lower budget levels, with allocations as low as 2.23%. This strategy, while cost-effective in the short term, poses significant risks during the response phase. In contrast, risk-aware models such as **wc**, **mc8**, and **mcr8** prioritize preparation, allocating a substantial share of the budget (approximately 75%-88%), thereby ensuring enhanced response capabilities.

At lower budgets, risk-averse models allocate around 57%-59% of the budget to preparation. As budgets increase, these allocations remain stable, averaging approximately 78% for models **wc**, **mc8**, and **mcr8**. The **jab** model, however, employs a different strategy, initially investing minimally in preparation but sharply increasing this share as the budget grows.

Risk-averse models underscore the importance of allocating a significant portion of the budget to preparation in uncertain or high-risk environments. This approach minimizes unforeseen complications and enhances responsiveness during the critical response phase. In contrast, the **jab** model's minimal preparation investment at lower budgets may prove detrimental in scenarios where effective anticipation is crucial.

The average preparation allocation of approximately 78% (as detailed in Table 11) provided by risk-averse models can serve as a benchmark for decision-makers in crisis management. This distribution strikes a balance, ensuring robust anticipation capabilities while reserving sufficient resources for the response phase.

In uncertain or high-stakes situations, the **wc**, **mc8**, and **mcr8** models are recommended due to their strong emphasis on risk mitigation and preparation. Conversely, the **jab** model may be more suitable in low-uncertainty environments where immediate risk is less pronounced.

This analysis highlights the critical role of risk consideration in disaster supply chain management. By evaluating the level of uncertainty and available resources, decision-makers can select the most appropriate model to optimize demand satisfaction while adhering to budget constraints.

# 8 Conclusion

In this research, we proposed risk-averse models for designing a blood supply chain to prepare for and respond to a humanitarian emergency caused by a disaster, such as an earthquake. In our study, we chose a two-stage model combining the preparation stage, where tactical decisions are made, and the response stage, where operational decisions are made. The objective is to minimize the unmet blood demand and, given the uncertainty about this demand, we proposed to use mean-CVaR a combination of risk measures: expectation (risk-neutral) and CVaR (risk-averse).

We proposed three models based on a discrete approximation of uncertainty. The first optimizes the
<figure><img src="image_6.png" /><figcaption>Figure 5: Percentage of the available budget used for the preparedness, by model and budget</figcaption></figure>

SAA approximation of mean-CVAR the second is a relaxation on the first (by relaxing binary variable indexed by scenario), and the third adopts a robust perspective and optimizes unmet demand in the worst-case scenario. We then evaluated the performance of these models on an out-of-sample set of scenarios, based on data from the literature and new, larger, instances. We computed an estimate of our risk measure mean-CVAR and its confidence interval and compared the performances of our models and the one from Jabbarzadeh et al. (2014). Numerical results show that the models using discrete mean-CVAR and worst-case objectives perform significantly better (in average 85, 67% and 76, 33% respectively) than the solution robustness and model robustness approaches (in average 9%) for the objective of the minimization of the unmet blood demand under uncertainty about this demand.

Our results also show that allocating the budget between the two phases affects the efficiency of the supply chain and that the preparation phase is crucial and clearly determines the system's performance in the response phase.

Despite their strong performance, the proposed optimization models exhibit limitations when addressing large instances. To address this, we introduced two decomposition methods that are better suited for large-scale problems due to their ability to divide the problem into more manageable subproblems. Numerical results show that these methods not only reduce the computational burden but also significantly enhance the solver's ability to explore feasible solutions efficiently.

Decomposition methods reduced the percentage of unsolved instances by at least a factor of six compared to monolithic models, decreasing from 25.2% to 0% and 38.87% to 6.2%. Additionally, they
<table><caption>Table 11: Average optimality gap and average resolution time per number of locations and number of scenarios</caption><thead><tr><th>Budget</th><th>jab</th><th>wc</th><th>mc8</th><th>mcr8</th></tr></thead><tbody><tr><td>10</td><td>2.23</td><td>59.31</td><td>57.35</td><td>57.35</td></tr><tr><td>20</td><td>8.35</td><td>76.97</td><td>75.95</td><td>75.95</td></tr><tr><td>30</td><td>43.74</td><td>84.73</td><td>82.74</td><td>82.09</td></tr><tr><td>40</td><td>35.08</td><td>83.56</td><td>85.55</td><td>85.55</td></tr><tr><td>50</td><td>31.14</td><td>88.52</td><td>86.87</td><td>86.88</td></tr><tr><td>avg.</td><td>24.13</td><td>78.62</td><td>77.69</td><td>77.56</td></tr></tbody></table>

led to a notable increase in the percentage of instances solved to optimality, from 35.13% to 82.2% and from 23.26% to 43.67%. However, a non-negligible percentage of instances still fail to converge to optimality, and some remain unsolvable due to memory constraints when binary variables indexed by scenarios are retained.

From an optimization standpoint, we observe that, for approximating mean-CVAR in this context, loosening the integrality constraints in the second stage proves to be a more practical approach than relying on a worst-case objective. This strategy not only accelerates the solution process by leveraging appropriate decomposition techniques but also enhances the overall quality of the solution, making it a more efficient alternative under computational constraints.

These findings highlight the need for more advanced resolution techniques, such as branch-and-cut algorithms, reformulations, dedicated heuristics, or alternative approximations like polyhedral uncertainty sets. The development of these techniques presents a valuable research direction to further improve the management of humanitarian emergencies caused by various types of disasters.

While our models and methods provide significant improvements in the design and management of blood supply chains for humanitarian emergencies, several limitations remain, as:

* *Resource Constraints:* Our models do not account for real-world limitations on the number of mobile collection units or human resources available during the response phase. Incorporating such constraints could further enhance the realism of the models.

* *Dynamic Demand Fluctuations:* We assume static uncertainty in demand. In practice, demand may fluctuate in real-time as more information becomes available during a disaster. Extending our models to handle dynamic or real-time demand updates would be valuable.

* *Transportation Challenges:* The current models do not explicitly account for disruptions in transportation, such as road closures or delays caused by disaster conditions. Modeling such disruptions could improve the robustness of the proposed supply chain.

* *Multi-Stakeholder Coordination:* Our models assume centralized decision-making. In reality, multiple stakeholders are often involved, which can lead to coordination issues and delays. Addressing these factors could improve the applicability of our approach in real-world settings.

Future work can explore these limitations by integrating additional real-world constraints, adopting stochastic or adaptive optimization techniques, and developing collaborative decision-making models involving multiple stakeholders.

## Availability of Source Code and Data

The source code used to implement the experiments described in this article is available as open source on GitHub but its address cannot be communicated in the submitted version for the purpose of the
Acknowledgements

This research was made possible by a grant from the ANR in the DESIDE project (ANR-20-CE40-0021).
The experiments presented in this document were carried out using the computing resources of the MCIA
(Mésocentre de Calcul Intensif Aquitain). Computer time for this study was provided by the computing
facilities of the MCIA (Mésocentre de Calcul Intensif Aquitain).

References

Afsari, N., SeyedShenava, S. J., and Shayeghi, H. (2024). A MILP Model Incorporated With the Risk Management Tool for Self-Healing Oriented Service Restoration. *Journal of Operation and Automation in Power Engineering*, 12(1):1–13.

Ahmadimanesh, M., Tavakoli, A., Pooya, A., and Dehghanian, F. (2020). Designing an optimal inventory management model for the blood supply chain: Synthesis of reusable simulation and neural network. *Medicine*, 99(29):e21208.

Ahmed, S. and Shapiro, A. (2002). The sample average approximation method for stochastic programs with integer recourse. Submitted for publication. *Science*, 12.

Artzner, P., Delbaen, F., Eber, J.-M., and Heath, D. (1999). Coherent measures of risk. *Mathematical finance*, 9(3):203–228.

Asadpour, M., Olsen, T. L., and Boyer, O. (2022). An updated review on blood supply chain quantitative models: A disaster perspective. *Transportation Research Part E: Logistics and Transportation Review*, 158:102583.

Beinsure (2023). 2023 Global Natural Catastrophes: 10 Costliest Insured Losses. https://beinsure.com/global-natural-catastrophes-report/.

Ben-Ameur, W. and Neto, J. (2007). Acceleration of cutting-plane and column generation algorithms: Applications to network design. *Networks*, 49(1):3–17. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/net.20137.

Benders, J. F. (1962). Partitioning procedures for solving mixed-variables programming problems. *Numerische Mathematik*, 4(1):238–252.

Birge, J. R. and Louveaux, F. (2011). *Introduction to Stochastic Programming*. Springer Publishing Company, Incorporated, 2nd edition.

Boonmee, C., Arimura, M., and Asada, T. (2017). Facility location optimization model for emergency humanitarian logistics. *International Journal of Disaster Risk Reduction*, 24:485–498.

Bruno, G., Diglio, A., Piccolo, C., and Cannavacciuolo, L. (2019). Territorial reorganization of regional blood management systems: Evidences from an Italian case study. *Omega*, 89:54–70.

Bushaj, S., Büyüktahtın, İ. E., and Haight, R. G. (2022). Risk-averse multi-stage stochastic optimization for surveillance and operations planning of a forest insect infestation. *European Journal of Operational Research*, 299(3):1094–1110.
Chaiwuttisak, P., Smith, H., Wu, Y., Potts, C., Sakuldamrongpanich, T., and Pathomsiri, S. (2016). Location of low-cost blood collection and distribution centres in Thailand. *Operations Research for Health Care*, 9:7–15.

Chen, L., Tang, H., Wu, J., Li, C., and Wang, Y. (2022). A robust optimization framework for energy management of CCHP users with integrated demand response in electricity market. *International Journal of Electrical Power & Energy Systems*, 141:108181.

Dolan, E. D. and Moré, J. J. (2002). Benchmarking optimization software with performance profiles. *Mathematical Programming*, 91(2):201–214.

Fahimnia, B., Jabbarzadeh, A., Ghavamifar, A., and Bell, M. (2017). Supply chain design for efficient and effective blood supply in disasters. *International Journal of Production Economics*, 183:700–709.

Fazli-Khalaf, M., Khalilpourazari, S., and Mohammadi, M. (2019). Mixed robust possibilistic flexible chance constraint optimization model for emergency blood supply chain network design. *Annals of Operations Research*, 283(1-2):1079–1109.

Fischetti, M., Ljubić, I., and Sinnl, M. (2017). Redesigning Benders Decomposition for Large-Scale Facility Location. *Management Science*, 63(7):2146–2162.

Ghatreh Samani, M. R., Torabi, S. A., and Hosseini-Motlagh, S.-M. (2018). Integrated blood supply chain planning for disaster relief. *International Journal of Disaster Risk Reduction*, 27:168–188.

Hemmelmayr, V., Doerner, K. F., Hartl, R. F., and Savelsbergh, M. W. P. (2010). Vendor managed inventory for environments with stochastic product usage. *European Journal of Operational Research*, 202(3):686–695.

Inuiguchi, M. and Ramík, J. (2000). Possibilistic linear programming: A brief review of fuzzy mathematical programming and a comparison with stochastic programming in portfolio selection problem. *Fuzzy Sets and Systems*, 111(1):3–28.

Jabbarzadeh, A., Fahimnia, B., and Seuring, S. (2014). Dynamic supply chain network design for the supply of blood in disasters: A robust model with real world application. *Transportation Research Part E: Logistics and Transportation Review*, 70:225–244.

Khalilpourazari, S. and Hashemi Doulabi, H. (2023). A flexible robust model for blood supply chain network design problem. *Annals of Operations Research*, 328(1):701–726.

Kim, S., Pasupathy, R., and Henderson, S. (2015). A Guide to Sample Average Approximation. In *Handbook of simulation optimization*, volume 216, pages 207–243. Springer.

Meneses, M., Santos, D., and Barbosa-Póvoa, A. (2023). Modelling the Blood Supply Chain. *European Journal of Operational Research*, 307(2):499–518.

Mulvey, S. M. and Vanderbei, R. J. (1995). Robust optimization of large-scale systems. *Operations Research*, 43(2):264.

Najafi, A., Pourakbari-Kasmaei, M., Jasinski, M., Lehtonen, M., and Leonowicz, Z. (2022). A max-min-max robust optimization model for multi-carrier energy systems integrated with power to gas storage system. *Journal of Energy Storage*, 48:103933.

Ni, W., Shu, J., and Song, M. (2018). Location and Emergency Inventory Pre-Positioning for Disaster Response Operations: Min-Max Robust Model and a Case Study of Yushu Earthquake. *Production & Operations Management*, 27(1):160–183.
Noyan, N. (2012). Risk-averse two-stage stochastic programming with an application to disaster management. *Computers & Operations Research*, 39(3):541-559.

Noyan, N., Meraklı, M., and Küçükyavuz, S. (2022). Two-stage stochastic programming under multivariate risk constraints with an application to humanitarian relief network design. *Mathematical Programming*, 191(1):7-45.

Pirabán, A., Guerrero, W. J., and Labadie, N. (2019). Survey on blood supply chain management: Models and methods. *Computers & Operations Research*, 112:104756.

Pisciella, P., Vespucci, M., Bertocchi, M., and Zigrino, S. (2016). A time consistent risk averse three-stage stochastic mixed integer optimization model for power generation capacity expansion. *Energy Economics*, 53:203-211.

Rahmaniani, R., Crainic, T. G., Gendreau, M., and Rei, W. (2017). The Benders decomposition algorithm: A literature review. *European Journal of Operational Research*, 259(3):801-817.

Rajendran, S. and Ravi Ravindran, A. (2019). Inventory management of platelets along blood supply chain to minimize wastage and shortage. *Computers & Industrial Engineering*, 130:714-730.

Rockafellar, R. T. and Uryasev, S. (2002). Conditional value-at-risk for general loss distributions. *Journal of Banking & Finance*, 26(7):1443-1471.

Salehi, F., Mahootchi, M., and Husseini, S. M. M. (2019). Developing a robust stochastic model for designing a blood supply chain network in a crisis: A possible earthquake in Tehran. *Annals of Operations Research*, 283(1):679-703.

Samani, M. R. G. and Hosseini-Motlagh, S.-M. (2019). An enhanced procedure for managing blood supply chain under disruptions and uncertainties. *Annals of Operations Research*, 283(1):1413-1462.

Schultz, R. and Tiedemann, S. (2006). Conditional Value-at-Risk in Stochastic Programs with Mixed-Integer Recourse. *Mathematical Programming*, 105(2):365-386.

Shapiro, A. (2003). Monte Carlo Sampling Methods. In *Handbooks in Operations Research and Management Science*, volume 10 of *Stochastic Programming*, pages 353-425. Elsevier.

Shih, H. and Rajendran, S. (2020). Stochastic Inventory Model for Minimizing Blood Shortage and Outdating in a Blood Supply Chain under Supply and Demand Uncertainty. *Journal of Healthcare Engineering*, 2020:1-14.

Slyke, R. M. V. and Wets, R. (1969). L-Shaped Linear Programs with Applications to Optimal Control and Stochastic Programming. *SIAM Journal on Applied Mathematics*, 17(4):638-663.

Tabatabaie, M., Ardalan, A., Abolghasemi, H., Naieni, K. H., Pourmalek, F., Ahmadi, B., and Shokouhi, M. (2010). Estimating blood transfusion requirements in preparation for a major earthquake: The tehran, iran study. *Prehospital and Disaster Medicine*, 25(3):246-252.

Tirkolaee, E. B., Golpîra, H., Javanmardan, A., and Maihami, R. (2023). A socio-economic optimization model for blood supply chain network design during the COVID-19 pandemic: An interactive possibilistic programming approach for a real case study. *Socio-Economic Planning Sciences*, 85:101439.

UNOCHA (2018). On-Site Operations Coordination Centre (OSOCC) Guidelines. https://www.insarag.org/guidance-notes/manuals/v-osocc/. Last accessed July 2024.
Van Denakker, T. A., Al-Riyami, A. Z., Feghali, R., Gammon, R., So-Osman, C., Crowe, E. P., Goel, R., Rai, H., Tobian, A. A., and Bloch, E. M. (2023). Managing blood supplies during natural disasters, humanitarian emergencies, and pandemics: lessons learned from COVID-19. *Expert Review of Hematology*, 16(7):501-514.

Yu, C.-S. and Li, H.-L. (2000). A robust optimization model for stochastic logistic problems. *International Journal of Production Economics*, 64(1):385-397.

Yu, G., Haskell, W. B., and Liu, Y. (2017). Resilient facility location against the risk of disruptions. *Transportation Research Part B: Methodological*, 104:82-105.

Zahiri, B. and Pishvaee, M. S. (2017). Blood supply chain network design considering blood group compatibility under uncertainty. *International Journal of Production Research*, 55(7):2013-2033.

Zahiri, B., Torabi, S. A., Mohammadi, M., and Aghabegloo, M. (2018). A multi-stage stochastic programming approach for blood supply chain planning. *Computers & Industrial Engineering*, 122:1-14.

Zeng, B. and Zhao, L. (2013). Solving two-stage robust optimization problems using a column-and-constraint generation method. *Operations Research Letters*, 41(5):457-461.

Zhang, D., Zhang, Y., Li, S., and Li, S. (2023). A novel min-max robust model for post-disaster relief kit assembly and distribution. *Expert Systems with Applications*, 214:119198.

Zhao, K., Xiang, C., and Chen, Y. (2021). Research on Blood Supply Chain Optimization under Supply Unreliable and Stochastic Demand. *Journal of Physics*.

# A Improvements of the Benders algorithm

This section details the algorithmic enhancements implemented to improve the convergence of the Benders algorithm presented in Section 5.1.

First, the relaxed master program ($RMP^B$) is initialized using the procedure described in section A.1. Then, model ($MP^B$) is solved using the *branch-and-Benders-cut* approach, i.e. solving the problem using a commercial branch-and-cut solver and adding the Benders optimality cuts within a lazy constraint callback.

## A.1 Initialization of the relaxed master program

Algorithm 3 details the procedure used to initialize the relaxed master program ($RMP^B$) at the beginning of the Benders algorithm, proposed in Fischetti et al. (2017) and applied to model ($MCVAR$).

At each iteration $k$, we solve the linear relaxation of ($RMP^B$) where only a subset (empty at the first iteration $k=0$) of the optimality cuts is included. This yields an optimal fractional solution ($x^{k*}, \bar{s}^{k*}, \eta^{k*}, z^{k*}$) of this relaxation. The simple in-out procedure uses a stabilizing point ($\tilde{x}, \tilde{s}, \tilde{\eta}, \tilde{z}$) and defines the separation point ($x^{\text{sep}}, \bar{s}^{\text{sep}}, \eta^{\text{sep}}, z^{\text{sep}}$) as ($x^{\text{sep}}, \bar{s}^{\text{sep}}, \eta^{\text{sep}}, z^{\text{sep}}$) = $\epsilon(x^{k*}, \bar{s}^{k*}, \eta^{k*}, z^{k*}) + (1-\epsilon)(\tilde{x}, \tilde{s}, \tilde{\eta}, \tilde{z})$, where $\epsilon \in [0, 1]$. The stabilization point is initially chosen as an arbitrary feasible solution of ($MP^B$) and is updated at each iteration.

Cuts are generated for subproblem solutions that are violated at the separation point. After five consecutive iterations without improvement of the value of the relaxed master program, the stabilization is stopped by setting $\epsilon = 1$ (which switches to the standard cutting plane approach). The procedure terminates when the stabilization is inactive, and either no new optimality cuts have been added at the last iteration or the value of the relaxed master program has not improved for five consecutive iterations.
To reduce the size of the initial master program, every fifth iteration and after the very last iteration, the optimality cuts that have been included in (RMP<sup>B</sup>) but are not active at the current solution (<b>x</b><sup>k*</sup>, <b>s</b><sup>k*</sup>, <b>η</b><sup>k*</sup>, <b>z</b><sup>k*</sup>) are removed from the model.

**Initial stabilizing point** The initial stabilizing point ($\tilde{\mathbf{x}}, \tilde{\mathbf{s}}, \tilde{\mathbf{\eta}}, \tilde{\mathbf{z}}$) is determined as a feasible solution to the linear relaxation of ($\overline{\text{MCVAR}}$), using the following simple heuristic. First, as long as the budget permits, fixed centers $\tilde{\mathbf{x}}$ are established, and the initial blood supply in all hospitals is set to zero ($\tilde{\mathbf{s}} = \mathbf{0}$). Next, we assume there is no flow in the network, meaning that for each scenario $\omega \in \Omega$, the total unmet demand equals the total demand. Then $\tilde{\mathbf{\eta}}$ is set as $\tilde{\mathbf{\eta}} := \text{VaR}_{\alpha, \omega \in \Omega}[\sum_{t \in \mathcal{T}^*} \sum_{h \in \mathcal{H}} \xi_{h,t}^\omega]$. From this value, we compute the positive difference $v^\omega$ between the total unmet demand of scenario $\omega \in \Omega$ and the $\tilde{\mathbf{\eta}}$. Finally, the estimation of the objective value $\tilde{z}^\omega$ of the subproblem for scenario $\omega \in \Omega$ is set as $\tilde{z}^\omega := (1-\lambda) \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \xi_{h,t}^\omega + \frac{\lambda}{1-\alpha} v^\omega$.

## A.2 Repairing heuristic for bunching

Algorithm 4 takes as input an optimal solution ($\hat{v}^\omega, \hat{y}^\omega, \hat{s}^\omega, \hat{q}^\omega, \hat{\theta}^\omega$) of problem ($SP_\omega^B(\bar{\mathbf{x}}, \bar{\mathbf{s}}, \bar{\mathbf{\eta}})$) and derives a feasible (hopefully optimal) solution for the separation problem ($SP_\omega^B(\bar{\mathbf{x}}, \bar{\mathbf{s}}, \bar{\mathbf{\eta}})$) of another scenario $\omega'$ such that $\xi_{h,t}^\omega \le \xi_{h,t}^\omega$ for all $h \in \mathcal{H}, t \in \mathcal{T}^*$ (i.e. with a smaller vector of demands). Its principle is to keep the same routing decisions ($\hat{y}^\omega := \hat{y}^\omega$ and adjust the flow variables in order to satisfy the constraints (7)-(8). Since the demands in $\omega'$ are lower than in $\omega$, one can easily show that the constructed solution ($\hat{v}^\omega, \hat{y}^\omega, \hat{s}^\omega, \hat{q}^\omega, \hat{\theta}^\omega$) is feasible for the subproblem $\omega'$. By construction, the constraints (7)-(8) are valid. As $\hat{q}^\omega \le \hat{q}^\omega$, the capacity constraints (5), (6) and (9) remain satisfied. Moreover, the economic cost of the construct solution cannot exceed the budget consumed by the solution of subproblem $\omega$, which uses more flow to satisfy more demand. Therefore, the budget constraint (2) also remains valid.

# B Performance profiles for different number of scenarios in the training phase

Figures 6 - 8 show the results for the instances with 20, 40, and 80 locations respectively and different numbers of scenarios (from 20 to 300).
<figure><img src="image_7.png" /><figcaption>Figure 6: Performance profiles for the instances with 20 locations</figcaption></figure>
<caption>**Algorithm 3:** Initialization of the Benders relaxed master program using in-out stabilization</caption>

**Data:** $\Omega$, $(\tilde{\mathbf{x}}, \tilde{\mathbf{s}}, \tilde{\eta}, \tilde{\mathbf{z}})$

**Result:** $(Q_D^\omega)_{\omega \in \Omega}$

1. $Q_D^\omega \leftarrow \emptyset$ \quad $\forall \omega \in \Omega$;
2. converged $\leftarrow$ false;
3. $\epsilon \leftarrow 0.2$;
4. old_objective $\leftarrow -\infty$;
5. $k \leftarrow 0$, no_impr $\leftarrow 0$;
6. **while** $\text{converged}$ **do**
7. &nbsp;&nbsp;&nbsp;&nbsp;new_cuts $\leftarrow 0$;
8. &nbsp;&nbsp;&nbsp;&nbsp;$(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \eta^{k*}, \mathbf{z}^{k*})$ $\leftarrow$ Solve (RMP) with $Q_D^\omega$;
9. &nbsp;&nbsp;&nbsp;&nbsp;new_objective $\leftarrow$ optimal value of (RMP) with $Q_D^\omega$;
10. &nbsp;&nbsp;&nbsp;&nbsp;$(\tilde{\mathbf{x}}, \tilde{\mathbf{s}}, \tilde{\eta}, \tilde{\mathbf{z}}) \leftarrow \frac{1}{2}[(\tilde{\mathbf{x}}, \tilde{\mathbf{s}}, \tilde{\eta}, \tilde{\mathbf{z}}) + (\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \eta^{k*}, \mathbf{z}^{k*})]$;
11. &nbsp;&nbsp;&nbsp;&nbsp;$(\mathbf{x}^{\text{sep}}, \bar{\mathbf{s}}^{\text{sep}}, \eta^{\text{sep}}, \mathbf{z}^{\text{sep}}) \leftarrow \epsilon(\mathbf{x}^{k*}, \bar{\mathbf{s}}^{k*}, \eta^{k*}, \mathbf{z}^{k*}) + (1 - \epsilon)(\tilde{\mathbf{x}}, \tilde{\mathbf{s}}, \tilde{\eta}, \tilde{\mathbf{z}})$;
12. &nbsp;&nbsp;&nbsp;&nbsp;**for** $\omega \in \Omega$ **do**
13. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$(\pi^{\omega*}, \zeta^{\omega*}, \chi^{\omega*}, \Lambda^{\omega*}, \nu^{\omega*})$ $\leftarrow$ Solve SP$^B(\mathbf{x}^{\text{sep}}, \bar{\mathbf{s}}^{\text{sep}}, \eta^{\text{sep}}, \xi^\omega)$;
14. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $z_\omega^{\text{sep}} < \pi^{\omega*} \left( \sum_{j \in \mathcal{J}} f_j x_j^{\text{sep}} + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h^{\text{sep}} - B \right) + \sum_{j \in \mathcal{J}} \sum_{t \in \mathcal{T}^*} \zeta_{j,t}^{\omega*} (x_j^{\text{sep}} - 1) - \nu^{\omega*} \eta^{\text{sep}} - \sum_{j \in \mathcal{J}} \sum_{t \in \mathcal{T}^*} \beta \chi_{j,t}^{\omega*} x_j^{\text{sep}} + \sum_{h \in \mathcal{H}} \Lambda_{h,1}^{\omega*} (\xi_{h,1}^\omega - \bar{s}_h^{\text{sep}})$ **then**
15. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$Q_D^\omega \leftarrow Q_D^\omega \cup (\pi^{\omega*}, \zeta^{\omega*}, \chi^{\omega*}, \Lambda^{\omega*}, \nu^{\omega*})$;
16. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;new_cuts $\leftarrow$ new_cuts + 1;
17. &nbsp;&nbsp;&nbsp;&nbsp;**end**
18. &nbsp;&nbsp;&nbsp;&nbsp;**end**
19. &nbsp;&nbsp;&nbsp;&nbsp;**if** $\epsilon = 1$ and (new_cuts = 0 or no_impr=5) **then**
20. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;converged $\leftarrow$ true;
21. &nbsp;&nbsp;&nbsp;&nbsp;**end**
22. &nbsp;&nbsp;&nbsp;&nbsp;**if** new_objective $\le$ old_objective **then**
23. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;no_impr $\leftarrow$ no_impr + 1;
24. &nbsp;&nbsp;&nbsp;&nbsp;**end**
25. &nbsp;&nbsp;&nbsp;&nbsp;**else**
26. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;no_impr $\leftarrow 0$;
27. &nbsp;&nbsp;&nbsp;&nbsp;**end**
28. &nbsp;&nbsp;&nbsp;&nbsp;**if** $\epsilon \ne 1$ and no_impr = 5 **then**
29. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;no_impr $\leftarrow 0$, $\epsilon \leftarrow 1$;
30. &nbsp;&nbsp;&nbsp;&nbsp;**end**
31. &nbsp;&nbsp;&nbsp;&nbsp;**if** $(k \equiv 0 \pmod 5)$ or converged **then**
32. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $\omega \in \Omega$ **do**
33. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $(\pi^{\omega*}, \zeta^{\omega*}, \chi^{\omega*}, \Lambda^{\omega*}, \nu^{\omega*}) \in Q_D^\omega$ **do**
34. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $z_\omega^{k*} > \pi^{\omega*} \left( \sum_{j \in \mathcal{J}} f_j x_j^{k*} + \sum_{h \in \mathcal{H}} \bar{c}_h \bar{s}_h^{k*} - B \right) + \sum_{j \in \mathcal{J}} \sum_{t \in \mathcal{T}^*} \zeta_{j,t}^{\omega*} (x_j^{k*} - 1) - \nu^{\omega*} \eta^{k*} - \sum_{j \in \mathcal{J}} \sum_{t \in \mathcal{T}^*} \beta \chi_{j,t}^{\omega*} x_j^{k*} + \sum_{h \in \mathcal{H}} \Lambda_{h,1}^{\omega*} (\xi_{h,1}^\omega - \bar{s}_h^{k*})$ **then**
35. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;remove $(\pi^{\omega*}, \zeta^{\omega*}, \chi^{\omega*}, \Lambda^{\omega*}, \nu^{\omega*})$ from $\in Q_D^\omega$;
36. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
37. &nbsp;&nbsp;&nbsp;&nbsp;**end**
<caption>**Algorithm 4:** Greedy algorithm for bunching</caption>

**Data:** $\eta$, $(\hat{v}^{\omega}, \hat{y}^{\omega}, \hat{s}^{\omega}, \hat{q}^{\omega}, \hat{\theta}^{\omega})$, $\omega$ and $\omega'$

**Result:** $(1 - \lambda) \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \hat{\theta}_{h,t}^{\omega'} + \frac{\lambda}{1-\alpha} \hat{v}^{\omega'}$

1 & $\hat{y}^{\omega'} \leftarrow \hat{y}^{\omega}$;
2 & $\hat{s}^{\omega'} \leftarrow \hat{s}^{\omega}$;
3 & for $h \in \mathcal{H}$ do
4 & & $\Delta_s \leftarrow 0$;
5 & & $\Delta \leftarrow 0$;
6 & & for $t = T$ to 2 do
7 & & & $\Delta_s \leftarrow \hat{s}_{h,t}^{\omega} - \hat{s}_{h,t}^{\omega'}$;
8 & & & $\Delta \leftarrow \xi_{h,t}^{\omega} - \xi_{h,t}^{\omega'} + \Delta_s$;
9 & & & $\hat{\theta}_{h,t}^{\omega'} \leftarrow \max(0, \hat{\theta}_{h,t}^{\omega} - \Delta)$;
10 & & & $\Delta \leftarrow \max(0, \Delta - \hat{\theta}_{h,t}^{\omega})$;
11 & & & $\hat{s}_{h,t-1}^{\omega'} \leftarrow \max(0, \hat{s}_{h,t-1}^{\omega} - \Delta)$;
12 & & & $\Delta \leftarrow \max(0, \Delta - \hat{s}_{h,t-1}^{\omega})$;
13 & & while $\Delta > 0$ do
14 & & & Choose $(i, j) \in \mathcal{E}$;
15 & & & & $\hat{q}_{i,j,h,t}^{\omega'} \leftarrow \max(0, \hat{q}_{i,j,h,t}^{\omega} - \Delta)$;
16 & & & & $\Delta \leftarrow \max(0, \Delta - \hat{q}_{i,j,h,t}^{\omega})$;
17 & & end
18 & end
19 & $\Delta_s \leftarrow \hat{s}_{h,1}^{\omega} - \hat{s}_{h,1}^{\omega'}$;
20 & $\Delta \leftarrow \xi_{h,t}^{\omega} - \xi_{h,t}^{\omega'} + \Delta_s$;
21 & $\hat{\theta}_{h,1}^{\omega'} \leftarrow \max(0, \hat{\theta}_{h,t}^{\omega} - \Delta)$;
22 & $\Delta \leftarrow \max(0, \Delta - \hat{\theta}_{h,1}^{\omega})$;
23 & while $\Delta > 0$ do
24 & & Choose $(i, j) \in \mathcal{E}$;
25 & & & $\hat{q}_{i,j,h,1}^{\omega'} \leftarrow \max(0, \hat{q}_{i,j,h,1}^{\omega} - \Delta)$;
26 & & & $\Delta \leftarrow \max(0, \Delta - \hat{q}_{i,j,h,1}^{\omega})$;
27 & & end
28 & end
29 & $\hat{v}^{\omega'} \leftarrow \max(0, \sum_{h \in \mathcal{H}} \sum_{t \in \mathcal{T}^*} \hat{\theta}_{h,t}^{\omega'} - \eta)$;
<figure><img src="image_8.png" /><figcaption>Figure 7: Performance profiles for the instances with 40 locations</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Performance profiles for the instances with 80 locations</figcaption></figure>