Bi-objective optimization for a multi-period COVID-19 vaccination planning problem.

This work investigates a new multi-period vaccination planning problem that simultaneously optimizes the total travel distance of vaccination recipients (service level) and the operational cost. An optimal plan determines, for each period, which vaccination sites to open, how many vaccination stations to launch at each site, how to assign recipients from different locations to opened sites, and the replenishment quantity of each site. We formulate this new problem as a bi-objective mixed-integer linear program (MILP). We first propose a weighted-sum and an ϵ -constraint methods, which rely on solving many single-objective MILPs and thus lose efficiency for practical-sized instances. To this end, we further develop a tailored genetic algorithm where an improved assignment strategy and a new dynamic programming method are designed to obtain good feasible solutions. Results from a case study indicate that our methods reduce the operational cost and the total travel distance by up to 9.3% and 36.6%, respectively. Managerial implications suggest enlarging the service capacity of vaccination sites can help improve the performance of the vaccination program. The enhanced performance of our heuristic is due to the newly proposed assignment strategy and dynamic programming method. Our findings demonstrate that vaccination programs during pandemics can significantly benefit from formal methods, drastically improving service levels and decreasing operational costs.

Introduction

The Coronavirus Disease 2019 (COVID-19) was first reported in late 2019. Since then, it has been declared a pandemic and spread worldwide, heavily threatening human lives in multiple waves. COVID-19 is highly contagious compared to previous pandemics [1], and it has dramatically hindered production in entire industries and exposed the fragility of global supply chains (SCs) [2]. It has a substantial impact on the global economy and has caused economic recession and widespread job losses [1]. The COVID-19 pandemic has led to disruptions and substantial losses in global SCs [3], where such disruptive impacts frequently yield ripple effects [4]. While most businesses were aware of the negative and catastrophic consequences of the ongoing pandemic, they lacked direction on how to simulate SC disruptions and their effects on performance in pandemic scenarios. Due to a lack of such guidance, delayed reactions and absence of awareness of the pandemic’s repercussions resulted in delays, high coordination efforts, and protracted shortage periods resulting from the late deployment of recovery activities [5]. Under such a disrupted scenario, different strategies and actions are required to address this challenge, including robust SC resilience strategies [6], [7], [8] and rendering global SCs more integrated and digitally ready to help improve the quality of response to epidemic-related interference [9]. In addition, optimization tools can help alleviate such problems [10]. In this context, we propose a massive vaccination decision support tool based on formal methods to optimize conflict objectives.

To break the transmission chain and mitigate these impacts, the World Health Organization (WHO) has issued specific health and social distancing guidelines [11], which have been implemented in most countries. To minimize physical contacts, governments closed many social places and even stopped public transit services at the early stage of the pandemic [12]. As infection rates peaked, authorities in numerous countries enforced a lockdown, such as in Italy, the US, and Australia [13], [14], [15]. Although reasonable control over the pandemic has been achieved through ongoing restrictions, infections still surge rapidly in new waves [16]. Moreover, COVID-19 is still circulating with mutations of the novel coronavirus (e.g., Delta, Kappa, and Omicron variants), generating a constant pandemic threat for upcoming seasons in various countries [17], [18].

Immunization with a safe and effective vaccine can reduce COVID-19-related illnesses, hospitalizations, and deaths, besides restoring societal functioning [19], [20]. To achieve mass vaccination, many governments attempt to fast-track the development of COVID-19 vaccines [21]. Various agencies are simultaneously working to develop and manufacture a vaccine to effectively cease or at least drastically decelerate the spread of the pandemic [22], resulting in various vaccine candidates being authorized [23]. More than 7.69 billion doses of vaccines will have been distributed worldwide by Nov 20, 2021 [24]. It is widely accepted that the risk of intermittent outbreaks of COVID-19 will continue until safe and effective vaccines become globally available and a vaccination program is successfully implemented [25]. Vaccinating a large population helps achieve herd immunity, which is an effective method to fight the virus [23], [26]. It is estimated that 70–85% of the population needs to be vaccinated to reach herd immunity [27].

However, the endgame of the COVID-19 pandemic is not vaccines, but vaccination [28]. As a result, the largest vaccination program in human history is currently in action, pushing the COVID-19 vaccine supply chain to its limits [23]. When launching these massive vaccination programs, many countries have opened online appointment channels, where residents must specify their personal information, including their addresses and expected vaccination dates, when making a reservation for vaccination. Given this reservation information of a large population, decision-makers must establish a vaccination plan subject to various constraints, e.g., limited budget and vaccine availability. In practice, many mass vaccination plans are made manually, leading to high operational costs and low service levels. An optimal vaccination plan involves the operation of multiple vaccination sites, the distribution of vaccines, and the allocation of recipients to vaccination sites. These decisions are tightly correlated, which indicates that a centralized and integrated optimization can save costs and improve service levels. Therefore, it is of great value to investigate how to establish efficient massive vaccination programs to fight COVID-19.

Although the importance of massive vaccination programs is well recognized, various challenges greatly restrict the efficiency and effectiveness of vaccination plans against COVID-19. Generally, the resources available to tackle a vaccination plan are limited, while opening a vaccination site uses expertise and resources, e.g., a facility, doctors, nurses, and other staff. From an operational cost perspective, a decision-maker would want to open only a few vaccination sites each day to save costs. Decision-makers thus face the challenge of allocating and scheduling limited resources with a high expected service level. Consequently, the impact on the regular operation of the healthcare system can also be mitigated. However, fewer opened vaccination sites may force the population to travel a long distance to get vaccinated, which leads to a low service level. As indicated by Li et al. [29], the convenience of receiving vaccination services would affect recipients decisions, especially travel distance and time. Moreover, an equitable optimization of cost and service level must also be taken into account [30]. Therefore, decision-makers face a multi-objective optimization problem when making a vaccination plan, in which the following decisions that are relevant and even conflicting must be simultaneously made to balance the operational cost and the service level. First, as vaccines are stored at a central warehouse, decision-makers must determine when to replenish each vaccination site. Second, they must determine the work schedule of each vaccination site, i.e., when to open each site. Third, for each opened site, to avoid work overload or low capacity utilization, the number of vaccination stations to be launched must be determined according to the number of recipients allocated, i.e., the vaccination capacity of the site. Finally, residents must be allocated to the opened sites to get vaccinated. Thus, the problem corresponds to a capacitated bi-objective multi-period vaccination planning problem (MVPP) with replenishment. This new problem generalizes the multi-period facility location problem (MFLP) since one has to determine the capacity and replenishment of each facility at each period. To the best of our knowledge, the capacitated bi-objective MFLP with replenishment has not been studied, which motivates us to develop new models and algorithms for this critical problem.

In this paper, we study the MVPP that simultaneously optimizes two objectives: (i) minimizing the total operational cost, including the fixed vaccination site opening cost, capacity selection cost, replenishment cost, and inventory cost; (ii) minimizing the total travel distance of the population (to improve recipients’ convenience, i.e., service level). The decisions to be made, for each period, are: 1) which vaccination sites to open; 2) how many vaccination stations to launch at each opened site; 3) the replenishment quantity at each site; and 4) how to assign the population to the opened sites.

We develop a mathematical programming model and design two algorithms to solve it. The first one is a weighted-sum method, and the second is based on an -constraint algorithm. Both methods provide a representative set of optimal Pareto solutions. To handle the practical case, we further develop a tailored genetic algorithm with new features designed specifically for the problem at hand.

The remainder of this paper is organized as follows. Section 2 overviews related literature. Section 3 describes the studied problem in detail and presents a MILP model. Section 4 is devoted to the proposed algorithms. Section 5 presents the computational study, the analysis of the results, and the sensitivity analysis. Finally, Section 6 concludes the study and indicates some research perspectives.

Literature review

