An equitable and accessible vaccine supply chain network in the epidemic outbreak of COVID-19 under uncertainty

Abstract

Vaccination is one of the most efficient ways to restrict and control the spread of epidemic outbreaks such as COVID-19. Due to the limited COVID-19 vaccine supply, an equitable and accessible plan should be prepared to cope with. This research focuses on designing a vaccine supply chain while aiming to achieve an equitable and accessible network. We present a novel mathematical formulation that helps to optimize vaccine distribution to inoculate people with various priority levels to achieve an equitable plan. The transshipment strategy is also incorporated into the model to enhance the accessibility of COVID-19 vaccine types between health facilities. The nature of COVID-19 is dynamic over time due to mutations, and the protection level of each vaccine type against this disease is not exact. Besides, complete information about the demand for different vaccine types is not available. Hence, we use Multi-Stage Stochastic Programming as a reliable strategy that is organized to manage stochastic data in a dynamic environment for the first time in the vaccine supply chain network. The scenarios in this approach are generated using a Monte Carlo simulation method, and then a forward scenario reduction technique is conducted to construct a suitable scenario tree. The practicality and capability of the model are shown in a real-life case of Iran. The results show that the performance of the Multi-Stage Stochastic Programming is significantly improved compared with the two-stage stochastic programming regarding the total cost of the vaccine supply chain and the number of the shortage units.

1 Introduction

1.1 Background and real motivational case

In the last decades, pandemic diseases such as the Ebola virus, Influenza virus, Severe Acute Respiratory Syndrome (SARS) and Middle East Respiratory Syndrome (MERS) have spread worldwide. Unfortunately, these pandemics have infected a large population of global and killed millions. Since late 2019, the world has been challenging by a new pandemic disease called Coronavirus or COVID-19 with its fast spread and raised death rate (Alam et al. 2021). With 233,092,663 infected cases and 4,769,789 deaths by September 28th, 2021 (www.worldometers.info), COVID-19 has become significant fear and stress worldwide. During this pandemic, several mutated versions of the Coronavirus have been detected in different countries. COVID-19 pandemic has been the firmest risk in the global supply chain during the last decade (Araz et al. 2020); therefore, it needs an appropriate plan to manage and prevent from infecting people to COVID-19.

Vaccination operations can play a crucial role in controlling and preventing the outbreak of COVID-19 (Yang et al. 2021). Simultaneously, as the COVID-19 outbreak, well-known companies in vaccine production, such as Pfizer, have started their studies on various aspects of COVID-19. They have tried to produce a particular vaccine and reach the results and samples of vaccines by the end of 2020 (www.pfizer.com). So far, several types of COVID-19 vaccines such as Pfizer, AstraZeneca, and Moderna have been produced according to standards approved by the World Health Organization (WHO). General vaccination needs to be supported by governmental regulators, equitably distributed to health facilities in different regions (Bollyky et al. 2020), and distributed among the people considering the priority policy (Choi et al. 2021). Since companies’ production capacity is limited and, on the other hand, countries’ demand is more than produced vaccines; therefore, covering all vaccine demands for a country may be delayed by a long time, leading to inevitable shortages in health facilities. Hence, some people who belong to a sensitive group, such as health facility personnel, patients, and older people, should be received the vaccine before other groups (Shim 2021). Therefore, prioritizing people based on medical status and closeness to infected patients is appropriate for distributing vaccines (Persad et al. 2020). Besides, considering equity policy in public distributions health networks is critical. Disregarding equity constraints is unacceptable and unfair during to COVID-19 outbreak (Enayati and Özaltın 2020; Mitchell et al. 2009).

Due to the high volatility of demand and effectiveness rate for each vaccine type during the COVID-19 outbreak, the demand for each type is an uncertain parameter. Hence, to respond to this uncertainty, Decision-Makers (DMs) should use reliable strategies. Uncertainty and time are two critical factors for DMs when faced with dynamic decisions, and they must be adopted appropriate decision strategies to respond to changes. Usually, when DMs access new recourses and information, these decisions work. One of these strategies is Multi-Stage Stochastic Programming (MSSP), which can help design the strategic decisions and respond to unknown prior (Zahiri et al. 2018). This approach has been used for design the supply chain network in a limited number of studies (Birge and Louveaux 2011; Fattahi et al. 2017; Zahiri et al. 2017). This approach has several decisions to be made throughout a path of uncertain events in a multi-period setting. At every decision stage, the best performance is investigated for observed previous occasions and anticipation of unknown future realization, while the Two-Stage Stochastic Programming (TSSP) approach is utilized only in a single-period setting (Zahiri et al. 2017). Finally, the MSSP aims to decrease forced risk via a scenario tree presenting data uncertainty and its potential influence on the sequence decisions (Zahiri et al. 2018).

Additionally, transshipment is one available chance to increase a health network’s performance (Dehghani et al. 2021; Samani et al. 2020) and consists of any vaccine transfers between health facilities. In other words, transshipment helps hospitals and health centers in allocating shortages further efficiently by utilizing nearby hospitals/health centers’ inventory. Notably, diverse conditions should be considered for each type of COVID-19 vaccine due to their temperature-sensitive. For example, the Moderna vaccine is more fragile compared to other types of vaccines and should minimize excessive movement during transshipment as possible. All vaccines must be packed in a cooler with specific degrees securely with bubble wrap to protect them (Gliatto et al. 2021) and moved safely by the appropriate vehicle. This paper shows how transshipment can be considered an inventory control system as a great strategy to rebalance the various types of vaccine inventory in the health centers and hospitals and eventually decrease the shortage. Against this, some health facilities like urban hospitals across the country may need special equipment such as ultra-cold freezers to store extra COVID-19 vaccines (AboulFotouh et al. 2021). Hence, different conditions for storing each type of vaccine should be considered to maintain.

According to the concerns mentioned above, in this study, a Mixed Integer Linear Programming (MILP) model is developed for designing a three-level distribution network of COVID-19 vaccine. This study is an idea and could potentially defeat the enormous challenge caused by COVID-19. The concerned model considers the choice of proper hub distribution centers, the supply capacity of logistic centers and hub distribution centers, allocation of health facilities to the opened hub distribution centers, and opened hub distribution centers to logistic centers. This problem aims to minimize the total cost of the network with respect to equitable distribution of all types of vaccines among regions, priority for different groups of people, and uncertain demand. Additionally, we use MSSP to cope with the uncertainty of input parameters in model formulation. In MSSP, we generate scenarios by a simulation approach and then obtain a tree scenario from a reduction method. The method makes the vaccine supply chain resilient to shortage and uncertainty.

1.2 Research questions

This paper presents a MILP mathematical model for designing an efficient vaccine supply chain under uncertainty to prevent and control COVID-19 transmission. Motivated by a real case in Tehran, a logistic center helps the regions imply the vaccination operations. This research aims to answer the following questions:

• How to consider an optimal distribution policy for covering the demands of people?

• How to inject vaccines for different groups of people based on their priorities?

• How to consider an equitable distribution among all regions?

• How can transshipment consider an appropriate mechanism to rebalance the vaccine accessibility in the vaccine supply chain network?

• Which parameter deals with uncertainty and how to cope with uncertainty?

The rest of this research is structured as follows: Sect. 2 includes a survey of the literature. In Sect. 3, the problem definition and the mathematical formulation are presented. The solution methods are elaborated in Sect. 4. Then, the computational results of the proposed model by a real case study are provided in Sect. 5. Eventually, concluding remarks and future research are gathered in the last Sect. 6.

2 Literature review

As discussed in the previous section, injecting vaccines for different groups of people based on their priorities, equitable distribution among all regions, and transshipment among health facilities for rebalancing the vaccine accessibility are vital during the COVID-19 outbreak. Since designing the COVID-19 vaccine supply chain network is seldom discussed in the literature, this paper focuses on supply chain design issues for the COVID-19 vaccine. To highlight the novelties of this study relative to the prior literature, in this section, two relevant research streams to the research are reviewed, applications of Operation Research (OR) during the COVID-19 outbreak and applications of the vaccine supply chain. Finally, the research gaps are recognized, and the contributions of this study are discussed.

2.1 Applications of OR during the COVID-19 outbreak

The COVID-19 is an innovative coronavirus spread from the Wuhan region in China in late 2019 (Choi 2020). The Coronavirus is transmitted human-to-human via the droplets generated from talking, coughing, and sneezing (WHO). Besides, necessary commodities like medical items and foods are required in the supply chain during the COVID-19 outbreak. The reduction creates a severe weakness in the global supply chain, causing uneven supply and demand requirements (Linton and Vakil 2020).

Several authors reviewed recent papers on this field’s characteristics, impacts of this outbreak, and emergency resources allocation during COVID-19. Some of these papers provide a strategic controller for this pandemic. For example, Govindan et al. (2020) divided citizens into four groups, normal, slightly sensitive, sensitive, and so sensitive. For this classification, they utilized some information based on the immune level of the body, pre-existing illnesses (like heart attack and diabetes), and age. They developed an effective support method based on physicians’ knowledge and fuzzy inference policy to help manage demand by a decision-making method for healthcare supplies, decrease fear and stress in society, and eliminate the outbreak chain of COVID-19. Ivanov (2020) presented a simulation approach for global supply chain disruption and investigated a set of sensitivity analyses for the COVID-19 scenario. Their results showed that facilities’ opening and closing timing at various levels in supply chains is the main factor for managing the epidemic outbreak effects on the network. Their method can help recognize the correct and incorrect risk reduction/preparation policies during pandemic outbreaks. Ivanov and Dolgui (2020) described that the presented OR approaches could help to cope with the impact of the COVID-19 outbreak on the forecast, discovery, containment, control, and removal steps following the classification of WHO. A multi-period stochastic programming model is proposed by Mehrotra et al. (2020) to distribute the ventilator inventory during the COVID-19 pandemic. They generated four different types of demand scenarios for a case study in the United States. Shokrani et al. (2020) studied creating alternative supply chains to produce medical devices like face shields (protective equipment) during COVID-19.

