Dynamic Planning of a Two-Dose Vaccination Campaign with Uncertain Supplies.

The ongoing COVID-19 pandemic has led public health authorities to face the unprecedented challenge of planning a global vaccination campaign, which for most protocols entails the administration of two doses, separated by a bounded but flexible time interval. The partial immunity already offered by the first dose and the high levels of uncertainty in the vaccine supplies have been characteristic of most of the vaccines used worldwide and made the planning of such a vaccination campaign extremely complex. Motivated by this compelling challenge, we propose a stochastic optimization framework for optimally scheduling a two-dose vaccination campaign in the presence of uncertain supplies, taking into account constraints on the time that should elapse between the two doses and on the capacity of the healthcare system. The proposed framework seeks to maximize the vaccination coverage, considering the different levels of immunization obtained with partial (one dose only) and complete vaccination (two doses). We cast the optimization problem as a convex second-order cone program, which can be efficiently solved through numerical techniques. We demonstrate the potential of our framework on a case study calibrated on the COVID-19 vaccination campaign in Italy. The proposed method shows good performance when unrolled in a sliding-horizon fashion, thereby offering a powerful tool to help public health authorities calibrate the vaccination campaign, pursuing a trade off between efficacy and the risk due to shortages in supply.

Introduction

In response to the COVID-19 health crisis, we witnessed an unprecedented mobilization of the scientific community. Within this joint effort, mathematical models and operational research (OR) have emerged as valuable tools to help assist public health authorities in their complex decisions during epidemic outbreaks (Nowzari, Preciado, Pappas, 2016, Zino, Cao, 2021). In particular, during the first stages of the COVID-19 pandemic, they have been supportive in forecasting the epidemic spread (Calafiore, Novara, Possieri, 2020, Parino, Zino, Porfiri, Rizzo, 2021b, Taylor, Taylor, 2021, Truszkowska, et al., 2021), including local outbreaks (Chang & Kaplan, 2021), understanding how to allocate patients to facilities (Hosseini-Motlagh, Samani, & Homaei, 2021), and helping to devise nonpharmaceutical interventions (Eryarsoy, Shahmanzari, Tanrisever, 2022, Fischetti, Fischetti, Stoustrup, 2021, Tsiligianni, Tsiligiannis, Tsiliyannis, 2022). With the availability of vaccines, new challenges have emerged concerning the planning of an unprecedented mass vaccination campaign. Also in this phase, mathematical models and OR have been key to help assist public health authorities (Choi, 2021).

Most of the OR literature regarding vaccination has focused on optimizing the vaccine supply chain (Duijzer, van Jaarsveld, & Dekker, 2018). For instance, in Lin, Zhao, & Lev (2022), the authors focus on matching uncertainties on the supply-side and on the demand-side; in Arifoğlu & Tang (2022), an incentive program for coordinating the supply chain is proposed for the seasonal Flu vaccine; and in De Boeck, Decouttere, Jónasson, & Vandaele (2021), supply chains in resource-limited settings are optimized using a spatial model. Also, from an epidemic modeling perspective, the design of optimal vaccination strategies is a largely debated topic and a whole body of literature has been developed toward encapsulating vaccines and their effect in epidemic models. In these models, vaccinations are incorporated by means of additional compartments in population or network models, or via game-theoretic mechanisms to capture individual-level decisions on whether to take the vaccine shot or not (Wang et al., 2016).

Despite these efforts, the 2021 COVID-19 vaccination campaign has laid bare unprecedented difficulties in planning a massive vaccine rollout that entails the administration of two vaccine doses in a short time, facing limited stocks and uncertain supplies, especially in the early phases of the campaign. Many of the works modeling the COVID-19 vaccination campaign mimic the effect of immunization by encapsulating a simplistic implementation of the vaccination process into the disease progression model. In Giordano (2021); Sinha, Kumar, & Chandra (2021), one-dose vaccination processes are used to investigate the prioritization of the recipients of the vaccine. Leveraging a coordinate descent algorithm, Bertsimas, Digalakis, Jacquillat, Li, & Previero (2022) study how to optimally locate mass vaccination centers. The problem of optimally scheduling vaccination appointments given a set of spatially-located vaccination sites is treated in Zhang, Li, Cao, & Wen (2022) by formulating a mixed-integer linear program and approximating the solution via heuristics.

