# Production-sharing of critical resources with dynamic demand under pandemic situation: The COVID-19 pandemic

Behnam Vahdani, Mehrdad Mohammadi, Simon Thevenin, Patrick Meyer, Alexandre Dolgui

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

## To cite this version:

Behnam Vahdani, Mehrdad Mohammadi, Simon Thevenin, Patrick Meyer, Alexandre Dolgui. Production-sharing of critical resources with dynamic demand under pandemic situation: The COVID-19 pandemic. Omega, 2023, 120, pp.102909. 10.1016/j.omega.2023.102909. hal-04232233

HAL Id: hal-04232233

https://hal.science/hal-04232233v1

Submitted on 11 Oct 2023

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Production-sharing of critical resources with dynamic demand under pandemic situation: The COVID-19 pandemic ☆

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

Behnam Vahdani<sup>a</sup>, Mehrdad Mohammadi<sup>a,b,*</sup>, Simon Thevenin<sup>c</sup>, Patrick Meyer<sup>a</sup>,
Alexandre Dolgui<sup>c</sup>

## ARTICLE INFO

**Article history:**

Received 28 September 2022

Accepted 26 May 2023

Available online 4 June 2023

**Keywords:**

OR in health services

Production-inventory

Resource sharing

Epidemiological model

Benders decomposition

## ABSTRACT

The COVID-19 virus's high transmissibility has resulted in the virus's rapid spread throughout the world, which has brought several repercussions, ranging from a lack of sanitary and medical products to the collapse of medical systems. Hence, governments attempt to re-plan the production of medical products and reallocate limited health resources to combat the pandemic. This paper addresses a multi-period production-inventory-sharing problem (PISP) to overcome such a circumstance, considering two consumable and reusable products. We introduce a new formulation to decide on production, inventory, delivery, and sharing quantities. The sharing will depend on net supply balance, allowable demand overload, unmet demand, and the reuse cycle of reusable products. Undeniably, the dynamic demand for products during pandemic situations must be reflected effectively in addressing the multi-period PISP. A bespoke compartmental susceptible-exposed-infectious-hospitalized-recovered-susceptible (SEIHRS) epidemiological model with a control policy is proposed, which also accounts for the influence of people's behavioral response as a result of the knowledge of adequate precautions. An accelerated Benders decomposition-based algorithm with tailored valid inequalities is offered to solve the model. Finally, we consider a realistic case study – the COVID-19 pandemic in France – to examine the computational proficiency of the decomposition method. The computational results reveal that the proposed decomposition method coupled with effective valid inequalities can solve large-sized test problems in a reasonable computational time and 9.88 times faster than the commercial Gurobi solver. Moreover, the sharing mechanism reduces the total cost of the system and the unmet demand on the average up to 32.98% and 20.96%, respectively.

© 2023 The Author(s). Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/)

## 1. Introduction

The world has been experiencing severe disasters ranging from natural (e.g., hurricanes, floods, tornadoes, pandemics, etc.) to man-made disasters (e.g., terror, error, etc.). Recently, the world has been tolerating a crucial pandemic due to the emergence of the SARS-CoV-2 virus (COVID-19, hereafter) over the past three

years. This has rapidly spread worldwide after initial exposure in China [1-3]. As of the beginning of the pandemic, all governments faced high demand for medical products/resources to control the virus's spread and cure patients in critical situations. Hence, countries worldwide agreed to increase the production of medical products and equipment ranging from medical masks to ventilators, developing health protocols, raising awareness through social media, and correcting misconceptions [4-6]. They were also aware that such an action plan necessitates re-planing medical products' production system and re-allocating/sharing limited health resources among different centers/regions to combat the pandemic as supply chain viability measures [7-12]. The success of such an action plan highly depends on two main problems: a precise enough estimation of the demand for medical products; and an efficient logistics network to produce, stock, distribute, and share such products [13-16].

<footnote><sup>a</sup> IMT Atlantique, Lab-STICC, UMR CNRS 6285, Brest F-29238, France</footnote>

<footnote><sup>b</sup> Department of Industrial Engineering and Innovation Sciences, Eindhoven University of Technology, Eindhoven 5600MB, the Netherlands</footnote>

<footnote><sup>c</sup> IMT Atlantique, LS2N-CNRS, La Chantrerie, 4, rue Alfred Kastler, Nantes cedex 3, F-44307, France</footnote>

<footnote>* Area: Production Management, Scheduling and Logistics. This manuscript was processed by Associate Editor Ivanov.</footnote>

<footnote>* Corresponding author at: Department of Industrial Engineering and Innovation Sciences, Eindhoven University of Technology, Eindhoven 5600MB, the Netherlands.</footnote>

<footnote>E-mail addresses: behnam.vahdani@imt-atlantique.fr (B. Vahdani), m.mohammadi1@tue.nl (M. Mohammadi), simon.thevenin@imt-atlantique.fr (S. Thevenin), patrick.meyer@imt-atlantique.fr (P. Meyer), alexandre.dolgui@imt-atlantique.fr (A. Dolgui).</footnote>
As far as the estimation of the demand is concerned, there are various forecasting procedures, which can provide a precise prediction of the progression of the pandemic, including agent-based models [17], metrological and meta-population models [18], compartmental epidemiological models [19–24], time-series methods [25,26], machine learning [27], and deep learning approaches [28]. Except for compartmental epidemiological models, the other approaches are highly dependent on data availability, while in the early stages of a pandemic, data is either unavailable or unreliable [29]. Furthermore, long-standing data cannot be provided due to the emergence of new variants of the virus. Also, they cannot reflect various government policies, including better treatment or public awareness programs [30]. In this regard, these methods cannot consider a specific target to guarantee the achievement of sketched control policies, such as reducing the number of deaths or system costs [31]. However, to implement control policies, integrating compartmental epidemiological models with the optimal control problem (OCP) is a powerful approach, which can cope with a significant number of the listed barriers [32–37]. More importantly, another advantage of epidemiological models is the possibility of computing an endemic equilibrium point by a dynamic optimization approach [38].

Despite all the advantages outlined for epidemiological models with OCP, researchers are less likely to utilize this approach due to the complexity of its implementation, especially in modeling the COVID-19 pandemic. What is more, although limited studies have applied this methodology, their findings cannot be used to plan logistics activities directly. Indeed, the considerations linked to their interaction with optimization models for logistics planning have not been taken into account because the presented models have only considered the aspect of epidemiology. As an illustration, we need a specific approach for identifying and estimating the demand for each piece of treatment equipment when trying to manage the supply, production, distribution, and sharing of items like ventilators, medical clothes, and ordinary and ICU beds.

Albeit a precise-enough estimation of demand aids in re-planning production and distribution decisions for medical products, it is nearly impossible to meet all demands due to limitations in proper infrastructures and facilities' production capacity. As a result, a sharing system between healthcare facilities can significantly impact the effective utilization of available inventory. In this regard, there are various approaches to provide a sharing mechanism, including queuing theory [39–41], simulation [42], agent-based algorithms [43], dynamic optimization [44], game theory [45], and data envelopment analysis [46]. In contrast, optimization models have been less employed due to their complexity. Although limited studies have been conducted to provide ventilator and patient-sharing systems among healthcare facilities with the emergence of the COVID-19 pandemic, a broad range of concerns has been neglected in these studies that should be taken into account to increase their effectiveness. As an illustration, neglecting supply and production limitations can lead to serious doubts about the feasibility of the system's results. Moreover, some medical equipment, such as ventilators and ICU beds, can be reused after an almost certain period as a lead time, which is about two weeks during the COVID-19 pandemic. So, this characteristic must be considered to update the inventory level in healthcare facilities, which may consequently affect the sharing decisions. What is more, a functional sharing system should provide the possibility that healthcare facilities can receive a certain amount of demand over their capacity. Such that they can fulfill this excess demand through service sharing instead of just accepting the demand based on a fixed capacity. In this regard, considering a broad range of factors are necessitated, including net supply balance, allowable demand overload, maximum sharing capacity, unmet demand,

and adjusting inventory level relying on a reuse cycle for reusable equipment.

According to the above-mentioned description, we address the two problems (i.e., precise-enough estimation of the dynamic demand and planning an efficient logistics network to produce, distribute, and share medical products) through a comprehensive two-stage framework. In the first stage, we suggest a SEIHRS model with the OCP for controlling the dynamicity of the demands of consumable and reusable products. Through the OCP, various control policies, including raising awareness via spreading information and improving normal and critical treatments for hospitalized people, are examined. The second stage deals with modeling a new multi-period production-inventory-sharing problem (PISP), where a broad range of decisions are made, including supply, production, distribution, inventory, and sharing.

The contributions of this study are fivefold. First, we introduce a comprehensive optimization framework to combat a pandemic. Second, we develop a tailored SEIHRS model with a control policy to handle the demand dynamicity. Third, we formulate the multi-period PISP under the equipment reuse cycle. Fourth, we offer an accelerated Benders decomposition algorithm to solve the problem. Fifth, we render new valid inequalities in the context of the multi-period PISP, which are based on the real demand, the inventory levels, and the unmet demand. The proposed optimization framework is finally applied to a realistic case study of the COVID-19 pandemic in France for the first time.

The remainder of the paper is organized as follows. Section 2 offers a comprehensive investigation concerning associated studies in the literature. Section 3 presents the description of the SEIHRS epidemiological model with the OCP. Section 4 proposes the multi-period PISP formulation. Section 5 presents the valid inequalities along with an accelerated Benders decomposition algorithm. Section 6 describes the case study and provides computational results and a set of sensitivity analyses. Section 7 provides managerial insights and implications. Finally, this study is concluded, and future research directions are proposed in Section 8.

## 2. Literature review

The current study concentrates on two main phases: 1) disease progression modeling and 2) logistics planning, especially focused on sharing/re-allocating mechanisms, to offer a comprehensive framework to combat the COVID-19 pandemic. In the following, the relevant studies of the literature are reviewed and discussed to better position the current work and its contributions compared to the literature. First, Table 1 classifies the relevant papers according to different criteria. Next, the studies of each phase are separately discussed in the following subsections.

### 2.1. Phase 1—Disease progression modeling under the COVID-19 pandemic

Considering the advantages and disadvantages listed in the previous section for disease progression modeling methods, only studies that have used epidemiology models are examined in this section. Although the COVID-19 pandemic just emerged three years ago, a significant number of studies have been conducted to predict transmission dynamics and formulate disease progression via epidemiology models. In this regard, three review articles have been presented [29,30,47,48], wherein epidemiological models such as susceptible-infected-recovered (SIR) and susceptible-exposed-infectious-recovered (SEIR) have been mostly employed. In this regard, the corresponding studies can be separated into two categories: the first category covers simply epidemiological models, while the second incorporates OCP into epidemiological models.
<table><caption>**Table 1**
Classification of the relevant papers*.</caption><thead><tr><th rowspan="3">Ref.</th><th colspan="2">Phase 1</th><th colspan="3">Phase 2</th><th colspan="5">Characteristics</th><th colspan="4">Sol.</th></tr><tr><th rowspan="2">EM-wo-OCP</th><th rowspan="2">EM-w-OCP</th><th rowspan="2">Sup.</th><th rowspan="2">Prod.</th><th rowspan="2">Dist./Shr.</th><th colspan="4">Demand</th><th rowspan="2">Resource Type</th><th rowspan="2">Inv. Cap.</th><th rowspan="2">Shr. Cap.</th><th rowspan="2">Re-use Cap.</th><th rowspan="2">Appr.</th></tr><tr><th>Identical</th><th>Distinct</th><th>Periodic</th><th>Overload</th><th>Unmet</th><th>Type</th></tr></thead><tbody><tr><td>[49]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[50]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[51]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[52]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[53]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[54]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[55]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[56]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[57]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[58]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[59]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[60]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[61]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[62]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[63]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[64]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[65]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[66]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[67]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[68]</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[69]</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[70]</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[71]</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[72]</td><td></td><td></td><td></td><td></td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[73]</td><td></td><td></td><td></td><td></td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[74]</td><td></td><td></td><td></td><td></td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[75]</td><td>✓</td><td></td><td></td><td>✓</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[76]</td><td>✓</td><td></td><td></td><td>✓</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>[77]</td><td></td><td></td><td></td><td>✓</td><td>✓</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>This study</td><td></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

In the first category, Anand et al. [49] proposed a modified SIR model (MSIRM) with a case study in India, where infected individuals were tested and quarantined. Cooper et al. [50] proposed an MSIRM with a dynamic susceptible population, where the proposed model was examined on different datasets, including China, South Korea, India, Australia, Italy, and the USA. Chen et al. [51] addressed an MSIRM with a case study in China, in which recovery and transmission rates were time-dependent. Abou-Ismail [52] applied SIR and SEIR models and offered a susceptible-quarantined-quarantined-confirmed (SUQC) model. Wangping et al. [53] extended the SIR model with a case study in Italy, where the transmission rate was time-dependent, and the Markov Chain Monte Carlo method was employed to estimate the basic reproductive number. Ifguis et al. [54] applied the SIR model for a case study in Morocco. Calafiore et al. [55] developed an MSIRM based on a case study in Italy, where their contributions focused on determining the proper number of infected people and enhancing identification and prediction processes. He et al. [56] proposed a modified SEIR model (MSEIRM) with a case study in China, where a particle swarm optimization algorithm was used to adjust the SEIR's parameters. Moreover, they examined the system's behavior regarding seasonality and stochastic infections. Moreover, in order to forecast the basic reproduction number and pandemic spread rate, Liao et al. [57] developed an integrated strategy comprising a SIR model and a polynomial regression algorithm as a machine learning method. Indeed, the SIR model's parameters were estimated using the machine learning algorithm. Similarly, Alanazi et al. [58] applied machine learning methods to estimate the parameters of a SIR model. Bagal et al. [59] applied the SIR model for analyzing various lockdown scenarios in India. Chen et al. [60] proposed a SIR model under an uncertain process by means of uncer-

tain differential equations, where the proposed model was investigated in a case study of China. Kyurkchiev et al. [61] proposed new modifications for SIR and SEIR models, wherein intervention polynomial factor was applied to consider various scenarios for the contagious disease deployment. Zisad et al. [62] incorporated a neural network model into an SEIR model to enhance the prediction accuracy, where the neural network was utilized to predict the quarantined population.

Alenezi et al. [63] proposed a time-dependent SIR model for a case study in Kuwait. Hassen et al. [64] proposed an MSIRM called SIR-Poisson to forecast the number of infected individuals following a Poisson distribution. The accuracy of the proposed model was examined in Maghreb's central regions, including Algeria, Tunisia, and Morocco. Prodanov [65] introduced an iterative approximation algorithm to estimate the parameters of the SIR model. Ala'raj et al. [78] applied ARIMA models to adjust the parameters of a modified susceptible-exposed-infected-recovered-deceased (SEIRD) model under time-dependent birth and death rates. What is more, Efimov and Ushirobira [66] proposed an MSEIRM by taking into account societal feedback on disease and confinement aspects as well as time-dependent parameters. Singh and Gupta [67] developed a generalized SIR model consisting of miscellaneous time-dependent patterns of emerging, peak amplitude, and dwindling. Also, they considered a logistic growth model to formulate infection spread. López and Rodo [68] proposed a time-dependent SEIR model for a case study in Spain and Italy, where infection spread during the latent period was considered under various proportions of containment.

In the second category, few studies have been conducted due to the complexity of reflecting OCP factors. Das and Samanta [69] suggested an MSIRM based on a fractional-order framework,

<footnote>*EM-w/wo-OCP: Epidemiological model with/without OCP; Sup.: Supply; Prod.: Production; Dist.: Distribution; Shr.: Sharing; Inv.: Inventory; Cap.: Capacity; Sol. Appr.: Solution approach.</footnote>
in which OCP was applied to investigate the effectiveness of a social distance policy, where lockdown was one of the policy's parts. Saha et al. [70] suggested a modified SEIRS model that accounts for the influence of information on adequate precautions on people's behavioral responses. They used OCP to minimize disease burden by implementing two control policies: boosting disease awareness through information and improving treatment for hospitalized patients. Saha and Samanta [71] suggested an MSEIRM in which OCP was used to assess the efficacy of social distancing and proper hygiene policy on raising awareness.

**Research gaps and questions** According to the literature, increasing awareness about disease symptoms by spreading information and enhancing treatments to hospitalized patients in different stages of the COVID-19 pandemic substantially impacts medical products demand, including hydro-alcoholic gels, masks, medical clothes, ordinary and ICU beds, and ventilators. On this subject, the demand for each of the mentioned medical products is associated with a specified part of the target population. As an illustration, all susceptible individuals must use masks and hydro-alcoholic gels to prevent the spread of the virus regarding public health protocols. Also, hospitalized individuals in general wards need ordinary beds, and hospitalized individuals in ICUs require ICU beds and ventilators. Moreover, all hospitalized individuals in general wards and ICUs need medical clothes. Therefore, reflecting these concerns to formulate a realistic disease progression model is necessitated for identifying and estimating the demand for each piece of treatment equipment. Needless to say, sketching a meaningful interaction between the estimation phase of demand for mentioned products and the planning phase of logistics activities, including production, distribution, and sharing, is undeniable. Regarding the explained research gaps, the main research questions are (1) How can we estimate demands for different medical products according to changes in disease progression? (2) How can we use estimated demand to plan logistics activities? (3) Does considering dynamic demand lead to more accurate planning for logistics activities?

## 2.2. Phase 2—Logistics planning: sharing/re-allocating mechanisms

Resource re-allocation or sharing is a reputable concept that has applications in diverse domains, including judicial services [79], closed-loop supply chain [80], bike-sharing [81], and humanitarian logistics operations [82]. However, there are a limited number of studies in the context of mathematical modeling of sharing mechanisms in healthcare systems. It is interesting to know that these limited studies have also been formed due to the emergence of the COVID-19 pandemic.

As the first papers dealing with the COVID-19 pandemic issues, Mehrotra et al. [72] proposed a ventilator sharing framework to combat the COVID-19 pandemic under demand uncertainty, where a central agency allocated ventilators to various regions under safety stock and inventory levels. In the sharing process, surplus ventilators were returned to the central agency from the regions and redistributed by the central agency with and without lead time according to the demands of other regions in which demand uncertainty was handled by stochastic programming. Parker et al. [73] proposed two sets of mathematical models to redistribute the COVID-19 patients as demands, and beds and nurses as resources among healthcare facilities regarding inventory level. Also, a robust optimization approach handled the demand uncertainty, and the objective functions minimized the total surge capacity and nurse overflow. Blanco et al. [74] proposed a single product-sharing system among healthcare facilities in the COVID-19 pandemic regarding inventory and capacity levels, where a stochastic programming approach handled the demand uncertainty. Also, they examined four objective functions that stand for different aspects of the unmet demand. Rozhkov et al. [76] proposed a simulation

model to analyze a production-inventory problem in a multi-layer supply chain system with perishable products during the COVID-19 pandemic, where two supply chain reactive strategies were considered, including order cancellation and excess inventory. On this subject, they modeled the pandemic dynamics and its disruptions by an agent-based forecast model as a separate system and analyzed its outcomes on the demand, supply, and capacities of the considered supply chain system. Li et al. [75] proposed an epidemiological model to formulate a coupled supply chain and disease network during the COVID-19 pandemic, where the proposed SEIHRS model worked alongside a production-inventory dynamic programming model. They considered four significant feedback functions to provide realistic circumstances of the COVID-19 pandemic and its disruptions, including disease-dependent capacity, disease-dependent demand, shortage-dependent transmission, and shortage-dependent fatality.

**Research gaps and questions** According to the literature, the supply and production of medical equipment have not been studied. Also, these plans are generally considered for one product type, while different products require different considerations. In this regard, some medical equipment, such as ventilators and ICU beds, can be reused after a certain period. So, this specification needs to be reflected to compute the inventory level of healthcare facilities and sharing decisions. Additionally, in real-world applications, sometimes healthcare facilities have a certain level of participation in the sharing process despite having excess capacity. Therefore, it seems necessary to consider a limited capacity for sharing process. What is more, a practical sharing system ought to allow healthcare facilities to accommodate some demand that exceeds their capacity so that they might share services to meet this surplus demand rather than simply accepting it based on a predetermined capacity. Regarding the explained research gaps, the main research questions are (1) How can we reduce health system costs and unmet demand by planning production, distribution, and sharing? (2) How does reusable products' reuse cycle affect production to sharing decisions? (3) How do consumable and reusable products plan from production to distribution in a multi-period healthcare system?

Finally, with the intention of demonstrating research gaps and our contributions better, Table 1 is provided.

## 3. Epidemiological model