COVID-19, which is an infectious disease resulting from a previous variant of the coronavirus [31], has created a new type of SC risk because it can render numerous components of an SC inefficient or unusable for an uncertain time duration [32]. Therefore, since the COVID-19 virus outbreak and the beginning of the associated pandemic, scholars have published their studies on various SC-related issues raised by COVID-19, such as the negative impact on supply chain efficiency and performance [4], viable supply chain model [7], [8], ripple effect in SCs [33], [34], and re-configurable SCs [35], [36]. However, few studies focus on the vaccination planning problem. Vaccines are essential for planning and deploying recovery operations in an epidemic. Next, we first review the vaccination planning problem. We then survey the related MFLP. Finally, we identify the research gaps based on our review.

Vaccination planning problem

Research on vaccination planning problems generally focuses on the location of vaccination sites and the vaccine distribution and inventory decisions. The location of the vaccine distribution center has been widely studied. Li et al. [29] propose a multi-objective mixed-integer non-linear programming model to help the CDC determine the locations of vaccination stations. The model jointly considers travel distance and operational cost within a single planning period. Rastegar et al. [37] present a MILP model for the location-inventory problem to provide equitable influenza vaccine distribution among different groups of people with varying priorities over multiple periods. Bertsimas et al. [38] present a bilinear and non-convex model to optimize the location of vaccination sites and subsequent vaccine allocation. Results suggest that the locations of vaccination sites can have a massive impact on the effectiveness of the vaccination campaign.

In terms of the vaccine distribution problem, Proano et al. [39] formulate a mixed-integer non-linear programming model and propose a constructive heuristic to determine the best combination of vaccines and their prices. Duijzer et al. [40] analyze the optimal allocation of a vaccine to avoid infection of the maximum number of people. They define a unique dose-optimal vaccination fraction that maximizes the health benefits per dose of vaccine in a population. Enayati et al. [41] study the optimal influenza vaccine distribution in a heterogeneous population to minimize the number of vaccine doses distributed to extinguish an emerging outbreak in its early stages. Moreover, an equity constraint was proposed to help public health authorities consider fairness when making vaccine distribution decisions. Boeck et al. [42] review the vaccine distribution chains in low- and middle-income countries. Yang et al. [43] focus on the vaccine distribution networks in these countries. They formulate the network design problem as a mixed-integer program (MIP) and present a new algorithm for typical problems that can hardly be solved using commercial MIP software.

Regarding the vaccine inventory, Jacobson et al. [44] propose a stochastic inventory model to capture production interruption and evaluate the impact of pediatric vaccine inventory level on vaccination coverage. Samii et al. [45] focus on an inventory rationing problem with a single vaccination period and two demand classes. They utilize service level and fill rate expressions to address the two interdependent problems of selecting an appropriate allocation mechanism and determining the optimal reserved quantity. Pambudi et al. [46] discussed how vaccine cold chain management and cold storage technology can help address the challenges of vaccination programs. Specifically, it examines several ways for conserving vaccines in liquid or frozen form to ensure that they do not get destroyed during distribution from manufacturing plants. Yarmand et al. [47] propose a two-phase vaccination policy to contain an epidemic where the first phase considers the early stages of the epidemic, and the second phase the middle of the outbreak after observing the initial outcome. The two-phase vaccination policy provides a lower infection rate and considerable cost savings. Guo and Cao [48] build a two-stage analytical framework to analyze vaccination hesitancy, showing that vaccination decisions are affected by reference point formation and updating. Wong et al. [49] explored Hong Kong adolescents attitudes towards the COVID-19 vaccination, pointing out that vaccine hesitancy is a key barrier to herd immunity, and that the vaccine’s safety and efficacy, as well as the risk of infection, were the main concerns. Macintyre et al. [50] analyze the impact of various COVID-19 vaccine strategies on COVID-19 case numbers and mortality under a limited supply scenario, as well as varied vaccine efficacy and speed of vaccination for mass vaccination. Martonosi et al. [51] address the COVID-19 pricing issue and use optimization and game theoretic approaches to model the vaccine market as a duopoly with two manufacturers, Pfizer-BioNTech and Moderna.

Research on the vaccine supply chain generally focuses on the facility location, distribution, and inventory of vaccines to minimize cost. To our knowledge, the capacitated multi-period vaccination site location with replenishment has not been addressed in this stream of literature.

Multi-period facility location problem

As mentioned, our problem is a generalization of the MFLP, which is an extension of the facility location problem (FLP). Generally, an FLP involves a location decision for optimally locating facilities and an allocation decision to allocate customers to facilities. Opening new facilities involve time and capital investment, and it is one of the most critical strategic decisions for any institution. Hence, the World Bank and various government projects have extensively used facility location models, such as for schools [52], or healthcare [53], [54]. Depending on the number of periods in the planning horizon in which location and allocation decisions are to be made, FLPs can be categorized as single-period and multi-period [55]. Most of the early work was done for the static or single-period case. Once the facilities have been optimally located, they are assumed fixed regardless of how demand and costs may change in future periods. In practice, however, the demand is time-varying, and other parameters such as operation and transportation costs may change over time. Therefore, it might be better to relocate the facilities, which promotes the development of multi-period location models, in which the optimal location of facilities is determined by minimizing the total relocation cost.

MFLPs have been widely studied after the initial works of [56], [57], [58]. Traditionally, the objective of an MFLP is to determine the spatial distribution of facilities at each period of a finite planning horizon to minimize the total fixed and variable costs for satisfying demands over time [59], [60]. Recently, [55], [61] discuss fundamental modeling aspects and address several variants of the MFLP. In the context of our vaccination planning problem, demand is highly dynamic over times, which makes the MFLP relevant. In most current studies, the service capacity of a facility is assumed to be fixed and given once opened. However, this may not be reasonable for highly dynamic demand environments, where decision-makers can flexibly determine the facility size to save cost and meet the dynamic demand. Some applications have dealt with facilities that vary their capacities over time [62], [63], [64]. In the forestry industry, this has been called modular capacities, where capacity can be added or removed from a facility in fixed-size blocks, called modules [65]. This problem has received some attention with mathematical programming [66], and heuristics based on Lagrangean relaxation [67] and evolutionary search [68]. In addition, current studies on the MFLP do not involve replenishment decisions, which is the case of the investigated problem.

There is a potential trade-off between the total operational cost and the service level that is often reflected by the travel distance to receive a vaccine. Therefore, decision-makers may expect to open fewer vaccination sites to save costs. However, they also hope to provide a high service level by locating more vaccination sites to facilitate access. This, in turn, incurs a higher operational cost. Thus, decision-makers face a multi-objective optimization problem. However, the cost and service level trade-off has rarely been studied in the above-reviewed literature.

Contributions

In summary, the investigated problem is new, and the current models and solution algorithms cannot be directly applied to solve it. To this end, we study the capacitated bi-objective MVPP with replenishment. Our paper makes the following contributions:

We present a new formulation for the bi-objective MVPP to simultaneously minimize the operational costs and travel distance (maximize service level) and discuss its practical application in the context of allocation of limited budget and vaccine availability to a vaccination planning platform.

We solve the location selection and population allocation to vaccination sites and make plans for a series of practical operations such as work scheduling, service capacity setting, vaccine supply, and storage. The solution of this novel bi-objective vaccination planning problem can provide a practical decision support tool to massive vaccination during the outbreak of epidemics like COVID-19 with detailed operational directives for successful applications.

We present a genetic algorithm with several refinements that exploit the problem’s underlying structure to very efficiently solve instances of practical sizes.

Problem description and formulation