Another critical issue concerns devising inventory policies to optimize the completion of the recommended two-dose administration. In Parino, Zino, Calafiore, & Rizzo (2021a), a model predictive control approach is used to optimally design two-dose vaccination rollouts in a deterministic framework. A similar problem is considered in Mak, Dai, & Tang (2021) using a deterministic compartmental model. Other studies assess the effectiveness of delaying the second vaccine dose (Moghadas, et al., 2021, Tuite, Fisman, Zhu, Salomon, 2021). While these works provide interesting insight, they rely on the limiting assumption that vaccine supplies are deterministic. To the best of our knowledge, the problem of a two-dose vaccine allocation with a stochastic supply is treated only in Shumsky, Smith, Hoen, & Gilbert (2021). Therein, the problem is formulated as a stochastic dynamic program, for which some insight into the structure of optimal solutions can be gained. However, the size of the state space and the complexity of the stochastic dynamic program make the computation of optimal solutions unfeasible, limiting the findings to heuristics.

To fill in this gap, here we present an effectively-solvable optimization framework for devising a two-dose vaccination campaign in the presence of uncertain supplies. First, we introduce a linear programming (LP) approach for dealing with the ideal, deterministic scenario in which the time profile of the vaccine supplies is assumed to be exactly known in advance. Such a simplified approach allows us to discuss the constraints, including the bounds on the interval between the administration of the first and the second dose of the vaccine, and to formulate an objective function that accounts for the effectiveness of the partial immunization achieved with a single dose of the vaccine (Polack, et al., 2020, Voysey, et al., 2021). Then, we discuss a more realistic scenario in which random uncertainty affects the vaccine supply, toward the definition of a stochastic optimization problem with chance constraints and affine decision variables to address such uncertainty in supply. This problem is then shown to be representable in the form of a convex second-order cone program (SOCP), and hence amenable to an efficient global solution (Lobo, Vandenberghe, Boyd, & Lebret, 1998).

We use the proposed optimization framework to investigate the optimal vaccination policies for the 2021 COVID-19 vaccination campaign. In particular, we calibrate the model on the Italian vaccination campaign, utilizing publicly available data to fit the uncertainty model for the supplies (Vaccini, 2021), and we consider the two most commonly used vaccines, which have different characteristics in terms of effectiveness and administration timing: Comirnaty by Pfizer-BioNTech and Vaxzevria by AstraZeneca (Polack, et al., 2020, Voysey, et al., 2021). Besides demonstrating the use of the proposed optimization framework in a real-world scenario, our findings allow us to contribute to the fast-growing field of vaccine distribution in OR by discussing some important properties of the optimal solution. In particular, we focus on the effectiveness of implementing sliding-horizon schemes, whereby the optimization problem is solved over a shorter time-window than the entire optimization horizon, the optimal solution is implemented, and then re-evaluated at fixed time-steps. Furthermore, we analyze the sensitivity of the solution to the uncertainty in the supply. Finally, we discuss how the proposed optimization framework can be utilized by public health authorities to trade-off efficacy of the vaccination campaign and the risk of setbacks due to shortages in the supply.

Our contribution is threefold. First, we propose a novel optimization framework to address the so far overlooked challenge of optimally devising a two-dose vaccination campaign in the presence of uncertain supplies. The flexibility of our optimization framework allows for its applicability to many real-world scenarios. Second, we provide a mathematical formulation of the optimization problem in the form of a SOCP. Such a formulation enables the use of efficient and computationally-effective solvers and algorithms to derive the optimal solution of the optimization problem. Third, we discuss a relevant application of the proposed framework to the 2021 COVID-19 vaccination campaign in Italy, illustrating how the framework allows for a straightforward calibration of the model parameters to real-world data and discussing some important properties of the optimal solution and of the real-world implementation of our OR method.

Problem setting

Notation

We denote the set of real, nonnegative real, nonnegative integer, and strictly positive integer numbers as , , , and , respectively. A vector is denoted with bold lowercase font, with th entry . A matrix is denoted with bold capital font, with denoting the th entry of the th row. Vector is a all-0 vector, except for a 1 in the th entry. The all-1 and all-0 matrices are denoted by and , respectively, and the identity matrix by . When needed, dimensions are denoted as one positive integer index (for vectors) and two indices (for matrices), which are the number of rows and columns, respectively. The superscript denotes the transpose operator. Inequalities for vectors or matrices are intended entrywise (i.e., means ). We denote by the probability of an event and by the expectation of a random variable (RV). We denote by the cumulative distribution function (CDF) of a Gaussian RV with zero mean and unit variance, and by its inverse (quantile function).

Problem setting

Due to the practical nature of the problem, which considers the scheduling of vaccination for large populations (e.g., entire countries), we approximate the population as a continuum of individuals, which allows us to avoid all the issues and limitations related to the optimization of integer variables. Let

be a discrete-time index;

denote the time-horizon of the vaccination campaign. The time-horizon is measured in time units of constant (fixed) duration. In the rest of this paper, we will consider one week as a time unit;

and denote the minimum and maximum interval between the two doses of the vaccine (here, in weeks), respectively (i.e., the second dose should be administered between the beginning of the th week and the beginning of the th week after the first one); is the duration of the time-window for the second dose;

