# Fair-split distribution of multi-dose vaccines with prioritized age groups and dynamic demand: The case study of COVID-19

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

► To cite this version:

Behnam Vahdani, Mehrdad Mohammadi, Simon Thevenin, Michel Gendreau, Alexandre Dolgui, et al.. Fair-split distribution of multi-dose vaccines with prioritized age groups and dynamic demand: The case study of COVID-19. European Journal of Operational Research, 2023, 310 (3), pp.1249-1272. 10.1016/j.ejor.2023.03.032. hal-04232211

HAL Id: hal-04232211

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

Submitted on 9 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Innovative Applications of O.R.

# Fair-split distribution of multi-dose vaccines with prioritized age groups and dynamic demand: The case study of COVID-19

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

Behnam Vahdaniª, Mehrdad Mohammadiª,b,*, Simon Theveninᶜ, Michel Gendreauᵈ, Alexandre Dolguijc, Patrick Meyerª

## ARTICLE INFO

**Article history:**

Received 4 April 2022

Accepted 25 March 2023

Available online 20 April 2023

**Keywords:**

Distribution

Vaccine distribution

Fairness

Dynamic demand

Benders decomposition

## ABSTRACT

The emergence of the SARS-CoV-2 virus and new viral variations with higher transmission and mortality rates have highlighted the urgency to accelerate vaccination to mitigate the morbidity and mortality of the COVID-19 pandemic. For this purpose, this paper formulates a new multi-vaccine, multi-depot location-inventory-routing problem for vaccine distribution. The proposed model addresses a wide variety of vaccination concerns: prioritizing age groups, fair distribution, multi-dose injection, dynamic demand, etc. To solve large-size instances of the model, we employ a Benders decomposition algorithm with a number of acceleration techniques. To monitor the dynamic demand of vaccines, we propose a new adjusted susceptible-infectious-recovered (SIR) epidemiological model, where infected individuals are tested and quarantined. The solution to the optimal control problem dynamically allocates the vaccine demand to reach the endemic equilibrium point. Finally, to illustrate the applicability and performance of the proposed model and solution approach, the paper reports extensive numerical experiments on a real case study of the vaccination campaign in France. The computational results show that the proposed Benders decomposition algorithm is 12 times faster, and its solutions are, on average, 16% better in terms of quality than the Gurobi solver under a limited CPU time. In terms of vaccination strategies, our results suggest that delaying the recommended time interval between doses of injection by a factor of 1.5 reduces the unmet demand up to 50%. Furthermore, we observed that the mortality is a convex function of fairness and an appropriate level of fairness should be adapted through the vaccination.

© 2023 The Authors. Published by Elsevier B.V.

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

## 1. Introduction

As of 20 March 2022, the ongoing SARS-CoV-2 virus (COVID-19, hereafter) pandemic has engendered approximately 427 million confirmed cases and 6.09 million fatalities worldwide (Johns Hopkins University & Medicine, 2021). Since the beginning of the pandemic, almost all countries around the world have adopted a variety of measures to restrict the propagation of COVID-19. These measures include individual precautions ranging from wearing face masks to social distancing, non-pharmaceutical measures

from governments such as quarantine and border controls, to name a few, and pharmaceutical propositions resorting to vaccination as the most effective way (Kaplan, 2020; Pagliusi et al., 2020). In this regard, the accelerated advancement and outset of some innovative COVID-19 vaccines is a remarkable and exceptional achievement of contemporary science, where Pfizer-BioNTech was the first vaccine to be conferred a temporary authorization for emergency uses in the United Kingdom (Shaw et al., 2021).

Although the rapid development of such effective vaccines against COVID-19 in a short time is among the greatest scientific achievements (Graham, 2020; Lurie et al., 2020), successfully vaccinating the global population encounters many challenges, from production to distribution, deployment, and even social acceptance. More precisely, two significant concerns need to be handled to guarantee the success of a vaccination plan in pandemic situations.

Firstly, coming up with an accurate estimation of the demand for vaccines, regardless of their types, which can immune people

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

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

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

<footnote>ᵈCIRRELT and Département de Mathématiques et Génie Industriel, Polytechnique Montréal, P.O. Box 6079, Station Centre-Ville, Montréal H3C 3A7, Canada</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), michel.gendreau@polymtl.ca (M. Gendreau), alexandre.dolgui@imt-atlantique.fr (A. Dolgui), patrick.meyer@imt-atlantique.fr (P. Meyer).</footnote>
against COVID-19 to confine its propagation. This challenge be-
comes even more critical knowing that the number of companies
producing COVID-19 vaccines is minimal, and almost all countries
have to meet their needs by purchasing vaccines from them. How-
ever, countries that can produce vaccines prioritize meeting their
own needs first (Alam et al., 2021). This issue forces other coun-
tries to be more accurate in estimating their demand when pur-
chasing vaccines. Nonetheless, the precision of such estimation in
the COVID-19 pandemic could be jeopardized due to many fac-
tors, including the emergence of new variants of the virus, non-
compliance with health protocols by some individuals, monopo-
lies in the distribution and consumption of vaccines by developed
countries, and insufficient infrastructure and distribution capacity
of vaccines (Alam et al., 2021; Burgos et al., 2021; Tyagi et al.,
2021; Wouters et al., 2021). Such factors lead to important fluctu-
ations and dynamicity in vaccine demand (Besiou & Van Wassen-
hove, 2021). Hence, providing a precise estimation of the progres-
sion of the pandemic allows more fruitful management of the lo-
gistics challenges. In this regard, in order to provide an accurate
estimation of disease propagation, a broad range of efficient ap-
proaches have been employed in the literature, including artifi-
cial intelligence approaches such as time series (Betcheva et al.,
2021; Qi et al., 2020), machine learning (Kushwaha et al., 2020),
and deep learning (Yang et al., 2020), agent-based approaches (Kai
et al., 2020), metrological and meta-population models (Ma et al.,
2021), and compartmental epidemiological models (Arfan et al.,
2021). From the technical point of view, apart from compartmental
epidemiological models, the other approaches hang on data acces-
sibility. Hence, these approaches almost become inapplicable since
data is unavailable and unreliable in the early phases of a pan-
demic (Gnanvi et al., 2021). Also, historical data cannot be ex-
poited due to the occurrence of different virus variants. Moreover,
these approaches cannot consider various government strategies,
such as quarantine and raising people's awareness (Biswas & Alfan-
dari, 2022; Nikolopoulos et al., 2021). Finally, they cannot optimize
the measures related to implementing the policies established by
the governments, such as reducing the number of deaths or min-
imizing system costs to control the propagation of a pandemic. In
this regard, the synergy between compartmental epidemiological
models and the optimal control problem (OCP) as a powerful dy-
namic optimization approach not only tackles the barriers men-
tioned above (Castilho, 2006) but also offers an endemic equilib-
rium point of the COVID-19 contagious influence (Gaff & Schae-
fer, 2009). However, due to the complexity of implementing these
methods, researchers have rarely used them.

Secondly, it is naive to believe that most of the vaccination
process has been completed by purchasing the estimated amount
of vaccine. In fact, distributing vaccines, especially COVID-19 vac-
cines with miscellaneous concerns such as minimal supplies, pri-
oritizing age groups, fair distribution, multi-dose injection, split
delivery of vaccines, vaccine loss during transit, and vial-opening,
is the second challenge that makes it a highly daunting process
(Abbasi et al., 2020; Chen et al., 2020; De Boeck et al., 2020;
Ivanov, 2020). Hence, this paper aims to develop an efficient vac-
cine distribution network to deliver a limited supply of vaccines
between different vaccination centers. Developing such a network
requires a broad range of decisions to be made (Duijzer et al.,
2018c; Manupati et al., 2021): 1) Location-allocation, where to lo-
cate vaccine distribution centers and allocate vaccination centers to
them (Ahmadzadeh & Vahdani, 2017; Mohammadi et al., 2014; Ni
et al., 2021; Niakan et al., 2015; Rahimi-Vahed et al., 2015; Vahdani
& Ahmadzadeh, 2019; Vahdani et al., 2017), 2) Inventory, where
(at which center) to hold inventory considering the perishability
of vaccines (Tang et al., 2022; Vahdani & Ahmadzadeh, 2021; Vah-
dani et al., 2011), and 3) Delivery-Routing, what amount of vaccines
to be delivered/routed between different centers using which ve-

hicles (Alam et al., 2021; Capelle et al., 2019; Gmira et al., 2021; Nabavi et al., 2022; Rahimi-Vahed et al., 2015; Rostami et al., 2021; Vidal et al., 2012; 2021).

Speaking of which, a careful examination of the COVID-19 pan-
demic reveals that additional factors should also be addressed in
new pandemics when planning vaccine distribution. This paper ac-
counts for four additional factors to design an efficient vaccination
network. These factors are classified into prioritizing the popula-
tion for vaccination, fair distribution of vaccine, multi-dose vaccina-
tion planning, and split delivery.

First, since the supply of COVID-19 vaccines has initially been
restricted, the world realized that to pilot a vaccination campaign
efficiently, governments all over the world should make a tough
decision on who to vaccinate first. In fact, health policy-makers
all around the world came to a global consensus on vaccination
strategies to prioritize the vaccination of front-line workers (e.g.,
healthcare workers) against the disease, then vulnerable popula-
tions such as older people and people with weak immune systems.
Accordingly, such prioritization decisions should be taken into ac-
count when designing a vaccine distribution network (Mohammadi
et al., 2022).

Second, the lack of COVID-19 vaccines raises the issue of a fair
distribution of the produced vaccines among or within countries.
In fact, vaccine distribution presents a slew of complicated and
contentious matters, including public health, finance, public per-
suasion, and diplomacy (DeRoo et al., 2020; Forni & Mantovani,
2021). Various global organizations, national administrators, and
vaccine manufacturers realized that ethics is vital in making de-
cisions. However, inadequate progress has been achieved in defin-
ing the fair distribution of vaccines globally and nationally. An "eq-
uitable vaccine distribution" is often praised without specifying a
structure or offering any suggestions. This paper investigates the
fairness in the distribution of vaccines among different vaccination
centers, each being responsible to vaccinate the population of their
corresponding region.

Third, another factor affecting the operation of vaccine distribu-
tion centers is multi-dose vials of vaccines recommended/imposed
by the manufacturers to achieve a sufficient level of immunity
against COVID-19. Such recommendation/necessity for multi-dose
vaccination encounters the distribution center with important
challenges. As a matter of fact, multi-dose vaccination imposes a
commitment to the distribution network to imperatively fulfill the
demand for the second dose, if the first dose has previously been
fulfilled. Another point to note is that individuals must receive the
second dose of a vaccine after a given time interval, depending on
the type of their first-dose vaccine (Silva-Cayetano et al., 2021). In
this regard, a number of recent studies have articulated that inject-
ing a second dose from another type of vaccine is possible if it is
compatible with the first dose, although the majority of evidence
has recommended that it is better if two doses have the same kind
(Sampath et al., 2021).

Fourth, the delivery of multi-dose vaccinations puts far more
pressure on the distribution network. This becomes even more
problematic when we know that COVID-19 vaccines must be trans-
ported in refrigerators or freezers, limiting the number of accessi-
ble vehicles. One of the most efficient strategies to cope with this
obstacle is to use the split delivery method (Haddad et al., 2018;
Mohammadi et al., 2020; Veysmoradi et al., 2018).

According to the above-mentioned description, we address the
two problems as a comprehensive two-stage framework. The first
stage concentrates on estimating the dynamic demand for vaccines
and assesses the number of vaccines that need to be provided to
control the pandemic. The second stage focuses on offering a ded-
icated multi-vaccine, multi-depot location-inventory-routing prob-
lem (MVMDLIRP) to distribute vaccines among vaccination centers,
wherein a broad range of characterful concerns of the COVID-19
vaccine distribution is reflected. The first stage offers an adjusted SIR model with the optimal control problem for controlling the dynamicity of demand for vaccines, and the resulting model provides the number of vaccines needed during the vaccination campaign. The second stage formulates a dedicated MVMDLIRP, where a broad range of decisions are considered, including supply, location-allocation, shipment, inventory, and routing concerning various concerns, including perishability, split delivery, capacity, population priorities, and multi-dose delivery. To the best of our knowledge, the MVMDLIRP has yet to be investigated, let alone the other reflected considerations in the current study that make it more realistic and sophisticated. Hence, the contributions of this study regarding the investigated literature in the next section are as follows:

* Introducing a comprehensive and exclusive optimization framework to distribute vaccines to control the pandemic situation,

* Proposing an adjusted SIR model with the optimal control problem to overcome the demand dynamicity and determine the number of vaccines that need to be supplied to control the pandemic, and

* Formulating the MVMDLIRP regarding a broad range of applicable concerns, including prioritizing age groups, fair distribution, multi-dose injection, split delivery of vaccines, vaccine fraction during transit and vial-opening, vaccine deterioration, and dynamic demand.

The remainder of the paper is organized as follows. Section 2 scrutinizes related literature. Section 3.1 provides the description of the SIR epidemiological model with the optimal control problem. Section 4 presents the MVMDLIRP formulation. Section 5 offers the details of the accelerated Benders decomposition algorithm. Section 6 describes the case study and provides numerical results. Next, Section 7 provides a set of sensitivity analyses and managerial insights. Finally, this study is concluded, and future research directions are proposed in Section 8.

## 2. Literature review

The current research focuses on two main phases: estimating the demand for vaccines through disease progression modeling via epidemiological models and the logistics planning of vaccine distribution. Hence, this section scrutinized the related literature of the two mentioned subjects to illustrate research gaps and our contributions.

### 2.1. Phase I—Epidemiological models at the service of vaccine distribution

Asano et al. (2008) suggested a modified Susceptible-Infected-Recovered (SIR) model with the optimal control problem for managing vaccine bait distribution to control the spread of rabies in raccoons, where the rate of vaccine bait distribution was considered as a control variable and determined the number of vaccines needed to control the disease. Nguyen & Carlson (2016) extended the conventional SIR model as a semi-Markovian process to reflect the real-time circumstances of an outbreak for controlling vaccine allocation. Büyüktahtakın et al. (2018) proposed an epidemics-logistics optimization model to control the Ebola virus. This model's decisions include allocating resources regarding capacity and budget restrictions. The objective function minimized the total number of infected people and deaths of infected individuals, where a modified SIR model was proposed to model the Ebola disease transmission. Duijzer et al. (2018b) considered the conventional SIR model to investigate two vaccination strategies, including early and later vaccines. They suggested that a hybrid strategy can diminish the number of infected individuals by more than 50% compared to the best absolute strategy. Duijzer et al.

(2018a) employed the conventional SIR model to introduce an analytical model for allocating the necessitated number of vaccines in order to control an epidemic, where the objective function maximized the health benefit, which was defined as the total number of people who escaped from infection with or without vaccination. Enayati & Özaltın (2020) proposed a modified susceptible-exposed-infected-recovered (SEIR) model to manage influenza vaccine distribution by minimizing the number of vaccines used to immunize people and control such an outbreak, where the Gini coefficient reflected the equity concern. Westerink-Duijzer et al. (2020) applied the conventional SIR model to formulate an influenza disease transmission as part of an analytical model to explore sharing vaccines as a redistribution problem, where health agencies might cooperate based on cooperative game theory. Gamchi et al. (2021) proposed a bi-objective mathematical model for influenza vaccine distribution, where the conventional SIR model with the optimal control problem, adapted from Alcaraz & Vargas-De-León (2012), was employed to prioritize regions and individuals for vaccination. Then a classical vehicle routing problem was used to distribute vaccines. The first objective function minimized the social costs of infected individuals before and after vaccination, and the second objective function minimized the fixed costs of vehicles used in the distribution process.

### 2.2. Phase II—Vaccine distribution problem

Chen et al. (2014) proposed a mathematical model to formulate an inventory-transshipment problem for vaccine distribution in developing countries. The objective function maximizes the number of immunized children and additional doses of vaccines. The authors consider two types of vaccines (frozen and refrigerated), vaccine fractions during transit and vial-opening, and deterioration of vaccine inventory have been the main contributions. Lim et al. (2022) proposed a multi-modal location-transshipment model to design a distribution network of vaccines, where the objective function minimizes the total costs of opening and operating facilities and transportation. The main feature is the consideration of a replenishment frequency, where two discrete options, including monthly and quarterly, were reckoned. Abbasi et al. (2020) proposed a transshipment-assignment problem to allocate COVID-19 vaccine to people, where the objective function minimized the infection risk of people based on the susceptibility rating of target groups, transshipment time, and the number of over-supplied vaccines.

They also considered the possibility of transmitting vaccines between vaccination centers, and different target groups to receive vaccines. Larissa et al. (2021) proposed a mathematical model to formulate a classical vehicle routing and scheduling problem, where the objective function minimized the transit time and vehicle and road penalties, which represented unreliable conditions that can jeopardize vaccine distribution. However, none of the specific considerations of vaccine distribution were considered in the suggested model. Yang et al. (2021) presented a multi-period location-routing model to distribute vaccines.

The considered problem accounts for vehicle and facility capacities, and the objective function minimized the costs of opening facilities and transportation. The authors applied the replenishment frequency suggested by Lim et al. (2022). Rastegar et al. (2021) proposed a multi-period inventory-location model to distribute the influenza vaccine during the COVID-19 pandemic. The model accounts for fairness concerns, where the objective function maximized the minimum delivery-to-demand ratio as an equity measurement. They considered multiple target groups for vaccination without priorities, but a minimum coverage rate for each group was considered as another equitable feature. Tavana et al. (2021) made a similar study to Rastegar et al. (2021), with the
<figure><img src="image_4.png" /><figcaption>Fig. 1. The SIR model.</figcaption></figure>

exception that they formulated the COVID-19 vaccine distribution. The main feature was the consideration of various vaccines, which required very cold and ultra-cold refrigeration for equipping distribution centers and their storage.

Chen et al. (2021) proposed an agent-based simulation model to allocate COVID-19 vaccines to people with the help of social contact network. They concluded that the proposed simulation model was much more effective than other simulations whose dialectics were based on the number of infected, hospitalized, or dead people. Tang et al. (2022) proposed a bi-objective location-inventory model for a vaccination planning problem, where the first objective minimized the total operational costs and the second one minimized the travel distance of people to vaccination centers. Also, they offered a genetic algorithm and dynamic programming method to solve the proposed model. Gilani & Sahebi (2022) considered a sustainable production-inventory-transshipment vaccine problem as a vaccine supply chain during the COVID-19 pandemic, where the decisions of manufacturers' capacity and capacity expansion, the locations of the packing centers of different vaccines and their distribution centers, and transshipment of vaccines were considered. They proposed a multi-objective mathematical model under uncertain supply to formulate the addressed problem. The first objective function minimized the total costs, the second one minimized the total environmental effects, and the third maximized the number of created jobs. Moreover, they offered a robust data-driven approach to tackle the challenges of the uncertain supply of vaccines. Finally, Table 1 summarizes the significant concerns of the investigated literature to demonstrate research gaps and our contributions.

### 3. The estimation of demand for COVID-19 vaccine

#### 3.1. Epidemiological model

In order to formulate the propagation of contagious disease and calculate the number of susceptible, infected, and recovered individuals in a population, the standard SIR model with the possibility of vaccination, as the well-known epidemiological model, is depicted in Fig. 1(a). As can be seen, susceptible individuals become either infected once exposed to the disease or recover if vaccinated, where the population size is fixed and predetermined. Notably, the standard model cannot consider the unique concerns of spreading COVID-19 because one of the crucial factors in model-

