COVID-19 vaccine distribution planning using a congested queuing system-A real case from Australia.

Crisis-induced vaccine supply chain management has recently drawn attention to the importance of immediate responses to a crisis (e.g., the COVID-19 pandemic). This study develops a queuing model for a crisis-induced vaccine supply chain to ensure efficient coordination and distribution of different COVID-19 vaccine types to people with various levels of vulnerability. We define a utility function for queues to study the changes in arrival rates related to the inventory level of vaccines, the efficiency of vaccines, and a risk aversion coefficient for vaccinees. A multi-period queuing model considering congestion in the vaccination process is proposed to minimise two contradictory objectives: (i) the expected average wait time of vaccinees and (ii) the total investment in the holding and ordering of vaccines. To develop the bi-objective non-linear programming model, the goal attainment algorithm and the non-dominated sorting genetic algorithm (NSGA-II) are employed for small- to large-scale problems. Several solution repairs are also implemented in the classic NSGA-II algorithm to improve its efficiency. Four standard performance metrics are used to investigate the algorithm. The non-parametric Friedman and Wilcoxon signed-rank tests are applied on several numerical examples to ensure the privilege of the improved algorithm. The NSGA-II algorithm surveys an authentic case study in Australia, and several scenarios are created to provide insights for an efficient vaccination program.

Introduction

The COVID-19 pandemic that was first reported in December 2019, has significantly influenced nearly all healthcare systems worldwide. It has resulted in more than 250 million confirmed cases and more than five million reported deaths across the world (World Health Organization, 2021). The socio-economic impacts of the COVID-19 pandemic (e.g., self-isolation, travel bans, transportation restrictions, and border shutdowns) have led to economic crisis and recession in advanced economic systems (Nicola et al., 2020). The new variants of COVID-19 with severe transmissibility have also driven healthcare systems to face intense financial and infrastructural challenges in a short time (Pescaroli et al., 2021), exposing healthcare supply chains to the risk of failure. According to Altmann et al. (2020), these challenges have motivated policymakers to speed up the process of mass vaccination. In this vein, key stakeholders in healthcare systems, like public hospitals and state and federal governments, must apply innovative supply chain strategies in the mass vaccination of their citizens to efficiently reduce the risks of infection and death. A vaccine supply chain consists of complex processes (e.g. vaccine allocation and vaccine distribution) that require coordination between key stakeholders to help achieve timely delivery and distribution of vaccines (Brown et al., 2014). According to Alam et al. (2021), a vaccine supply requires a ”development phase and a fulfilment phase when it is approved by the healthcare authorities”. In other words, a vaccine must follow a more rigorous procedure than other products (e.g. textile, sport, automotive) to become commercialised and distributed. While many companies are able to accelerate production and distribution even in difficult situations (e.g. the supply of personal protective equipment during the COVID-19 pandemic), vaccine production and distribution in emergencies still requires maintaining sensible R&D activities, coordinating clinical trials, and authorising safe and transparent mechanisms that impact the entire supply chain network. It is, therefore, imperative for policymakers to understand how current healthcare systems must be redeveloped or upgraded to the so-called “crisis-induced supply chain” format to respond effectively and efficiently to the high vaccination demand with limited supply. We consider the crisis-induced supply chain a supply chain that aims to respond efficiently and effectively to emergencies (e.g., the COVID-19 pandemic). These chains are often organised in a temporary fashion but adopt similar supply chain strategies (e.g., responsive or efficient approaches) (Dasaklis et al., 2017). The emphasis of these supply chains is on immediate action and reaction to a crisis. Supply chain coordination in such situations appears to be more challenging than in typical situations, as the level of uncertainty, is high. Uncertainty in a crisis-induced supply chain is originated from the unfamiliarity of stakeholders with the crisis-induced scenarios (e.g., the COVID-19 pandemic) and its unexpected consequences (e.g., blockages across major supply chain networks, shortage in supply and trade wars). In addition, procurement, transportation and distribution in emergencies require a great deal of international collaboration between public and private sectors to ease the equality and accessibility in vaccination. In this vein, the COVID-19 vaccine distribution planning may be incompatible with other vaccine distribution networks in some circumstances (e.g., the lack of inventory accessibility) and the increased price of vaccines due to the R&D expenditures), necessitating the development of and planning for new resources to improve the supply chain performance.

In addition, the uneven geographic distribution of hospitals and their branches in metropolitan compared to suburban areas and the rate of vulnerability to COVID-19 among different age groups have caused serious challenges in terms of ensuring equity in vaccine distribution (Bettampadi et al., 2021, Lee et al., 2012) as well as managing service costs (e.g., vaccine ordering and holding costs) (Hovav and Tsadikovich, 2015). For instance, the ratio of Inner Melbourne’s population ( citizens) to Inner Melbourne’s public hospitals ( hospitals) is greater than the remaining population of Melbourne ( citizens) to the remaining public hospitals ( hospitals) (Australian Bureau of Statistics, 2021). Furthermore, more than 15% of Melbourne’s population include citizens aged 65 and above (Australian Bureau of Statistics, 2021), who are at the greatest risk of death from a COVID-19 infection (Daoust, 2020). According to Mallapaty (2020), more than 80% of COVID-19 deaths occur in people aged 65 and above. These statistics reveal the complexity of the COVID-19 vaccine supply. In this vein, an integrated approach must be developed and implemented with multiple supply scenarios to help adopt a crisis-induced supply chain network in the short term (Chowdhury et al., 2021, Jahani et al., 2018). On the one hand, when more hospitals have vaccine allocation and distribution, this increases the pace of vaccination and, eventually, citizen satisfaction. On the other hand, when too many hospitals are equipped for vaccination, this disperses the targeted population and significantly increases the vaccine supply chain costs (i.e., ordering and holding costs). In addition, the operational costs rely on the number of open hospitals, their locations, their distances from the main healthcare network centre, and the targeted population in a specific living area (Li et al., 2021). This situation may become even more complex when the crisis-induced vaccine supply chains must cope with varying types of vaccines. This is because vaccines of varying brands have different characteristics, which require different production, transportation, storage, and administration strategies. For example, while both Pfizer–Biontech (known as Pfizer) and AstraZeneca–Oxford (known as AstraZeneca) are COVID-19 vaccines with two doses that must be given in different time windows, Johnson & Johnson’s Janssen (known as J&J) vaccine is comprised of only one dose. The differences in transportation and storage, efficiency, and costs also vary depending on which vaccine has priority (Abbasi et al., 2020, Tavana et al., 2021). Therefore, it is imperative to develop distribution mechanisms to respond effectively and efficiently to such challenges and mitigate likely disparities.

Despite all efforts to develop and implement healthcare supply chain strategies for a crisis-induced vaccine supply chain (e.g. see Abbasi et al., 2020, Chen et al., 2020, Ferranna et al., 2021) in prior studies, these analyses did not fully consider the multifaceted aspects of vaccine distribution mechanisms to mitigate the risk of failure in crisis-induced vaccine supply chains. For example, despite some studies have only focused on equity in vaccine distribution, based on variables such as the level of vulnerability, accessibility and availability Abbasi et al., 2020, Chen et al., 2020, Ferranna et al., 2021, the important of other constraints such as wait time and product holding and ordering costs are often overlooked due to the complexity of planning systems. More importantly, vaccines are very sensitive to the temperature change, requiring a high-level coordination across the chain to prevent vaccine failure. This situation is even more critical in a crisis-induced environment where the production and distribution processes strongly impact human life. In addition, the crisis-induced vaccine supply chains are more expensive than typical vaccine supply chains, as crisis-induced environments may create unexpected circumstances (e.g. inequality in distribution, a dramatic increase in the vaccine price, increased rate of vaccine failure). According to Alam et al. (2021), these unexpected circumstances are rooted in challenges such as the lack of vaccine monitoring bodies, the difficulty in tracking vaccinated population, inappropriate coordination with local organisations, the lack of proper storage systems and failure in vaccine marketisation. In this vein, millions of COVID-19 vaccines are being distributed to hundreds of local public hospitals on a daily or weekly basis, requiring public hospitals to consider a wide variety of factors in prioritising vaccine demands among distinct groups of citizens (e.g., vulnerable communities) in different time windows based on the type of vaccine. Therefore, a robust distribution and storage planning model is required to conceptualise the effective and responsive vaccination process at the operational level. The model should also ensure transparency in vaccine distribution by minimising the waiting lists of vaccine demands and reducing the expenses of vaccine ordering and holding.

This study aims to highlight the role of queuing theory in managing vaccine supply chains in emerging situations by considering two opposing yet interlinked factors in the success of crisis-induced vaccine supply chains: (i) minimising the wait time of vaccination, which will increase the effectiveness and responsiveness of such supply chains and (ii) minimising the total ordering and holding costs of vaccines and vaccination. This study thus develops a mathematical model emphasising the crisis-induced vaccine supply in healthcare supply chains to ensure effective and efficient vaccine distribution and storage at the operational level. We refer to the ordering costs as fixed costs (e.g., vaccine price) and variable costs (e.g. transportation) while holding costs are a combination of fixed and variable costs (e.g., vaccine failure, equipment, repair, and maintenance) that hospitals (or local authorities) must spend for vaccination. This study focuses on the greater Melbourne area as a case-study. It tests the proposed model based on real data from public hospitals and demographic information of citizens of the greater Melbourne. The findings of this study can be generalised to all countries affected by the COVID-19 pandemic (and other potential pandemics in the future), as this study makes several recommendations on how policymakers and managers in public hospitals must work together to establish a vaccine supply chain in times of crisis.

This paper is structured as follows: Section 2 describes the crisis-induced vaccine supply chain concept and later presents prior studies that adopted mathematical modelling techniques to develop optimised vaccine distribution networks. This section then discusses the theoretical framework of this study by highlighting the key role of the queuing theory in developing COVID-19 vaccine distribution planning using a congested queuing system. Section 3 describes the research problem and defines assumptions, properties, and the supply chain optimisation model. We present the solution methodology in Section 4, describing the non-dominated sorting genetic algorithm (NSGA-II) and solutions to the model. This is followed by Section 5, which tests the proposed model in an authentic case study. Section 5 also conducts a sensitivity analysis and performance evaluation to determine how the model and the proposed solution approaches perform in different scenarios. Finally, Section 6 concludes the main findings and provides recommendations for future policymakers and healthcare managers in the current and potential crisis-induced events, followed by recommendations for future studies.

Literature review

Crisis-induced vaccine supply chain

A conventional vaccine supply chain includes a complex set of processes, equipment, logistics, and transportation that involves vaccine production, procurement and prompt distribution (Brown et al., 2014, Shittu et al., 2016). Traditionally, manufacturers, wholesalers, or distributors supplied vaccines to care service providers (e.g., hospitals, clinics, and pharmacies) in a standardised chain to minimise vaccination wait time and speed (Li et al., 2018). A good example is the distribution mechanism of the influenza vaccine that occurs seasonally and often follows routinised planning for ordering, storage, and distribution. However, crisis-induced vaccine supply chains are developed in emergencies and require immediate planning for vaccine distribution; these may be more fragile to disruptions (e.g., shortages in supply, storage, and access restrictions). Indeed, prior studies often focused on disruptions in everyday situations (Chen et al., 2020, Lee et al., 2016, MacIntyre et al., 2021), overlooking the importance of disruptions in emergencies. That is, crisis-induced vaccine supply chains can experience a wide variety of disruptions not only in the development and production phases, but also in distribution mechanisms where the vaccination process may be influenced by ambiguity in vaccine efficiency, volatility in demand and complexity in pricing (Abbasi et al., 2020, Chaturvedi and Chakravarty, 2021, Chen et al., 2021). Other studies also emphasised lead-time disruptions in vaccine logistics and transportation networks and called for further research to examine these disruptions in crisis-induced scenarios (Lin et al., 2020, Esmizadeh et al., 2021).