be the number of individuals in the population who are eligible for receiving the vaccine but have not received the first dose at the beginning of week ;

be a -dimensional vector where the component denotes the number of individuals at the beginning of week who have received the first dose weeks before, but have not received the second dose yet;

denote the number of individuals who, at the beginning of week , have received the first dose weeks before or more, but have not received the second dose yet. For these individuals, the second dose is delayed by more than the maximum amount of time admitted;

be the number of individuals who have already received the second dose of the vaccine at the beginning of week ;

be the number of first doses administered to individuals counted by the variable during week ;

be a -dimensional vector, where the th component , for , denotes the number of second doses administered during week to individuals who have received the first dose weeks before (i.e., those counted by ); the entry denotes the number of second doses administered to individuals who have received the first dose since more than weeks (i.e., those counted by );

be the (cumulative) vaccine supply received up to and including the beginning of week ;

be the total capacity of vaccines that the healthcare system is able to administer during week .

The dynamics of the vaccination procedure is described by the following -dimensional system of linear recurrence equations:which are initialized at the beginning of the vaccination campaign with , where is the total number of individuals in the population that should be vaccinated, for , , and .

The linear dynamics in Eq. (1) can be conveniently cast aswhere the state variable is the -dimensional vector , the decision variables are the scalar and the -dimensional vector , and , , and are defined as

Optimization problem

An LP approach for the deterministic scenario

Before introducing the stochastic framework for the uncertain supply case, we first discuss a “nominal” scenario in which supplies are deterministic.

We assume that the cumulative supply is a deterministic sequence, entirely known at time .

Under this assumption, the problem of designing an optimal two-dose vaccination strategy can be cast as an LP, hence its global optimal solution can be computed with extreme numerical efficiency (Boyd & Vandenberghe, 2004). Our purpose is to devise an optimal vaccination strategy in terms of first and second doses, i.e., designing and , over a given time horizon , so that an appropriate objective function involving the number of vaccinated individuals over the considered time horizon is maximized, while satisfying suitable constraints on the vaccination process. We start by discussing the constraints.

Nonnegativity. The state variable and the optimization variables should be nonnegative (entrywise), i.e., Note that in Eq. (3), nonnegativity is required also at time , to ensure that the control implemented at time yields a feasible state.

Stock. The cumulative number of doses administered until week cannot exceed the cumulative vaccine supply received up to the beginning of that week, which can be conveniently formulated as

Capacity. The total number of doses administered during week , cannot exceed the weekly capacity, i.e.,

Delay. Second doses should be performed in the interval between weeks and weeks after the first dose. We may allow for a fraction of the population who received the first dose to wait past the limit weeks for receiving the second one. If excess wait is not allowed, then is set in the model. This yields the following constraints:

For the sake of readability, we define three vectors , , and that gather the state and decision variables over the entire duration of the vaccination campaign, respectively. We define the following linear objective function to be maximized:The first summand of Eq. (9) carries a positive contribution and accounts for the number of fully vaccinated individuals during the entire duration of the vaccination campaign. The second summand, on the other hand, entails a negative contribution and accounts for the population that is still waiting to start the vaccination process. The parameter weights the relative contribution of the two terms: for , the optimizer maximizes the number of fully vaccinated individuals; for , instead, the key objective tends to be the minimization of the population that has still to receive the first dose. Note that, even though the expression in Eq. (9) does not directly depend on the decision variables, such a dependence is indirectly shaped by the evolution of the state , which is determined by Eq. (2) and thus depends on the decision variables and . To sum up, under Assumption 1, the vaccination problem can be formulated as an LP, in which the linear objective function is maximized subject to the linear constraints in Eqs. (2)–(8).

Eq. (9) can be re-written in order to explicitly consider the number of individuals with one and two vaccine doses. By rewriting the number of unvaccinated individuals as the difference between the total population and the number of individuals vaccinated with at least one dose, we obtainThe first term accounts for the number of fully vaccinated individuals; the second one for the number of individuals with a single dose; while the last one (i.e., ) is a constant that can be obviously dropped in the process of maximization. Hence, if we conveniently set the ratio equal to the ratio between the effectiveness of a single dose of vaccine and the one of two doses, then maximizing the function implies the maximization of the vaccination coverage over time, in terms of the sum of the total number of individuals vaccinated with a single dose multiplied by the effectiveness of a single dose and the total number of individuals with two doses multiplied by the effectiveness of two doses, up to some constants that play no role in the maximization.

Dealing with supply uncertainty