Gray (2020) studied the different shipping risks during COVID-19 outbreak. They enumerated these risks involve labor issues, regulatory closures, and decreased marine container movements of shipping services. Paul and Chowdhury (2020) investigated several policies for high demand and necessary commodities like toilet paper during COVID-19 outbreak and suggested a strategy to deal with disruption. Queiroz et al. (2020) presented a framework for supply chain management during COVID-19 outbreak considering six perspectives. Analytics mechanisms for forecasting the growth rate and planning of COVID-19 are proposed by Nikolopoulos et al. (2021). Their forecast policy is based on statistical, machine, epidemiological, and a new hybrid forecasting system by closest clustering and neighbors. Murray (2020) predicted the stress on ICUs, hospitals, and ventilators and considered the overall demand beyond the healthcare system’s current capacity. The author investigated the effect of COVID-19 on hospitals and the number of deaths in the United States and Europe. The researcher also suggested measures to increase the supply of necessary services and commodities temporarily. Onar et al. (2020) developed a decision method for determining the inventory of healthcare facilities under multiple scenarios of pandemics like COVID-19. They presented a new spherical regret-based multi-criteria decision-making approach and utilized it for assessing the regrets of not keeping stock in the healthcare facilities.

In the waste disposal scopes, Kargar et al. (2020) presented a linear programming model for reverse logistics network design with three objective functions during COVID-19. They evaluated the risk associated with the shipping and medicine of infectious medical waste and the number of uncollected wastes in medical waste generation centers. Tirkolaee et al. (2021) developed a new MILP model for the sustainable multi-trip location-routing problem. They considered time windows for medical waste management during COVID-19 to define hospitals’ service priority with various disruptions. They also considered minimizing transportation time, environmental disturbance, and the violation of service priorities required on nearby disposal places as objective functions. Besides, they utilized a fuzzy chance-constrained programming method for the model to face uncertainty. Saadat et al. (2020) addressed medical waste disposal in the environment as a hugely challenging problem and evaluated COVID-19 environmental aspects. In another paper, Wang et al. (2020) proposed several procedures to disinfect hospital waste/wastewater during the COVID-19 outbreak in China. Then, Manupati et al. (2021) suggested a framework for choosing the best technique for healthcare waste disposal according to socio-technical and triple bottom line aspects during COVID-19.

A review of literature on this section, a substantial proportion of the papers has worked on OR techniques to examine emergency resources allocation and their economic aspects during COVID-19. OR technics can design an efficient vaccine supply chain network.

2.2 Applications of the vaccine supply chain

Scholars have been interested in the vaccine supply chain from several aspects due to vaccines’ critical purpose in preventing infectious diseases outbreaks. It has been shown by Yang et al. (2021) that vaccination is an efficient way to prevent infectious disease transmission. On-time vaccine distribution plays a significant role in the performance and effectiveness of the disaster response. Mainly, modeling a supply chain for proper vaccine distribution affects the health network’s performance during the COVID-19 pandemic (AlTheeb and Murray 2017). Yang et al. (2021) designed a vaccine distribution network and formulated a Mixed Integer Programming (MIP) model for low and middle-income countries. They considered the costs of the hub facility and round-trip transportation as the objective function. They also proposed a novel algorithm for general problems with large size and examined their algorithm applying data obtained from four various countries of Africa. Some researchers have investigated the different aspects of the vaccine supply chain. Lin et al. (2020) focused on the cold supply chain of vaccines sensitive to temperature. They studied the distributor’s aims to move vaccines by a cold chain or non-cold chain. Dai et al. (2020) considered a supply chain of vaccine transportation with a distributor and a retailer. They proposed a model of decision-making time-delay problem that the decision-makers do not make a decision immediately. The findings proved that the system would go toward chaos or miss stability if the decision delay passed a specific threshold. Chandra and Kumar (2019) tried to recognize and investigate the vaccine supply chain’s critical issues for immunizing kids in developing countries.

The priority of various people groups is the most useful policy in discounting the pandemics’ impacts on vaccines’ distribution (Uscher-Pines et al. 2006). This policy plays an essential role in reducing the total social cost related to an infected person (the medicine drugs, hospitalization cost, and other disease costs) (Torabi 2018). Healthcare providers, vital services providers, and people with a high risk of outbreaks are priority groups for vaccination recommended by WHO. In reality, vaccine prioritization aims to decrease death, decrease social shocks, and financial disasters (Straetemans et al. 2007). There are some papers in which priority classes are considered by simulation and mathematical models under various objectives. Lee et al. (2010) presented a simulation based on the computer of group allocation and prioritization during the Influenza outbreak in 2009. Pregnant women, kids under six months old, workers of the health system, people with 6 months–24 years old, and people with 25–64 years are classified as five priority groups by authors. Hu and Zhang (2013) established a novel algorithm for setting the number of vaccine doses allotted to various age groups to reduce infected people during pandemics. Hovav and Herbon (2017) introduced a new vaccination distribution model with the high-risk people priority to minimize non-medical and medical costs. Gamchi et al. (2021) presented a vehicle routing problem for distributing vaccines to manage the outbreak of pandemics. They considered two objective functions to minimize the used vehicle cost and the social cost by assuming the various priority groups under the Susceptible-Infected-Recovered (SIR) epidemic model. Besides, equitable distribution is critical in the vaccine supply chain network during an outbreak. Enayati and Özaltın (2020) developed the optimal influenza vaccine distribution with multiple subgroups of people. They utilized a compartmental model for influenza transmission and presented a mathematical model to minimize the number of vaccines to control the outbreak effectively. The authors also considered an equity policy to make vaccine distribution decisions.

To the best of our knowledge and reviewed papers, there is no study yet dealing with an efficient network for distributing the vaccines of COVID-19 that have taken into account transshipment between health facilities to increase accessibility and equitable distribution. However, we seek to apply OR techniques to design an efficient COVID-19 vaccine supply chain network to fill the gap.

2.3 Contributions and research gaps

Regarding literature, no paper has been done to focus on the equitable and accessible COVID-19 vaccination service for these new challenges that the world is facing. This research connects the gap between study and application and helps make these relevant decisions through a mathematical programming model. According to the mentioned issues, the purpose of this study is to develop a MILP model for designing a vaccine supply chain network with a multi-stage planning horizon under uncertain demand. In summary, the main contributions for this research explain as follows:

• Proposing a novel mathematical model for designing a multi-period vaccine supply chain network considering national logistic centers, potential hub centers, hospitals, and health centers.

• Formulating an equity strategy for distributing the different types of COVID-19 vaccine in the regions.

• Classifying the priority groups of people for injecting the COVID-19 vaccine.

• Proposing a transshipment strategy as a powerful mechanism in an inventory control system to access the COVID-19 vaccine in the health centers and hospitals and eventually reduce the shortage costs.

• Applying the multi-stage stochastic programming to deal with uncertain demand for different types of COVID-19 vaccines.

• Using a scenario generation method by a Monte Carlo simulation method and a forward scenario reduction method to convert scenario fans into a scenario tree.

• Applying the proposed model for a real case of Tehran in Iran during the COVID-19 outbreak to address the model’s practicality.

3 Problem description and mathematical formulation

In this paper, a new mathematical model for a vaccine supply chain network is developed. To do so, we present a network consists of three levels, national logistic centers, potential hub centers, and health facilities (hospitals and health centers) in regions. In the first echelon, vaccines are driven from the national logistic center and then equitability transferred to the opened hubs. In the second echelon, vaccine units moved from the opened hub to health facilities based on their demands. Finally, in the third echelon, vaccines are injected into people in the health facilities, and then the leftover vaccines can be stored or be transshipped from a health facility to other health facilities in this network. Notably, vaccines are specific medicines that are fragile and very sensitive to temperature, unlike most drugs. Therefore, a portable vaccine refrigerator or qualified pack out should be used for transport vaccines to keep them safe. A qualified pack out includes a container and supplies specifically designed for use when packing each type of vaccine for transport. Also, a suitable pack out is “qualified” by laboratory testing under controlled conditions to ensure it could achieve and maintain desired temperatures for a set amount of time (www.cdc.gov/vaccines/covid-19/clinical-considerations/homebound-persons.html). Hence, these considerations for transporting each type of vaccine to the centers should be used. Besides, prioritizing vaccines will contribute to an efficient and effective immunization plan as the primary goal of vaccination is to control pandemic faces with infectious diseases. Another important decision is equitable allocations concerning distributing vaccines to all districts. Mentioned policies are described as follows:

• 1) Equity policy: Each person is arranged inside a particular group of risk priorities based on sensitivity risk. People are allocated to have different groups based on age group, job type, and previous diseases by the segmentation system. The combinations of the criteria are used to classify priority groups into High-risk and Low-risk groups. The proposed model allows the prioritizing of people by assigning them to pre-determined risk clusters and next allocating them to the nearest health facilities.

• 2) Priority strategy: Another subject concerns the equity policy of vaccine distribution. Proposing a fair and equitable plan to distribute different types of the vaccine among people, which have been stored in the logistic centers or hub distribution centers, is an undeniable truth and vital measure. The equity level could be set by a maximal ratio to determine the fittest allocation of vaccines according to the available supply.

• 3) Transshipment policy: Lateral transshipment consists of every inventory movement between hospitals and health centers in regions. This approach conventionally adjusts inventory in health facilities by reallocating. This means transshipment can be applied as a useful way to balance the current lack between the demand of the present/future and the inventory of different vaccine types among health centers and hospitals.

This research aims to minimize the total cost, and Fig. 1 demonstrates a schematic overview of the concerned network.

Fig. 1

An overview of the concerned network

The following assumptions are considered in this paper:

• Optimal locations of hub distribution centers are selected from the candidate sites.

• The different types of vaccines are distributed equitably among the regions based on their demand.

• Due to the high volatility of demand for different types of vaccines during COVID-19 outbreak, the demand is uncertain.

• Two priority categories for injection are considered that include the high-risk and low-risk group of people.

• The transshipment strategy is considered between health facilities to prevent shortages.