This section develops a SEIHRS epidemiological model with OCP (see Fig. 1) to estimate the demands of consumable and reusable products. As a new extension, hospitalized patients are separated into two groups. The first group represents the patients in general wards who need regular hospital beds and typical treatments. The second group comprises patients under severe conditions who require ICU beds, ventilators, and critical treatments. In addition, as the physical condition of the patients in the first group deteriorates, they are transferred to the second group. More importantly, in retrospect, people show different behavior towards minor and intense symptoms of corona disease. For example, most people do not pay attention to minor symptoms and postpone visiting physicians, while this is usually the opposite in relation to intense symptoms. In addition, as people become more aware of the consequences of delaying the start of needed treatment, this attitude usually changes. Furthermore, during the COVID-19 pandemic, due to the insufficient knowledge of scientists about the origin of this disease, many of the treatments used have been based on their experience with similar diseases in the past. Therefore, with the passage of time, conducting numerous clinical trials, increasing awareness of the complications of the disease, and how patients respond to the prescribed drugs, the treatments required in minor and in-
<figure><img src="image_4.png" /><figcaption>Fig. 1. The proposed SEIHRS model.</figcaption></figure>

tense conditions have been adjusted. Consequently, various control policies are considered to reflect such circumstances.

## 3.1. Main assumptions and notations

The main assumptions to develop the SEIHRS epidemiological model are as follows: 1) As symptomatic infected population density rises, so does the accumulation of information, but over time it reaches saturation; 2) The recovery condition is not permanent and a number of recovered individuals return to the susceptible ones; 3) This model does not take into account the diseases' natural histories, such as the incubation and latency periods; 4) Every member of the population engages in social interaction; 5) Every member of the population has the same health traits, such as immune response, immunity, etc.; and 6) All afflicted individuals are contagious and disseminate the illness among the susceptible individuals.

The list of the notations of the proposed mathematical model is provided as Table 2.

The susceptible people become exposed once they encounter infected individuals who are asymptomatic, symptomatic, or hospitalized under normal or critical treatments via the term $(\beta_1A + \beta_2I + \beta_3HN + \beta_4HC)S$.

The exposed individuals are then divided into either asymptomatic or symptomatic segments regarding probabilistic $\rho$ and $(1 - \rho)$, respectively if any disease symptoms have been recognized in the infected individuals [70]. Similarly, symptomatic individuals are separated into either normal or critical hospitalized people with probabilities $(1 - c)$ and $c$, respectively, if they need supervised care in the hospital. Furthermore, the recovery condition is not permanent and a number of recovered individuals return to the susceptible segment, further with a constant rate $\xi$.

Although raising awareness persuades susceptible individuals to change their behavior and protect themselves from contagion, every person would not become cautious adequately forever due to various reasons ranging from financial matters to personality characteristics. Hence, a proportion of susceptible individuals adopt precautionary measures and alter their conventional habits, where these alterations depend not only on the number of susceptible individuals but also on the awareness density. Furthermore, the rate of information dissemination is solely determined by the number of symptomatic people. In this regard, let $u_1k$ be the rate at which the susceptible individuals' behavior changes to prevent disease spread as a result of appropriate disease prevention proceedings (e.g., social distancing, cleanliness, and self-isolation). In addition, $\frac{\partial I}{\partial t}$ is a saturation function that represents how information grows among symptomatic infected individuals [70,83].

## 3.2. The proposed SEIHRS model

The developed SEIHRS model for COVID-19 is a system dynamics model, where susceptible individuals acquire essential information concerning the disease and recommendations to avoid it. The infected individuals diverge into asymptomatic and symptomatic segments; however, due to the latency in the appearance of the symptoms, the possibility of converting asymptomatic individuals to symptomatic ones is considered. Additionally, symptomatic individuals are divided into minor and intense symptoms, which require normal and critical treatments. Furthermore, as the physical state of the patients in the first group deteriorates, they are transferred to the second group. Also, a number of hospitalized people can recover with necessary care, although they can become susceptible again due to the nature of the virus. Hence, the proposed system seeks to examine how to control the spread of the disease by investigating factors such as raising awareness about precautionary measures and offering better treatments. In this regard, the proposed SEIHRS model of Fig. 1 is mathematically formulated as the system of Eq. (1):

$$ \frac{dS}{dt} = \Lambda - (\beta_1A + \beta_2I + \beta_3HN + \beta_4HC)S - dS + \xi R - u_1kSZ \quad S(0) > 0, $$

$$ \frac{dE}{dt} = (\beta_1A + \beta_2I + \beta_3HN + \beta_4HC)S - (k+d)E \quad E(0) \geq 0, $$

$$ \frac{dA}{dt} = \kappa \rho E - \eta A - (d + \delta_1)A \quad A(0) \geq 0, $$

$$ \frac{dI}{dt} = \kappa(1 - \rho)E + \eta A - \omega I - (d + \delta_2)I \quad I(0) \geq 0, $$

$$ \frac{dHN}{dt} = \omega(1 - c)I - (d + \delta_3)HN - \gamma_1HN - \theta HN \quad HN(0) \geq 0, \quad (1) $$

$$ \frac{dHC}{dt} = \omega c I - (d + \delta_4)HC - \gamma_2HC + \theta HN \quad HC(0) \geq 0, $$

$$ \frac{dR}{dt} = \gamma_1 HN + \gamma_2 HC - \xi R - dR + u_1 kSZ \quad R(0) \geq 0, $$

$$ \frac{dZ}{dt} = \frac{pI}{1 + qI} - a_0 Z \quad Z(0) \geq 0 $$

This system represents the equilibrium of input and output flows at each state, which depends on time and simulates its fluctuation and dynamicity. In this regard, the first equation of the system (1) displays the number of susceptible people at time t. It calculates the difference between the total population, and the number of recovered individuals who return to susceptible ones, as the
<table><caption>**Table 2**
Notations of the proposed SEIHRS model.</caption><thead><tr><th>Notation</th><th>Description</th></tr></thead><tbody><tr><td><strong>Parameters</strong></td><td></td></tr><tr><td>Δ</td><td>Total population at the beginning of the pandemic</td></tr><tr><td>β<sub>1</sub>/β<sub>2</sub></td><td>Disease contagion rates per contact by asymptomatically/symptomatically infected individuals</td></tr><tr><td>β<sub>3</sub>/β<sub>4</sub></td><td>Disease contagion rates per contact by individuals hospitalized under normal/critical treatments</td></tr><tr><td>d</td><td>Natural death rate of individuals</td></tr><tr><td>δ<sub>1</sub>/δ<sub>2</sub></td><td>Particular death rate among asymptomatically/symptomatically infected individuals</td></tr><tr><td>δ<sub>3</sub>/δ<sub>4</sub></td><td>Particular death rate among individuals hospitalized under normal/critical treatments</td></tr><tr><td>ρ</td><td>Probability of being asymptomatically infected among exposed individuals</td></tr><tr><td>c</td><td>Probability of being hospitalized under critical treatments among symptomatically infected individuals</td></tr><tr><td>κ, η, ω, γ<sub>1</sub>, γ<sub>2</sub></td><td>Conversion rates from exposure to infected, asymptomatic to symptomatic, symptomatic to hospitalized, and normal and critical hospitalized to recovered individuals</td></tr><tr><td>ξ</td><td>Return rate of recovered individuals to be susceptible again</td></tr><tr><td>a<sub>0</sub></td><td>Depreciation rate of information in the course of time due to the physical fading of memory regarding sequels</td></tr><tr><td>k</td><td>Information interaction rate</td></tr><tr><td>p, q</td><td>Information growth rate and the unresponsiveness rate to information</td></tr><tr><td><strong>Variables</strong></td><td></td></tr><tr><td><em>State variables</em></td><td></td></tr><tr><td>N(t)</td><td>Total population at time t</td></tr><tr><td>S(t)</td><td>Number of susceptible individuals at time t</td></tr><tr><td>E(t)</td><td>Number of exposed individuals at time t</td></tr><tr><td>A(t), I(t)</td><td>Number of asymptomatically/symptomatically infected individuals at time t</td></tr><tr><td>HN(t), HC(t)</td><td>Number of hospitalized individuals under normal/critical treatments at time t</td></tr><tr><td>R(t)</td><td>Number of recovered individuals at time t</td></tr><tr><td>Z(t)</td><td>Awareness density because of the spread of information provided by social media platforms and educational campaigns concerning awareness programs to illustrate precaution measures (Z(t) = 0 if lack of symptomatic infected individuals)</td></tr><tr><td><em>Control variables</em></td><td></td></tr><tr><td>u<sub>1</sub></td><td>Response intensity to the information</td></tr><tr><td>u<sub>2</sub></td><td>Awareness intensity of normal symptoms</td></tr><tr><td>u<sub>3</sub></td><td>Awareness intensity of intense symptoms</td></tr><tr><td>u<sub>4</sub></td><td>Response intensity to normal treatments</td></tr><tr><td>u<sub>5</sub></td><td>Response intensity to critical treatments</td></tr></tbody></table>

input flows; the number of susceptible individuals who become exposed, the number of susceptible individuals who die, and the number of susceptible individuals who follow health protocols, as the output flows. The second equation of the system (1) indicates the number of exposed people at time t. It calculates the difference between the number of susceptible individuals who become exposed, as the input flow; the number of exposed individuals who die, and the number of exposed individuals who become infected, as the output flows. The third equation of the system (1) indicates the number of asymptomatically infected people at time t. It computes the difference between the number of exposed individuals who become asymptomatically infected, as the input flow; the number of asymptomatically infected people who die naturally and unnaturally, and the number of asymptomatically infected people who become symptomatically infected, as the output flows. The fourth equation of the system (1) indicates the number of symptomatically infected people at time t. It calculates the difference between the number of exposed individuals who become symptomatically infected, and the number of asymptomatically infected people who become symptomatically infected, as the input flows; the number of symptomatically infected people who die naturally and unnaturally, and the number of symptomatically infected people who become hospitalized, as the output flows.

The fifth equation of the system (1) indicates the number of hospitalized individuals under normal treatments at time t. It calculates the difference between the number of symptomatically infected people who become hospitalized under normal treatments, as the input flow; the number of hospitalized individuals under normal treatments who die naturally and unnaturally, the number of hospitalized individuals under normal treatments who become recovered, and the number of hospitalized individuals under normal treatments who transfer to ICUs, as the output flows. The sixth equation of the system (1) indicates the number of hospitalized individuals under critical treatment at time t. It calculates the difference between the number of symptomatically infected people who

become hospitalized under critical treatments, and the number of hospitalized individuals under normal treatments who transfer to ICUs, as the input flows; the number of hospitalized individuals under critical treatments who die naturally and unnaturally, and the number of hospitalized individuals under critical treatments who become recovered, as the output flows. The seventh equation of the system (1) indicates the number of recovered individuals at time t. It calculates the difference between the number of hospitalized individuals under normal and critical treatments who become recovered, and the number of susceptible individuals who follow health protocols, as the input flows; the number of recovered individuals who die naturally, and the number of recovered individuals who return to susceptible, as the output flows. Finally, the eighth equation of the system (1) indicates the awareness density at time t. It calculates the difference between the information growth, as the input flow, and the information depreciation, as the output flow.

For the calculation of the basic reproduction number $R_0$ and the endemic equilibrium point, interested readers are referred to Appendix A. Furthermore, the stability analysis of the $R_0$ equilibria is provided in Appendix B.

### 3.3. The optimal control problem – OCP

As discussed earlier, the main goal of the proposed SEIHRS model is to examine how to control the propagation of the disease by implementing a number of policies as control variables. Needless to say, managers are always looking for feedback from the system after the implementation of these interventions so that they can have precise control over them. In addition, it should not be forgotten that the implementation of these interventions always imposes costs on the system. This matter becomes more important when in the conditions of a pandemic; the available funds must be spent more obsessively. Therefore, an approach should be adopted that guarantees the achievement of the goal in addition to con-
sidering the mentioned limitations. However, the proposed SEIHRS model alone cannot provide the possibility of achieving the listed considerations. In this regard, OCP can be incorporated into the proposed system of equations (1) to determine control variables, satisfy the defined constraints of the devised system, and at the same time minimize the total costs of the system as a performance criterion.

In this regard, we considered two main policies, (1) raising awareness between susceptible and symptomatic infected people concerning minor and intense symptoms via spreading information, when the susceptible may observe precaution measures, and the symptomatic infected can seek treatment without disregarding their symptoms [70], and (2) improving normal and critical treatments for hospitalized individuals in general and ICU wards due to the increase of knowledge about virus behavior. In this regard and referring to Fig. 1, symptomatic infected individuals are divided into two categories. The first category visits medical centers with the mildest symptoms and receives normal treatment. The second category requires critical treatment due to ignoring the slightest symptoms, which exacerbates their health condition. This, in turn, leads to saturating general and ICU wards in medical centers. Therefore, let $\frac{\varepsilon_1 u_2(t)I}{1+\zeta_1I}$ and $\frac{\varepsilon_2 u_3(t)I}{1+\zeta_2I}$ be saturation functions for general and ICU wards, where $\varepsilon_1$ and $\varepsilon_2$ are acceptance rates of minor and severe symptoms by infected individuals with intensities $u_2$ and $u_3$ and saturation constants $\zeta_1^{-1}$ and $\zeta_2^{-1}$. Due to the unknown of COVID-19, treatments will be developed over time as more knowledge becomes available. This has an impact on disease progression and mortality rates. However, the ratio of patients to existing equipment is quite high, which leads to saturation in the system. In light of this reality, let $\frac{\varepsilon_3 u_4(t)HN}{1+\zeta_3HN}$ and $\frac{\varepsilon_4 u_5(t)HC}{1+\zeta_4HC}$ be saturation treatment functions, where $\varepsilon_3$ and $\varepsilon_4$ are normal and critical treatment rates with intensities $u_4$ and $u_5$ and saturation constants $\zeta_3^{-1}$ and $\zeta_4^{-1}$.

In this OCP, five control variables $u_i(t)$, $i = 1, \dots, 5$ ($u_i(t) \in [0, 1]$) are defined. In control variable $u_1(t)$, values 0 and 1 signify no reaction and the complete response of susceptible individuals to information, respectively. For control variables $u_2(t)$ and $u_3(t)$, value 0 signifies when a person is unaware of his or her normal or intense symptoms and does not seek medical advice, respectively; and value 1 represents complete awareness of symptoms. Finally, control variables $u_4(t)$ and $u_5(t)$, values 0 and 1 signify no response and complete response to normal and critical treatments in general and ICU wards, respectively. Hence, in order to determine the optimal response to minor and intense symptoms and optimal normal and critical treatments, an objective function is considered, which minimizes the system's total costs including the cost of spreading awareness among susceptible and symptomatic infected individuals and the cost of normal and critical treatments in hospitals. The region of the control interventions can be defined as (2):

$$ \cap = \{u_1(t), u_2(t), u_3(t), u_4(t), u_5(t) \mid u_i(t) \in [0, 1]^5, i = 1, \dots, 5; t \in [0, T_f]\} \quad (2) $$

where $T_f$ is the final time to implement the control policies, and $u_i(t)$ for $i = 1, \dots, 5$ are bounded and measurable variables. Hence, with the intention of formulating the explained control policies, the following control problem is proposed with initial conditions $S(0) > 0$, $E(0) \ge 0$, $A(0) \ge 0$, $I(0) \ge 0$, $HN(0) \ge 0$, $HC(0) \ge 0$, $R(0) \ge 0$ and $Z(0) \ge 0$.

$$ \begin{aligned} & J[u_1(t), u_2(t), u_3(t), u_4(t), u_5(t)] \\ & = \int_0^{T_f} [w_1 I(t) + w_2 HN(t) + w_3 HC(t) + w_4 u_1^2(t) \\ & \qquad + w_5 u_2^2(t) + w_6 u_3^2(t) + w_7 u_4^2(t) + w_8 u_5^2(t)] dt \end{aligned} \quad (3) $$

s.t.:

$$ \begin{aligned} \frac{dS}{dt} &= \Lambda - (\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S - dS + \xi R - u_1(t)kSZ, \\ \frac{dE}{dt} &= (\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S - (k + d)E, \\ \frac{dA}{dt} &= \kappa \rho E - (\eta + d + \delta_1)A, \\ \frac{dI}{dt} &= \kappa (1 - \rho)E + \eta A - (\omega + d + \delta_2)I - \frac{\varepsilon_1 u_2(t)I}{1 + \zeta_1 I} - \frac{\varepsilon_2 u_3(t)I}{1 + \zeta_2 I}, \\ \frac{dHN}{dt} &= \omega(1 - c)I - (d + \delta_3 + \gamma_1 + \theta)HN + \frac{\varepsilon_1 u_2(t)I}{1 + \zeta_1 I} - \frac{\varepsilon_3 u_4(t)HN}{1 + \zeta_3 HN}, \\ \frac{dHC}{dt} &= \omega cI - (d + \delta_4 + \gamma_2)HC + \theta HN + \frac{\varepsilon_2 u_3(t)I}{1 + \zeta_2 I} - \frac{\varepsilon_4 u_5(t)HC}{1 + \zeta_4 HC}, \\ \frac{dR}{dt} &= \gamma_1 HN + \gamma_2 HC - \xi R - dR + u_1(t)kSZ + \frac{\varepsilon_3 u_4(t)HN}{1 + \zeta_3 HN} + \frac{\varepsilon_4 u_5(t)HC}{1 + \zeta_4 HC}, \\ \frac{dZ}{dt} &= \frac{pl}{1 + qI} - a_0 Z \end{aligned} \quad (4) $$

The first term in the objective function (3) calculates the cost when symptomatic infected individuals consult a medical professional concerning normal and intense symptoms and seek necessary medical care without disregarding symptoms. The second and the third terms calculate the costs of hospitalized people in general and ICU wards, respectively. The fourth term calculates the cost of spreading awareness among susceptible individuals. The fifth and sixth terms compute the productivity loss due to mild or intense illness, respectively. Finally, the seventh and eighth terms calculate the opportunity loss, such as productivity loss, resulting from hospitalization to general and ICU wards, respectively. Remarkably, the square of the control variable reveals the intensity of side effects of awareness programs and treatment policies [31,32,38,84,85]. The positive parameters $w_i$, $i = 1, \dots, 8$ are importance weights and balance the units of the integrand.

**Theorem 3.** The control system (3) and (4) has optimal control variables such that $J(u_1^*, u_2^*, u_3^*, u_4^*, u_5^*) = \min J(u_1, u_2, u_3, u_4, u_5)$.

**Proof.** A mathematical proof has been provided in Appendix C. □

Finally, Appendix D provides the way to obtain optimal control variables. Noteworthy, as far as the pandemic is not over, $S(t)$, $A(t)$, $I(t)$, $HN(t)$, and $HC(t)$ help to estimate the demand of medical products corresponding to each segment of the population. For instance, considering the COVID-19 pandemic, the demands for hydro-alcoholic gels and masks over time can be estimated by $S(t)$; $HN(t) + HC(t)$ can represent the demand for medical clothes over time; the demand for ordinary beds over time can be estimated by $HN(t)$; and finally, the demands for ICU beds and ventilators over time can be estimated by $HC(t)$. More importantly, at the endemic equilibrium point, the demand for the above-mentioned products and equipment can be estimated using obtained optimal state variables ($S^*, HN^* + HC^*, HN^*$, and $HC^*$) resulting from the optimal system.

## 4. The proposed PISP model

This section formulates a logistics network that includes a set of suppliers and manufacturers that provide medical products and equipment, as well as a set of demand points/hospitals (hereafter, we only use the term hospitals as the demand and cure points of the population). The manufacturers and hospitals are clustered in a given number of regions, in such a way that the hospitals of a given region receive medical products and equipment from the manufacturers located in the same region. A set of medical
products and equipment are considered, which are divided into consumable and reusable products. Manufacturers produce medical products and equipment using raw materials provided by suppliers and distribute them only among the hospitals of their corresponding region. Due to demand fluctuations and capacity limitations in manufacturers and suppliers, hospitals might encounter shortages. Hence, an applicable sharing system among hospitals in each region is developed, wherein hospitals with surplus products can transfer the extra portion to the hospitals with extra demands. In addition, hospitals can receive a higher percentage of their patients' actual demands and respond to them through the sharing system. More importantly, a lead time is defined for reusable products according to their occupancy times in ICU wards. This is indeed in line with real-world practice, where these products become free and reusable over a certain period of time, and they will be reconsidered as available inventory.

Before entering the details of the proposed mathematical model, assumptions and necessary notations are described.

## 4.1. Assumptions and notations

The main assumptions of the defined problem are as follows: 1) The capacity of suppliers, manufacturers, and hospitals is limited; 2) The sharing capacity of any hospital, and the total sharing capacity over all regions are limited; 3) ICU beds and ventilators can be reused after a certain period as the occupancy time; 4) The shortage of medical products in hospitals is allowable, which is reflected by unmet demand; 5) Hospitals with surplus inventory can transfer the extra portion to hospitals with extra demands; 6) Hospital's unmet demand is met using the net supply balance and sharing from other hospitals; 7) Hospitals with extra inventory cannot receive products from the sharing process; 8) The sharing approach enables hospitals to get a greater proportion of their patients' actual demands and address them; and 9) The demand for medical products in hospitals is uncertain and has a dynamic pattern, which is estimated by the developed SEIHRS model with OCP.