In addition, prior studies only adopted two protection tactics against disruptions: reactive and proactive tactics (Sazvar et al., 2021, Sinha et al., 2021). While the former is used to protect supply chains from potential disruptions, the latter reduces post-disruption losses (Marcello et al., 2014). These studies mentioned above addressed protection tactics in the context of inventory hedging, contingent sourcing, supplier diversification, and demand switching. The inventory hedging tactic is described as excess inventory carried by the nodes to hedge against network disruptions and supply chain risks. This is often not ideal for perishable products, as its adoption may lead to significant losses under a disruption-free scenario (Marcello et al., 2014, Jahani et al., 2017). Other contexts may not be ideal to crisis-induced vaccine supply chains for the following reasons (Seib et al., 2012):

Demand switching (reactive tactic) is not appropriate for crisis-induced vaccine supply chains since it contradicts the policymakers’ objectives to ensure COVID-19 herd immunity. Demand switching occurs when customers change their preferences based on their perception of a product’s characteristics (e.g., its quality), which in turn impacts decisions on the product supply (e.g., import, inventory, distribution) (Hsieh, 2011). In a crisis with a high-level supply–demand uncertainty, customers prefer to have access to a wide variety of product choices. In contrast, policymakers prefer to limit or edit the customers’ choices to manage high-demanding products better and facilitate inventory and distribution (Altmann et al., 2020). Therefore, although demand switching can increase customer satisfaction in normal situations, paying attention toward mitigating the risk of death (e.g. through herd immunity) rather than promoting demand switching (e.g. switching from AstraZeneca to Pfizer) may be preferred.

Working in co-operation with various suppliers ensures a crisis-induced supply chain strategy through supplier diversification (proactive tactic), although corresponding to tight supplier-related regulatory requirements may increase the complexity of the vaccine supply chain and vaccine expenses.

The contingency sourcing tactic (reactive tactic) that is often suitable for vaccine supply chains may not benefit a crisis-induced vaccine supply chain since it may result in quality assurance difficulties.

Therefore, this study adopts a set of multi-objective optimisation algorithms to acquire optimal solutions and optimise wait time, holding costs, and ordering costs simultaneously in a crisis-induced vaccine supply chain. To our knowledge, this has not been investigated in prior literature.

There is currently a vast and well-theorised body of knowledge on vaccine supply chains (See Table 1), most of which were motivated by cost-effectiveness and efficiency of vaccine procurement, allocation, and distribution to vaccinate a large population. Many studies applied simulation models to better respond to uncertainties and complexities arising from real-world problems in vaccine supply chains (i.e., vaccine production, allocation, and distribution) (Lee et al., 2011, Lee et al., 2012, Seib et al., 2012, Davila-Payan et al., 2014, Marcello et al., 2014, Kazaz et al., 2016, Shittu et al., 2016). In this vein, Brown et al. (2014) and Lee et al. (2016) applied event-driven simulation techniques for resource allocation to investigate how a newly developed vaccine can be incorporated into an existing vaccine supply chain. However, these two studies relied on the capacity of an existing supply chain, and elements like disruptions in production and distribution, especially in emergencies, were overlooked. Other techniques of the simulation were used for various interventions in inventory level, demand, and supply (Lee et al., 2012, Seib et al., 2012, Davila-Payan et al., 2014, Govindan et al., 2020) or predicting the impact of newly developed vaccines in supply chain performance (Lee et al., 2011, Davila-Payan et al., 2014, Marcello et al., 2014).

Table 1

A summary of recent studies on vaccine supply chains.

Author(s)	Research objectives	Features								COVID-19	Real case-study	Research methodologies
		Ordering/holding costs	Transportation costs	Lead time disruption	Priority	Congestion	Production disruption		Multi-period
Brandeau et al. (2003)	To determine the optimal allocation of limited resources for epidemic control among multiple segments of a country’s population.	–	✓	✓	–	–	–	–	✓	–	–	Mathematical model (differential equations)
Chick et al. (2008)	To determine the best scenario for distribution.	–	–	✓	–	–	–	✓	–	–	–	Mathematical model (game theory)
Arora et al. (2010)	To propose a resource allocation approach for optimising regional aid during public health emergencies.	–	✓	–	✓	–	–	–	–	–	USA	Mathematical model (quadratic optimisation)
Lee et al. (2011)	To develop a computational model for vaccine supply chains.	–	–	✓	–	–	✓	✓	✓	–	Thailand	Simulation methods
Abrahams and Ragsdale (2012)	To propose a decision support system for scheduling in travel vaccine administration.	✓	✓	–	–	–	–	✓	✓	–	–	Binary integer programming, genetic algorithm
Lee et al. (2012)	To propose an optimal age-specific vaccination strategy against pandemic influenza.	–	–	✓	✓	–	–	–	–	–	–	Simulation model
Seib et al. (2012)	To examine vaccine providers’ willingness to respond.	–	–	✓	–	–	✓	✓	✓	–	–	Statistical method (regression)
Brown et al. (2014)	To explore the impact of cost on vaccine availability.	–	✓	–	–	–	–	–	–	–	Benin	Computational model (HERMES)
Davila-Payan et al. (2014)	To examine the impact of delivery and system factors in H1N1 vaccination rate.	–	–	–	✓	–	✓	✓	✓	–	USA	Statistical method (regression)
Marcello et al. (2014)	To propose a distribution of pandemic influenza vaccine.	–	–	–	✓	–	–	–	✓	–	USA	Simulation model
Kazaz et al. (2016)	To develop a model of supply chain and investigate the impact of various interventions.	–	✓	✓	–	–	✓	–	–	–	–	Mathematical and simulation model
Lee et al. (2016)	To redesign the vaccine supply chain to improve supply chain efficiency.	–	✓	–	–	–	–	✓	✓	–	Mozambique	Computational model (HERMES)
Shittu et al. (2016)	To develop a simulation model to explore the effects of supply and demand on storage capacity requirements.	–	✓	–	–	–	–	–	✓	–	Nigeria	Simulation model
Li et al. (2018)	To examine vaccine inventory levels in vaccine distribution systems .	–	✓	–	–	–	–	✓	✓	–	–	Simulation method
Wichapa and Khokhajaikiat (2018)	To propose a hybrid multi-objective location problem model.	–	✓	–	–	–	–	–	✓	–	Thailand	Goal programming and genetic algorithm, fuzzy AHP
Abbasi et al. (2020)	To optimise vaccine allocation on the downstream part of a COVID-19 vaccine supply chain.	–	✓	✓	✓	–	–	–	–	✓	Australia	Mathematical model (mixed integer programming)
Chen et al. (2020)	To propose an allocation model for the COVID-19 vaccine distribution networks.	–	–	–	✓	–	✓	–	✓	✓	USA	Mathematical model (game theory)
Govindan et al. (2020)	To develop a practical decision support system for COVID-19 healthcare supply chain.	–	–	–	✓	–	✓	–	✓	✓	WHO	Mathematical model (fuzzy inference system)
Lin et al. (2020)	To develop a conceptual model for transportation decisions in a cold vaccine supply chain.	–	✓	✓	–	–	✓	–	–	–	China	Mathematical model (game theory)
Buhat et al. (2021)	To determine an optimal and equitable allocation of COVID-19 vaccines while minimising deaths and satisfying the priority groups for immediate vaccination.	–	–	–	–	✓	–	✓	–	✓	Philippines	Mathematical model (linear programming)
Chen et al. (2021)	To examine how a limited number of vaccine doses can be strategically distributed to individuals to reduce the overall burden of the pandemic.	–	–	–	–	✓	–	✓	–	✓	–	Mathematical model (SEIR model)
Chaturvedi and Chakravarty (2021)	To overcome the epidemic by considering vaccination rate effects on the dynamics of COVID-19 control.	–	–	–	–	–	–	–	–	✓	India–Brazil–USA	Simulation method mathematical model (SIR model)
Foy et al. (2021)	To assess age-specific vaccine allocation strategies.	–	–	–	–	✓	–	–	–	✓	India	Mathematical models (age-structured model+ expanded SEIR)
Jahn et al. (2021)	To provide evidence-based guidance to the authorities to minimise COVID-19-related hospitalisations and deaths.	–	–	–	–	✓	–	–	–	✓	Austria	Simulation method (agent-based simulation method)
MacIntyre et al. (2021)	To respond to limited supply (age-targeted or ring vaccination) and mass vaccination.	–	–	–	–	–	–	–	–	✓	Australia	Mathematical model (age-structured model)
Minoza et al. (2021)	To formulate a multi-objective linear programming model to optimise vaccine distribution.	–	–	–	–	–	–	✓	✓	✓	Philippines	Mathematical model (multi-objective linear programming)
Roy et al. (2021)	To incorporate epidemiological factors, like population density, susceptible count and infected ratio, and transportation costs, and disseminate vaccines among zones.	–	✓	–	–	–	–	✓	–	✓	USA	Mathematical model (linear programming)
Sazvar et al. (2021)	To understand the application of capacity planning in terms of redundancy and design a supply chain network that is resilient toward the demand side.	–	–	–	✓	–	✓	✓	✓	–	Iran	Mathematical model (robust fuzzy optimisation)
Shim (2021)	To determine optimal vaccine allocation for minimising infections, deaths, and years of life lost while accounting for population factors	–	–	–	–	–	–	–	–	✓	South Korea	Mathematical model (age-structured dynamic model)
Sinha et al. (2021)	To identify conditions under which a strategic inventory reserve policy cannot be practically implemented to meet service level targets	✓	✓	✓	–	–	–	–	✓	✓	India	Mathematical model (game theory)
Tavana et al. (2021)	To develop equitable COVID-19 vaccine distribution in developing countries.	✓	✓	–	–	–	–	–	✓	✓	India	Mathematical model (mixed-integer linear programming)
This study	To propose a vaccine distribution model for the COVID-19 pandemic	✓	✓	✓	✓	✓	✓	✓	✓	✓	Australia	Mathematical model (non-linear programming), goal attainment, NSGA-II

Over the past decade, vaccine supply chains have experienced a growing level of complexity and uncertainty in their production, distribution, and allocation systems, requiring supply chain practitioners and scholars to develop and implement hybrid vaccine production and distribution methods. For instance, Kazaz et al. (2016) developed a hybrid mathematical and simulation model to examine the impact of disruptions in a vaccine supply chain. Wichapa and Khokhajaikiat (2018) proposed a hybrid model of vaccine distribution using a goal programming technique, combined with the genetic algorithms technique and fuzzy AHP to locate the best spots for temporary hubs in a vaccine supply chain. Abrahams and Ragsdale (2012) applied a hybrid model of binary integer programming and genetic algorithm to design a decision support system for patient scheduling in travel vaccine administration. However, none of these studies tested their models in an authentic case study to deal with real issues in a vaccine supply chain (e.g., holding and ordering costs) for the healthcare ecosystem, especially in emergencies (e.g., the COVID-19 Pandemic).

Recently, many studies have applied mathematical modelling techniques to address the issue of vaccine distribution for immediate mass vaccination. MacIntyre et al. (2021) used a deterministic mathematical model of epidemic response with limited supply in Australia. They concluded that an age-targeted vaccination strategy plays a crucial role in mass vaccination where the capacity of supply is limited, although their results were concluded based on prior research. Some other scholars focused on modelling vaccine supply at the country level. Foy et al. (2021) applied a simulated susceptible–exposed–infectious–recovered (SEIR) mathematical model to evaluate various vaccine allocation strategies for different age groups in India. They found that speeding the vaccination process for aged people can significantly reduce the population’s mortality rate. They also recommended that future research consider the types of vaccines as a predominant characteristic of a mass vaccination distribution mechanism. In this vein, Moghadas et al. (2021) compared Pfizer–Biontech and Moderna vaccines and highlighted the significant and positive role of vaccine type in deploying a vaccine distribution mechanism. Using SEIR mathematical modelling, Ferranna et al. (2021) demonstrated that essential workers (e.g. caregivers, nurses and practitioners) must be prioritised over aged people for COVID-19 vaccination to help significantly reduce infection-fatality risk in the United States of America (USA). Chaturvedi and Chakravarty (2021) focused on an analytically simulated susceptible–infectious–recovered model to predict the value of each dose of a COVID-19 vaccine in India, the USA and Brazil, and identified parameters that can contribute to ending the pandemic. They achieved comparable results for all three countries, but they did not provide any recommendations to help healthcare ecosystems address problems associated with the socio-economic status of citizens (e.g., access to care, investment in vaccines).