• The initial inventory at the health facilities at the beginning of the first period is equal to zero.

• Each type of vaccine is transported between all nodes by its special cold chain equipment.

3.1 Model formulation

This section presents a new mathematical MILP model for designing a vaccine supply chain network. In the beginning, the sets, parameters, and variables employed in the concerned model are defined in Table 1.

Table 1 Nomenclature, list of sets, parameters, binary and positive variables

In the following, the mathematical model with respect to mentioned notifications is formulated. It is notable that the “minimize” function is abbreviated by “$Min$” for the objective function.

$$Min \sum_{h}{cw}_{h}{w}_{h}+\sum_{v}\sum_{l}\sum_{h}\sum_{t}2\times{cd}_{vlh}^{^{\prime}}{f}_{lh}{u}_{lht}^{^{\prime}}+\sum_{v}\sum_{h}\sum_{c}\sum_{t}2\times{cd}_{vhc}{f}_{hc}{u}_{hct}+\sum_{v}\sum_{{t}_{1}}\sum_{c}\sum_{{c}^{^{\prime}}\ne c}2\times{cz}_{vc{c}^{^{\prime}}}{f}_{c{c}^{^{\prime}}}{z}_{vc{c}^{^{\prime}}t}+\sum_{v}\sum_{c}\sum_{t}{ck}_{vc}{k}_{vct}+\sum_{v}\sum_{c}\sum_{t}{ck}_{vc}^{^{\prime}}{k}_{vct}^{^{\prime}}+\sum_{v}\sum_{c}\sum_{t}{cs}_{vc} {In}_{vct}$$

(1)

The objective function (1) minimizes the total cost that includes seven components: fixed cost of opening hub distribution centers ($\sum_{h}{cw}_{h}{w}_{h}$), traveling cost from the national logistic center to hub distribution centers ($\sum_{v}\sum_{l}\sum_{h}\sum_{t}2\times{cd}_{vlh}^{^{\prime}}{f}_{lh}{u}_{lht}^{^{\prime}}$), traveling cost from hub distribution centers to health facilities ($\sum_{v}\sum_{h}\sum_{c}\sum_{t}2\times{cd}_{vhc}{f}_{hc}{u}_{hct}$), transshipment cost between health facilities ($\sum_{v}\sum_{{t}_{1}}\sum_{c}\sum_{{c}^{^{\prime}}\ne c}2\times{cz}_{vc{c}^{^{\prime}}}{f}_{c{c}^{^{\prime}}}{z}_{vc{c}^{^{\prime}}t}$), penalty cost of shortage in demand for high-risk ($\sum_{v}\sum_{c}\sum_{t}{ck}_{vc}{k}_{vct}$), low-risk group ($\sum_{v}\sum_{c}\sum_{t}{ck}_{vc}^{^{\prime}}{k}_{vct}^{^{\prime}}$), and vaccine storing cost in health facilities ($\sum_{v}\sum_{c}\sum_{t}{cs}_{vc} {In}_{vct}$), respectively.

Assignment constraints:

$$\sum_{l}{u}_{lht}^{^{\prime}}\le 1 \forall h\in H,t\in T$$

(2)

$$\sum_{h}{u}_{hct}\le 1 \forall c\in C,t\in T$$

(3)

$${u}_{hct}\le {w}_{h} \forall h\in H, c\in C,t\in T$$

(4)

$${x}_{vlht}^{^{\prime}}\le M {u}_{lht}^{^{\prime}} \forall v\in V, l\in L,h\in H, t\in T$$

(5)

$${x}_{vhct}\le M{u}_{hct } \forall v\in V,h\in H,c\in C, t\in T$$

(6)

Constraint (2) guarantees that each hub distribution center can be assigned to one national logistic center. Constraint (3) ensure that each health facility can be assigned to one opened hub distribution center. Constraint (4) determines if a hub distribution center is opened, a health facility can be allocated to that hub distribution center. Constraint (5) signifies that the different vaccine types can be transferred from the national logistic center to the hub distribution center if they are connected. Constraint (6) assures that the different vaccine types can be transferred from the hub distribution center to health facilities if there is a connection between them.

Inventory constraints:

$${I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{v{c}^{^{\prime}}c{t}_{1}}={e}_{vc{t}_{1}}+{e}_{vc{t}_{1}}^{^{\prime}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{v{cc}^{^{\prime}}{t}_{1}}+{In}_{vc{t}_{1}} \forall v\in V, c\in C$$

(7)

$${I}_{vct}={e}_{vct}+{e}_{vct}^{^{\prime}}+{In}_{vct} \forall v\in V, c\in C, t\in T,t\ne {t}_{1}$$

(8)

$${I}_{vc{t}_{1}}={q}_{vc{t}_{1}} \forall v\in V, c\in C$$

(9)

Constraints (7)-(9) update the inventory of each vaccine type for each period in health facilities. Constraint (7) sets the balance of vaccine units at the first period considering transshipment between health facilities, while Constraint (8) establishes the same balance for the remaining periods without considering transshipment. The transshipment can be done at the beginning of the time horizon. Constraint (9) indicates that the initial inventory of each type of vaccine is known beforehand in the first period.

Demand constraints for different vaccine types in health facilities:

$${e}_{vct}+{k}_{vct}={\tilde{d }}_{vct} \forall v\in V, c\in C, t\in T$$

(10)

$${e}_{vct}^{^{\prime}}+{k}_{vct}^{^{\prime}}={\tilde{d }}_{vct}^{^{\prime}} \forall v\in V, c\in C, t\in T$$

(11)

In this problem, all demands of different groups may not be fulfilled. For this purpose, Constraint (10) indicates that the demands for each vaccine type will be satisfied by fulfillment demands for the high-risk group in health facilities; otherwise, they will deal with a shortage. Constraint (11) also indicates that demand for each vaccine type will be satisfied by fulfillment demands for the low-risk group in health facilities; otherwise, they will deal with a shortage.

Order constraints:

$${y}_{vct}=\sum_{h}{x}_{vhc(t+1)} \forall v\in V, c\in C, t\in T,t\ne \left|T\right|$$

(12)

$${In}_{vct}+\sum_{h}{x}_{vhc(t+1)}={I}_{vc(t+1)} \forall v\in V, c\in C, t\in T$$

(13)

$${To}_{vc}-{I}_{vct}={y}_{vct} \forall v\in V, c\in C, t\in T,t\ne 1$$

(14)

Constraints (12)–(14) update the order quantity of each vaccine type in health facilities at each period. Constraint (12) states that the order of each vaccine type for each health facility at a given period arrives at the end of that period (it is used at the beginning of the next period), which is assumed that the lead time is one period. Constraint (13) illustrates that the total inventory of each vaccine type at the end of a given period and the received vaccine units from hub distribution centers at the next period are used to fulfill the initial inventory at the next period. The level of target inventory in a health facility is calculated by Constraint (14).

Flow and capacity constraints:

$$\sum_{l}{x}_{vlht}^{^{\prime}}=\sum_{c}{x}_{vhct }\forall v\in V, h\in H, t\in T$$

(15)

$$\sum_{h}\sum_{t}{x}_{vlht}^{^{\prime}}\le {S}_{vl}^{^{\prime}} \forall v\in V, l\in L$$

(16)

$$\sum_{c}{x}_{vhct}\le {S}_{vh} \forall v\in V,h\in H, t\in T$$

(17)

$${In}_{vct}\le {S}_{vc}^{\prime\prime} \forall v\in V,c\in C, t\in T$$

(18)

Constraint (15) balances the flows in hub distribution centers and health facilities. In other words, the number of transferred vaccine units from national logistic centers is equal to the number of transferred vaccine units from hub distribution centers to health facilities. Constraints (16) and (17) are capacity constraints limiting each vaccine type’s outcoming supply from national logistic and hub distribution centers. Constraint (18) bounds the capacity limitation of health facilities for each type of vaccine. Two decision variables are defined: ${x}_{vlht}^{^{\prime}}$, which determines the number of each vaccine-type in each national logistic center that should be transferred to the hub distribution center; and ${x}_{vhct}$, which determines the number of vaccine types in each hub distribution center that should be transferred to health facilities.

Priority constraints:

The allocation of available COVID-19 vaccines for injecting to the high-risk group is of great importance. The impact of prioritization on the vaccination of people cannot be overlooked for controlling the COVID-19 outbreak; therefore, it is important in this model. The following constraints can formulate this policy:

$${e}_{vc{t}_{1}}=min\left\{{\tilde{d }}_{vc{t}_{1}},{I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}\right\} \forall v\in V, c\in C$$

(19)

$${e}_{vc{t}_{1}}^{^{\prime}}=min\left\{{\tilde{d }}_{vct}^{^{\prime}},{I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}\right\} \forall v\in V, c\in C$$

(20)

$${e}_{vct}=min\left\{{\tilde{d }}_{vct},{I}_{vct}\right\} \forall v\in V, c\in C, t\in T\backslash t={t}_{1}$$

(21)

$${e}_{vct}^{^{\prime}}=min\left\{{\tilde{d }}_{vct}^{^{\prime}},{I}_{vct}-{e}_{vct}\right\} \forall v\in V, c\in C, t\in T\backslash t={t}_{1}$$

(22)

Constraints (19)-(22) determine that the available vaccine units should be injected into the high-risk group of people in the first step and then serviced to other groups based on their priority. Constraint (19) enforces that the available vaccine units fulfill the demand of the high-risk group in each health facility (${I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}$) during the first period. It means that if the available vaccine units are equal to (or less than) the high-risk group’s demand, no vaccine units will remain for another group, and it faces a shortage. Otherwise, if the available vaccine unit is more than the demand of the high-risk group, the demand of the low-risk group will be fulfilled with the remains of the available vaccine units (${I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}$) at the first period in Constraint (20). Also, if the remains of the available vaccine units (${I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}$) are more than the demand of the low-risk group; in this case, none of the groups deal with shortages. In this case, the remain vaccine units are stored as inventory or shared among health facilities via a transshipment mechanism. If the vaccine units are insufficient to fully satisfy the high-risk group’s demand, the high-risk group’s remaining demand and the low-risk group’s demands deal with a shortage by Constraint (19). Constraints (21) and (22) enforce the same log for other periods ($t\in T\backslash t={t}_{1}$) without transshipment decision variables between health facilities.

Note that prioritization constraints are in the form of a non-linear model with respect to the decision variables. Thus, for linearization, Constraints (19)–(22) can be replaced by the following Constraints (23)–(34):

$${e}_{vct}\le {\tilde{d }}_{vct} \forall v\in V, c\in C, t\in T$$

(23)

$${e}_{vct}^{^{\prime}}\le {\tilde{d }}_{vct}^{^{\prime}} \forall v\in V, c\in C, t\in T$$

(24)

$${e}_{vc{t}_{1}}\le {I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}} \forall v\in V, c\in C, t\in T$$

(25)

$${e}_{vc{t}_{1}}^{^{\prime}}\le {I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}} \forall v\in V, c\in C$$