The list of the notations of the proposed mathematical model is provided as Table 3.

## 4.2. The PISP formulation

This section proposes the multi-period PISP formulation as a mixed-integer linear programming model using the notations described in Table 3, where the demand for consumable and reusable products are estimated from the proposed SEIHRS model (4).

### 4.2.1. Objective function

The objective function of the multi-period PISP formulation is presented as (5) which minimizes the total costs of the system.

$$ \min Z = \sum_{j \in J} \sum_{t \in T} OC_j w_{jt} + \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} SC_{pm} y_{pm} + \sum_{p \in P} \sum_{j \in J} \sum_{m \in K} \sum_{t \in T} TC_{pjm} u_{pjmt} \\ + \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} PC_{pm} x_{pmt} + \sum_{p \in P} \sum_{m \in K} \sum_{m' \in H} \sum_{t \in T} TC_{pmm'} q_{mm'}^{pt} \\ + \sum_{p \in P} \sum_{m \in H} \sum_{m' \in H} \sum_{t \in T} TC_{pmm'} s_{mm'}^{pt} + \sum_{p \in P} \sum_{m \in K \cup H} \sum_{t \in T} I_{pm} I_{pmt} + CU \times S_{max} \quad (5) $$

The objective function (5) includes eight different terms including: 1) fixed ordering costs for suppliers, 2) fixed production setup costs at manufacturers, 3) variable transportation costs between suppliers and manufacturers, 4) variable production costs at manufacturers, 5) variable transportation costs between manufacturers and hospitals, 6) variable sharing costs among hospitals, 7) variable inventory holding costs at manufacturers and hospitals, and 8) unmet demand costs, respectively.

### 4.2.2. Constraints

In the following, the constraints of the proposed multi-period PISP model are explained. The first set of constraints, i.e. (6)–(12), formulate the amount of unmet demand as well as the amount of raw materials and products among suppliers, manufacturers, and hospitals.

$$ S_{max} \geq S_{pmt} \quad \forall p \in P, m \in H, t \in T \quad (6) $$

$$ \sum_{m \in K} u_{pjmt} \leq U_{pj} w_{jt} \quad \forall j \in J, p \in P, t \in T \quad (7) $$

$$ \sum_{j \in J} u_{pjmt} = \gamma_p x_{pmt} \quad \forall p \in P, m \in K, t \in T \quad (8) $$

$$ x_{pmt} \leq Q_{pm} y_{pmt} \quad \forall p \in P, m \in K, t \in T \quad (9) $$

$$ \sum_{m' \in H_i(m)} q_{mm'}^{pt} \leq x_{pm, t-1} + I_{pm, t-1} \quad \forall i \in I, p \in P, m \in K_i(m), t \in T \quad (10) $$

$$ g_{pmt} + \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{mm'}^{pt})^{+} \leq \sum_{m' \in K_i(m)} q_{m'm}^{pt} \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (11) $$

$$ \sum_{m \in K_i(m)} q_{mm'}^{pt} \leq C_{pm'}^{max} \quad \forall i \in I, p \in P, m' \in H_i(m), t \in T \quad (12) $$

In this regard, constraint (6) calculates the maximum unmet demand among all hospitals. Constraint (7) represents supplier capacities and ensures that only suppliers that have already set up can supply raw materials. Constraint (8) represents the balance of raw materials quantity transferred from suppliers to manufacturers and the production quantity in manufacturers. Next, constraint (9) guarantees the production capacity of manufacturers for each product type at each period. Furthermore, constraint (10) ensures that the production quantity transferred from each manufacturer to hospitals at the beginning of each period can not exceed the production quantity of the previous period and the inventory at the beginning of the preceding period. In addition, constraint (11) ensures that supplied and shared products via each hospital during each period cannot exceed the amount of product delivered to that hospital. Finally, constraint (12) guarantees that the delivery quantity of any type of product from all manufacturers to each hospital in each period cannot exceed the hospital's maximum supply capacity.

The second set of constraints, i.e. (13)–(17), monitor the inventory level of different products at manufacturers and hospitals.

$$ I_{pmt} = I_{pm, t-1} + x_{pm, t-1} - \sum_{m \in H_i(m)} q_{mm'}^{pt} \quad \forall i \in I, p \in P, m \in K_i(m), t \in T \quad (13) $$

$$ I_{pmt} \leq I_{pm}^{max} \quad \forall p \in P, m \in M, t \in T \quad (14) $$

$$ I_{pmt} = I_{pm, t-1} + \sum_{m' \in K_i(m)} q_{m'm}^{pt} - g_{pm, t-1} + \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{mm'}^{p,t-1})^{-} \\ - \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{mm'}^{p,t-1})^{+} \quad \forall i \in I, p \in P', m \in H_i(m), t \in T \quad (15) $$
<table><caption>**Table 3**
Notations of the proposed multi-period PISP model.</caption><thead><tr><th>Notation</th><th>Description</th></tr></thead><tbody><tr><td><strong>Sets & Indices</strong></td><td></td></tr><tr><td>m, m' &isin; M</td><td>Set of nodes, including manufacturers and hospitals in all regions</td></tr><tr><td>K & H</td><td>Set of manufacturers and hospitals in all regions</td></tr><tr><td>K<sub>i</sub>(m)</td><td>Subsets of manufacturers in region i; K = &cup;<sub>i</sub>K<sub>i</sub>(m)</td></tr><tr><td>H<sub>i</sub>(m)</td><td>Subsets of hospitals in region i; H = &cup;<sub>i</sub>H<sub>i</sub>(m)</td></tr><tr><td>j &isin; J</td><td>set of suppliers of raw materials</td></tr><tr><td>p &isin; P</td><td>Set of medical products, including two sets of consumable P' (hydro-alcoholic gel, mask, medical clothes, and ordinary beds) and reusable P'' (ICU beds and ventilator)</td></tr><tr><td>t &isin; T</td><td>Set of time periods</td></tr><tr><td><strong>Parameters</strong></td><td></td></tr><tr><td>L</td><td>Lead-time periods as the occupancy time, at which the reusable products are reconsidered as available inventory</td></tr><tr><td>&gamma;<sub>p</sub></td><td>Conversion rate of raw materials to manufacture product p</td></tr><tr><td>D<sub>pmt</sub></td><td>Demand for product p in hospital m at period t</td></tr><tr><td>U<sub>pj</sub></td><td>Capacity limitation for product p at supplier j</td></tr><tr><td>Q<sub>pmt</sub></td><td>Capacity limitation for product p at manufacturer m</td></tr><tr><td>OC<sub>j</sub></td><td>Ordering cost of product p at supplier j</td></tr><tr><td>SC<sub>pmt</sub></td><td>Production cost of product p at manufacturer m</td></tr><tr><td>PC<sub>pmt</sub></td><td>Unit production cost of product p at manufacturer m</td></tr><tr><td>TC<sub>pjm</sub></td><td>Transportation cost per unit of product p from supplier j to manufacturer m</td></tr><tr><td>CU</td><td>Unit cost of unmet demand at hospitals</td></tr><tr><td>I<sup>0</sup><sub>pmt</sub></td><td>Initial inventory of product p at manufacturer/hospital m at the beginning of the planning horizon</td></tr><tr><td>IC<sub>pmt</sub></td><td>Unit inventory holding cost of product p at manufacturer/hospital m</td></tr><tr><td>I<sup>max</sup><sub>pmt</sub></td><td>Maximum inventory capacity for product p at manufacturer/hospital m</td></tr><tr><td>G<sup>max</sup><sub>pmt</sub></td><td>Maximum sharing capacity for product p at hospital m</td></tr><tr><td>TQ<sub>p</sub></td><td>Total sharing capacity for product p over all regions</td></tr><tr><td>&mu;<sub>p</sub></td><td>Allowed demand surplus (%) for product p at hospitals that can be met through the sharing system</td></tr><tr><td>TC<sub>pmm'</sub></td><td>Unit transportation cost of product p shared among hospitals m and m' during the sharing process</td></tr><tr><td><strong>Decision variables</strong></td><td></td></tr><tr><td>w<sub>jt</sub></td><td>1 if an order is placed to supplier j at period t; 0 otherwise</td></tr><tr><td>y<sub>pmt</sub></td><td>1 if product p is set up for production in manufacturer m at period t</td></tr><tr><td>b<sub>pmt</sub></td><td>Net service balance of product p in hospital m at the beginning of period t before transferring the products in the sharing process; b<sub>pmt</sub> = b<sup>+</sup><sub>pmt</sub> + b<sup>-</sup><sub>pmt</sub>, where b<sup>+</sup><sub>pmt</sub> and b<sup>-</sup><sub>pmt</sub> are positive and negative parts of net service balance b<sub>pmt</sub> (i.e., b<sup>+</sup><sub>pmt</sub> = max{0, b<sub>pmt</sub>} and b<sup>-</sup><sub>pmt</sub> = min{0, b<sub>pmt</sub>})</td></tr><tr><td>v<sub>m</sub></td><td>1 if b<sub>pmt</sub> &gt; 0; 0 otherwise</td></tr><tr><td>x<sub>pmt</sub></td><td>Production quantity of product p in manufacturer m in period t</td></tr><tr><td>u<sub>pjm</sub></td><td>Quantity of raw material for product p transferred from supplier j to manufacturer m in period t</td></tr><tr><td>q<sup>pt</sup><sub>mm'</sub></td><td>Quantity of product p transferred from manufacturer m to hospital m' at the beginning of period t</td></tr><tr><td>l<sub>pmt</sub></td><td>Inventory of product p in node m at the end of period t</td></tr><tr><td>g<sub>pmt</sub></td><td>Available quantity of product p to be supplied by hospital m in period t to satisfy its internal demand and external sharing</td></tr><tr><td>s<sup>pt</sup><sub>mm'</sub></td><td>Quantity of product p shared from hospital m to hospital m' in period t; s<sup>pt</sup><sub>mm'</sub> = 0 and s<sup>pt</sup><sub>mm'</sub> = -s<sup>pt</sup><sub>m'm</sub>, where (s<sup>pt</sup><sub>m'm</sub>)<sup>+</sup> and (s<sup>pt</sup><sub>m'm</sub>)<sup>-</sup> are the positive and negative parts of shared products (i.e., (s<sup>pt</sup><sub>m'm</sub>)<sup>+</sup> = max{0, s<sup>pt</sup><sub>m'm</sub>} and (s<sup>pt</sup><sub>m'm</sub>)<sup>-</sup> = min{0, s<sup>pt</sup><sub>m'm</sub>} for m &ne; m')</td></tr><tr><td>S<sub>pmt</sub></td><td>Amount of unmet demand of product p in hospital m at the end of period t</td></tr><tr><td>S<sub>max</sub></td><td>Maximum unmet demand among all hospitals</td></tr></tbody></table>

$$I_{pmt} = I_{pm,t-1} + \sum_{m' \in K_i(m)} q_{m'm}^{pt} - g_{pm,t-1} + \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{m'm}^{p,t-1}) - \\ - \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{m'm}^{p,t-1})^{+} \quad \forall i \in I, p \in P'', m \in H_i(m), t \in T; 1 < t \le L \quad (16)$$

$$I_{pmt} = I_{pm,t-1} + \sum_{m' \in K_i(m)} q_{m'm}^{pt} - g_{pm,t-1} + g_{pm,t-L} + \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{m'm}^{p,t-1})^{-} \\ - \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{m'm}^{p,t-1})^{+} \quad \forall i \in I, p \in P'', m \in H_i(m), t \in T; L < t \quad (17)$$

In this regard, constraint (13) represents product inventory in each manufacturer at the beginning of each period. Constraint (14) represents the storage capacity of facilitates, including manufacturers and hospitals, for each product type and period. Constraint (15) represents product inventory in each hospital at the beginning of each period for the consumable products. Constraint (16) represents product inventory in each hospital at the beginning

of periods 2 to L for the reusable products. Constraint (17) represents product inventory in each hospital from the beginning of period L-onwards for reusable products regarding the reuse cycle.

The rest of the constraints, i.e. (18)-(28), are defined to form the sharing system.

$$D_{pmt} \le (1 + \mu_p)g_{pmt} \quad \forall p \in P, m \in H_i(m), t \in T \quad (18)$$

$$\sum_{m \in H} g_{pmt} \le TQ_p \quad \forall p \in P, t \in T \quad (19)$$

$$g_{pmt} \le I_{pmt} \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (20)$$

$$b_{pmt} = g_{pmt} - D_{pmt} \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (21)$$