In the approach discussed in the previous section, we computed the entire optimal vaccination schedule over the horizon , based on an assumed perfect knowledge at of the entire supply function over (Assumption 1). However, this assumption may be quite unrealistic in practice, since on the one hand there is uncertainty on the future delivery of vaccines, and on the other hand, there is no necessity of devising a fully “here-and-now” schedule at for the whole foreseeable future, since as time unrolls we may want to exploit the opportunity of suitably adjusting the schedule on the basis of the actual supply. In the remaining of this section, we discuss how to tackle both these issues, and we eventually formalize the problem of designing an optimal vaccination strategy in the scenario of uncertain supplies as a second-order cone programming (SOCP) problem (Lobo et al., 1998).

In the following, we assume that the (cumulative) supply at the generic time-instant is known up to a given level of uncertainty, and we assume that the level of uncertainty (i.e., its standard deviation) may vary with time in a fashion known a priori. This setting is summarized in the following assumption.

We assume that the cumulative supply is in the formwhere is the given nominal (or expected) time profile of the cumulative supply, is the given time profile of the standard deviation of the uncertainty, and , , is a sequence of independent and identically distributed (IID) RVs with zero mean and unit variance, representing the uncertainty on the cumulative supply at time .

The uncertainty model in Assumption 2 is able to capture the presence of delays (or earliness) in the actual supply with respect to the scheduled one and relies on the assumption that suppliers, when not able to comply with an order, will try to accommodate it in the next delivery. Hence, a delay (or earlyness) at time does not affect the expected future supplies.

At this stage, we let the decision variables related to the first and second dose to react to changes in the cumulative supply. Specifically, we assume vaccination policies to take an affine form, i.e.,where and represent the nominal vaccination schedules, and and represent the adjustment gains, whose purpose is to correct the nominal decisions and as a consequence to deviations of the cumulative supply from its nominal (or expected) value. Note that the policies in Eq. (12) are fully causal, since the decision variables at depend only on the uncertain supply realization up to , while the future supplies remain uncertain and are suitably treated via probabilistic constraints, as discussed next. The idea of linear policies, or affine-adjustable variables, has been introduced in the context of robust optimization in Ben-Tal, Goryashko, Guslitzer, & Nemirovski (2004) and has been adopted in several uncertain decision problems, see, e.g., Calafiore (2009).

The decision variables of the optimization problem will now be , and , , , rather than , . The optimization problem will therefore provide us with the parameters to plug into Eq. (12) and, from these equations at each , based on the observation of the actual outcome of , we shall decide for the actual number of first and second dose vaccinations to be delivered during day .

Under Assumption 2, substituting Eqs. (11) and (12) into the dynamics in Eq. (2), we can determine the recursions that govern the evolution of the stochastic system. Specifically, the state variable of the system evolves according to the stochastic recursionSince the state is an RV, it might be convenient and useful for the derivations in the following to explicitly write two different expressions for the temporal evolution of its expected value and for the shift of the RV with respect to it. Specifically, if we define the expected value and the shift , we have thatandwhere the vector is defined asFinally, we derive a recursion equation for the covariance matrix of , defined as , for which it holds

Due to the stochasticity of the system, the constraints should be expressed in terms of probabilities (Calafiore, Ghaoui, 2006, Prékopa, 1995). To this aim, we fix the risk parameter , which represents the probability that a constraint is violated, and we cast the constraints as follows.

Nonnegativity. Since the state and the number of doses performed are stochastic, we re-formulate Eqs. (3)–(5) in a probabilistic fashion as

Stock. Similarly, we impose that

Capacity. We set

Delay. The constraints can be stated on the expected states, i.e.,

Finally, the objective function in Eq. (9) is now an RV, since it depends on the realization of the stochastic process . The approach we take now is to maximize its expected value, yielding the following expectation objectivewhere must satisfy the recursion equation in Eq. (14), and the constraints discussed above. In general, maximizing Eq. (24) might be a difficult task, due to the stochastic nature of the system and, ultimately, the nature of the uncertainties . In the following, we will show that, under some reasonable assumptions on the nature of the uncertainties , the optimization problem can be cast as a SOCP, and thus easily solved.

Convex second-order cone problem formulation

From now on, we make the following assumptions about the uncertainties.

We assume that Assumption 2 holds, and that , for is a sequence of IID Gaussian RVs with zero mean and unit variance.

Under Assumption 3, Eq. (13) yields a sequence of Gaussian random vectors with expected value and covariance matrix  (Ross, 2019). Using this property, we derive the following explicit expressions for the constraints.

Nonnegativity (state variable). From Eq. (18), we observe that the generic th component of the state variable is a Gaussian-distributed RV with expected value and variance . We have therefore that Thus, if and only if (iff) , which in turn holds iff . Hence, being , the generic th constraint in Eq. (18) is equivalent toNotice that we have from Eq. (17) that , and that we can write the following recursion , . Hence, we can write where . It follows that Eq. (25) is expressed as a second-order cone constraint , where, by Eq. (16), is an affine function of the variables , . Hence, the (probabilistic) nonnegativity constraint on the state variables in Eq. (18) can be written as