Some other studies identified the challenges of a vaccine distribution network, implying the importance of queuing mechanism and wait time calculation among different age groups. However, many of these studies did not fully or partially consider cost factors a vital characteristic of an efficient vaccine distribution mechanism. For example, studies by Jahn et al. (2021) in Austria and Shim (2021) in South Korea suggested that the efficiency of a vaccine distribution mechanism through a dynamic aged-based population model is critical for limited vaccine supply. Neither of these two studies considered the wait time of vaccine delivery and total investment vaccine costs as two essential parameters of vaccine distribution to minimise the rate of hospitalisation and mortality while vaccination efforts are in progress. In contrast, Chen et al. (2021), who applied an agent-based modelling approach to vaccine allocation, found that “[the] individuals’ degree (number of social contacts) and total social proximity time is significantly more effective than the currently used age-based allocation strategy in terms of number of infections, hospitalisations, and deaths”.

Cost factors were also reported to be a challenging element of COVID-19 vaccine distribution networks. While some studies relating to COVID-19 vaccine supply management modelled transportation expenses in their analyses (e.g., see Abbasi et al., 2020, Roy et al., 2021, Tavana et al., 2021), Buhat et al. (2021) suggested future studies must consider real scenario-based transportation costs in their analyses to warrant accuracy in their modelling. Many studies reported that decisions about holding and ordering costs play a crucial role in analysing the efficiency of a vaccine supply chain (e.g., see Abrahams and Ragsdale, 2012, Brown et al., 2014, Abbasi et al., 2020, Buhat et al., 2021), as the proposed objective functions are based on controlling or minimising the related costs in a vaccine supply chain. However, the literature has paid less attention to the mutual effects of these two components in vaccine supply chains (e.g. see Abrahams and Ragsdale, 2012). One reason is that these two factors (wait time and ordering and holding costs) are contradictory and cannot be easily rationalised in a crisis-induced vaccine supply chain where vaccinee satisfaction, the minimisation of total costs, and the pace of vaccination must be considered together.

Our review of the most relevant studies indicates that only a few studies developed and implemented an efficient distribution network of vaccines in an emergency — one that ensured the vaccine would be effectively delivered to the neediest individuals while considering both holding and ordering costs. We also found that none of the previous studies considered the congestion in the vaccine supply chain that stemmed from real-world endogenous and exogenous parameters. Indeed, many studies applied the mathematical models in real case studies associated with scenarios like influenza and H1N1, and only some of these studies characterised multi-objective parameters in examining efficiency in a crisis-induced vaccine supply chain (e.g., see Brandeau et al., 2003, Wichapa and Khokhajaikiat, 2018, Minoza et al., 2021).

To bridge these theoretical gaps, this study develops a conceptual model of a crisis-induced vaccine supply chain deployed in an authentic case study to ensure efficient coordination and distribution of different COVID-19 vaccine types to people with varying vulnerability levels to the COVID-19 pandemic. In addition, we theorise a queuing model for handling congestion in vaccine distribution networks to understand better how to (i) manage the congestion in the process of mass vaccination and (ii) minimise the expected average wait time and total investment in holding and ordering costs of vaccines simultaneously.

Thus, the following contributions are presented:

Designing a vaccine distribution mechanism to reduce the total investment vaccine cost and minimise the average expected wait time of vaccinees in an authentic case study.

Modelling congestion in a real-world vaccine distribution network and associating the wait times of vaccinees with the efficiency of vaccines and a risk aversion coefficient.

Improving the NSGA-II algorithm by adding several solution repair steps to the classic NSGA-II algorithm to ensure the algorithm’s efficiency through several statistical tests.

Designing the new chromosome structure to determine the number of servers (i.e. nurses and practitioners) and the number of orders during a specific time planning horizon.

Theoretical framework

As explained previously, prior studies considered either wait time or vaccine expenses as constraints in the vaccine supply chain and distribution networks. The reason is that prior studies often overlooked the adversarial relationships between wait time and vaccine holding and ordering costs. This study thus adopts the queuing theory to address these two constraints simultaneously and help policymakers make appropriate decisions in a crisis-induced vaccine supply chain.

The queuing theory is the mathematical modelling of variables and constraints that contribute to waiting lines (Mohammadi et al., 2014). Based on this theory, a queuing system includes an arrival process of entities that need to be served (Singer and Donoso, 2008). In this regard, if the system only explains the status quo, the queuing system will be a descriptive model. In contrast, when it intends to facilitate decision-making about a given situation, the queuing system will be a prescriptive model (Li and Zhang, 2015). Like other complex systems, in healthcare ecosystems, prescriptive modelling helps optimise the course of actions that pursue efficiency in the system (Boulton et al., 2016). Emergency service provision has been a topic of interest for many studies that adopted the queuing theory (e.g. see Sayarshad et al., 2020, Geroliminis et al., 2009, Smith, 1991). However, the application of this theory in the vaccine supply chain, especially in a crisis, is not well articulated in prior research. One reason is that prior studies appear not to consider a vaccine supply chain an emergency. In other words, prior studies did not fully consider efficiency from the lens of crisis. It is noteworthy that the wait time of vaccination, type of vaccines, vaccinees’ demographic information, and a vaccine’s holding and ordering costs are simultaneously critical in emergencies.

Using the queuing theory, this study assumes that distribution networks in a crisis-induced vaccine supply chain perform in a queuing system in which the vaccinee’s characteristics, the vaccine’s characteristics, and the distribution network’s characteristics simultaneously play a vital role in enhancing efficiency in the supply chain. For example, vaccinees of different ages may express a specific interest in a particular vaccine type (e.g., Pfizer vs. AstraZeneca) and expect to receive the vaccination service during a specific time frame. At the same time, a vaccine distribution network (e.g., public hospitals) may have to prioritise vaccinees under specific conditions (e.g., vulnerability, accessibility and affordability). These constraints can complicate the supply chain, but adopting the queuing theory can help optimise queuing systems and identify the best scenarios even if the constraints play an opposing role in the mathematical model.

Problem description and mathematical model

Problem definition

In this study, we assume a conceptual multi-period model for the COVID-19 vaccine allocation. We develop four time windows for COVID-19 vaccination in the greater Melbourne area. We understand that all hospitals are responsible for administrating all vaccine types. Vaccinees must receive a specific vaccine type from one hospital. This model considers all hospitals in a queuing mechanism. Based on a queuing mechanism, we assumed that each vaccinee’s demand has a Poisson distribution. Therefore, we identify the vaccinee’s demand based on the total of the Poisson distributions. We also consider that each hospital owns an  queuing model to serve vaccinees. Our two objective functions (the minimisation of the average of expected wait times and the minimisation of holding and ordering costs of vaccines) contradict each other; vaccinees intend to reduce wait times, while hospitals seek lower costs in ordering and holding vaccines. The final, best solution is a trade-off between the two primary needs of hospitals and vaccinees.

A centralised booking system (CBS) is therefore required to facilitate the administration tasks (i.e., appointments booking and scheduling) for mass vaccination (Abbasi et al., 2020). This CBS would be responsible for checking the eligibility of vaccinees and ensuring that all eligible vaccinees can receive a certain type of vaccine in a specific hospital. The system has a self-service-oriented workflow mechanism with a web-portal interface that enables potential vaccine recipients to enter their personal information and other relevant details (e.g., age, location, existing health conditions, gender, number of vulnerable family members, and daily transport mode). The key functions of this web portal are to allocate vaccinees to the most appropriate hospital based on susceptibility to COVID-19 infection and to identify the most suitable hospital based on its location in relation to the suburbs. According to the information provided through the CBS, each vaccinee is classified into a specific risk priority group. We propose a segmentation mechanism based on recommendations provided by the US government and research published on COVID-19 and influenza vaccine distribution networks (Kazaz et al., 2016, Uscher-Pines et al., 2006, Chen et al., 2020). As COVID-19 vaccines are limited in an emergency, the CBS enables the initial filtering and prioritising of vaccinees by assigning them to predetermined risk clusters and then distributing them to the nearest hospitals for vaccination. Fig. 1 exhibits the structure of the aforementioned vaccine network.

Fig. 1

Proposed conceptual vaccine allocation network.

Assumptions

In this paper, the demand of each hospital is estimated based on the priorities of age groups and their population in each suburb. The supply of vaccines to different hospitals is based on their operational capacity and vaccine availability. We base the allocation of each vaccine type in each time window on the estimated demand and supply. Based on the problem definition, the assumptions of the queuing model are summarised as follows:

Hospitals provide immobile services, which means each vaccinee moves to a specific hospital for vaccination.

The queuing system in each hospital follows an independent exponential service time.

Each vaccinee reaches a hospital at an independent Poisson arrival rate.

The remaining vaccines in each time window, who did not receive the vaccination, can be allocated to the next time window. In other words, the vaccines can be stored in hospitals, based on the hospital’s storage capacity and the vaccine’s characteristics, to be used in the next time window.

In each hospital, the expected service (vaccination) times of all vaccine types in any time window are equal.

Nomenclature

The notations employed in the mathematical model of this study are as follows:

Indices

Parameters

Variables

M/M/C queuing model

An  queuing model is mainly used to analyse service stations in hospitals with more than one server (e.g., a nurse). This model assumes that vaccinees’ arrival rate follows a Poisson distribution and develops an exponentially distributed service time. The number of servers and the services provided by each server are independent of each other. It is also assumed that vaccinees form only one queue for each server and those at the beginning of the queue receive the service (i.e., vaccination) when the servers are released. As soon as the vaccinee is in the queuing system to receive the service, servers will begin vaccination. In this model, the rate of vaccinees leaving the system will be different from the service rate. If the number of vaccinees in a hospital, , is less than the number of servers for the hospital, the number of vaccinees leaving the hospital is , because the time between two consecutive departures is the minimum of  exponential random variables; i.e., the service time for  of the  servers who are working is an exponential random variable with parameter . However, when the number of vaccinees in the system is more than , the vaccinee exit rate is equal to  because the time between two consecutive outputs in this case will be the minimum of  exponential random variables. Consequently, the service rate can be written as follows:

The arrival rate of vaccinees to each hospital is not only related to the expected demand and inventory level of vaccines; it is also related to the efficiency of vaccines and the individual preferences of the vaccinees. The arrival rate of vaccinees is considered an essential parameter that often significantly affects hospitals’ utilisation or productivity rate. Although quantifying individual preference might be difficult for each vaccinee, it is essential to contemplate individual preference in the arrival rate of vaccinees. The utility function is a suitable tool to consider the vaccinee’s preference regarding the decision made under each vaccinee’s uncertainty. The concept of utility function is adopted in some queuing models. In Mohajeryami et al. (2015), the impact of time-based demand response programs on calculating incentive payments has been investigated considering the vaccinee’s utility function. In Ma et al. (2014), vaccinee and utility company viewpoints have been modelled using utility and cost functions to determine electricity consumption. Cicek and Delic (2015) define the total welfare of the society resulting from electricity consumption in terms of all vaccinees’ utility functions and power production costs. Our study adopts the utility function to examine changes in arrival rates related to vaccinees’ perceptions of the vaccine efficacy. Each vaccinee behaves independently based on his/her preference. The different arrival rates of vaccinees can be modelled as utility function, which measures the arrival rate for each vaccinee as a function of the efficiency of the selected vaccine.

This study assumes that the utility to a hospital, , is a function of the vaccine efficacy, ), and a risk aversion coefficient, . The risk aversion coefficient represents different vaccinees’ various risk aversion behaviours (Xie et al., 2021). in this vein, the utility function is a twice-differentiable function with the following properties (Samadi et al., 2012):

Marginal utility is positive which means that as the vaccine efficacy grows, the utility increases which can be written as follows:

Marginal utility is a non-increasing function which means that the marginal utility for the efficiency of vaccines decreases as the efficiency increases, expressed as:

The marginal utility is zero when there is no efficacy expected from a specific vaccine. This means that the arrival rate should be calculated disregarding the vaccine characteristics and type, thus:

The utility function is non-decreasing with the risk aversion coefficient. This property means that the larger is, the greater the utility. This can be expressed as:

Based on the above properties, the utility function is an ascending and concave function that is gradually saturated with the vaccine efficacy. While different functions can satisfy the above properties, this study considers only the exponential utility function to study the arrival rate of vaccinees. The exponential function is a popular utility function used to describe individuals’ risk preferences (Niromandfam et al., 2020). The exponential utility function can be expressed as follows:  Moreover, it is assumed that the arrival rate of vaccinees is related to the inventory level, as well. When the inventory level of vaccines increases, the arrival rate will also increase. Similarly, when the inventory level decreases, the arrival rate to hospitals will be reduced. Therefore, in this study, it is assumed that the arrival rate of vaccinees to the queuing system follows a Poisson distribution (with parameter ) and can be calculated as follows:

According to Khodemani-Yazdi et al. (2019) and considering an  queuing system with congestion in the vaccination procedure of each hospital, the expected wait time of vaccinees for each hospital  during a specified period  is formulated by Eq. (8). This is considered in the first objective function (12). The first term, , calculates the average service time related to the hospital, and the second calculates the average wait in the queue before service. As an essential part of Eq. (8), the idleness probability (having no vaccinee in the queue) in each hospital  is calculated in Eq. (9).

Knowing the probability distribution function (PDF) and cumulative distribution function (CDF) for vaccinees’ wait time in each queue will allow us to calculate the probability that vaccinees will wait in a queue for more than a specified time (Kumar and Sharma, 2018). These functions are formulated in Eqs.  (10), (11). By examining the values of probabilities, the systems (i.e., hospitals) can alter service strategies like changes in the number of servers (i.e., nurses and practitioners) or service rates to manage and analyse the queuing system performance.

Supply chain optimisation model

Using the definitions and assumptions mentioned above, we formulate our conceptual multi-period vaccine allocation model in Eqs. (12), (13), as the objective functions of the optimisation model, subject to Constraints (8), (9), and (14) to (19). Eq. (12) controls the congestion in the system by minimising the expected average wait times of vaccinees in all hospitals and time windows. Eq. (13) minimises the aggregate ordering and holding costs of vaccines for all hospitals during each time window. According to Azzi et al. (2014), ordering costs are usually included in cost pools (e.g., purchase, transportation, contracts with suppliers, and the expediting orders) related to the number of products (here vaccines) in each time window. The holding costs include the inventory service costs and storage space costs. We present more details on the calculation of holding and ordering costs in Section 5.

Eq. (14) is the inventory balance constraint defined for each hospital and vaccine type over each time window. Eq. (15) calculates the total demand of any vaccine type for each hospital in the entire planning time horizon. It considers the weight of priority groups for each hospital to illuminate the diverse range of vaccinees in the suburbs who are allocated to a particular hospital. The total demand is also based on the number of doses defined for each vaccine type. Although the demand in each time window can be estimated by considering an equal distribution (see ), policymakers and hospital managers can also expect any distribution given that . Constraint (16) sets an upper bound for limiting the inventory level to the storage capacity of each hospital based on the vaccine type. Constraint (16) represents the maximum number of vaccines that can be supplied during each time window. The number of servers in each hospital and each time window is also limited to the maximum storage capacity of hospitals using Constraint (18). Finally, Constraint (19) indicates the boundary of decision variables.

The total expected wait time expression (i.e.,  calculated in Eq. (8)) leads to the non-linearity of the model in the first objective function. Thus, the suggested model is based on a bi-objective non-linear programming model.

Solution methodology

Solving a non-linear programming optimisation problem is computationally challenging (Gholizadeh et al., 2021). As the two objective functions of the model are opposing yet interlinked, the “non-inferior” solutions would help improve both objective functions without sacrificing one over the other. This requires adopting a set of Pareto-optimal solutions to help address the multi-objective optimisation problem. To achieve this, we assumed the two opposing yet interlinked objective functions as a multi-objective optimisation problem and developed algorithms that pursue the optimal values of these objective functions. We first introduced the goal attainment algorithm as a classical optimisation. According to Kiresuk et al. (2014), the goal attainment algorithm is one of the efficient methods for solving bi-objective small- to medium-sized optimisation problems. In relation to our NP-Hard problem set and especially with large settings, we expect the goal attainment algorithm to not reach an optimal solution as an exact optimisation method. So, we propose the NSGA-II and the modified NSGA-II (here in this study) as meta-heuristic algorithms to solve the large-scale problems. We chose the NSGA-II as the optimiser algorithm due to its stochastic nature, and high efficiency and performance in complex functions. In addition, several other advantages of this optimiser algorithm are listed as follows (Subashini and Bhuvaneswari, 2012):

It uses the non-dominated sorting techniques to provide the solution as close to the Pareto-optimal solution as possible.

It uses the crowding distance techniques to provide insensibility to the non-convexity of Pareto-optimal front in our problem set.

It also uses the elitist techniques to preserve the best solution for the current population in the next generation.

The following subsections explain the techniques and algorithms developed for this study.

Goal attainment algorithm

The goal attainment algorithm is one of the best techniques to address multi-objective problems proposed by Kiresuk et al. (2014). This technique utilises a group of design goals, ) linked to a group of objectives, ). The formulation of the problem allows the objectives to be under- or over-achievement, making the goals of the initial design imprecise (Kargar et al., 2020). The relative degree of goal under- or over-achievement is managed by a weighting coefficient vector,  = (), and is illustrated as a standard optimisation problem by the use of the following formulation:

Goal attainment minimises a slack variable (), defined as the maximum over  of the expressions in Eq. (20). The term () indicates an element of slackness in the problem. The vector of weighting, , gives the designer the strength of indicating a measure for relative trade-offs between the objectives. For example, setting the vector of weighting  equal to the initial goals illustrates that the same percentage of goal under- or over-achievement, , is obtained. We can also limit the design by setting a specific weighting factor equal to zero (i.e. ). The goal attainment technique makes the intuitive interpretation of the design problem convenient by converting it to standard optimisation method.

Non-dominated sorting genetic algorithm (NSGA-II)

Concerning the complicated nature of our multi-objective non-linear model, an evolutionary algorithm is to be more beneficial for the proposed solution, especially for large-scale problems (e.g., mass vaccination). Asghari et al. (2020) and Karimi-Mamaghan et al. (2020) revealed that the NSGA-II is an excellent evolutionary algorithm and a powerful yet simple method to solve a bi-objective problem. The NSGA-II method holds a broad application to resolve theoretical and practical models (Chaleshtori et al., 2019). Algorithm 1 introduces the steps of the NSGA-II approach in which a population consisting of competing individuals is generated first. The algorithm then ranks and sorts each individual depending on non-domination levels. Following the implementation of evolutionary operations, the algorithm produces a novel offspring pool. It then gathers the offspring and parents prior to portioning the novel pool into fronts, followed by the algorithm that runs the niching technique by determining a crowding distance for each individual. The NSGA-II uses this distance in selecting operators to keep different fronts by ensuring that the individual maintains a certain crowding distance apart from others. We followed steps of the NSGA-II algorithm described in Coello et al. (2007). The steps are presented as follows.

Initialisation

This stage explains how the parameters of NSGA-II (i.e., the maximum number of iterations, population size, crossover rate, and mutation rate) are initialised. No universal constant parameters for NSGA-II exist, as parameters influence the NSGA-II efficiency. In general, parameters are required to be tuned into a particular solution. Thus, NSGA-II parameters should be set as if the highest level of utilisation is obtained. In doing so, NSGA-II must find the best solution in the early steps of its process. To improve the NSGA-II response, the highest probable selectiveness, a constrained initial population size, and a high crossover and mutation rate are utilised. In this paper, we apply the Taguchi approach (explained in Section 4.5) to adjust the parameters. Table A.5 summarises the adjusted parameters of NSGA-II as employed in this study. For the next stage, we use a random generation methodology to produce the initial population.

Table A.5

Tuned values of the parameters in the NSGA-II and the modified algorithm.

Algorithm	MaxIt	nPop	pCrossover	pMutation
NSGA-II	300	50	0.5	0.1
Modified NSGA-II	100	50	0.5	0.3

Encoding and decoding processes

Darmawan et al. (2021) declared that defining the required chromosome is the critical part of an evolutionary algorithm. To code our problem with the proposed evolutionary algorithm, a type of chromosome is introduced; the first part of the chromosome is designed to determine the number of vaccine doses to order, and the second part determines the number of servers (i.e., nurses or practitioners for vaccination) at each hospital within a specific time window. The entire coding process occurs in two phases: coding and decoding. In relation to the coding phase, the chromosome structure is deliberated as a matrix called  with a size of , as shown in Fig. 2. Each gene in this matrix is a real random number between zero and one. The second part of the chromosome is a matrix with the size of . Each element of this matrix is an integer number between one and the maximum number of servers in each hospital.

Fig. 2

Structure of the chromosome defined in the NSGA-II algorithm.

In the second phase, the chromosome is decoded. When the definition of the chromosome and its decoding techniques are carried out, the model’s constraints must be decoded according to Eq. (21). To undertake the decoding phase in the first part of the chromosome, each element in the row of the matrix  is transformed into the minimum and the maximum amount of vaccine  during each time window, , , respectively. In this study, the minimum amount of available vaccines in each time window is considered to be zero (). The second part of the proposed chromosome does not require the decoding process and can determine the number of parallel servers in each hospital.

Cost function evaluation and constraint handling

To meet limitations when dealing with chromosome non-feasibility, the model needs to develop an eminent penalty strategy. Constraints may hold different scales. To equalise the scale, therefore, all limitations are normalised. Thus, the final quantity of the violation of limitations can be measured as the total of the normalised violation of limitations. In the proposed problem, the penalty equals zero once a chromosome is feasible, and it may become a positive value in the resulting solution, even if one of the limitations is not met. For a general formulation of the limitations (), the penalty of the chromosome is specifically defined as:

Where , , and  demonstrate a large positive number, the constraint, and the assigned penalty for the chromosome x, respectively. Concerning infeasible chromosomes, the defined penalty is represented by  in Eq. (23) to achieve the related cost function.  It is noteworthy that the penalty function is applied for all constraints with inequality and not for equal equations. Eq. (14) is an equilibrium constraint (controlling constraint). It determines the levels of vaccine inventory within the planning time horizon for hospitals. Indeed, the inventory levels should be adjusted by the number of vaccine orders as endogenous variables during the optimisation procedure. Moreover, Eq. (15) describes how policymakers and hospital managers can estimate the demand for vaccine doses during the planning time horizon. The estimated parameter of  is therefore applied in Eq. (14) to indicate the amount of inventory and number of vaccine orders.

The selection operator

Up to this point, no novel solutions exist; only good solutions are repeated in the mating pool, and the bad solutions are removed from the parent population. Using the Roulette Wheel Selection method, we adopted and applied several techniques to select the mating pool in this study. As a result, two chromosomes are chosen randomly each time, and more favourable ones would be in the mating pool.

The crossover operator

After generating the members by the selection operator, the crossover operator is applied to them. The two chromosomes derived from the parent population are chosen according to the Roulette Wheel method each time, and the use of this operator creates two offspring chromosomes. In this study, the arithmetic crossover operator is utilised. In the arithmetic crossover, a linear combination of the two chromosomes is defined, on the one hand. On the other hand, the two chromosomes are selected according to the Roulette Wheel method for the crossover operator, i.e., () which produces two offspring chromosomes (). The following formulas define the proposed crossover operator:

Where  is a random number in the range of .

For the second part of the chromosome, a single point crossover is applied. In this crossover, one point is randomly selected from the parent chromosome, and then the offspring is created by relocating the end of the parent chromosome.

The mutation operator

The mutation operator employs offspring created in the previous stage, leading to population diversity. In this study, the Gaussian mutation operator introduced by Deb and Deb (2014) is applied. The steps of the Gaussian mutation operator are as follows:

To create a random number from a standard normal distribution ()

To use Eq. (26) to create offspring  for the parent

Where  is a set in this study.

For the second part of the chromosome, a permutation crossover operator (also known as reverse crossover) is applied. In this crossover, two points are randomly selected, and then the values of the points between the selected points are reversed.