$$S_{pmt} = -b_{pmt} - \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} s_{m'm}^{pt} \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (22)$$

$$b_{pmt}^{+} \le TQ_p v_m \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (23)$$

$$b_{pmt}^{-} \ge TQ_p(v_m - 1) \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (24)$$
$$ \sum_{\substack{m' \in H_i(m) \\ m' \neq m}} (s_{mm'}^{pt})^{+} \leq b_{pmt}^{+} \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \quad (25) $$

$$ b_{pmt}^{+} \geq (s_{mm'}^{pt})^{+} \quad \forall i \in I, p \in P, m, m' \in H_i(m), t \in T; m \neq m' \quad (26) $$

$$ b_{pmt}^{-} \leq (s_{mm'}^{pt})^{-} \quad \forall i \in I, p \in P, m, m' \in H_i(m), t \in T; m \neq m' \quad (27) $$

$$ (s_{mm'}^{pt})^{+} \leq G_{pm}^{max} \quad \forall i \in I, p \in P, m, m' \in H_i(m), t \in T; m \neq m' \quad (28) $$

For this aim, constraint (18) ensures that each hospital's demand for each product type in each period is met by the amount of hospital supply and the allowable demand overload ratio. Constraint (19) guarantees that all hospitals' total available supply for any product cannot exceed the total supply capacity of that hospital. Constraint (20) guarantees that the supply quantity of each product type in each period via each hospital cannot exceed the hospital inventory at the beginning of the period. Constraint (21) represents the net supply balance of each product in each hospital before sharing. Constraint (22) states that each hospital's unmet demand for any product is met using the net supply balance and product sharing from other hospitals. Constraints (23) and (24) impose that only one of $b_{pmt}^{+}$ and $b_{pmt}^{-}$ becomes nonzero. Constraint (25) guarantees that the total product transferred from each hospital to other hospitals should not exceed that hospital's maximum transferable capacity. Constraint (26) ensures that each hospital with extra supply can offer an outgoing capacity transfer. Constraint (27) ensures that a hospital with extra supply can not receive an incoming capacity transfer. Constraint (28) prevents product sharing across regions. Finally, constraints (29) determine the types of decision variables.

$$ \begin{aligned} & u_{p,jmt}, x_{pmt}, q_{mm'}^{pt}, I_{pmt}, g_{pmt}, b_{pmt}^{+}, (s_{mm'}^{pt})^{+}, S_{pmt}, S_{max} \geq 0 \\ & s_{mm'}^{pt}, b_{pmt} \in \mathbb{R} \\ & b_{pmt}^{-}, (s_{mm'}^{pt})^{-} \leq 0 \\ & v_m, w_{jt}, y_{pmt} \in \{0, 1\} \end{aligned} \quad (29) $$

## 5. Benders decomposition

This section offers an accelerated Benders decomposition (ABD) algorithm [86] to solve the proposed PISP model. Benders decomposition divides the original problem into a master problem and several sub-problems such that each of which is generally easier to solve than the original problem. The sub-problems' variables are estimated by applying linear programming duality. The remaining variables are in the master problem, as well as an artificial variable that describes a lower bound (upper bound) on the sub-problems' objective function for a minimization (maximization) problem. Afterward, a cutting plane algorithm solves the resultant model, where the values of the master problem's variables are first determined. Then, the sub-problems are solved with these variables determined at each iteration. A feasibility cut is added to the master problem if the sub-problems are infeasible and unbounded; else, an optimality cut is added. If sub-problems are feasible, an upper bound can be obtained, and if the optimal solution is obtained by solving the master problem, a lower bound can be obtained. The process is repeated until an optimal solution is detected or the optimality gap falls below a specified threshold. What is more, in order to enhance Benders decomposition's performance, a number of acceleration methods are applied.

### 5.1. Valid inequalities

This section introduces a set of new valid inequalities that can be applied to strengthen the linear programming relaxation of the PISP formulation. The first class of inequalities depends on the products' type (i.e., consumable and reusable). Three factors are considered in the valid inequalities for consumable products (i.e., Inequality (30)): inventory level, unmet demand quantity, and aggregated demand at hospitals. In addition to these factors, the valid inequalities of reusable products (i.e., Inequalities (31) and (32)) encompass the released inventory level due to the reuse cycle.

$$ \sum_{m \in H_i(m)} I_{pmt} + \sum_{m \in H_i(m)} S_{pmt} \geq \sum_{m \in H_i(m)} \sum_{h \leq t' \leq t} D_{pmt'} \left( 1 - \sum_{m \in K_i(m)} \sum_{h \leq t' \leq t} y_{pmt'} \right) \quad \forall i \in I, p \in P', t, h \in T; 1 \leq h \leq t \quad (30) $$

$$ \sum_{m \in H_i(m)} I_{pmt} + \sum_{m \in H_i(m)} S_{pmt} \geq \sum_{m \in H_i(m)} \sum_{h \leq t' \leq t} D_{pmt'} \left( 1 - \sum_{m \in K_i(m)} \sum_{h \leq t' \leq t} y_{pmt'} \right) \quad \forall i \in I, p \in P'', t, h \in T; 1 \leq t \leq L, 1 \leq h \leq t \quad (31) $$

$$ \begin{aligned} & \sum_{m \in H_i(m)} I_{pmt} + \sum_{m \in H_i(m)} S_{pmt} + \sum_{m \in H_i(m)} D_{pmt, t-L} \geq \\ & \qquad \sum_{m \in H_i(m)} \sum_{h \leq t' \leq t} D_{pmt'} \left( 1 - \sum_{m \in K_i(m)} \sum_{h \leq t' \leq t} y_{pmt'} \right) \\ & \quad \forall i \in I, p \in P'', t, h \in T; L < t, L < h \leq t \end{aligned} \quad (32) $$

Inequality (30) validates that the aggregated demand between periods $h$ to $t$ is typically fulfilled either from available inventory at hospitals through direct supply or sharing or from backlog if the required inventory is not available. Inequalities (31) and (32) represent a similar concept, except that after period $L$, the inventory level of reusable products is increased by the amount of demand at period $t-L$, which can be employed to fulfill demands through direct supply or sharing.

Moreover, since the primary master problem in this study only contains logical constraints upon $y_{pmt}$ and a limited number of optimality cuts, the master problem contributes to a supply chain with a restricted number of equipped suppliers and manufacturers in the initial iterations of the Benders algorithm. As a result, the lower bounds for early master problems are negligible. More importantly, only a tiny portion of demands can be fulfilled in sub-problems with few equipped facilities, imposing a high cost of unmet demand on the problem. In order to reduce this repercussion, the second class of inequalities is in charge of this task. In this regard, since suppliers provide raw materials to manufacturers and manufacturers make finished products, we must set up at least one supplier and manufacturer, which is guaranteed by inequalities (33) and (34). Finally, inequality (35) validates that the unmet demand of product $p$ in hospital $m$ during each period cannot exceed the product's demand in that hospital.

$$ \sum_{j \in J} \sum_{t \in T} w_{jt} \geq 1 \quad (33) $$

$$ \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} y_{pmt} \geq 1 \quad (34) $$

$$ S_{pmt} \leq D_{pmt} \quad \forall p \in P, m \in H, t \in T \quad (35) $$
## 5.2. Knapsack inequalities

The knapsack inequality is added to the master problem to accelerate the branch-and-bound process in the solver. Indeed, the convergence speed of Benders decomposition can be improved by adding the following cuts, which assist progressive solvers like Gurobi in deriving a range of valid inequalities [87]. Consequently, the following cuts (36) are added in iteration n + 1 of the master problem.

$$ LB^n \leq \sum_{j \in J} \sum_{t \in T} OC_j w_{jt} + \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} SC_{pm} y_{pmt} + \zeta \quad (36) $$

where $LB^n$ signifies the best specified lower bound found so far, and $\zeta$ is an additional variable representing the total cost except ordering and setup costs. Finally, the PISP Benders reformulation is provided in Appendix F.

# 6. Computational experiments

This section presents the computational results of the devised framework to deal with pandemic situations under dynamic demand. For this aim, Section 6.1 designs a set of experiments based on a real case study, the COVID-19 pandemic in France. Next, Section 6.2 provides a comparative analysis between the proposed ABD algorithm and the Gurobi solver without (GRB) and with the valid inequalities (GRB_VI). The comparison between GRB and GRB_VI helps to investigate how much the valid inequalities have been effective in reducing the computational time of the Gurobi solver. Finally, a set of comprehensive sensitivity analyses are presented in Section 6.3 to investigate the behavior of the proposed PISP model to any changes in its input parameters.

## 6.1. Experimental design and case study

This section designs a set of experiments for evaluating the devised framework's performance in dealing with the COVID-19 pandemic in France, where the demands of consumable and reusable products are determined using the proposed SEIHRS model. The experiments are derived from a real case study, the COVID-19 pandemic in France. In this case study, a set of six medical products and equipment are considered, which are divided into consumable and reusable products, including masks, hydro-alcoholic gel, medical clothes, and ordinary beds as consumable products, and ICU beds and ventilators as reusable ones. Among 13 metropolitan and 5 overseas regions of France, this case includes 12 regions, and each region has a given number of factories, hospitals, and population as depicted in Fig. 2. Furthermore, there exist 12 suppliers to supply raw materials throughout France. For each region, we nominate a set of geographically dispersed hospitals throughout the region, and each may also integrate a set of other local/smaller hospitals. In the case of integrating multiple hospitals, the capacity of the representative hospital will be the sum of the capacity of integrated hospitals. Each main hospital or a group of hospitals is responsible for fulfilling the demand of a populated area under their coverage. Accordingly, the dynamic demand of products at each hospital (i.e., $D_{pmt}$) is obtained by running the proposed system dynamics model (4) for that area, where the total population and death rates are the mean values for covered populated areas, and the population is the sum of the populations in those areas. Moreover, we assume that the disease does not spread among different regions and the proposed SEIHRS model is executed separately for each region. Detailed information on all other parameters in the SEIHRS model and PISP formulation is provided in

Figure 3 illustrates the outcomes of the SEIHRS model in determining the state variables for two distinct hospitals that belong to two different regions. It is worth mentioning that these

results only represent the first peak of the disease. As can be observed, the dynamicity of the state variables in the two hospitals are significantly different, and they do not follow a specific distribution. These issues demonstrate the necessity of employing epidemiological-based system dynamics models (i.e., the SEIHRS model) to determine the demand. As a result, the policy-makers in each region should plan separately based on their production and sharing capacities. Needless to say, this knowledge leads to better utilization of the capacity. What is more, national policy-makers may also arrange to transfer excess capacity to other regions.

Remarkably, since the basic reproduction number $R_0$ is obtained greater than 1 for the case study, the endemic equilibrium point of the pandemic cannot be obtained by solving the system of Eq. (A.6). Hence, it can be realized that the optimal values of the state variables cannot be used to estimate the demand for masks and hydro-alcoholic gels ($S^*$), ordinary beds ($HN^*$), ICU beds and ventilators ($HC^*$), and medical clothes ($HN^* + HC^*$). Therefore, the cumulative value of each of the state variables ($S(t), HN(t), HC(t), HN(t) + HC(t)$) in a week represents the demand for medical products in each period of the PISP model.

To test the performance of the proposed ABD algorithm, a set of 25 test problems are derived from the case study, wherein the numbers of regions, suppliers, manufacturers, and hospitals are fixed, but the period increases from 4 weeks to 52 weeks. Since the SEIHRS model is run for each hospital for a time horizon of 365 days and it is also impossible and illogical to plan the PISP daily, each period of the PISP model represents a week and the demand of each period is the cumulative demand of the whole week.

The PISP model and the ABD algorithm were coded in Python 3 using the Gurobi library, and all experiments were done on a server containing four Intel XEON processors with 5 GB of RAM running at 2.3 GHz. Furthermore, two stopping criteria are considered when executing the accelerated Benders decomposition algorithm and the Gurobi solver: 1) a gap of 1% and 2) a CPU time of 7200s. For the ABD algorithm, the first criterion is the gap (%) between the obtained lower-bound and upper-bound at each iteration of the algorithm; and for the Gurobi solver, it is the gap (%) between the best-found solution and the current obtained solution.

## 6.2. Numerical results

This section provides a comparative analysis of the performance of the proposed ABD algorithm with the GRB and GRB_VI. Table 4 shows the results of this comparison for 25 test problems (i.e., column "$|T|$" as the number of periods), in terms of both objective function values (i.e., columns "Obj. Values") and computational time (i.e., columns "Time (s)"). In columns "Obj. Values", the values are proportional to the objective function value of the first test problem, which has been shown as "X" in Table 4.

As can be seen in Table 4, all three methods have been able to obtain optimal solutions for test problems up to 26 time periods except the test problem with 20 periods; however, this has occurred in less computational efforts for the ABD algorithm compared to both GRB and GRB_VI. Furthermore, the GRB_VI has been also faster compared to the GRB in the mentioned test problems. In addition, the ABD algorithm has been able to obtain a gap of less than 1% for all test problems before reaching the maximum allowable CPU time of 7200s, even for larger test problems. It has also happened for the GRB_VI for almost all test problems except the last four test problems. Indeed, the proposed valid inequalities help the mathematical model to be solved faster, and they obtain better results under limited computational efforts. However, they are still ineffective when solving large-sized test problems for $|T| \ge 46$.

On the other hand, the GRB has been unable to obtain the optimal solution for test problems with 20 and $\ge 28$ periods, among which for test problems with $|T| \ge 42$ a feasible solution has not
<figure><img src="image_5.png" /><figcaption>Fig. 2. The case study information on each region.</figcaption></figure>

been even obtained within the limited computational effort of 7200s. It is also clear that there is a considerable difference between the computational time of GRB and GRB_VI compared to the ABD algorithm, which indicates the effectiveness of the proposed ABD algorithm in solving the PISP.

More importantly, to demonstrate the efficiency of the sharing mechanism, two critical evaluation criteria have been examined for all test problems, namely the total cost of the system and the maximum unmet demand. These evaluations have been reported in columns “Obj.” and “S<sub>max</sub>” under “Sharing impacts (%),” respectively, where positive and negative values, respectively, show the increase and the decrease percentages in the corresponding values. For this aim, the ABD algorithm has been used to solve all test problems under two options: 1) sharing and 2) without sharing. In the second option, the parameter G<sub>max</sub> is fixed to 0 to prohibit the hospitals from sharing products. As can be seen, the sharing mechanism reduces the total system costs and unmet demand in all test problems except for test problems with a smaller number of periods. For instance, the total cost of the system is increased by 15.71% and 2.35% for test problems |T| = 4 and |T| = 6,

respectively. The unmet demand is also increased for test problems |T| = {4, 6, 8, 10}. The reason for these exceptional increases in the total cost of the system goes back to the fact that in primer periods, the demand for medical products increases sharply up to period t = 10 (when the peak of the pandemic occurs in most of the populated areas), consequently, the maximum unmet demand increases as well. In such situations, since the planning horizon is short, the system increases the inventory levels with the hope to fulfill the demand in all possible ways (i.e., direct or sharing fulfillment), while this increase is not that significant to absorb the effect of the unmet demand cost in the objective function. Therefore, the maximum unmet demand and the total cost of the system increase. However, when a higher number of periods are considered, the system allows for a significant increase in the inventory levels even in primary periods to fulfill the demand, as well as the fact that the extra inventory is allowed to be shared in further periods. In addition, this phenomenon indirectly decreases the total cost where the system prefers a higher transfer quantity of products in the lower frequency of supplier and manufacturer setups. Indeed, the system's increase in transportation costs
<figure><img src="image_6.png" /><figcaption>Fig. 3. The outcomes of the SEIHRS model showing the dynamic demands.</figcaption></figure>

<table><caption>**Table 4**
ABD vs. GRB & GRB_VI and the impact of sharing.</caption><thead><tr><th rowspan="2">|T|</th><th colspan="3">Obj. values</th><th colspan="3">Time (s)</th><th colspan="2">Sharing impact (%)</th></tr><tr><th>ABD</th><th>GRB</th><th>GRB_VI</th><th>ABD</th><th>GRB</th><th>GRB_VI</th><th>Obj.</th><th>S<sub>max</sub></th></tr></thead><tbody><tr><td>4</td><td>1.0X</td><td>1.0X</td><td>1.0X</td><td>9</td><td>28</td><td>18</td><td>15.71</td><td>25.99</td></tr><tr><td>6</td><td>2.0X</td><td>2.0X</td><td>2.0X</td><td>15</td><td>56</td><td>50</td><td>2.35</td><td>31.81</td></tr><tr><td>8</td><td>2.8X</td><td>2.8X</td><td>2.8X</td><td>74</td><td>1586</td><td>1220</td><td>-10.83</td><td>14.03</td></tr><tr><td>10</td><td>3.8X</td><td>3.8X</td><td>3.8X</td><td>56</td><td>750</td><td>623</td><td>-21.95</td><td>9.64</td></tr><tr><td>12</td><td>4.4X</td><td>4.4X</td><td>4.4X</td><td>129</td><td>2313</td><td>2005</td><td>-29.83</td><td>-9.77</td></tr><tr><td>14</td><td>5.3X</td><td>5.3X</td><td>5.3X</td><td>178</td><td>1548</td><td>1248</td><td>-35.08</td><td>-21.06</td></tr><tr><td>16</td><td>6.3X</td><td>6.3X</td><td>6.3X</td><td>412</td><td>5217</td><td>4350</td><td>-38.32</td><td>-29.07</td></tr><tr><td>18</td><td>6.9X</td><td>6.9X</td><td>6.9X</td><td>463</td><td>3560</td><td>3019</td><td>-40.61</td><td>-32.58</td></tr><tr><td>20</td><td>8.1X</td><td>33.4X</td><td>8.1X</td><td>558</td><td>7200*</td><td>2750</td><td>-42.02</td><td>-35.16</td></tr><tr><td>22</td><td>8.7X</td><td>8.7X</td><td>8.7X</td><td>623</td><td>5661</td><td>2572</td><td>-42.79</td><td>-35.31</td></tr><tr><td>24</td><td>9.3X</td><td>9.3X</td><td>9.3X</td><td>1131</td><td>6183</td><td>3334</td><td>-43.12</td><td>-35.20</td></tr><tr><td>26</td><td>10.3X</td><td>10.3X</td><td>10.3X</td><td>1228</td><td>6607</td><td>2873</td><td>-43.19</td><td>-36.44</td></tr><tr><td>28</td><td>11.5X</td><td>42.6X</td><td>11.5X</td><td>2002</td><td>7200*</td><td>6800</td><td>-42.94</td><td>-35.95</td></tr><tr><td>30</td><td>25.4X</td><td>44.8X</td><td>25.4X</td><td>2315</td><td>7200*</td><td>665</td><td>-42.50</td><td>-35.35</td></tr><tr><td>32</td><td>13.1X</td><td>47.3X</td><td>13.1X</td><td>2505</td><td>7200*</td><td>1316</td><td>-41.74</td><td>-34.68</td></tr><tr><td>34</td><td>27.9X</td><td>49.6X</td><td>27.9X</td><td>2580</td><td>7200*</td><td>1553</td><td>-40.95</td><td>-30.01</td></tr><tr><td>36</td><td>29.5X</td><td>52.0X</td><td>29.5X</td><td>2692</td><td>7200*</td><td>1661</td><td>-40.15</td><td>-29.22</td></tr><tr><td>38</td><td>35.1X</td><td>54.3X</td><td>35.1X</td><td>2858</td><td>7200*</td><td>1662</td><td>-39.26</td><td>-28.42</td></tr><tr><td>40</td><td>17.7X</td><td>56.7X</td><td>17.7X</td><td>3034</td><td>7200*</td><td>2044</td><td>-38.33</td><td>-27.61</td></tr><tr><td>42</td><td>59.1X</td><td>**</td><td>101.6X</td><td>3669</td><td>**</td><td>7200*</td><td>-37.36</td><td>-26.81</td></tr><tr><td>44</td><td>61.6X</td><td>**</td><td>104.5X</td><td>5413</td><td>**</td><td>7200*</td><td>-36.32</td><td>-26.03</td></tr><tr><td>46</td><td>64.1X</td><td>**</td><td>**</td><td>5665</td><td>**</td><td>**</td><td>-35.31</td><td>-25.28</td></tr><tr><td>48</td><td>66.6X</td><td>**</td><td>**</td><td>6029</td><td>**</td><td>**</td><td>-34.23</td><td>-24.55</td></tr><tr><td>50</td><td>69.0X</td><td>**</td><td>**</td><td>6743</td><td>**</td><td>**</td><td>-33.32</td><td>-23.85</td></tr><tr><td>52</td><td>71.5X</td><td>**</td><td>**</td><td>7126</td><td>**</td><td>**</td><td>-32.33</td><td>-23.2</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>Average:</td><td>-32.98</td><td>-20.96</td></tr></tbody></table>

<footnote>* Time limit of 7200 s reached.</footnote>

<footnote>** No feasible solution after 7200 s.</footnote>
<table><caption>**Table 5**
Comparison between ABD & GRB & GRB_VI.</caption><thead><tr><th rowspan="2">|T|</th><th colspan="3">Obj. values (%)</th><th colspan="3">Time (ratio)</th></tr><tr><th>ABD vs. GRB</th><th>ABD vs. GRB_VI</th><th>GRB_VI vs. GRB</th><th>ABD vs. GRB</th><th>ABD vs. GRB_VI</th><th>GRB_VI vs. GRB</th></tr></thead><tbody><tr><td>4</td><td>0.00</td><td>0.00</td><td>0.00</td><td>3.16</td><td>2.02</td><td>1.56</td></tr><tr><td>6</td><td>0.00</td><td>0.00</td><td>0.00</td><td>3.73</td><td>3.33</td><td>1.12</td></tr><tr><td>8</td><td>0.00</td><td>0.00</td><td>0.00</td><td>21.43</td><td>16.49</td><td>1.30</td></tr><tr><td>10</td><td>0.00</td><td>0.00</td><td>0.00</td><td>13.40</td><td>11.13</td><td>1.20</td></tr><tr><td>12</td><td>0.00</td><td>0.00</td><td>0.00</td><td>17.93</td><td>15.54</td><td>1.15</td></tr><tr><td>14</td><td>0.00</td><td>0.00</td><td>0.00</td><td>8.70</td><td>7.01</td><td>1.24</td></tr><tr><td>16</td><td>0.00</td><td>0.00</td><td>0.00</td><td>12.66</td><td>10.56</td><td>1.20</td></tr><tr><td>18</td><td>0.00</td><td>0.00</td><td>0.00</td><td>7.69</td><td>6.52</td><td>1.18</td></tr><tr><td>20</td><td>-75.91</td><td>0.00</td><td>-75.91</td><td>*</td><td>4.93</td><td>*</td></tr><tr><td>22</td><td>0.00</td><td>0.00</td><td>0.00</td><td>9.09</td><td>4.13</td><td>2.20</td></tr><tr><td>24</td><td>0.00</td><td>0.00</td><td>0.00</td><td>5.47</td><td>2.95</td><td>1.85</td></tr><tr><td>26</td><td>0.00</td><td>0.00</td><td>0.00</td><td>5.38</td><td>2.34</td><td>2.30</td></tr><tr><td>28</td><td>-72.89</td><td>0.00</td><td>-72.89</td><td>*</td><td>3.40</td><td>*</td></tr><tr><td>30</td><td>-43.34</td><td>0.00</td><td>-43.34</td><td>*</td><td>13.29</td><td>*</td></tr><tr><td>32</td><td>-72.25</td><td>0.00</td><td>-72.25</td><td>*</td><td>22.73</td><td>*</td></tr><tr><td>34</td><td>-43.88</td><td>0.00</td><td>-43.88</td><td>*</td><td>24.79</td><td>*</td></tr><tr><td>36</td><td>-43.28</td><td>0.00</td><td>-43.28</td><td>*</td><td>25.75</td><td>*</td></tr><tr><td>38</td><td>-35.45</td><td>0.00</td><td>-35.45</td><td>*</td><td>24.69</td><td>*</td></tr><tr><td>40</td><td>-68.81</td><td>0.00</td><td>-68.81</td><td>*</td><td>28.60</td><td>*</td></tr><tr><td>42</td><td>**</td><td>-41.82</td><td>*</td><td>**</td><td>*</td><td>*</td></tr><tr><td>44</td><td>**</td><td>-41.04</td><td>*</td><td>**</td><td>*</td><td>*</td></tr><tr><td>46</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td></tr><tr><td>48</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td></tr><tr><td>50</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td></tr><tr><td>52</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td><td>**</td></tr></tbody><tfoot><tr><td colspan="3">Average (in case of optimality of GRB):</td><td>9.88</td><td>7.46</td><td>1.48</td></tr><tr><td colspan="3">Average (in case of optimality of GRB_VI):</td><td>N/A</td><td>12.12</td><td>N/A</td></tr></tfoot></table>

is more tolerable/affordable compared to the increase in setup costs.

Regardless of the system's total cost increase, the sharing mechanism has led to a significant decrease of 32.98% and 20.96% on average in the total cost and the maximum unmet demand, respectively. Furthermore, these savings in the majority of the test problems are significantly more than the average reported savings. These maximum savings belong to test problems with $|T| \approx 24$. After a detailed investigation, we figured out that the maximum savings happen in test problems wherein the system experiences the pandemic's peak demands. This, in turn, demonstrates the effectiveness and outstanding capabilities of the proposed framework to handle the sharp increase in demands whenever required.

Moreover, Table 5 provides a pairwise comparative analysis between the proposed ABD algorithm, the GRB, and the GRB_VI, in terms of both the objective function gap (%) and the computational time ratio. In columns "Obj. Values", negative values indicate that the first algorithm has obtained a lower objective function compared to the second algorithm (i.e., first algorithm vs. the second algorithm) under a limited computational time of 7200s. In column "Time (ratio)", the values represent how many times the first algorithm is faster than the second algorithm. The column "ABD vs. GRB" of "Obj. Values" reveals that, under the limited computational time of 7200s, the proposed ABD algorithm can obtain much better solutions in terms of objective function values compared to GRB with an average of -57% (i.e., the solutions of ABD algorithm are 57% better than the solutions of GRB without valid inequalities). Looking at column "ABD vs. GRB_VI" in "Obj. Values" figures out that both ABD and GRB_VI obtain the optimal solutions for test problems with $|T| = \{4-40\}$, but with a much less computational time for the ABD algorithm. For two test problems $|T| = \{42, 44\}$, the GRB_VI has not been able to obtain the optimal solution, but its best-found solutions compared to the optimal solutions of ABD are in gaps of 41.82% and 41.04%, respectively.

Moreover, column "GRB_VI vs. GRB" in "Obj. Values" shows that the valid inequalities are indeed effective and help the GRB to obtain even optimal solutions for test problems $|T| = \{20, 28-40\}$ under the limited computational time of 7200s. The columns "Time (ratio)" present interesting results on the speed of the three methods. As can be seen, the proposed ABD algorithm is always faster than GRB_VI, and GRB_VI is also faster than GRB. It is worth mentioning that the paired comparison between the two methods in terms of speed has been reported whenever both methods have reached the optimal solution under the limited computational time of 7200 s. In summary, for the test problems for which the GRB has obtained the optimal solution (i.e., $|T| = \{4-18, 22-26\}$), it was observed that the ABD is on average 9.88 and 7.46 times faster than GRB and GRB_IV, respectively. Considering the same test problems, the GRB_VI is on average, 1.48 times faster than GRB. It reveals that the valid inequalities accelerate the GRB on the average up to 1.48 times faster when obtaining the optimal solutions. For another set of test problems (i.e., $|T| = \{4-40\}$) where the GRB_VI has obtained the optimal solution, the proposed ABD algorithm is, on average, 12.12 times faster than GRB_VI.

### 6.3. Sensitivity analyses

This section renders a broad range of sensitivity analyses to validate and illustrate the capabilities of the developed model and the solution approach. This investigation encompasses three main categories: 1) sensitivity of the objective function value, the inventory, and the sharing of products to input parameters, 2) inter-relation between the total demand, inventory, and sharing, and 3) capacity utilization. These three categories are explained in the following subsections, respectively.

#### 6.3.1. Sensitivity to input parameters