Nonnegativity (optimization variables). Eq. (19) yields . Hence iffSimilar, from Eq. (20) we derive the explicit constraint

Stock. We rewrite the stock constraints in Eq. (21) aswhere . We observe that is the linear combination of independent Gaussian RVs. Hence, also is a Gaussian RV with and , where . The constraints in Eq. (29) are thus rewritten as , where we set for brevity . The above probability is greater than or equal to iff , whence we obtain the explicit constraint:

Capacity. The capacity constraints in Eq. (22) are expressed, for all , as , which, proceeding as in the previous cases, is equivalent to

Delay. Finally, we note that the expression in Eq. (23) is already formulated in terms of expected states.

If we stack all the , , , and variables in four vectors , , , and , we can formulate the optimization problem as:s.t.:   Eqs. (14), (26), (27), (28), (30), (31), (23).

Problem (32) has variables — state variables and decision variables— and constraints, and involves the optimization of a linear objective function subject to linear constraints and the intersection of second-order cones. Hence, it is a SOCP, for which the global optimal solution can be computed in polynomial time via interior point methods (Lobo et al., 1998), and standard and efficient algorithms have been developed, see, e.g., Mosek (Andersen, Roos, & Terlaky, 2003).

Despite the fact that the optimal policy in Eq. (12) is reactive and takes into account the uncertainty, one should consider that, as time unrolls and vaccines are supplied, part of the uncertainty is revealed. Hence, in practical applications, it makes sense to deploy Problem (32) within a sliding-horizon framework.

Fixing a simulation horizon and a re-evaluation step , every time steps starting from , one could utilize Problem (32) to compute the optimal vaccination policy for the following time steps and apply it for the first steps. Then, after time-steps, one can re-evaluate the policy by utilizing again Problem (32), with updated initial conditions and uncertainty model (informed by the time-series of supplied observed in the real world), and apply it for the following time-steps. This process can be repeated until the end of the time-horizon of the vaccination campaign. Similar sliding-horizon schemes have been successfully used in the literature of optimization under uncertainty, sometimes under the name of folding horizon. See, e.g., Ben-Tal, Golany, Nemirovski, & Vial (2005); Cohen, Golany, & Shtub (2016); Dabadghao, Marandi, & Roy (2021); De Ruiter, Brekelmans, & den Hertog (2016); Gorissen, Yaníkoğlu, & den Hertog (2015).

Particular combinations of the parameters and input values may result in unfeasible optimization problems. Specifically, the nonnegativity of the optimization variables in Eqs. (19) and (20) and the stock constraint in Eq. (21) may not be simultaneously verified for some values of . This can typically happen in the first weeks, in the presence of small supplies with high uncertainty. A simple way to address this issue in practice is to fix and for all and solve the Problem (32) starting from week , where : Problem (32) is feasible for all .

Application to 2021 COVID-19 vaccination campaign in Italy

To demonstrate the proposed optimization framework against real-world data, we use it to investigate the 2021 COVID-19 vaccination campaign in Italy.

Vaccine parameters

We consider Comirnaty and Vaxzevria, which have been the two mostly used vaccines in the first phases of the vaccination campaign.

Comirnaty. In Italy, the two doses are administered at least weeks apart and at most weeks apart. The efficacy during the interval between the two doses was estimated in clinical trials as 52%, while both doses yield a 95% protection (Polack, 2020). Hence, following Remark 1, we set .

Vaxzevria. The Italian Medicines Agency indicates that the administration of the second dose should be ideally performed during the twelfth week and, in any case, at least ten weeks after the first dose. Consequently, we set and . Clinical data confirm a vaccine efficacy of after the first dose and 82% after two doses (Voysey, 2021), yielding .

Expected supply and uncertainty model

We use the official data of the daily supply of COVID-19 vaccines in Italy from January 4, 2021 to July 19, 2021 (Vaccini, 2021). We set the discrete time unit equal to one week and we aggregate the data accordingly.

Based on these data, we construct a model for the expected cumulative supply . The unfolding of the supply delivered (see Fig. 1 A, B) suggests that the size of weekly supplies increases linearly with time. Hence, we suppose a quadratic trend for the expected cumulative supply . We model the data by using a weighted least-square technique (Chatterjee & Simonoff, 2013). In particular, since the th entry of the cumulative supply vector is the sum of vaccine deliveries, each one characterized by its level of uncertainty, in fitting the weighted least-square model, we associate with the th data-point of the Italian cumulative supply the weight . Using the weighted least-square estimator we obtain for Comirnaty and for Vaxzevria, where Comirnaty is administrated starting from the week of January 19, while the administration of Vaxzevria starts from the week of February 22. Also, we use the fitted model to estimate the uncertainty associated with the predictions (i.e., the model prediction errors or forecast error), whose values are reported in Table S1 in the Supplementary Material. We use these estimations as , and we finally obtain the supply in the form discussed in Eq. (11), under Assumption 3, i.e., assuming that , for is a sequence of IID Gaussian RV with zero mean and unit variance. The resulting fitting with the correspondent 95% prediction intervals for the two vaccine supplies are reported in Fig. 1.