ing the spread of this virus is considering the test results on infected people. So, according to the test results, people must be divided into two groups: quarantined and unquarantined (Anand et al., 2020). Consequently, the adjusted version of the SIR model with the possibility of vaccination is offered to model the spread of COVID-19. The configuration of this model is depicted in Fig. 1(b), and its notations are presented in Table 2.

As can be seen, the susceptible individuals (S) become either infected (I) once exposed to the disease via the term $(\beta_1 I + \beta_2 Q + \beta_3 U)S$, or they get recovered (R) if vaccinated via term $uS$. Unlike the classical SIR model, the adjusted SIR model considers that infected individuals become either quarantined (Q) or unquarantined (U) via rates $k\tau$ and $k(1-\tau)$, respectively. In fact, the Q compartment consists of the infected individuals whose test has been positive. These individuals are quarantined/isolated from the rest of the population, and they no longer jeopardize the health of other susceptible individuals. On the other hand, the U compartment represents a group of individuals who are infected but have not still been tested. Contrary to the Q individuals, the U individuals pose the risk of infection to the susceptible individuals. More importantly, these two compartments are correlated with the testing factor $\tau$. Indeed, the higher the level of this factor, the higher the number of individuals who become quarantined. Hence, there are fewer unquarantined people and less risk of infection for susceptible individuals. Finally, the quarantined and unquarantined individuals are differently removed from the system and become recovered (R). It is indeed assumed that the vaccinated individuals will never get infected again.

The justified SIR model of Fig. 1(b) is mathematically formulated as the system of Eq. (1). Indeed, the system of Eq. (1) represents the balance of input and output time-dependent flows at each state and models its fluctuation and dynamicity.

$$ \frac{dS}{dt} = \Lambda - (\beta_1 I + \beta_2 Q + \beta_3 U)S - dS - uS \quad S(0) > 0, $$

$$ \frac{dI}{dt} = (\beta_1 I + \beta_2 Q + \beta_3 U)S - (d + \delta_1 + k)I \quad I(0) \geq 0, $$

$$ \frac{dQ}{dt} = k\tau I - (d + \delta_2)Q - v_1 Q \quad Q(0) \geq 0, $$

$$ \frac{dU}{dt} = k(1 - \tau)I - (d + \delta_3)U - v_2 U \quad U(0) \geq 0, $$

$$ \frac{dR}{dt} = v_1 Q + v_2 U - dR + uS \quad R(0) \geq 0, \quad (1) $$
<table><caption>**Table 1**
Classification of the relevant papers.</caption><thead><tr><th rowspan="3">Ref.</th><th colspan="5">Phase 1</th><th colspan="5">Phase 2</th><th rowspan="3">Solution Method<sup>c</sup></th></tr><tr><th colspan="4">Epidemiological model<sup>a</sup></th><th rowspan="2">Disease</th><th rowspan="2">Cooperation<sup>a</sup></th><th colspan="5">Vaccine supply chain</th></tr><tr><th>wOCP</th><th>wOCP</th><th>Adjusted</th><th>Location</th><th>Allocation</th><th>Transshipment</th><th>Inventory</th><th>Routing</th><th>Specification<sup>b</sup></th><th>Uncertainty</th></tr></thead><tbody><tr><td>Asano et al. (2008)</td><td></td><td></td><td></td><td>Rabies in raccoons</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Nguyen & Carlson (2016)</td><td>✓</td><td>✓</td><td>✓</td><td>General</td><td>SMP</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Büyüktahtakın et al. (2018)</td><td>✓</td><td>✓</td><td>✓</td><td>Ebola</td><td>OR</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Duijzer et al. (2018b)</td><td>✓</td><td></td><td></td><td>Influenza</td><td>OR</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Duijzer et al. (2018a)</td><td>✓</td><td></td><td></td><td>Influenza</td><td>AM</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Enayati & Özaltın (2020)</td><td>✓</td><td></td><td>✓</td><td>Influenza</td><td>OR</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Westerink-Duijzer et al. (2020)</td><td>✓</td><td></td><td></td><td>Influenza</td><td>AM</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Gamchi et al. (2021)</td><td></td><td>✓</td><td></td><td>General</td><td>OR</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Chen et al. (2014)</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>Lim et al. (2022)</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>Abbasi et al. (2020)</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>Larissa et al. (2021)</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>Yang et al. (2021)</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>Rastegar et al. (2021)</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>Tavana et al. (2021)</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>Chen et al. (2021)</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>Tang et al. (2022)</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>Gilani & Sahebi (2022)</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>Thul & Powell (2023)</td><td></td><td></td><td></td><td>COVID-10</td><td>SMP</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Yin et al. (2023)</td><td>✓</td><td>✓</td><td>✓</td><td>COVID-19</td><td>OR</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Current study</td><td></td><td>✓</td><td>✓</td><td>COVID-19</td><td>OR</td><td></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td></td><td></td></tr></tbody></table>

<footnote><sup>a</sup> woOCP: without OCP; wOCP: with OCP; SMP: Semi-Markovian Process; OR: Operations Research; AM: Analytical Model.</footnote>

<footnote><sup>b</sup> A: Multi-period; B: Multiple vaccines; C: Fairness; D: Age groups; E: Prioritized age groups; F: Split delivery; G: Multi-dose; H: Cumulative demand in the second dose; I: Time lag between injection of doses; J: Vaccine fraction; K: Vaccine deterioration; L: Unmet demand; M: Storage capacity; N: Shipping capacity.</footnote>

<footnote><sup>c</sup> MH: Meta-heuritic.</footnote>
<table><caption>**Table 2**
Notations of the adjusted SIR model.</caption><thead><tr><th>Notation</th><th>Description</th></tr></thead><tbody><tr><td>Λ</td><td>Total population at the beginning of the pandemic</td></tr><tr><td>d</td><td>Natural death rate of individuals</td></tr><tr><td>N(t)</td><td>Total population at time t</td></tr><tr><td>S</td><td>Number of susceptible individuals</td></tr><tr><td>I</td><td>Number of infected individuals</td></tr><tr><td>Q</td><td>Number of quarantined individuals</td></tr><tr><td>U</td><td>Number of unquarantined individuals</td></tr><tr><td>R</td><td>Number of recovered individuals</td></tr><tr><td>β<sub>1</sub>/β<sub>2</sub>/β<sub>3</sub></td><td>Disease contagion rates per contact by infected/quarantined/unquarantined individuals</td></tr><tr><td>u</td><td>Proportion of susceptible individuals who are vaccinated</td></tr><tr><td>k</td><td>Information interaction rate through which individuals admit to taking a test</td></tr><tr><td>τ</td><td>The probability of a test to be positive</td></tr><tr><td>ν/ν<sub>1</sub>/ν<sub>2</sub></td><td>Progression rate by which infected/quarantined/unquarantined individuals get recovered</td></tr><tr><td>δ<sub>1</sub>/δ<sub>2</sub>/δ<sub>3</sub></td><td>Particular death rate among infected/quarantined/unquarantined individuals</td></tr></tbody></table>

By solving the isoclines, the basic reproduction ratio $\mathcal{R}_0$ and the endemic equilibrium point can be calculated. In this regard, two equilibrium points are considered as the disease-free equilibrium point (i.e., $E_0(S_0, 0, 0, 0) = (\frac{\Lambda}{d}, 0, 0, 0, 0)$) and the endemic equilibrium point (i.e., $E^*(S^*, I^*, Q^*, U^*, R^*)$, where $\mathcal{R}_0 < 1$).

The basic reproduction ratio $\mathcal{R}_0$ (Wallinga et al., 2010) is defined as the number of new infections caused by a single infectious individual in a completely susceptible population (Duijzer et al., 2018b). Accordingly, in the case of $\mathcal{R}_0 \ge 1$, the pandemic is anticipated to progress. Next, the method proposed by (Van den Driessche & Watmough, 2002) is used to calculate $\mathcal{R}_0$. For this purpose, let $x = (I, Q, U)$ and constitute $\frac{dx}{dt}, \mathcal{S}(x)$, and $v(x)$, where $\mathcal{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} = \mathcal{S}(x) - v(x) \qquad (2) $$

$$ \mathcal{S}(x) = \begin{pmatrix} (\beta_1 I + \beta_2 Q + \beta_3 U)S \\ 0 \\ 0 \end{pmatrix} \qquad (3) $$

$$ v(x) = \begin{pmatrix} (d + \delta_1 + k)I \\ -k\tau I + (d + \delta_2 + \nu_1)Q \\ -k(1 - \tau)I + (d + \delta_3 + \nu_2)U \end{pmatrix} \qquad (4) $$

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

$$ F = (D\mathcal{S}(x))(E_0) = \begin{pmatrix} \beta_1 S_0 & \beta_2 S_0 & \beta_3 S_0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \qquad (5) $$

$$ V = (Dv(x))(E_0) = \begin{pmatrix} (d + \delta_1 + k) & 0 & 0 \\ -k\tau & (d + \delta_2 + \nu_1) & 0 \\ -k(1 - \tau) & 0 & (d + \delta_3 + \nu_2) \end{pmatrix} \qquad (6) $$

Finally, the basic reproduction number is equal to the spectral radius of $FV^{-1}$ (i.e., $\mathcal{R}_0 = \rho(FV^{-1})$, where $\rho(A)$ is the spectral radius of the matrix A). The endemic equilibrium point is then obtained by solving the system of Eq. (7).

$$ \begin{aligned} \Lambda - (\beta_1 I^* + \beta_2 Q^* + \beta_3 U^*)S^* - dS^* - uS^* &= 0, \\ (\beta_1 I^* + \beta_2 Q^* + \beta_3 U^*)S^* - (d + \delta_1 + k)I^* &= 0, \\ k\tau I^* - (d + \delta_2)Q^* - \nu_1 Q^* &= 0, \\ k(1 - \tau)I^* - (d + \delta_3)U^* - \nu_2 U^* &= 0, \\ \nu_1 Q^* + \nu_2 U^* - dR^* + uS^* &= 0, \end{aligned} \qquad (7) $$

## 3.2. The optimal control problem (OCP)

At this stage, we investigate how intervention measures can mitigate the population's disease burden. For this aim, the system of Eq. (1) is revised to reflect the impact of the control variable $u(t)$ over time. In this regard, the aim is to minimize the cost incurred by infected and quarantined individuals and the immunization program in order to maximize the number of recovered individuals using the possible minimal control variable $u(t)$. The optimal value of the control variable $u(t)$ (i.e., $u^*(t)$) can be found through the optimal control problem (8) and (9). Indeed, this value makes the system follow a trajectory state variable that minimizes the performance measure $J$, where positive parameters $w_1, w_2$, and $w_3$ are weight constants balancing the units of the integrand.

$$ J[u(t)] = \int_{0}^{T_f} [w_1 I(t) + w_2 Q(t) + w_3 u^2(t)] dt \qquad (8) $$

s.t.:

$$ \frac{dS}{dt} = \Lambda - (\beta_1 I + \beta_2 Q + \beta_3 U)S - dS - u(t)S. $$

$$ \frac{dI}{dt} = (\beta_1 I + \beta_2 Q + \beta_3 U)S - (d + \delta_1 + k)I, $$

$$ \frac{dQ}{dt} = k\tau I - (d + \delta_2)Q - \nu_1 Q, $$

$$ \frac{dU}{dt} = k(1 - \tau)I - (d + \delta_3)U - \nu_2 U, $$

$$ \frac{dR}{dt} = \nu_1 Q + \nu_2 U - dR + u(t)S \qquad (9) $$

with initial conditions $S(0) > 0, I(0) \ge 0, Q(0) \ge 0, U(0) \ge 0, R(0) \ge 0$. The objective function (8) denotes the total incurred cost, where the first term calculates the cost when infected individuals consult medical professionals concerning symptoms and intend to take a test. The second term also represents the cost imposed on the government for quarantined individuals (e.g., social and home-care aids). The third term computes the cost of the immunization program. Noteworthy, the square of the control variable $u(t)$ reflects the severity of vaccination side effects (Laarbi et al., 2013). So, the region for the control intervention $u(t) \in [0, 1]$ is given as $\cap = \{u(t) | u(t) \in [0, 1], t \in [0, T_f]\}$, where $T_f$ is the final time up to which the control policy is executed, and also $u(t)$ is a measurable and bounded function. The optimal control intervention $u^*(t)$ exists in $\cap$ that minimizes the cost function $J$.

**Theorem 1.** The optimal control intervention $u^*$ in $\cap$ of the optimal control problem (8) and (9) exists such that $J(u^*) = \min[J(u)]$.

**Proof.** A mathematical proof is provided in Appendix A. □
In what follows, in order to obtain the optimal control variable $u^*$, the Hamiltonian function is formulated by introducing adjoint variables $\lambda = (\lambda_1, \lambda_2, \dots, \lambda_5) \in \mathbb{R}^5$, and minimized by applying Pontryagin's Maximum Principle (Pontryagin, 1987) (please see more explanations in Appendix B).

**Theorem 2.** If $u^*$ is the optimal control variable and $S^*, I^*, Q^*, U^*, R^*$ are optimal state variables of the optimal control problem (8) and (9), there exist then adjoint variables $\lambda = (\lambda_1, \lambda_2, \dots, \lambda_5) \in \mathbb{R}^5$ that satisfy the canonical system of Eq. (10):

$$ 
\begin{aligned}
\frac{d\lambda_1}{dt} &= \lambda_1(\beta_1 I + \beta_2 Q + \beta_3 U + d + u) \\
&\quad - \lambda_2(\beta_1 I + \beta_2 Q + \beta_3 U) - \lambda_5(u), \\
\frac{d\lambda_2}{dt} &= -w_1 + \lambda_1(\beta_1 S) - \lambda_2(\beta_1 S - d - \delta_1 - k) \\
&\quad - \lambda_3(k\tau) - \lambda_4(k(1 - \tau)), \\
\frac{d\lambda_3}{dt} &= -w_2 + \lambda_1(\beta_2 S) - \lambda_2(\beta_2 S) + \lambda_3(d + \delta_2 + v_1) - \lambda_5(v_1), \\
\frac{d\lambda_4}{dt} &= \lambda_1(\beta_3 S) - \lambda_2(\beta_3 S) + \lambda_4(d + \delta_3 + v_2) - \lambda_5(v_2), \\
\frac{d\lambda_5}{dt} &= \lambda_5(d),
\end{aligned}
\quad (10) $$

With the transversality conditions $\lambda_i(T_f) = 0, i = 1, \dots, 5$, the optimal control variable $u^*$ is given as:

$$ u^* = \min \left\{ \max \left\{ 0, \frac{(\lambda_1 - \lambda_5)S}{2w_3} \right\}, 1 \right\}. \quad (11) $$

**Proof.** A mathematical proof is provided in Appendix C. □

Finally, we can establish the optimal system by replacing the obtained optimal control variable (11) at the optimal control problem (8) and (9) with optimal state variables. Noteworthy, as far as the pandemic is not over, $uS(t)$ helps estimate the demand for vaccines for the population that should be provided. More importantly, at the endemic equilibrium point (i.e., $\mathcal{R}_0 < 1$), the demand for the vaccines can be estimated as $u*S^*$ resulting from the optimal system. Accordingly, the government should increase the value of $u$ until it reaches its optimal value $u^*$ to enter the endemic equilibrium point.

## 4. Mathematical formulation

This section formulates the MVMDLIRP to distribute purchased vaccines from a central storage hub to a set of vaccination centers through a set of intermediate distribution centers/depots. In this network, the location of the central storage hub and vaccination center is fixed, while the locations of the distribution centers need to be determined. Both distribution and vaccination centers are able to hold an inventory of vaccines in different periods. Moreover, two groups of vaccines are available that differ in their physicochemical properties, including mRNA vaccines (e.g., Pfizer-BioNTech and Moderna) and Adenovirus vector vaccines (e.g., AstraZeneca and Sputnik V). These vaccines require different storage modes while the former should be stored in ultra-cold freezers between [−80 °C, −15 °C], the latter should be stored in refrigerators between [2 °C, 8 °C]. Therefore, different inventory holding costs are applied to each type of vaccine. It is also considered that a fraction of the vaccines might perish/deteriorate from one period to the next and be lost. To make the problem more realistic, the population is divided into different classes based on their age, and each class has a particular priority to be vaccinated.

According to the description, in order to provide an applicable plan for vaccine distribution, three main categories of decisions should be made simultaneously: design, distribution, and inventory.

The design decision focuses on locating distribution centers from a set of potential locations. Notably, these potential places have specialized storage facilities and usually only need to make small changes in the deployment of equipment. Therefore, when the pandemic spread is controlled, we need to close these centers or reduce their numbers. So, this can be performed with minimal energy, time, and cost, and these facilities can return to their previous activities. Therefore, this level of decision should be reviewed periodically according to the volume of demand for vaccination. Therefore, in the presented model, we have periodically reviewed the decisions related to this level, confirming that the location decision in the current problem does not belong to strategic level decisions. Next, decisions related to the distribution of vaccines should be made. In this regard, the vaccine distribution should be performed from opened distribution centers as intermediate nodes between the central storage hub and vaccination centers. In fact, the inventory levels in opened distribution centers and vaccination centers are the decision criteria for supplying vaccines from the central storage hub and sending them to vaccination centers.

Additionally, the delivery of vaccines from opened distribution centers to vaccination centers is performed via a limited number of refrigerated trucks, and each truck has a limited capacity. At each distribution center, the trucks are loaded with appropriate vaccines, and each truck is responsible for delivering a number of vaccines to different vaccination centers. Consequently, the shipping of vaccines from distribution centers to vaccination centers is addressed as a vehicle routing problem. What is more, due to the dynamic nature of demand, resource limitations, and capacities of trucks and vaccination centers, in most cases, it is impossible to meet the vaccine demand of vaccination centers in one visit. Therefore, in order to overcome this problem, split delivery is intended for the distribution of vaccines. Thus, several trucks from different distribution centers can meet the demand for a vaccination center.

Moreover, since the amount of vaccines produced is much lower than the current demand, a fair distribution system is considered between vaccination centers. For the sake of fairness, a service gap level is defined, which guarantees that the difference between the ratio of satisfied demand to the demand of two different vaccination centers cannot exceed the specified service gap. Therefore, we can guarantee that all vaccination centers can supply the vaccine demand of the people assigned to them within an acceptable boundary. Also, a predetermined percentage of each vaccination center's total demand is guaranteed to be met. In this regard, the shortage of vaccines is also considered by unmet demand that happens when the demand for vaccines is more significant than the supply power. Additionally, the cost of unmet demand is correspondent by morbidity and mortality costs that an unvaccinated individual imposes on the health system. Notably, we have considered this cost to be the same for all target age groups because, on the one hand, this issue is of equal importance to all residents in a community, and on the other hand, it expresses a kind of fairness in the distribution of vaccines.

More importantly, COVID-19 vaccines should be injected in two doses, which should be done in a predefined time interval. Although this time has a lower bound according to the assessment of the World Health Organization and vaccine-producing companies, this lag is not necessarily fixed. In fact, the inventory of purchased vaccines, the severity of the disease's spread, and the government's policies in the vaccination process can change the time interval. In addition, as previously discussed, the second dose should be the same as the first dose, so each vaccination center must respond to two types of demand in each period. The first type, which we call new demand, is related to people who want to inject their first dose. The second type of demand is related to people who want to inject their second dose, and the type of vaccine required for
the second injection must be the same as the first dose. Therefore, each period's cumulative demand must be calculated for each distribution center. The main assumptions of the described problem are as follows:

* The locations of the central storage hub and vaccination centers are predefined;

* The capacities of the central storage hub, distribution centers, vaccination centers, and vehicles are limited;