This section provides a broad sensitivity analysis of three main criteria including the objective function value (i.e., the total cost of

<footnote>* Time limit of 7200 s reached.</footnote>

<footnote>** No feasible solution after 7200 s.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** The impact of the cost of unmet demand on total cost, inventory, and sharing.</figcaption></figure>

the system), the total amounts of inventory, and sharing with respect to changes of certain input parameters including the cost of unmet demand (CU), the demand of products (D), maximum sharing capacity ($G_{max}$), inventory holding cost (ICH), and allowable demand overload ratio ($\mu$).

Figure 4 shows how unmet demand cost affects the three mentioned criteria. Figure 4(a) indicates that if CU increases, the total costs also increase. In fact, with a 35% of increase in the cost of unmet demand, it was observed that supply, manufacturing, distribution, and sharing expenses outweigh the overall cost of unmet demands. However, by increasing this cost from 35% to 60%, the costs of unmet demands surpass the system's other expenses. Finally, by increasing this cost by more than 60%, the system reaches its maximum capacity and achieves a state of equilibrium. At this stage, a slight decrease can be seen in the system's total cost. The reason for this decrease goes back to the fact that after increasing the cost of unmet demand higher than 60%, the last term of the objective function that calculates the total cost of the unmet demand (i.e., $CU \times S_{max}$) becomes the dominant term of the objective function. Accordingly, the model puts a higher effort into minimizing this term. Therefore, a sudden fall happens in the maximum number of unmet demands ($S_{max}$). This decrease is significant

and, overall, leads to a decrease in the total cost of the system.

In terms of the changes in the total inventory of products in the system, Fig. 4(b) illustrates that as the cost of unmet demand increases, the system seeks to minimize mortality by increasing the production of hospital wards-related items such as medical clothes, ordinary beds, ICU beds, and ventilators. As a result, these items' inventory levels increase. Simultaneously, individuals are attempting to take better care of themselves by observing public health, which is related to a rise in the usage of health products, such as masks and hydro-alcoholic gel, which will boost consumption and reduce the remaining inventory of these products. In connection with the changes in the total quantity of product sharing, Fig. 4(c) articulates that the sharing of all products has increased with the increase in the cost of unmet demand. This increase is entirely meaningful for ventilators and medical clothing. Indeed, the system decides to transfer the increased level of product inventory (seen in Fig. 4(b)) via sharing mechanism to fulfill the demand as much as possible to absorb the increase of unmet demand cost. Therefore, the system has attempted to fulfill demand by utilizing all available capacities. In addition, the considerable increase in the sharing of medical clothes is due to the demand for this
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** The impact of demand on total cost, inventory, and sharing.</figcaption></figure>

product being impacted by hospitalized individuals in general and ICU wards.

Figure 5 shows how changes in demand D affect the three criteria. As can be seen in Fig. 5(a), if D increases, the total costs also increase. This trend is logical due to rising supply, production, distribution, and sharing expenses across various facilities. Also, Fig. 5(b) illustrates that the inventory of masks and hydro-alcoholic gels has reduced due to increased consumption of these products in response to rising demands. Other products – which are the essential products for curing patients – are facing an increase in their total inventory, which reflects the system's behavior when the demand increases for such essential products. Indeed, the system produces more and stores more to fulfill the demand for essential products. For the reusable products (i.e., ICU beds and ventilators), the inventory also experiences a special increase, which is the return of these products to the utilization cycle. Finally, it might be claimed from Fig. 5(c) that there is an increasing-then-decreasing oscillation cycle in terms of the total sharing of products. Indeed, initially, the system strives to compensate for its lack of capacity by sharing, but once there is no more capacity to share, the system begins to boost product supply and production, resulting in less sharing.

Figure 6 shows how changes in the maximum sharing capacity $G_{max}$ affect the three concerned criteria. As can be seen in Fig. 6(a), raising the maximum sharing capacity lowers the total cost of the system. This trend has two primary explanations. The first is that as sharing capacity grows, more sharing is allowed through the system (see Fig. 6(c)). Consequently, the supplies and the productions are set up with lower frequency to effectively use the current capacity, mainly from the sharing mechanism. Another reason is that certain hospitals serve as intermediate distribution centers. In reality, due to the geographic structure of each region, hospitals in each region may either acquire products from manufacturers or distribute their excess products in the established sharing mechanism to other hospitals. This issue becomes even more pressing when this procedure lowers the cost of shipping items by eliminating the need to order products from manufacturers and allowing products to be moved between hospitals at significantly cheaper transportation costs (i.e., due to lower distance).

Moreover, as the sharing capacity rises in Fig. 6(b), the total inventory level also climbs because the system can now store more products and share them through the system. Also, as can be observed, the ventilator and ICU beds have greater inventory levels than other products being impacted by two reasons; first, these
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The impact of maximum sharing capacity on total cost, inventory, and sharing.</figcaption></figure>

products are expensive to be produced, so the system attempts to store these products as much as affordable by the system; second; these products may be reused after a predefined lead time. When looking at Fig. 6(c), it is evident that as the capacity for sharing grows; more products are allowed to be shared since sharing is typically less expensive compared to manufacturing, especially for expensive products (e.g., beds and ventilators). In this regard, the bigger the capacity, the less the influence on the sharing of low-cost products (e.g., masks, hydro-alcoholic gels, and medical cloths) because the cost of manufacturing is affordable compared to the cost of sharing for these products.

Figure 7 shows how changes in the inventory holding cost ICH impact the three mentioned criteria. As Fig. 7(a) depicts, the total cost of the system rises monotonically with the increase of ICH. This is because the cost of holding products in hospitals and the objective function have a direct linear connection. What is more in Fig. 7(b), when inventory holding costs rise, more expensive products (e.g., medical clothes, beds, and ventilators) will have lower inventory levels since the ICH is a fraction of the production cost PC (see Table G.8).

Furthermore, for products with a lower production cost, the total inventory level rises. This pattern may be seen in hospitals that

have greater holding costs. Indeed, they would prefer not to keep inventory in the system and instead use the sharing mechanism to satisfy their demand. On the other hand, this behavior enables other hospitals with lower holding costs to keep additional inventory to satisfy their demand. Finally, Fig. 7(a) illustrates that raising the cost of inventory in the hospital lowers product sharing. This is due to hospitals' unwillingness to keep inventory to participate in the sharing process and preferring to get products only as much of their demand directly from the manufacturers.

Figure 8 illustrates the impact of allowable demand overload ratio $\mu$ on the three mentioned criteria. Based on Fig. 8(a), as the $\mu$ increases and consequently more and more amounts of products are allowed to be stored extra than the actual demand, the possibility of fulfilling demand surpluses through the sharing mechanism increases accordingly. As a result, the supply and production setups happen with a lower frequency, and their fixed setup costs decrease consequently. Furthermore, the transportation cost of direct transferring of products from manufacturers to hospitals decreases. All these result in lower total costs of the system. Moreover, Fig. 8(b) reveals that hospitals are more likely to employ the sharing mechanism and to decrease the total inventory of products as the overload ratio rises. As a result, hospitals strive to maintain
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** The impact of inventory holding cost at hospitals on total cost, inventory, and sharing.</figcaption></figure>

<table><caption>**Table 6**

Impact of disease parameters on the demand for products.</caption><thead><tr><th>Demand for ...</th><th>β<sub>1</sub></th><th>β<sub>2</sub></th><th>β<sub>3</sub></th><th>β<sub>4</sub></th><th>ζ</th><th>κ</th><th>η (If ...)</th><th>ρ (If ...)</th><th>ω (If ...)</th><th>c & θ (If ...)</th><th>η<sub>1</sub> & η<sub>2</sub></th></tr></thead><tbody><tr><td>Demand for ...</td><td>β<sub>1</sub></td><td>β<sub>2</sub></td><td>β<sub>3</sub></td><td>β<sub>4</sub></td><td>ζ</td><td>κ</td><td>η (If ...)</td><td>ρ (If ...)</td><td>ω (If ...)</td><td>c & θ (If ...)</td><td>η<sub>1</sub> & η<sub>2</sub></td></tr><tr><td>Masks & Hydro-alcoholic gel</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>—</td><td>— (If β<sub>2</sub> > β<sub>1</sub>)<br/>+ (β<sub>2</sub> > β<sub>1</sub>)</td><td>— (β<sub>1</sub> > β<sub>2</sub>)<br/>+ (β<sub>1</sub> < β<sub>2</sub>)</td><td>— (β<sub>3</sub>, β<sub>4</sub> > β<sub>2</sub>)<br/>+ (β<sub>3</sub>, β<sub>4</sub> < β<sub>2</sub>)</td><td>— (β<sub>4</sub> > β<sub>3</sub>)<br/>+ (β<sub>4</sub> < β<sub>3</sub>)</td><td>+</td></tr><tr><td>Medical clothes</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>+</td><td>—</td><td>—</td></tr><tr><td>Ordinary bed</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>+</td><td>—</td><td>—</td></tr><tr><td>ICU bed & Ventilator</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>+</td><td>—</td><td>—</td></tr></tbody></table>

their inventory levels to a minimum. In addition, when the allowable demand overload ratio rises in Fig. 8(c), the total sharing for all products intends to be increased. Remarkably, this trend is particularly apparent in consumable products because these products cannot be reused. As a result, since the inventory level of these products is derived from both production and sharing amounts, the excess demand is met entirely through the sharing system.

### 6.3.2. Inter-relation between cost, inventory, and sharing

This section seeks the inter-relation between the total demand and the total inventory and sharing of the products in the system. Indeed, the goal is to investigate how different hospitals manage

the inventory level and the sharing of products based on their particular amount of demand. In this regard, Fig. 9(a) depicts the quantity of total demand as well as the proportional level of total inventory and total sharing to the total demand over a period of 3 months (i.e., $|T| = 12$), wherein the results of the even periods have been depicted for the sake of simplicity. As can be observed in Fig. 9(a)(a), there is a modest decrease in the total demand as the pandemic progresses, which is logical since we have only a single peak of demand at the early stages of the pandemic in our case study, and as long as the pandemic progresses, the demand decreases due to the control of the pandemic. In addition, due to a drop in the number of susceptible individuals after the pandemic
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** The impact of allowable demand overload ratio at hospitals on total cost, inventory, and sharing.</figcaption></figure>

peak, the inventory-to-demand ratio in Fig. 9(a)(b) for masks and hydro-alcoholic gels is falling. However, this ratio is higher for ICU beds and ventilators due to the rise in hospitalized individuals in general and ICU wards. As a result, the hospitals have attempted to fulfill the demand by storing more products. Furthermore, the sharing-to-demand ratio in Fig. 9(a)(c) indicates that the hospital is able to engage in the product-sharing procedure owing to the required inventory level.

In order to show how different hospitals manage their inventory and sharing, three different hospitals have been selected that behave differently in managing their inventory and sharing. The hospital of Fig. 9(b) refuses to engage in the sharing mechanism since no input (Fig. 9(b)(c)) and output (Fig. 9(b)(d)) sharing happens in this hospital. Consequently, the inventory level (Fig. 9(b)(b)) of the hospital remains at its minimum for the majority of the products. This indicates that the demand and the supply of products from manufacturers are in an implicit equilibrium. In a second hospital depicted in Fig. 9(c), the situation is different. In fact, this hospital faces an inventory shortage owing to high demand and the possibility of receiving surplus demand. In this situation, the hospital fulfills its demand (Fig. 9(c)(a)) simultaneously from its inventory (Fig. 9(c)(b)) as well as the input sharing

(Fig. 9(c)(c)) coming from other hospitals, and it shares no products with others (Fig. 9(c)(d)). Finally, another hospital is available with excess products as Fig. 9(d). This hospital fulfills its demand (Fig. 9(d)(a)) only from its own inventory (Fig. 9(d)(b)) and shares the excess products with other hospitals (Fig. 9(d)(c)).

### 6.3.3. Supply and production capacity utilization

An interesting investigation in this study has been done in Fig. 10 to see how much of the supply and production capacities have been utilized over the first three months (i.e., $|T| = 12$) of the pandemic.

For this curiosity, Fig. 10(a) depicts the percentage of suppliers that contributed to the fulfillment of the raw materials required to produce the medical products. As can be seen, from the early periods of the pandemic, at least 50% of the suppliers have been used, and this percentage goes up to full utilization of 100% from the sixth week of the pandemic (i.e., the moment that the pandemic approaches its peak). Furthermore, Fig. 10(b) and (c) depict the percentage of manufacturers that contributed to the production of medical products and their capacity percentage utilization, respectively. Knowing that manufacturers are involved in the system with an initial inventory, the whole number and the whole
<figure><img src="image_12.png" /><figcaption>(a) The ratios of total inventory and sharing to demand</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) The ratios of inventory and sharing to demand in a hospital with inbound sharing</figcaption></figure>

<caption>(d) The ratios of inventory and sharing to demand in a hospital with outbound sharing</caption>

<caption>**Fig. 9.** Inter-relation between cost, inventory, and sharing.</caption>

capacity of the manufacturers are not used in the early weeks of the pandemic; however, the involvement of the manufacturers increases as far as the pandemic progresses. Furthermore, Fig. 10(c) shows that the manufacturers do not always utilize their maximum production capacity for certain products (e.g., masks, hydro-alcoholic gels, medical clothes, and ordinary beds). However, they produce a higher rate of the more essential products (i.e., ICU beds and ventilators).

## 7. Managerial insights

Based on the findings of Section 6.3, this section provides a list of managerial insights for the health policy-makers to better plan for possible future health pandemics.

The first managerial implication of the proposed model relates to its ability to act as both business and mortality models. The behavior of the model changes by changing the unmet demand cost. Although the proposed model has a cost objective function, manipulating the unmet demand costs switches the proposed model

from a cost-oriented model to a more mortality-oriented model. In fact, there is a turning point for this switch when increasing the unmet demand, the proposed model goes from cost-oriented planning to mortality-oriented planning. For the current case study, the switch happens when increasing the unmet demand cost up to 35%. Managers should also notice that increasing the unmet demand will also affect the level of inventory and sharing. Therefore, requirements should be met when switching the behavior of the model. In terms of the inventory of products, two different phenomena happen. As the cost of unmet demand increases, the system seeks to minimize mortality by increasing the production of hospital ward-related items such as medical clothes, ordinary beds, ICU beds, and ventilators. As a result, these items' inventory levels increase. On the other hand, when experiencing a high unmet demand cost, individuals are attempting to take better care of themselves by observing public health, which is related to a rise in the usage of health products, such as masks and hydro-alcoholic gel, which will boost consumption and reduce the remaining inventory of these products.
<figure><img src="image_14.png" /><figcaption>Fig. 10. The capacity utilization of suppliers and manufacturers.</figcaption></figure>

Another important observation is the impact of demand varia-
tion on the sharing mechanism. We observed that the amount of
sharing is a concave function of the demand. The results show that
there are concave cycles in the total sharing amount of products
when changing the demand. Initially, when increasing the demand,
the system attempts to mitigate the risk of capacity shortage in
some hospitals through sharing mechanism, but once the maxi-
mum sharing capacity is reached, the system begins to satisfy the
demand by producing and supplying new products, which results
in less sharing. This phenomenon highlights the impact of shar-
ing mechanisms between hospitals to better respond to the surging
product demand. As a matter of fact, sharing mechanisms provide
both cost and responsiveness advantages. In terms of cost, sharing
of products between neighboring hospitals reduces the transporta-
tion costs of transferring newly produced products from manufac-
turers to geographically far hospitals. On the other hand, this shar-
ing mechanism helps to increase the responsiveness and agility of
hospitals to their demands. In fact, due to usually less distance be-
tween neighboring hospitals compared to regional manufacturers,
some hospitals can act as intermediate suppliers. In this situation,
some hospitals should be allowed to store more products to re-
respond to the shortage of farther hospitals more rapidly.

It is worth mentioning that the demand is not an indepen-
dent parameter but a function of disease parameters in the epi-
demiological model. Any changes in these parameters would ei-
ther increase or decrease the demand for products, under spe-
cific conditions. Table 6 shows how increasing the disease pa-
rameters increase ("+") or decreases ("-") the demand for prod-
ucts under specific conditions. Accordingly, any changes in dis-
ease parameters can be mapped to the changes in the objective
function value, inventory, and sharing of products through the de-
mand intermediate parameter. By manipulating the disease param-

eters, health policy-makers can change product demand and at-
tempt better production-inventory-sharing planning.

We observed that as long as the inventory holding cost in-
creases, the inventory level of less critical products (e.g., mask,
hydro-alcoholic gel) increases while the inventory level of critical
products (e.g., beds and ventilators) decreases. It was revealed if
the unit inventory holding cost exceeds the unit of manufacturing
setup cost, the inventory of critical products decreases. In this sit-
uation, the system also discourages the sharing mechanism (due
to the lack of inventory), consequently leading to a higher level
of unmet demand. Therefore, there is an important interplay be-
tween the inventory holding cost and the amount of product shar-
ing. Accordingly, among different cost elements in the system, it is
recommended to keep the inventory holding cost at its minimum
level to encourage the sharing of products between hospitals. With
this observation together with the role of intermediate supply of
some hospitals, policy-makers may provide the geographically cen-
tral hospitals with a higher inventory capacity and less inventory
holding cost, and the farther hospitals (from the manufacturers)
with less inventory holding capacity. This allows for storing more
products in central hospitals with a lower inventory holding cost,
and these hospitals share the excess product with farther ones. In
this case, the distribution system goes from an isolated structure,
where hospitals order and fulfill separately, to a more hierarchical
and collaborative system, where hospitals order and fulfill collabo-
ratively.

8. Conclusions

The COVID-19 pandemic has posed a variety of challenges to
human civilizations, emphasizing the necessity for a decision-
making framework to address them. In this study, we have pre-
sented a novel comprehensive framework that includes three phases of demand controlling, production-distribution-sharing planning, and solution methodology, wherein various specific concerns of the health pandemic have been reflected, particularly the COVID-19 pandemic. In the first phase, we proposed a SEIHRS epidemiological model with the optimal control problem for handling the dynamicity of the demands of consumable and reusable products. In the second phase, we have introduced and formulated a multi-period production-inventory-sharing problem in terms of a mixed-integer linear programming model. Finally, an accelerated Benders decomposition algorithm with a set of tailored valid inequalities has been proposed to solve the proposed model.

Although the framework has been partially adjusted to deal with the COVID-19 pandemic, it can surely be applied to deal with subsequent pandemics as well. The second phase, in particular, is a contrivance for planning supply, manufacturing, distribution, and sharing decisions that may be implemented directly or with minimal adjustments in a variety of populated areas. The computational results show that the proposed accelerated Benders decomposition algorithm is efficient in handling large-sized test problems. In this regard, it was observed that the proposed decomposition method coupled with effective valid inequalities could solve large-sized test problems in a reasonable computational time and 9.88 times faster than the Gurobi solver. In addition, the effectiveness of the valid inequalities was also reported that helping the proposed model to be solved 1.48 times faster using the Gurobi solver. Moreover, the sharing mechanism reduces the total cost of the system and the unmet demand, on average, up to 32.98% and 20.96%, respectively. Furthermore, the two parameters of maximum sharing capacity and allowable demand overload ratio are two valuable factors at hand for managers to regulate inventory and sharing decisions in the proposed sharing mechanism, which significantly influences the system's total cost and unmet demand.

Some ideas for expanding the current study on the prediction and planning stages can be explored. The application of machine learning algorithms to estimate the parameters of the SEIHRS model during the prediction phase can improve demand prediction accuracy. The proposed SEIHRS model is a single-region-based model wherein the evolution of the pandemic in each region is modeled in isolation and no interaction between regions is considered. This assumption is similar to the situation when all regions are under a quarantine restriction and no inter-regional displacement or immigration happens. Therefore, as a new research direction, we propose to develop a multi-regional epidemiological model with inter-regional interactions [8].

In the planning phase, although we have planned the production, distribution, and sharing of the products according to the existing geographical borders for various regions, redistricting models may be applied to define regions, which also improves the efficacy of the sharing mechanism. Furthermore, the cost of unmet demand for all products is treated the same in this study, despite the fact that there is a substantial difference between the lack of ventilators and the absence of masks. As a result, using a cost function rather than a parameter can help improve the planning results.

Last but not least, developing an integrated optimization framework would be another research direction that is worthy of investigation. The current work considers that the disease parameters are constant during the planning horizon and no vaccination takes place. Therefore, a new framework can be developed to incorporate these two aspects. To do so, the proposed optimization model can be integrated with the SEIHRS model in an interactive way to develop a myopic model. In this model, the SEIHRS model and the optimization model interact in each period, and a vaccination rate can be added to the SEIHRS model. Furthermore, from one period to another, disease parameters can also vary.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## CRediT authorship contribution statement