Fig. 1

Italian cumulative supply for the two types of vaccines (blue) and the fitted model , with the confidence intervals (orange and gray, respectively).

Finally, we set the maximum weekly capacity , i.e., we assume that the healthcare system can accommodate a number of vaccine doses that is slightly greater than the expected weekly supply.

Analysis of the optimal solutions

We compute the optimal vaccination rollout using the optimization framework described in Section 3.3, with the expected supplies and uncertainty model derived for the 2021 COVID-19 vaccination campaign in Italy. We use the actual time-series of the Italian supply from Vaccini (2021) as the realization of the process. Unless stated otherwise, we set .

First, we compare the performance of the optimal solution with the baseline performance of the solution of the LP proposed in Section 3.1, computed assuming that the exact time-series of the supplies are known a priori. This solution, even though unrealistic since it ignores the uncertainty, serves as a benchmark for evaluating the performance of the other approaches and the effect of uncertainty on the vaccination campaign. With respect to the LP solution, with both vaccines we observe that the presence of uncertainty may slow down a vaccination campaign. In particular, we observe a reduction of the total number of doses administrated for the Comirnaty vaccination campaign. For Vaxzevria, which is characterized by higher levels of uncertainty, the performance are even worse, as we record a drop in the number of administered doses.

We next discuss a sliding-horizon implementation of the proposed approach, following Remark 3. We compare different solutions, obtained using different sliding-horizon policies (all with ): (a) no re-evaluation, which is the numerical solution of SOCP computed over the entire horizon; (b) weekly re-evaluation, in which the sliding horizon moves on a weekly basis (): every week the uncertainty model is re-estimated following the procedure described in Sec. 4.2, appropriately shifting the time index; then, the optimal solution is re-computed by solving the SOCP over a time horizon and it is applied for week, before the following re-evaluation; (c) monthly re-evaluation, in which the uncertainty model is updated and the strategy is re-evaluated on a 4-week basis (); and (d) bimonthly re-evaluation, in which the uncertainty model is updated and the strategy is re-evaluated every weeks. In other words, once an optimal solution is found, it will be applied for the next weeks before being re-evaluated. Together with the re-evaluation of the optimal solution based on the increased data on the time-series of the supplies, a key advantage of the sliding-horizon technique resides in the possibility of utilizing such data to also re-calibrate the uncertainty model at each iteration, thereby reducing the level of uncertainty in the supply as the vaccination campaign proceeds.

In Fig. 2 A–D, for both vaccines we report the optimal solutions for the cases of no re-evaluation, weekly, and monthly re-evaluation. We notice that the solutions for the entire horizon result in a smooth trend of first and second doses, while the optimal solutions with the sliding horizon have more discontinuous trends, as they are able to re-evaluate the campaign in response to changes in the supply. Fig. 2E, F shows the temporal evolution of the number of vaccinated population under the different policies, illustrating how the implementation of a sliding horizon could be beneficial. We analyze in detail the different cases in Table 1 . We compare the strategy that uses the sliding horizon with the case with no re-evaluation. The comparison is presented as percentage variation of i) the: average utility function during the optimization , which is related to the average coverage of the population during the vaccination campaign, as discussed in Remark 1, and ii) the total number of doses of vaccine administered during the entire campaign. The results show the beneficial effect of implementing a sliding-horizon technique and highlight a nontrivial trade-off between performance and re-evaluation periods, whereby more frequent re-evaluations may result in a performance increase.

Fig. 2

Optimal solution with no re-evaluation (blue), weekly (orange), and monthly re-evaluation (green) for the two vaccines. Panels A–B show the time-series of first doses; C–D of second doses; and E–F the evolution of the number of fully vaccinated individuals.

Table 1

Performance with sliding horizon, with respect to the scenario with no re-evaluation.

	Comirnaty		Vaxzevria
	ΔJ	ΔV	ΔJ	ΔV
re-evaluation T_r=1 week	2.95%	6.30%	0.26%	7.74%
re-evaluation T_r=4 weeks	1.71%	6.64%	0.15%	5.54%
re-evaluation T_r=8 weeks	1.52%	6.68%	0.13%	3.52%