This section first elaborates on the studied MVPP with replenishment and then proposes a MILP model. Consider a set of vaccination sites, each site located at (). A set of recipients are to be vaccinated within the planning horizon . Each recipient has a specific location () and an appointment date indicating the day the recipient expects to get vaccinated. The travel distance between recipient and vaccination site is . Each vaccination site opened on day has a fixed setup cost and a variable cost depending on the number of vaccination stations deployed. A fixed replenishment cost is incurred whenever vaccination site is replenished. The unit inventory cost at vaccination site is . Each vaccination site has a maximum replenishment quantity , inventory capacity , and a maximum number of stations that can be deployed. The maximum number of recipients that a vaccination station can serve each day (service capacity) is denoted by . We assume that replenishments occur at the beginning of the day before the vaccination service starts. The objective is to simultaneously minimize the total operational cost and the total travel distance of all recipients by optimizing: 1) the selection of the opened sites for each day; 2) the replenishment date and quantity for each site; 3) the number of stations launched at each opened site; 4) the inventory quantity at each site at each day, and; 5) the allocation of recipients to the opened sites at each day.

We model the problem with the following binary variables: is equal to 1 if recipient is served at vaccination site ; is equal to 1 if site is opened on day ; is equal to 1 if site is replenished on day . We also use variables to represent the replenishment quantity of vaccine at on day ; the number of vaccination stations at site on day ; the inventory quantity at site at the end of day . Based on the above description, we define the following notation.

We next present a bi-objective MILP model () for the studied MVPP as follows.

The objective function (1) minimizes the total operational cost, which contains the fixed setup cost for opening vaccination sites, variable cost for launching vaccination stations, fixed replenishment cost, and the inventory cost. The objective function (2) minimizes the travel distance of all recipients to get vaccinated. Constraints (3) guarantee that a recipient is served by exactly one vaccination site. Constraints (4) limit the maximum number of recipients that a vaccination site can serve per day. Constraints (5) restrict the maximum number of vaccination stations that can be launched at an opened site. Constraints (6) correspond to the inventory flow balance at each vaccination site, and constraints (7) limit their maximum inventory. Constraints (8) indicate that the maximum replenishment quantity at a vaccination site on each day cannot exceed a given threshold. Constraints (9)–(12) define the domains of variables.

The considered MVPP contains a capacitated MFLP, which generalizes the uncapacitated MFLP. The MVPP is an NP-hard problem since the latter is known to be NP-hard [60].

Solution methods

To solve the considered MVPP, we first apply the widely-used weighted-sum (WS) approach and -constraint (EC) method. The WS method transforms the original bi-objective problem into a single objective one by associating each objective function with a weighting coefficient. It then solves a series of single-objective MVPPs by varying the objectives’ weights. The WS method is extensively used because it is simple to understand and easy to implement [69], [70]. However, one of the disadvantages of the WS is its inability to find specific Pareto optimal solutions in the case of a non-convex objective space [71], [72]. The EC method is often preferred over the WS method in solving bi-objective optimization problems because the EC can find solutions also in the non-convex regions [73]. In general, the EC optimizes one objective and transforms the other one into a constraint bounded by . It then solves a series of single-objective problems by varying the value of . Given proper increments of , the EC is guaranteed to find the entire set of Pareto optimal solutions for a general multi-objective problem [74]. Moreover, [75] shows that an optimal solution identified via the EC is guaranteed to be Pareto optimal if the constraints representing the objectives are binding and the solution is feasible. Hence, the EC is more often utilized against the WS in recent research on multi-objective location-allocation problems [76].

In this paper, the WS and EC methods are developed based on the proposed MILP model. Our preliminary results show that they work well for small-sized instances, but they have difficulty in solving large-sized instances due to the NP-hardness of the studied problem. To this end, we further develop a non-dominated sorting genetic algorithm (NSGA-II) to provide a set of approximate Pareto solutions for large-sized instances. NSGA-II is one of the most popular multi-objective optimization algorithms with three special characteristics: fast non-dominated sorting approach, fast crowded distance estimation procedure, and simple crowded comparison operator [77]. Besides, NSGA-II is the most widely applied multi-objective evolutionary algorithm as observed in the existing literature [78], [79], [80], [81], [82]. We next present the WS, EC, and NSGA-II algorithms in sequence.

The weighted-sum approach

The WS method combines the two objectives to form a single-objective optimization problem. The idea is to associate a weighting coefficient with each objective function and minimize the weighted sum of the two objectives. In detail, we define two parameters and which represent the preference of decision-makers on each objective, where , and . Since the two objectives and have different scales, we normalize them to eliminate the effects of dimensions [83]. The normalized objectives are calculated by the following equation: [84], [85]. We follow this literature and normalize the objectives into the interval [0,1]. Then, we denote the single-objective MVPP with the WS method as S-MVPP (WS), which is shown as follows:The values of , , , and are obtained by exactly solving the following mono-objective problems: where is the vector of variables and is the solution space of .

We can obtain a set of Pareto solutions by solving a series of S-MVPP (WS) with different combinations of and values. Initially, the value of is set to . Then the value of is increased by at each iteration. Algorithm 1 reports the main steps of the WS method.

Algorithm 1

Weighted-sum method.

The -constraint approach

We next apply the EC method to solve the bi-objective MVPP. The EC method is widely used in solving multi-objective optimization problems [86]. Its basic idea is to transform the original multi-objective problem into a single-objective one, in which only one objective is directly optimized. Other objectives are transformed into constraints. Compared to the WS method, the EC has three typical characteristics [87]. First, it does not need to normalize different objectives with different units or scales. Second, it can flexibly obtain the expected number of Pareto solutions and manage to find non-extreme solutions. Third, it can provide a good trade-off between solution quality and computation time. For the studied MVPP, we transform it into a single-objective MVPP by optimizing the second objective, with the first one as a constraint. We denote the single-objective MVPP with an EC as S-MVPP (), which is shown as follows:

Constraint (19) represents the EC that restricts the value of the first objective. The S-MVPP () is then solved to minimize the total travel distance (second objective) with the total operational cost (first objective) being bounded by . A set of Pareto solutions can be obtained by solving a series of S-MVPP () with different values of , whose value is bounded by an upper limit and a lower limit . Then, the S-MVPP is solved iteratively by defining a step size , where is the expected number of Pareto solutions. Initially, takes the value of . Then, it is reduced by the step size at each iteration. In the iteration, where , let the value of the first objective be represented by . Then is set to be for the iteration. We solve a series of S-MVPP () as long as the value of . Finally, a set of Pareto solutions are obtained. The detailed steps of the EC method are shown in Algorithm 2 . Note that the selection of the principal objective may affect the method’s performance. We also tested setting the first objective as the principal one, but the results were inferior to those of the presented setting. The computational results of Section 5.3.1 support this choice.

Algorithm 2

-constraint method.

Non-Dominated Sorting Genetic Algorithm II (NSGA-II) approach

The WS and EC methods must solve a series of MILP models during the solution process. Our results show that both methods lose efficiency in solving practical-sized instances. To this end, we tailor the NSGA-II to tackle large-sized instances. The NSGA-II has been successfully applied in the literature to solve multi-objective optimization problems [82]. The detailed process of NSGA-II can be found in [78]. It can quickly obtain a set of approximate Pareto solutions through crossover and mutation operators based on an initial population. In this process, the generation of the initial population plays a key role, and the performance of NSGA-II is strongly dependent on its quality. For our problem, it is non-trivial to construct a good feasible solution. In particular, it is challenging to appropriately assign many recipients to vaccination sites, given the various limitations, e.g., vaccination capacity, inventory capacity, and replenishment capacity. To tackle this difficulty, we first design a new assignment strategy. Each recipient gets a score based on the distance difference between their nearest and second nearest vaccination sites. We then apply a dynamic programming method to determine the replenishment time and quantity, and the inventory quantity. With the above two novel features, our tailored NSGA-II obtains good quality Pareto front for the studied bi-objective MVPP. The flow chart of NSGA-II is shown in Fig. 1 , and dummyTXdummy- we describe each step in detail next.

Fig. 1

Flow chart of NSGA-II.