Stopping criteria

The definition of the stop condition is a group of criteria meeting a feasible solution. Different strategies can be used to stop algorithms. In the current study, the suggested NSGA-II stops when the number of iterations reaches the maximum threshold. In other words, we assume the maximum iteration () for our proposed NSGA-II algorithm.

Modified NSGA-II

The classic NSGA-II, described in the previous section, is an efficient multi-objective genetic algorithm for solving multi-objective optimisation problems. The main goal in multi-objective optimisation problems is to find Pareto-optimal solutions. For this purpose, NSGA-II utilises the two basic concepts of the non-dominated sorting and crowding distance. To solve the constrained multi-objective mathematical model, we propose a crucial modification to the solution repairs algorithm. We subsequently prove that the improved algorithm provides better performance for our proposed model.

Solution repairs

In this method, constraints are tackled so that the obtained solutions are within feasible space and so the algorithm can search in a direction that will not produce infeasible solutions. In our optimisation model, Constraints (16), (17) are contemplated in the initial population generation, crossover, and mutation operations stages in the NSGA-II algorithm to repair the infeasible solutions. The solutions are repaired to prevent the movement of solutions that do not satisfy these constraints into the next steps of the algorithm. The three stages of solution repairing, determined by lines with ** in the pseudo-code of Algorithm 1, are conducted during the search process to prevent the entry of infeasible solutions to the next steps.

In the proposed solution repair mechanism, no repair is made if the answer is feasible due to Constraints (16), (17). However, if the answer is infeasible, the repair mechanism is applied to the infeasible solution. This corrective approach prevents the production of infeasible solutions and thus makes the NSGA-II algorithm seek optimum solutions in the feasible space. More specifically, if Constraint (16) is violated, the repair mechanism is applied to each row of the first part of the proposed chromosome (as shown in Fig. 2). If Constraint (17) is violated, the repair mechanism is applied to each column of the first part of the chromosome. In Eqs. (27), (28), we introduce a set of feasible responses and allocate them to the rows and columns of the chromosome, respectively.

The right-hand side terms of Eqs. (27), (28) provide feasible solutions for Constraints (16), (17) , respectively.

Considering that each term of the summation is positive, each term is less than or equal  and satisfies Constraint (16).

Considering that each term of the summation is positive, the term is less than or equal to  and satisfies Constraint (17).

Performance evaluation metrics

Four standard performance metrics are applied to evaluate the performance of the proposed algorithms:

Spacing metric (SP): This measures the standard deviation of the distances between the solutions in the Pareto front. Eqs. (31) to (33) set the metric, where  is the cardinality of the Pareto front,  is a minimal distance of the solution  from the set of remaining solutions,  is the mean of , and  is the cardinality of the objective functions.

Number of Pareto solutions (NPS): This measures the size of the Pareto front and is described as the cardinality of the group of all non-dominated solutions.

Mean Ideal Distance (MID): This measures the convergence rate of the Pareto front toward a specific point (0,0). Eqs. (34), (35) set the metric, where  is the cardinality of Pareto solutions, and  is the cardinality of the objective functions.

Diversity metric (DM): It calculates the relative distance among solutions in Pareto front according to Eq. (36).

Parameter tuning

To obtain valuable solutions, the parameters of both algorithms are calibrated in this Subsection using the Taguchi method (Sadeghi et al., 2011, Sels and Vanhoucke, 2012, Zarandi et al., 2013). Roy (2010) proposed factorial designs to investigate the impact of several factors on an average response. The response is the fitness value of a solution, whereas factors are the parameters of the solution algorithms. Taguchi studied fractional factorial experiments (FFEs) to reduce a large number of experiments in the factorial designs (Roy, 2010). In the Taguchi method, the factors are divided into noise factors () and controllable factors (), where only  can be directly controlled in the experiments. Taguchi proposed a procedure to control  for reducing the variation around the target regarding orthogonal arrays. The design that is less disrupted by the noise factors is robust. There are two ways to perform an analysis of results in the Taguchi method: (i) variance analysis for experiments with a single replicate; and (ii) the use of the signal to evaluate noise ratio () for experiments with multiple replications. The analysis with multiple replications has the better performance, so the  is applied to analyse the solutions. For more information about the Taguchi method, see (Taguchi et al., 2005).

A real case-study

An authentic case study based on the healthcare ecosystem of greater Melbourne, Australia, is undertaken to validate the proposed crisis-induced vaccine supply chain model. We perform several sensitivity analyses for the investigation of objective functions’ sensitivity to the input parameters. Finally, we code multiple numerical examples and calculate the four efficiency metrics to evaluate and compare the performance of the suggested solution methodologies. The proposed model was run in MATLAB2020b software by a computer with Corei7 CPU and 8 GB RAM.

Health networks and their public hospitals across the greater Melbourne are responsible for carrying out COVID-19 mass vaccinations in four different time windows to vaccinate Melbournians until mid-2022 (Moore, 2021). We identified 39 public hospitals (: seven large (over 500 beds), ten medium (250–500 beds), eight small (100–250 beds), and 14 very small (fewer than 100 beds) hospitals) across the greater Melbourne (Department of Health, 2021b). Fig. 3 shows the local hospital networks across Melbourne. We loaded Australia’s shape file for suburbs in Python (using Shapefile and Seaborn modules) and detected Melbourne’s suburbs for our case. According to the State Revenue Office (2021), the greater Melbourne includes 284 suburbs with approximately five million residents. Based on reports published by each local hospital network in 2019 (Department of Health, 2021c), we identified all those suburbs covered by their local hospital networks and estimated the number of people who had to refer to each hospital for vaccination. For instance, the Eastern Health network includes the Box Hill hospital (large), Maroondah hospital (medium), Angliss hospital (small), and Wantirna hospital (very small). As of August 2021, the Australian government has ordered two vaccine types for Victorians (the Pfizer and AstraZeneca vaccines), it has also expressed its willingness to order the Moderna vaccine as a backup.

Fig. 3

Melbourne regional area and its suburbs and hospitals — colours show the population range for each suburb.

We also included the Johnson & Johnson vaccine in the model to compare double-dose vaccines (Pfizer, AstraZeneca, and Moderna) with the single-dose type vaccine Johnson & Johnson (J&J). We present the details of these four vaccine types () in Table 2. The Australian government has undertaken the vaccination effort based on priority age groups to better respond to the crisis. According to the Department of Health (2021a), people aged 60 and over, critical and high-risk workers (e.g. caregivers, nurses, quarantine workers), and people with chronic diseases or underlying medical conditions are at the first layer of priority. The other priority layers belongs to people aged 50–59, 40–49 and 16–39 years respectively. According to the Centres for Disease Control and Prevention (2021), the death rate among people aged 60+ and 50–59 years is 95 times and 35 times higher than other age groups, respectively. Therefore, the inclusion of age priority in our model will help make the model more accurate.

Table 2

Ordering and holding costs per vaccine dose per year in AUD*.

Vaccine	Transportation cost				Total ordering
	Vaccine price per dose	Flight cost (Boeing 777) per dose per 24 h	Average cost from airport to hospital (30 km)		costs per dose
Moderna	44.00	0.32	0.05		44.367
Pfizer	27.00	0.32	0.05		27.367
AsZe	5.50	0.32	0.05		5.867
J&J	13.50	0.32	0.05		13.867
Vaccine	Holding costs				Total holding
	Vaccine Failure per year	Costs of # of freezers per dose per year	Maintenance	Real state costs per dose per year	costs per dose
Moderna	0.53	0.08	0.06	0.24	0.908
Pfizer	0.32	0.36	1.32	0.24	1.920
AsZe	0.07	0.08	0.06	0.24	0.380
J&J	0.16	0.08	0.06	0.24	0.380

To calculate the ordering and holding costs, we adopted the techniques proposed by Li et al. (2021) and Hovav and Tsadikovich (2015). To calculate ordering cost, we considered the price of each vaccine and the average transportation costs from the USA and UK to Melbourne Airport and from the airport to each hospital. The vaccine prices are fixed costs but vary depending on the vaccine type per dose. Transportation costs were calculated based on air shipping costs and road transport costs. We obtained the air shipping costs from The International Air Transport Association (2021) and found that each flight takes approximately 24 h from the USA and 21 h from the UK to the Melbourne Airport. Assuming that vaccines would be transported using Boeing 777-300ER aeroplanes, each flight can transport 1,000,000 doses to Melbourne (Beoing, 2021).

We also identified the leasing cost of each flight, which is approximately AUD13,000.1. We considered this exchange rate per hour. To transport 1,000,000 doses of the four different vaccine types, we developed four equal time windows to provide COVID-19 vaccinations in 12 different scenarios (see Table 3). For example, in Scenarios 1, 5 and 9, the transportation costs are different from Scenarios 2, 6, and 10, as these scenarios present different percentages in vaccine delivery and distribution. In Scenarios 1, 5, and 9, for instance, the government orders 300,000 Pfizer, 600,000 AstraZeneca and 100,000 Johnson & Johnson vaccines. As Pfizer vaccines require special conditions, compared to other types, we calculated the ordering costs per dose to better manage the flexibility of vaccine costs in our scenarios. To calculate road transport costs, we used the average fixed costs from Melbourne Airport to hospitals per hour per truck (4.5-tonne). According to the Transport Industry Council (2021), the fixed cost of each truck (4.5-tonne) per hour per kilometre in 2021 is AUD11.56, and the variable cost is 35 cents per kilometre. We also contacted our network in some hospitals to identify how these vaccines must be transported to hospitals. We assume that hospitals use the Haier  °C freezer (959-litre), which costs AUD25,500 and can store up to 3000 Pfizer vaccines. Each truck (4.5-tonne) can carry up to five Haier  °C freezers (959-litre). For other vaccine types, we understand that only a typical refrigerator temperature is required; therefore, any type of refrigerator can be applied (we assumed the Laboratory Refrigerator BYC-5L666, with a storage capacity of 2500 vaccines, would be appropriate). In Scenarios 1, 5, and 9, for instance, 20 trucks must be employed for 300,000 Pfizer vaccines and 35 trucks for 700,000 AstraZeneca and Johnson & Johnson vaccines. The average distance between Melbourne Airport and the hospitals was calculated to be 30 kilometres. We did not include the fixed and variable costs or insurance expenses during flights, as it would be the vaccine producer’s responsibility to warrant the delivery of vaccines.

Table 3

Objective functions values for the obtained Pareto solutions in the simulated real case scenarios.

Scenario #	Time windows				Vaccine types				Pareto solutions		Average values
	t_1	t_2	t_3	t_4	Pfizer	AsZe	J&J	Moderna	z_1 (minutes)	z_2 (mAUD)	z_1 (minutes)	z_2 (mAUD)
1	20%	30%	40%	10%	30%	60%	10%	0%	23.323600	8.250	23.32365	8.210
									23.323650	8.210
									23.323658	8.160
2					60%	30%	10%	0%	23.155200	9.900	23.1552	9.790
									23.155220	9.680
3					20%	20%	60%	0%	23.831400	9.380	23.8314	9.377
									23.831430	9.370
4					70%	0%	30%	0%	23.071000	14.340	23.0751	14.318
									23.075000	14.300
5	20%	50%	20%	10%	30%	60%	10%	0%	23.040000	7.360	23.0910	7.330
									23.040670	7.340
									23.141800	7.330
									23.141830	7.290
6					60%	30%	10%	0%	23.101200	13.480	24.9700	11.123
									24.970000	10.900
7					20%	20%	60%	0%	23.462000	9.370	23.4600	8.547
									23.831000	8.240
8					70%	0%	30%	0%	23.590000	15.880	23.5900	15.870
									23.597000	15.870
9	20%	10%	50%	20%	30%	60%	10%	0%	23.776000	10.398	23.7700	10.390
									23.776200	10.390
10					60%	30%	10%	0%	24.170000	11.690	24.1700	11.680
									24.175000	11.660
11					20%	20%	60%	0%	23.383000	6.932	23.3830	6.930
									23.383900	6.931
12					70%	0%	30%	0%	23.240000	9.750	23.2400	9.640
									23.244000	9.520