(26)

$${\tilde{d }}_{vct}-{e}_{vct}\le M {g}_{vct} \forall v\in V, c\in C, t\in T$$

(27)

$${\tilde{d }}_{vct}^{^{\prime}}-{e}_{vct}^{^{\prime}}\le M {g}_{vct}^{^{\prime}} \forall v\in V, c\in C, t\in T$$

(28)

$${I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}\le M\left(1-{g}_{vc{t}_{1}}\right) \forall v\in V, c\in C$$

(29)

$${I}_{vc{t}_{1}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}-{e}_{vc{t}_{1}}^{^{\prime}}\le M\left(1-{g}_{vc{t}_{1}}^{^{\prime}}\right) \forall v\in V, c\in C$$

(30)

$${e}_{vct}\le {I}_{vct} \forall v\in V, c\in C, t\in T,t\ne 1$$

(31)

$${e}_{vct}^{^{\prime}}\le {I}_{vct}-{e}_{vct} \forall v\in V, c\in C, t\in T,t\ne 1$$

(32)

$${I}_{vct}-{e}_{vct}\le M\left(1-{g}_{vct}\right) \forall v\in V, c\in C, t\in T,t\ne 1$$

(33)

$${I}_{vct}-{e}_{vct}-{e}_{vct}^{^{\prime}}\le M\left(1-{g}_{vct}^{^{\prime}}\right) \forall v\in V, c\in C, t\in T,t\ne 1$$

(34)

Constraints (23) and (24) restrict the injected vaccine units to high-risk and low-risk groups of people, respectively, in each health facility that should be less than the demand of each group. Constraints (25) and (26) restrict the injected vaccine units to high-risk and low-risk groups of people in each health facility to the vaccine units in inventory and transshipment. According to Constraints (11) and (12), we define the shortage variables for high-risk and low-risk groups by ${k}_{vct}$ and ${k}_{vct}^{^{\prime}}$ that is equal to ${\tilde{d }}_{vct}-{e}_{vct}$ and ${\tilde{d }}_{vct}^{^{\prime}}-{e}_{vct}^{^{\prime}}$, respectively. Therefore, Constraints (27) and (28) indicate that shortages for both groups are positive values (or zero). Constraint (29) guarantees that the high-risk group’s demand is satisfied by the available vaccine units in each health facility at the first period. Constraint (30) illustrates that the low-risk group’s demand will be fulfilled with the remains of the available vaccine units. Constraint (31) restricts the injected vaccine units to the high-risk group of people in each health facility to its inventory at the beginning of each period. Constraint (32) also restricts the injected vaccine units to the low-risk group of people in each health facility to the remainder vaccine units. The injected vaccine units to high-risk groups are less or equal to inventory in each health facility shown in Constraint (33). Constraint (34) limits the number of injected vaccine units to the low-risk group of people.

Equity constraints:

In this pandemic, equity in the distribution of the COVID-19 vaccine is one of the main goals of our research. This assumption is well-considered in the following constraints, which tend to hold the minimum level vaccine distribution for all regions based on their total demands.

$$\frac{\sum_{l}\sum_{t}{x}_{vlht}^{^{\prime}}}{\sum_{c\in r}\sum_{t}\left({\tilde{d }}_{vct}+{\tilde{d }}_{vct}^{^{\prime}}\right) {u}_{hct}}\ge {\alpha }^{min} \forall v\in V, h\in H, r\in R$$

(35)

Constraint (35) ensures that the transferred vaccine units to a given hub distribution center must be at least a percentage of health facilities’ demands in the district assigned to that hub distribution center. The numerators represent the number of received vaccines to each hub distribution center in all periods. Simultaneously, the denominators reflect the total demands of high and low-risk groups in the health facilities of a given district allocated to a given hub distribution center during all periods.

4 Methodology

Uncertainty is an essential factor that affects the design of the supply chain network. The demand for each type of vaccine for the high-risk and low-risk groups in respect to the high volatility during the COVID-19 outbreak is considered an imprecise parameter. This section presents the MSSP method for solving the stochastic problem, and this method requires a decision to be made facing an uncertain future. After a while, the uncertainty is resolved, and it is possible to take any action to respond the realization of stochastic parameters. The MSSP method involves decision-making across all periods and therefore is an appropriate method for the proposed model. The applied methodology is explained in Sect. 4.1. Next, scenario generation and scenario reduction for demand scenarios are presented in Sects. 4.2 and 4.3. Figure 2 shows the algorithm of the MSSP approach.

Fig. 2

The algorithm of the MSSP approach

4.1 Multi-stage stochastic programming

In our problem, uncertainty and time play a key role, and therefore a decision model must be proposed to let the user employ a decision strategy that can reply to events as they display. Hence, a MSSP approach (which was first pointed out by Dupačová 1995) is presented to lecture a multi-period optimization model with uncertain demand. If each period and scenario are denoted by $t\in T$ and $s\in S$, respectively, the stochastic demand parameter indicates by ${\mathbb {D}}_{t}^{s}$ at period $t$ under scenario $s$. The relation between periods and stages of the MSSP approach is $\left|T\right|={\mathbb {N}}-1$. In the MSSP, much emphasis is disposed on the decision to be made in the current day, given current sources, later uncertainties, and probable future actions. A scenario tree presents uncertainty in stochastic programming, and it can be addressed in two forms, node-based and scenario-based. The compact form of the MSSP approach is as follows:

$$Min \sum_{s=1}^{\left|\mathrm{S}\right|}{p}_{s}f({\mathcal{x}}_{1}^{s},{\mathcal{x}}_{2}^{s},\dots ,{\mathcal{x}}_{\mathbb{N}}^{s};{\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})$$

(36)

Subject to:

$${\mathcal{x}}_{\mathbbm {n}}^{s}\in {\mathcal{X}}_{\mathbbm {n}}^{s} \forall {\mathbbm {n}}\in {\mathbb{N}}, s\in S$$

(37)

$${\mathcal{x}}_{{\mathbbm {n}}_{1}}^{s}={\mathcal{x}}_{{\mathbbm {n}}_{1}}^{{s}^{^{\prime}}} \forall \left(s,{s}^{^{\prime}}\right)\in S$$

(38)