Step 1: Chromosome Encoding. Each individual (solution) is represented by a chromosome that comprises genes whose number is equal to the number of variables in the solution. For our problem, each chromosome consists of three parts. The first part represents the value of variables which takes a value of 1 if the corresponding vaccination site is opened. The second part corresponds to the recipient allocation decisions, denoted as , indicating the assignment of recipients to a site. The third part indicates the replenish quantity . The number of genes is . Fig. 2 shows a chromosome with =6 recipients, =3 vaccination sites, and =2 days. Vaccination sites 1 and 3 open on the first day, and sites 2 and 3 open on the second day. The six recipients are allocated to vaccination sites 3, 2, 1, 3, 3, and 1. Site 1 receives two doses of the vaccine on the first day, site 2 receives one dose on the second day, and site 3 receives three doses on the first day.

Fig. 2

Chromosome illustration.

Step 2: Population initialization. A set of individuals constitutes the initial population. We develop a new heuristic to generate the initial population, as shown in Algorithm 3 . The new heuristics consists of three phases. In the first one, the opened vaccination sites are determined. Then, the set of recipients that each opened vaccination site serves on each day is determined using an improved assignment strategy. The third step calculates the replenishment quantity using a dynamic programming method. Next, we elaborate on each phase to show how each individual is generated. This process is repeatedly performed to generate solutions corresponding to the initial population, where the random elements of Phase 1 ensure diversification.

Algorithm 3

Chromosome initialization.

1) Determine the opened vaccination sites on each day (). Let be the number of recipients to be vaccinated at day , be the service capacity for vaccination site , which is the minimum between the vaccination capacity and replenishment capacity. In Algorithm 4 , we generate two extreme solutions, one where all vaccination sites open each day (lines 1–2), and another where the fewest vaccination sites are opened (lines 3–12).

Algorithm 4

Determine the opened vaccination sites.

2) Allocate recipients to opened vaccination sites (). Let be the number of recipients assigned to site at day . In Algorithm 5 , for each recipient that is not allocated, we calculate the distance difference between its nearest and second nearest opened vaccination sites with remaining capacity. Next, for each day and each opened vaccination site , we define a set with all the unallocated recipients whose nearest site is and the expected vaccination date is . We also define another set that contains all the unallocated recipients whose nearest and second nearest sites are and , where , and whose expected vaccination date is . Let , where all recipients in are sorted in non-increasing order of . Then we compute , and assign the first recipients to . The values of the first recipients in the set is set to , and . The rationale of such an assignment rule is to prioritize recipients assigned to a site that is generally far from their location if not assigned to their nearest site. This assignment strategy provides better performance than sorting all recipients in non-decreasing order of their nearest sites. The computational results presented in Section 5.3.2 support the advantage of our strategy.

Algorithm 5

Assign recipients with an improved strategy.

3) Determine the replenishment variables , , and inventory quantity variables . Once th site opening variables and recipient assignment variables are determined, the remaining problem is to compute the optimal replenishment and inventory at each opened site . We recall that is the number of recipients to be served at the vaccination site on the day , is the fixed replenishment cost at site , and is the unit inventory cost at site . Then, for each site we compute the optimal value of the replenishment date variables , replenishment quantity variables , and inventory variables by solving the following MILP model:

It is obvious that the above discrete economic ordering problem at each opened site satisfies the Wagner-Whitin condition [88]. That is to say, the condition is satisfied, indicating the replenishment takes place only when the inventory is completely consumed. To this end, we define parameter as the replenishment and inventory costs incurred between periods and if all vaccines needed during this interval are received at time .

The first term denotes the replenishment cost at period . The second term indicates the inventory costs from periods . The inventory quantity at period is equal to the sum of vaccines needed during periods to .

Let be the minimum cost from the first day to day when the inventory is zero at the end of day . We can use the following state transition function:where =0. The objective value can be obtained by calculating .

The detailed implementation process of the proposed dynamic programming (DP) is shown in Algorithm 6 . This procedure provides the values for variables , , and . Finally, variables can be computed as .

Algorithm 6

Determine , , and with a DP method.

Step 3: Perform the non-dominated sorting to rank the initial population and create Pareto fronts. Specifically, each individual gets a fitness value equal to its non-domination level. The first front is a completely non-dominant set, and other fronts are only dominated by the individuals in precedent fronts. In addition, the crowding distance for each individual is calculated, which is then used to measure its closeness to neighbors. Generally, a larger average crowding distance denotes a better population diversity.

Step 4: Parents are selected using a tournament selection operator based on their fitness value and crowding distance. Individuals with better fitness or larger crowding distance are more likely to be selected.

Step 5: Generate offsprings by crossover and mutation operations. Specifically, we select the first part (i.e., ) as the basic chromosome because once changes, and change deterministically. The crossover and mutation operations are carried out on the basic chromosome. After each crossover and mutation operation, Algorithms 5 and 6 are used to generate new solutions (offspring population) based on the new values of variables . The offspring is then added to the population, and a selection operation is performed based on non-domination to select the individuals of the next generation. Finally, elitism is ensured by selecting the best individuals.

Step 6: Stopping criterion. If the maximum number of iterations is reached, stop and output a set of non-dominated Pareto solutions; otherwise, return to Step 4.

Computational study

In this section, we evaluate the performance of the proposed algorithms by conducting computational experiments. We first present a real-world case study based on the vaccination program in Tongzhou District, Beijing City, China. We then perform numerical experiments on 64 random instances with up to 200,000 recipients, 50 vaccination sites, and 10 days. The random instances and results are available online at https://www.dmu-yantongli.com/instances.

The proposed weighted-sum and -constraint algorithms are coded in C++ linked with CPLEX 12.10, and the NSGA-II method is implemented in MATLAB R2020b. All runs are performed on a PC with Core i5 at 1.6 GHz and 8 GB RAM. The parameters used in the algorithms are stated as follows. For the weighted-sum method, the value of is set to 0.05. For the -constraint method, and is set to 20. For both methods, the total time limit and that for a single iteration are set to 3600 and 600 seconds, respectively, for the small- and medium-sized instances. The two parameters are set to 7200 and 1200 seconds, respectively, for large-sized instances. In terms of the NSGA-II, the population size , the number of generations , and the time limit are set to 600, 200, and 3600 seconds, respectively, for the small- and medium-sized instances. These three parameters are set to 300, 100, and 7200, for large-sized instances.

Case study of Beijing Tongzhou District

This section presents a case study from the Tongzhou district of Beijing City, China, to validate the proposed model and solution method. The case corresponds to a vaccination program covering 12 days and 202,595 recipients. The Center for Disease Control of Tongzhou District is responsible for the mass vaccination of all the residents in the district. The vaccination program consists of 17 vaccination sites, and the vaccination capacity for each site per day is listed in Table 1 , along with their locations, visualized in Fig. 3 .

Table 1

Information of vaccination sites in Tongzhou district.

No.	Name	Latitude	Longitude	Capacity/d	# of recipients
1	Yangzhuang Street	39.91143	116.63804	840	5930
2	Beiyuan Street	39.91414	116.65996	780	8111
3	People Mall	39.91425	116.67745	2960	31546
4	Tongyun Street	39.91899	116.70526	700	4225
5	Jiukeshu Street	39.90094	116.65871	1120	9487
6	Yuqiao Street	39.89856	116.69962	1120	12340
7	Yongshun Town	39.95379	116.68951	1400	14200
8	Songzhuang Town	39.95223	116.73489	1680	18195
9	Liyuan Town	39.87097	116.65589	1950	18249
10	Zhangjiawan Town	39.86281	116.74197	960	10215
11	Lucheng Town	39.86046	116.82250	840	9080
12	Taihu Town	39.82754	116.64454	1560	13190
13	Xiji Town	39.80478	116.83826	840	8830
14	Huoxian Town	39.78210	116.80376	720	7290
15	Majuqiao Town	39.74870	116.57271	1680	13583
16	Yujiawuxiang Town	39.72415	116.71123	1120	7432
17	Yongledian Town	39.72156	116.80356	1120	10692
Total				21390	202595