We used the vaccine failure rate, storage, and maintenance equipment costs per dose to calculate holding cost. We adopted the vaccine failure rate per vaccine from the study by Li et al. (2021, p. 3) that calculates this rate as “the average cost of vaccines multiplied by the average power failure rate of storage equipment”. Based on the historical data and our expert network in hospitals, we calculated this failure at 1.4% per year. For example, the failure rate of a dose of a Pfizer vaccine is AUD44.00 × 1.2% = AUD0.53 per year. We also calculated the required number of Haier  °C freezers and their maintenance costs.

We simplified the vaccine allocation among hospitals by giving them equal allocation based on their number of beds. The vaccine allocation in large hospitals was 47% of total vaccines in Victoria per week. In medium hospitals, it was 33%, in small hospitals 15%, and in very small hospitals, 5%. Assuming that all hospitals at the same level will receive the same order, we can generate the number of doses each hospital will receive per week. For example, the Alfred, as a large-size hospital, will receive  doses per week. We understand that each Haier  °C freezer can store up to 3000 doses of the Pfizer vaccine, requiring a total of 100 Haier  °C freezers for all hospitals. Each Haier  °C freezer costs AUD25,500, requiring hospitals to pay a total of AUD2,550,000 to purchase these freezers. We did not calculate these expenses for other vaccines, as we assume our hospitals are well-equipped to store other types of vaccines that need regular refrigeration temperatures. However, freezer and refrigerator variable costs (maintenance and electricity) are estimated at AUD1068 per year, requiring AUD1068 × 100  AUD106,800 for all Haier  °C freezers and up to AUD148,000 per year for refrigerators. We also obtained the required information for the average rents per square metre per year for commercial buildings in different hospital addresses from Real Estate Victoria (2021) a figure that was treated equally for all vaccine types. All ordering and holding costs are reported in Table 2.

Each hospital is considered an  queuing system. , or  in our model, is the number of parallel servers (i.e. nurses or practitioners ready for vaccination) in each hospital in a given time window. We assumed a maximum of ten numbers of parallel servers in each hospital. Fig. 4 demonstrates the number of vaccinations during the last 100 days (by the time of this research) in Victoria. The original time series (the orange chart) is decomposed to its trend and seasonality, aimed at finding the daily arrival rate of vaccinees to the predetermined hospitals. We used the classical decomposition approach with an additive decomposition model. It is noteworthy that trend, seasonal, and remainder components are added to create the time series. According to the average number of vaccinees acquired in each cycle of the seasonal trend and efficiency rates of the vaccines reported in Pilishvili et al. (2021), we estimate the initial inter-arrival rate of vaccinees in every hospital as six people per hour (), which is reasonable in the vaccination process. According to Eq. (7), the inter-arrival rate can change based on the alteration of vaccine inventory during the planning time horizon. To estimate the risk aversion coefficient used in  Eq. (7), we assume an initial value () with a heuristic approach. We also implement other tests in using this coefficient in the model. The associated service time is approximately three people per hour (). We consider four time windows () as the time horizon. Each time window is assumed as three months to help complete mass vaccination in one year (mid-2022). The purpose of the model is to minimise the expected wait time of vaccinees and minimise the aggregate holding and ordering costs of vaccines during the planning time horizon.

Fig. 4

Number of people who received at least one dose of a COVID-19 vaccine in Victoria during the past 106 days (12 Apr 2021–27 Jul 2021) and its decomposed trend and seasonality (COVID LIVE, 2021).

As mentioned previously, we considered several different scenarios for vaccine supply to include conditions in a crisis-induced vaccine supply chain (see Table 3). The expected demand is declared in the percentages of vaccine orders in each time window for each vaccine type (). As the Pfizer and AstraZeneca vaccines are reported more effectively than the Johnson & Johnson vaccine at the time this study was conducted (Centers for Disease Control and Prevention, 2021), in most demand scenarios, we assume that the demand for double-dose vaccines would be higher than the demand for single-dose vaccines, as reported by the Department of Health (2021a). This helps us estimate each hospital’s vaccine storage capacity. We also assume that the maximum supply of each vaccine for each hospital during each window () is 1.5 times greater than the hospital’s expected demand. We then apply a modified NSGA-II algorithm to find near-optimal solutions as we dealt with a large-scale case study at the state level (see Section 5.3 for problem sizes).

The levels of parameters for both the NSGA-II and its modified version are shown in Table A.1. We utilise the Taguchi  orthogonal array to run the experiments and adjust the parameters.

Table A.1

NSGA-II and modified NSGA-II parameters.

Parameters	Value
	Level 1	Level 2	Level 3
Maximum iteration(MaxIt)	100	200	300
Number of population(nPop)	20	30	50
Crossover probability(pCrossover)	0.2	0.5	0.8
Mutation probability(pMutation)	0.1	0.3	0.4

Four different responses, each representing a specific solution quality, are calculated according to the metrics introduced in Section 4.4. When a minimisation problem is involved, MID and SP with lower values better explain the efficiency of the algorithms. In contrast, NPS and DM with greater values better explain the efficiency of the algorithms. Considering that each parameter level of the L9 forms orthogonal arrays, each problem identified in our real case study (39 public hospitals, three types of vaccines, and four time windows as the time horizon) is run five times. The aggregation of the obtained metrics for the five runs are normalised in Table A.2, Table A.3 for the NSGA-II and the modified algorithm, respectively. These normalised values are used in the Taguchi method as the response values of experiments based on different combinations of the parameter levels. Since a solution with the highest performance is desired, the aim is to find the maximum , calculated by:  Where,  is the response in  replication of the Taguchi method and  is the number of replications in the experiments. Table A.2, Table A.3 show the experimental results of NSGA-II and modified NSGA-II under different scenarios of the parameter combinations, respectively. “1”, “2”, and “3” refer to the first, second, and third levels of the parameters. Table A.2, Table A.3 also report . In addition, Fig. A.1, Fig. A.2 exhibit the best values of  (the highest average values). We obtained the optimal parameter values in both algorithms as shown in Table A.5. We implement the proposed Taguchi method by Minitab 15.1 software. After adjusting the parameters, we summarise the results (the Pareto solutions and their averages) in Table 3. The obtained solutions are feasible in both algorithms because none of the models’ constraints has been violated.

Table A.2

Experimental results of NSGA-II.

MaxIt	nPop	pCrossover	pMutation	Sum_1	Sum_2	Sum_3	Sum_4	Sum_5	S/N
1	1	1	1	4	1.2256	1.3333	1.5654	1.1302	3.0422
1	2	2	2	2.0401	1.1119	4	1.3333	1.3333	3.4937
1	3	3	3	4	1.3333	2.7970	1.5272	1.3243	4.5548
2	1	2	3	4	1.3333	1.2834	1.4385	1.7947	3.9236
2	2	3	1	1.500	1.500	1.500	1.4372	4	4.2482
2	3	1	2	3.6667	1.2199	1.5415	0.9491	1.2487	2.3378
3	1	3	2	1.3333	1.6129	1.3333	3.6667	1.1105	3.1971
3	2	1	3	1.7943	1.0701	1.3333	1.3333	3.6667	3.2179
3	3	2	1	1.7575	1.8303	1.2819	1.7005	4	4.8436

Table A.3

Experimental results of modified NSGA-II.

MaxIt	nPop	pCrossover	pMutation	Sum_1	Sum_2	Sum_3	Sum_4	Sum_5	S/N
1	1	1	1	1.1949	1.7863	4	1.6359	2.6314	4.9638
1	2	2	2	2.7776	1.8078	2.3615	1.5782	4	6.6600
1	3	3	3	2.1645	3.9112	1.4420	4	1.2238	5.5280
2	1	2	3	1.3597	1.6786	4	1.8536	1.3764	4.4924
2	2	3	1	1.1062	2.0326	3.6667	2.2777	1.1612	3.8344
2	3	1	2	2.1116	1.7519	4	1.3515	1.3333	4.6278
3	1	3	2	1.3333	1.3333	3.6667	1.9990	1.3333	3.9530
3	2	1	3	4	1.3333	1.0269	1.3709	1.6841	3.0389
3	3	2	1	1.3333	1.7106	2.8741	1.8239	3.6667	5.5274

Fig. A.1

Average  for different levels of parameters in the NSGA-II model.

Fig. A.2

Average  for different levels of parameters in the modified NSGA-II model.

As evident in Table 3, the average wait time of vaccinees in all hospitals and time windows () present a reasonable time with a tight variation range (23.09 to 24.17 min). Our queuing model is expected to control congestion in each time window and allocate sufficient vaccines for meeting hospital demands. In the case of the system’s total cost (), however, the demand violations regarding time windows and vaccine types will cause significantly different values (6.931 to 15.883 million Australian Dollars). This unveils the importance of the second objective function in our proposed model. Given that the Pareto solutions from Scenario 11 lead to lower average aggregated ordering and holding costs, we choose this option as the optimal scenario. Scenario 11 presents the best vaccination plan for scheduling the highest vaccine demand (50%) in the third time window using the single-dose vaccine, Johnson & Johnson, as the primary type (60%). Furthermore, comparing the results of scenarios that use a higher proportion of double-dose vaccines (i.e., Pfizer and AstraZeneca), in each bundle of time window allocations (e.g., Scenarios 1, 2, and 4 in the 20-30-40%-10% bundle or Scenarios 5, 6, and 8 in the 20-50-20%-10% bundle) intuitively demonstrates the privilege of the cheaper double-dose vaccine type (i.e., AstraZeneca). According to Usher (2020), governments prefer to purchase COVID-19 vaccines from their international allies due to the high accessibility (i.e., Pfizer or AstraZeneca imported from the USA and UK, respectively (see Table 2)). Therefore, the trade-off between these scenarios will benefit policymakers in optimal vaccine ordering during each time window. For instance, comparing Scenarios 2, 4, 6, 8, and 10, we found Scenario 2 as the best allocation strategy with lower total cost. Each combination of vaccine types (e.g., 60% of Pfizer, 30% of AstraZeneca, 10% of Johnson & Johnson, and 0% of Moderna) should be allocated 20-30-40%–10% over the entire planning time horizon. As evident in Table 3, the Pareto set obtained for each scenario reveals a correlation between two objective functions, meaning that an increase in the first objective function may cause a decrease in the second objective function. This trend confirms the opposing yet interlinked behaviour of two objective functions.

Best scenario analysis

We consider Scenario 11 (the best option introduced in Table 3) for our following analyses. The model outputs suggest the optimal order and inventory amounts for each vaccine type in each hospital over a specific time window. Fig. B.3 graphically shows the accumulative number of orders for each suburb of Melbourne during each time window, considering the first Pareto solution. Fig. 5 illustrates the order amount and inventory levels for a large hospital (Box Hill) compared to a small hospital (Williamstown). It is recommended to order almost the same amount of vaccines in all time windows but to use the inventories remaining from the previous ones. In this plan, the highest demand (50%) is estimated in the third time window, and the inventory level of the fourth time window will be used in accordance with its planned demand (i.e., 20-10-50%–20%). The figures also illustrate that although the demand for one vaccine type can be higher than others, the model may plan to order lower amounts and use the remaining inventory to allocate vaccines optimally. This is an example of inventory management for the vaccines using their particular specifications.

Fig. B.3

Optimal ordering levels for the Melbourne’s suburbs in the best scenario.

Fig. 5

Optimal inventory and ordering levels for large (Box Hill) and small (Williamstown) hospitals for the first Pareto front in the best scenario.

The other significant finding of the model is associated with the number of parallel servers required in each hospital over each time window. Fig. 6 shows the number of parallel servers and their utilisation rate for both hospitals mentioned above (i.e., Box Hill and Williamstown). The utilisation rate refers to the probability that servers (i.e., caregivers, nurses, and practitioners) would not be vaccinated prior to their service. This probability is shown as . Fig. 6 shows that the maximum number of servers (9 servers in our case) is only required in the third time window in the Box Hill hospital (a large hospital). The corresponding findings in the figure also demonstrate although we use the maximum number of servers in this time window, the utilisation rate of the hospital may not be 100%. This can be justified by the formula of the utilisation rate related to the expected arrival rate of vaccinees to a specific hospital () and the expected service rate of the hospital (). The highest utilisation rates will be concluded by near eight parallel servers; see the third time window for the Williamstown hospital (a small hospital), shown in Fig. 6. It is noteworthy that our objective is not to maximise the utilisation rates but to minimise the wait times of vaccinees.