In particular, the use of a sliding horizon yields an increase in the total number of doses administered. For Comirnaty —which is characterized by lower levels of uncertainty— such an increase is approximately by and it seems not to be strongly influenced by the re-evaluation frequency. However, we observe that shorter re-evaluation periods can be beneficial toward increasing the vaccination coverage during the campaign, with increases from for bimonthly re-evaluation to for weekly re-evaluation. On the other hand, the higher levels of uncertainty in the Vaxzevria supply seem to further increase the benefits of frequent re-evaluations, as also the total number of doses administered seems to be strongly affected by the re-evaluation policy. In fact, the implementation of a sliding-horizon technique increases the doses administered by in the case of 1 week re-valuation.

Summarizing, being able to set up flexible vaccination rollouts that can be re-evaluated in order to accommodate potential changes in the expected supply by utilizing the proposed optimization framework with a sliding horizon can bring to a considerable improvement of the performance of the vaccination campaign. In addition, the frequency of re-evaluation seems key in determining the increase in the performance, especially when the supply is highly uncertain.

We solved the problem using Mosek (Andersen et al., 2003) implemented with CVX in MatLab R2021a. The most computationally demanding problem (Comirnaty) involves 667 variables and 809 constraints. It was solved in approximately 24 seconds on a 2.9 Gigahertz 6-Core Intel Core i9.

Sensitivity with respect to the uncertainty

In the previous analysis, we considered the level of uncertainty estimated from Italian data (Vaccini, 2021). Here, we extend our analysis by considering the effect of different levels of uncertainty on the performance of our optimization framework. We use as a baseline the level of uncertainty estimated from the Italian data, and, for the two sliding-horizon policies described in the above (weekly and monthly), we consider different scenarios in which the level of uncertainty is increased by two, four or ten times, or decreased by the same amounts —for Vaxzevria, we omit to report the results with tenfold uncertainty because, due to the high baseline uncertainty, the problem is often infeasible. Table 2 reports the variation of the average utility function and total number of doses administered with respect to the baseline scenario.

Table 2

Effect of the level of uncertainty on the performance of the optimization framework.

	Comirnaty				Vaxzevria
	weekly		monthly		weekly		monthly
	ΔJ	ΔV	ΔJ	ΔV	ΔJ	ΔV	ΔJ	ΔV
10σ	-10.32%	-21.04%	-8.05%	-8.26%	-	-	-	-
4σ	-3.17%	-7.37%	-2.89%	-0.97%	-1.66%	-20.27%	-1.70%	-21.90%
2σ	-0.95%	-2.50%	-1.03%	-0.16%	-0.44%	-5.89%	-0.55%	-5.95%
0.5σ	0.53%	1.28%	0.64%	2.10%	0.26%	2.12%	0.31%	1.63%
0.25σ	0.82%	1.87%	0.96%	3.21%	0.40%	3.25%	0.45%	1.27%
0.1σ	0.99%	2.19%	1.14%	3.79%	0.49%	4.01%	0.54%	1.25%

Predictably, increasing the uncertainty negatively affects the performance, while a decrease has a positive impact. Interestingly, the effect of such an impact is highly nonlinear and is affected by the characteristics of the vaccine and the re-evaluation policy implemented. In particular, it seems that our optimization framework is able to deal with increased levels of uncertainty without big performance losses up to a critical threshold. Once this critical threshold is reached, the performance quickly drops (as can be observed with a tenfold increase in the uncertainty for Comirnaty and with a fourfold increase for Vaxzevria, which has already a higher baseline level of uncertainty). Besides depending on the characteristics of the vaccine, the critical point seems also to be influenced by the re-evaluation policy. Specifically, policies with less frequent re-evaluations seem to be more resilient than policies with frequent re-evaluations. These observations should be carefully considered in view of our findings in Section 4.3. In fact, in the previous section, we showed that frequent re-evaluations of the optimal policy may be beneficial in the presence of higher levels of uncertainty. However, it seems that above a critical threshold for the uncertainty —fortunately, several times larger than the one observed in real data— the optimization framework fails in dealing with the uncertainty. A possible explanation could be that, at each re-evaluation, the uncertainty model for the expected supply is re-computed. In the presence of an extremely high uncertainty, this may lead to overfitting the data, and the advantage of having frequent re-evaluation is therefore lost, thus decreasing the performance. To sum up, our findings suggest that frequent re-evaluation of the vaccination policy may be a double-edged sword when the uncertainty level is very high and that the accurate estimation of the uncertainty model (based on historical data and on information given by the suppliers) is key toward optimally design the sliding horizon policy. However, our estimation of the uncertainty model for the COVID-19 vaccination campaign in Italy suggests that real-world scenarios have levels of uncertainty that are below such a critical threshold, and frequent re-evaluations are thus often preferable.

Sensitivity with respect to the risk parameter