* The shortage of vaccines is allowed, which is considered by unmet demand;

* COVID-19 vaccines must be injected twice with a predefined time interval;

* The second dose of vaccines must be injected the same as the first dose for each individual;

* People are classified into a number of classes based on their ages;

* The satisfaction of demand for vaccines is based on the priority of each class of people;

* Each vaccine center can be visited by vehicles more than once, which is implied by the split delivery;

* The demand for vaccines is uncertain and has a dynamic pattern, which is estimated by the proposed adjusted SIR model with OCP;

* Vaccines can become unusable due to fractions during shipping or opening and deterioration;

* The unit cost of unmet demand is an estimation of the cost that an unvaccinated individual imposes on the health system (e.g., test, treatment, hospitalization, and death costs);

* The government obliges the people to receive the total dose of vaccines in order to receive social services.

## 4.1. Notations

Table 3 lists the notations to formulate the MVMDLIRP.

## 4.2. The MVMDLIRP formulation

This section formulates the MVMDLIRP in terms of a mixed-integer non-linear programming model using the notations described in Section 4.1, where the demand for vaccines is estimated from the proposed modified SIR model (9). The non-linear terms are linearized afterward.

### 4.2.1. Objective function

The objective function of the MVMDLIRP formulation is presented as (12) which minimizes the total costs of the system.

$$ \min Z = \sum_{j \in J} \sum_{t \in T} f_j Z_{jt} + \sum_{p \in P} \sum_{j \in J} \sum_{t \in T} c_{pj} U_{pjt} \\ + \sum_{p \in P} \sum_{j \in J} \sum_{i \in I} \sum_{t \in T} c_{spjit} T_{pjit} + \sum_{p \in P} \sum_{j \in J} \sum_{t \in T} h_{pjt}^{DC} I_{pjt} \\ + \sum_{p \in P} \sum_{i \in I} \sum_{t \in T} h_{pi}^{VC} I'_{pit} + \sum_{j \in M} \sum_{i \in M} \sum_{v \in V} \sum_{t \in T} \pi_v \Delta_{ji} X_{jivt} \\ + cu \sum_{p \in P} \sum_{k \in K} \sum_{i \in I} \sum_{t \in T} w_{pkit} \left( \frac{R_{pkit}}{d_{pkit}} \right) \quad (12) $$

The objective function (12) includes seven different terms, including 1) fixed cost of opening distribution centers, 2) variable transportation cost of shipping vaccines from the central storage hub to distribution centers, 3) variable transportation cost of shipping vaccines from distribution centers to vaccination centers, 4) variable inventory holding cost of vaccines in distribution centers, 5) variable inventory holding cost of vaccines in vaccination centers, 6) variable transportation cost (per distance) of vehicles, and

7) unmet demand cost. The last term also guarantees the prioritization of different age groups for vaccination. In fact, this term seeks to minimize the total remaining fraction of unmet demand, wherein the fraction of unmet demand of a given vaccination center is weighted by its related priority score.

### 4.2.2. Constraints

In the following, the constraint body of the proposed MVMDLIRP formulation is explained. The first set of constraints, i.e., (13)–(17), determine different quantities in the model from distributed vaccines to unmet demands.

$$ \sum_{p \in P} (1 - \gamma_{pjt}) U_{pjt} \geq \sum_{p \in P} \sum_{i \in I} T_{pjit} Y_{jit} \quad \forall j \in J, t \in T \quad (13) $$

$$ \sum_{k \in K} G_{pkit} \leq \sum_{j \in J} (1 - \mu_{pjit}) T_{pjit} Y_{jit} \quad \forall p \in P, i \in I, t \in T \quad (14) $$

$$ \beta_k d_{pkit} \leq G_{pkit} \leq d_{pkit} \quad \forall p \in P, k \in K, i \in I, t \in T, t \leq L_p \quad (15) $$

$$ \beta_k d_{pkit} + G_{pkit'} \leq G_{pkit} \leq d_{pkit} + G_{pkit'} \quad \forall p \in P, k \in K, i \in I, t \in T, t > L_p, t' = t - L_p, \quad (16) $$

$$ R_{pkit} = R_{pkit, t-1} + d_{pkit} - G_{pkit} \quad \forall p \in P, k \in K, i \in I, t \in T \quad (17) $$

Constraint (13) guarantees that the number of vaccines distributed from each distribution center should be less than the number of arrival vaccines to that center, incorporating the shipping lost rate. Similarly, constraint (14) holds the same condition on arrival and distributed vaccines at each vaccination center. Constraints (15) and (16) determine the amount of demand that needs to be fulfilled for the first and the second dose of each vaccine, respectively, such that the fulfilled demand of a vaccine for its first dose in period t should be re-fulfilled in period t + L. Finally, constraint (17) monitors the unmet demands at each period based on the actual demand and the amount of demand that each vaccination center fulfills.

By involving inventory and shipping lost rates, the set of constraints (18) and (19) preserve the inventory level of vaccines in distribution and vaccination centers, respectively.

$$ I_{pjt}^{DC} = (1 - \alpha_{pjt, t-1}) I_{pjt, t-1}^{DC} + (1 - \gamma_{pjt}) U_{pjt} - \sum_{i \in I} T_{pjit} \quad \forall p \in P, j \in J, t \in T \quad (18) $$

$$ I_{pit}^{VC} = (1 - \lambda_{pit, t-1}) I_{pit, t-1}^{VC} + \sum_{j \in J} (1 - \mu_{pjit}) T_{pjit} - \frac{\sum_{k \in K} G_{pkit}}{1 - \phi_{pit}} \quad \forall p \in P, i \in I, t \in T \quad (19) $$

The set of constraints (20) and (24) ensure the capacity restrictions such that (20) represents the supply capacity for the number of vaccines transferred from the central storage hub; constraints (21) limit the number of vaccines delivered to each vaccination center; constraints (22) ensure the capacity of vehicles when shipping vaccines; and finally, constraints (23) and (24) limit the level of inventory at distribution and vaccination centers, respectively.

$$ \sum_{j \in J} U_{pjt} \leq Q_p^{CH} \quad \forall p \in P, t \in T \quad (20) $$

$$ \sum_{j \in J} (1 - \mu_{pjit}) T_{pjit} \leq \bar{Q}_{pi}^{VC} \quad \forall p \in P, i \in I, t \in T \quad (21) $$

$$ \sum_{i \in I} \sum_{p \in P} T_{pjit} Z K_{jvt} \leq Q_v^V \quad \forall j \in J, v \in V, t \in T \quad (22) $$
<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 colspan="2"><strong>Sets & Indices</strong></td></tr><tr><td>j &isin; J</td><td>Set of distribution centers</td></tr><tr><td>i, i' &isin; I</td><td>Set of vaccination centers</td></tr><tr><td>m &isin; M</td><td>Set of entire nodes in the network; M = I &cup; J</td></tr><tr><td>p &isin; P</td><td>Set of vaccines</td></tr><tr><td>k &isin; K</td><td>Set of population classes</td></tr><tr><td>v &isin; V</td><td>Set of shipping vehicles</td></tr><tr><td>t, t' &isin; T</td><td>Set of time periods</td></tr><tr><td colspan="2"><strong>Parameters</strong></td></tr><tr><td>w<sub>pkit</sub>/d<sub>pkit</sub></td><td>Priority/demand of each population class k at each vaccination center i for receiving vaccine p in period t</td></tr><tr><td>&beta;<sub>k</sub></td><td>Minimum percentage of the demand of each population class k that should be fulfilled</td></tr><tr><td>&delta;</td><td>Service level gap at which a fair distribution of vaccine among different vaccination centers should be established; &delta; &ge; 1. This is also called the maximum ratio among the proportions of fulfilled demand at all vaccination centers</td></tr><tr><td>L<sub>p</sub></td><td>Time interval between two doses of vaccines. In fact, It is imposed to the system that each population class must receive two doses of the same vaccine. Therefore, the demand level in the network is indeed doubled</td></tr><tr><td>f<sub>j</sub></td><td>Fixed cost of establishing the distribution center j</td></tr><tr><td>&gamma;<sub>pjt</sub></td><td>Lost rate of vaccines while shipping the vaccine p from the central storage hub to distribution center j in period t</td></tr><tr><td>&mu;<sub>pjt</sub></td><td>Lost rate of vaccines while shipping the vaccine p from distribution center j to vaccination center i in period t</td></tr><tr><td>c<sub>pj</sub></td><td>Unit shipment cost of vaccine p from the central storage hub to distribution center j</td></tr><tr><td>c<sub>pji</sub></td><td>Unit shipment cost of vaccine p from distribution center j to vaccination center i</td></tr><tr><td>&Delta;<sub>ij</sub></td><td>Distance between each pair of nodes i and j in the network (i, j &isin; M)</td></tr><tr><td>&pi;<sub>v</sub></td><td>Variable (per distance unit) cost of vehicle v</td></tr><tr><td>I<sup>DC</sup><sub>p,j,0</sub>/I<sup>DC</sup><sub>p,j,1,0</sub></td><td>Initial inventory of vaccine p hold at distribution center j/vaccination center i</td></tr><tr><td>I<sup>DC</sup><sub>pjt</sub>/I<sup>DC</sup><sub>pjt,0</sub>/I<sup>DC</sup><sub>pjt,1,0</sub></td><td>Unit inventory holding cost of vaccine p at distribution center j/vaccination center i in period t</td></tr><tr><td>&alpha;<sub>pjt</sub></td><td>Perishing (deterioration) rate of vaccine p at distribution center j at the end of period t</td></tr><tr><td>&lambda;<sub>pji</sub>/&phi;<sub>pit</sub></td><td>Perishing/deterioration rate of vaccine p at vaccination center i at the end of period t</td></tr><tr><td>Q<sup>CH</sup><sub>p</sub></td><td>Storage capacity constraint of vaccine p at the central storage hub</td></tr><tr><td>Q<sup>VC</sup><sub>pj</sub>/Q<sup>VC</sup><sub>pi</sub></td><td>Storage capacity of vaccine p at distribution center j/vaccination center i</td></tr><tr><td>Q<sup>V</sup><sub>v</sub></td><td>Shipping capacity of vehicle v</td></tr><tr><td>Q<sup>VC</sup><sub>pi</sub></td><td>Arrival capacity of vaccine p at vaccination center i</td></tr><tr><td>cu</td><td>Unit cost of unmet demand</td></tr><tr><td>M</td><td>A big enough number</td></tr><tr><td colspan="2"><strong>Decision variables</strong></td></tr><tr><td>I<sup>DC</sup><sub>pjt</sub>/I<sup>VC</sup><sub>pjt</sub></td><td>Inventory level of vaccine p at distribution center j/vaccination center i at the end of period t</td></tr><tr><td>U<sub>pjt</sub></td><td>Amount of vaccine p shipped from the central storage hub to distribution center j in period t</td></tr><tr><td>T<sub>pjt</sub></td><td>Amount of vaccine p shipped from distribution center j to vaccination center i in period t</td></tr><tr><td>G<sub>pkit</sub></td><td>Amount of fulfilled demand of population class k for vaccine p from vaccine center i in period t</td></tr><tr><td>R<sub>pkit</sub></td><td>Amount of demand of population class k for vaccine p left unmet from vaccine center i in period t; R<sub>pkit,0</sub> = 0</td></tr><tr><td>Z<sub>jt</sub></td><td>1 if distribution center j is opened/utilized in period t; 0 otherwise</td></tr><tr><td>Y<sub>jit</sub></td><td>1 if vaccination center i is assigned to distribution center j in period t; 0 otherwise</td></tr><tr><td>X<sub>jitv</sub></td><td>1 if node j is on the route of vehicle v before node i at time period t; 0 otherwise</td></tr><tr><td>ZK<sub>jvt</sub></td><td>1 if node j is on the route of vehicle v at time period t; 0 otherwise</td></tr><tr><td>UU<sub>jitv</sub></td><td>The subtour elimination variable</td></tr></tbody></table>

$$I_{pjt}^{DC} \le Q_{pj}^{DC} Z_{jt} \quad \forall p \in P, j \in J, t \in T \quad (23)$$

$$I_{pit}^{VC} \le Q_{pi}^{VC} \quad \forall p \in P, i \in I, t \in T \quad (24)$$

In order to address the fairness in distributing the vaccines among different vaccination centers, constraints (25) and (26) are introduced such that the former determines the service gap, i.e., the maximal ratio between fulfilled demand proportions to the whole demand for all vaccines among all population classes. The latter also signifies that all vaccination centers that are served by the same distribution center will receive the same proportions of their demand.

$$\frac{\sum_{p \in P} \sum_{k \in K} G_{pkit}}{\sum_{p \in P} \sum_{k \in K} d_{pkit}} \le \frac{\sum_{p \in P} \sum_{k \in K} G_{pkit't}}{\sum_{p \in P} \sum_{k \in K} d_{pkit't}} \cdot \delta \quad \forall i, i' \in I; i \neq i', t \in T \quad (25)$$