Fig. 6

Number of servers and utilisation rate for large (Box Hill) and small (Williamstown) hospitals considering the best scenario for demands.

Sensitivity analysis

As mentioned before, minimising the expected wait time for the vaccination process is considered the first desirable goal for vaccinees. However, the minimisation of vaccine holding and ordering costs is desirable for hospitals and is considered the second research objective. We understand that the duration of the vaccination process depends on the availability of vaccines, vaccine types, and the number of servers (i.e., nurses or practitioners). Therefore, it can be inferred that by allocating the optimal number of servers for mass vaccination, more vaccines (and vaccine types) would be available. This would result in more vaccinees receiving vaccines during a specified time window. Thus, vaccinees can expect that the wait time for their vaccination would be significantly reduced. Simultaneously, the proposed model provides the optimal number of servers (i.e. nurses and practitioners) to reduce hospital vaccine holding and ordering costs. In order to numerically demonstrate the contradiction of objective functions, we solve the problem for a convex combination of both objective functions by minimising  for various  values selected from the interval  (Eisenhart, 2003, Ida, 1997, Rangan and Poolla, 1997). The obtained results, shown in Fig. 7, confirm that the objective functions are conflicting yet interlinked, and the selection of separate objective functions is necessary for our problem setting (Chaleshtori et al., 2020).

Fig. 7

Demonstration of the contradiction within the objective functions using .

We recognise the effect of our proposed model’s most significant parameter in providing managerial insights. In the previous section, we observed the changes of the second objective function, while the first one was stable. As evident in Eq. (8), the first objective (the expected wait time of vaccinees) is closely related to the expected service rate of vaccinees to hospitals () (see Fig. 8). We conduct all experiments for the best scenario detected in our real case problem (see Table 3). All findings in this section confirm that this parameter significantly influences the costs associated with the system. However, it can substantially change the vaccinees’ expected wait times. Both charts shown in Fig. 8 validate the inverse relation between the expected service rate of vaccinees and their wait times. Comparing the findings in Fig. 8, Fig. 8 demonstrates that the increased service rate is most readily in large hospitals (e.g., Box Hill hospital). Therefore, it is inferred that the  parameter must be determined more carefully for large hospitals than small or very small hospitals.

Fig. 8

Expected service rate of hospitals for large (Box Hill) and small (Williamstown) hospitals in the best scenario.

We also observe a non-monotonic trend in both objective functions. As discussed previously, the first objective function is non-linear and the contradiction between two objective functions leads to this non-monotonic trend in both objective functions.

Performance evaluation

To verify the efficiency of the proposed solution approaches, we implement the model in multiple numerical experiments. The input data for these ten problem sizes are developed according to the real case data and reported in Table A.4. The corresponding outcomes, listed in Table 4, show that our proposed modified NSGA-II algorithm outperforms the classic NSGA-II and goal attainment algorithms. In all instances, the DM, SP and NPS metrics of the modified NSGA-II are better than the corresponding values for the two other algorithms. In the MID measure, the modified NSGA-II algorithm also outperforms the other algorithms. Moreover, we can conclude that both the classic and modified NSGA-II algorithms are efficient tools for solving small- to large-sized problems. The goal attainment approach could not solve large size (#7–10) problems in less than one hour. This is expected, as our problem is NP-hard, and exact methods cannot reach an optimal solution for large-scale problems in a reasonable computational time.

Table A.4

Inputs parameters for numerical examples.

Problem	Size of the problem	Number of	Assigned	Demand	Ordering	Holding	Maximum	Maximum
number	(K-T-F)	servers (e.g. nurses)	server		cost	cost	capacity	available vaccine
1	(5-2-3)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
2	(5-3-5)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
3	(5-4-3)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
4	(6-3-2)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
5	(6-3-10)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
6	(4-9-12)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
7	(5-12-14)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
8	(6-18-15)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
9	(6-24-20)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)
10	(6-24-40)	U(1,10)	U(0,1)	U(4E+02,8E+04)	U(5,30)	U(1,5)	U(12E+04,3E+05)	U(4E+04,95E+04)

Table 4

Performance metrics for numerical examples using the proposed algorithms.

Problem	NPS			DM			SP			MID
number	mod. NSGA-II	NSGA-II	Goal-att.	mod. NSGA-II	NSGA-II	Goal-att.	mod. NSGA-II	NSGA-II	Goal-att.	mod. NSGA-II	NSGA-II	Goal-att.
1	4	5	2	1.03E＋04	8.12E＋03	7.33E＋03	3.83E＋03	4.57E＋03	4.49E＋03	40.22	74.72	76.58
2	4	4	8	3.78E＋04	3.29E＋04	3.04E＋04	3.097E＋03	3.28E＋03	3.44E＋03	21.37	41.61	42.5
3	6	6	4	1.78E＋04	1.68E＋04	1.68E＋04	4.21E＋03	5.05E＋03	5.08E＋03	22.81	35.67	36.55
4	8	7	7	3.91E＋04	3.39E＋04	3.06E＋04	3.57E＋03	4.74E＋03	4.63E＋03	23.89	36.48	36.64
5	4	4	7	1.22E＋04	1.2E＋04	1.06E＋04	946.40	1.79E＋03	1.78E＋03	39.72	51.18	52.47
6	8	2	6	9.71E＋03	8.9E＋03	7.65E＋03	1.27E＋03	2.15E＋03	2.22E＋03	45.85	87.5	87.26
7	3	2	–	4.08E＋04	4.02E＋04	–	3.64E＋03	6.17E＋03	–	35.54	53.19	–
8	8	6	–	1E+04	8.24E＋04	–	4.04E＋03	4.79E＋03	–	28.53	50.59	–
9	3	3	–	2.27E＋04	1.9E＋04	–	2.79E＋03	5.54E＋03	–	41.38	56.38	–
10	2	2	–	2.39E＋04	2.05E＋04	–	1.15E＋03	2.003E＋03	–	30.65	38.84	–

Considering that the results may be inconsistent after changing numerical examples, we use several statistical hypothesis tests to ensure the reliability of previous comparisons. Thus, the implementations of the solution algorithms concerning several statistics are compared by employing the non-parametric Friedman test (Fadeyi, 2021). The results obtained from the SPSS software ( test, degree of freedom (), and significant level ()) are summarised in Table 6. The results illustrate significantly varied levels of performance between the algorithms in the (), (), and () metrics.

Table 6

Statistical outputs obtained from the non-parametric Friedman test, applied on the results of the numerical examples.

Performance	Test statistics
metric	χ^2	df	Sig.
NPS	1.625	2	0.444
DM	12	2	0.002
SP	9	2	0.011
MID	10.33	2	0.006

Besides the general validation of the performance metrics, gained from the Friedman test, we apply a Post-Hoc test called Wilcoxon signed-rank to distinguish which algorithms are significantly Superior to the others (Lin et al., 2021). Table 7 summarises the results.

Table 7

Statistical outputs obtained from the Wilcoxon signed-rank test, applied on each pair of results in the numerical examples.

Wilcoxon signed-rank	Performance	Test pairs
Statistics	metric	Mod. NSGA-II vs. Cls. NSGA-II	Mod. NSGA-II vs. Goal att.	Cls. NSGA-II vs. Goal att.
Z	NSP	−1.841	−0.106	−1.089
Sig.
Z	SP	−2.666	−2.201	−1.05
Sig.
Z	DM	−1.784	−2.201	−2.201
Sig.
Z	MID	−2.803	−2.207	−1.782
Sig.

The results of Table 7 indicate that the modified NSGA-II performs better than the other two algorithms regarding (), (), and () metrics. However, the red values of the table show that the performance of the specified metric cannot be significantly better than that within the pair of algorithms. This weakens the results gained from Table 4, Table 5 regarding the privilege of the classic NSGA-II over the goal attainment approach. Despite this uncertainty in the classic algorithm, all tests approve that the proposed modified NSGA-II algorithm surpasses the other algorithms and can be used without failure.

Table 5

Computational times (seconds) for solving the proposed algorithms.

Problem #	mod. NSGA-II	NSGA-II	Goal-att.
1	57	58	94
2	102	100	220
3	136	135	248
4	400	410	680
5	522	529	935
6	736	745	1672
7	898	987	≥3600
8	991	1002	≥3600
9	1102	1143	≥3600
10	1196	1225	≥3600

Analysis of the queuing system

Based on the exponential utility function introduced in Eq. (6), the arrival rate at hospitals is in contrast to vaccinees’ risk aversion coefficient that decreases (or increases) by increasing (or decreasing) the efficiency of vaccines. In this study, we consider the vaccinee risk aversion behaviour to be an uncertain parameter because social studies are needed to identify this coefficient. Thus, the probabilistic nature of the vaccinee’s behaviour is modelled via its PDF as a uniform distribution, with a range of , where  is the upper bound of the distribution. We also understand that the vaccine efficacy varies depending on the vaccine type (95%, 95%, 70%, and 72%, for Moderna, Pfizer, AsZe, and J&J vaccines respectively) (Lopez Bernal et al., 2021). As per our assumptions, the arrival rate to a hospital follows a Poisson distribution. Therefore, as per Eq. (6),  follows a Poisson probability distribution with a parameter equal to the sum of the vaccines’ utility function. To satisfy the stability condition of the  queuing system, we adjust the service rate of each hospital, , so that the utilisation rate of the hospital, , is always less than one. It is presented as follows:

To plot the PDF and CDF of wait times for the queues of a hospital, besides the above-mentioned setup, we consider seven hours for a working day and the average number of servers in the hospital (five servers). Fig. 9(a) illustrates the wait time PDF. Increasing the risk aversion coefficient of vaccinees to a specific value (the peak for each PDF trend) increases the expected wait time of vaccinees. For values greater than the peak, increasing the risk aversion coefficient of vaccinees leads to a decrease in the expected probability of the wait time in the queue for receiving the vaccine service. Using this plot, The crisis-induced vaccine supply chain designers can find the peak values as the best wait times (regarding the assumed risk aversion) for setting up the queuing system. Moreover, it can be inferred from the positive skew of the PDFs that the average queue wait time is greater than the median wait time. This means that the wait time for half of the vaccinees is shorter than the average wait time for receiving the vaccines.

Fig. 9

Probability distribution function (PDF) and cumulative distribution function (CDF) of wait times in a queue of vaccines concerning various risk aversion coefficients.

Fig. 9(b) shows the CDFs of wait time of vaccinees in a queue of vaccines concerning various boundaries () for the risk aversion coefficient. The crisis-induced vaccine supply chain designers can employ these CDFs to calculate the probability of an expected wait time. In other words, the designers can determine how probable a wait time is for a hospital’s queues considering a specific risk aversion for vaccinees. This will give them insight in determining proper values for the parameters of the queuing system aimed at controlling any potential congestion in the queues of a hospital.

Discussion and conclusion

The healthcare systems have been experiencing significant burdens since the emergence of COVID-19. The number of hospitalised infected people has significantly increased, resulting in a shortage in healthcare staff to administrate COVID-19 vaccines in a timely manner. We argue that responses to such emerging situations require the development of crisis-induced supply chains that increase the pace of supply and distribution and reduce ordering and storage costs. Simultaneously, governments in general and public health authorities in particular must take proactive actions in such supply chains to ensure that vaccines are equally available to citizens. Vaccines must be administered promptly and with respect to the rate of vulnerability to COVID-19. The timely distribution and administration of COVID-19 vaccines have two important outcomes: (i) reducing the ordering and holding costs of such vaccines that require special conditions for storage and (ii) increasing the public resilience against the COVID-19 pandemic and boosting vaccinee satisfaction, which has proven to have positive social and economic consequences (Pescaroli et al., 2021). To achieve this, healthcare systems must respond effectively and efficiently to vaccine demand.