Fig. 3

The location of vaccination sites in Tongzhou District.

We obtained the vaccination data for 12 days from Mar 19, 2021, to Mar 30, 2021. The total number of vaccinated recipients is 202,595. We estimate the location of each recipient based on the total number of residents in each community and the vaccination rate. In particular, let the total number of residents in Tongzhou district be , and the number of permanent residents in each community be , where and is the number of communities. We then calculate the vaccination rate of the Tongzhou district, which is denoted by . Next, we estimate the vaccinated population in each community by . The obtained number of recipients at each community is shown in the last column of Table 1. We generate a random coordinate within the corresponding administrative subdistrict for each vaccinated recipient.

The average working time is 8 h/d, and each vaccination takes about five minutes, indicating that the service capacity of a vaccination station is 100/d. The maximum number of vaccination stations that can be launched at each vaccination site varies from 7–30. We assume that the expected service day for recipients is uniformly distributed along the planning horizon. The fixed setup cost for vaccination site is related to its location, human and other resources, and falls in the interval [2000, 3000]. The fixed replenishment cost at vaccination site is mainly dependent on the transportation cost of vaccines and lies in the interval [1000, 2000]. The cost of a vaccination station at vaccination site per day is set to 600, which corresponds to the salary for a medical professional. The unit inventory cost at the vaccination site per day is randomly generated from the interval [0.2, 0.5]. The maximum replenishment and inventory capacity at vaccination site on each day are randomly generated as , where is randomly generated from 1, 2, 3. A summary of these values is presented in Table 2 .

Table 2

Values for the model parameters.

PAR.	Description	Value
|J|	# of recipients	202595
|K|	# of sites	17
|T|	# of days	12
f_k	Setup cost	[2380 2795 2418 2688 2272 2872 2823 2966 2402 2692 2553 2183 2007 2553 2398 2724 2668]
e_k	Replenishment cost	[1989 1358 1910 1622 1080 1334 1518 1067 1767 1539 1622 1942 1771 1070 1827 1694 1047]
h_k	Inventory holding cost	[0.30 0.23 0.49 0.33 0.38 0.29 0.47 0.45 0.42 0.43 0.25 0.22 0.33 0.26 0.44 0.33 0.31]
m_k	Maximum number of vaccination stations	[9 8 30 7 12 12 14 17 20 10 9 16 9 8 17 12 12]
O_k	Maximum replenishment quantity	[16883 16883 33766 33766 33766 50649 50649 16883 50649 50649 50649 33766 50649 33766 50649 33766 16883 ]
V_k	Maximum inventory quantity	[16883 33766 50649 50649 50649 33766 50649 16883 33766 33766 33766 16883 50649 50649 16883 33766 50649]

Results and analysis

We solve the case study using the proposed NSGA-II method. A time limit of 2 hours is imposed. The size of the initial population and the number of iterations are set to 300 and 100, respectively. The obtained Pareto solutions (blue triangles) are shown in Fig. 4 .

Fig. 4

Pareto front and solutions obtained by experience-based methods.

Decision-makers can select a solution from the obtained Pareto solutions based on their preference. In this case, we apply the fuzzy logic decision method to help decision-makers choose their preferred solution [89]. We consider three combinations of weights and for the two objectives, each representing a specific preference on the two objectives. In scenario 1, the weight coefficients and are set to 0.2 and 0.8, respectively, indicating that the decision-maker assigns a higher weight to the travel distance and tends to focus on the service quality. Scenario 2 allocates equal weights to both objectives (, ), indicating that the decision-maker believes operational cost and service quality are equally important. In Scenario 3, and are set to 0.8 and 0.2, which indicates that the decision-maker views the operational cost as a more critical criterion and tends to make a cost-effective plan. We apply the fuzzy logic decision method to select the preferred solutions for the three scenarios. The selected solutions for the three scenarios are shown in red “+" from top-down in Fig. 4, from which we can see that the method selects three distinct solutions for the three scenarios.

We then compare the obtained Pareto front to the solutions obtained by two sequential heuristics adopted in practice. The heuristics used in practice consist of scenarios where decision-makers make decisions sequentially based on their experiences. They follow two greedy rules: 1) each recipient is assigned to the vaccination site in his administrative subdistrict (denoted by heuristic solution 1); 2) each recipient is assigned to its nearest vaccination site (denoted by heuristic solution 2). The obtained solutions of the two heuristics are presented in Fig. 4. We see from Fig. 4 that for both heuristic solutions, we can find better solutions in the Pareto front that dominate them. In detail, we find the following.

First, heuristic solution 1 is dominated by all preferred solutions, particularly by the preferred solution with = 0.5. The reduced operational cost and travel distance can be as high as 1.67E+05 and 2.66E+08, respectively. Therefore, the proposed model can significantly reduce the operational costs (with an average of 9.3%) and the total travel distance (with an average of 36.6%) when compared with real-world practice.

Second, regarding heuristic solution 2, it is dominated by the preferred solution with = 0.2, with savings on operational cost and travel distance of 3190 and 2.63E+07, respectively.

These results indicate that our integrated decision model and solution methods are effective and outperform the sequential heuristics that are adopted in practice. More specifically, our model shows superior ability to improve resource utilization, which is especially valuable for carrying out a mass vaccination campaign. Besides, the service level is improved, and the time to reach mass vaccination sites is reduced, which helps mitigate the risks of virus spread. Thus, our tools can help decision-makers reduce operational costs and improve service levels during mass vaccination.

The preferred solutions corresponding to the three scenarios clearly show the trade-offs between the two optimized objectives, i.e., the operational cost and travel distance. To better understand the cost components of the three solutions, we visualize in Fig. 5 their objective function values, the number of opened vaccination sites, the number of launched vaccination stations, and the capacity utilization rate. The capacity utilization rate is calculated as .

Fig. 5

Results obtained from the trade-off analysis.

The analysis of these detailed solutions allows us to derive the following insights.

1) As the weight increases from 0.2 to 0.8, the operational cost decreases, and the travel distance increases. This is expected since a larger weight for the operational cost objective might force the model to find a cost-effective vaccination plan. This leads to a situation where recipients must travel a longer distance to get vaccinated. Besides, there is an apparent conflict between these two objectives: promoting one objective is achieved at the expense of deteriorating the other.

2) Figs. 5(b) and 5(c) show that when the value of increases, the total number of opened vaccination sites and stations decreases, and the proportion of closed vaccination sites increases. Specifically, when a larger weight is allocated to the travel distance objective (i.e., and ), almost all the vaccination sites are opened for vaccination each day. When a larger weight is allocated to the operational cost objective (i.e., and ), the total number of opened vaccination sites decreases from 200 to 151, increasing the proportion of closed vaccination stations from 2% to 27%. Similarly, when increases from 0.2 to 0.8, the number of vaccination stations used decreases from 2086 to 2041, as shown in Fig. 5(c).

3) In Fig. 5(d), we see that as the weight of the operational cost objective increases from 0.2 to 0.8, the service capacity utilization rate increases from 78.50% to 96.64%. This is aligned with the previous finding, when most people are assigned to fewer sites, increasing utilization.

In summary, this trade-off analysis provides managerial insights on the bi-objective MVPP. The proposed decision framework can better balance the workload at each vaccination site owing to the scientific location, assignment, and service setting. Resource wasting can be eliminated, which dramatically helps improve the vaccination operation performance. Our method provides a set of Pareto solutions from which decision-makers can choose flexibly based on their preference and constraints. This information is beneficial for decision-makers when they face such a complex problem in which the operational cost and the service level must be simultaneously optimized. Accordingly, the situation shows that to fight the COVID-19 outbreak, one can rely on advanced solution methodologies. Our paper contributes with an efficient procedure to help policymakers.