$$\frac{\sum_{p \in P} \sum_{k \in K} G_{pkit}}{\sum_{p \in P} \sum_{k \in K} d_{pkit}} \le \frac{\sum_{p \in P} \sum_{k \in K} G_{pkit't}}{\sum_{p \in P} \sum_{k \in K} d_{pkit't}} + (2 - Y_{jit} - Y_{ji't}) \quad \forall i, i' \in I; i \neq i', j \in J, t \in T \quad (26)$$

The remaining set of constraints (27)–(42) are to design the distribution network. Constraint (27) ensures that only opened/established distribution centers can distribute vaccines to vaccination centers. Constraints (28) and (29), respectively, guarantee that each vaccination center is assigned to at least one distribution center at each period as well as more than a single vaccination center is assigned to each distribution center at each period. Constraints (30) and (31) ensure that shipping vehicles serve only established distribution centers and only the established links between them and vaccination centers. Constraint (32) forces that each vehicle can travel through the link from node i to node j if and only if it is decided that node j be on the path of the vehicle v. Constraint (33) guarantees that each vehicle is assigned to at most a single distribution center. Constraint (34) indicates that each node should be on the path of at least one vehicle.

$$Y_{jit} \le Z_{jt} \quad \forall j \in J, i \in I, t \in T \quad (27)$$

$$\sum_{j \in J} Y_{jit} \ge 1 \quad \forall i \in I, t \in T \quad (28)$$

$$\sum_{i \in I} Y_{jit} \ge 1 \quad \forall j \in J, t \in T \quad (29)$$

$$Z_{jt} \ge ZK_{jvt} \quad \forall j \in J, v \in V, t \in T \quad (30)$$
$$Y_{jit} \geq X_{jivt} \quad \forall j \in J, i \in I, v \in V, t \in T \quad (31)$$

$$ZK_{jvt} \geq X_{jivt} \quad \forall j \in M, i \in I, v \in V, t \in T \quad (32)$$

$$\sum_{j \in J} ZK_{jvt} \leq 1 \quad \forall v \in V, t \in T \quad (33)$$

$$\sum_{v \in V} ZK_{jvt} \geq 1 \quad \forall j \in M, t \in T \quad (34)$$

Constraint (35) states that a vaccination center can be assigned to a distribution center if both centers are on the same route. Constraint (36) allows split delivery and ensures that each vaccination center is visited at least once. That is to say, the number of vaccines required by each vaccination center can be fulfilled not as a whole but through smaller deliveries by different vehicles from distribution centers. Constraint (37) ensures that each vehicle should be dispatched from at most a single distribution center. Constraint (38) forces that each distribution center dispatches at least one vehicle. In addition, constraint (39) ensures that at most one vehicle could be assigned for each route. The connectivity condition for vehicles is guaranteed by constraint (40), and constraint (41) indicates that distribution centers can deliver vaccines only to the assigned vaccination centers. Constraint (42) guarantees sub-tour elimination. Finally, constraint (43) determines the types of decision variables.

$$\sum_{j \in M} X_{ijvt} + \sum_{j \in M} X_{j'jvt} - Y_{j'it} \leq 1 \quad \forall i \in I, j' \in J, v \in V, t \in T \quad (35)$$

$$\sum_{j \in M} \sum_{v \in V} X_{jivt} \geq 1 \quad \forall i \in I, t \in T \quad (36)$$

$$\sum_{j \in J} \sum_{i \in I} X_{jivt} \leq 1 \quad \forall v \in V, t \in T \quad (37)$$

$$\sum_{v \in V} \sum_{i \in I} X_{jivt} \geq 1 \quad \forall j \in J, t \in T \quad (38)$$

$$\sum_{v \in V} X_{jivt} \leq 1 \quad \forall i, j \in M: i \neq j, t \in T \quad (39)$$

$$\sum_{i \in M} X_{jivt} - \sum_{i \in M} X_{ijvt} = 0 \quad \forall j \in M, v \in V, t \in T \quad (40)$$

$$T_{pjit} \leq M Y_{jit} \quad \forall p \in P, j \in J, i \in I, t \in T \quad (41)$$

$$UU_{jvt} - UU_{ivt} + |I|X_{jivt} \leq |I| - 1 \quad \forall i, j \in I, v \in V, t \in T \quad (42)$$

$$\begin{cases} U_{pjt}, T_{pjit}, G_{pkit}, R_{pkit}, I_{pjt}, I'_{pit}, UU_{jvt} \geq 0 \\ Z_{jt}, Y_{jit}, ZK_{jvt}, X_{jivt} \in \{0, 1\} \end{cases} \quad \forall p \in P, j \in J, i \in I, v \in V, t \in T \quad (43)$$

The proposed formulation of this section is a non-linear programming model due to the multiplication $T_{pjit}Y_{jit}$ in constraints (13) and (14) as well as the multiplication $T_{pjit}ZK_{jvt}$ in constraint (22). The linearization has been provided in Appendix D.

The employment of SIR models in a vaccine distribution network is twofold: Apriori or Interactive. In an apriori way (Gamchi et al., 2021), the SIR model is used to estimate the dynamic demand with a given vaccination strategy (i.e., fixed rate of vaccination (i.e., u in 1)). In this approach, once the demand is estimated for the whole planning horizon, the demand is given as an input parameter to the optimization model. Furthermore, the decisions in one period do not affect the estimated demand of further

periods, except the unmet demand accumulated over periods. In an interactive way (Bertsimas et al., 2022), the SIR model interacts with an optimization model in each period. In this way, the vaccination rate u is a decision variable, and the optimization model determines it according to the technical constraints of the distribution network. In this way, the decisions in each period directly affect the estimated demand for further periods, and no unmet demand is considered. Accordingly, the SIR model is updated at the end of each period based on the output of the optimization model in that period, and it estimates the demand for the next period. The new demand is then given to the optimization model for further vaccination decisions. These optimization models, which are executed for each period, are called myopic models. In this paper, we have employed the apriori approach with a given vaccination rate u, where the SIR model is executed once throughout the whole planning horizon.

## 5. Benders decomposition - BD

This section employs a Benders decomposition algorithm (Benders, 1962) to solve the proposed model for the MVMDLIRP. Benders (1962) proposed the BD algorithm to deal with problems with complicated variables, where it decentralizes the structure of a problem to provide an easier formation of the problem to reduce the computational burden (Abou-Ismail, 2020; Rahmaniani et al., 2018). From then until now, it has been gaining prestige as one of the most powerful exact algorithms to solve a broad range of NP-hard problems, such as vehicle routing (Corrèa et al., 2007), facility location (Boland et al., 2016), logistics network design (Cordeau et al., 2006), transportation (Coelho & Laporte, 2014), and inventory vehicle routing (Alkaabneh et al., 2020). What is more, many studies in the literature have used the BD algorithm to solve the vehicle routing problem and its derivatives, such as location-routing (Çalık et al., 2021), location-inventory-routing (Zheng et al., 2019), and production-inventory-routing (Cordeau et al., 2015) problems. As a result, several efficient valid inequalities have been presented for these problems that could limit the solution space and improve the convergence speed of the BD. Hence, since the proposed model in the current research is one of the derivatives of the vehicle routing problem, the BD algorithm is employed to solve the proposed model.

### 5.1. BD in general

The idea behind the BD algorithm is to divide the original problem into a master problem and a set of sub-problems, with the hope that the decomposed problem be easier to be solved (Abou-Ismail, 2020; Rahmaniani et al., 2018). The master problem incorporates a part of the original model with only integer variables. On the other hand, sub-problems are formed by applying programming duality over the rest of the original model, knowing that the value of the integer variables is given. Along with the integer variables in the master problem, an artificial variable is considered that describes a lower bound (upper bound) on the sub-problems' objective function for a minimization (maximization) problem. Through the BD algorithm (see Fig. 2), master and sub-problems are solved iteratively, such that the master problem is solved first and the values of the integer variables are determined. Next, the sub-problems are solved for the given value of the integer variables from the master problem. Finally, a feasibility cut (based on the values of the sub-problems' variables) is added to the master problem if some sub-problems are infeasible or 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. This procedure is repeated for further iterations until a
<figure><img src="image_5.png" /><figcaption>Fig. 2. BD's flowchart.</figcaption></figure>

stopping criterion is met (i.e., the maximum number of iterations or a threshold for the gap between the obtained lower and upper bounds). Afterward, a set of accelerators, in terms of inequalities, are employed to enhance BD's performance, the so-called accelerated BD (ABD) algorithm. In the following subsection, a set of valid inequalities are introduced.

## 5.2. Valid inequalities

The first set of valid inequalities (44)–(46) have been widely used in the literature for inventory and location routing problems (Coelho & Laporte, 2014; Darvish et al., 2019). These inequalities are employed to empower the linear programming relaxation of the MVMDLIRP formulation.

$$ \sum_{v \in V} \sum_{t' = t_1}^{t_2} ZK_{livt'} \geq \frac{\sum_{t' = t_1}^{t_2} \sum_{p \in P} \sum_{k \in K} d_{kit'} - \sum_{p \in P} \sum_{k \in K} d_{kit}}{\sum_{t' = t_1}^{t_2} \sum_{p \in P} \sum_{k \in K} d_{kit'}} \quad \forall i \in I, t_1, t_2 \in T; t_2 \geq t_1 \quad (44) $$

$$ X_{livt} = 0 \quad \forall i \in M, v \in V, t \in T \quad (45) $$

$$ 2 \cdot ZK_{jvt} \leq \sum_{i \in I} X_{jivt} + \sum_{i \in I} X_{ijvt} \quad \forall j \in J, v \in V, t \in T \quad (46) $$

The idea behind the valid inequality (44) relates to whether the inventory held at each vaccination center at each period is sufficient to fulfill future demands. More precisely, if the inventory held in period $t_1$ by vaccination center $i$ is sufficient to fulfill its demand for periods $[t_1, t_2]$, then no delivery of vaccines to vaccination center $i$ is required. Indeed, if $\sum_{p \in P} I_{pit_1}^{VC} \geq \sum_{t' = t_1}^{t_2} \sum_{p \in P} \sum_{k \in K} d_{kit'}$,

then $\sum_{v \in V} \sum_{t' = t_1}^{t_2} ZK_{livt'} \geq 0$. On the other hand, if the inventory is not sufficient to fulfill future demands, then a delivery must take place. Moreover, the valid inequality (45) eliminates links between a node and itself. Finally, the valid inequality (46) forces that one route can start and end at a distribution center in each period if the distribution center is opened/established.

Another type of valid inequality, called knapsack inequality (Santoso et al., 2005), is used in the proposed ABD algorithm. The knapsack inequality is added to the master problem with the aim of accelerating the branch-and-bound process of the solver. Indeed, the convergence speed of Benders decomposition can be improved by adding the knapsack inequality, which assists progressive solvers like Gurobi in deriving a range of valid inequalities.

$$ LB^n \leq \sum_{j \in J} \sum_{t \in T} f_j Z_{jt} + \sum_{j \in J} \sum_{i \in M} \sum_{v \in V} \sum_{t \in T} \pi_v \Delta_{ij} X_{jivt} + \zeta \quad (47) $$

where $LB^n$ represents the best lower bound found so far (i.e., the objective function value (the best lower bound found by the solver) in the previous iteration in case of feasibility (infeasibility)), and $\zeta$ is an additional variable representing the total cost except fixed opening cost distribution centers as well as the variable shipping cost of the vehicles.

The Benders reformulation of the MVMDLIRP is provided in Appendix E.

# 6. Computational experiments

This section presents the computational results of the designed vaccine distribution network under demand dynamicity. In the following, Section 6.1 designs a set of experiments based on a real case study, the vaccination campaign during the COVID-19 pandemic in France. Next, Section 6.2 presents a comparative analysis between the proposed ABD algorithm and the Gurobi solver
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** The case study representing different regions with corresponding population indices.</figcaption></figure>

(GRB). Finally, a set of comprehensive sensitivity analyses are conducted in Section 7 to investigate the behavior of the proposed MVMDLIRP model to any changes of its input parameters. These sensitivity analyses help to draw valuable insights for health policymakers (Fig. 3).

## 6.1. Experimental design and case study

This section describes a real case study (see Fig. 3) on the vaccination campaign in France. France encompasses a total number of 18 regions (i.e., 13 metropolitan and 5 overseas regions) and is divided into 101 departments (i.e., 96 metropolitan and 5 overseas departments). Among all, a number of 12 metropolitan regions (see Fig. 3(a) with the total eligible-for-vaccination population of each region in parentheses) are considered to be the potential location for establishing distribution centers (i.e., $|I| = 12$), and a number of 80 metropolitan departments are considered as the vaccination centers (i.e., $|I| = 80$) or the areas that have a particular demand for vaccines. The metropolitan regions in France possess different properties (e.g., demographic, educational, social, economic, etc.), and these disparities among regions make vaccination a complex work for France's health policymakers. A summary of these properties has been provided per region in Table F.1 in Appendix F. In this case study, the population in each area is divided into four age-based classes (i.e., $|K| = 4$; ages: "18–49", "50–64", "65–74", "> 75"), and each area includes a different number of individuals per each age-based class (see Fig. 3(b)). The central storage hub of the purchased vaccines is located in the capital (i.e., region de-France). A total number of 24 refrigerated trucks are considered (i.e., $|V| = 24$) to ship the vaccines from distribution centers to vaccination centers.

The dynamic demand for vaccines in each populated area is obtained by running the proposed system dynamics model (9) for that area, where the birth and death rates as well as the population are different for each area. Detailed information on all other parameters of the adjusted SIR model and the MVMDLIRP model is provided in Table 4. Except for the birth and death rates as well as the population in each set of populated areas, other parameters of the adjusted SIR model are considered similar for the whole of France (Angeli et al., 2021; OPECST, 2021).

In Table 4, the unit shipment costs $c$ and $cs$ are calculated based on the rent price of refrigerated trucks divided by their holding capacity (per pallet of 10k vaccines). Depending on the type and capacity of the truck, this cost varied from 10 to 40 euros. In addition, the variable (per distance) shipment cost $\pi$ represents the vehicles' traveling cost (e.g., fuel). The unit holding costs $h^{DC/VC}$

are also proportional to the shipment cost of vaccines, $c$, since the shipment cost is indeed the cost of utilizing mobile refrigerators to hold vaccines. Furthermore, the unit unmet demand cost is calculated as the expected cost that an unvaccinated individual imposes on the health system (e.g., treatments, medicines, hospitalization, test costs, etc.).

Figure 4 illustrates the outcomes of the adjusted SIR model in determining the state variables for two distinct populated areas, the demand points of vaccines. As can be observed, the dynamicity of the state variables in the two areas 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 adjusted SIR model) to determine the demand for vaccines. Remarkably, since the pandemic is not over at the time of this research, the state variable $uS$ (Fig. 4f) is used to estimate the demand for vaccines in each populated area.

The performance of the proposed ABD algorithm is tested through a set of 36 test problems with a different number of vaccination centers and time periods. All test problems except the last one are indeed smaller parts of the case study with fewer vaccination centers or fewer time periods. Since the adjusted SIR model is run for each populated area for a time horizon of 365 days, and it is also impossible and illogical to plan the MVMDLIRP daily, each period of the MVMDLIRP model represents a week and the demand for each period is the cumulative demand of the whole week.

Both the MVMDLIRP model and the ABD algorithm were coded in Python 3.8 using the Gurobi 9.1.2 library, and all experiments were done on a server containing four Intel XEON processors with 5 gigabytes of RAM memory running at 2.3 gigahertz. Furthermore, two stopping criteria were considered, when executing the ABD algorithm and GRB, as 1) a gap of 1% and 2) a CPU time of 7200 seconds. For the ABD algorithm, the first criterion is the gap (%) between the obtained lower-bound and upper-bound at each algorithmic iteration; and for GRB, it is the gap (%) between the best-found solution and the current obtained solution.

## 6.2. Numerical results

This section conducts a comparative analysis of the performance of the proposed ABD algorithm with valid inequalities (ABD-w-VI) and without valid inequalities (ABD-wo-VI) compared to GRB. Table 5 shows the results of this comparison for 36 test problems with a different number of vaccination centers (i.e., column $|I|$) and different periods (i.e., column $|T|$) in terms of 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 value of the test problem #1. The goal of this comparison is to evaluate the performance of the proposed ABD algorithm compared to GRB and evaluate the benefits of valid inequalities for the proposed ABD algorithm.

As can be seen in Table 5, both ABD-w-VI and ABD-wo-VI algorithms have been able to obtain solutions with a gap of less than 1% for all test problems before reaching the maximum allowable CPU time of 7200 seconds, even for larger test problems. However, GRB has been unable to find the optimal solution for test problems ($|I|, |T| \in \{(60, 20-24), (70, 14-18), (80, 8-12)\}$) under a computational effort of 7200 seconds. Furthermore, GRB has not even been able to find a feasible solution for larger test problems ($|I|, |T| \in \{(70, 20-24), (80, 14-24)\}$) under the limited computational time. Moreover, columns "Gap" compare the performance of the ABD-wo-VI algorithm and GRB with the ABD-w-VI algorithm in terms of objective function gap (%) and the computational time ratio. Regarding the objective function gap, we only compare the performance of GRB, since the ABD-wo-VI has reached the same solutions as those of the ABD-w-VI algorithm. In terms of objective function gap, the ABD-w-VI algorithm has obtained much better
<figure><img src="image_7.png" /><figcaption>Fig. 4. The outcomes of the dynamic optimization model.</figcaption></figure>

<table><caption>**Table 4**
Input data.</caption><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>80</td><td>λ</td><td>[0.05, 0.1]</td><td>L (week)</td><td>[2, 6]</td></tr><tr><td>J</td><td>12</td><td>Q<sup>CH</sup></td><td>[2M, 5M]</td><td>cu (€)</td><td>2000</td></tr><tr><td>P</td><td>3</td><td>Q<sup>DC</sup></td><td>[0.5M, 1M]</td><td>β<sub>1</sub></td><td>3.70E-6</td></tr><tr><td>K</td><td>4</td><td>Q<sup>VC</sup></td><td>[2E + 4, 10E + 4]</td><td>β<sub>2</sub></td><td>1.48E-5</td></tr><tr><td>V</td><td>24</td><td>Q<sup>VC</sup></td><td>[5E + 4, 2E + 5]</td><td>β<sub>3</sub></td><td>7.00E-4</td></tr><tr><td>T</td><td>24</td><td>Q<sup>V</sup></td><td>[0.3M, 0.5M]</td><td>k</td><td>0.05</td></tr><tr><td>w</td><td>[0,0.5]</td><td>f (€)</td><td>[0.2M, 0.5M]</td><td>δ<sub>1</sub></td><td>[1E-4, 1E-3]</td></tr><tr><td>β</td><td>[0.1,0.25]</td><td>c (€)</td><td>[10, 40]</td><td>δ<sub>2</sub></td><td>[2E-4, 2E-3]</td></tr><tr><td>δ</td><td>1.5</td><td>cs (€)</td><td>[5%, 20%] × c</td><td>δ<sub>3</sub></td><td>[5E-4, 5E-3]</td></tr><tr><td>γ</td><td>[0.05, 0.1]</td><td>h<sup>DC</sup> (€)</td><td>[5%, 10%] × c</td><td>τ</td><td>0.2</td></tr><tr><td>μ</td><td>[0.05, 0.1]</td><td>h<sup>VC</sup> (€)</td><td>[5%, 10%] × c</td><td>ν<sub>1</sub></td><td>0.2</td></tr><tr><td>α</td><td>[0.05, 0.1]</td><td>π (€)</td><td>1% × c</td><td>ν<sub>2</sub></td><td>0.05</td></tr></tbody></table>

ter solutions than GRB with a mean gap of 16% over test problems for which GRB has reached at least a feasible solution. In terms of computational time, a paired comparative analysis has only been done over the test problems for which both corresponding methods have obtained the optimal solution with a gap of less than 1% within the computational time of 7200 seconds.

It can be observed that the ABD-w-VI algorithm is faster than both the ABD-wo-VI algorithm and GRB for all the corresponding test problems with mean ratios of 1.75 and 12, indicating that the ABD-w-VI algorithm is, on average, 1.752 and 12 times faster than the ABD-wo-VI algorithm and GRB for solving the MVMDLIRP model. To better show the difference between algorithms, Fig. 5 illustrates the comparison between three solution techniques (i.e., ABD-w/wo-VI and GRB) in terms of both objective function values (Fig. 5(a)) and the computational time (Fig. 5(b)).

In order to explore the detail of results, Table 6 reports different information for all test problems, including the number of opened distribution centers (“OCD” in %), the unmet demand (“UMD” in %) for each age-based class, and the contribution of each vaccine to the vaccination of each age-based class (“CVG” in %).

Looking at column “ODC (%)", we observe that the larger the size of the problem (i.e., the number of vaccination centers and

planning periods), the higher the number of distribution centers to fulfill the demand. Overall, about 82% of the total number of distribution centers have been used to fulfill the demand in the distribution network. The reason for the value “0” for instances with $|T| = 2$ (i.e., two weeks) is that the initial inventory of vaccines in the vaccination centers is sufficient to fulfill the demand for such a short planning horizon. When solving the original case study with $|I| = 80$ for planning horizons of $|T| \ge 8$ (i.e., greater than 8 weeks), the whole capacity of the distribution network is utilized by opening 100% of the distribution centers.

The results of columns “UMD (%)” indicate that, in overall, 6.4%, 9.4%, 13.8%, and 23.1% of the demand have remained unmet for age-based classes “> 75”, “65–74”, “50–64”, and “18–49”, respectively. Comparing age-based classes shows that a lower percentage of the population remains unvaccinated, from elders to younger individuals. This comes from the prioritization of the age-based classes in the objective function. More importantly, it can be seen that the unmet demand of elders remains fixed when increasing the number of vaccination centers (or indirectly increasing the potential demand for vaccination); however, the unmet demand of younger individuals increases. The reason is that no matter how much the population is, the system put its all efforts into vac-
<table><caption>**Table 5**
ABD vs. GRB.</caption><thead><tr><th rowspan="2">#</th><th rowspan="2">|I|</th><th rowspan="2">|T|</th><th colspan="3">Obj. values</th><th colspan="3">Time (seconds)</th><th colspan="3">Gap (... vs. ABD-w-VI)</th></tr><tr><th>ABD-w-VI</th><th>ABD-wo-VI</th><th>GRB</th><th>ABD-w-VI</th><th>ABD-wo-VI</th><th>GRB</th><th>Obj. (%)</th><th>Time (ratio)</th></tr></thead><tbody><tr><td>1</td><td>60</td><td>2</td><td>1.000X</td><td>1.000X</td><td>1.000X</td><td>57</td><td>62</td><td>130</td><td>0</td><td>1.09</td><td>2.29</td></tr><tr><td>2</td><td>60</td><td>4</td><td>2.066X</td><td>2.066X</td><td>2.066X</td><td>75</td><td>83</td><td>393</td><td>0</td><td>1.1</td><td>5.28</td></tr><tr><td>3</td><td>60</td><td>6</td><td>3.152X</td><td>3.152X</td><td>3.152X</td><td>88</td><td>97</td><td>749</td><td>0</td><td>1.1</td><td>8.48</td></tr><tr><td>4</td><td>60</td><td>8</td><td>4.251X</td><td>4.251X</td><td>4.251X</td><td>104</td><td>118</td><td>1818</td><td>0</td><td>1.13</td><td>17.52</td></tr><tr><td>5</td><td>60</td><td>10</td><td>5.361X</td><td>5.361X</td><td>5.361X</td><td>132</td><td>168</td><td>2560</td><td>0</td><td>1.27</td><td>19.37</td></tr><tr><td>6</td><td>60</td><td>12</td><td>6.481X</td><td>6.481X</td><td>6.481X</td><td>190</td><td>314</td><td>3496</td><td>0</td><td>1.65</td><td>18.45</td></tr><tr><td>7</td><td>60</td><td>14</td><td>7.610X</td><td>7.610X</td><td>7.610X</td><td>206</td><td>340</td><td>4386</td><td>0</td><td>1.65</td><td>21.30</td></tr><tr><td>8</td><td>60</td><td>16</td><td>8.752X</td><td>8.752X</td><td>8.752X</td><td>282</td><td>488</td><td>5976</td><td>0</td><td>1.73</td><td>21.21</td></tr><tr><td>9</td><td>60</td><td>18</td><td>9.906X</td><td>9.906X</td><td>9.906X</td><td>347</td><td>618</td><td>6666</td><td>0</td><td>1.78</td><td>19.22</td></tr><tr><td>10</td><td>60</td><td>20</td><td>11.075X</td><td>11.075X</td><td>14.790X</td><td>389</td><td>700</td><td>7200</td><td>33.54</td><td>1.8</td><td>*</td></tr><tr><td>11</td><td>60</td><td>22</td><td>12.258X</td><td>12.258X</td><td>17.430X</td><td>470</td><td>879</td><td>7200</td><td>42.19</td><td>1.87</td><td>*</td></tr><tr><td>12</td><td>60</td><td>24</td><td>13.457X</td><td>13.457X</td><td>20.761X</td><td>556</td><td>1073</td><td>7200</td><td>54.27</td><td>1.93</td><td>*</td></tr><tr><td>13</td><td>70</td><td>2</td><td>1.005X</td><td>1.005X</td><td>1.005X</td><td>46</td><td>58</td><td>193</td><td>0</td><td>1.27</td><td>4.20</td></tr><tr><td>14</td><td>70</td><td>4</td><td>2.075X</td><td>2.075X</td><td>2.075X</td><td>66</td><td>89.1</td><td>531</td><td>0</td><td>1.35</td><td>8.02</td></tr><tr><td>15</td><td>70</td><td>6</td><td>3.158X</td><td>3.158X</td><td>3.158X</td><td>139</td><td>204</td><td>1154</td><td>0</td><td>1.47</td><td>8.29</td></tr><tr><td>16</td><td>70</td><td>8</td><td>4.257X</td><td>4.257X</td><td>4.257X</td><td>187</td><td>275</td><td>2283</td><td>0</td><td>1.47</td><td>12.23</td></tr><tr><td>17</td><td>70</td><td>10</td><td>5.367X</td><td>5.367X</td><td>5.367X</td><td>191</td><td>302</td><td>3608</td><td>0</td><td>1.58</td><td>18.92</td></tr><tr><td>18</td><td>70</td><td>12</td><td>6.489X</td><td>6.489X</td><td>6.489X</td><td>378</td><td>631</td><td>6967</td><td>0</td><td>1.67</td><td>18.43</td></tr><tr><td>19</td><td>70</td><td>14</td><td>7.624X</td><td>7.624X</td><td>10.636X</td><td>422</td><td>709</td><td>7200</td><td>39.51</td><td>1.68</td><td>*</td></tr><tr><td>20</td><td>70</td><td>16</td><td>8.772X</td><td>8.772X</td><td>12.837X</td><td>487</td><td>891</td><td>7200</td><td>46.34</td><td>1.83</td><td>*</td></tr><tr><td>21</td><td>70</td><td>18</td><td>9.934X</td><td>9.934X</td><td>14.631X</td><td>516</td><td>970</td><td>7200</td><td>47.29</td><td>1.88</td><td>*</td></tr><tr><td>22</td><td>70</td><td>20</td><td>11.112X</td><td>11.112X</td><td>**</td><td>537</td><td>1208</td><td>**</td><td>**</td><td>2.25</td><td>**</td></tr><tr><td>23</td><td>70</td><td>22</td><td>12.308X</td><td>12.308X</td><td>**</td><td>685</td><td>1569</td><td>**</td><td>**</td><td>2.29</td><td>**</td></tr><tr><td>24</td><td>70</td><td>24</td><td>13.467X</td><td>13.467X</td><td>**</td><td>727</td><td>1687</td><td>**</td><td>**</td><td>2.32</td><td>**</td></tr><tr><td>25</td><td>80</td><td>2</td><td>1.007X</td><td>1.007X</td><td>1.007X</td><td>56</td><td>77</td><td>296</td><td>0</td><td>1.37</td><td>5.29</td></tr><tr><td>26</td><td>80</td><td>4</td><td>2.080X</td><td>2.080X</td><td>2.080X</td><td>121</td><td>190</td><td>1855</td><td>0</td><td>1.57</td><td>15.32</td></tr><tr><td>27</td><td>80</td><td>6</td><td>3.168X</td><td>3.168X</td><td>3.168X</td><td>193</td><td>326</td><td>4558</td><td>0</td><td>1.69</td><td>23.60</td></tr><tr><td>28</td><td>80</td><td>8</td><td>4.275X</td><td>4.275X</td><td>6.648X</td><td>283</td><td>501</td><td>7200</td><td>55.49</td><td>1.77</td><td>*</td></tr><tr><td>29</td><td>80</td><td>10</td><td>5.416X</td><td>5.416X</td><td>8.570X</td><td>385</td><td>689</td><td>7200</td><td>58.24</td><td>1.79</td><td>*</td></tr><tr><td>30</td><td>80</td><td>12</td><td>6.497X</td><td>6.497X</td><td>9.704X</td><td>482</td><td>877</td><td>7200</td><td>49.35</td><td>1.82</td><td>*</td></tr><tr><td>31</td><td>80</td><td>14</td><td>7.638X</td><td>7.638X</td><td>**</td><td>598</td><td>1244</td><td>**</td><td>**</td><td>2.08</td><td>**</td></tr><tr><td>32</td><td>80</td><td>16</td><td>8.786X</td><td>8.786X</td><td>**</td><td>711</td><td>1536</td><td>**</td><td>**</td><td>2.16</td><td>**</td></tr><tr><td>33</td><td>80</td><td>18</td><td>9.945X</td><td>9.945X</td><td>**</td><td>1026</td><td>2319</td><td>**</td><td>**</td><td>2.26</td><td>**</td></tr><tr><td>34</td><td>80</td><td>20</td><td>11.118X</td><td>11.118X</td><td>**</td><td>1359</td><td>3126</td><td>**</td><td>**</td><td>2.3</td><td>**</td></tr><tr><td>35</td><td>80</td><td>22</td><td>12.314X</td><td>12.314X</td><td>**</td><td>1958</td><td>4836</td><td>**</td><td>**</td><td>2.47</td><td>**</td></tr><tr><td>36</td><td>80</td><td>24</td><td>13.474X</td><td>13.474X</td><td>**</td><td>2842</td><td>7133</td><td>**</td><td>**</td><td>2.51</td><td>**</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Comparison between ABD and GRB.</figcaption></figure>