Therefore, in this paper, a bi-objective non-linear programming optimisation model was proposed to design optimal distribution planning for ordering and holding different types of COVID-19 vaccines. Several planning time windows were determined to cover the growing demand of COVID-19 vaccines with limited supply, especially in situations where certain age groups have priority of vaccination over others. Therefore, we addressed two opposing yet interlinked objectives in this study. The first objective of the model focuses on the minimisation of the average wait time of vaccinees in hospitals, while the second objective function minimises the total ordering and holding costs of vaccines.

The first objective provides policymakers with a set of directions to better manage congestion in the predetermined hospitals using a queuing system. We also defined a utility function to study potential changes in arrival rates. Indeed, different vaccinees may have different perspectives toward vaccine types and their efficacy. The utility function was also related to a risk aversion coefficient. We assumed that the arrival rate to a hospital is related to a Poisson probability distribution with a parameter equal to the sum of the vaccines’ utility function. We also proposed two optimisation algorithms (i.e., goal attainment and NSGA-II) to solve the opposing yet interlinked nature of the objective functions of the proposed model (i.e., contradictions between arrival rate of vaccinees to a hospital versus the ordering and holding costs of a vaccine type).

We proposed some modifications in the NSGA-II algorithm as solution repairs to introduce a modified algorithm. The performance of the proposed algorithms was investigated by four metrics. We then compared the results to find a better solution approach. The results demonstrated that the heuristic approaches (the classic and modified NSGA-II algorithms) perform better than the goal-attainment technique in solving large problems. Moreover, the modified heuristic algorithm outperformed the classic method, meaning that our new algorithm can also solve NP-hard optimisation problems with reasonable computational time.

Based on the results obtained by the non-parametric Friedman and Wilcoxon signed-rank tests, the performance reported by the modified NSGA-II was significantly higher than that in other algorithms involving all performance metrics. Consequently, we proposed the modified NSGA-II as a specific algorithm for finding optimal solutions in a crisis-induced vaccine supply chain. In addition, we discussed how to expect a proper wait time for a hospital queue based on the corresponding PDF and CDF diagrams and various risk aversion levels.

The proposed model was implemented using an authentic case study in the greater Melbourne to show how to manage congestion and costs in the ordering and storing of vaccines in the proposed hospitals. For a better understanding of the input parameters of the queuing model and how to estimate them, several recommendations were provided. The findings offer a way to coordinate a crisis-induced vaccine supply chain through optimising the number of required hospitals, the minimum expected wait times for the vaccination of each vaccinee, and the maximum number of orders and inventories in each hospital during a specific time window. We also conducted a sensitivity analysis on the main parameter to evaluate the model’s performance and facilitate the provision of managerial insights. The findings validated the model and determined the importance of the expected service rate of medical facilities — especially for large hospitals.

Implications for policymakers and healthcare managers

This study proposed and validated a practical model for vaccine distribution in a crisis-induced supply chain. As of August 2021, many countries (e.g. Australia) were unable to produce COVID-19 vaccines locally due to the unavailability of production technology, the lack of R&D investments, and the need for immediate action regarding the COVID-19 pandemic. Therefore, healthcare systems in these countries are reinforced by the public to quickly import the best available vaccines. This demand has also caused international competition in purchasing and distributing internationally recognised vaccine types among developing economies. The shortage in the COVID-19 vaccine supply and the high cost of transportation, storage, and administration in crisis-induced vaccine supply chains cannot be easily resolved unless all the issues mentioned above are addressed. Thus, the findings of this study enable policymakers and hospital managers in several ways.

First, the proposed method and model can help policymakers quickly and responsively develop optimal operational decisions for vaccine supply chain distribution in emergencies under reservations arising from limited supply, diversity of vaccine types, unfamiliarity with unexpected consequences, and disease vulnerability. As evident in prior studies (Tavana et al., 2021, Ferranna et al., 2021), the prioritisation of certain groups over others (e.g., based on age and critical and high-risk jobs) for vaccination will significantly reduce the number of infected and dead cases but does not necessarily guarantee herd immunisation at the national level. We recommend that policymakers consider two opposing yet interlinked constraints in mass vaccination to succeed in herd immunisation in a short time: (i) wait time of vaccination and (ii) total ordering and holding costs of vaccines. At the national level, vaccination requires policymakers to maximise the utilisation of healthcare resources, which has proven to increase pressure on healthcare ecosystems both financially and psychologically (Chen et al., 2021, Foy et al., 2021). The findings suggest the best scenario for the optimal number of orders and amount of inventory for each vaccine type at each hospital, indicating that our model can help hospital managers make optimal operational decisions for vaccine ordering, storage, and distribution. However, we recommend that hospital managers not only consider the demographic information of their covered regions (e.g., the population of each suburb, the diversity of age groups) but also assess their vaccine distribution and administration capacity to better coordinate the vaccine supply and reduce ordering and storage costs of the vaccines. Vaccine supply chain coordination in emergencies requires a partnership between governments to better manage vaccine orders and inventories. Our findings help local and federal governments coordinate vaccine distribution locally to reduce unnecessary vaccine ordering and storage, allowing other local governments to meet ordering requirements. More particularly, our model facilitates the process of vaccine ordering and storage for state governments by enabling them to take great responsibility in vaccine ordering and distribution on their own, minimising the centralised role of federal governments in the process of vaccine distribution. Indeed, although we understand the benefits of centralised structures such as (i) having the capacity to mount a quick response to a crisis, (ii) preventing local competition over resources, and (iii) increasing equity in vaccination (Desson et al., 2020), it may – to some extent – increase ordering and holding costs if federal governments show inflexibility in the local decision-making processes. In other words, enabling state governments allows a flexible response for public hospitals to manage vaccine inventory in accordance with their resources and local orders that can in turn reduce resource mobilisation across states (e.g. vaccines and nurses transportation) to meet the required demand. Therefore, the findings of this study can encourage state governments to develop partnership with local public hospitals and perform responsively toward the growing demand for vaccination in their service regions. It in turn reduces the burden on the federal government and enables state governments to better recognise vaccine eligibility and prioritisation. In this view, local public hospitals can better connect vaccine demand and supply by leveraging local call centres, online registration and vaccine administration that can occur in a short period of time. In addition, as healthcare networks and main public hospitals can more quickly make decisions about their vaccine inventory and storage capacity, the likelihood of supply disruption in the crisis induced vaccine supply chain will be minimal.

Our model can also be applied in middle-to-low economies, where healthcare ecosystems may not be equipped appropriately to respond to the high demand for vaccination, ordering and holding costs of vaccines on a large scale due to financial restraints. Our queuing model adopted a four-stage protocol of mass vaccination. The vaccination protocols may vary depending on vaccine supply and demand. These protocols may prioritise vaccinees based on their level of vulnerability to the COVID-19 pandemic. We applied these four stages based on the COVID-19 vaccine administration protocol proposed and implemented by the Australian government, although our model can be flexibly applied to other vaccine administration processes (e.g. 3- or 5-stage protocols). The policymakers can also select the best vaccination scenario (e.g., the best proportions of vaccine types) that will minimise the vaccination wait time and maximise the capacity of vaccine ordering and holding in their healthcare ecosystems.

As mentioned previously, the logistics of a crisis-induced vaccine supply chain is with a high-level uncertainty due to the nature of the emergency. It requires immediate actions toward international agreements, transportation, storage and administration. In this vein, the frequent change in the definition of ’full vaccination’ (shifting from 2- to 3- or 4-dose vaccination) has put much pressure on the federal governments to maintain a centralised high-level standard of vaccination. Our model can help routinise daily vaccine demands and supply in local areas to develop the vaccination national plan in synchronicity, providing a relatively practical toolkit at both strategic and operational decision-making levels to enhance vaccine distribution processes in emergencies.

Future research directions

Although this proposed model performs well in vaccine distribution planning, especially in emergencies, the model has limitations that future research can investigate. First, in the proposed model, we only considered a capacity for the servers, and there was no restriction on the daily entrance of vaccinees to hospitals. In reality, however, a capacity limit (especially for the booking system) can be demonstrated to make the model more accurate. Other types of queuing models, like  can be implemented to model this problem ( represents the size of the calling population in the CBS). In addition, considering the behaviour of vaccinees in queues like bulking, reneging, or jockeying, recalculating the average wait time would be an interesting topic for future research.

Moreover, the inclusion of some uncertainties may extend the model and improve the accuracy of findings. For instance, vaccine producers may fail to fulfil orders (supply uncertainty) due to the international competition in emergencies or vacinees may refuse to accept the second/third/fourth dose of a special vaccine (demand uncertainty) due to vaccine side effects. In addition, the efficacy rate of vaccine types is also a great concern in many developed economies, as it may fail the national vaccination plan for these countries and cause massive pressure on economy. Future research is encouraged to develop a specific ratio for each vaccine type to include the efficacy rate in the model.

Policymakers can allocate financial incentives to allow citizens to select a vaccine type from a diverse range of vaccines based on their efficacy rate, price, and number of doses. Our model can be applied in situations where these incentives are essential. However, we direct future research to consider policy objectives in developing the model to make the analyses more precise, helping policymakers better decide which vaccine type would be cost-beneficial for citizens of different socio-economic statuses. For example, our analyses did not indicate which vaccine administration protocol (i.e., 3-, 4-, or 5-stage protocol) minimises the holding costs or wait times. This can direct future studies to consider the role of government’s initiatives and interventions in vaccine distribution in emergencies. Future research can examine the implementation of the model with other optimisation algorithms and compare the performance of the obtained solutions with the proposed algorithms.

We used an inventory model for our supply chain design considering the main costs of the chain (i.e., holding and ordering costs). Another holistic view can integrate the spoilage and shortage costs of vaccines in the model. Quantifying these costs and minimising the resultant total cost for any selected crisis-induced supply scenario will enhance the performance of the vaccine supply chain.

We did not analyse the public acceptance of vaccines in our model. Indeed, these analyses require employing other research designs in which surveys on the vaccinee’s experience would be the central component of data collection and analysis. Our review indicated that the public acceptance of vaccines plays a crucial role in succeeding the vaccine administration and minimising the ordering and holding costs. Therefore, we encourage future research to develop and examine the vaccinee’s experience using a survey-based research design to help understand how vaccinees adopt new vaccines especially in emergencies.

CRediT authorship contribution statement

Hamed Jahani: Concept, Design, Analysis, Writing – review & editing. Amir Eshaghi Chaleshtori: Concept, Design, Analysis, Writing – review & editing. Seyed Mohammad Sadegh Khaksar: Concept, Design, Analysis, Writing – review & editing. Abdollah Aghaie: Concept, Design, Analysis, Writing – review & editing. Jiuh-Biing Sheu: Concept, Design, Analysis, Writing – review & editing.

Declaration of Competing Interest

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

k	Vaccine type, k=1,2,…,K.
t	Time window, t=1,2,…,T.
f	hospital, f=1,2,…,F.
j	Priority group, j=1,2,…,J.

s_k	The number of doses defined for each COVID-19 vaccine k, s_k=1,2.
δ_f	Percentage of severs (i.e. nurses, caregivers and practitioners) assigned to undertake vaccination in hospital f.
w_jf	The weight of priority group j in hospital f.
n_kf	The number of estimated vaccinees interested in receiving vaccine k in the CBS developed for hospital f.
D_ktf	The amount of vaccine k that a vaccinee demands in hospital f over time window t.
a_ktf	Ordering cost of vaccine k for hospital f at time window t.
b_ktf	Holding cost of vaccine k for hospital f at time window t.
μ_f	The service rate of hospital f; =1/E[Service time].
Ikfmax	Maximum capacity of hospital f for holding costs of vaccine k.
Xktmax	The maximum amount of vaccine k that can be supplied at time window t.

c_tf	The number of parallel servers in hospital f over time window t.
W_tf	The expected wait times of vaccinees in hospital f at time window t.
X_ktf	The number of vaccine k orders in hospital f over time window t.
I_ktf	The amount of inventory of vaccine k in hospital f at the end of time window t.
λ_f	The arrival rate of vaccinees to hospital f; =1/E[Inter arrival time], where E[.] is the expectation function.