One of the key parameters in our framework is the risk parameter , i.e., the probability of violating the constraints in the stochastic optimization framework. Practically, setting the value of means to determine the probability with which we allow the computed optimal solution to be unfeasible due to the stochasticity of the supply. Thus, large values of correspond to allowing higher risks of unfeasible plans that could result in problematic online steerings of the vaccination rollout. On the other hand, smaller values of would provide more conservative solutions, which are often safer but may yield unsatisfactory performance. Intuitively, policy makers should set the value of to trade off between potentially high performance and the risk of implementing problematic steerings. In the simulations above, we proposed optimal solutions using . Here, we explore the effect of different values of , spanning from very safe implementations with , to extremely risky scenarios, in which we set .

The results are reported in Table 3 , where we compare the performances with the one obtained in the safest scenario, (). Predictably, we observe decreasing performance as long as we approach safer scenarios. Interestingly, while increasing toward more risky strategies has a strong initial impact on the performance, such an effect seems to saturate, becoming smaller and smaller as we move to more risky strategies. In fact, while increasing from 0.005 to 0.1 increases the number of doses administered for Comirnaty with a weekly evaluation by , a further increase of the parameter by 0.1 would increase the doses by just a further . This suggests that policy makers can determine an optimal trade-off value for such a parameter, corresponding to the point at which the increase in the performance starts becoming not sufficiently large to justify the corresponding increases in the risk. Another interesting observation can be found when comparing the solutions with different sliding-horizon policies. In fact, the optimal solutions with monthly re-evaluation seem to be slightly more affected by changes in than those revised on a weekly basis. Hence, opting for vaccination strategies that are more risky could be more beneficial when policy makers cannot re-evaluate the strategy with high frequency. However, in the presence of high levels of uncertainty (i.e., with Vaxzevria) and with less frequent re-evaluations, the increase in the performance associated with the implementation of risky strategies seems to quickly saturate, suggesting that conservative scenarios may be preferred when the uncertainty in the supply is large. From the policy makers, point of view, this could constitute an interesting approach to evaluate the trade-off between the maximization of the vaccination and the risk of an infeasible vaccination plan.

Table 3

Effect of the risk parameter on the performance of the optimization framework.

	Comirnaty				Vaxzevria
	weekly		monthly		weekly		monthly
	ΔJ	ΔV	ΔJ	ΔV	ΔJ	ΔV	ΔJ	ΔV
ε=0.4	1.41%	3.55%	1.77%	4.04%	0.72%	7.41%	0.85%	5.84%
ε=0.3	1.23%	3.18%	1.58%	3.39%	0.63%	6.47%	0.75%	5.89%
ε=0.2	1.01%	2.75%	1.33%	2.45%	0.52%	5.39%	0.65%	6.22%
ε=0.1	0.71%	2.05%	1.00%	1.21%	0.36%	3.89%	0.49%	5.84%
ε=0.05	0.49%	1.46%	0.71%	0.46%	0.25%	2.77%	0.34%	4.53%
ε=0.01	0.14%	0.39%	0.10%	-0.19%	0.07%	0.77%	0.09%	1.34%

Conclusion

We proposed a novel optimization framework for dealing with the optimal design of a two-dose vaccination rollout in the presence of uncertain vaccine supply. We cast the optimization problem as a numerically-solvable SOCP, with an objective function that takes into account the vaccine efficacy and the partial immunity gained after a single dose. The framework incorporates realistic constraints on the administration timing of the vaccine —minimum and maximum period that should elapse between the doses, limited capacity of the healthcare system, and tunable parameters that model the uncertainty in the supply.

We demonstrated the proposed optimization framework on a case study calibrated on the COVID-19 vaccination campaign in Italy. The results illustrated the applicability of the proposed framework to design vaccination rollout, in real-world scenarios, thanks to its formulation as a SOCP. In particular, we showed that a sliding-horizon implementation could improve the performance, even in the presence of high levels of uncertainty. This case study also showed how the framework might be used by public health authorities to accurately design a vaccination rollout, by setting the re-evaluation frequency of the rollout and trading off performances and risk due to uncertainties in the supplies.

The flexibility of our framework paves the way for several extensions. First, further features may be incorporated (e.g., additional booster doses, or the combination of different types of vaccine) through additional state variables. Second, we considered a uniform population, where the vaccination of each individual provides the same contribution to the objective function. Partitioning the population in classes, depending on the individuals’ risk —and thus on the benefit of vaccinating them— may be a suitable strategy to use our framework for devising a vaccination rollout not only in terms of scheduling first and second doses, but also in terms of their prioritization. All these extensions could be implemented at the cost of slightly increasing the computational complexity, without changing the very fabric of the optimization framework. Finally, we modeled uncertainty by means of Gaussian noise in the cumulative supply, which captures delays (or earliness) of committed weekly supplies. Further sources of uncertainty, such as the cancellation of supply may be modeled by including additional features to the uncertainty model, such as a Bernoulli process for missing deliveries. These extensions, however, would call for a different formulation of the stochastic optimization problem and are envisaged for our future research.