Sensitivity analysis

To analyze the impact of service capacity at each vaccination site, we expand it by increasing the number of potential stations. The average value of for the above Pareto front is 13, which is set as a basis. We then consider three scenarios where the maximum number of vaccination stations is set to 20, 25, and 30. These scenarios are then solved by the same NSGA-II method, and the obtained Pareto fronts for the three scenarios are presented in Fig. 6 .

Fig. 6

Effect of changing service capacity.

From the results shown in Fig. 6, we see that as the service capacity increases from 13 to 30, the obtained Pareto front becomes better and that the larger difference happens concerning the cost, not distance. This indicates that decision-makers can significantly reduce operational costs and improve service levels by enlarging the service capacity limits.

Results analysis on randomly generated instances

We generate 64 instances in three groups. The first one contains 20 small-sized instances with up to 1000 recipients, 10 vaccination sites, and 10 days, i.e., , , . The second group contains 20 medium-sized instances with up to 10,000 recipients, 10 vaccination sites, and 10 days, i.e., , , . The third group contains 24 large-sized instances with up to 200,000 recipients, 50 vaccination sites, and 10 days, i.e., , , . The coordinates of the vaccination sites and of the recipients are randomly generated from [1, 200]. The distance between vaccination site and recipient is calculated as . The expected service day for recipient is randomly generated from the interval [1, ]. Finally, we generate an instance for each combination of , and , totaling 64 instances with up to 200,000 recipients, 50 vaccination sites, and 10 scheduling days.

These instances are solved with each of the three methods.

Performance indicators

We define several indicators for evaluating the performance of different solution methods for multi-objective optimization problems. To this end, we first define a reference Pareto solution set which contains all solutions provided by the three developed methods, and dominated solutions are excluded. Then we present four performance metrics, i.e., the number of non-dominated solutions, hypervolume ratio, average -dominance, and computation time, widely adopted in evaluating multi-objective optimization solutions [89], [90]. The first and last indicators are trivial, while the second and third indicators are introduced by [91].

The number of non-dominated solutions indicates the number of Pareto solutions in the approximate Pareto solution set excluding those dominated by the solutions in the referenced set . This indicator can be viewed as the contribution of the current algorithm to the reference set. A large number of non-dominated solutions indicates a better performance.

The hypervolume ratio represents the ability of coverage of a Pareto solution set. Let be the hypervolume of the reference set, which is calculated based on the nadir point . It represents the objective space covered by the set . The hypervolume ratio of set can be calculated as . A large hypervolume ratio indicates the solution set has better coverage, which suggests that the corresponding algorithm is better.

The -dominance indicator measures the average distance between the reference solution set and the approximate solution set . The -dominance value of a solution in the approximate solution set is calculated as follows:

The average -dominance indicator for each approximate solution set is calculated as:

The average -dominance value is always greater than or equal to 1, and a smaller value indicates the approximate solution set is closer to the reference solution set, i.e., a better performance of the corresponding algorithm.

The comparison results on the 64 randomly generated instances are presented in Tables 3 , 4 , and 5 , where “Weighted-sum”, “-constraint”, and “NSGA-II” represent the three proposed solution methods, respectively. Column “#” denotes the sequence number of the instances; “#Nd” indicates the number of non-dominated solutions in the obtained approximate Pareto solution set; “Hv” indicates the value of hypervolume ratio; “” represents the value of average -dominance, and “Time (s)” is the computational time in seconds. Recall that WS refers to the weighted-sum method and EC to the -constraint algorithm.

Table 3

Comparison results for small-sized instances.

#	|J|	|K|	|T|	Weighted-sum				ϵ-constraint				NSGA-II
				#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)
1			5	13	0.958	1.023	58	14	0.977	1.011	24	8	0.999	1.001	818
2		5	10	16	0.967	1.012	360	9	0.968	1.009	336	68	0.997	1.001	908
3			5	16	0.977	1.016	70	18	0.977	1.014	51	77	0.998	1.001	673
4	200	10	10	18	0.991	1.008	1411	20	0.978	1.008	628	117	0.980	1.005	1514
5			5	13	0.947	1.029	215	17	0.968	1.012	267	16	0.996	1.003	973
6		5	10	13	0.930	1.025	1388	19	0.981	1.009	2767	91	0.994	1.002	1646
7			5	16	0.981	1.022	1286	13	0.977	1.017	1070	145	0.998	1.001	1477
8	400	10	10	19	0.990	1.014	3591	18	0.978	1.010	813	195	0.987	1.003	2258
9			5	12	0.944	1.019	174	12	0.969	1.009	248	19	0.997	1.001	1662
10		5	10	8	0.769	1.052	3642	12	0.996	1.005	2833	40	0.981	1.004	2085
11			5	15	0.977	1.024	3252	17	0.978	1.015	1472	144	0.996	1.002	2313
12	600	10	10	16	0.967	1.016	3275	21	0.982	1.009	3178	167	0.990	1.003	2966
13			5	10	0.963	1.014	244	13	0.974	1.006	756	14	0.992	1.002	1740
14		5	10	15	0.978	1.026	2796	15	0.989	1.008	2737	91	0.994	1.001	1640
15			5	17	0.957	1.031	3969	16	0.972	1.018	2861	171	0.999	1.001	1961
16	800	10	10	12	0.874	1.068	4131	21	0.989	1.010	3644	165	0.995	1.003	3268
17			5	15	0.989	1.014	1956	18	0.988	1.006	3367	6	0.913	1.019	2393
18		5	10	14	0.927	1.023	3615	17	0.942	1.017	4123	119	0.995	1.002	2829
19			5	15	0.975	1.020	4108	15	0.966	1.015	2959	79	0.988	1.008	1668
20	1000	10	10	12	0.926	1.038	4027	16	0.871	1.092	4126	266	0.997	1.002	3612
Average				14	0.949	1.025	2178	16	0.971	1.015	1913	100	0.989	1.003	1920

Table 4

Comparison results for medium-sized instances.

#	|J|	|K|	|T|	Weighted-sum				ϵ-constraint				NSGA-II
				#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)
21			5	15	0.977	1.016	3586	10	0.839	1.066	3901	15	0.983	1.007	3608
22		5	10	10	0.844	1.035	4147	6	0.388	1.232	3973	242	0.999	1.001	3621
23			5	8	0.718	1.177	3604	13	0.901	1.046	4196	100	0.999	1.001	3620
24	2000	10	10	11	0.825	1.117	4180	12	0.663	1.406	3831	191	0.997	1.002	3615
25			5	7	0.417	1.118	3610	12	0.708	1.046	3760	27	0.994	1.002	3613
26		5	10	9	0.770	1.105	3932	10	0.565	1.169	4055	69	0.996	1.002	3626
27			5	8	0.789	1.093	4193	13	0.844	1.084	3732	122	0.999	1.001	3614
28	3000	10	10	11	0.849	1.082	4148	10	0.558	1.443	3913	235	0.999	1.001	3619
29			5	6	0.516	1.094	3888	9	0.588	1.100	4012	27	0.981	1.005	3602
30		5	10	7	0.839	1.035	3685	13	0.754	1.065	4138	54	0.998	1.003	3624
31			5	8	0.782	1.075	3733	13	0.788	1.074	4060	154	0.999	1.002	3606
32	4000	10	10	7	0.761	1.098	3791	11	0.583	1.224	3733	175	0.998	1.001	3618
33			5	7	0.844	1.023	3755	15	0.973	1.009	4117	16	0.963	1.006	3609
34		5	10	7	0.533	1.117	3931	7	0.417	1.182	4115	96	0.997	1.002	3643
35			5	9	0.668	1.103	4074	9	0.603	1.196	3644	66	0.995	1.003	3628
36	5000	10	10	7	0.678	1.216	3662	8	0.431	1.583	3707	170	0.996	1.002	3615
37			5	6	0.984	1.004	3687	10	0.609	1.010	4149	0	0.839	1.009	3643
38		5	10	9	0.533	1.054	4036	8	0.623	1.062	3892	124	0.860	1.008	3708
39			5	8	0.788	1.028	3691	10	0.514	1.101	4135	26	0.996	1.002	3691
40	10,000	10	10	7	0.792	1.068	3615	6	0.308	1.338	4100	47	0.991	1.006	3690
Average				8	0.745	1.083	3847	10	0.633	1.172	3958	98	0.979	1.003	3631