**Behnam Vahdani:** Formal analysis, Methodology, Project administration, Software, Validation, Supervision, Writing – review & editing.
**Mehrdad Mohammadi:** Formal analysis, Methodology, Project administration, Software, Validation, Supervision, Writing – review & editing.
**Simon Thevenin:** Formal analysis, Project administration, Validation, Supervision, Writing – review & editing.
**Patrick Meyer:** Formal analysis, Project administration, Validation, Writing – review & editing, Supervision.
**Alexandre Dolgui:** Formal analysis, Project administration, Validation, Supervision, Writing – review & editing.

## Data availability

Data will be made available on request.

## Appendix A. Basic reproduction number

In what follows, the basic reproduction number $\mathcal{R}_0$ and the endemic equilibrium point are calculated by solving the isoclines, where two equilibrium points can be considered including disease-free equilibrium point $E_0(S_0, 0, 0, 0, 0, 0, 0) = (\frac{\Lambda}{d}, 0, 0, 0, 0, 0, 0)$ and endemic equilibrium point $E^*(S^*, E^*, A^*, I^*, HN^*, HC^*, R^*, Z^*)$. The basic reproduction number $\mathcal{R}_0$ is an indicator for the study of contagious disease dynamics. If $\mathcal{R}_0 \ge 1$, the outbreak is anticipated to continue, and the outbreak is going to end otherwise (i.e., $\mathcal{R}_0 < 1$). In order to calculate $\mathcal{R}_0$, we applied the proposed method by Van den Driessche and Watmough [88]. For this purpose, let $x = (E, A, I, HN, HC)$ and constitute $\frac{dx}{dt}, \mathfrak{S}(x)$, and $v(x)$, where $\mathfrak{S}(x)$ refers to the segments in which new infection terms are introduced, and $v(x)$ refers to the remainder of the segments. Accordingly, we have:

$$ \frac{dx}{dt} = \mathfrak{S}(x) - v(x) \qquad (A.1) $$