<footnote>* Time limit of 7200 seconds reached **Average**: ≈ 16 ≈ 1.75 ≈ 12.</footnote>

<footnote>* No feasible solution after 7200 seconds.</footnote>
<table><caption>**Table 6**
Detailed results.</caption><thead><tr><th rowspan="2">#</th><th rowspan="2">|I|</th><th rowspan="2">|T|</th><th rowspan="2">ODC (%)</th><th colspan="4">UMD (%)</th><th colspan="12">CVG (%)</th></tr><tr><th>"18-49"</th><th>"50-64"</th><th>"65-74"</th><th>">75"</th><th colspan="3">"18-49"</th><th colspan="3">"50-64"</th><th colspan="3">"65-74"</th><th colspan="3">">75"</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>AZ</th><th>MO</th><th>Pf</th><th>AZ</th><th>MO</th><th>Pf</th><th>AZ</th><th>MO</th><th>Pf</th><th>AZ</th><th>MO</th><th>Pf</th></tr></thead><tbody><tr><td>1</td><td>60</td><td>2</td><td>0</td><td>25.0</td><td>23.9</td><td>23.0</td><td>21.5</td><td>7.9</td><td>13.9</td><td>78.2</td><td>7.9</td><td>13.9</td><td>78.2</td><td>8.3</td><td>13.8</td><td>77.9</td><td>8.9</td><td>13.7</td><td>77.4</td></tr><tr><td>2</td><td>60</td><td>4</td><td>50</td><td>35.2</td><td>31.0</td><td>26.9</td><td>22.1</td><td>8.1</td><td>16.1</td><td>75.8</td><td>6.3</td><td>15.8</td><td>77.9</td><td>7.8</td><td>15.3</td><td>76.8</td><td>8.8</td><td>15.5</td><td>75.7</td></tr><tr><td>3</td><td>60</td><td>6</td><td>67</td><td>34.1</td><td>27.4</td><td>20.9</td><td>14.9</td><td>7.8</td><td>21.9</td><td>70.4</td><td>6.0</td><td>19.7</td><td>74.3</td><td>7.1</td><td>19.4</td><td>73.5</td><td>7.5</td><td>18.9</td><td>73.6</td></tr><tr><td>4</td><td>60</td><td>8</td><td>75</td><td>31.7</td><td>22.3</td><td>15.3</td><td>9.1</td><td>10.6</td><td>19.8</td><td>69.6</td><td>7.4</td><td>19.8</td><td>72.7</td><td>6.1</td><td>19.2</td><td>74.8</td><td>6.0</td><td>17.9</td><td>76.0</td></tr><tr><td>5</td><td>60</td><td>10</td><td>80</td><td>25.0</td><td>16.1</td><td>8.8</td><td>2.2</td><td>11.1</td><td>18.6</td><td>70.3</td><td>6.8</td><td>20.2</td><td>72.9</td><td>5.2</td><td>18.9</td><td>75.9</td><td>5.2</td><td>16.8</td><td>78.0</td></tr><tr><td>6</td><td>60</td><td>12</td><td>83</td><td>21.9</td><td>12.1</td><td>4.1</td><td>0.0</td><td>12.0</td><td>18.2</td><td>69.8</td><td>5.9</td><td>20.6</td><td>73.5</td><td>4.7</td><td>18.5</td><td>76.9</td><td>5.0</td><td>16.3</td><td>78.7</td></tr><tr><td>7</td><td>60</td><td>14</td><td>86</td><td>18.7</td><td>9.2</td><td>0.0</td><td>0.0</td><td>13.2</td><td>17.6</td><td>69.2</td><td>4.9</td><td>21.5</td><td>73.6</td><td>4.3</td><td>18.6</td><td>77.1</td><td>5.1</td><td>16.1</td><td>78.7</td></tr><tr><td>8</td><td>60</td><td>16</td><td>88</td><td>15.7</td><td>4.2</td><td>0.0</td><td>0.0</td><td>14.2</td><td>17.5</td><td>68.3</td><td>4.4</td><td>22.4</td><td>73.2</td><td>4.4</td><td>19.2</td><td>76.4</td><td>5.4</td><td>15.8</td><td>78.8</td></tr><tr><td>9</td><td>60</td><td>18</td><td>89</td><td>11.9</td><td>0.7</td><td>0.0</td><td>0.0</td><td>15.1</td><td>17.0</td><td>67.8</td><td>3.6</td><td>23.7</td><td>72.7</td><td>4.5</td><td>19.2</td><td>76.4</td><td>5.8</td><td>15.3</td><td>78.8</td></tr><tr><td>10</td><td>60</td><td>20</td><td>90</td><td>11.9</td><td>0.0</td><td>0.0</td><td>0.0</td><td>16.2</td><td>16.4</td><td>67.4</td><td>3.4</td><td>23.9</td><td>72.7</td><td>4.5</td><td>19.0</td><td>76.5</td><td>6.2</td><td>15.4</td><td>78.4</td></tr><tr><td>11</td><td>60</td><td>22</td><td>91</td><td>10.1</td><td>0.0</td><td>0.0</td><td>0.0</td><td>17.7</td><td>14.8</td><td>67.5</td><td>3.4</td><td>23.7</td><td>72.8</td><td>4.7</td><td>18.4</td><td>76.8</td><td>6.5</td><td>15.0</td><td>78.5</td></tr><tr><td>12</td><td>60</td><td>24</td><td>92</td><td>2.6</td><td>0.0</td><td>0.0</td><td>0.0</td><td>19.2</td><td>14.0</td><td>66.9</td><td>3.4</td><td>23.8</td><td>72.8</td><td>4.9</td><td>18.6</td><td>76.5</td><td>6.7</td><td>15.0</td><td>78.3</td></tr><tr><td>13</td><td>70</td><td>2</td><td>0</td><td>25.1</td><td>24.1</td><td>23.2</td><td>21.6</td><td>7.9</td><td>13.9</td><td>78.2</td><td>7.9</td><td>13.9</td><td>78.2</td><td>8.3</td><td>13.8</td><td>77.8</td><td>8.8</td><td>13.7</td><td>77.5</td></tr><tr><td>14</td><td>70</td><td>4</td><td>60</td><td>35.5</td><td>31.4</td><td>27.4</td><td>22.6</td><td>8.4</td><td>12.6</td><td>79.0</td><td>6.5</td><td>13.3</td><td>80.1</td><td>8.0</td><td>12.8</td><td>79.2</td><td>9.3</td><td>12.2</td><td>78.5</td></tr><tr><td>15</td><td>70</td><td>6</td><td>72</td><td>36.9</td><td>29.9</td><td>23.2</td><td>16.4</td><td>8.7</td><td>12.0</td><td>79.2</td><td>6.4</td><td>13.3</td><td>80.3</td><td>7.8</td><td>12.1</td><td>80.1</td><td>8.2</td><td>12.3</td><td>79.5</td></tr><tr><td>16</td><td>70</td><td>8</td><td>84</td><td>33.2</td><td>25.0</td><td>16.5</td><td>9.2</td><td>8.6</td><td>17.1</td><td>74.3</td><td>5.7</td><td>15.1</td><td>79.2</td><td>6.7</td><td>15.1</td><td>78.3</td><td>7.1</td><td>15.2</td><td>77.7</td></tr><tr><td>17</td><td>70</td><td>10</td><td>93</td><td>30.7</td><td>19.2</td><td>10.3</td><td>2.3</td><td>8.9</td><td>18.1</td><td>72.9</td><td>5.3</td><td>17.7</td><td>77.0</td><td>5.9</td><td>16.7</td><td>77.4</td><td>6.6</td><td>15.8</td><td>77.6</td></tr><tr><td>18</td><td>70</td><td>12</td><td>96</td><td>26.3</td><td>14.0</td><td>4.4</td><td>0.0</td><td>9.5</td><td>19.1</td><td>71.4</td><td>5.0</td><td>18.2</td><td>76.8</td><td>5.5</td><td>17.2</td><td>77.3</td><td>6.9</td><td>16.4</td><td>76.7</td></tr><tr><td>19</td><td>70</td><td>14</td><td>98</td><td>22.0</td><td>9.4</td><td>0.0</td><td>0.0</td><td>10.6</td><td>20.0</td><td>69.4</td><td>4.6</td><td>18.7</td><td>76.7</td><td>5.5</td><td>17.9</td><td>76.6</td><td>7.3</td><td>15.6</td><td>77.1</td></tr><tr><td>20</td><td>70</td><td>16</td><td>100</td><td>18.5</td><td>4.6</td><td>0.0</td><td>0.0</td><td>12.1</td><td>19.6</td><td>68.3</td><td>4.3</td><td>19.4</td><td>76.4</td><td>5.8</td><td>17.8</td><td>76.4</td><td>7.7</td><td>15.2</td><td>77.1</td></tr><tr><td>21</td><td>70</td><td>18</td><td>100</td><td>14.7</td><td>0.6</td><td>0.0</td><td>0.0</td><td>15.5</td><td>17.6</td><td>66.9</td><td>4.1</td><td>19.7</td><td>76.2</td><td>6.1</td><td>17.5</td><td>76.3</td><td>7.9</td><td>14.5</td><td>77.7</td></tr><tr><td>22</td><td>70</td><td>20</td><td>100</td><td>12.7</td><td>0.0</td><td>0.0</td><td>0.0</td><td>17.5</td><td>16.3</td><td>66.2</td><td>4.1</td><td>20.0</td><td>75.9</td><td>6.2</td><td>17.3</td><td>76.4</td><td>7.9</td><td>14.1</td><td>78.0</td></tr><tr><td>23</td><td>70</td><td>22</td><td>100</td><td>11.5</td><td>0.0</td><td>0.0</td><td>0.0</td><td>19.1</td><td>15.2</td><td>65.7</td><td>4.3</td><td>20.1</td><td>75.7</td><td>6.7</td><td>17.4</td><td>75.9</td><td>7.9</td><td>13.9</td><td>78.2</td></tr><tr><td>24</td><td>70</td><td>24</td><td>100</td><td>8.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>21.1</td><td>14.2</td><td>64.7</td><td>5.3</td><td>19.1</td><td>75.7</td><td>8.7</td><td>19.4</td><td>71.9</td><td>7.9</td><td>12.9</td><td>79.2</td></tr><tr><td>25</td><td>80</td><td>2</td><td>0</td><td>27.4</td><td>26.3</td><td>24.7</td><td>22.9</td><td>7.9</td><td>13.9</td><td>78.2</td><td>7.9</td><td>13.9</td><td>78.2</td><td>8.4</td><td>13.8</td><td>77.8</td><td>9.0</td><td>13.7</td><td>77.3</td></tr><tr><td>26</td><td>80</td><td>4</td><td>75</td><td>38.5</td><td>34.0</td><td>28.9</td><td>23.6</td><td>8.4</td><td>12.5</td><td>79.0</td><td>6.7</td><td>13.3</td><td>80.0</td><td>8.4</td><td>12.6</td><td>79.0</td><td>9.6</td><td>12.3</td><td>78.0</td></tr><tr><td>27</td><td>80</td><td>6</td><td>88</td><td>37.5</td><td>31.9</td><td>24.0</td><td>16.8</td><td>9.0</td><td>10.9</td><td>80.1</td><td>6.7</td><td>12.7</td><td>80.6</td><td>8.4</td><td>11.6</td><td>80.0</td><td>8.5</td><td>12.5</td><td>79.0</td></tr><tr><td>28</td><td>80</td><td>8</td><td>100</td><td>36.2</td><td>29.0</td><td>19.5</td><td>12.2</td><td>11.0</td><td>21.8</td><td>67.1</td><td>7.5</td><td>19.9</td><td>72.7</td><td>6.3</td><td>19.6</td><td>74.1</td><td>6.6</td><td>19.4</td><td>74.0</td></tr><tr><td>29</td><td>80</td><td>10</td><td>100</td><td>34.7</td><td>25.2</td><td>15.3</td><td>8.3</td><td>11.1</td><td>20.1</td><td>68.8</td><td>6.6</td><td>20.5</td><td>72.9</td><td>5.4</td><td>19.5</td><td>75.1</td><td>5.6</td><td>18.4</td><td>76.0</td></tr><tr><td>30</td><td>80</td><td>12</td><td>100</td><td>29.3</td><td>19.0</td><td>11.4</td><td>5.0</td><td>12.4</td><td>19.0</td><td>68.6</td><td>5.8</td><td>20.6</td><td>73.6</td><td>4.8</td><td>18.7</td><td>76.6</td><td>5.4</td><td>17.4</td><td>77.3</td></tr><tr><td>31</td><td>80</td><td>14</td><td>100</td><td>26.0</td><td>12.4</td><td>8.0</td><td>1.0</td><td>13.2</td><td>18.6</td><td>68.2</td><td>4.9</td><td>20.8</td><td>74.3</td><td>4.4</td><td>18.5</td><td>77.1</td><td>5.7</td><td>17.5</td><td>76.8</td></tr><tr><td>32</td><td>80</td><td>16</td><td>100</td><td>23.5</td><td>9.6</td><td>2.0</td><td>0.0</td><td>14.4</td><td>17.9</td><td>67.7</td><td>4.4</td><td>21.7</td><td>73.9</td><td>4.5</td><td>19.0</td><td>76.6</td><td>5.9</td><td>16.2</td><td>77.9</td></tr><tr><td>33</td><td>80</td><td>18</td><td>100</td><td>19.7</td><td>2.6</td><td>0.0</td><td>0.0</td><td>15.4</td><td>17.4</td><td>67.2</td><td>3.6</td><td>23.0</td><td>73.4</td><td>4.5</td><td>19.2</td><td>76.3</td><td>6.2</td><td>16.3</td><td>77.5</td></tr><tr><td>34</td><td>80</td><td>20</td><td>100</td><td>16.7</td><td>1.0</td><td>0.0</td><td>0.0</td><td>16.3</td><td>16.5</td><td>67.2</td><td>3.4</td><td>23.4</td><td>73.2</td><td>4.6</td><td>18.8</td><td>76.6</td><td>6.3</td><td>15.6</td><td>78.1</td></tr><tr><td>35</td><td>80</td><td>22</td><td>100</td><td>14.5</td><td>0.0</td><td>0.0</td><td>0.0</td><td>17.8</td><td>15.2</td><td>67.1</td><td>3.4</td><td>23.7</td><td>72.9</td><td>4.7</td><td>18.4</td><td>76.9</td><td>6.6</td><td>15.3</td><td>78.1</td></tr><tr><td>36</td><td>80</td><td>24</td><td>100</td><td>11.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>19.0</td><td>13.7</td><td>67.3</td><td>3.5</td><td>23.5</td><td>73.1</td><td>5.0</td><td>18.0</td><td>77.0</td><td>6.7</td><td>14.8</td><td>78.4</td></tr><tr><td>Average:</td><td></td><td></td><td></td><td>23.2</td><td>13.8</td><td>9.4</td><td>6.4</td><td>12.7</td><td>16.6</td><td>70.7</td><td>5.3</td><td>19.3</td><td>75.4</td><td>6.0</td><td>17.2</td><td>76.7</td><td>7.0</td><td>15.4</td><td>77.6</td></tr></tbody></table>