Table 5

Comparison results for large-sized instances.

#	|J|	|K|	|T|	Weighted-sum				ϵ-constraint				NSGA-II
				#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)
41			5	4	0.402	1.404	7326	1	0.000	1.856	7200	129	0.999	1.000	5548
42		30	10	4	0.519	1.374	7289	1	0.000	2.740	7200	138	0.998	1.000	5578
43			5	5	0.380	1.521	7449	1	0.000	2.562	7200	148	0.998	1.001	4926
44	20,000	50	10	5	0.481	1.568	7407	1	0.000	3.882	7200	148	0.997	1.001	5663
45			5	5	0.423	1.270	7467	1	0.000	1.623	7200	113	0.998	1.001	6494
46		30	10	5	0.465	1.357	7324	1	0.000	2.127	7200	143	0.998	1.001	6395
47			5	3	0.313	1.455	7428	1	0.000	2.079	7200	142	0.997	1.000	6900
48	30,000	50	10	5	0.448	1.365	7424	1	0.000	2.812	7200	120	0.995	1.001	7218
49			5	4	0.388	1.197	7716	1	0.000	1.426	7200	132	0.999	1.000	7240
50		30	10	4	0.647	1.259	7528	1	0.000	1.980	7200	125	0.997	1.001	7286
51			5	3	0.232	1.407	7945	1	0.000	1.835	7200	123	0.999	1.000	7218
52	40,000	50	10	4	0.299	1.462	7669	1	0.000	2.407	7200	141	0.998	1.001	7255
53			5	5	0.384	1.141	7688	0	0.000	3.021	7200	67	0.995	1.001	7291
54		30	10	3	0.284	1.347	7483	1	0.000	1.788	7200	129	0.999	1.000	7319
55			5	5	0.695	1.188	8101	0	0.000	3.906	7200	97	0.996	1.001	7275
56	50,000	50	10	3	0.261	1.524	7882	–	–	–	7200	139	0.998	1.000	7279
57			5	1	0.190	1.136	7262	–	–	–	7200	53	1.000	1.000	7587
58		30	10	3	0.402	1.126	7404	–	–	–	7200	62	0.998	1.000	7353
59			5	3	0.246	1.193	7663	–	–	–	7200	45	1.000	1.000	7782
60	100,000	50	10	2	0.295	1.256	7670	–	–	–	7200	95	0.996	1.000	7308
61		30		0	0.212	1.187	8245	–	–	–	7200	23	1.000	1.000	7755
62	150,000	50	10	–	–	–	7200	–	–	–	7200	35	1.000	1.000	7587
63		30		–	–	–	7200	–	–	–	7200	11	1.000	1.000	8963
64	200,000	50	10	–	–	–	7200	–	–	–	7200	19	1.000	1.000	13858
Average				3	0.332	1.156	7540	1	0.000	1.502	7200	99	0.998	1.000	7295

Results for small-sized instances

Table 3 reports the results for small-sized instances with up to 1000 recipients. We see that the tailored WS and EC methods obtain on average 14 and 16 non-dominated Pareto solutions, significantly fewer than the 100 non-dominated solutions provided by NSGA-II. This implies that NSGA-II contributes the most to the reference Pareto solution set. The values of the hypervolume ratio obtained by WS and EC methods are 0.949 and 0.971, respectively. In contrast, this value reaches 0.989 for NSGA-II, indicating a better coverage in the objective space of NSGA-II, and again a better performance when compared to the exact methods. Similarly, WS and EC methods yield a large value of -dominance, while the -dominance value of NSGA-II is exceptionally close to one, meaning the Pareto solution provided by NSGA-II is close to the reference set. In terms of the computation time, NSGA-II and EC take almost the same time. NSGA-II obtains the best solution quality within one hour for all the instances. The above comparisons demonstrate that NSGA-II has the best performance regarding the three quality indicators.

Further observing the computational results obtained by WS and EC methods, we can see that EC obtains a larger value of hypervolume ratio than WS, and a smaller value of -dominance. Additionally, with less computation time, the number of non-dominated Pareto solutions obtained by the EC method is greater than that of the WS method. This indicates that the Pareto solutions provided by the EC method cover a larger space and are closer to the reference Pareto solution front, implying a better performance of the EC method against WS. To sum up, for these small-sized instances, although CPLEX can optimally solve each single-objective MVPP (), the NSGA-II method still performs the best for solving the set of small-sized instances under the bi-criteria analysis.

Results for medium-sized instances

Table 4 compares results for medium-sized instances with up to 10,000 recipients. We see that the tailored WS, EC, and NSGA-II obtain on average 8, 10, and 98 non-dominated Pareto solutions, respectively. This implies that NSGA-II still contributes the most to the reference Pareto front. Regarding the hypervolume ratio indicator, the WS method provides a value of 0.745 compared to only 0.633 of the EC method. The NSGA-II method still performs the best on this indicator, with a value of 0.979. Regarding the -dominance indicator, the values obtained by WS, EC, and NSGA-II are 1.083, 1.172, and 1.003, respectively, which indicate that the WS and NSGA-II methods perform better than EC. In addition, the -dominance value of WS and NSGA-II methods are close to one, meaning the Pareto solution set provided by these two methods is close to the reference set. Still, our heuristic performs one order of magnitude better than the WS method. In terms of the computation time, NSGA-II is the fastest, and EC takes the longest time.

The above comparisons demonstrate that the NSGA-II method has the best performance regarding the four performance indicators for medium-sized instances.

Results for large-sized instances

Table 5 presents the computational results for the large-sized instances with up to 200,000 recipients and 50 vaccination sites. As the number of recipients and vaccination sites increases, we can see that the performance of EC degrades drastically. For example, it provides on average 13 non-dominated Pareto solutions for instances 1–40 of Tables 3 and 4 that have fewer recipients and vaccination sites, while it only finds at most one and often none for instances 41–64 of Table 5.

The same conclusion can be drawn by observing other performance indicators. Though performing slightly better than the EC, the WS method has difficulties in finding non-dominated Pareto solutions for instances with more than 100,000 recipients. NSGA-II shows its high stability in consistently finding a large number of non-dominated Pareto solutions, a large value of hypervolume ratio, and a small value of average -dominance. NSGA-II can find a set of non-dominated Pareto solutions for each instance. It significantly outperforms the other two methods since the values of hypervolume ratio and -dominance for these instances are very close or equal to one. This is because the Pareto solutions in the reference set are mostly or all contributed by the NSGA-II. The results in Table 5 clearly show the superiority of NSGA-II.

Further observing the computational results of Table 5, WS and EC find on average three and one non-dominated Pareto solutions, respectively, which are significantly fewer than the 99 non-dominated solutions provided by NSGA-II. It can be observed that EC cannot find any non-dominated Pareto solution for 11 out of the 24 instances, and the WS method loses its effectiveness when the number of recipients reaches 150,000 with 50 vaccination sites. The values of the hypervolume ratio obtained by WS and EC are 0.332 and 0.000, respectively. This value reaches 0.998 for NSGA-II, indicating a better coverage in the objective space of NSGA-II. The EC yields a much larger value of -dominance than the two other methods, which indicates that the Pareto solutions of EC are far from the reference set. The -dominance value of NSGA-II is much closer to 1.000 than that of WS, meaning the Pareto solution set provided by NSGA-II is closer to the reference set. The computation time of the EC method is sometimes equal to the time limit of 7200 seconds because it fails to find a feasible solution when solving S-MVPP (). For the WS and EC methods, the total computation time of a run may exceed the time limit of 7200 seconds since we check the stopping criterion after each solution of the single-objective problem.