$${\mathcal{x}}_{\mathbbm {n}}^{s}={\mathcal{x}}_{\mathbbm {n}}^{{s}^{^{\prime}}} \forall {\mathbbm {n}}\in {\mathbb{N}}\backslash \{{\mathbbm {n}}={\mathbbm {n}}_{1}\},(s,{s}^{^{\prime}})\in S :({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(39)

Eqs. (38) and (39) are the non-anticipativity constraints in stage 1 (${\mathbbm {n}}_{1}$) and other stages (${\mathbbm {n}}\in {\mathbb {N}}\backslash \{{\mathbbm {n}}={\mathbbm {n}}_{1}\}$), respectively, in the MSSP model. Where {${\omega }_{1},{\omega }_{2},\dots ,{\omega }_{{\mathbb {N}}-1}$} is a random sequence of parameters in stages and {${\omega }_{1}^{\mathrm{s}},{\omega }_{2}^{\mathrm{s}},\dots ,{\omega }_{{\mathbb {N}}-1}^{\mathrm{s}}$} is a realization of the stochastic parameter under scenario $s$ (${\omega }_{\mathbbm {n}}^{s}=({\mathbb {D}}_{t}^{s})$ when ${\mathbbm {n}}=t$). The probability of each scenario $s$ is ${p}_{\mathrm{s}}$ ${(\mathcal{x}}_{1}^{s},{\mathcal{x}}_{2}^{s},\dots ,{\mathcal{x}}_{\mathbb {N}}^{s})$ are decision vectors in stages $\{\mathrm{1,2},\dots , {\mathbb {N}}\}$ under scenario $\mathrm{s},$ and it is required decisions in each stage that have to be made before realization of stochastic parameters. Also, ${\mathcal{X}}_{\mathbbm {n}}^{\mathrm{s}}$ is the feasible region for decisions in stage ${\mathbbm {n}}$ under scenario $s$. For example, a usual scenario tree is illustrated in Fig. 3, a node-based form of a three-stage stochastic program. This scenario tree has ten nodes {${n}_{1},\dots , {n}_{10}$}, ${\mathbb{N}}=3$ stages {${\mathbbm {n}}_{1},{\mathbbm {n}}_{2},{\mathbbm {n}}_{3}$}, $\left|T\right|$=2 periods {${t}_{1},{t}_{2}$} and six scenarios {${s}_{1},\dots , {s}_{6}$}.

Fig. 3

An illustration of a scenario tree

The non-anticipativity constraints for the proposed model (see Sect. 3) are presented as follow:

$${u}_{lh{t}_{1}}{{^{\prime s}}}={u}_{lh{t}_{1}}{{^{\prime {s}}}^{^{\prime}}} \forall l,h, (s,{s}^{^{\prime}})\in S$$

(40)

$${u}_{hc{t}_{1}}^{s}={u}_{hc{t}_{1}}^{{s}^{^{\prime}}} \forall h,c, (s,{s}^{^{\prime}})\in S$$

(41)

$${u}_{lht}{{^{\prime s}}}={u}_{lht}{{^{\prime {s}}}^{^{\prime}}} \forall l,h,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(42)

$${u}_{hct}^{s}={u}_{hct}^{{s}^{^{\prime}}} \forall h,c,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S :({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(43)

The decision variables ${u}_{lht}{{^{\prime s}}}$, which determines that a hub distribution center is allocated to a national logistic center, and ${u}_{lht}^{s}$, which determines that a health facility is allocated to a hub distribution center, are binary variables. At each stage of MSSP, these allocation decisions should be made before identifying the uncertainty realization of stochastic parameters. Several studies considered the non-anticipativity equations for decision variables that have to be done at each stage after the uncertain parameters are exposed. Therefore, the non-anticipativity equations for several decision variables in the proposed model can be declared as follows:

$${x}_{vlht}{{^{\prime s}}}={x}_{vlht}{{^{\prime {s}}}^{^{\prime}}} \forall v,l,h,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(44)

$${x}_{vhct}^{s}={x}_{vhct}^{{s}^{^{\prime}}} \forall v,h,c,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S:({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(45)

$${{y}_{vct}^{s}= y}_{vct }^{{s}^{^{\prime}}}\forall v,l,h,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(46)

$${k}_{vct}^{s}={k}_{vct}^{{s}^{^{\prime}}} \forall v,h,c,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S:({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(47)

$${k}_{vct}{{^{\prime s}}}={k}_{vct}{{^{\prime {s}}}^{^{\prime}}} \forall v,l,h,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(48)

$${I}_{vct}^{s}={I}_{vct}^{{s}^{^{\prime}}} \forall v,h,c,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S:({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(49)

$${In}_{vct}^{s}={In}_{vct}^{{s}^{^{\prime}}} \forall v,l,h,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(50)

$${e}_{vct}^{s}={e}_{vct}^{{s}^{^{\prime}}} \forall v,h,c,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S:({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(51)

$${e}_{vct}{{^{\prime s}}}={e}_{vct}{{^{\prime {s}}}^{^{\prime}}} \forall v,l,h,t \backslash \{{t}_{1}\}, (s,{s}^{^{\prime}})\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(52)

Constraint (44) illustrates that the number of vaccine-type $v$ transferred from national logistic center $l$ to hub distribution center $h$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (45) states that the number of vaccine-type $v$ transferred from hub distribution center $h$ to health facility $c$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (46) represents that the order quantity of vaccine-type $v$ for health facility $c$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (47) assures that shortage of vaccine-type $v$ for the high-risk group in health facility $c$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (48) assures that shortage of vaccine-type $v$ for the low-risk group in health facility $c$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (49) guarantees that inventory level of vaccine-type $v$ in health facility $c$ at the beginning of period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (50) shows that inventory level of vaccine-type $v$ in health facility $c$ at the end of period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (51) determines that the number of vaccine-type $v$ fulfilled the demand of high-risk group in health facility $c$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. Constraint (52) guarantees that the number of vaccine-type $v$ fulfilled the demand of low-risk group in health facility $c$ at period $t \backslash \{{t}_{1}\}$ are equal for each pair scenarios $s, {s}^{^{\prime}}$. The node-based form of the MSSP method for the proposed model can be referred to in Appendix A.

4.2 Scenario generation for stochastic demand parameter

Generating a scenario tree is one of the most critical challenges in a multi-stage stochastic program that describes a multivariate stochastic method for the stochastic parameter. The scheme of scenario generation is going to explain in detail in this section. The discretization method for the continuous probability density function through a discrete scenario set is named scenario generation. The Monte Carlo approach is used to simulate the stochastic parameter for generating discrete scenarios. The demands of regions are assumed that follow a stochastic method. A modified first-order autoregressive model first introduced by Sodhi (2005) should be implemented to manage high time-variable demands. In this case, the demand for health facility $c$ under scenario $s$ at period $t$ can be calculated as follows:

$${\tilde{d }}_{vct}^{s}-{\overline{d}}_{vct}={\alpha }_{vc} \left({\tilde{d }}_{vc\left(t-1\right)}^{s}-{\overline{d}}_{vc\left(t-1\right)}\right)+{{\beta }_{vct}^{s}}^{N} \forall v,c,t \backslash \{{t}_{1}\}, s\in S: ({\omega }_{1}^{s},{\omega }_{2}^{s},\dots ,{\omega }_{{\mathbb{N}}-1}^{s})=({\omega }_{1}^{{s}^{^{\prime}}},{\omega }_{2}^{{s}^{^{\prime}}},\dots ,{\omega }_{{\mathbb{N}}-1}^{{s}^{^{\prime}}})$$

(53)

where ${\overline{d}}_{vct}$ is a time-depend pre-determined parameter applied to seize the seasonality of demands over the planning horizon, $\alpha$ is the autoregressive parameter and ${{\beta }_{vct}^{s}}^{N}$ is the error term with $N\left(0,{\sigma }^{2}\right)$. The next section will clarify how to reduce scenarios number and turn them toward a scenario tree.

4.3 Scenario reduction

The large number of scenarios achieved by the scenario generation method in Sect. 4.2 will result in an unsolvable optimization problem in a scenario-based stochastic program. Therefore, decreasing the number of scenarios is significant and is first presented by Dupačová et al. (2003). Suppose $\mathcal{M}$ is an $r$-dimensional stochastic distribution function $\omega ={\left\{{\omega }_{t}\right\}}_{t=1}^{\left|T\right|}$. When $\mathcal{M}$ has limited support, discrete scenarios such as $supp(\mathcal{M})=\left\{{\omega }^{1},{\omega }^{2}, \dots ,{\omega }^{\left|s\right|}\right\}, {\omega }^{s}={\left\{{\omega }_{t}^{s}\right\}}_{t=1}^{\left|T\right|}, s=\mathrm{1,2},\dots , \left|S\right|$ can be expressed through $\left|S\right|$. Besides, the corresponding scenario probabilities are $m_{s}$ and ${\sum }_{s=1}^{\left|S\right|}m_{s}=1$. Furthermore, let $\widehat{\mathcal{M}}$ is a distribution function with the equal dimension of another stochastic method, $\widehat{\omega }={\left\{{\widehat{\omega }}_{t}\right\}}_{t=1}^{\left|T\right|}$. Then, $supp(\widehat{\mathcal{M}})=\left\{{\widehat{\omega }}^{1},{\widehat{\omega }}^{2}, \dots ,{\widehat{\omega }}^{\left|{s}^{^{\prime}}\right|}\right\}$ with $\left|{S}^{^{\prime}}\right|$ scenarios and the related probabilities ${\widehat{\mathcal{m}}}_{{s}^{^{\prime}}}$, ${s}^{^{\prime}}=1,\dots ,\left|{S}^{^{\prime}}\right|$ and ${\sum }_{{s}^{^{\prime}}=1}^{\left|{S}^{^{\prime}}\right|}{\widehat{\mathcal{m}}}_{{s}^{^{\prime}}}=1$. The distance of the Kantorovich approach (${K}^{n}$) between $\mathcal{M}$ and $\widehat{\mathcal{M}}$ is the optimal solution that is presented as follows:

$${K}^{\mathcal{n}}\left(\mathcal{M},\widehat{\mathcal{M}}\right)=inf\left\{\sum_{s=1}^{\left|S\right|}\sum_{{s}^{^{\prime}}=1}^{\left|{S}^{^{\prime}}\right|}{\delta }_{s{s}^{^{\prime}}} {dis}_{\left|T\right|}\left({\omega }^{s},{\widehat{\omega }}^{{s}^{^{\prime}}}\right) :{\delta }_{s{s}^{^{\prime}}}\ge 0, \sum_{s=1}^{\left|S\right|}{\delta }_{s{s}^{^{\prime}}}={\widehat{\mathcal{m}}}_{{s}^{^{\prime}}}, \sum_{{s}^{^{\prime}}=1}^{\left|{S}^{^{\prime}}\right|}{\delta }_{s{s}^{^{\prime}}}={\mathcal{m}}_{s}, s=1,\dots ,\left|S\right|, {s}^{^{\prime}}=1,\dots ,\left|{S}^{^{\prime}}\right|\right\}$$

(54)

When ${dis}_{t}\left({\omega }^{s},{\widehat{\omega }}^{{s}^{^{\prime}}}\right)$ is equal to $\sum_{\mathcal{n}=1}^{t}\Vert {\omega }_{\mathcal{n}}^{s}-{\widehat{\omega }}_{\mathcal{n}}^{{s}^{^{\prime}}}\Vert , t=1,\dots ,\left|T\right|$ and measures the distance between two scenarios at period $t$. In the following, if we consider $\widehat{\mathcal{M}}$ that is the distribution of the reduced probability of $\omega ={\left\{{\omega }_{t}\right\}}_{t=1}^{\left|T\right|}$, discrete support of $\widehat{\mathcal{M}}$ consists of scenarios ${\omega }^{{s}^{^{\prime}}}={\left\{{\omega }_{t}^{{s}^{^{\prime}}}\right\}}_{t=1}^{\left|T\right|}$; and ${s}^{^{\prime}}\in \left\{1,\dots ,\left|s\right|\right\}\backslash DS$ when $DS$ is the deleted scenarios set. For a constant set $DS\subset \left\{1,\dots ,\left|S\right|\right\}$, the minimum distance between $\mathcal{M}$ and $\widehat{\mathcal{M}}$ can be calculated as follows that presented by Dupačová et al. (2003):

$${K}^{\mathcal{n}}\left(\mathcal{M},\widehat{\mathcal{M}}\right)=\sum_{s\in DS}{\mathcal{m}}_{s}\underset{{s}^{^{\prime}}\in \mathit{DS}}{\underset{\prod}{\mathit{min}}} {dis}_{\left|T\right|}({\omega }^{s},{\omega }^{{s}^{^{\prime}}})$$

(55)

The probabilities ${\widehat{\mathcal{m}}}_{{s}^{^{\prime}}}$ for scenarios of $\widehat{\mathcal{M}}$, i.e., the protected scenarios ${\omega }^{{s}^{{\prime}}}$; ${s}^{{\prime}}\notin DS$, are provided by ${\widehat{\mathcal{m}}}_{{s}^{^{\prime}}}={\mathcal{m}}_{{s}^{{\prime}}}+{\sum }_{s\in {DS}_{{s}^{^{\prime}}}}{\mathcal{m}}_{s}$ when ${DS}_{{s}^{^{\prime}}}=\left\{s\in DS: {s}^{^{\prime}}={s}^{^{\prime}}(s)\right\}$, and ${s}^{^{\prime}}(s)$ is an index set selection of nearest scenarios to ${\omega }^{s}$, $\forall s\in DS$, i.e. ${s}^{^{\prime}}\left(s\right)\in$ arg min ${dis}_{\left|T\right|}\left({\omega }^{s},{\omega }^{{s}^{^{\prime}}}\right)$,$\forall s\in DS$. For optimizing selection set $DS$ with determined cardinality $\#DS$, the following reduction problem can be solved as follows:

$$Min\left\{\sum_{s\in DS}{\mathcal{m}}_{s}\underset{{s}^{{\prime}}\notin \mathit{DS}}{\underset{\prod}{\mathit{min}}}{dis}_{\left|T\right|}\left({\omega }^{s},{\omega }^{{s}^{^{\prime}}}\right):DS\subset \left\{1,\dots ,\left|S\right|\right\}, \#DS=\left|S\right|-NS\right\},$$

(56)

When $NS=\left|s\right|-\#DS$ is the protected scenario number after reduction. Moreover, the above relation is similar to a problem of set covering proven by Dupačová et al. (2003). Even for unusual cases when one scenario is deleted ($\#DS=1$) and one scenario is kept ($\#DS=\left|S\right|-1$), this problem can be solvable. Two heuristics named forward and backward reduction are presented by Dupačová et al. (2003). The one scenario optimal selection repeats recursively till the $NS$ scenario is reached in the selection algorithm of the forward scenario. Against the one scenario, optimal deletion repeats recursively till the deletion of the $\left|S\right|-NS$ scenario, in the reduction algorithm of backward scenario. In the following, the construction of the scenario tree is defined in detail.

The generated scenario fan is turned to a representative scenario tree by Heitsch and Romisch (2005) to reduce the scenario numbers. Each scenario $s=\mathrm{1,2},\dots ,\left|S\right|$ can be considered as ${\omega }^{s}=({\omega }_{0}^{s}, {\omega }_{1}^{s},\dots ,{\omega }_{\left|T\right|}^{s})$ with probability ${\mathcal{m}}_{s}$ for a scenario fan with probability distribution $\mathcal{M}$ linked to the multivariate scenarios. As all scenarios stay the same for the beginning node (${\omega }_{0}^{1}={\omega }_{0}^{2},=\dots ={\omega }_{0}^{\left|s\right|}={\omega }_{0}^{*}$); therefore, the whole number of nodes is $\left|S\right|\left|T\right|+1$. Also, the scenario tree with ${\mathcal{M}}_{\theta }$ probability distribution has an origin node ${\omega }_{0}^{*}$. Besides, the node numbers in the scenario tree are also smaller rather than the scenario fan and ${K}^{\mathcal{n}}(\mathcal{M},{\mathcal{M}}_{\theta })$. We used the forward scenario reduction procedure to perform for period $t=\mathrm{1,2},\dots ,\left|T\right|$ and scenarios serial clustering to turn the scenario fan into a scenario tree. ${\theta }_{t}$ is considered for completing the scenario reduction procedure under the status ${\sum }_{t=1}^{\left|T\right|}{\theta }_{t}\le \theta$ to become a scenario tree with ${K}^{\mathcal{n}}(\mathcal{M},{\mathcal{M}}_{\theta })\le \theta$, at period $t$. It means that a maximal reduction procedure is applied at period $t$ with $\sum_{s\in DS}m_{s}\underset{{s}^{{\prime}}\notin {DS}}{\underset{\prod}{\mathit{min}}}{dis}_{t}({\omega }^{s},{\omega }^{{s}^{{\prime}}})\le {\theta }_{t}$. Furthermore, the distance among every two scenarios is obtained through ${dis}_{t}\left({\omega }^{s},{\omega }^{{s}^{{\prime}}}\right)=\sum_{\mathcal{n}=1}^{t}\Vert {\omega }_{\mathcal{n}}^{s}-{\omega }_{\mathcal{n}}^{{s}^{{\prime}}}\Vert$ at period $t$. Usually, $\theta ={\theta }_{rel}{\theta }_{max}$ is applied to calculate $\theta$ (Heitsch and Romisch 2005). ${\theta }_{rel}$ is a constant rate between 0 and 1 that describes a measure for the reduction value in the first scenario fan and ${\theta }_{max}$ is the greatest distance between the probability distribution of the first scenario fan and one of its scenarios with a probability unit. For generating a scenario tree, ${\theta }_{t}$ is obtained at period $t$ as follows:

$$\theta _{t} = \frac{\theta }{{\left| T \right| + 1}}\left[ {\frac{1}{2} + \bar{\rho }\left( {1 - \frac{t}{{\left| T \right| + 1}}} \right)} \right],\quad t = 1, \ldots ,\left| T \right|$$

(57)

When $\overline{\rho }$ is a constant value between zero and one that is considered equal to one in this research. Therefore, ${\sum }_{t=1}^{\left|T\right|}{\theta }_{t}$ is equal to $\frac{\theta }{\left|T\right|+1}(\frac{\left|T\right|}{2}+\overline{\rho }\frac{\left|T\right|}{2})$ and on the other, ${\sum }_{t=1}^{\left|T\right|}{\theta }_{t}\le \theta$ is due to $\overline{\rho }\in \left[\mathrm{0,1}\right]$. As a result, ${\sum }_{t=1}^{\left|T\right|}{\theta }_{t}=\frac{\left|T\right|}{\left|T\right|+1}\theta$ because $\overline{\rho }$ is considered equal to one.

5 Case study

As the capital of Iran, Tehran is one of the most populated cities in which healthcare providers face a growing demand for different health services. A real case study for the southwestern region (region 10) of Tehran is applied to evaluate the proposed model’s performance and feasibility. Municipality of Tehran segmented region10 to 17 districts as showed in Fig. 4. Each district is considered a demand node. Also, 10 candidate locations for establishing hub distribution centers have been selected in region 10 (Fig. 5). A national logistic center is established in Tehran (out of region 10) that transfers vaccine units to hub distribution centers in all regions of Tehran (Fig. 5). The population of demand nodes is obtained from Mousazadeh et al. (2018) is reported in Table 2.

Fig. 4

The boundary of the 17 demand nodes in region 10 of Tehran (from Google Maps)

Fig. 5

The location of all centers and hospitals in region 10 of Tehran (from Google Maps)

Table 2 The population of each demand node in region 10 of Tehran (Mousazadeh et al. 2018)

Besides, geographic dispersion of the logistic center, 10 candidates hub distribution centers, 21 health centers, and 10 hospitals located in region 10 are shown in Fig. 5. Also, the name of hospitals in demand nodes is recorded in this figure. As can be seen, at least one health center is there in each region, and some demand nodes have more than one hospital. For example, three hospitals are located in demand node 1, while there is no hospital in demand node 3.

There is a compelling indication in researches to show that two of the critical relevant factors in quantifying people’s susceptibility are ‘age and medical history’. Therefore, in the case study, we identify two priority groups:

• 1) High-risk group: Infected patients to COVID-19 are quickly transmitting the virus to others, or care operators are in an environment with high transmission potential that has severe risk and those who are an elevated chance, compared to others, of receiving serve disease COVID-19. They include older people, people with special medical situations, healthcare operators, other care operators, including injury support operators and quarantine operators. Based on existing documents in the Statistical Center of Iran (SCI) and Iran Ministry of Health and Medical Education (IMHME, http://corona.behdasht.gov.ir), 15% of Tehran’s population is in the high-risk group. Therefore, the high-risk group demands in a given district are estimated at around 15% of that district’s population.

• 2) Low-risk group: This group includes the rest of the people. According to IMHME and SCI’s documents, 85% of Tehran’s population is low-risk. Therefore, the demands of the low-risk group in a given district are estimated at around 85% of its population.

In addition, due to various conditions for transporting and storing COVID-19 vaccines, the type of vaccine is considered as a set in the proposed model. Therefore, the traveling and storing costs are different for each type. In this case study, three types of vaccines, AstraZeneca, Sinopharm, and Barakat, are injected in region 10 of Tehran. According to Table 3, cold chain requirements for these three types are in the same conditions; hence, the traveling cost (2 US dollars) and storing cost (200 US dollars per 1 year) are considered equal based on the National Committee on Combating Coronavirus of Iran and www.temparmour.com.

Table 3 The injected vaccine-types in region 10 of Tehran (www.icao.int)

As mentioned earlier, we convert a scenario fan of the demand parameter to a scenario tree. Therefore, we generate a scenario fan containing 170 scenarios by the Monte Carlo simulation method and according to Eq. (53) in Sect. 4.2. Then, we reduce the number of scenarios by the scenario reduction method in Sect. 4.3 and convert them to the scenario tree. Finally, eight scenarios remain in the scenario tree, as shown in Fig. 6. The constant rate (${\theta }_{rel}$) is considered 0.65 in implementing the reduction method. The probabilities of 22 nodes in this scenario tree are estimated based on available data obtained from health facilities of Tehran’s region 10 and the COVID-19 experts’ opinions (Table 4).

Fig. 6

The scenario tree of the case study

Table 4 The probabilities of the scenario tree’ nodes

The fixed costs are considered for establishing hub distribution centers to calculate them with other costs. Therefore, the values of ${cw}_{h}$ for a week are set as [535.25, 561.52, 521.20, 596.83, 446.34, 554.53, 541.05, 363.80, 375.22, and 460.45 US dollars] for hub distribution center 1 to 10, respectively, based on the opinions of experts and the reports of the Iran Ministry of Roads and Urban Development (https://tn.ai/2109911). Furthermore, the transporting cost between two nodes is calculated by multiplying the distance between them (via the roadway network in Google Maps) in the traveling cost per kilometer for a vaccine unit in 2 for the vehicle’s round trip. The demand is considered an uncertain parameter based on scenarios with underlying conditions estimated by IMHME reports. Due to the variety popularity of each vaccine type between people and its risk of shortage, the shortage cost for each type and group of people is different. Accordingly, the shortage penalty costs are assumed for high-risk: 35, 30, 25 and low-risk: 20, 15, 10 US dollars for AstraZeneca, Sinopharm, and Barakat, respectively, based on experts’ opinion in IMHME. In the following sections, the proposed model is implemented with GAMS 23.5.2 by CPLEX solver, and all sensitivity analyses are solved on a Core i5 2.70 GHz system with RAM 8 GB.

5.1 The result of the case study

In this section, we report the general results of the case study. The total cost of the vaccine supply chain network is 283172 US dollars associated with the establishment, transportation, transshipment, shortage, and storing costs. Usually, the high demand for each vaccine type in each period faces a shortage. In this model, the demand is considered as an uncertain parameter based on the scenario. Moreover, the geographic dispersion of the assignments and established hub distribution centers is provided in Fig. 7. As shown in Fig. 7a, two hub distribution centers 1 and 3 are opened throughout region 10 of Tehran. As can be seen, the scattering of the established hub distribution centers is gathered in the focal districts of region 10, where almost the districts’ population is higher compared to other districts. We can say that the number of demands for vaccine units is high around these centers, and the gathering of hub distribution centers in this area is because of this. It should be mentioned that hub distribution centers are not opened in some districts. Nonetheless, the lack of hub distribution centers in a district does not determine that the vaccine unit can not be distributed to that district. The logic is that the distribution services can be provided by established hub distribution centers in neighboring districts regarding their coverage space. The allocation pattern of each hospital and health center to hub distribution centers is provided in Fig. 7b and Fig. 7c, respectively. The number of allocated health centers and hospitals to each established hub distribution center is shown in Table 5.

Fig. 7

Assignment of (a) hub distribution centers to logistic center, (b) hospitals to hub distribution centers, and (c) health centers to hub distribution centers

Table 5 The number of allocated health facilities to the opened hub distribution centers

The summaries of the numerical results of the case study are presented in Tables 6, 7, 8, 9, and 10. Table 6 shows the number of transferred vaccine types distributed in hub distribution centers under each scenario. As inserted in this table, an area with an average of around 50,000 people is covered by each hub center during the time horizon under each scenario. Table 7 illustrates the number of transferred vaccine units to hospitals and health centers during the time horizon under each scenario. According to the results in this table, an average of around 31,000 vaccine units in 10 hospitals and 69,000 vaccine units in 21 health centers under each scenario will be injected into people during the time horizon, which will prevent a multitude of people from being affected COVID-19 infection. Notably, the number of injected vaccine units in each health center and each hospital is different. In other words, the capacities of supply are used as much as possible to avoid allocating health facilities to far centers.

Table 6 Number of different vaccine types transferred to opened hub distribution centers under each scenario

Table 7 Number of different vaccine types transferred to health facilities under each scenario

Table 8 Number of different vaccine types injected to two groups under each scenario

Table 9 The number of injected vaccine units to the high-risk group of people in each demand node under scenario 1

Table 10 Number of injected vaccine units to the low-risk group of people in each demand node under scenario 1

The vaccine units in health centers and hospitals are injected to meet the high-risk group of people, and the remaining are used for people with low-risk, and the remaining are stored or transshipped to other health facilities. Consequently, the number of demands for vaccine units in the next periods affects the number of required vaccine units. Table 8 reports the number of injected vaccine units to two groups of people under each scenario.

The number of injected vaccine units to a high-risk group of people in each demand node (district) in each period is presented in Table 9. Due to the one-day lead time, vaccine units are not injected in the first period. As can be seen in this table, the district with more demand has more injected units; for example, 2905 units of vaccine are injected in demand node 1. This table shows that demand nodes 1 and 7 have the highest and lowest ratio for vaccine injection, respectively, to the high-risk group in region 10. The remainders of the injected vaccine units are into the low-risk group of people in each demand node in each period under scenario 1 reported in Table 10. As mentioned, vaccine units are not injected in the first period due to the one-day lead time. This table shows that demand nodes 1 and 9 have the highest and lowest ratio for vaccine injection to the low-risk group in region 10.

In this paper, the proposed model is solved via the CPLEX solver in GAMS 23.5.2 and motivated by a real case on region 10 of Teheran. The model runs in reasonable times (00:02:14) with zero Gap, and sensitivity analyses are completed in less than 6 min. Hence, applying algorithms for solving the proposed model concerning the vaccine distribution network and data is not essential. The proposed model can be used for large-size case studies than region 10. Therefore, seven problem instances with various settings and sizes are generated to investigate the limit of the proposed model. The sets of these numerical examples, including the number of all nodes, national distribution centers, potential hub distribution centers, hospitals, health centers, and districts, are shown in Table 11. The computational time and Gap for solve each instance are also illustrated in Table 11. The maximum legal time is set to 10 h to solve the model that is not found any response for Instance 7 in less than 10 h. These results demonstrate that even though the computational time increases by expanding the size of problem occurrences, the small and medium sizes of the problem are solved by the presented model.

Table 11 Problem instances

5.2 Impact of equity in the distribution of vaccine units

It is more satisfactory to inject the available vaccines among all people fairly, even though it might not be enough for everyone at the current instance. In this study, equity among regions is modeled through a set of constraints that enforce a minimum supply level for each opened hub distribution center. Table 12 reports a comparison of the number of vaccine units transferred to opened hub distribution centers considering two states when the equity is assessed and when the equity is ignored in our model. The minimum percentage of total demand is considered 0.65. In other words, 65 percent of the total demand of each area should be fulfilled to obtain an equitable distribution network.

Table 12 The effect of equity policy on the number of transferred vaccine units to hub distribution centers under scenario

For example, the number of vaccine units transferred to the hub distribution center 1 in scenario 3 equals 53,006 by ignoring equity strategy, while this number is decreased to 52,178 by considering equity. This strategy has 828 units reduction, leading to increases in other hub centers by applying this strategy. Consequently, when equity is ignored, unfulfilled demand is much in some areas while others are minimal. The variation range of vaccine units in each hub distribution center under each scenario is calculated in the last line of Table 12.

In Fig. 8, the effect of the equity policy on vaccine units transferred to each opened hub distribution center is investigated. This figure shows that the number of vaccine units moved to hub 1 increases as the number of vaccine units transferred to hub center 3 decreases by considering equity strategy. In other words, the fairness and equity policy entail distributing vaccine units to areas based on their demand; consequently, the shortage of vaccine units may decrease in some areas.

Fig. 8

Comparison of transferred vaccine units to hub distribution centers with equity and without equity policy under each scenario

5.3 Impact of priority setting

The proposed model can prioritize sensitive people, whereby the priority group first receives vaccines from the health facilities. Then, the remainders of vaccines are injected into the second priority group of people. In Table 13, a comparison between the vaccine units injected to high-risk and low-risk groups of people is considered in two states: with priority and without priority setting. This table shows that the number of assigned vaccine units to the high-risk group grows by using the priority strategy. Consequently, the number of injections to the low-risk group of people decreases. For example, without the priority policy, the number of vaccine units injected into the high-risk group is 17679 and into the low-risk group is 81964 in scenario 4. In contrast, by considering priority setting, these numbers have 153 unit changes. These judgments are shown in more detail in Fig. 9.

Table 13 Number of injected vaccine units to groups of people considering their priority under each scenario

Fig. 9

Comparison of injected vaccine units to the high-risk and low-risk group of people concerning priority policy under each scenario (s)

5.4 Performance of MSSP approach

In this section, we report the result of the proposed model according to MSSP and compare it with the two-stage stochastic programming (TSSP). Huang and Ahmed (2009) presented a standard, called the relative value of MSSP (RVM), for comparing the importance of MSSP and TSSP models in a stochastic optimization problem as follows:

$$RVM=\left(\frac{{OF}^{TSSP}-{OF}^{MSSP}}{{OF}^{TSSP}}\right)\times 100$$

(58)

where ${OF}^{TSSP}$ and ${OF}^{MSSP}$ are the values of objective functions for the TSSP and MSSP models, respectively. In the case study for modeling the TSSP, the non-anticipativity constraints presented in Sect. 4.1 should be omitted from the model, and constraint (59) is added. In the TSSP, the location decisions of hub distribution centers, allocation decisions of opened hub distribution centers to logistic centers, and health facilities to hub distribution centers for each period should be decided at the beginning of the plan. In other words, these decisions are not dependent on uncertain parameters.

$${\mathcal{x}}^{s}={\mathcal{x}}^{{s}^{^{\prime}}} \forall (s,{s}^{^{\prime}})\in S$$

(59)

The optimal objective values of the TSSP and MSSP models are reported for the case study in Table 14. Moreover, the RVM value is presented according to Eq. (58) and is equal to 1.059. Table 14 shows that the objective function value is reduced meaningfully in the MSSP model compared with the TSSP model. Therefore, it can mean that we obtain a significant improvement by handling the MSSP model.

Table 14 Comparison between TSSP and MSSP approach

6 Conclusion

In the following, the summary of this paper, the main findings, the research limitations, and suggestions for future study are presented.

6.1 Summary of this study

This study designs a vaccine supply chain to obtain an equitable and accessible network. To do so, a novel mathematical formulation is presented to optimize the distributed vaccine units to areas and the injected vaccine units into people with various priority levels. We model uncertainty in demand due to the dynamic nature of COVID-19 outbreak through the MSSP approach. The Monte Carlo simulation method generates the scenarios in this approach. Then, a forward scenario reduction method is applied to create a suitable scenario tree. A real case of Iran is applied to show the practicality and capability of the presented model and approach, but this study’s result is not confined to just this case. The main findings of the presented model are stated in the following.

6.2 Main findings of this study

In this paper, the managerial insights have resulted as follows:

• 1) The multi-stage stochastic programming performance is significantly improved compared with the two-stage approach concerning the total cost of the vaccine supply chain network.

• 2) Lateral transshipment makes the vaccine supply chain flexible to the shortage. This policy uses inventory in health facilities to respond to shortages in other health facilities; either shortage is due to lack of enough supply of vaccine units or uncertain demand.

• 3) Distributing vaccine units through the national logistic center’s equity strategy can be applied to distribute vaccine units based on the demands of each district.

• 4) Considering the priority policy for different groups is useful and leads to timely vaccination of priority groups in the community.

6.3 Suggestions for future study

Although the proposed model for the vaccine supply chain network is examined for region 10 of Tehran, this research has several restrictions that could be fixed in future studies. To minimizing the number of unvaccinated people, mobile centers could be considered. Mobile centers could be assigned to hub distribution centers and reduced the risk of unfulfilled demand in regions. Also, uncertainty in other parameters could be applied in the model. Other health objective functions could also be combined alongside. Regarding the different maximum shelf-life of each vaccine type, we can assume that some of the vaccine units at health facilities are not injected during their shelf-life. Hence, those vaccines must be discarded and added the cost of expiry to the objective function.

On the other hand, the proposed model is complex and cannot be solved in a reasonable computational time for large-size instances. To do so, the number of continuous variables, binary variables, and constraints is increased. This study proposed a MILP model in which the number of binary variables is the main indicator for the complexity of the model (Alemany et al. 2018; Viana and Pedroso 2013). In other words, when the size of the network is increased, the model complexity and the number of binary variables increase, presenting a meta-heuristics or heuristics algorithm could be very significant.

Appendix A

See Table 15.

Table 15 Nomenclature of the multi-stage stochastic model

$$Min TC= \sum_{h}{cw}_{h}{w}_{h}+\sum_{nt\left(n,t\right)}{p}^{n}\left(\sum_{v}\sum_{l}\sum_{h}2\times{cd}_{vlh}^{^{\prime}}{f}_{lh}{u}_{lht}{{^{\prime n}}}+\sum_{v}\sum_{h}\sum_{c}2\times{cd}_{vhc}{f}_{hc}{u}_{hct}^{n}\right)+\sum_{v}\sum_{{t}_{1}}\sum_{c}\sum_{{c}^{^{\prime}}\ne c}2\times{cz}_{vc{c}^{^{\prime}}}{f}_{c{c}^{^{\prime}}}{z}_{vc{c}^{^{\prime}}t}+\sum_{nt\left(n,t\right)}{p}^{n}\left(\sum_{v}\sum_{c}{ck}_{vc}{k}_{vct}^{n}+\sum_{v}\sum_{c}{ck}_{vc}^{^{\prime}}{k}_{vct}{{^{\prime n}}}+\sum_{v}\sum_{c}{cs}_{vc} {In}_{vct}^{n}\right)$$

(A1)

$$\sum_{l}{u}_{lht}{{^{\prime n}}}\le 1 \forall h\in H ,nt({n}_{1},t)$$

(A2)

$$\sum_{h}{u}_{hct}^{n}\le 1 \forall c\in C,nt({n}_{1},t)$$

(A3)

$${u}_{hct}^{n}\le {w}_{h} \forall h\in H, c\in C,nt({n}_{1},t)$$

(A4)

$${x}_{vlht}{{^{\prime n}}}\le M{ u}_{lht}{{^{\prime {n}_{1}}}} \forall v\in V,l\in L,h\in H,nt(n,t)\backslash {n}_{1}$$

(A5)

$${x}_{vhct}^{n}\le M {u}_{hct}^{{n}_{1}} \forall v\in V,h\in H,c\in C,nt(n,t)\backslash {n}_{1}$$

(A6)

$${I}_{vc{t}_{1}}^{n}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}={e}_{vc{t}_{1}}^{n}+{e}_{vc{t}_{1}}{{^{\prime n}}}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}+{In}_{vc{t}_{1}}^{n} \forall v\in V,c\in C,nt({n}_{1},{t}_{1})$$

(A7)

$${I}_{vct}^{n}={e}_{vct}^{n}+{e}_{vct}{{^{\prime n}}}+{In}_{vct}^{n} \forall v\in V,c\in C,nt\left(n,t\right)\backslash {t}_{1}$$

(A8)

$${I}_{vc{t}_{1}}^{n}={q}_{cv{t}_{1}} \forall v\in V,c\in C,nt(n,{t}_{1})$$

(A9)

$${e}_{vct}^{n}+{k}_{vct}^{n}={d}_{vct}^{n} \forall v\in V,c\in C,nt\left(n,t\right)$$

(A10)

$${e}_{vct}{{^{\prime n}}}+{k}_{vct}{{^{\prime n}}}={d}_{vct}{{^{\prime n}}} \forall v\in V,c\in C, nt\left(n,t\right)$$

(A11)

$${y}_{vct}^{n}=\sum_{h}{x}_{vhc(t+1)}^{n} \forall v\in V,c\in C,nt\left(n,t\right)\backslash t=\left|T\right|$$

(A12)

$${In}_{vct}^{n}+\sum_{h}{x}_{vhc(t+1)}^{n}={I}_{vc(t+1)}^{n} \forall v\in V,c\in C,nt\left(n,t\right)\backslash t=\left|T\right|$$

(A13)

$${To}_{vc}-{I}_{vct}^{n}={y}_{vct}^{n} \forall v\in V,c\in C,nt\left(n,t\right)\backslash {t}_{1}$$

(A14)

$$\sum_{l}{x}_{vlht}{{^{\prime n}}}=\sum_{c}{x}_{vhct}^{n} \forall v\in V,h\in H, nt\left(n,t\right)\backslash {n}_{1}$$

(A15)

$$\sum_{h}\sum_{nt\left(n,t\right)}{x}_{vlht}{{^{\prime n}}}\le {S}_{vl}^{^{\prime}} \forall v\in V,l\in L,n\in N$$

(A16)

$$\sum_{c}{x}_{vlht}{{^{\prime n}}}\le {S}_{vh} \forall v\in V,h\in H, nt\left(n,t\right)$$

(A17)

$${In}_{vct}^{n}\le {S}_{vc}^{\prime\prime} \forall v\in V,c\in C, nt\left(n,t\right)$$

(A18)

$${e}_{vct}^{n}\le {d}_{vct}^{n} \forall v\in V,c\in C, nt\left(n,t\right)$$

(A19)

$${e}_{vc{t}_{1}}^{n}\le {I}_{vc{t}_{1}}^{n}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}} \forall v\in V,c\in C, nt\left(n,{t}_{1}\right)$$

(A20)

$${d}_{vct}^{n}-{e}_{vct}^{n}\le M {g}_{vct}^{n} \forall v\in V,c\in C, nt\left(n,t\right)$$

(A21)

$${I}_{vc{t}_{1}}^{n}+\sum_{{c}^{^{\prime}}\ne c}{z}_{v{c}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}^{n}\le M\left(1-{g}_{vc{t}_{1}}^{n}\right) \forall v\in V,c\in C, nt\left(n,{t}_{1}\right)$$

(A22)

$${e}_{vct}^{n}\le {I}_{vct}^{n} \forall v\in V,c\in C, nt\left(n,t\right)\backslash {t}_{1}$$

(A23)

$${I}_{vct}^{n}-{e}_{vct}^{n}\le M\left(1-{g}_{vct}^{n}\right) \forall v\in V,c\in C,nt\left(n,t\right)\backslash {t}_{1}$$

(A24)

$${e}_{vct}{{^{\prime n}}}\le {d}_{vct}{{^{\prime n}}} \forall v\in V,c\in C,nt\left(n,t\right)$$

(A25)

$${e}_{vc{t}_{1}}{{^{\prime n}}}\le {I}_{vc{t}_{1}}^{n}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}^{n} \forall v\in V,c\in C, nt\left(n,{t}_{1}\right)$$

(A26)

$${d}_{vct}{{^{\prime n}}}-{e}_{vct}{{^{\prime n}}}\le M {g}_{vct}{{^{\prime n}}} \forall v\in V,c\in C,nt\left(n,t\right)$$

(A27)

$${I}_{vc{t}_{1}}^{n}+\sum_{{c}^{^{\prime}}\ne c}{z}_{{vc}^{^{\prime}}c{t}_{1}}-\sum_{{c}^{^{\prime}}\ne c}{z}_{{vcc}^{^{\prime}}{t}_{1}}-{e}_{vc{t}_{1}}^{n}-{e}_{vc{t}_{1}}{{^{\prime n}}}\le M\left(1-{g}_{vc{t}_{1}}{{^{\prime n}}}\right) \forall v\in V,c\in C, nt\left(n,{t}_{1}\right)$$

(A28)

$${e}_{vct}{{^{\prime n}}}\le {I}_{vct}^{n}-{e}_{vct}^{n} \forall v\in V,c\in C,nt\left(n,t\right)\backslash {t}_{1}$$

(A29)

$${I}_{vct}^{n}-{e}_{vct}^{n}-{e}_{vct}{{^{\prime n}}}\le M\left(1-{g}_{vct}{{^{\prime n}}}\right) \forall v\in V,c\in C, nt\left(n,t\right)\backslash {t}_{1}$$

(A30)

$$\frac{\sum_{l}\sum_{t, nt\left(n,t\right)}{x}_{vlht}{{^{\prime n}}}}{\sum_{c\in r}\sum_{t, nt\left(n,t\right)}\left({d}_{vct}^{n}+{d}_{vct}{{^{\prime n}}}\right){u}_{hct}^{n}}\ge {\alpha }^{min} \forall v\in V,h\in H,r\in R, n\in N$$

(A31)