cinating the elders as many as possible. The rest of the results in columns “CVG (%)” show that Pfizer (Pf) is always higher than Moderna (MO) and AstraZeneca (AZ) vaccines. The main reason is the higher quantity of Pfizer vaccines purchased and stored by the French government. These results of Table 6 are helpful to extract some rules for the government when designing a vaccine distribution network. For instance, when planning to vaccinate eligible individuals through a set of $|I| = 80$ vaccination centers for a horizon of $|T| = 24$ weeks, the optimal vaccination strategy would be opening the whole distribution centers and allocating, on average, 74% (average over 67.3, 73.1, 77.0, and 78.4), 17.5% (average over 13.7, 23.5, 18.0, and 14.8), and 8.5% (average over 19.0, 3.5, 5.0, and 6.7) of the available capacity (purchased or stored) of vaccines Pfizer, Moderna, and AstraZeneca vaccines to vaccinate the population.

## 7. Sensitivity analyses and managerial insights

### 7.1. Sensitivity analyses

This section provides a comprehensive sensitivity analysis of three main criteria, including the total cost of the system, the total inventory of different vaccines, and the total unmet demand of each age-based class with respect to changes of certain input parameters, including the minimum percentage of demand to be met ($\beta$), the maximum supply of vaccines at vaccination centers ($\bar{Q}^{VC}$), inventory holding capacity at vaccination centers ($Q^{VC}$), inventory holding cost at vaccination centers ($h^{VC}$), the demand of vaccines

($d$), the time interval between two doses injection ($L$), fairness service level gap ($\delta$), and unmet cost ($cu$). Figures 6–13 illustrate these sensitivities, wherein the value “1.00” on x-axis represents the original level of the corresponding parameter in the case study and other values on x-axis are proportional to the original value.

Figure 6 illustrates how altering parameter $\beta$ affects the three mentioned criteria. Precisely, Fig. 6(a) shows that increasing $\beta$ forces the system to vaccinate more individuals, and consequently, the corresponding costs also increase linearly. Figure 6(b) illustrates that the level of inventory for all three vaccines decreases since a larger demand should be fulfilled, and consequently, a larger portion of the inventory should be spent. Since a higher amount of inventory belongs to the Pfizer vaccine, its decrease happens with a higher slope. Finally, Fig. 6(c) illustrates that imposing the system to fulfill more demands causes more unmet demand for all age-based classes. The reason goes back to the limited availability of vaccines in the distribution system. Indeed, the system’s capacity to fulfill extra demand is limited, and the extra demand remains unmet.

Figure 7 investigates whether changing the supply capacity of vaccines in vaccination centers impacts the concerned three criteria. Figure 7(a) shows that the total cost of the system decreases as more vaccines are allowed to be supplied in the vaccination center. The potential reason for this decrease is twofold, the reduction in the total inventory holding cost or the reduction in the unmet demand cost. As Fig. 7(b) illustrates, the level of inventory increases, while Fig. 7(c) shows that the unmet demand decreases.
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Impact of minimum percentage of demand to be met on total cost, inventory, and unmet demand.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Impact of the maximum supply of vaccines on total cost, inventory, and unmet demand.</figcaption></figure>

What can be inferred is that the reduction of the unmet demand cost outweighs the increase of the inventory holding cost. Consequently, the total cost of the system decreases. The trends in Fig. 7(b) and (c) were expected since any increase in the supply capacity of the vaccines permits to store of more vaccines (i.e., an increase of the inventory level), and much larger storage helps to vaccinate more individuals (i.e., less unmet demand). What can be extracted more from Fig. 7(c) is that the unmet demand of age-based classes with lower priority (i.e., "18–49") always remains higher than the unmet demand of high-priority classes (i.e., "> 75").

Moreover, Fig. 8 investigates whether any change in the storage capacity of vaccination centers affects the performance of the distribution system. As can be seen, by increasing the storage capacity, all three concerned criteria show a decreasing trend.

Importantly, significant decreases (e.g., 80% for age-based class "18–49") happen in the unmet demand for all age-based classes

when the capacity of the vaccination centers increases up to 40% (from 0.6 to 1.0 on x-axis of Fig. 8(c)). Indeed, the increase of the capacity could be either increasing the capacity of the current vaccination centers by adding more vaccination lines or even opening new vaccination centers. It can also be observed that increasing the current capacity of the vaccination campaign in France up to 80% (point 1.8 on x-axis) will cause important reductions in the unmet demand in all age-based classes; however, any increase higher than 80% will not affect the unmet demands. Similar to Fig. 7(c), the unmet demand of age-based class "18–49" always remains higher than other classes; however, a significant decrease happens on the unmet of this class when increasing the capacity.

Since the vaccines require special modes of storage (i.e., ultra-cold freezers/refrigerators), the inventory holding cost plays a vital role in the distribution centers' performance. In this regard, Fig. 9 illustrates how the system reacts to the increase in the inventory holding cost. As can be seen in Fig. 9(a), the total cost
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Impact of inventory holding capacity on total cost, inventory, and unmet demand.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Impact of inventory holding cost on total cost, inventory, and unmet demand.</figcaption></figure>

of the system increases with the increase of the holding cost. When this parameter increases, the reaction of the system is to absorb such an increase by decreasing the inventory level (Fig. 9(b)), with the hope to reduce the total inventory holding cost. On the other hand, any reduction in the inventory level of vaccines signifies less capability of the system to vaccinate the population. Accordingly, the unmet demand increases as illustrated in Fig. 9(c). In such a situation, the increase of the total unmet demand cost prevails to the reduction of the total holding inventory cost; hence, the system’s total cost increases. In addition, Fig. 9(b) reveals that by increasing the inventory holding cost, the system puts its effort into reducing the share of Pfizer and Moderna since they require ultra-cold freezers and possess a higher holding cost. However, the inventory level of the AstraZeneca vaccine remains stable since it costs less in terms of inventory holding.

Figure 10 depicts the impact of demand alteration on the performance of the system and illustrates how the increase in demand affects the cost of the system as well as the levels of inventory and unmet demand. Figure 10(a) represents the reaction of the objective function to the increase in demand. As can be seen, this reaction consists of two phases: a decreasing trend followed by an increasing trend in the total cost of the system. In the decreasing phase, the demand for vaccines is small, and it is mostly fulfilled via the initial inventory of the vaccines in both distribution and vaccination centers. During the time that the system consumes the initial inventory in this phase, no extra vaccines are neither distributed nor stored in the network. Hence, both inventory holding and transportation costs are reduced, leading to a reduction in the system’s total cost. On the other hand, once the system runs out of the initial inventory, the system starts distributing and storing new vaccines. This phenomenon then increases the total cost of the sys-
<figure><img src="image_13.png" /><figcaption>Fig. 10. Impact of demand on total cost, inventory, and unmet demand.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Impact of the time interval between two doses injection on total cost, inventory, and unmet demand.</figcaption></figure>

In Fig. 10(b), the decrease in the inventory level of Pfizer and Moderna is higher than AstraZeneca, since a more significant initial inventory has been considered for the two former vaccines. Finally, Fig. 10(c) shows an expected trend for the unmet demand, as any increase in the demand for vaccines increases the potential unmet demand due to the limited availability of the vaccines in the system.

A more interesting result has been obtained when investigating the impact of the time interval between two doses of injection of the vaccines on the performance of the system, as depicted in Fig. 11. Figure 11(a) illustrates that if the time interval is reduced to half of their recommended delay, the unmet demand increases by a factor of 2.5 as shown in Fig. 11(c). Indeed, when the system focuses on a complete vaccination with a short time interval, more individuals are completely vaccinated, but a considerable part of the population remains un-vaccinated (i.e., unmet demand).

Furthermore, if the recommended time interval increases by a factor of 1.5, the unmet demand decreases up to 50%. Increasing

the time interval signifies that more individuals receive at least one dose of vaccines. On the other hand, this increase will not always reduce the unmet demand since, in further periods, the accumulated demand for the second dose of vaccines in early periods is summed up with new demands for the first dose; hence, augmented demand for vaccines is imposed on the system. Therefore, the unmet demand starts to be increasing when the system is disabled to absorb such an augmented increase in the demand. Accordingly, Fig. 11(a) and (b), collectively, show that the optimal performance of the system happens when delaying the time interval between two doses of injection of vaccines with a factor of 1.5.

Figure 12 depicts the impact of the fairness service level on the unmet demand. Lower values of this service level guarantee fair access for the population among different regions to vaccines. On the contrary, higher values of this service level allow a higher dispersion of vaccine allocation among different regions. The former situation should provide a higher equity and higher satisfaction among the whole population; however, the latter puts more effort into al-
<figure><img src="image_15.png" /><figcaption>Fig. 12. Impact of fairness service level gap on the unmet demand.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Impact of unmet cost on the objective function and the amount of unmet demand.</figcaption></figure>

locating vaccines to high-populated regions with a higher risk of infection due to higher social contacts. Figure 12 illustrates that increasing fairness (i.e., decreasing $δ$) is not always the best strategy in vaccination campaigns to stop pandemics. In fact, an optimum level of fairness leads to the minimum level of unmet demand. More explanations have been provided through managerial insights in Section 7.2.

In another sensitivity analysis, we investigate the impact of unmet demand cost $cu$ on the objective function (12) (Fig. 13(a)), the total distribution cost (i.e., the objective function (12)), except the last term), and the total amount of unmet demand (i.e., the last part of the objective function (12)). As can be seen from Fig. 13(a), increasing the unit of unmet demand cost increases the objective function, since the network aims to vaccinate as many individuals as possible to alleviate the impact of unmet demand. Figure 13(b) shows, in detail, how increasing the unit unmet demand cost increases (decreases) the total distribution cost (amount of unmet demand).

## 7.2. Discussion and managerial insights

In this section, some of the more interesting results are discussed, and a sort of managerial insights are provided.

We observed that if the recommended time interval increases by a specific factor, the unmet demand decreases significantly. It means that vaccinating a higher number of individuals, even with

a single dose, remarkably decreases the death rate. However, it implies vaccines with a high level of efficacy for their first-dose injection. Therefore, the availability of vaccines also plays a vital role in piloting a vaccination campaign. In reality, depending on the availability of vaccines, several countries have recommended an increase in time intervals between two injections (ECD, 2021). In this regard, by June 2021, 16 EU countries increased the time interval to provide more individuals with their first vaccination dose. Romania applied for such extension only in particular cases Dascalu et al. (2021); however, seven other EU countries (Iceland, Latvia, Lithuania, Malta, Slovakia, Slovenia, and Spain) did not apply such a strategy. Shortly after, seven countries removed this extension in September 2021 (Austria, Belgium, Czechia, Finland, Ireland, Portugal, and Spain). The possibility of such an increase in the interval between two injections introduced the concept of fractional dosing. Actually, due to the limited availability of vaccines, especially at the beginning of the pandemic, fractional dosing could be a crucial decision to make a trade-off between reducing the number of potential deaths (by vaccinating individuals with complete doses) or decelerating the pandemic (by providing a higher number of individuals with the first dose).

Scientists have always debated the impact of equity and fairness in pandemics (Enayati & Özaltın, 2020; Mohammadi et al., 2022). Although a higher level of fairness in vaccine allocation provides more social satisfaction among the population, the fair distribution of vaccines is not always the best barrier strategy to stop
pandemics. In this paper, we observed that the unmet demand is a convex function of fairness with a global minimum. Therefore, low and high levels of fairness do not represent the optimum allocation strategy. In fact, a low level of fairness puts more focus on the vaccination of populated areas with a higher risk of infection while ignoring less populated areas where the virus may transmit without sufficient barriers (e.g., vaccination). On the contrary, a high level of fairness guarantees fair access to the population to vaccines. However, there would be an imbalance between high and low-populated areas regarding barrier actions. Therefore, there should always be a trade-off between efforts for fair access to vaccination and attempts to control/stop the pandemic in high-risk (populated) areas. The focus on populated areas reveals the importance of the level of social contact between individuals during the progress of the pandemic. Therefore, from an optimization viewpoint, the effort of health policymakers should not only be vaccinating high-risk individuals to decrease deaths but also vaccinating the population with higher social contacts (i.e., populated areas). Although the former may directly decrease the number of deaths among the population, the latter directly breaks/decreases the circulation of the virus in society and indirectly reduces the number of deaths among the population. In this regard, one way to control social contact could be imposing quarantine restrictions in some regions while vaccinating others.

An important practicality of the proposed model for health policymakers can be highlighted from the results shown in Fig. 13(b). In fact, parameter *cu* takes the role of an interplay between two viewpoints, cost-oriented or mortality-oriented vaccination strategies. As a matter of fact, increasing this parameter drives the proposed model from a cost-oriented network toward a more mortality-oriented network. Accordingly, although the objective function of the proposed model is apparently a business objective, adjusting the unit unmet demand cost can satisfy mortality-related concerns. In this regard, a possibility could be separating the two parts of the objective function (i.e., distribution cost and amount of unmet demand) and optimizing them through a bi-objective model. The gain of such a model could be generating a set of non-dominated solutions and providing health policymakers with a portfolio of different vaccination strategies.

## 8. Conclusions

The most effective intervention strategy to control the COVID-19 pandemic, which requires careful preparation due to its availability, is vaccination, the principal hope of today's communities. In this regard, we proposed a novel exhaustive scheme that contains three phases of demand regulating, location-inventory-routing planning, and solution approach, wherein a broad range of concerns of the COVID-19 vaccine procurement and distribution have been considered. In the first stage, we proposed a modified SIR epidemiological model for handling the dynamicity of the vaccine demands, where the optimal control problem was involved in determining the number of vaccines that need to be purchased to control the pandemic situation. In the second stage, we have presented and formulated a multi-period, multi-vaccine, multi-depot location-inventory-routing problem in the form of a mixed-integer programming model, where various realistic concerns of the COVID-19 vaccine distribution, ranging from multi-dose vaccine regimen to prioritizing age groups, have been reflected. Finally, we offered an accelerated Benders decomposition algorithm to solve sizable test problems in the third stage. Although the primary aim of the current research was to cope with the concerns of the COVID-19 vaccine procurement and distribution, the soundness outcomes revealed that it could indeed be applied to cope with various vaccine distribution problems. In fact, due to the complexity of the considerations that had to be reflected

in the distribution planning of the COVID-19 vaccines, the proposed optimization model factored in almost all essential concerns in the vaccine distribution. As a result, it could be fruitful to distribute existing or forthcoming vaccines. The computational results demonstrated that the proposed framework could provide a practical plan for the preparation and distribution of vaccines. In this regard, we provided operational planning to prioritize vaccination target groups, and we were also able to schedule a time interval between injections of different vaccines if possible to control the pandemic better.

For future research, we can mention a number of interesting cases that can be considered in both levels of the presented framework. Regarding the first level of the framework (i.e., the adjusted SIR model), considering the efficacy of various vaccines would be interesting to be considered in the model. Various facets of this problem make it worthwhile to study. If immunizations are only moderately effective, the booster dose should be administered sooner. Furthermore, if the vaccine's efficiency is excellent, the option of delaying the second dosage of the vaccine rises, allowing a more significant percentage of the target population to get vaccinated. The second involves anti-vaccination individuals. These people increase the possibility of infection and raise mistrust among others who have not yet been vaccinated. Regarding the second phase of the framework (i.e., vaccine distribution network), the proposed optimization model can be integrated with the adjusted SIR model in an interactive way to develop a myopic model. In this model, the SIR model and the optimization model interact in each period, and the vaccination rate *u* in Section 1 is considered as a decision variable, and the vaccination decisions in each period directly affect the estimated demand of the further periods.

## Appendix A. Proof of Theorem 1

**Proof.** The optimal control intervention exists when the following conditions are satisfied:

* The solution space of system dynamics (1) with control variable $u \cap n \neq \emptyset$.

* The mentioned set $n$ is closed, and convex, and the state system is represented with a linear function of the control variable where coefficients depend on time and also on state variables.

* Integrand of (8): $L$ is convex on $n$ and $L(S, I, Q, UQ, R, u) \geq f(u)$, where $f(u)$ is continuous, and $\|u\|^{-1}f(u) \to \infty$, when $\|u\| \to \infty$. Indeed, $\|\cdot\|$ indicates the $L(0, T_f)$ norm.

From the system dynamics (1), the total population $N = S + I + Q + UQ + R$. Hence, $\frac{dN}{dt} = \Lambda - dN - \delta_1I - \delta_2Q - \delta_3UQ \leq \Lambda - dN$, and $0 < N(t) \leq N(0)e^{-dt} + \frac{\Lambda}{d}(1 - e^{-dt})$, where $N(0) = S(0) + I(0) + Q(0) + UQ(0) + R(0)$. Also, when $t \to \infty$, we have $0 < N(t) \leq \frac{\Lambda}{d}$.

For the control variable in $n$, the solution of the system dynamics (1) is bounded and right-hand side functions are also locally Lipschitzian. The Picard Lindelof theorem demonstrates that the first condition is met (Coddington & Levinson, 1955). The control set $n$ is closed and convex by definition. Again, all the equations of the system dynamics (1) are written as linear equations in $u$ where state variables depend on the coefficients, and hence the second condition is also satisfied. Moreover, the quadratic nature of the control variable guarantees the convex property of integrand $L(S, I, Q, UQ, R, u)$. Let $f(u) = w_3u_2$, next $L(S, I, Q, UQ, R, u) \geq f(u)$. Here, $f$ is continuous and $\|u\|^{-1}f(u) \to \infty$, when $\|u\| \to \infty$. Hence, the third condition is also met. Therefore, it is concluded that there is the control variable $u^*$ with the condition $J(u^*) = \min[J(u)]$ (Fleming et al., 2012; Gaff & Schaefer, 2009). □
## Appendix B. Pontryagin's maximum principle (Pontryagin, 1987)

By Pontryagin's Maximum Principle, we have derived here the necessary conditions for optimal control functions for the optimal control problem (8) and (9) (Fleming et al., 2012; Pontryagin, 1987). Let us define the Hamiltonian function as (B.1) and (B.2):

$$ \bar{H}(S, I, Q, UQ, R, u, \lambda) = L(S, I, Q, UQ, R, u) + \lambda_1 \frac{dS}{dt} + \lambda_2 \frac{dI}{dt} \\ + \lambda_3 \frac{dQ}{dt} + \lambda_4 \frac{dUQ}{dt} + \lambda_5 \frac{dR}{dt} \quad (B.1) $$

Hence,

$$ \begin{aligned} \bar{H} = & w_1 I(t) + w_2 Q(t) + w_3 UQ(t) \\ & + \lambda_1 [\Lambda - (\beta_1 I + \beta_2 Q + \beta_3 UQ)S - dS - u(t)S] \\ & + \lambda_2 [(\beta_1 I + \beta_2 Q + \beta_3 UQ)S - (d + \delta_1 + k)I] \\ & + \lambda_3 [k\tau I - (d + \delta_2)Q - v_1 Q] \\ & + \lambda_4 [k(1 - \tau)I - (d + \delta_3)UQ - v_2 UQ] \\ & + \lambda_5 [v_1 Q + v_2 UQ - dR + u(t)S] \end{aligned} \quad (B.2) $$

## Appendix C. Proof of Theorem 2

**Proof.** Let $u^*$ be optimal control variable and $S^*, I^*, Q^*, UQ^*, R^*$ be optimal state variables of system dynamics (1) that minimize the objective function (8). With the help of Pontryagin's Maximum Principle and the defined adjoint variables $\lambda_i$, $i = 1, \dots, 5$, we have $\frac{\partial \lambda_1}{\partial t} = -\frac{\partial \bar{H}}{\partial S}$, $\frac{\partial \lambda_2}{\partial t} = -\frac{\partial \bar{H}}{\partial I}$, $\frac{\partial \lambda_3}{\partial t} = -\frac{\partial \bar{H}}{\partial Q}$, $\frac{\partial \lambda_4}{\partial t} = -\frac{\partial \bar{H}}{\partial UQ}$, $\frac{\partial \lambda_5}{\partial t} = -\frac{\partial \bar{H}}{\partial R}$, which lead to obtain Eq. (10). In regard to the transversality conditions $\lambda_i(T_f) = 0$, $i = 1, \dots, 5$ and optimality conditions, we have:

$$ \begin{aligned} \frac{\partial \bar{H}}{\partial u} |_{u=u^*} &= 0 \\ \frac{\partial \bar{H}}{\partial u} |_{u=u^*} &= 2w_3 u^* - \lambda_1 S + \lambda_5 S = 0 \\ u^* &= \frac{(\lambda_1 - \lambda_5)S^*}{2w_3} \end{aligned} \quad (C.1) $$

Now from these findings along with the characteristics of control set $\square$, we have:

$$ u^* = \begin{cases} 0 & \text{if } \frac{(\lambda_1 - \lambda_5)S^*}{2w_3} < 0 \\ \frac{(\lambda_1 - \lambda_5)S^*}{2w_3} & \text{if } 0 \le \frac{(\lambda_1 - \lambda_5)S^*}{2w_3} \le 1 \\ 1 & \text{if } \frac{(\lambda_1 - \lambda_5)S^*}{2w_3} > 1 \end{cases} \quad (C.2) $$

which is equivalent to (11). □

## Appendix D. Linearization

With the help of auxiliary variables $TY_{pjit}$ and $TZ_{pjivt}$ ($TY_{pjit}$, $TZ_{pjivt} \ge 0$), these two multiplications are linearized by the set of constraints (D.1)-(D.3) and constraints (D.4)-(D.6), respectively.

$$ TY_{pjit} \le T_{pjit} \quad \forall p \in P, j \in J, i \in I, t \in T \quad (D.1) $$

$$ TY_{pjit} \le MY_{jit} \quad \forall p \in P, j \in J, i \in I, t \in T \quad (D.2) $$

$$ TY_{pjit} \ge T_{pjit} - M(1 - Y_{jit}) \quad \forall p \in P, j \in J, i \in I, t \in T \quad (D.3) $$

$$ TZ_{pjivt} \le T_{pjit} \quad \forall p \in P, j \in J, i \in I, v \in V, t \in T \quad (D.4) $$

$$ TZ_{pjivt} \le MZK_{jvt} \quad \forall p \in P, j \in J, i \in I, v \in V, t \in T \quad (D.5) $$

$$ TZ_{pjivt} \ge T_{pjit} - M(1 - ZK_{jvt}) \quad \forall p \in P, j \in J, i \in I, v \in V, t \in T \quad (D.6) $$

Finally, the constraints (13), (14), and (22) are replaced by the following constraints (13), (14), and (22), respectively.

$$ \sum_{p \in P} (1 - \gamma_{pjt}) U_{pjt} \ge \sum_{p \in P} \sum_{i \in I} TY_{pjit} \quad \forall j \in J, t \in T \quad (D.7) $$

$$ \sum_{k \in K} G_{pkit} \le \sum_{j \in J} (1 - \mu_{pjit}) TY_{pjit} \quad \forall p \in P, i \in I, t \in T \quad (D.8) $$

$$ \sum_{i \in I} \sum_{p \in P} TZ_{pjivt} \le Q_v^V \quad \forall j \in J, v \in V, t \in T \quad (D.9) $$

## Appendix E. The MVMDLIRP Benders reformulation

Let $\bar{Z}_{jt}$, $\bar{Y}_{jit}$, $\bar{Z}_{Kjvt}$, and $\bar{X}_{jivt}$ be the vectors of fixed $Z_{jt}$, $Y_{jit}$, $ZK_{jvt}$, and $X_{jivt}$ variables, respectively. The Benders sub-problem can be written as follows:

$$ \min Z = \sum_{p \in P} \sum_{j \in J} c_{pj} U_{pjt} + \sum_{p \in P} \sum_{j \in J} \sum_{i \in I} \sum_{t \in T} c_{spjt} T_{pjit} + \sum_{p \in P} \sum_{j \in J} \sum_{t \in T} h_{pjit}^{DC} I_{pjit} \\ + \sum_{p \in P} \sum_{i \in I} \sum_{t \in T} h_{pi}^{VC} I'_{pkit} + cu \sum_{p \in P} \sum_{k \in K} \sum_{i \in I} \sum_{t \in T} w_{pkit} \left( \frac{R_{pkit}}{d_{pkit}} \right) \quad (E.1) $$

s.t.: Constraints (15)-(21), (24), (25), (D.1), (D.4), (D.7)-(D.9) and,

$$ I_{pjt}^{DC} \le Q_{pj}^{DC} \bar{Z}_{jt} \quad \forall p \in P, j \in J, t \in T \quad (E.2) $$