Table 5 shows that WS and EC methods lose their efficiency and efficacy in solving such large-sized instances, and NSGA-II again performs the best by providing more non-dominated solutions with very good quality.

Performance analysis

In this section, we study the performance of some critical elements of our algorithm that support their choices and pertinence. In Section 5.3.1, we evaluate the choice of the primary objective function of the EC method. In Section 5.3.2, we assess the value of our improved assignment strategy in the design of the initial population. Finally, in Section 5.3.3, we study the efficiency of the dynamic programming which is key to obtaining a solution at each iteration of our heuristic.

Performance analysis of the -constraint method

The selection of the objective to be optimized plays a critical role in the EC method. We denote the version where the first objective (second objective) is selected as the principal objective as -constraint- (-constraint-). We then solve the same 40 instances with the two methods. The results are shown in Table 6 . These indicators are calculated based on the reference set as introduced in Section 5.2.1.

Table 6

Comparison results of different objective functions used as -constraint.

#	|J|	ϵ-constraint-f_1				ϵ-constraint-f_2
		#Nd	Hv	e	Time (s)	#Nd	Hv	e	Time (s)
1–4	200	19.00	0.990	1.006	735	20.00	0.986	1.010	260
5–8	400	19.00	0.971	1.025	2410	19.25	0.987	1.015	1229
9–12	600	13.00	0.660	1.109	3094	18.75	0.990	1.008	1933
13–16	800	10.75	0.714	1.183	3287	19.25	0.987	1.013	2500
17–20	1000	7.75	0.454	1.245	3773	17.00	0.949	1.057	3644
21–24	2000	5.25	0.332	1.315	3886	12.50	0.720	1.190	3975
25–28	4000	6.50	0.347	1.213	3741	11.25	0.716	1.140	3865
29–32	6000	4.50	0.307	1.258	3669	11.50	0.797	1.121	3986
33–36	8000	4.75	0.155	1.219	3669	10.00	0.662	1.180	3896
36–40	10000	3.25	0.286	1.123	3704	10.50	0.759	1.060	4047
Average		9.38	0.522	1.170	3197	15.00	0.855	1.079	2933

It can be seen that -constraint- yields a much larger number of non-dominated solutions, larger hypervolume ratio, and smaller average -dominance, which indicates its better performance. Moreover, the -constraint- is also faster than -constraint-. These results strongly support our choice in the design of this algorithm.

Performance analysis of the improved assignment strategy

We further analyze the excellent performance of the proposed NSGA-II method. As stated, our NSGA-II method differs from existing ones by employing an improved recipient assignment strategy to vaccination sites and an effective dynamic programming method to obtain the optimal inventory and replenishment quantity at each vaccination site. We first highlight the improved recipient assignment strategy by comparing it with a basic strategy where each recipient is assigned to its nearest vaccination site. Four instances are selected and the detailed comparison results are shown in Fig. 7 . We can clearly see that the obtained Pareto front with the improved assignment strategy dominates that of the basic one.

Fig. 7

The Effect of different strategies for assigning the recipients.

Effect of the dynamic programming method

Next, we show the effect of the new dynamic programming (DP) method used to establish optimal replenishment and inventory plans given a recipient assignment. The same four instances are selected. We compare the NSGA-II with our DP strategy against a random strategy (shown in Algorithm 7) where the replenishment quantity is randomly generated in the feasible interval. The computational results are shown in Fig. 8 . The results show that our DP scheme significantly improves the solution quality.

Algorithm 7

Determine the replenishment quantity without the dynamic programming method.

Fig. 8

Effect of a random replenishment strategy instead of our DP.

In summary, the excellent performance of the tailored NSGA-II method is due to the developed new recipient assignment strategy and the DP.

Conclusion

We have investigated the integrated bi-objective multi-period capacitated vaccination planning problem with replenishment. The aim is to obtain a cost-effective vaccination plan and provide a high-quality service level. In particular, the studied problem simultaneously optimizes the total operational cost and travel distance. The first objective represents an economic criterion, and the second one can be viewed as a service quality measure. The studied problem is a generalization of the multi-period facility location problem. It further considers replenishment and capacity decisions at each vaccination site. We have proposed a bi-objective MILP model for the studied new problem and three solution methods (i.e., weighted-sum, -constraint, and NSGA-II) to obtain approximate Pareto solutions. The weighted-sum and -constraint methods are developed based on the proposed MILP. The tailored NSGA-II method uses a problem structure-based initial solution generation scheme to provide approximate Pareto solutions.

To evaluate the performance of the proposed model and algorithms, we have presented a real-world case study and numerical experiments on 64 random instances. Results on the real-world case study indicate that: 1) the proposed model and algorithms can provide a set of Pareto solutions to the studied bi-objective MVPP; 2) the Pareto solutions obtained by the developed methods dominate those obtained by experience-based greedy methods and a real-world solution; and 3) the proposed decision framework provides decision-makers with the flexibility to select preferred solutions from a set of Pareto solutions based on their preferences. In particular, our solutions decrease the operational costs by 9.3% on average and decrease the total distance by 36.6%. Results on randomly generated instances with up to 200,000 recipients, 50 vaccination sites, and 10 days show that the proposed model and algorithms effectively solve the studied bi-objective MVPP. The results consistently indicate that the proposed NSGA-II method significantly outperforms the weighted-sum and -constraint. In particular, the NSGA-II method can provide a large number of Pareto solutions for all large-sized instances.

The mass vaccination planning problem is a challenge for both researchers and practitioners. In the current problem setting, vaccine replenishment is assumed to be shipped directly. The replenishment quantity of a vaccination site may be small, which motivates us to include vehicle routing in the studied problem to further reduce operational costs. In this case, the problem structure changes. We must propose new models for the problem and develop efficient algorithms to solve practical-sized instances. Furthermore, despite the appointment information, the demand may vary due to no-shows or walk-in recipients. Therefore, future studies may focus on a stochastic version considering demand uncertainty. Besides, the equity issue that considers the vaccination distribution or the balance between efficiency and equity deserves further study.

CRediT authorship contribution statement

Lianhua Tang: Writing – original draft, Methodology, Software, Validation, Investigation, Data curation. Yantong Li: Conceptualization, Formal analysis, Methodology, Writing – review & editing, Visualization, Supervision. Danyu Bai: Funding acquisition, Project administration. Tao Liu: Resources, Project administration, Writing – review & editing. Leandro C. Coelho: Writing – review & editing, Visualization.

Sets:
J	set of vaccination recipients, indexed by j;
K	set of vaccination sites, indexed by k;
T	planning horizon, indexed by t.
Parameters:
c_jk	distance between recipient j and vaccination site k;
α_jt	equal to 1 if recipient j needs to be vaccinated on day t;
e_k	fixed cost for replenishing vaccination site k;
f_k	opening cost for vaccination site k;
g_k	unit cost for launching a vaccination station at vaccination site k per day;
h_k	unit inventory holding cost at vaccination site k per day;
m_k	maximum number of vaccination stations at vaccination site k each day;
O_k	maximum replenishment quantity at vaccination site k;
Q	maximum number of recipients that can be served by a station per day;
V_k	maximum inventory capacity at vaccination site k.
Decision variables:
x_jk	equal to 1 if recipient j is served at vaccination site k, and 0 otherwise;
v_kt	equal to 1 if vaccination site k is opened on day t;
w_kt	equal to 1 if vaccination site k is replenished on day t;
z_kt	replenishment quantity of vaccine at site k on day t;
u_kt	number of launched vaccination stations at vaccination site k on day t;
I_kt	inventory quantity at vaccination site k at the end of day t.