$$ \mathfrak{S}(x) = \begin{pmatrix} (\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \qquad (A.2) $$

$$ v(x) = \begin{pmatrix} (k+d)E \\ -k\rho E + (\eta + d + \delta_1)A \\ -k(1-\rho)E - \eta A + (\omega + d + \delta_2)I \\ -\omega(1-c)I + (d + \delta_3 + \gamma_1 + \theta)HN \\ -\omega cI + (d + \delta_4 + \gamma_2)HC - \theta HN \end{pmatrix} \qquad (A.3) $$

In what follows, the equivalent linearized matrices of $\mathfrak{S}(x)$ and $v(x)$ at disease-free equilibrium $E_0 = (S_0, 0, 0, 0, 0, 0, 0)$ can be constructed as (A.4) and (A.5), respectively.

$$ F = (D\mathfrak{S}(x))(E_0) = \begin{pmatrix} 0 & \beta_1 S_0 & \beta_2 S_0 & \beta_3 S_0 & \beta_4 S_0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \qquad (A.4) $$
$$V = (D\nu(x))(E_0) = \begin{pmatrix} k+d & 0 & 0 & 0 & 0 \\ -\kappa\rho & \eta+d+\delta_1 & 0 & 0 & 0 \\ -\kappa(1-\rho) & -\eta & \omega+d+\delta_2 & 0 & 0 \\ 0 & 0 & -\omega(1-c) & d+\delta_3+\gamma_1+\theta & 0 \\ 0 & 0 & -\omega cI & -\theta & d+\delta_4+\gamma_2 \end{pmatrix} \quad (A.5)$$

Consequently, regarding matrices F and V, the basic reproduction number is equal to the special vector of $FV^{-1}$. Moreover, the endemic equilibrium point is obtained by solving the system of Eq. (A.6).

$$\Lambda - (\beta_1 A^* + \beta_2 I^* + \beta_3 HN^* + \beta_4 HC^*) S^* - dS^* + \xi R^* - u_1 k S^* Z^* = 0,$$

$$(\beta_1 A^* + \beta_2 I^* + \beta_3 HN^* + \beta_4 HC^*) S^* - (k+d)E^* = 0,$$

$$\kappa \rho E^* - (\eta + d + \delta_1) A^* = 0,$$

$$\kappa(1-\rho)E^* + \eta A^* - (\omega + d + \delta_2)I^* = 0,$$

$$\omega(1-c)I^* - (d + \delta_3 + \gamma_1 + \theta)HN^* = 0,$$

$$\omega c I^* - (d + \delta_4 + \gamma_2)HC^* + \theta HN^* = 0,$$

$$\gamma_1 HN^* + \gamma_2 HC^* - \xi R^* - dR^* + u_1 k S^* Z^* = 0,$$

$$\frac{pI^*}{1+qI^*} - a_0 Z^* = 0$$

(A.6)

After solving the system of Eq. (A.6), a unique equilibrium point is obtained when the basic reproduction number is greater than one ($R_0 > 1$). However, values smaller than one for the basic reproduction number ($R_0 < 1$) cause the non-occurrence of an endemic equilibrium. Therefore, the following statement can be deduced. The system of equations (1) has a disease-free equilibrium point $E_0(S_0, 0, 0, 0, 0, 0, 0, 0)$ for any values of parameters. And for $R_0 > 1$, the system has a unique endemic equilibrium $E^*(S^*, E^*, A^*, I^*, HN^*, HC^*, R^*, Z^*)$.

## Appendix B. Stability analysis

Local and global stability conditions of the equilibria are analyzed as follows. For this aim, we analyze the local and global stability conditions for the disease-free equilibrium point. Let $\alpha_0 = \kappa + d$, $\alpha_1 = \eta + d + \delta_1$, $\alpha_2 = \omega + d + \delta_2$, $\alpha_3 = d + \delta_3 + \gamma_1 + \theta$, $\alpha_4 = d + \delta_4 + \gamma_2$. The Jacobian matrix $\tilde{J}$ of system (1) is considered as follows:

$$\tilde{J} = \begin{pmatrix} \alpha_{11} & 0 & -\beta_1 S & -\beta_2 S & -\beta_3 S & -\beta_4 S & \zeta & -u_1 \kappa S \\ \beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC & -\alpha_0 & \beta_1 S & \beta_2 S & \beta_3 S & \beta_4 S & 0 & 0 \\ 0 & \kappa \rho & -\alpha_1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \kappa (1-\rho) & \eta & -\alpha_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \omega(1-c) & -\alpha_3 & 0 & 0 & 0 \\ 0 & 0 & 0 & \omega c & \theta & -\alpha_3 & 0 & 0 \\ u_1 k Z & 0 & 0 & 0 & \gamma_1 & \gamma_2 & -(\xi+d) & u_1 k S \\ 0 & 0 & 0 & \frac{p}{(1+qI)^2} & 0 & 0 & 0 & -a_0 \end{pmatrix} \quad (B.1)$$

where $\alpha_{11} = -(\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC) - d - u_1 k Z$.

### B1. Local stability of $E_0$

Jacobian matrix at the disease-free equilibrium point $E_0 = (\Delta/d, 0, 0, 0, 0, 0, 0, 0)$, which is shown by $\tilde{J}E_0$ is constructed and solved with Python to obtain eigenvalues of $\tilde{J}E_0$, where $S_0 = \Delta/d$. In this regard, we have eight eigenvalues, where $\lambda_1 = -d$, $\lambda_2 = -(\xi + d)$, $\lambda_3 = -a_0$, and other five eigenvalues can be obtained from the roots of the following equation:

$$\lambda^5 + Q_1 \lambda^4 + Q_2 \lambda^3 + Q_3 \lambda^2 + Q_4 \lambda + Q_5 = 0 \quad (B.2)$$

Regarding the values of parameters in Table G.8, by solving the above equation with Python, we can conclude that $Q_5$ is positive for $R_0 < 1$. So for $R_0 < 1$, the characteristic equation has roots with negative real parts only when $Q_2, Q_3, Q_4 > 0$, and we have the following theorem.

**Theorem 1.** Disease-free equilibrium point $E_0$ of system (1) is locally asymptotically stable for $R_0 < 1$ when $Q_2, Q_3, Q_4 > 0$.

### B2. Global stability of $E_0$

**Theorem 2.** Disease-free equilibrium point $E_0$ of the system (1) is globally asymptotically stable when $p < d$ and $u_1 k \Delta < a_0 d$ hold.

**Proof.** Let us consider the Lyapunov function as follows:

$$V_1 = (S - S_0 - S_0 \ln(\frac{S}{S_0})) + E + A + I + HN + HC + R + Z \quad (B.3)$$
where $V_1$ is a positive definite function for all $(S, E, A, I, HN, HC, R, Z)$ other than the disease-free equilibrium point. Time derivative of $V_1$ along the solution of the system (1) is calculated as follows:

$$ 
\begin{aligned}
\frac{dV_1}{dt} &= (1 - \frac{S_0}{S})[\Lambda - (\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S - dS + \xi R - u_1(t)kSZ] \\
&\quad + (\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S - d(E + A + I + HN + HC + R) - \delta_1 A - \delta_2 I \\
&\quad - \delta_3 HN - \delta_4 HC - \xi R + u_1 kSZ + \frac{pl}{1+ql} - a_0 Z \\
&= \Delta - d(E + A + I + HN + HC + R) - \delta_1 A - \delta_2 I - \delta_3 HN - \delta_4 HC - \frac{\Delta S_0}{S} + dS_0 - \frac{\xi S_0 R}{S} + u_1 kS_0 Z \\
&\quad + \frac{pl}{1+ql} - a_0 Z \leq (\frac{\Delta}{S} - d)(S - S_0) + (p - d)I + \frac{1}{d}(u_1 k\Delta - a_0 d)Z \\
&= -\frac{d}{S}(S - S_0)^2 + (p - d)I + \frac{1}{d}(u_1 k\Delta - a_0 d)Z
\end{aligned}
\quad (B.4) $$

Therefore, $\frac{dV_1}{dt} < 0$ when $p < d$ and $u_1 k\Delta < a_0 d$ hold. Also, $\frac{dV_1}{dt} = 0$ when $S = S_0$ and $E = A = I = HN = HC = R = Z = 0$. Regarding LaSalle's invariance principle [89], $E_0$ is globally asymptotically stable when $R_0 < 1$, and when the considered parametric limitations are satisfied. Similarly, the local and global stability of $E^*$ can be examined. $\square$

## Appendix C. Proof of Theorem 3

**Proof.** The optimal control interventions exist when the following conditions are satisfied:

* The solution space of system dynamics (1) with control variables $u_i, i = 1, \dots, 5$ in $\cap \neq \emptyset$.

* The state system is a linear function of control variables, whose coefficients depend on state variables and time, where set $\cap$ is convex and closed.

* Integrand of (3): $L$ is convex on $\cap$ and $L(S, E, A, I, HN, HC, R, u_1, u_2, u_3, u_4, u_5) \geq f(u_1, u_2, u_3, u_4, u_5)$, where $f(u_1, u_2, u_3, u_4, u_5)$ is continuous, and $\|(u_1, u_2, u_3, u_4, u_5)\|^{-1}f(u_1, u_2, u_3, u_4, u_5) \to \infty$, when $\|(u_1, u_2, u_3, u_4, u_5)\| \to \infty$. Indeed, $\|\cdot\|$ indicates the $L^5(0, T_f)$ norm.

Regarding system dynamics (1), the total population is $N = S + E + A + I + HN + HC + R$. Hence, $\frac{dN}{dt} = \Lambda - dN - \delta_1 A - \delta_2 I - \delta_3 HN - \delta_4 HC \leq \Lambda - dN$, and $0 < N(t) \leq N(0)e^{-dt} + \frac{\Lambda}{d}(1 - e^{-dt})$, where $N(0) = S(0) + E(0) + A(0) + I(0) + HN(0) + HC(0) + R(0)$. Also, when $t \to \infty$, we have $0 < N(t) \leq \frac{\Lambda}{d}$. In addition, $\frac{dZ}{dt} = \frac{pl}{1+ql} - a_0 Z \leq pl - a_0 Z$. So, $\frac{dZ}{dt} + a_0 Z \leq \frac{p\Lambda}{d}$, and $0 < Z(t) \leq Z(0)e^{-a_0 t} + \frac{p\Lambda}{a_0 d}(1 - e^{-a_0 t})$. Also, when $t \to \infty$, we have $0 < Z(t) \leq \frac{p\Lambda}{a_0 d}$.

The solution of system dynamics (1) is bounded for each of the control variables in $\cap$, and the right-hand side functions are also locally Lipschitzian. The Picard Lindelöf theorem demonstrates that the first condition is met [90].

Regarding the definition, the control set $\cap$ is convex and closed. The second condition is likewise met because all of the equations in (1) are represented as linear equations in $u_i^*, i = 1, \dots, 5$, whose state variables are dependent on coefficients. Furthermore, the convex property of integrand $L(S, E, A, I, HN, HC, R, u_1, u_2, u_3, u_4, u_5)$ is guaranteed by the quadratic nature of all control variables. Also, we have $L(S, E, A, I, HN, HC, R, u_1, u_2, u_3, u_4, u_5) = w_1 I + w_2 HN + w_3 HC + w_4 u_1^2 + w_5 u_2^2 + w_6 u_3^2 + w_7 u_4^2 + w_8 u_5^2 \geq w_4 u_1^2 + w_5 u_2^2 + w_6 u_3^2 + w_7 u_4^2 + w_8 u_5^2$.

Let $\bar{u} = \min(w_4, w_5, w_6, w_7, w_8) > 0$ and $f(u_1, u_2, u_3, u_4, u_5) = \bar{u}(u_1^2 + u_2^2 + u_3^2 + u_4^2 + u_5^2)$. Hence,

$L(S, E, A, I, HN, HC, R, Z, u_1, u_2, u_3, u_4, u_5) \geq f(u_1, u_2, u_3, u_4, u_5)$

In this regard, $f$ is continuous and $\|(u_1, u_2, u_3, u_4, u_5)\|^{-1}f(u_1, u_2, u_3, u_4, u_5) \to \infty$, when $\|(u_1, u_2, u_3, u_4, u_5)\| \to \infty$. Hence, the third condition is also met. Therefore, we can conclude that there are control variables $u_i^*, i = 1, \dots, 5$ with the condition $J(u_1^*, u_2^*, u_3^*, u_4^*, u_5^*) = \min[J(u_1, u_2, u_3, u_4, u_5)]$ [38,91]. $\square$

## Appendix D. Optimal control variables

In what follows, in order to obtain optimal control variables, the Hamiltonian function is formulated by introducing adjoint variables $\lambda = (\lambda_1, \lambda_2, \dots, \lambda_8) \in \mathbb{R}^8$, and minimized by applying Pontryagin's Maximum Principle [92], which has been explained in Appendix E.

**Theorem 4.** If $S^*, E^*, A^*, I^*, HN^*, HC^*, R^*, Z^*$ are optimal state variables, and $u_i^*, i = 1, \dots, 5$ are optimal control variables of system (3) and (4), there are adjoint variables $\lambda = (\lambda_1, \lambda_2, \dots, \lambda_8) \in \mathbb{R}^8$ that satisfy the canonical system of Eq. (D.1):

$$ \frac{d\lambda_1}{dt} = \lambda_1(\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC + d + u_1 kZ) - \lambda_2(\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC) - \lambda_7 u_1 kZ, $$

$$ \frac{d\lambda_2}{dt} = \lambda_2(k+d) - \lambda_3(\kappa\rho) - \lambda_4(\kappa(1-\rho)), $$

$$ \frac{d\lambda_3}{dt} = \lambda_1(\beta_1 S) - \lambda_2(\beta_1 S) + \lambda_3(\eta + d + \delta_1) - \lambda_4(\eta), $$

$$ \frac{d\lambda_4}{dt} = -w_1 + \lambda_1(\beta_2 S) - \lambda_2(\beta_2 S) + \lambda_4\left(\omega + d + \delta_2 + \frac{\varepsilon_1 u_2}{(1 + \zeta_1 I)^2} + \frac{\varepsilon_2 u_3}{(1 + \zeta_2 I)^2}\right) - $$
$$ \lambda_5 \left( \omega(1-c) + \frac{\varepsilon_1 u_2}{(1 + \zeta_1 I)^2} \right) - \lambda_6 \left( \omega c + \frac{\varepsilon_2 u_3}{(1 + \zeta_2 I)^2} \right) - \lambda_8 \left( \frac{p}{(1 + qI)^2} \right), $$

$$ \frac{d\lambda_5}{dt} = -w_2 + \lambda_1(\beta_3 S) - \lambda_2(\beta_3 S) + \lambda_5 \left( d + \delta_3 + \gamma_1 + \theta + \frac{\varepsilon_3 u_4}{(1 + \zeta_3 HN)^2} \right) - \qquad (D.1) $$

$$ \lambda_6(\theta) - \lambda_7 \left( \gamma_1 + \frac{\varepsilon_3 u_4}{(1 + \zeta_3 HN)^2} \right), $$

$$ \frac{d\lambda_6}{dt} = -w_3 + \lambda_1(\beta_4 S) - \lambda_2(\beta_4 S) + \lambda_6 \left( d + \delta_4 + \gamma_2 + \frac{\varepsilon_4 u_5}{(1 + \zeta_4 HC)^2} \right) - \lambda_7 \left( \gamma_2 + \frac{\varepsilon_4 u_5}{(1 + \zeta_4 HC)^2} \right), $$

$$ \frac{d\lambda_7}{dt} = -\lambda_1(\xi) + \lambda_7(\xi + d), $$

$$ \frac{d\lambda_8}{dt} = \lambda_1(u_1 k S) - \lambda_7(u_1 k S) + \lambda_8(a_0) $$

Regarding the transversality conditions $\lambda_i(T_f) = 0, i = 1, \dots, 8$, optimal control variables $u_i^*$, $i = 1, \dots, 5$ are obtained as:

$$ u_1^* = \min \left\{ \max \left\{ 0, \frac{kS^*Z^*(\lambda_1 - \lambda_7)}{2w_4} \right\}, 1 \right\}, \qquad (D.2) $$

$$ u_2^* = \min \left\{ \max \left\{ 0, \frac{\varepsilon_1 I^*}{1 + \zeta_1 I^*} (\lambda_4 - \lambda_5) \right\}, 1 \right\}, \qquad (D.3) $$

$$ u_3^* = \min \left\{ \max \left\{ 0, \frac{\varepsilon_2 I^*}{1 + \zeta_2 I^*} (\lambda_4 - \lambda_6) \right\}, 1 \right\}, \qquad (D.4) $$

$$ u_4^* = \min \left\{ \max \left\{ 0, \frac{\varepsilon_3 HN^*}{1 + \zeta_3 HN^*} (\lambda_5 - \lambda_7) \right\}, 1 \right\}, \qquad (D.5) $$

$$ u_5^* = \min \left\{ \max \left\{ 0, \frac{\varepsilon_4 HC^*}{1 + \zeta_4 HC^*} (\lambda_6 - \lambda_7) \right\}, 1 \right\}, \qquad (D.6) $$

**Proof.** Let $S^*, E^*, A^*, I^*, HN^*, HC^*, R^*, Z^*$ be optimal state variables, and $u_i^*$, $i = 1, \dots, 5$ be optimal control variables of system (1) that minimize objective function (3). With the help of Pontryagin's Maximum Principle and defined adjoint variables $\lambda_i, i = 1, \dots, 8$, we have $\frac{d\lambda_1}{dt} = -\frac{\partial \tilde{H}}{\partial S}, \frac{d\lambda_2}{dt} = -\frac{\partial \tilde{H}}{\partial E}, \frac{d\lambda_3}{dt} = -\frac{\partial \tilde{H}}{\partial A}, \frac{d\lambda_4}{dt} = -\frac{\partial \tilde{H}}{\partial I}, \frac{d\lambda_5}{dt} = -\frac{\partial \tilde{H}}{\partial HN}, \frac{d\lambda_6}{dt} = -\frac{\partial \tilde{H}}{\partial HC}, \frac{d\lambda_7}{dt} = -\frac{\partial \tilde{H}}{\partial R}, \frac{d\lambda_8}{dt} = -\frac{\partial \tilde{H}}{\partial Z}$, which lead to obtain Eq. (D.1). In regard to the transversality conditions $\lambda_i(T_f) = 0, i = 1, \dots, 8$ and optimality conditions, we have:

$$ \frac{\partial \tilde{H}}{\partial u_1} |_{u_1} = u_1^* = 0 $$ $$ \frac{\partial \tilde{H}}{\partial u_1} |_{u_1} = u_1^* = 2w_4 u_1^* - \lambda_1 kSZ + \lambda_7 kSZ = 0 \qquad (D.7) $$

$$ u_1^* = \frac{kS^*Z^*(\lambda_1 - \lambda_7)}{2w_4} $$

$$ \frac{\partial \tilde{H}}{\partial u_2} |_{u_2} = u_2^* = 0 $$ $$ \frac{\partial \tilde{H}}{\partial u_2} |_{u_2} = u_2^* = 2w_5 u_2^* + \lambda_4(-\varepsilon_1 I/1 + \zeta_1 I) + \lambda_5(\varepsilon_1 I/1 + \zeta_1 I) = 0 \qquad (D.8) $$

$$ u_2^* = \frac{\left( \frac{\varepsilon_1 I^*}{1 + \zeta_1 I^*} \right) (\lambda_4 - \lambda_5)}{2w_5} $$

$$ \frac{\partial \tilde{H}}{\partial u_3} |_{u_3} = u_3^* = 0 $$ $$ \frac{\partial \tilde{H}}{\partial u_3} |_{u_3} = u_3^* = 2w_6 u_3^* + \lambda_4(-\varepsilon_2 I/1 + \zeta_2 I) + \lambda_6(\varepsilon_2 I/1 + \zeta_2 I) = 0 \qquad (D.9) $$
$$u_3^* = \frac{\left( \frac{\varepsilon_2 I^*}{1 + \zeta_2 I^*} \right) (\lambda_4 - \lambda_6)}{2w_6}$$

$$\frac{\partial \bar{H}}{\partial u_4} |_{u_4} = u_4^* = 0$$

$$\frac{\partial \bar{H}}{\partial u_4} |_{u_4} = u_4^* = 2w_7 u_4^* + \lambda_5 (-\varepsilon_3 HN/1 + \zeta_3 HN) + \lambda_7 (\varepsilon_3 HN/1 + \zeta_3 HN) = 0 \quad (D.10)$$

$$u_4^* = \frac{\left( \frac{\varepsilon_3 HN^*}{1 + \zeta_3 HN^*} \right) (\lambda_5 - \lambda_7)}{2w_7}$$

$$\frac{\partial \bar{H}}{\partial u_5} |_{u_5} = u_5^* = 0$$

$$\frac{\partial \bar{H}}{\partial u_5} |_{u_5} = u_5^* = 2w_8 u_5^* + \lambda_6 (-\varepsilon_4 HC/1 + \zeta_4 HC) + \lambda_7 (\varepsilon_4 HC/1 + \zeta_4 HC) = 0 \quad (D.11)$$

$$u_5^* = \frac{\left( \frac{\varepsilon_4 HC^*}{1 + \zeta_4 HC^*} \right) (\lambda_6 - \lambda_7)}{2w_8}$$

Now, we can deduce the following optimal control variables from these observations and the properties of control set $\mathcal{N}$:

$$u_1^* = \begin{cases} 0 & \text{if } \frac{kS^*Z^*(\lambda_1 - \lambda_7)}{2w_4} < 0 \\ [5pt] \frac{kS^*Z^*(\lambda_1 - \lambda_7)}{2w_4} & \text{if } 0 \le \frac{kS^*Z^*(\lambda_1 - \lambda_7)}{2w_4} \le 1 \\ [5pt] 1 & \text{if } \frac{kS^*Z^*(\lambda_1 - \lambda_7)}{2w_4} > 1 \end{cases} \quad (D.12)$$

$$u_2^* = \begin{cases} 0 & \text{if } \frac{\left( \frac{\varepsilon_1 I^*}{1 + \zeta_1 I^*} \right) (\lambda_4 - \lambda_5)}{2w_5} < 0 \\ [5pt] \frac{\left( \frac{\varepsilon_1 I^*}{1 + \zeta_1 I^*} \right) (\lambda_4 - \lambda_5)}{2w_5} & \text{if } 0 \le \frac{\left( \frac{\varepsilon_1 I^*}{1 + \zeta_1 I^*} \right) (\lambda_4 - \lambda_5)}{2w_5} \le 1 \\ [5pt] 1 & \text{if } \frac{\left( \frac{\varepsilon_1 I^*}{1 + \zeta_1 I^*} \right) (\lambda_4 - \lambda_5)}{2w_5} > 1 \end{cases} \quad (D.13)$$

$$u_3^* = \begin{cases} 0 & \text{if } \frac{\left( \frac{\varepsilon_2 I^*}{1 + \zeta_2 I^*} \right) (\lambda_4 - \lambda_6)}{2w_6} < 0 \\ [5pt] \frac{\left( \frac{\varepsilon_2 I^*}{1 + \zeta_2 I^*} \right) (\lambda_4 - \lambda_6)}{2w_6} & \text{if } 0 \le \frac{\left( \frac{\varepsilon_2 I^*}{1 + \zeta_2 I^*} \right) (\lambda_4 - \lambda_6)}{2w_6} \le 1 \\ [5pt] 1 & \text{if } \frac{\left( \frac{\varepsilon_2 I^*}{1 + \zeta_2 I^*} \right) (\lambda_4 - \lambda_6)}{2w_6} > 1 \end{cases} \quad (D.14)$$

$$u_4^* = \begin{cases} 0 & \text{if } \frac{\left( \frac{\varepsilon_3 HN^*}{1 + \zeta_3 HN^*} \right) (\lambda_5 - \lambda_7)}{2w_7} < 0 \\ [5pt] \frac{\left( \frac{\varepsilon_3 HN^*}{1 + \zeta_3 HN^*} \right) (\lambda_5 - \lambda_7)}{2w_7} & \text{if } 0 \le \frac{\left( \frac{\varepsilon_3 HN^*}{1 + \zeta_3 HN^*} \right) (\lambda_5 - \lambda_7)}{2w_7} \le 1 \\ [5pt] 1 & \text{if } \frac{\left( \frac{\varepsilon_3 HN^*}{1 + \zeta_3 HN^*} \right) (\lambda_5 - \lambda_7)}{2w_7} > 1 \end{cases} \quad (D.15)$$
$$ u_5^* = \begin{cases} 0 & \text{if } \frac{(\frac{\varepsilon_4 HC^*}{1 + \zeta_4 HC^*}) (\lambda_6 - \lambda_7)}{2w_8} < 0 \\ [5pt] \frac{(\frac{\varepsilon_4 HC^*}{1 + \zeta_4 HC^*}) (\lambda_6 - \lambda_7)}{2w_8} & \text{if } 0 \le \frac{(\frac{\varepsilon_4 HC^*}{1 + \zeta_4 HC^*}) (\lambda_6 - \lambda_7)}{2w_8} \le 1 \\ [5pt] 1 & \text{if } \frac{(\frac{\varepsilon_4 HC^*}{1 + \zeta_4 HC^*}) (\lambda_6 - \lambda_7)}{2w_8} > 1 \end{cases} \quad (D.16) $$

which is corresponding to (D.2)-(D.6). □

As mentioned earlier, Pontryagin's Maximum Principle helps to obtain optimal controls $(u_1^*, u_2^*, u_3^*, u_4^*, u_5^*)$. In this regard, the impacts of the application of one or all control variables to obtain the minimal cost have been examined one by one. So, the control system in Eqs. (3) and (4) is solved with respect to the mentioned parameters. On this subject, a number of combinations of control variables are reflected, including $[u_1, (u_1, u_2), (u_1, u_3), (u_1, u_4), (u_1, u_5), (u_1, u_2, u_3, u_4, u_5)]$. The computational simulation for all combinations is executed by Python. Next, the Forward-backward sweep method is employed to find the optimal control variables, where the optimal and adjoint state systems are solved by forward and backward in time, respectively. Afterward, in order to update the optimal controls by Hamiltonian for the optimality of the system, the steepest descent method is employed. This procedure carries on until the convergence [93].

## Appendix E. Pontryagin's maximum principle

**Proof.** We derived the necessary conditions for the optimal control problem for system (3) and (4) by using Pontryagin's Maximum Principle [91,92]. Let $\tilde{H}$ be the Hamiltonian function, which is defined as (E.1) and (E.2):

$$ \begin{aligned} \tilde{H}(S, E, A, I, HN, HC, R, Z, u_1, u_2, u_3, u_4, u_5, \lambda) \\ = L(S, E, A, I, HN, HC, R, Z, u_1, u_2, u_3, u_4, u_5) \\ + \lambda_1 \frac{dS}{dt} + \lambda_2 \frac{dE}{dt} + \lambda_3 \frac{dA}{dt} + \lambda_4 \frac{dI}{dt} + \lambda_5 \frac{dHN}{dt} + \lambda_6 \frac{dHC}{dt} + \lambda_7 \frac{dR}{dt} + \lambda_8 \frac{dZ}{dt} \end{aligned} \quad (E.1) $$

Hence,

$$ \begin{aligned} \tilde{H} ={}& w_1 I + w_2 HN + w_3 HC + w_4 u_1^2 + w_5 u_2^2 + w_6 u_3^2 + w_7 u_4^2 + w_8 u_5^2 \\ & + \lambda_1 [\Lambda - (\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S - dS + \xi R - u_1(t)kSZ] \\ & + \lambda_2 [(\beta_1 A + \beta_2 I + \beta_3 HN + \beta_4 HC)S - (k+d)E] \\ & + \lambda_3 [\kappa \rho E - \eta A - (d+\delta_1)A] \\ & + \lambda_4 [\kappa (1-\rho)E + \eta A - \omega I - (d+\delta_2)I - \frac{\varepsilon_1 u_2(t)I}{1+\zeta_1 I} - \frac{\varepsilon_2 u_3(t)I}{1+\zeta_2 I}] \\ & + \lambda_5 [\omega(1-c)I - (d+\delta_3)HN - \gamma_1 HN - \theta HN + \frac{\varepsilon_1 u_2(t)I}{1+\zeta_1 I} - \frac{\varepsilon_3 u_4(t)HN}{1+\zeta_3 HN}] \\ & + \lambda_6 [\omega c I - (d+\delta_4)HC - \gamma_2 HC + \theta HN + \frac{\varepsilon_2 u_3(t)I}{1+\zeta_2 I} - \frac{\varepsilon_4 u_5(t)HC}{1+\zeta_4 HC}] \\ & + \lambda_7 [\gamma_1 HN + \gamma_2 HC - \xi R - dR + u_1 kSZ + \frac{\varepsilon_3 u_4(t)HN}{1+\zeta_3 HN} + \frac{\varepsilon_4 u_5(t)HC}{1+\zeta_4 HC}] \\ & + \lambda_8 [\frac{pI}{1+qI} - a_0 Z] \end{aligned} \quad (E.2) $$

## Appendix F. The PISP Benders reformulation

Let $\tilde{w}_{jt}$, $\tilde{y}_{pmt}$, and $\tilde{v}_m$ be the vectors of fixed $w_{jt}$, $y_{pmt}$, and $v_m$ variables, respectively. Based on the notation of Table F.7, the Benders sub-problem can be written as follows:

$$ \begin{aligned} \min Z ={}& \sum_{p \in P} \sum_{j \in J} \sum_{m \in K} \sum_{t \in T} TC_{pjm} u_{pjmt} + \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} PC_{pm} x_{pmt} + \sum_{p \in P} \sum_{m \in K} \sum_{m' \in H} \sum_{t \in T} TC_{pmm'} q_{mm'}^{pt} \\ & + \sum_{p \in P} \sum_{m \in H} \sum_{m' \in H} \sum_{t \in T} TC_{pmm'} s_{mm'}^{pt} + \sum_{p \in P} \sum_{m \in K \cup H} \sum_{t \in T} IC_{pm} I_{pmt} + CU \times S_{max} \end{aligned} \quad (F.1) $$

s.t.: Constraints (6), (8), (10)-(22), (25)-(29), and

$$ \sum_{m \in K} u_{pjmt} \le U_{pj} \tilde{w}_{jt} \quad \forall j \in J, p \in P, t \in T \quad (F.2) $$
<table><caption>**Table F.7**

Vectors of the dual variables associated with the Benders sub-problem constraints.</caption><thead><tr><th>Not.</th><th>Desc.</th><th>Not.</th><th>Desc.</th></tr></thead><tbody><tr><td>α</td><td>{α<sub>pmt</sub> ≥ 0 | p ∈ P, m ∈ H, t ∈ T}</td><td>β</td><td>{β<sub>jpt</sub> ≥ 0 | j ∈ J, p ∈ P, t ∈ T}</td></tr><tr><td>Φ</td><td>{Φ<sub>pmt</sub> ∈ ℝ | p ∈ P, m ∈ H, t ∈ T}</td><td>η</td><td>{η<sub>pmt</sub> ≥ 0 | p ∈ P, m ∈ K, t ∈ T}</td></tr><tr><td>λ<sub>1</sub></td><td>{λ<sub>ipm,1</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ K<sub>i</sub>(m)}</td><td>λ</td><td>{λ<sub>ipm</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ K<sub>i</sub>(m), t ∈ T; t > 1}</td></tr><tr><td>τ</td><td>{τ<sub>pmt</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td><td>π<sub>1</sub></td><td>{π<sub>ipm,1</sub> ∈ ℝ | p ∈ P, i ∈ I, m ∈ K<sub>i</sub>(m)}</td></tr><tr><td>π</td><td>{π<sub>ipmt</sub> ∈ ℝ | p ∈ P, i ∈ I, m ∈ K<sub>i</sub>(m), t ∈ T; t > 1}</td><td>ε</td><td>{ε ∈ pmt ≥ 0 | p ∈ P, m ∈ M, t ∈ T}</td></tr><tr><td>φ<sub>1</sub></td><td>{φ<sub>ipm,1</sub> ∈ ℝ | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m)}</td><td>φ</td><td>{φ<sub>ipmt</sub> ∈ ℝ | p ∈ P', i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T; t > 1}</td></tr><tr><td>φ'</td><td>{φ'<sub>ipmt</sub> ∈ ℝ | p ∈ P'', i ∈ I, m ∈ H<sub>i</sub>(m); 1 < t ≤ L}</td><td>φ''</td><td>{φ''<sub>ipmt</sub> ∈ ℝ | p ∈ P'', i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T; t > L}</td></tr><tr><td>Ω</td><td>{Ω<sub>pmt</sub> ≥ 0 | p ∈ P, m ∈ H, t ∈ T}</td><td>ω</td><td>{ω<sub>pt</sub> ≥ 0 | p ∈ P, t ∈ T}</td></tr><tr><td>χ</td><td>{χ<sub>ipmt</sub> ≥ 0 | p ∈ P, i ∈ I, φ' =, m ∈ H<sub>i</sub>(m), t ∈ T}</td><td>Υ</td><td>{Υ<sub>ipmt</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td></tr><tr><td>W</td><td>{W<sub>ipmt</sub> ∈ ℝ | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td><td>ℝ</td><td>{ℝ<sub>ipmt</sub> ∈ ℝ | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td></tr><tr><td>V</td><td>{V<sub>ipmt</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td><td>V'</td><td>{V'<sub>ipmt</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td></tr><tr><td>Z</td><td>{Z<sub>ipmt</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td><td>Z'</td><td>{Z'<sub>ipmt</sub> ≥ 0 | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td></tr><tr><td>Z''</td><td>{Z''<sub>ipmt</sub> ∈ ℝ | p ∈ P, i ∈ I, m ∈ H<sub>i</sub>(m), t ∈ T}</td><td>N</td><td>{N<sub>ipmm't</sub> ≥ 0 | p ∈ P, i ∈ I, m, m' ∈ H<sub>i</sub>(m), t ∈ T; m ≠ m'}</td></tr><tr><td>N'</td><td>{N'<sub>ipmm't</sub> ≥ 0 | p ∈ P, i ∈ I, m, m' ∈ H<sub>i</sub>(m), t ∈ T; m ≠ m'}</td><td>M</td><td>{M<sub>ipmm't</sub> ∈ ℝ | p ∈ P, i ∈ I, m, m' ∈ H<sub>i</sub>(m), t ∈ T; m ≠ m'}</td></tr><tr><td>M'</td><td>{M'<sub>ipmm't</sub> ∈ ℝ | p ∈ P, i ∈ I, m, m' ∈ H<sub>i</sub>(m), t ∈ T; m ≠ m'}</td><td>M''</td><td>{M''<sub>ipmm't</sub> ≥ 0 | p ∈ P, i ∈ I, m, m' ∈ H<sub>i</sub>(m), t ∈ T; m ≠ m'}</td></tr><tr><td>H</td><td>{H<sub>ipmm't</sub> ≥ 0 | p ∈ P, i ∈ I, m, m' ∈ H<sub>i</sub>(m), t ∈ T; m ≠ m'}</td><td></td><td></td></tr></tbody></table>

$$x_{pmt} \leq Q_{pm} \bar{y}_{pmt} \quad \forall p \in P, m \in K, t \in T \qquad (F.3)$$

$$b_{pmt}^{+} \leq T Q_p \bar{v}_m \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \qquad (F.4)$$

$$b_{pmt}^{-} \geq T Q_p (\bar{v}_m - 1) \quad \forall i \in I, p \in P, m \in H_i(m), t \in T \qquad (F.5)$$

The Benders master problem can be written as follows:

$$\min Z = \sum_{j \in J} \sum_{t \in T} OC_j W_{jt} + \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} SC_{pm} y_{pmt} + \zeta \qquad (F.6)$$

s.t.: Constraints (30)–(35) and (36), and

$$\begin{align}
& - \sum_{p \in P} \sum_{m \in H} \sum_{t \in T} D_{pmt} \alpha_{pmt} - \sum_{j \in J} \sum_{m \in H} \sum_{t \in T} U_{pj} W_{jt} \beta_{jpt} - \sum_{p \in P} \sum_{m \in K} \sum_{t \in T} Q_{pm} y_{pmt} \eta_{pmt} \nonumber \\
& + \sum_{i \in I} \sum_{m \in K_i(m)} \sum_{p \in P} I_{pm}^0 (\pi_{ipm,1} - \lambda_{ipm,1}) - \sum_{p \in P} \sum_{m \in M} \sum_{t \in T} I_{pm}^{\max} \epsilon_{pmt} \nonumber \\
& + \sum_{i \in I} \sum_{p \in P} \sum_{m \in H_i(m)} I_{pm}^0 \phi_{ipm,1} + \sum_{p \in P} \sum_{m \in H} \sum_{t \in T} D_{pmt} \Omega_{pmt} - \sum_{p \in P} \sum_{t \in T} T Q_p \omega_{pt} \nonumber \\
& - \sum_{i \in I} \sum_{p \in P} \sum_{m' \in H_i(m)} \sum_{t \in T} G_{pm'}^{\max} \chi_{ipm't} + \sum_{i \in I} \sum_{p \in P} \sum_{m' \in H_i(m)} \sum_{t \in T} D_{pmt} W_{ipmt} \nonumber \\
& - \sum_{i \in I} \sum_{p \in P} \sum_{m' \in H_i(m)} \sum_{t \in T} T Q_p \nu_m Z_{ipmt} + \sum_{i \in I} \sum_{p \in P} \sum_{m' \in H_i(m)} \sum_{t \in T} T Q_p (\nu_m - 1) Z'_{ipmt} \nonumber \\
& - \sum_{i \in I} \sum_{p \in P} \sum_{m \in H_i(m)} \sum_{t \in T} G_{pm}^{\max} \mathcal{H}_{ipmm't} \leq \zeta \nonumber \\
& m' \neq m \nonumber
\end{align}$$

$$\forall (\alpha, \beta, \Phi, \eta, \lambda_1, \lambda, \tau, \pi_1, \pi, \epsilon, \phi_1, \phi, \phi', \phi'', \Omega, \omega, \chi, \mathcal{Y}, \mathcal{W}, \mathcal{R}, \mathcal{V}, \mathcal{Z}', \mathcal{Z}, \mathcal{Z}', \mathcal{Z}'', \mathcal{N}, \mathcal{N}', \mathcal{M}, \mathcal{M}', \mathcal{M}'', \mathcal{H}) \in P_{\nabla} \qquad (F.8)$$

where $\nabla$ signifies the polyhedron specified by the constraints of the problem, and $P_{\nabla}$ implicates the set of extreme points of $\nabla$. Constraint (F.7) represents the optimality cut which can be generated consecutively in accordance with Bender's sub-problem solution and the vectors of the dual variables.

## Appendix G. Details of the case study's data

Detailed information on all other parameters in the SEIHRS model and PISP formulation is provided in Table G.8. The majority of the parameters are based on real data except for the setup costs and the supply and production capacities which have been estimated based on expert opinion. Except for the death rates as well as the population in each set of populated areas, other parameters of the SEIHRS model are considered similar for the whole of France [52,94,95].
<caption>**Table G.8**</caption>

Input data.

<table><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>I</td><td>12</td><td>γ</td><td>[0.1,0.5]</td><td>ω</td><td>0.02</td></tr><tr><td>J</td><td>12</td><td>OC (ε)</td><td>[15,000, 20,000]</td><td>γ<sub>1</sub></td><td>0.5</td></tr><tr><td>M</td><td>233</td><td>L (week)</td><td>2</td><td>γ<sub>2</sub></td><td>0.3</td></tr><tr><td>K</td><td>33</td><td>β<sub>1</sub></td><td>3.70E-11</td><td>p</td><td>0.01</td></tr><tr><td>H</td><td>200</td><td>β<sub>2</sub></td><td>1.48E-10</td><td>q</td><td>1</td></tr><tr><td>P</td><td>6</td><td>β<sub>3</sub></td><td>7.00E-12</td><td>a<sub>0</sub></td><td>0.06</td></tr><tr><td>PC (ε)</td><td>[0.5, 20, 000]</td><td>β<sub>4</sub></td><td>7.00E-12</td><td>c</td><td>0.2</td></tr><tr><td>SC (ε)</td><td>[500, 5M]</td><td>ζ</td><td>3.00E-10</td><td>ε<sub>1</sub></td><td>0.1</td></tr><tr><td>TC (ε /100km)</td><td>[5%, 10%] × PC</td><td>κ</td><td>0.6</td><td>ε<sub>2</sub></td><td>0.1</td></tr><tr><td>IC (ε)</td><td>[5%, 20%] × PC</td><td>k</td><td>0.002</td><td>ε<sub>3</sub></td><td>0.9</td></tr><tr><td>Q(#)</td><td>[100, 1M]</td><td>φ</td><td>0.1</td><td>ε<sub>4</sub></td><td>0.9</td></tr><tr><td>I<sup>max</sup> (m<sub>me</sub>E<sub>K</sub>) (#)</td><td>[20, 10M]</td><td>η</td><td>0.15</td><td>ζ<sub>1</sub></td><td>0.05</td></tr><tr><td>I<sup>max</sup> (m<sub>me</sub>H) (#)</td><td>[4, 10,000]</td><td>δ<sub>1</sub></td><td>1.00E-04</td><td>ζ<sub>2</sub></td><td>0.05</td></tr><tr><td>U(#)</td><td>[3000, 30, 000]</td><td>δ<sub>2</sub></td><td>5.00E-04</td><td>ζ<sub>3</sub></td><td>0.01</td></tr><tr><td>G<sup>max</sup> (#)</td><td>[400, 20, 000]</td><td>δ<sub>3</sub></td><td>2.00E-04</td><td>ζ<sub>4</sub></td><td>0.01</td></tr><tr><td>I<sup>0</sup> (#)</td><td>[10%, 20%] × I<sup>max</sup></td><td>δ<sub>4</sub></td><td>6.00E-04</td><td></td><td></td></tr><tr><td>μ (%)</td><td>[10,30]</td><td>θ</td><td>0.2</td><td></td><td></td></tr></tbody></table>

## References

[1] Ciotti M, Ciccozzi M, Terrinoni A, Jiang W-C, Wang C-B, Bernardini S. The COVID-19 pandemic. Crit Rev Clin Lab Sci 2020;57(6):365-88.

[2] Gereffi G. What does the COVID-19 pandemic teach us about global value chains? The case of medical supplies. J Int Bus Policy 2020;3(3):287-301.

[3] Sokat KY, Altay N. Serving vulnerable populations under the threat of epidemics and pandemics. J Humanit Logist Supply Chain Manag 2021;11(2):176-97.

[4] Al-Dmour H, Salman A, Abuhashesh M, Al-Dmour R, et al. Influence of social media platforms on public health protection against the COVID-19 pandemic via the mediating effects of public health awareness and behavioral changes: integrated model. J Med Internet Res 2020;22(8):e19996.

[5] Gonçalves P, Ferrari P, Crivelli L, Albanese E. Informing hospital capacity and social distancing policies during the covid-19 pandemic. Available at SSRN 38927652021;.

[6] Brusset X, Jebali A, La Torre D. Production optimisation in a pandemic context. Int J Prod Res 2023;61(5):1642-63.

[7] Ekici A, Keskinocak P, Swann JL. Modeling influenza pandemic and planning food distribution. Manuf Serv Oper Manag 2014;16(1):11-27.

[8] Long EF, Nohdurfurt E, Spinler S. Spatial resource allocation for emerging epidemics: a comparison of greedy, myopic, and dynamic policies. Manuf Serv Oper Manag 2018;20(2):181-98.

[9] Ivanov D, Dolgui A. The shortage economy and its implications for supply chain and operations management. Int J Prod Res 2022;60(24):7141-54.

[10] Ivanov D, Keskin BB. Post-pandemic adaptation and development of supply chain viability theory. Omega 2023;116:102806.

[11] Ivanov D, Dolgui A., Blackhurst J.V., Choi T.-M.. Toward supply chain viability theory: from lessons learned through COVID-19 pandemic to viable ecosystems. 2023.

[12] Li S, He Y, Huang H, Lin J, Ivanov D. Supply chain hoarding and contingent sourcing strategies in anticipation of price hikes and product shortages. IISE Trans 2023;0(0):1-27. just-accepted

[13] Queiroz MM, Ivanov D, Dolgui A, Wamba SF. Impacts of epidemic outbreaks on supply chains: mapping a research agenda amid the COVID-19 pandemic through a structured literature review. Ann Oper Res 2020;319(1):1159-96.

[14] Ivanov D, Dolgui A. Viability of intertwined supply networks: extending the supply chain resilience angles towards survivability, a position paper motivated by COVID-19 outbreak. Int J Prod Res 2020;58(10):2904-15.

[15] Choi T-M, Sethi SP. Innovative service operations for survivals of SMEs under COVID-19: two cases in hong kong. IEEE Eng Manag Rev 2021;49(1):50-4.

[16] Vahdani B, Mohammadi M, Thevenin S, Gendreau M, Dolgui A, Meyer P. Fair-split distribution of multi-dose vaccines with prioritized age groups and dynamic demand: the case study of COVID-19. Eur J Oper Res 2023;310(3):1249-72.

[17] Kai D., Goldstein G.-P., Morgunov A., Nangalia V., Rotkirch A.. Universal masking is urgent in the COVID-19 pandemic: SEIR and agent based models, empirical validation, policy recommendations. arXiv preprint arXiv: 2004135532020;1-19.

[18] Ma Y, Pei S, Shaman J, Dubrow R, Chen K. Role of meteorological factors in the transmission of SARS-CoV-2 in the United States. Nat Commun 2021;12(1):1-9.

[19] Goenka A., Liu L., Nguyen M.-H.. Modeling optimal quarantines under infectious disease related mortality2020;

[20] Goenka A, Liu L, Nguyen M-H. SIR economic epidemiological models with disease induced mortality. J Math Econ 2021;93:102476.

[21] Arfan M, Shah K, Abdeljawad T, Mlaiki N, Ullah A. A Caputo power law model predicting the spread of the COVID-19 outbreak in pakistan. Alex. Eng J 2021;60(1):447-56.

[22] Brusset X, Davari M, Kinra A, Torre DL. Modelling COVID-19 ripple effect and global supply chain productivity impacts using a reaction-diffusion time-space sis model. In: Advances in production management systems. Artificial intelligence for sustainable and resilient production systems: IFIP WG 5.7 international conference, APMS 2021, Nantes, France, September 5-9, 2021, proceedings, Part IV. Springer; 2021. p. 3-12.

[23] La Torre D, Liuzzi D, Marsiglio S. Geographical heterogeneities and externalities in an epidemiological-macroeconomic framework. J Public Econ Theory 2022;24(5):1154-81.

[24] La Torre D, Liuzzi D, Maggistro R, Marsiglio S. Mobility choices and strategic interactions in a two-group macroeconomic-epidemiological model. Dyn Games Appl 2022;12:110-32.

[25] Qi H, Xiao S, Shi R, Ward MP, Chen Y, Tu W, et al. COVID-19 transmission in Mainland China is associated with temperature and humidity: a time-series analysis. Sci Total Environ 2020;728:138778.

[26] Betcheva L., Erhun F., Feylessoufi A., Fryers P., Gonçalves P., Jiang H., Kattuman P., Pape T., Pari A., Scholtes S., et al. An adaptive research approach to COVID-19 forecasting for regional health systems in england. Available at SSRN 36952582021;.

[27] Kushwaha S, Bahl S, Bagha AK, Kulwinder Singh P, Javaid M, Haleem A, et al. Significant applications of machine learning for COVID-19 pandemic. J Ind Integr Manag 2020;5(4):453-79.

[28] Yang Z, Zeng Z, Wang K, Wong S-S, Liang W, Zanin M, et al. Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions. J Thorac Dis 2020;12(3):165.

[29] Gnanvi J, Salako KV, Kotanmi B, Kakaï RG. On the reliability of predictions on COVID-19dynamics: a systematic and critical review of modelling techniques. Infect Dis Model 2021;6:258-72.

[30] Nikolopoulos K, Punia S, Schäfers A, Tsinopoulos C, Vasilakis C. Forecasting and planning during a pandemic: COVID-19 growth rates, supply chain disruptions, and governmental decisions. Eur J Oper Res 2021;290(1):99-115.

[31] Kassa SM, Ouhinou A. The impact of self-protective measures in the optimal interventions for controlling infectious diseases of human population. J Math Biol 2015;70(1):213-36.

[32] Castilho C. Optimal control of an epidemic through educational campaigns. Electron J Differ Equ 2006;2006(125):1-11.

[33] Boucekkine R, Hritonenko N, Yatsenko Y. On the optimal control of the vintage capital growth model with endogenous labour supply. Math Econ Lett 2013;1(1):3-7.

[34] Boucekkine R., Camacho C., Fabbri G.. On the optimal control of some parabolic partial differential equations arising in economics2013. halshs-00831042;.

[35] Dolgui A, Ivanov D, Sethi SP, Sokolov B. Scheduling in production, supply chain and industry 4.0 systems by optimal control: fundamentals, state-of-the-art and applications. Int J Prod Res 2019;57(2):411-32.

[36] La Torre D, Malik T, Marsiglio S. Optimal control of prevention and treatment in a basic macroeconomic-epidemiological model. Math Soc Sci 2020;108:100-8.

[37] El Ouardighi F, Khmelnitsky E, Sethi SP. Epidemic control with endogenous treatment capability under popular discontent and social fatigue. Prod Oper Manag 2022;31(4):1734-52.

[38] Gaff H, Schaefer E. Optimal control applied to vaccination and treatment strategies for various epidemiological models. Math Biosci Eng 2009;6(3):469.

[39] Karsten F, Slikker M, Van Houtum G.-J. Resource pooling and cost allocation among independent service providers. Oper Res 2015;63(2):476-88.

[40] Gurvich I, Van Mieghem JA. Collaboration and multitasking in networks: prioritization and achievable capacity. Manag Sci 2018;64(5):2390-406.
[41] Carew S, Nagarajan M, Shechter S, Arneja J, Skarsgard E. Dynamic capacity allocation for elective surgeries: reducing urgency-weighted wait times. Manuf Serv Oper Manag 2021;23(2):407-24.

[42] Struben J, Lee BH, Bingham CB. Collective action problems and resource allocation during market formation. Strategy Sci 2020;5(3):245-70.

[43] Yildirim G, Tatar Y. Simplified agent-based resource sharing approach for WS-N-WSN interaction in IoT/CPS projects. IEEE Access 2018;6:78077-91.

[44] Yin X., Buyuktaktakin I.E., Patel B.P.. COVID-19: optimal allocation of ventilator supply under uncertainty and risk. Available at SSRN 38011832021;

[45] Zhao H, Deshpande V, Ryan JK. Inventory sharing and rationing in decentralized dealer networks. Manag Sci 2005;51(4):531-47.

[46] Pachkova EV. Restricted reallocation of resources. Eur J Oper Res 2009;196(3):1049-57.

[47] Rahimi I, Chen F, Gandomi AH. A review on COVID-19 forecasting models. Neural Comput Appl 2021;1-11. doi:10.1007/s00521-020-05626-8.

[48] Kaplan EH. Om forum-COVID-19 scratch models to support local decisions. Manuf Serv Oper Manag 2020;22(4):645-55.

[49] Anand N, Sabarinath A, Geetha S, Somanath S. Predicting the spread of COVID-19 using sir sir model augmented to incorporate quarantine and testing. Trans Indian Natl Acad Eng 2020;5(2):141-8.

[50] Cooper I, Mondal A, Antonopoulos CG. A SIR model assumption for the spread of COVID-19 in different communities. Chaos, Solitons Fractals 2020;139:110057.

[51] Chen Y-C, Lu P-E, Chang C-S, Liu T-H. A time-dependent sir model for COVID-19 with undetectable infected persons. IEEE Trans Netw Sci Eng 2020;7(4):3279-94.

[52] Abou-Ismail A. Compartmental models of the COVID-19 pandemic for physicians and physician-scientists. SN Compr Clin Med 2020;2:852-8.

[53] Wangping J, Ke H, Yang S, Wenzhe C, Shengshu W, Shanshan Y, et al. Extended SIR prediction of the epidemics trend of COVID-19 in Italy and compared with Hunan, China. Front Med 2020;7:169.

[54] Ifguis O, El Ghozlani M, Ammou F, Moutcine A, Abdellah Z. Simulation of the final size of the evolution curve of coronavirus disease in morocco using the SIR model. J Environ Public Health 2020;2020.

[55] Calafiore GC, Novara C, Possieri C. A modified SIR model for the COVID-19 connection in Italy. In: 2020 59th IEEE conference on decision and control (CDC). IEEE; 2020. p. 3889-94.

[56] He S, Peng Y, Sun K. SEIR modeling of the COVID-19 and its dynamics. Nonlinear Dyn 2020;101(3):1667-80.

[57] Liao Z, Lan P, Liao Z, Zhang Y, Liu S. TW-SIR: time-window based SIR for COVID-19forecasts. Sci Rep 2020;10(1):1-15.

[58] Alanazi SA, Kamruzzaman M, Alruwaili M, Alshammari N, Alqahtani SA, Karime A. Measuring and preventing COVID-19 using the SIR model and machine learning in smart health care. J Healthc Eng 2020;2020:1-12.

[59] Bagal DK, Rath A, Barua A, Patnaik D. Estimating the parameters of susceptible-infected-recovered model of COVID-19 cases in India during lockdown periods. Chaos, Solitons Fractals 2020;140:110154.

[60] Chen X, Li J, Xiao C, Yang P. Numerical solution and parameter estimation for uncertain SIR model with application to COVID-19. Fuzzy Optim Decis Mak 2021;20(2):189-208.

[61] Kyurkchiev N, Kyurkchiev V, Iliev A, Rahnev A. A new modifications of the SIR/SEIR models with "intervention polynomial factor", methodological aspects. Int J Differ Equ Appl 2021;20(1):15-30.

[62] Zisad SN, Hossain MS, Hossain MS, Andersson K. An integrated neural network and SEIR model to predict COVID-19. Algorithms 2021;14(3):94.

[63] Alenezi MN, Al-Anzi FS, Alabdulrazzaq H. Building a sensible sir estimation model for COVID-19outspread in kuwait. Alex Eng J 2021;60(3):3161-75.

[64] Hassen HB, Elaud A, Salah NB, Masmoudi A. A sir-poisson model for COVID-19: evolution and transmission inference in the maghreb central regions. Arabian J Sci Eng 2021;46(1):93-102.

[65] Prodanov D. Analytical parameter estimation of the sir epidemic model. applications to the COVID-19 pandemic. Entropy 2021;23(1):59.

[66] Efimov D, Ushirobira R. On an interval prediction of COVID-19 development based on a SEIR epidemic model. Annu Rev Control 2021;51:477-87.

[67] Singh P, Gupta A. Generalized sir (GSIR) epidemic model: an improved framework for the predictive monitoring of COVID-19 pandemic. ISA Trans 2021. doi:10.1016/j.isatra.2021.02.016. in press

[68] López L, Rodo X. A modified SEIR model to predict the COVID-19outbreak in Spain and Italy: simulating control scenarios and multi-scale epidemics. Results Phys 2021;21:103746.

[69] Das M, Samanta G. Optimal control of fractional order COVID-19 epidemic spreading in Japan and India 2020. Biophys Rev Lett 2020;15(04):207-36.

[70] Saha S, Samanta G, Nieto JJ. Epidemic model of COVID-19 outbreak by inducing behavioural response in population. Nonlinear Dyn 2020;102(1):455-87.

[71] Saha S, Samanta G. Modelling the role of optimal social distancing on disease prevalence of COVID-19 epidemic. Int J Dyn Control 2021;9(3):1053-77.

[72] Mehrotra S, Rahimian H, Barah M, Luo F, Schantz K. A model of supply-chain decisions for resource sharing with an application to ventilator allocation to combat COVID-19. Naval Res Logist (NRL) 2020;67(5):303-20.

[73] Parker F., Sawczuk H., Ganjkhanloo F., Ahmadi F., Ghobadi K.. Optimal resource and demand redistribution for healthcare systems under stress from COVID-19. arXiv preprint arXiv:201103528 2020.

[74] Blanco V., Gázquez R., Leal M.. Reallocating and sharing health equipments in sanitary emergency situations: the COVID-19 case in Spain. arXiv preprint arXiv:2012020622020.

[75] Li X, Ghadami A, Drake JM, Rohani P, Epureanu BI. Mathematical model of the feedback between global supply chain disruption and COVID-19 dynamics. Sci Rep 2021;11(1):15450.

[76] Rozhkov M, Ivanov D, Blackhurst J, Nair A. Adapting supply chain operations in anticipation of and during the COVID-19 pandemic. Omega 2022;110:102635.

[77] Mohammadi M, Dehghan M, Pirayesh A, Dolgui A. Bi-objective optimization of a stochastic resilient vaccine distribution network in the context of the COVID-19 pandemic. Omega 2022;113:102725.

[78] Ala'raj M, Majdalawieh M, Nizamuddin N. Modeling and forecasting of COVID-19 using a hybrid dynamic model based on SEIRD with ARIMA corrections. Infect Dis Model 2021;6:98-111.

[79] Ko J, Nazarian E, Nam Y, Guo Y. Integrated redistricting, location-allocation and service sharing with intra-district service transfer to reduce demand overload and its disparity. Comput Environ Urban Syst 2015;54:132-43.

[80] Saedinia R, Vahdani B, Etebari F, Nadjafi BA. Robust gasoline closed loop supply chain design with redistricting, service sharing and intra-district service transfer. Transp Res Part E 2019;123:121-41.

[81] Neumann-Saavedra BA, Crainic TG, Gendron B, Mattfeld DC, Römer M. Integrating resource management in service network design for bike-sharing systems. Transp Sci 2020;54(5):1251-71.

[82] Vahdani B, Veysmoradi D, Mousavi S, Amiri M. Planning for relief distribution, victim evacuation, redistricting and service sharing under uncertainty. Socio-Economic Plan Sci 2021;101158.

[83] Buonomo B, d'Onofrio A, Lacitignola D. Globally stable endemicity for infectious diseases with information-related changes in contact patterns. Appl Math Lett 2012;25(7):1056-60.

[84] Behncke H. Optimal control of deterministic epidemics. Optim Control Appl Methods 2000;21(6):269-85.

[85] Joshi HR, Lenhart S, Li MY, Wang L. Optimal control methods applied to disease models. Contemp Math 2006;410:187-208.

[86] Benders JF. Partitioning procedures for solving mixed-variables programming problems. Numer Math 1962;4(1):238-52.

[87] Santoso T, Ahmed S, Goetschalckx M, Shapiro A. A stochastic programming approach for supply chain network design under uncertainty. Eur J Oper Res 2005;167(1):96-115.

[88] Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci 2002;180(1-2):29-48.

[89] La Salle JP. The stability of dynamical systems. SIAM; 1976.

[90] Coddington EA, Levinson N. Theory of ordinary differential equations. Tata McGraw-Hill Education; 1955.

[91] Fleming W, Rishel R, Marchuk G, Balakrishnan A, Borovkov A, Makarov V, et al. Applications of mathematics. Determ Stoch Optim Control 1975;1:1-212.

[92] Pontryagin LS. Mathematical theory of optimal processes. CRC Press; 1987.

[93] Kirk DE. Optimal control theory: an introduction. Courier Corporation; 2004.

[94] Salje H, Kiem CT, Lefrancq N, Courtejoie N, Bosetti P, Paireau J, Andronico A, Hozé N, Richet J, Dubost C-L, et al. Estimating the burden of SARS-CoV-2 in France. Science 2020;369(6500):208-11.

[95] Kiem T, Lefrancq N, Courtejoie N, Bosetti P, Paireau J, Andronico A, et al. Erratum: estimating the burden of SARS-CoV-2 in France. Science 2020;368(6498).