$$ \frac{\sum_{p \in P} \sum_{k \in K} G_{pkit}}{\sum_{p \in P} \sum_{k \in K} d_{pkit}} \le \frac{\sum_{p \in P} \sum_{k \in K} G_{pkit't}}{\sum_{p \in P} \sum_{k \in K} d_{pkit't}} + (2 - \bar{Y}_{jit} - \bar{Y}_{ji't}) \quad \forall i, i' \in I; i \ne i', j \in J, t \in T \quad (E.3) $$

$$ T_{pjit} \le M\bar{Y}_{jit} \quad \forall p \in P, j \in J, i \in I, t \in T \quad (E.4) $$

$$ UU_{jvt} - UU_{ivt} + |I|\bar{X}_{jivt} \le |I| - 1 \quad \forall i, j \in I, v \in V, t \in T \quad (E.5) $$

$$ TY_{pjit} \le M\bar{Y}_{jit} \quad \forall p \in P, j \in J, i \in I, t \in T \quad (E.6) $$
$$TY_{pjit} \geq T_{pjit} - \mathcal{M}(1 - \bar{Y}_{jit}) \quad \forall p \in P, j \in J, i \in I, t \in T \quad (E.7)$$

$$TZ_{pjivt} \leq \mathcal{M}\bar{Z}_{jvt} \quad \forall p \in P, j \in J, i \in I, v \in V, t \in T \quad (E.8)$$

$$TZ_{pjivt} \geq T_{pjit} - \mathcal{M}(1 - \bar{Z}_{jvt}) \quad \forall p \in P, j \in J, i \in I, v \in V, t \in T \quad (E.9)$$

Let $\mathcal{A} = \{A_{pt} \geq 0 | p \in P, t \in T\}$, $\mathcal{B} = \{B_{jt} \geq 0 | j \in J, t \in T\}$, $\mathcal{C} = \{C_{pit} \geq 0 | p \in P, i \in I, t \in T\}$, $\mathcal{D} = \{D_{pjit} \geq 0 | p \in P, i \in I, j \in J, t \in T\}$, $\mathcal{E} = \{\mathcal{E}_{pjit} \geq 0 | p \in P, i \in I, j \in J, t \in T\}$, $\mathcal{F} = \{\mathcal{F}_{pjit} \geq 0 | p \in P, i \in I, j \in J, t \in T\}$, $\mathcal{G} = \{G_{pjit} \in \mathbb{R} | p \in P, j \in J, t \in T\}$, $\mathcal{H} = \{\mathcal{H}_{pjit} \geq 0 | p \in P, j \in J, t \in T\}$, $\mathcal{J} = \{\mathcal{J}_{pkit} \geq 0 | p \in P, i \in I, k \in K, t \in T; t \leq L_p\}$, $\mathcal{K} = \{\mathcal{K}_{pkit} \geq 0 | p \in P, i \in I, k \in K, t \in T; t \leq L\}$, $\mathcal{L} = \{\mathcal{L}_{pkit} \geq 0 | p \in P, i \in I, k \in K, t \in T; t > L_p\}$, $\mathcal{N} = \{\mathcal{N}_{ii't} \geq 0 | i, i' \in I, t \in T; i! = i'\}$, $\mathcal{O} = \{\mathcal{O}_{ii'jt} \geq 0 | i, i' \in I, j \in J, t \in T; i! = i'\}$, $\mathcal{P} = \{\mathcal{P}_{pit} \geq 0 | p \in P, i \in I, t \in T\}$, $\mathcal{Q} = \{\mathcal{Q}_{pit} \in \mathbb{R} | p \in P, i \in I, t \in T\}$, $\mathcal{R} = \{\mathcal{R}_{pit} \geq 0 | p \in P, i \in I, t \in T\}$, $\mathcal{S} = \{S_{pki} \in \mathbb{R} | p \in P, i \in I, k \in K\}$, $\mathcal{T} = \{T_{pkit} \in \mathbb{R} | p \in P, i \in I, k \in K, t \in T\}$, $\mathcal{U} = \{U_{pjit} \geq 0 | j \in J, v \in V, t \in T\}$, $\mathcal{V} = \{V_{pjit} \geq 0 | p \in P, i \in I, j \in J, v \in V, t \in T\}$, $\mathcal{X} = \{X_{pjivt} \geq 0 | p \in P, i \in I, j \in J, v \in V, t \in T\}$, $\mathcal{Y} = \{Y_{pjivt} \geq 0 | i \in I, j \in J, v \in V, t \in T\}$, $\mathcal{Z} = \{Z_{pjit} \geq 0 | p \in P, i \in I, j \in J, t \in T\}$ be the vectors of the dual variables associated with the Benders sub-problem constraints. The Benders master problem can be written as follows:

$$\min Z = \sum_{j \in J} \sum_{t \in T} f_j Z_{jt} + \sum_{j \in J} \sum_{i \in M} \sum_{v \in V} \sum_{t \in T} \pi_v \Delta_{ij} X_{pjivt} + \zeta \quad (E.10)$$

s.t.: Constraints (27)–(40), (44–47), and

$$\begin{aligned} & - \sum_{p \in P} \sum_{t \in T} Q_p^{CH} A_{pt} - \sum_{p \in P} \sum_{i \in I} \sum_{j \in J} \sum_{t \in T} \mathcal{M}_{jit} Y_{jit} \mathcal{E}_{pjit} \\ & - \sum_{p \in P} \sum_{i \in I} \sum_{j \in J} \sum_{t \in T} \mathcal{M}_{jit} (1 - Y_{jit}) \mathcal{F}_{pjit} \\ & - \sum_{p \in P} \sum_{j \in J} \sum_{t \in T} Q_{pj}^{DC} Z_{jt} \mathcal{H}_{pjit} + \sum_{p \in P} \sum_{i \in I} \sum_{k \in K} \sum_{t \in T} \sum_{t \leq L_p} \beta_k d_{pkit} \mathcal{J}_{pkit} \\ & - \sum_{p \in P} \sum_{i \in I} \sum_{k \in K} \sum_{t \in T} d_{pkit} \mathcal{K}_{pkit} \end{aligned}$$

$$\begin{aligned} & + \sum_{p \in P} \sum_{i \in I} \sum_{k \in K} \sum_{t > L_p} \beta_k d_{pkit} \mathcal{L}_{pkit} - \sum_{p \in P} \sum_{i \in I} \sum_{k \in K} \sum_{t > L_p} d_{pkit} \mathcal{M}_{pkit} \\ & - \sum_{i \in I} \sum_{i' \in I} \sum_{j \in J} \sum_{t \in T} (2 - Y_{jit} - Y_{ji't}) \mathcal{O}_{ii'jt} \\ & - \sum_{p \in P} \sum_{i \in I} \sum_{t \in T} Q_{pi}^{VC} \mathcal{P}_{pit} - \sum_{p \in P} \sum_{i \in I} \sum_{t \in T} \bar{Q}_{pi}^{VC} \mathcal{R}_{pit} \\ & + \sum_{p \in P} \sum_{i \in I} \sum_{k \in K} d_{pki,1} S_{pki} + \sum_{p \in P} \sum_{i \in I} \sum_{k \in K} \sum_{t \in T \setminus \{1\}} d_{pkit} \mathcal{T}_{pkit} \\ & - \sum_{j \in J} \sum_{v \in V} \sum_{t \in T} Q_v^{VC} U_{jvt} - \sum_{p \in P} \sum_{i \in I} \sum_{j \in J} \sum_{v \in V} \sum_{t \in T} \mathcal{M}_{jvvt} ZK_{jvvt} W_{pjivt} \\ & - \sum_{p \in P} \sum_{i \in I} \sum_{j \in J} \sum_{v \in V} \sum_{t \in T} \mathcal{M}_{jvvt} (1 - ZK_{jvvt}) X_{pjivt} \\ & + \sum_{i \in I} \sum_{j \in J} \sum_{v \in V} \sum_{t \in T} (-|I| + 1 + |X_{pjivt}|) Y_{ijvt} \\ & - \sum_{p \in P} \sum_{i \in I} \sum_{j \in J} \sum_{t \in T} \mathcal{M}_{jit} Y_{jit} Z_{pjit} \leq \zeta \end{aligned} \quad (E.11)$$

$$\zeta \geq 0 \quad (E.12)$$

$\forall \mathcal{A}, \mathcal{B}, \mathcal{C}, \mathcal{D}, \mathcal{E}, \mathcal{F}, \mathcal{G}, \mathcal{H}, \mathcal{J}, \mathcal{K}, \mathcal{L}, \mathcal{M}, \mathcal{N}, \mathcal{O}, \mathcal{P}, \mathcal{Q}, \mathcal{R}, \mathcal{S}, \mathcal{T}, \mathcal{U}, \mathcal{V}, \mathcal{W}, \mathcal{X}, \mathcal{Y}, \mathcal{Z} \in P_{\nabla} \text{ where } \nabla \text{ signifies the polyhedron specified by the constraints of the problem, and } P_{\nabla} \text{ implicates the set of extreme points of } \nabla$. Constraint (E.11) represents the optimality cut which can be generated consecutively in accordance with Benders sub-problem solution and the vectors of the dual variables. In addition, $\mathcal{M}_{jit}$ and $\mathcal{M}_{jvvt}$ are obtained via Eqs. (E.13) and (E.14), respectively.

$$\mathcal{M}_{jit} = \min \{ Q_p^{CH}, Q_{pj}^{DC}, Q_{pi}^{VC}, \sum_{p \in P} \sum_{k \in K} \sum_{l=t}^{T} d_{pkil} \} \quad \forall i \in I, j \in J, t \in T \quad (E.13)$$

$$\mathcal{M}_{jvvt} = \min \{ Q_{pj}^{DC}, Q_v^{VC}, \sum_{p \in P} \sum_{k \in K} \sum_{i \in I} \sum_{l=t}^{T} d_{pkil} \} \quad \forall j \in J, v \in V, t \in T \quad (E.14)$$

## Appendix F. Properties of French metropolitan properties

<caption>Table F.1.</caption>

<table><caption>**Table F.1**

Properties per region in France (Insee.fr, 2022).</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">Population</th><th rowspan="2">Education (%)<sup>a</sup></th><th rowspan="2">GDP (€)<sup>b</sup></th><th rowspan="2">Poverty (%)<sup>c</sup></th></tr><tr><th>Density (ratio)</th><th>Male (%)</th><th>Female (%)</th></tr></thead><tbody><tr><td>Auvergne-Rhne-Alpes</td><td>1,92</td><td>48,64</td><td>51,36</td><td>79,65</td><td>33,205</td><td>12,7</td></tr><tr><td>Bourgogne-Franche-Comté</td><td>1,00</td><td>48,66</td><td>51,34</td><td>76,24</td><td>27,026</td><td>12,9</td></tr><tr><td>Bretagne</td><td>2,05</td><td>48,60</td><td>51,40</td><td>81,70</td><td>29,500</td><td>10,9</td></tr><tr><td>Centre-Val de Loire</td><td>1,12</td><td>48,50</td><td>51,50</td><td>76,25</td><td>27,859</td><td>13,1</td></tr><tr><td>Grand Est</td><td>1,64</td><td>48,75</td><td>51,25</td><td>76,78</td><td>27,529</td><td>14,8</td></tr><tr><td>Hauts-de-France</td><td>3,20</td><td>48,45</td><td>51,55</td><td>74,31</td><td>26,948</td><td>18</td></tr><tr><td>Ile-de-France</td><td>17,10</td><td>48,26</td><td>51,74</td><td>81,27</td><td>57,600</td><td>15,6</td></tr><tr><td>Normandie</td><td>1,88</td><td>48,36</td><td>51,64</td><td>74,62</td><td>27,168</td><td>13,5</td></tr><tr><td>Nouvelle-Aquitaine</td><td>1,19</td><td>48,07</td><td>51,93</td><td>79,05</td><td>28,467</td><td>13,6</td></tr><tr><td>Occitanie</td><td>1,36</td><td>48,34</td><td>51,66</td><td>79,61</td><td>28,157</td><td>17,2</td></tr><tr><td>Pays de la Loire</td><td>1,97</td><td>48,71</td><td>51,29</td><td>78,85</td><td>30,258</td><td>10,8</td></tr><tr><td>Provence-Alpes-Cte d'Azur</td><td>2,71</td><td>47,76</td><td>52,24</td><td>78,64</td><td>31,580</td><td>17,3</td></tr></tbody></table>

<footnote><sup>a</sup> Rate of total regional population with a high school education and higher.</footnote>

<footnote><sup>b</sup> Gross Domestic Product of each individual.</footnote>

<footnote><sup>c</sup> Rate of total regional population with a salary less than poverty threshold.</footnote>
## References

Abbasi, B., Fadaki, M., Kokshagina, O., Saeed, N., & Chhetri, P. (2020). Modeling vaccine allocations in the COVID-19 pandemic: A case study in Australia. Available at SSRN 3744520.

Abou-Ismail, A. (2020). Compartmental models of the COVID-19 pandemic for physicians and physician-scientists. *SN Comprehensive Clinical Medicine, 2*, 852–858.

Ahmadzadeh, E., & Vahdani, B. (2017). A location-inventory-pricing model in a closed loop supply chain network with correlated demands and shortages under a periodic review system. *Computers and Chemical Engineering, 101*, 148–166.

Alam, S. T., Ahmed, S., Ali, S. M., Sarker, S., Kabir, G., et al., (2021). Challenges to COVID-19 vaccine supply chain: Implications for sustainable development goals. *International Journal of Production Economics, 239*, 108193.

Alcaraz, G. G., & Vargas-De-León, C. (2012). Modeling control strategies for influenza a H1N1 epidemics: SIR models. *Revista Mexicana de Física, 58*(1), 37–43.

Alkaabneh, F., Diabat, A., & Gao, H. O. (2020). Benders decomposition for the inventory vehicle routing problem with perishable products and environmental costs. *Computers and Operations Research, 113*, 104751.

Anand, N., Sabarinath, A., Geetha, S., & Somanath, S. (2020). Predicting the spread of COVID-19 using SIR SIR model augmented to incorporate quarantine and testing. *Transactions of the Indian National Academy of Engineering, 5(2)*, 141–148.

Angeli, M., Neofotistos, G., Mattheakis, M., & Kaxiras, E. (2021). Modeling the effect of the vaccination campaign on the COVID-19 pandemic. *Chaos, Solitons and Fractals, 154*, 111621.

Arfan, M., Shah, K., Abdeljawad, T., Mlaiki, N., & Ullah, A. (2021). A Caputo power law model predicting the spread of the COVID-19 outbreak in Pakistan. *Alexandria Engineering Journal, 60*(1), 447–456.

Asano, E., Gross, L. J., Lenhart, S., & Real, L. A. (2008). Optimal control of vaccine distribution in a rabies metapopulation model. *Mathematical Biosciences and Engineering, 5*(2), 219.

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

Bertsimas, D., Digalakis Jr, V., Jacquillat, A., Li, M. L., & Previero, A. (2022). Where to locate COVID-19 mass vaccination facilities? *Naval Research Logistics (NRL), 69*(2), 179–200.

Besiou, M., & Van Wassenhove, L. N. (2021). System dynamics for humanitarian operations revisited. *Journal of Humanitarian Logistics and Supply Chain Management, 11*(4), 599–608.

Betcheva, L., Erhun, F., Feylessoufi, A., Fryers, P., Gonçalves, P., Jiang, H., Kattuman, P., Pape, T., Pari, A., Scholtes, S. et al. (2021). An adaptive research approach to COVID-19 forecasting for regional health systems in England. Available at SSRN 3695258.

Biswas, D., & Alfandari, L. (2022). Designing an optimal sequence of non-pharmaceutical interventions for controlling COVID-19. *European Journal of Operational Research, 303*(3), 1372–1391.

Boland, N., Fischetti, M., Monaci, M., & Savelsbergh, M. (2016). Proximity Benders: A decomposition heuristic for stochastic programs. *Journal of Heuristics, 22*(2), 181–198.

Burgos, R. M., Badowski, M. E., Drwiega, E., Ghassemi, S., Griffith, N., Herald, F., Johnson, M., Smith, R. O., & Michienzi, S. M. (2021). The race to a COVID-19 vaccine: Opportunities and challenges in development and distribution. *Drugs Context, 10* 2020-12-2. https://doi.org/10.7573/dic.2020-12-2. PMID: 33643421; PMCID: PMC7888064.

Büyüktahtın, I. E., des Bordes, E., & Kıbış, E. Y. (2018). A new epidemics–logistics model: Insights into controlling the Ebola virus disease in west Africa. *European Journal of Operational Research, 265*(3), 1046–1063.

Capelle, T., Cortés, C. E., Gendreau, M., Rey, P. A., & Rousseau, L.-M. (2019). A column generation approach for location-routing problems with pickup and delivery. *European Journal of Operational Research, 272*(1), 121–131.

Castilho, C. (2006). Optimal control of an epidemic through educational campaigns. *Electronic Journal of differential equations, 2006*(125), 1–11.

Çalik, H., Oulamara, A., Prodhon, C., & Salhi, S. (2021). The electric location-routing problem with heterogeneous fleet: Formulation and Benders decomposition approach. *Computers and Operations Research, 131*, 105251.

Chen, J., Hoops, S., Marathe, A., Mortveit, H., Lewis, B., Venkatramanan, S., Hadadan, A., Bhattacharya, P., Adiga, A., Vullikanti, A. et al. (2021). Prioritizing allocation of COVID-19 vaccines based on social contacts increases vaccination effectiveness. *MedRxiv*.

Chen, S.-I., Norman, B. A., Rajgopal, J., Assi, T. M., Lee, B. Y., & Brown, S. T. (2014). A planning model for the WHO-EPI vaccine distribution network in developing countries. *IIE Transactions, 46*(8), 853–865.

Chen, X., Li, M., Simchi-Levi, D., & Zhao, T. (2020). Allocation of COVID-19 vaccines under limited supply. Available at SSRN 3678986.

Coddington, E. A., & Levinson, N. (1955). *Theory of ordinary differential equations*. Tata McGraw-Hill Education.

Coelho, L. C., & Laporte, G. (2014). Improved solutions for inventory-routing problems through valid inequalities and input ordering. *International Journal of Production Economics, 155*, 391–397.

Cordeau, J.-F., Laganà, D., Musmanno, R., & Vocaturo, F. (2015). A decomposition-based heuristic for the multiple-product inventory-routing problem. *Computers and Operations Research, 55*, 153–166.

Cordeau, J.-F., Pasin, F., & Solomon, M. M. (2006). An integrated model for logistics network design. *Annals of Operations Research, 144*(1), 59–82.

Corréa, A. I., Langevin, A., & Rousseau, L.-M. (2007). Scheduling and routing of automated guided vehicles: A hybrid approach. *Computers and Operations Research, 34*(6), 1688–1707.

Darvish, M., Archetti, C., Coelho, L. C., & Speranza, M. G. (2019). Flexible two-echelon location routing problem. *European Journal of Operational Research, 277*(3), 1124–1136.

Dascalu, S., Geambasu, O., Covaciu, O., Chereches, R. M., Diaconu, G., Dumitra, G. G., ... Popovici, E. D. (2021). Prospects of COVID-19 vaccination in Romania: Challenges and potential solutions. *Frontiers in Public Health, 9*, 1–5.

De Boeck, K., Decouttere, C., & Vandaele, N. (2020). Vaccine distribution chains in low-and middle-income countries: A literature review. *Omega, 97*, 102097.

DeRoo, S. S., Pudalov, N. J., & Fu, L. Y. (2020). Planning for a COVID-19 vaccination program. *JAMA, 323*(24), 2458–2459.

Van den Driessche, P., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences, 180*(1-2), 29-48.

Duijzer, L. E., van Jaarsveld, W., & Dekker, R. (2018a). The benefits of combining early aspecific vaccination with later specific vaccination. *European Journal of Operational Research, 271*(2), 606–619.

Duijzer, L. E., van Jaarsveld, W. L., Wallinga, J., & Dekker, R. (2018b). Dose-optimimal vaccine allocation over multiple populations. *Production and Operations Management, 27*(1), 143–159.

Duijzer, L. E., van Jaarsveld, W., & Dekker, R. (2018c). Literature review: The vaccine supply chain. *European Journal of Operational Research, 268*(1), 174–192.

(2021). Overview of the implementation of COVID-19 vaccination strategies and deployment plans in the EU/EEA. https://www.ecdc.europa.eu/en/publications-data/overview-implementation-covid-19-vaccination-strategies-and-deployment-plans/.

Enayati, S., & Özaltın, O. Y. (2020). Optimal influenza vaccine distribution with equity. *European Journal of Operational Research, 283*(2), 714–725.

Fleming, W. H., & Rishel, R. W. (2012). Deterministic and stochastic optimal control (Vol. 1). Springer Science & Business Media.

Forni, G., & Mantovani, A. (2021). COVID-19 vaccines: Where we stand and challenges ahead. *Cell Death and Differentiation, 28*(2), 626–639.

Gaff, H., & Schaefer, E. (2009). Optimal control applied to vaccination and treatment strategies for various epidemiological models. *Mathematical Biosciences and Engineering, 6*(3), 469.

Gamchi, N. S., Torabi, S. A., & Jolai, F. (2021). A novel vehicle routing problem for vaccine distribution using SIR epidemic model. *OR Spectrum, 43*(1), 155–188.

Gilani, H., & Sahebi, H. (2022). A data-driven robust optimization model by cutting hyperplanes on vaccine access uncertainty in COVID-19 vaccine supply chain. *Omega, 110*, 102637.

Gmira, M., Gendreau, M., Lodi, A., & Potvin, J.-Y. (2021). Tabu search for the time-dependent vehicle routing problem with time windows on a road network. *European Journal of Operational Research, 288*(1), 129–140.

Gnanvi, J. E., Salako, K. V., Kotanmi, G. B., & Kakaï, R. G. (2021). On the reliability of predictions on COVID-19 dynamics: A systematic and critical review of modelling techniques. *Infectious Disease Modelling, 6*, 258–272.

Graham, B. S. (2020). Rapid COVID-19 vaccine development. *Science, 368*(6494), 945–946.

Haddad, M. N., Martinelli, R., Vidal, T., Martins, S., Ochi, L. S., Souza, M. J. F., & Hartl, R. (2018). Large neighborhood-based metaheuristic and branch-and-price for the pickup and delivery problem with split loads. *European Journal of Operational Research, 270*(3), 1014–1027.

Inseef.fr (2022). Institut national de la statistique et des Études Économiques. https://www.insee.fr/fr/acueil.

Ivanov, D. (2020). Predicting the impacts of epidemic outbreaks on global supply chains: A simulation-based analysis on the coronavirus outbreak (COVID-19/SARS-CoV-2) case. *Transportation Research Part E: Logistics and Transportation Review, 136*, 101922.

Johns Hopkins University & Medicine (2021). Coronavirus COVID-19 global cases by the center for systems science and engineering (CSSE) at Johns Hopkins University (JHU). https://coronavirus.jhu.edu/map.html.

Kai, D., Goldstein, G.-P., Morgunov, A., Nangalia, V., & Rotkirch, A. (2020). Universal masking is urgent in the COVID-19 pandemic: SEIR and agent based models, empirical validation, policy recommendations. *arXiv preprint arXiv:2004.13553*

Planal, E. H. (2020). Om forum–COVID-19 scratch models to support local decisions. *Manufacturing and Service Operations Management, 22*(4), 645–655.

Kushwaha, S., Bahl, S., Bagha, A. K., Parmar, K. S., Javaid, M., Haleem, A., & Singh, R. P. (2020). Significant applications of machine learning for COVID-19 pandemic. *Journal of Industrial Integration and Management, 5*(04), 453–479.

Laarbi, H., Abta, A., Rachik, M., Labriji, E., Bouyaghroumni, J., & Labriji, E. (2013). Stability analysis and optimal vaccination strategies for an SIR epidemic model with a nonlinear incidence rate. *International Journal of Nonlinear science, 16*(4), 323–333.

Larissa, P. G. P., Zabinsky, Z. B., Zameer, M., Chu, Y., Muteia, M. M., Resende, M. G. C., ... Draiva, A., et al., (2021). A light-touch routing optimization tool (roOT) for vaccine and medical supply distribution in mozambique. *International Transactions in Operational Research, 28*(5), 2334–2358.

Lim, J., Norman, B. A., & Rajgopal, J. (2019). Redesign of vaccine distribution networks. *International Transactions in Operational Research, 29*(1), 200–225.

Lurie, N., Saville, M., Hatchett, R., & Halton, J. (2020). Developing COVID-19 vaccines at pandemic speed. *New England Journal of Medicine, 382*(21), 1969–1973.

Ma, Y., Pei, S., Shaman, J., Dubrow, R., & Chen, K. (2021). Role of meteorological factors in the transmission of SARS-CoV-2 in the United States. *Nature Communications, 12*(1), 1–9.
Manupati, V. K., Schoenherr, T., Subramanian, N., Ramkumar, M., Soni, B., & Panigrahi, S. (2021). A multi-echelon dynamic cold chain for managing vaccine distribution. *Transportation Research Part E: Logistics and Transportation Review, 156*, 102542.

Mohammadi, M., Dehbari, S., & Vahdani, B. (2014). Design of a bi-objective reliable healthcare network with finite capacity queue under service covering uncertainty. *Transportation Research Part E: Logistics and Transportation Review, 72*, 15–41.

Mohammadi, M., Dehghan, M., Pirayesh, A., & Dolgui, A. (2022). Bi-objective optimization of a stochastic resilient vaccine distribution network in the context of the COVID-19 pandemic. *Omega, 113*, 102725.

Mohammadi, S., Darestani, S. A., Vahdani, B., & Alinezhad, A. (2020). A robust neu-trosophic fuzzy-based approach to integrate reliable facility location and routing decisions for disaster relief under fairness and aftershocks concerns. *Computers and Industrial Engineering, 148*, 106734.

Nabavi, S. M., Vahdani, B., Nadjafi, B. A., & Adibi, M. A. (2022). Synchronizing vic-tim evacuation and debris removal: A data-driven robust prediction approach. *European Journal of Operational Research, 300*(2), 689–712.

Nguyen, C., & Carlson, J. M. (2016). Optimizing real-time vaccine allocation in a stochastic SIR model. *PLoS One, 11*(4), e0152950.

Ni, W., Shu, J., Song, M., Xu, D., & Zhang, K. (2021). A branch-and-price algorithm for facility location with general facility cost functions. *INFORMS Journal on Computing, 33*(1), 86–104.

Niakan, F., Vahdani, B., & Mohammadi, M. (2015). A multi-objective optimization model for hub network design under uncertainty: An inexact rough-interval fuzzy approach. *Engineering Optimization, 47*(12), 1670–1688.

Nikolopoulos, K., Punia, S., Schäfers, A., Tsinopoulos, C., & Vasilakis, C. (2021). Forecasting and planning during a pandemic: COVID-19 growth rates, supply chain disruptions, and governmental decisions. *European journal of operational research, 290*(1), 99–115.

OPECST (2021). Épidémie de COVID-19—point sur la modélisation épidémiologique pour estimer l’ampleur et le devenir de l’épidémie de COVID-19. https://www.senat.fr/espace_presse/actualites/202003/suivi_covid_senat/opecst.html.

Pagliusi, S., Jarrett, S., Hayman, B., Kreysa, U., Prasad, S. D., Reers, M., ... Baek, Y. O., et al., (2020). Emerging manufacturers engagements in the COVID-19 vaccine research, development and supply. *Vaccine, 38*(34), 5418–5423.

Pontryagin, L. S. (1987). *Mathematical theory of optimal processes*. CRC Press.

Qi, H., Xiao, S., Shi, R., Ward, M. P., Chen, Y., Tu, W., ... Zhang, Z. (2020). Covid-19 transmission in mainland China is associated with temperature and humidity: A time-series analysis. *Science of the Total Environment, 728*, 138778.

Rahimi-Vahed, A., Crainic, T. G., Gendreau, M., & Rei, W. (2015). Fleet-sizing for multi-depot and periodic vehicle routing problems using a modular heuristic algorithm. *Computers and Operations Research, 53*, 9–23.

Rahmanian, R., Crainic, T. G., Gendreau, M., & Rei, W. (2018). Accelerating the Benders decomposition method: Application to stochastic network design problems. *SIAM Journal on Optimization, 28*(1), 875–903.

Rastegar, M., Tavana, M., Meraj, A., & Mina, H. (2021). An inventory-location optimization model for equitable influenza vaccine distribution in developing countries during the COVID-19 pandemic. *Vaccine, 39*(3), 495–504.

Rostami, B., Desaulniers, G., Errico, F., & Lodi, A. (2021). Branch-price-and-cut algorithms for the vehicle routing problem with stochastic and correlated travel times. *Operations Research, 69*(2), 436–455.

Sampath, V., Rabinowitz, G., Shah, M., Jain, S., Diamant, Z., Jesenak, M., Vieths, S., Agache, I., Akdis, M., Barber, D., et al., (2021). Vaccines and allergic reactions: The past, the current COVID-19 pandemic, and future perspectives. *Allergy, 76*(6), 1640–1660.

Santoso, T., Ahmed, S., Goetschalckx, M., & Shapiro, A. (2005). A stochastic programming approach for supply chain network design under uncertainty. *European Journal of Operational Research, 167*(1), 96–115.

Shaw, R. H., Stuart, A., Greenland, M., Liu, X., Van-Tam, J. S. N., & Snape, M. D. (2021). Heterologous prime-boost COVID-19 vaccination: initial reactogenicity data. *The Lancet, 397*(10289), 2043–2046.

Silva-Cayetano, A., Foster, W. S., Innocentin, S., Belij-Rammerstorfer, S., Spencer, A. J., Burton, O. T., ... Conceicao, C., et al., (2021). A booster dose enhances immunogenicity of the COVID-19 vaccine candidate ChAdOx1 nCoV-19 in aged mice. *Med, 2*(3), 243–262.

Tang, L., Li, Y., Bai, D., Liu, T., & Coelho, L. C. (2022). Bi-objective optimization for a multi-period COVID-19 vaccination planning problem. *Omega, 110*, 102617.

Tang, Z., Jiao, Y., & Ravi, R. (2022). Combinatorial heuristics for inventory routing problems. *INFORMS Journal on Computing, 34*(1), 370–384.

Tavana, M., Govindan, K., Nasr, A. K., et al., (2021). A mathematical programming approach for equitable COVID-19 vaccine distribution in developing countries. *Ann Oper Res*. https://doi.org/10.1007/s10479-021-04130-z.

Thul, L., & Powell, W. (2023). Stochastic optimization for vaccine and testing kit allocation for the COVID-19 pandemic. *European Journal of Operational Research, 304*(1), 325–338.

Tyagi, K., Ghosh, A., Nair, D., Dutta, K., Bhandari, P. S., Ansari, I. A., & Misra, A. (2021). Breakthrough COVID19 infections after vaccinations in healthcare and other workers in a chronic care medical facility in New Delhi, India. *Diabetes and Metabolic Syndrome: Clinical Research and Reviews, 15*(3), 1007–1008.

Vahdani, B., & Ahmadzadeh, E. (2019). Designing a realistic ICT closed loop supply chain network with integrated decisions under uncertain demand and lead time. *Knowledge-Based Systems, 179*, 34–54.

Vahdani, B., & Ahmadzadeh, E. (2021). Incorporating price-dependent demands into a multi-echelon closed-loop network considering the lost sales and backorders: A case study of wireless network. *Networks and Spatial Economics, 21*, 639–680.

Vahdani, B., Soltani, M., Yazdani, M., & Mousavi, S. M. (2017). A three level joint location-inventory problem with correlated demand, shortages and periodic review system: Robust meta-heuristics. *Computers and Industrial Engineering, 109*, 113–129.

Vahdani, B., Zandieh, M., & Roshanaei, V. (2011). A hybrid multi-stage predictive model for supply chain network collapse recovery analysis: A practical framework for effective supply chain network continuity management. *International Journal of Production Research, 49*(7), 2035–2060.

Yeysmoradi, D., Vahdani, B., Farhadi Sartangi, M., & Mousavi, S. M. (2018). Multi-objective open location-routing model for relief distribution networks with split delivery and multi-mode transportation under uncertainty. *Scientia Iranica, 25*(6), 3635–3653.

Vidal, T., Crainic, T. G., Gendreau, M., Lahrichi, N., & Rei, W. (2012). A hybrid genetic algorithm for multidepot and periodic vehicle routing problems. *Operations Research, 60*(3), 611–624.

Vidal, T., Martinelli, R., Pham, T. A., & Hà, M. H. (2021). Arc routing with time-dependent travel times and paths. *Transportation Science, 55*(3), 706–724.

Wallinga, J., van Boven, M., & Lipsitch, M. (2010). Optimizing infectious disease interventions during an emerging epidemic. *Proceedings of the National Academy of Sciences, 107*(2), 923–928.

Westerink-Duijzer, L. E., Schlicher, L. P. J., & Musegaas, M. (2020). Core allocations for cooperation problems in vaccination. *Production and Operations Management, 29*(7), 1720–1737.

Wouters, O. J., Shadlen, K. C., Salcher-Konrad, M., Pollard, A. J., Larson, H. J., Teer-awattananon, Y., & Jit, M. (2021). Challenges in ensuring global access to COVID-19 vaccines: Production, affordability, allocation, and deployment. *The Lancet, 397*(10278), 1023–1034.

Yang, Y., Bidkhori, H., & Rajgopal, J. (2021). Optimizing vaccine distribution networks in low and middle-income countries. *Omega, 99*, 102197.

Yang, Z., Zeng, Z., Wang, K., Wong, S.-S., Liang, W., Zanin, M., ... Mai, Z., et al., (2020). Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions. *Journal of Thoracic Disease, 12*(3), 165.

Yin, X., Büyüktahtakin, I. E., & Patel, B. P. (2023). COVID-19: Data-driven optimal allocation of ventilator supply under uncertainty and risk. *European Journal of Operational Research, 304*(1), 255–275.

Zheng, X., Yin, M., & Zhang, Y. (2019). Integrated optimization of location, inventory and routing in supply chain network design. *Transportation Research Part B: Methodological, 121*, 1–20.