key: cord-0722570-bs9tsuhb
authors: Tang, Lianhua; Li, Yantong; Bai, Danyu; Liu, Tao; Coelho, Leandro C.
title: Bi-objective Optimization for a Multi-period COVID-19 Vaccination Planning Problem
date: 2022-02-16
journal: Omega
DOI: 10.1016/j.omega.2022.102617
sha: fae7add1b96b647b5346fddc82eb30f442afccd8
doc_id: 722570
cord_uid: bs9tsuhb

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 [Formula: see text]-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 enlarge the service capacity of vaccination sites to 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.

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 [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.

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.

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. [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. [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. [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, [39] formulate a mixed-integer non-linear programming model and propose a constructive heuristic to determine the best combination of vaccines and their prices. [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. [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. [42] review the vaccine distribution chains in low-and middle-income countries. [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, [44] propose a stochastic inventory model to capture production interruption and evaluate the impact of pediatric vaccine inventory level on vaccination coverage. [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. [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. [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. [48] build a two-stage analytical framework to analyze vaccination hesitancy, showing that vaccination decisions are affected by reference point formation and updating. [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. [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.

[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.

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 singleperiod 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 timevarying, 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 multiperiod 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.

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:

1. 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.

2. 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.

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

This section first elaborates on the studied MVPP with replenishment and then proposes a MILP model. Consider a set K = {1, ..., |K|} of vaccination sites, each site k ∈ K located at (X k , Y k ). A set J = {1, ..., |J|} of recipients are to be vaccinated within the planning horizon T = {1, ..., |T |}. Each recipient j ∈ J has a specific location (X j , Y j ) and an appointment date α jt ∈ {0, 1} indicating the day the recipient expects to get vaccinated. The travel distance between recipient j and vaccination site k is c jk . Each vaccination site k opened on day t has a fixed setup cost f k and a variable cost g k depending on the number of vaccination stations deployed. A fixed replenishment cost e k is incurred whenever vaccination site k is replenished. The unit inventory cost at vaccination site k is h k . Each vaccination site has a maximum replenishment quantity O k , inventory capacity V k , and a maximum number m k of stations that can be deployed. The maximum number of recipients that a vaccination station can serve each day (service capacity) is denoted by Q. 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: x jk is equal to 1 if recipient j is served at vaccination site k; v kt is equal to 1 if site k is opened on day t; w kt is equal to 1 if site k is replenished on day t. We also use variables z kt to represent the replenishment quantity of vaccine at k on day t; u kt the number of vaccination stations at site k on day t; I kt the inventory quantity at site k at the end of day t. Based on the above description, we define the following notation. 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; We next present a bi-objective MILP model (P) 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] .

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 nonconvex 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 multiobjective 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 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 w 1 and w 2 which represent the preference of decision-makers on each objective, where w 1 + w 2 = 1, and w 1 ≥ 0, w 2 ≥ 0. Since the two objectives f 1 and f 2 have different scales, we normalize them to eliminate the effects of dimensions [83] . The normalized objectives f i (i = 1, 2) 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 f I 1 , f I 2 , f N 1 , and f N 2 are obtained by exactly solving the following mono-objective problems:

where x is the vector of variables and χ is the solution space of x.

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

Algorithm 1 Weighted-sum method 1: Solve (14)- (15) and (16)- (17) to obtain the values of

Solve S-MVPP(WS), and obtain f 1 and f 2 . 5:

Update set P = P ∪ {(f 1 , f 2 )}. 6:

Set w1 = w1 + ω. 7: end while 8: Remove the dominated points from P and return P .

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 multiobjective 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: (12), and to

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 f N 1 and a lower limit f I 1 . Then, the S-MVPP( ) is solved iteratively by defining a step size ∆ =

where h is the expected number of Pareto solutions. Initially, takes the value of f N 1 − ∆. Then, it is reduced by the step size ∆ at each iteration. In the s th iteration, where s = {1, 2, , h}, let the value of the first objective be represented by f s 1 . Then is set to be f s 1 −∆ for the (s+1) th iteration. We solve a series of S-MVPP( ) as long as the value of > f I 1 . 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. 

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 Figure 1 , and we describe each step in detail next.

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 v kt which takes a value of 1 if the corresponding vaccination site is opened. The second part corresponds to the recipient allocation decisions, denoted as Y j , indicating the assignment of recipients to a site.

The third part indicates the replenish quantity z kt . The number of genes is 2|K||T | + |J|. Figure 2 shows a chromosome with |J|=6 recipients, |K|=3 vaccination sites, and |T |=2 days. Vaccination Step 2: Population initialization. A set C of P s 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 Assign recipients (Yj (x jk )) to opened vaccination sites. Algorithm 5 5:

Determine the replenishment date (w kt ), quantity (z kt ), and inventory (I kt ). Algorithm 6 6:

Add the individual to C. Set i = i + 1. 2) Allocate recipients to opened vaccination sites (Y j ). Let L kt be the number of recipients assigned to site k at day t. In Algorithm 5, for each recipient that is not allocated, we calculate the distance difference D j between its nearest and second nearest opened vaccination sites with remaining capacity. Next, for each day t and each opened vaccination site k, we define a set U 1 kt with all the unallocated recipients whose nearest site is k and the expected vaccination date is t. We also define another set U 2 kt that contains all the unallocated recipients whose nearest and second nearest sites are l and k, where l < k, and whose expected vaccination date is t. Let

where all recipients in U kt are sorted in non-increasing order of D j . Then we compute a k = min(|U kt |, R k − L kt ), and assign the first a k recipients to k. The values Y j of the first a k recipients in the set U kt is set to k, and L kt = L kt + a k . 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.

3) Determine the replenishment variables z kt , w kt , and inventory quantity variables I kt . Once th site opening variables v kt and recipient assignment variables x jk are determined, the remaining problem is to compute the optimal replenishment and inventory at each opened site k. We recall that L kt is the number of recipients to be served at the vaccination site k on the day t, e k is the 2: if i = 1, then 3:

Generate an extreme solution, i.e., v kt =1 for all k and all t. 4: else 5:

if i = 2, then 6:

Generate another extreme solution, i.e., 7:

for t = 1, 2..., |T |, do 8:

for k = 1, 2..., |K|, do 9:

Randomly generate v kt = 0 or 1. 10:

if fixed replenishment cost at site k, and h k is the unit inventory cost at site k. Then, for each site k we compute the optimal value of the replenishment date variables w kt , replenishment quantity variables z kt , and inventory variables I kt by solving the following MILP model:

s.t. (7), (8) , and to

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

The first term e k w k,t 1 denotes the replenishment cost at period t 1 . The second term indicates the inventory costs from periods s ∈ {t 1 , · · · , t 2 − 1}. The inventory quantity at period τ is equal to the sum of vaccines needed during periods s + 1 to t 2 .

Let f t be the minimum cost from the first day to day t when the inventory is zero at the end Algorithm 5 Assign recipients with an improved strategy 1: Input: αjt, c jk , R k , and initialize Yj = 0 (∀j), L kt = 0 (∀k, t). if Yj = 0 then 6:

Denote K the set of available vaccination sites that opened on the day t|αjt=1. 7:

Sort the sites (k ∈ K) in non-increasing order of c jk , and denote k1 and k2 the two nearest sites. 8:

Calculate Dj = c j,k 2 − c j,k 1 ; if there is only one vaccination site k1 in set K, denote Dj = c j,k 1 . 9:

end if 10:

end for 11:

for t = 1, 2..., |T |, do 12:

for k = 1, 2..., |K|, do 13:

if v kt = 1 and R k > L kt , then 14:

Obtain recipients set U 1 kt whose k1 = k, αjt = 1 and Yj = 0. 15:

Obtain recipients set U 2 kt whose k1 < k, k2 = k, αjt = 1 and Yj = 0. 16:

Denote U kt = U 1 kt ∪ U 2 kt , and set a number a k = min(|U kt |, R k − L kt ). 17:

Assign the first a k recipients with largest Dj to site k, and denote it as B kt . 18:

Denote Yj = k (j ∈ B kt ), and obtain L kt = L kt + a k .

end if 20:

end for 21:

end for 22: end while 23: Return Yj and L kt . of day t. We can use the following state transition function:

where f 0 =0. The objective value can be obtained by calculating f |T | .

The detailed implementation process of the proposed DP is shown in Algorithm 6. This procedure provides the values for variables w kt , z kt , and I kt . Finally, variables u kt can be computed as

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., v kt ) as the basic chromosome because once v kt changes, Y j and z kt 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 v kt . 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 P s 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.

Determine w kt , z kt , and I kt with a DP method 1: Input: L kt , O k , V k , e k , h k , and initialize w kt =0, z kt =0, and I kt =0 (∀k, t).

2: for k = 1, 2..., |K|, do 3:

(1) Calculate the replenishment and inventory costs c(t1, t2) from periods t1 to t2.

for t1 = 1, 2..., |T |, do 5:

for t2 = t1, ..., |T |, do

L kτ , and c(t1, t2) = 0 Denote η(k, t2) = min{δ(:, t2)}, and the index as ζ(t2).

end for 27:

{/ * δ(t1, t2) denotes the costs until t2 when the last replenishment date is t1; ζ(t2) represents the latest optimal replenishment date before time t2; η(k, t2) is the minimum cost for site k until time t2. * /} 28:

(3) Determine the replenishment dates, and replenishment and inventory quantities.

for t = 2, 3..., |T |, do 30: Γ(t, 2)=Γ(t − 1, 1) − 1, and Γ(t, 1)=ζ(Γ(t, 2)); where Γ(1, 1)=ζ(|T |) 31:

z k,Γ(t,1) =α(Γ(t, 1),Γ(t, 2)), w k,Γ(t,1) =1; compute I kt using constraints (21) 32:

if Γ(t, 1) = 1, then 33:

break {/ * All the replenishment dates have been found * /} 34: end if 35: end for 36:

{/ * Γ(t, 1) represents the replenishment date, when the replenishment quantity can meet the vaccines needed from period Γ(t, 1) to Γ(t, 2). * /} 37: end for 38: Return w kt , z kt , and I kt .

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 h 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 P s , the number of generations M , 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.

This section presents a case study from the Tongzhou district of Beijing City, China, to validate Table 1 , along with their locations, visualized in Figure 3 . 

We obtained the vaccination data for 12 days from Mar 19, 2021, to Mar 30, 2021. The total number of vaccinated recipients is 202,596. 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 T P , and the number of permanent residents in each community be P i , where N i=1 P i = T P and N is the number of communities. We then calculate the vaccination rate of the Tongzhou district, which is denoted by r = |J|/T P . Next, we estimate the vaccinated population in each community by V i = rP i . 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. |T | , where ρ is randomly generated from 1, 2, 3. A summary of these values is presented in Table 2 .

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 Figure 4 . 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 w 1 and w 2 for the two objectives, each representing a specific preference on the two objectives. In scenario 1, the weight coefficients w 1 and w 2 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 (w 1 = 0.5, w 2 = 0.5), indicating that the decision-maker believes operational cost and service quality are equally important. In Scenario 3, w 1 and w 2 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 Figure 4 , from which we can see that the method selects three distinct solutions for the three scenarios. We then compare the obtained Pareto front P to the solutions obtained by two sequential heuristics adopted in practice. The heuristics used in practice consist of scenarios where decisionmakers 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 Figure 4 . We see from Figure 4 that for both heuristic solutions, we can find better solutions in the Pareto front P 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 w 1 = 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 w 1 = 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 Figure 5 The analysis of these detailed solutions allows us to derive the following insights. 1) As the weight w 1 increases from 0.2 to 0.8, the operational cost decreases, and the travel distance increases. This is expected since a larger weight w 1 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. 3) In Figure 5(d) , we see that as the weight w 1 of the operational cost objective increases from 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.

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 m k for the above Pareto front P 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 Figure 6 . From the results shown in Figure 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. coordinates (X k , Y k ) of the vaccination sites k ∈ K and (X j , Y j ) of the recipients j ∈ J are randomly generated from [1, 200] . The distance between vaccination site k and recipient j is calculated as

The expected service day for recipient j ∈ J is randomly generated from the interval [1, |T |]. Finally, we generate an instance for each combination of |J|, |K| and |T |, 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.

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 P * 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 e-dominance, and computation time, widely adopted in evaluating multiobjective 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 P excluding those dominated by the solutions in the referenced set P * . 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 H(P * ) be the hypervolume of the reference set, which is calculated based on the nadir point (f N 1 , f N 2 ). It represents the objective space covered by the set P * . The hypervolume ratio of set P can be calculated as H v (P ) = H(P ) H(P * ) . A large hypervolume ratio indicates the solution set has better coverage, which suggests that the corresponding algorithm is better.

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

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

The average e-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 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. 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 nondominated Pareto solution for 11 out of the 24 instances, and the WS method loses its effectiveness Table 5 shows that WS and EC methods lose their efficiency and efficacy in solving such largesized instances, and NSGA-II again performs the best by providing more non-dominated solutions with very good quality.

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.

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-f 1 ( -constraint-f 2 ). 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 P * as introduced in Section 5.2.1. It can be seen that -constraint-f 2 yields a much larger number of non-dominated solutions, larger hypervolume ratio, and smaller average e-dominance, which indicates its better performance.

Moreover, the -constraint-f 2 is also faster than -constraint-f 1 . These results strongly support our choice in the design of this algorithm.

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 Figure 7 . We can clearly see that the obtained Pareto front with the improved assignment strategy dominates that of the basic one.

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 Figure 8 . The results show that our DP scheme significantly improves the solution quality. In summary, the excellent performance of the tailored NSGA-II method is due to the developed new recipient assignment strategy and the DP.

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 realworld case study and numerical experiments on 64 random instances. Results on the real-world case 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.

This work was partly supported by grant 61873173 from the National Natural Science Foundation of China (NSFC), and grant 2019-00094 from the Canadian Natural Sciences and Engineering

Research Council (NSERC). This support is greatly acknowledged. We thank the editor and two anonymous referees for their valuable comments that helped improve this paper.

Appendix: Basic strategy of NSGA-II

The random strategy to establish replenishment and inventory plans is shown in Algorithm 7.

Algorithm 7 Determine the replenishment quantity without the dynamic programming method 1: for k = 1, 2..., |K|, do 2:

for t 1 = 1, 2..., |T |, do 3:

Randomly generate z kt from [A, B], where A = max{0, L kt − I k,t−1 + max{0, L k,t+1 − O k }, which represents the replenishment quantity of vaccines should not less than the demands. B = min{ |T | τ =t L kτ −I k,t−1 , (V kt +L kt −I k,t−1 ), O k }, which represent the limits of maximum replenishment quantity and inventory quantity.

Calculate the inventory quantity I kt by constraints (6).

Verity whether z kt and I kt exceed the limits of replenishment and inventory quantity. If yes, go back to step 3; otherwise, get the values of z kt and I kt . 

Comparison with past pandemics

Modeling vaccine allocations in the COVID-19 pandemic: A case study in Australia

Data analytics for operational risk management

Predicting the impacts of epidemic outbreaks on global supply chains: A simulation-based analysis on the coronavirus outbreak (COVID-19/SARS-CoV-2) case

Impacts of epidemic outbreaks on supply chains: mapping a research agenda amid the COVID-19 pandemic through a structured literature review

A review of the existing and emerging topics in the supply chain risk management literature

Viable supply chain model: integrating agility, resilience and sustainability perspectiveslessons from and thinking beyond the COVID-19 pandemic

Viability of intertwined supply networks: extending the supply chain resilience angles towards survivability. A position paper motivated by COVID-19 outbreak

Innovative bring-service-near-your-home operations under corona-virus (COVID-19/sars-CoV-2) outbreak: Can logistics become the messiah?

OR/MS models for supply chain disruptions: a review

WHO, Coronavirus disease COVID-19 advice for the public

Heightened levels of maladaptive daydreaming are associated with COVID-19 lockdown, pre-existing psychiatric diagnoses, and intensified psychological dysfunctions: A multi-country study

10 Italian towns in lockdown over coronavirus fears

Curfew as US cities shut down in coronavirus fight

Australia's third-largest city of Brisbane to enter Covid lockdown

A decision support system for demand management in healthcare supply chains considering the epidemic outbreaks: A case study of coronavirus disease 2019 (COVID-19)

WHO, Tracking SARS-CoV-2 variants

Preparedness of countries to face COVID-19 pandemic crisis: Strategic positioning and underlying structural factors to support strategies of prevention of pandemic threats

Center for Disease Control and Prevention, COVID19 vaccination program interim playbook for jurisdiction operations

Optimal levels of vaccination to reduce COVID-19 infected individuals and deaths: A global analysis

SARS-CoV-2 vaccines in development

Towards effective COVID19 vaccines: Updates, perspectives and challenges (review)

Optimal planning of the COVID-19 vaccine supply chain

Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period

Attaining herd immunity to a new infectious disease through multi-stage policies incentivising voluntary vaccination

Herd immunity: Understanding COVID19

Transforming COVID-19 vaccines into vaccination: Challenges and opportunities for management scientists

Locate vaccination stations considering travel distance, operational cost, and work schedule

On the fair optimization of cost and customer service level in a supply chain under disruption risks

Digital technology and COVID-19

A decision support system for demand management in healthcare supply chains considering the epidemic outbreaks: A case study of coronavirus disease 2019 (COVID-19)

OR-methods for coping with the ripple effect in supply chains during COVID-19 pandemic: Managerial insights and research implications

Stochastic optimization of supply chain resilience under ripple effect: A COVID-19 pandemic related study

On the risk-averse selection of resilient multi-tier supply portfolio

Reconfigurable supply chain: the x-network

An inventory-location optimization model for equitable influenza vaccine distribution in developing countries during the COVID19 pandemic

Where to locate COVID19 mass vaccination facilities?

Making combination vaccines more accessible to low-income countries: The antigen bundle pricing problem

Dose-optimal vaccine allocation over multiple populations

Optimal influenza vaccine distribution with equity

Vaccine distribution chains in low-and middle-income countries: A literature review

Optimizing vaccine distribution networks in low and middle-income countries

An analysis of the pediatric vaccine supply shortage problem

Reservation and allocation policies for influenza vaccines

Vaccine cold chain management and cold storage technology to address the challenges of vaccination programs

Optimal two-phase vaccine allocation to geographically different regions under uncertainty

Can reference points explain vaccine hesitancy? A new perspective on their formation and updating

Adolescents attitudes to the COVID-19 vaccination

Modelling of COVID19 vaccination strategies and herd immunity, in scenarios of limited and full vaccine supply

Pricing the COVID-19 vaccine: A mathematical approach

On solving complex multi-period location models using simulated annealing

Ambulance location and relocation models

Modeling the budget-constrained dynamic uncapacitated facility locationnetwork design problem and solving it via two efficient heuristics: A case study of health care

Multi-Period Facility Location, Location Science

Multi-dimensional location problems

A comparative study of approaches to dynamic location problems

A dual-based procedure for dynamic facility location

Multi-period capacitated facility location under delayed demand satisfaction

Multiperiod capacitated location models, Discrete Location Theory

Facility location dynamics: An overview of classifications and applications

A stochastic multi-period capacitated multiple allocation hub location problem: Formulation and inequalities

Optimizing dynamic facility location-allocation for agricultural machinery maintenance using Benders decomposition

Robust facility location under demand uncertainty and facility disruptions

A multi-period facility location problem with modular capacity adjustments and flexible demand fulfillment

Dynamic facility location with generalized modular capacities

Lagrangian heuristics for large-scale dynamic facility location with generalized modular capacities

Heuristics for the dynamic facility location problem with modular capacities

Adaptive weighted-sum method for bi-objective optimization: Pareto front generation, Structural and Multidisciplinary Optimization

The weighted sum method for multi-objective optimization: new insights

A robust counterpart approach to the bi-objective emergency medical service design problem

A multiobjective, maximal conditional covering location problem applied to the relocation of hierarchical emergency response facilities

Evaluation of multi-objective optimization approaches for solving green supply chain design problems

Multiobjective decision making: theory and methodology

Robust, multi-objective optimization for the military medical evacuation location-allocation problem

Overview of NSGA-II for optimizing machining process parameters

A fast and elitist multiobjective genetic algorithm: NSGA-II

Multi-objective optimization of surface grinding process using NSGA II

Multi-objective optimization of cutting parameters with improved NSGA-II

Solving multi-objective parallel machine scheduling problem by a modified NSGA-II

A NSGA-II based memetic algorithm for multiobjective parallel flowshop scheduling problem

Microscopic optimization model and algorithm for integrating train timetabling and track maintenance task scheduling

Multiobjective optimisation of production, distribution and capacity planning of global supply chains in the process industry

Corrigendum to multiobjective optimisation of production, distribution and capacity planning of global supply chains in the process industry

An exact -constraint method for bi-objective combinatorial optimization problems: application to the traveling salesman problem with profits

Effective implementation of the -constraint method in multi-objective mathematical programming problems

Dynamic version of the economic lot size model

Integrated production inventory routing planning for intelligent food logistics systems

Service-oriented bi-objective robust collectiondisassembly problem with equipment selection

A critical survey of performance indices for multi-objective optimisation

Conceptualization, Model analysis

Danyu Bai: Funding acquisition

Writing-Review & Editing, Visualization