Designing an optimal sequence of non-pharmaceutical interventions for controlling COVID-19.

The COVID-19 pandemic has had an unprecedented impact on global health and the economy since its inception in December, 2019 in Wuhan, China. Non-pharmaceutical interventions (NPI) like lockdowns and curfews have been deployed by affected countries for controlling the spread of infections. In this paper, we develop a Mixed Integer Non-Linear Programming (MINLP) epidemic model for computing the optimal sequence of NPIs over a planning horizon, considering shortages in doctors and hospital beds, under three different lockdown scenarios. We analyse two strategies - centralised (homogeneous decisions at the national level) and decentralised (decisions differentiated across regions), for two objectives separately - minimization of infections and deaths, using actual pandemic data of France. We linearize the quadratic constraints and objective functions in the MINLP model and convert it to a Mixed Integer Linear Programming (MILP) model. A major result that we show analytically is that under the epidemic model used, the optimal sequence of NPIs always follows a decreasing severity pattern. Using this property, we further simplify the MILP model into an Integer Linear Programming (ILP) model, reducing computational time up to 99%. Our numerical results show that a decentralised strategy is more effective in controlling infections for a given severity budget, yielding up to 20% lesser infections, 15% lesser deaths and 60% lesser shortages in healthcare resources. These results hold without considering logistics aspects and for a given level of compliance of the population.

Introduction

Influenza outbreaks have always posed significant challenges to societies because of their adverse impact on public health and economic development, resulting in widespread loss of lives and livelihoods. The novel coronavirus 2019 (SARS - CoV-2) outbreak originated in Wuhan, China in December, 2019 and was declared a pandemic by WHO on 11th March, 2020 after 114 countries were affected with the virus (WHO, 2020a). Coronaviruses (CoV) belong to a large family of viruses that cause a wide range of illnesses from common cold to severe respiratory disorders and other diseases (Cui, Li, & Shi, 2019). The novel coronavirus (nCoV) is a virus strain that has been newly identified in humans and has been subsequently named the ‘COVID-19 virus’ (WHO, 2020b). COVID-19 has infected more than 266 million people across 213 countries as of 6th December, 2021, with more than 5.2 million reported deaths (Worldometers, 2021) with a projected economic loss of 3.94 trillion USD Chammah, Ralph (2021). The last time the world witnessed a pandemic of such large proportions was back in 1918. The 1918 influenza pandemic, better known as ‘The Spanish Flu’ was caused by an H1N1 virus with genes of avian origin. The deadly virus outbreak originated in USA during the spring of 1918 and resulted in approximately 500 million infections, 50 million deaths across the world (CDC, 2019), lowering the average life expectancy in the United States by 12 years (Jester, Uyeki, & Jernigan, 2018). Other notable infectious disease outbreaks which have emerged in the last century are severe acute respiratory syndrome (SARS), Middle East respiratory syndrome (MERS), Ebola, Human Immunodeficiency Virus (HIV).

In this paper, we study the policymaker’s problem of managing an epidemic outbreak, balancing its sanitary and economic impacts, by factoring in finite healthcare capacity and economic output. In the absence of a vaccine, it is imperative for affected countries to implement containment measures to check the epidemic spread. These interventions which do not involve any kind of immunization or preventive medical action are called non-pharmaceutical interventions (NPIs). CDC (2020) defines NPI as ”actions, apart from getting vaccinated and taking medicine, that people and communities can take to help slow the spread of illnesses like pandemic influenza”. ECDC (2020) clearly mentions that until safe and effective vaccines are approved, NPIs are the only public health tools to slow the growth of a pandemic. Medicines and therapeutic treatments are curative solutions to treat people who have already been infected whereas vaccines are preventive solutions offering immunity against the virus (De Keersmaecker & Cassidy, 2021). NPIs can range from mild measures like self-isolation of symptomatic individuals to more severe ones like travel restrictions, ban on public gatherings, school closures, service closures, with the most stringent one being a complete lockdown. Each of these containment measures come with varying economic and social costs (Brzezinski, Kecht, Van Dijcke, 2020, Charpentier, Elie, Laurière, Tran, 2020, Perra, 2021, Sadique, Adams, Edmunds, 2008). There are immediate direct costs like lost revenue and profits and far-reaching indirect effects like increased unemployment and suicide rates. The stringency of these measures has an impact on the infection rate and helps in controlling the rapid multiplication of the virus. The basic reproduction number, , pronounced “R naught”, also called the basic reproduction ratio, is an epidemiological metric used to describe the contagiousness of infected individuals (Delamater, Street, Leslie, Yang, & Jacobsen, 2019). is one of the most fundamental metrics for the study of infectious disease dynamics (Pellis, Ball, & Trapman, 2012). Dietz, K. (1993), defines as “the number of secondary cases one case would produce in a completely susceptible population.” When the for an outbreak is greater than 1, it leads to an epidemic because of the exponential growth of the pathogen in the population (Diekmann, Heesterbeek, & Metz, 1990). The objective of implementing NPIs is to reduce the as much as possible. The trade-off lies in the fact that the most draconian NPIs with the highest reduction in are also the most expensive to execute and sustain.

The timing of implementing NPIs plays a crucial role in epidemic decision making. Pei, Kandula and Shaman (2020) reported that if a lockdown had been imposed a week earlier in USA, there would have been 36,000 fewer deaths due to COVID-19. Most countries have faced the extremely difficult dilemma of safeguarding lives versus livelihoods for controlling the pandemic. This excruciating dilemma is real and unavoidable (Loayza & Pennings, 2020). Different nations have adopted varying strategies for tackling this challenging situation. While countries like South Korea and Taiwan have resorted to extensive testing and contact tracing for isolating the infected, European countries, namely France, Germany, Spain and Italy have imposed restrictive lockdowns for checking the spread. On the other hand, countries like Brazil, Sweden and USA have been late to react to the epidemic situation allowing infections to increase rapidly. Another important consideration for policymaking has been the implementation of decentralised (localised) and centralised containment measures. Lockdowns are not only economically expensive, but lead to social fatigue (Ouardighi, Khmelnitsky, & Sethi, 2020) and political instability (BBC, Ward), as has been observed during the multiple waves of the COVID-19 pandemic. During the second wave, UK (Shaw, 2021) and Germany (Kleczkowski, 2021) have enforced localised restrictions in order to allow the least affected regions to mobilize the economy. Hence, the comparative effectiveness of decentralised and centralised measures will serve as an important insight for policy-making, for controlling a pandemic.

In this paper, we propose a mixed integer linear programming (MILP) formulation for modelling the progression and control of the epidemic. We develop optimization models to address two strategies - centralised and decentralised, and compare the effectiveness of both in controlling the epidemic. The models combine the epidemiological dynamics with operations research for devising an optimal policy for managing the COVID-19 pandemic. The research questions that we explore are the following:

What is the optimal sequence of non-pharmaceutical interventions (NPIs) for controlling an epidemic under a given severity budget, considering economic constraints and healthcare resource shortages?

What is the benefit of a decentralised strategy (differentiated decisions across regions) in terms of reduction of infections and deaths?

Based on the study of extant literature, the unique contributions of this paper are the following:

This paper is the first to implement a Mixed Integer Linear Programming (MILP) model for scheduling an optimal sequence of non-pharmaceutical interventions (NPIs) for epidemic control. We also consider various restrictions on lockdowns implemented based on social fatigue and economic viability. Contrary to other papers evaluating a restricted subset of scenarios or policies, we consider all potential sequences of NPIs over the planning horizon.

We compare the decentralised and centralised strategies for epidemic control and evaluate their effectiveness in checking the growth of infections and deaths under multiple scenarios, using the data of the 13 French regions.

We factor in the phenomenon of the infection spread affecting healthcare professionals as a part of the epidemiological model and consider shortages in healthcare capacity in terms of doctors, Intensive Care Unit (ICU) beds and regular beds (non-ICU). This is an important consideration related to COVID-19 as it has been observed in many countries at the peak of their infection spread that a high number of doctors and medical workers have succumbed to the infection and been inactive from service. On April 2, 2020, Spain was reported to have around 15,000 medical workers infected with COVID-19, making up 14% of the confirmed cases then (Nugent, 2020). They were self-isolating and unable to render their services.

We show that the optimal sequence of NPIs necessarily follows a decreasing severity pattern, i.e., the most severe measures like lockdowns are scheduled first in the optimal sequence for a given severity budget. Based on this property, we simplify the MILP to an Integer Linear Program (ILP), which is computationally more efficient and accessible to the policymaker.

We consider a crisis situation where no vaccine is available yet, which can be a long period as evidenced by the COVID-19 case; however we propose a possible extension of the model for the case with vaccines available (Online Supplement).

Naturally, our decision model has some limitations. In addition to the assumptions of the epidemic model itself detailed in Section 3.1, we can identify the following main limitations:

Epidemiological data are constant over the planning horizon in the model, whereas the knowledge on virus transmissibility and fatality may evolve over the period, given the high number of new infections each week and additional information that can be collected from them. This limitation is inherent to the fact that the government cannot change decisions every week like with dynamic models, but needs to announce a plan over several weeks, to enable the population, companies, restaurants, schools etc. to anticipate and schedule their activities in the everyday life. This was the case in France where the government’s lockdown plan was in general announced for 6 to 10 consecutive weeks.

The severity budget is a mean (or total) budget over the population and not an identical budget per person, so this mean budget can lead to different severity levels across individuals or regions (more severity for regions in a bad situation with strained healthcare resources is compensated by less severity for less risky regions, for a given total budget). This assumption plays a role in the decentralized policy resulting in lesser infections and deaths than a centralized one. However, allowing discrepancy in the severity of measures may pose compliance issues and feelings of unfairness among the population, although differentiation of severity level has already been implemented in several countries (eg.,lock-downs only applied to large cities or most crowded regions). We deal with the compliance issue in the Online Supplement of this paper.

The remainder of the paper is organized as follows. In Section 2, we review the literature in two parts, we first look at papers with optimization models addressing operational challenges related to epidemic management and then analyse recent COVID-19 papers on containment strategies specifically. In Section 3, we introduce the epidemic compartmental model and notations. In Section 4, we explain the Time-based optimization models for the decentralised and centralised strategies. In Section 5, we discuss the Sequence-based optimization model. In Section 6, we report the numerical results and the relative performance of both strategies considering data from France. In Section 7, we discuss the managerial implications of the strategies to inform policy-making. In Section 8, we summarize our findings and propose future extensions of our work.

Literature review

Here, we review the literature in two different sections; first, we focus on epidemic logistics papers in the field of Operations Research, encompassing various operational issues faced during an epidemic outbreak and second, we discuss papers exclusive to COVID-19 control.

Epidemic logistics models

We refer to a couple of key papers which provide an exhaustive review of the different types of problems and methodologies implemented in epidemic logistics research. Dasaklis, Pappis and Rachaniotis (2012) review the literature for epidemic containment and characterize three major streams of research: first, where pharmaceutical containment measures like vaccines are planned and executed, second, where non-pharmaceutical interventions like school closure, lockdowns are imposed, third, where a combination of both are deployed. The authors also classify the research streams into pre-epidemic, post-epidemic and integrated categories while looking at various logistical issues encompassing different solution approaches like mathematical programming, game theory, queueing theory, data simulation etc. for stochastic and deterministic problems. Dimitrov and Meyers (2010) provide a detailed comparison of various mathematical approaches for analysing epidemic spreads. They infer that compartmental models are less complex computationally, compared to contact network models and agent-based simulations, although the latter are more accurate in capturing the epidemic progression.

Operational problems like facility location layout, distribution and logistics of medical aid during an epidemic outbreak have been widely covered in this area of research for both influenza epidemics and man-made epidemics (bio-terrorism). Here, we discuss some of the related papers. Jia, Ordonez and Dessouky (2007) analyse a maximal covering facility location problem of a large scale emergency response for medical supplies distribution by implementing different optimization solution approaches. Anparasan and Lejeune (2019) determine the optimal number, size and location of medical facilities and analyse the deployment of resources to attend medically ill patients, based on the cholera outbreak of Haiti in 2010. Wang, Wang, and Zeng (2009) present a multi-objective stochastic optimization model using genetic algorithms and Monte Carlo simulations for constructing an emergency epidemic logistics network. Liu and Liang (2013) propose a three level dynamic optimization model for allocating medical resources in epidemic controlling, using the epidemic diffusion model for forecasting demand, whereas Ekici, Keskinocak and Swann (2014) develop a food distribution network model using Mixed Integer Linear Programming during an influenza pandemic. Different types of objective functions and problem structures have also been studied in epidemic logistics. For example, He and Liu (2015) address a novel epidemic control objective function, where they minimize psychological suffering for infected individuals who remain undiagnosed or untreated because of relief supply shortage, using an epidemic optimization model. Liu and Zhang (2016) model medical resource allocation in the light of an influenza outbreak (time discretized SEIR model) and divide their analysis in different phases for forecasting, planning and execution. Yin and Büyüktahtakın (2021) propose a multi-stage stochastic optimization model for allocation of hospital resources under uncertainty, for the Ebola epidemic. Other types of epidemic outbreak include planned bio-terrorist attacks where a pathogen is deliberately released in a population, for eg. smallpox (Dasaklis, Rachaniotis, Pappis, 2017, Kaplan, Craft, Wein, 2002), anthrax (Craft, Wein, Wilkins, 2005, Wanying, Alain, Angel, 2016), influenza (Longini, Halloran, Nizam, & Yang, 2004). However, there have been very few papers analysing non-pharmaceutical interventions for controlling an epidemic. Yaesoubi and Cohen (2011) develop a Markov Decision Process framework for dynamic policymaking for handling an influenza outbreak considering social distancing and vaccines as the two major types of interventions. The paper extends the SEIR epidemiological model with transition probabilities and defines actions as NPIs along with a reward function which captures the health benefits and economic costs owing to infections and interventions.

The non-linearity of the epidemic compartmental model poses challenges for mathematical programming approaches and has been approximated linearly for tackling epidemic logistics problems. Büyüktahtakın, des Bordes and Kıbış (2018) devise a Mixed Integer Linear Programming (MILP) model to study the spatial spread of the Ebola epidemic and derive the optimal number, time and location of healthcare facilities under budget constraints. This paper relaxes the bi-linear transmission dynamics of susceptible and infected individuals and uses a linear epidemic equation avoiding the computational complexities arising out of non-linearity. Liu, Xu, Cao and Zhang (2020) use a similar relaxation of the epidemic equations as Büyüktahtakın et al. (2018) and propose a Mixed Integer Linear Programming model for an H1N1 epidemic where the decision variables are the number of healthcare facilities to open and close for minimizing fatalities, subject to budget constraints. Other linear epidemiological-logistics models can be found in Yin and Büyüktahtakın (2021) in the context of multi-stage stochastic MILP, and Yin, Büyüktahtakın and Patel (2021) for ventilator allocation.

We build on the simplified epidemic equations as introduced by Büyüktahtakın et al. (2018) and furthered by Liu et al. (2020), because it enables tractability in the MILP approach. We show in Appendix B that the linear epidemic model is a very good approximation of the non-linear model. Although we borrow this linear approximation of Buyuktahtakin’s model, a lot of the other attributes are different and unique. We model the implementation of NPIs for controlling the epidemic and compute the gain of decentralisation over a centralised containment strategy, instead of modelling facility location decisions. We consider differing transmission rates for the different regions based on population density and assume a common national budget for managing the economic costs of implementing NPIs for all regions combined. There are other papers in the past which have linearized the epidemic equations for ensuring tractability of the models. Reveller, Lynn and Feldmann (1969) study the resource allocation problem for minimizing infections for a tuberculosis outbreak by implementing an integer programming model. They linearize the epidemic equations by assuming an average number of active infected people during the planning horizon. Zaric and Brandeau (2001) model interventions for controlling an epidemic spread and maximizing quality-adjusted life years of the population under budget constraints. The authors argue that the non-linear functions lead to increased computational complexity and hence propose approximations and heuristics to tackle them.

COVID‐19 Containment

We now look at some of the influential papers which have analysed the impact of various containment strategies in the context of the COVID-19 pandemic. Perra (2021) provides an exhaustive review of recent studies on COVID-19 with a focus on the application of non-pharmaceutical interventions, analysing 348 papers across different methodologies (optimal control, agent based models, statistical analysis, mathematical programming etc.) encompassing both analytical and empirical work based on different countries for implementing and lifting NPIs. The author highlights the challenges of limited information with respect to NPIs and the need to develop more robust theories for capturing its long terms effects and behavioral dimensions. Ferguson et al. (2020) carry out a comparative study between a counterfactual, unmitigated model, mitigation and suppression for the population of UK for different scenarios. They consider the implementation of five NPI measures individually and in combination, but do not consider their sequences. For the mitigation approach, the healthcare capacity is overwhelmed multiple times. For suppression, the combination of all NPIs at the same time yields the highest impact on transmission reduction. Walker et al. (2020) from the same Imperial College team extend the analysis to contrast mitigation and suppression strategies for Low, Middle and High Income countries. Loayza and Pennings (2020) argue that owing to a bigger unorganized sector in low income countries, the compliance to stringent lockdowns will be lower in comparison to higher income countries. Loayza, N. (2020) contends that indiscriminate lockdowns can be counter-productive leading to displacement of labourers and further aggravation of the spread. Also, when compliance to lockdowns is low, it results in low marginal gains and eventually leads to second and third waves of infections. Barnett-Howell and Mobarak (2020) point out that the economic value generated by comparable social distancing measures is approximately 240 times larger for the United States, 70 times larger for Germany compared to the value created in Pakistan or Nigeria, according to their analysis.

Eichenbaum, Rebelo and Trabandt (2020) analyse the interaction between epidemic modelling and economic decisions using an agent based modelling technique for various containment scenarios. They look at “smart containment” where the infection condition of the entire population is known to the policymaker with the help of extensive testing and hence enables the maximization of social welfare of the population by using a tailored intervention policy. Alvarez, Argente and Lippi (2020) analyse the policymaker’s problem to handle the trade-off between output costs of lockdown policies and fatality costs because of infections by combining the SIR epidemiology framework in an optimal control model. The paper looks at various degrees of lockdown from 60% of the population being confined to 10% of the population with a parameterized effectiveness of the lockdown policy, and the subsequent economic impact in terms of GDP contraction. They argue that if the lockdown is not very effective and healthcare capacity is not ramped up, the severity and duration of the lockdown should be lowered. Using optimal control, Ouardighi et al. (2020) compare the effectiveness of mobility restrictions and secured social interactions for epidemic control, factoring in ill-preparedness of the policymaker and population social fatigue; Caulkins et al. (2021) model the optimal intensity of lockdown, whereas Charpentier et al. (2020) model the optimal lockdown level, level of testing efforts for Type 1 (antigen) and Type 2 (antibody), by considering the combined impact of the epidemic on healthcare, economic and social costs.

Mehrotra, Rahimian, Barah, Luo and Schantz (2020) propose a multi-period mixed integer linear programming model for allocation of ventilators to different regions from the centre, based on stochastic demand for the COVID-19 pandemic. Yin et al. (2021) develop a data-driven stochastic multi-stage optimization model for allocation of ventilators in different regions under asymptomatic uncertainty and risk, for COVID-19. Acemoglu, Chernozhukov, Werning and Whinston (2020) study a Multi-Group SIR model for three different demographic groups- young, middle-aged and old and devise a customized lockdown strategy for minimizing economic damage and infections. The authors suggest a differential lockdown policy with targeted shielding of the elderly and testing, tracing, isolation of all symptomatic patients for improving the economic and health impacts of lockdowns. Di Domenico, Pullano, Sabbatini, Boëlle and Colizza (2020) propose a stochastic age structured epidemic transmission model for capturing the progression of the epidemic in the Ile-de-France region in France, evaluating possible exit strategies after the lockdown. Davies et al. (2020) carry out a similar study for United Kingdom, and analyse various NPIs for controlling the epidemic across 186 counties. They compare local and national control strategies and consider various triggers for initiating lockdowns based on healthcare capacity.

Note that in the above papers, when government policies are evaluated for controlling Covid-19 (lockdowns or other NPIs), only a few scenarios are generated, evaluated (often with optimal control techniques) and compared, but no optimization approach has been proposed so far for designing the best sequence (or optimal mix) of NPIs over a planning horizon among the huge set of all potential sequences. We did not find either a comparison of decentralised vs centralised strategies for epidemic control. In the next section we present our epidemic and decision models.

Epidemic model and notations

Epidemic model

In this section, we introduce the compartmental model for computing the epidemic spread of COVID-19 in a given population, based on the classical model of Kermack and McKendrick (1932). A large number of recent articles (more than 40 papers with studies across 15 countries) have adopted the SIR compartmental model or its extensions for analysing the control of COVID-19 using NPIs (Perra, 2021). We divide the entire population into the compartments of Susceptible, Infected, Critical, Recovered and Dead for a given time period . The Susceptible compartment comprises of individuals who are vulnerable to the virus but have not been infected yet; Infected represents the active infections in the population at time ; Critical represents the number of individuals who require any kind of hospitalization; Recovered constitutes the number of people who have been cured, and Dead denotes the fatalities owing to COVID-19 specifically. We consider a discrete-time and deterministic epidemic model. We assume each time unit to be equivalent to one week accounting for weekly reviews and situational assessments from the policymaker’s perspective.

Figure 1 depicts the compartmental model schematically. The individuals in the susceptible compartment transition into the Infected compartment at a rate of . The infected population comprises of non-severe (mild or asymptomatic) and severe infections, where the latter require hospitalization. We denote by , the proportion of severe infections. Severe infections enter the Critical compartment at a rate of per week, from symptom onset. Non-severe infections recover at a rate of per week. From the Infected compartment, people transition into either the Critical compartment for hospitalisation or to the Recovered compartment directly, if non-severe. Individuals in the Critical compartment who are hospitalized, recover or die at a rate of . We assume that when a patient requires hospitalisation (Critical), he is admitted to a regular bed or an ICU bed (subject to availability). is the proportion of the hospitalised individuals who recover and is the proportion of people who die. For the schematic model, we represent the rates and proportions for ICU and non-ICU patients to be the same. However, in the numerical analysis we consider different rates based on empirical observations. Dead and Recovered are terminal compartments with no exit. Conversely, the Susceptible compartment has no inflow from other compartments.

Fig. 1

Epidemic Compartmental Model.

The assumptions for the epidemic compartmental model are the following:

We assume the transmission between Susceptible and Infected to be dependent only on the Infected population of the previous period for linearizing the epidemic state equation (Büyüktahtakın et al., 2018). In Appendix B we show empirically that this assumption of a linear epidemic model is not restrictive compared to the non-linear epidemic model.

The infection transmission rate is a function of the population density of the region. Hence, densely populated clusters will have a higher rate of spread than sparsely populated regions. (Gerritse, 2020) empirically shows that when population density doubles the transmission rate increases by 0.06 points.

When an infected individual recovers, he is no longer susceptible to the virus, since the proportion of reinfection cases with respect to the total is negligible (less than 1%) according to recent findings (Hall et al., 2020). Another very recent empirical study by Flacco et al. (2021) further strengthens this finding, reporting a reinfection rate across multiple age groups and cases with and without co-morbidities.

Patients requiring hospitalisation avail either a regular bed or an Intensive Care Unit (ICU) bed, and not both (Di Domenico et al., 2020).

A shortage of regular beds and doctors delays the recovery process and a shortage of ICU beds leads to deaths. Lin (2021) reports that COVID-19 mortality rates are strongly associated with ICU bed availability from an empirical study and we base our assumption on this.

Deaths are accounted for only the infections which require hospitalization (Di Domenico et al., 2020).

A given proportion of severe symptom infections self-isolate. Based on a recent study in the UK (Smith et al., 2021), it was reported that only 18% of COVID-19 patients self-isolated themselves.

Model notations

In Table 1 , we summarize the various NPI levels the policymaker can implement for controlling the pandemic and in Table 2 we describe the epidemic model parameters. In Table 3 , we describe the optimization model notations and constant parameters, and in Table 4 , we define the primary and auxiliary decision variables for the model.

Table 1

NPI List.

Parameter	Description
Self-Isolation (I)	Isolation or Quarantine of vulnerable/ symptomatic people for 7 days
Travel Restrictions (T)	Travel restrictions within and across regions based on distance.
School Closure (S)	Closure of schools and universities - primary, secondary, UG, PG.
Public Gathering Ban (G)	Restriction on gatherings for public events >50 people
Full Lockdown (L)	Full service closure, curfew, mass restriction of movement
NPI Combinations
I	Level 1 NPIs
I+T, I+G, I+S	Level 2 NPIs
I+T+S, I+T+G, I+S+G	Level 3 NPIs
I+T+S+G.	Level 4 NPIs
L	Level 5 indicates Lockdown

Table 2

Epidemic Parameters.

Parameter	Description
β_k	Weekly transmission rate of infections per infected person in region k
βkH	Weekly rate at which doctors are infected in region k
ρ	Percentage of infected people who become Critical (need hospitalisation)
μ	Weekly rate at which a non-severe infection recovers
ϕ	Weekly rate at which an Infected person enters the Critical state
γ	Weekly Rate at which a Critical infection leaves the hospital
λ	Percentage of hospitalised people who recover
r^H	Percentage of infected doctors who become re-infected

Table 3

Optimization Model Notations.

Parameter	Description
Indices
t	Index for time period (weeks), t∈{1,…,T}
i	Index for NPI level, i∈N={1,2,3,4,5}
k	Index for region, k∈{1,…,K}
r	Index for medical resource,r∈Res={1,2,3},1 = Doctors,2 = ICU,3 = Regular (Beds)
Other parameters
τ¯	Average time to quarantine for doctors after being infected.
c_i	Severity cost of implementing NPI i per period and per person.
α¯	Maximum reduction in transmission rate, α¯=Max{α_i|i∈N}
m	Number of patients attended by one doctor
b_2,k	Number of ICU beds in region k
b_3,k	Number of Regular beds in region k
θ	Proportion of Critical patients requiring ICU beds
η_rkt	Proportion of resource r allocated for non-COVID cases in period t, region k
l_0	Maximum number of changeovers of NPIs
l_1	Maximum number of weeks of lockdown over the T weeks
l_2	Maximum number of consecutive weeks of lockdown
l2′	Minimum number of weeks between two successive blocks of lockdowns
B	Upper bound on the average severity of NPIs per individual, per period (budget).
V_kt	Theoretical Lower bound of Infections for period t, region k
I_max	Theoretical Upper bound of Infections for all periods.
ϵ	% of Critical patients who self-isolate.
P	Total population of the country
p_k	Population of region k
I_k1,C_k1,H_k1	Epidemic state values for period 1 for region k, ∀k∈{1,…,K}

Table 4

Decision variables for MINLP and MILP models.

Variables	Description
PRIMARY
x_kit	=1 if NPI level i is selected for week t and region k, 0 otherwise
AUXILIARY
Epidemic
I_kt	Number of Infected people at week t in region k
C_kt	Number of Critical people at week t in region k
H_kt	Number of active hospital doctors at week t in region k
Shortage
s_rkt	Number of units of resource r (beds/doctors) in shortage at week t in region k
z_0,kt	=1 if the number of untreated patients due to shortage of doctors is less than the number of untreated patients due to shortage of regular beds, 0 otherwise, at week t in region k.
z_rkt	=1 if resource r demand exceeds supply, 0 otherwise, at week t and region k
Q_0,kt	Max between number of untreated patients due to shortages in doctors and non-ICU beds at week t in region k
Q_rkt	Min between number of patients needing (demand) and actually availing (capacity) resource r at week t in region k
Linearization of quadratic terms
X_kit	Replacement for x_kit(1−α_i)(1−ρϵ)I_kt
Y_0,kt	Replacement for H_ktz_1,kt
Y_rkt	Replacement for C_ktz_rkt
U_rkt	Replacement for s_rktz_0,kt,∀r∈{1,3}
Changeover
W_kit	=1, if NPI level i is not selected in week t−1 and selected in week t, 0 otherwise

As indicated in Table 1, we group the NPIs in different levels, from 1 to 5. For each level (for Level 2,3,4 since Level 1 and 5 have only one combination), we assign a single combination corresponding to the best reduction factor based on our estimations and other empirical studies (Ferguson, Laydon, Nedjati Gilani, Imai, Ainslie, Baguelin, Bhatia, Boonyasiri, Cucunuba Perez, Cuomo-Dannenburg, et al., Haug, Geyrhofer, Londei, Dervic, Desvars-Larrive, Loreto, Pinior, Thurner, Klimek, 2020, Lin, 2021). Level 5 corresponds to a lockdown, which is the most stringent containment measure possible.1

Before delving into the model further, we take a closer look at the variables used, as listed in Table 4. The primary decision variable is the NPI level, . The epidemic state variables defined by the compartmental model are the auxiliary variables. Out of the state variables, we consider only Infected and Critical (Dead can be expressed as a function of Critical), since the rest do not directly impact the policy maker’s decisions. We also model the epidemic impact on doctor availability with variable . The set of variables captures the shortage in capacity of medical resource (doctors, ICU, non-ICU beds) at week and region , in terms of the respective resource () units. For example, indicates the shortage of doctors and is the number of COVID-19 patients untreated due to lack of doctors, where is the number of patients one doctor treats. Similarly, and indicate the shortages in ICU and non-ICU beds, respectively, where one bed corresponds to one patient, for both. The combination of the and variables are used to ensure that is exactly equal to the shortage in capacity of medical resource .

For analysing the optimal policy making for containing the spread of COVID-19, we first introduce a Mixed Integer Non Linear Programming (MINLP) optimization model. We linearize the bilinear terms in the MINLP model subsequently to convert it into a Mixed Integer Linear Programming (MILP) model. We further simplify it using model properties into an Integer Linear Programming (ILP) model and compare the computational efficiency of both in the numerical section. We refer to the MINLP and MILP models as ‘time-based’ because the NPI decisions are associated with each week, whereas we refer to the ILP model as ‘sequence-based’ because entire sequences of NPIs spanning the planning horizon are decided.

Time‐based optimization models

Decentralised MINLP model

Here, we define the Time-based Mixed Integer Non-Linear Programming (MINLP) model for the decentralised strategy, with the aim to minimize the total number of unique infections:

Subject to:

In our model, we capture the shortage in healthcare resources and evaluate its impact on the epidemiological state variables. The maximum between untreated patients due to lack of doctors and regular beds gets added to the number of critical patients and the untreated patients due to lack of ICU beds gets added to the number of dead patients, influencing the NPI level selected by the policymaker subsequently.

The objective function in (1) minimizes the sum of unique infections over the weeks. In constraint (2) we impose that exactly one NPI can be chosen for each region in every week . In constraint (3) we ensure that the average severity of NPIs over the whole population and the horizon of weeks stays within the national severity budget . This severity budget is an indicator of the maximum economic loss that the policymaker is willing to afford for safeguarding the population’s health. The metric can be also be interpreted as labour hours lost due to containment measures. We assume that the economic impact of the contact reduction is a convex, quadratic function of as analysed by (Charpentier et al., 2020), i.e., . Constraints (4)–(6) imply that the shortage due to unserviced healthcare demand of medical resources in terms of doctors, ICU Beds, regular beds is defined by (Resource demand – Resource supply, 0). To ensure that the binary variables truly indicate if the demand of resource exceeds supply, we introduce another set of continuous variables defined as the minimum between demand and supply for resource . We add three constraints to achieve this. The first constraint computes = Capacity + (1-) Demand. The second and third constraints imply is less than the capacity and demand, respectively. Hence, when demand is greater than capacity, for to satisfy the constraints of being lesser than both capacity and demand, , so that exactly equals capacity. In this way, we ensure that the variable indicates when demand is greater than capacity, resulting in equalling the minimum value between capacity and demand of resource . For example, let us assume that the demand for doctors, and the availability (capacity) of doctors is . Since we have demand greater than capacity, because of the constraints with . This results in a shortage of doctors. However, if the availability of doctors was greater than demand, then and shortage would be 0. The same applies for all resources . A similar type of constraint has been discussed in Büyüktahtakın et al. (2018). Constraints (7) to (9), (10) to (12), (13) to (15) are the system of constraints for medical resources’ demand and supply for all resources . Constraints (16) to (18) compute the maximum between the number of patients untreated due to a shortage in doctors, and due to a shortage in regular (non-ICU) beds. Constraints (17) and (18) imply that is at least equal to the number of untreated patients due to the lack of doctors and due to regular beds, respectively. Hence, when the number of untreated patients due to doctors is greater than the number of untreated patients due to a lack of regular beds (non-ICU), we have , so that exactly equals the untreated patients due to doctors. Since shortages due to a lack of either of them leads to untreated patients, instead of adding them we take the maximum of the two. For example, if we assume one doctor treats 10 patients and there is a shortage of 2 doctors along with a shortage of 30 non-ICU beds, the number of untreated patients due to lack of doctors would be , and the total number of untreated patients would be .

Constraints (19) to (21) are the epidemic state equations based on the Kermack and McKendrick (1932) compartmental model. The choice of NPI level for time period and region , reduces the new infections by a reduction factor for level . A proportion of people with severe infections self-isolate and the rest are subjected to the NPI (self-isolation is common to all NPIs). We introduce the parameter for capturing the proportion of non-severe infections who need hospitalization. Hence, a proportion of non-severe infections recover directly (We test for ). The number of patients who need doctors and regular beds and remain untreated due to a shortage of either or both () are added back to the Critical (C) compartment, indicating that their recovery is prolonged, whereas the number of patients who need ICU support and cannot avail it due to a shortage in ICU beds are added to the Dead compartment, indicating their demise. Constraint (21) captures the state equation for hospital doctors for time period and region . In constraint (21), we assume that the infections in doctors are directly impacted by the number of actual hospitalizations in regular and ICU beds and the number of patients attended per doctor. We consider the minimum of demand and capacity to be precise, as captured by variables and respectively. We define as the rate of infection for doctors. Based on an empirical study by Abo-Leyah et al. (2021), healthcare workers have been observed to be three times more likely to contract the infection than the general population. Although healthcare workers are expected to take precautionary measures like wearing personal protective equipment and taking other medicinal safeguards, their risk of contracting the virus is high due to continuous exposure to infected patients. The number of recovered doctors from the last period is added back to the current period’s level of available doctors, in proportion to the effective recovery rate . Based on the empirical study by Dimeglio et al. (2021), we also account for the reinfection of doctors due to repeated exposure. The authors observed in their study that of the doctors were reinfected after a period of 167 days of being infected for the first time, indicating that the proportion of reinfection among medical professionals is very low. Hence, we add a reinfection parameter by which, doctors already infected once in period (period of recurrence observed was around 5.5 months, hence we estimate weeks) get infected again and are removed from the current period’s level of available doctors. Constraints (22), (23) and (24) are integrality and non-negativity constraints on variables. Note that we do not include the state equations corresponding to Susceptible, Recovered and Dead in the model as they do not have constraints involving them (we model deaths as a separate objective function subsequently). They can be computed ex-post based on the NPI decisions and infections resulting from it, as follows: (hospitalised patients).

We also test another objective function of minimizing the total number of deaths caused by the pandemic. We compute deaths due to the epidemiological progression, and due to ICU bed shortages captured by variable . The time horizon considered is weeks because the impact of the NPI in week () impacts the Infected state in , Critical infections of week , which in turn affects the deaths in week . The objective function becomes:

For the centralised model, we replace the NPI decision variables with because of homogeneous decisions across regions, the assignment constraint (2) is no more for each region and the sum over disappears in the budget constraint (3). The rest of the constraints are identical with the decentralised model. We still compute the Infected, Critical and the corresponding shortage variables at a regional level because of differentiated infection rates and capacities across regions.

MILP reformulation

In order to convert the former MINLP model into a Mixed Integer Linear Programming (MILP) model, we linearize the quadratic parts of the model in the following way. We define four sets of auxiliary variables for linearizing the products of , , respectively.

Linearization of :

Linearization of :

Linearization of :

Linearization of :

We use the linearization method proposed by Glover (1975), as implemented in Büyüktahtakın et al. (2018), where we add auxiliary variables and constraints. We introduce constraints (26) to (29) for linearization of the bi-linear terms in the objective function. We define and as the lower bound and upper bound, respectively, for the corresponding continuous variable . The rationale behind defining the lower bound as time dependent is to make it tighter and improve the computational time. Hence, we derive the lower bound for total active infected people at week by considering that the strongest NPI (Level 5) has been applied in all periods leading up to . We introduce the auxiliary variable as a replacement of . We define the lower bound for by . We add constraints (30) to (33) for linearizing the bilinear terms corresponding to available doctors in region and time , given by . Finally, we implement constraints (34) to (37) for linearizing the product of and and constraints (38) to (41) for linearizing the product of shortage variables and the indicator variable . We assign the lower bound for as 0 and the upper bound as for each time period and region .

Additional constraints on NPI sequences

In this subsection we introduce three scenarios S0, S1 and S2, of additional policy constraints in order to make the containment policy enforceable and sustainable from an implementation point of view. In the last two years, we have seen various containment plans being implemented by different countries for COVID-19. For example, after the second lockdown there was a phased exit with relaxations in the severity of measures every days for France (Thompson, 2021), whereas in UK a five week gap was announced between each stage of relaxing restrictions (Manning, 2021). We capture these aspects of policymaking in our model in Scenario . The other phenomenon that we consider is resistance to lockdowns, which has been observed in different parts of the world like the US, some countries in Europe and Brazil (Ward, 2020). With a fourth wave impending in Europe, there has been a resurgence of protests due to the announcement of new lockdown measures (BBC, 2021). Hence, we introduce different restrictions on the number of weeks of lockdown imposed, in Scenarios and , to account for this.

Scenario S0:

Restricted changeovers. No more than changeovers across weeks:

Two week rule: If an NPI is selected, it should be on for at least two consecutive weeks, whenever applied.

Scenario S1: No more than weeks of lockdown

Scenario S2: Maximum weeks of consecutive lockdowns, Minimum gap of weeks in between phases

Scenario is the default model, where we consider restrictions on NPIs in general. In constraints (42) to (44), we put a limit on the number of changeovers in NPIs across the planning horizon. We impose this condition in order to make the policy practically sustainable. Switching between NPI modes from 1 to 5 rapidly will be chaotic to enforce and populations will find it difficult to comply with. In constraints (42) and (43) we record the number of changeovers between NPIs using binary variables . In constraint (44), we put an upper bound on the number of changeovers. Using Constraint (45) and (46) we ensure that every NPI, whenever it is introduced, stays on for two weeks at least.

For Scenario , constraint (47) ensures that the total number of lockdowns is less than , as set by the policymaker. Constraints (48) and (49) together represent Scenario , which implies that there cannot be more than successive weeks of lockdowns (Constraint (48)), with a gap of at least weeks between two phases of lockdowns. We test the centralised and decentralised strategy for all scenarios , and separately.

Sequence‐based optimization model

In this section, we present a simplification of the MILP model based on nice properties of optimal solutions stated in Propositions 1 and 2, that are corollaries of the following lemmas.

For a given feasible solution, if for some region the NPIs at two periodsand, with, are swapped, then the infection level in periodstays the same.

Note that Lemma 1’s result is about the level of infections at a given week due to the inventory balance equations, and not the cumulative number of new infections until that week, which is not the same with a NPI switch.

To illustrate this we discuss a simple example. In Fig. 3 , we show two NPI sequences for and , where the first is the optimal one and the second is a modification of the optimal sequence where the NPIs in week (NPI level 5) and week (NPI level 2) have been swapped. From the graph on the left in Fig. 2 we can observe that the infection level in week is the same for both the optimal and swapped solution showing that Lemma 1 is satisfied. In the figure on the right we compare the cumulative infections for the optimal and swapped solutions and show that it is clearly higher for the swapped solution. This gives us the intuition behind Lemma 2. We also use the result of Lemma 1 to derive Lemma 2 (See Appendix A)

Fig. 3

NPI sequences for optimal and swapped solutions.

Fig. 2

Infection level (left) and cumulative infections (right) for the optimal and swapped solution.

If a solution satisfiesandwithand(i.e NPIis less effective than), then making the NPI swapandwill strictly improve the objective value (for both minimization of infections and deaths).

The proofs of the lemmas can be found in the Appendix. We now introduce the two propositions (corollaries of Lemma 2) that enable to turn the time-based model into a tractable sequence-based model.

For Scenarioand, the optimal NPI sequence is always ordered from higher severity to lower severity (for both minimization of infections and deaths).

For Scenario, the optimal NPI sequence always follows an ordered sequence from higher to lower severity for every sub-sequence containing a lockdown, and for non-lockdown NPIs across sub-sequences (for both minimization of infections and deaths).

Using Propositions 1 and 2, one can switch to a tractable ‘sequence-based’ model, where a pre-processing can easily compute the set of all sequences satisfying the property of descending severity, together with the budget constraint and the constraints of scenarios and . Doing so, the state variables and shortage variables become constants in the model. The cost of each sequence is given by where binary parameter is equal to 1 if NPI level is assigned to week in sequence , 0 otherwise. represent the state values corresponding to sequence at week and region . , and represent the min(demand,capacity) values and shortage values for healthcare resources, respectively. As a result of this simplification, all the linearization variables of Section 4.2 are discarded.

We discuss only the decentralised strategy, because this simplification turns the centralised optimization problem into an enumeration problem (generate all sequences in , evaluate the objective value of each and pick the best one). The sequence-based decentralised formulation is still an NP-hard Integer Linear Programming model, though much simpler than the earlier MILP formulation. For this ILP model we define a new decision variable: and the sequence-based model can be written as follows:

Minimization of Infections

Minimization of Deaths

Subject to:

The objective function (50) minimizes the sum of infections for the planning horizon, whereas objective (51) minimizes deaths. Constraint (52) implies that for each region , only one sequence is assigned from the feasible set . Similarly, Constraint (53) is the budget constraint, weighted by population, factoring in the NPI cost for each sequence . Constraint (54) defines the variable domain. We demonstrate the advantage of the sequence-based formulation over the time-based formulation in terms of computational efficiency in the next section.

Numerical analysis

In this section we discuss a case study based on the COVID-19 spread in the metropolitan region of France in Europe. We compare the containment models for both decentralised and centralised strategies by considering various lockdown constraints separately as scenarios (, , , as defined earlier). For the decentralised strategy, we study the epidemic spread in the 13 metropolitan regions of France. For our analysis, we have incorporated data from French national archives and databases for demographics and healthcare capacities. We consider population density, healthcare capacity in terms of regular beds, ICU beds, general practitioner for all 13 regions. Based on Salje et al. (2020), and Di Domenico et al. (2020), we calibrate the various epidemic parameters pertaining to the regions and at a national level. Naturally, the decentralised model can be implemented at a finer granularity than the region level, e.g. at French “départements” or city level, however we chose a division of the country into regions since we had available data for the 13 regions, and the benefit of decentralisation is already shown with no ambiguity at this level. The gain of more decentralization would be further improved if one applied it at a very fine granularity like cities or districts (since we would give even more flexibility for solutions), naturally with higher implementation difficulties.

Calibration of parameters

For quantifying the impact of various NPIs on the reduction of infection transmission and , we resort to social contact matrices. Social contact matrices divide the population into age groups and define the per day contact between different age groups. Several empirical studies have been carried out in the past to record actual social contact patterns in different countries. Prem, Cook and Jit (2017) developed social contact matrices for 152 countries based on contact surveys and demographic data. Béraud et al. (2015) carried out one of the first large population contact based surveys for France and computed social contact matrices for modelling infection spread in the population.

Based on these studies we consider contact matrices for the French population for age groups (0–19, 20–39, 40–59, 60 -75, 75+) and modes of contacts such as home, school, work, other locations and calibrate the contacts for different NPIs. For consideration of reduction in contacts for different NPIs our assumptions are closely related to Di Domenico et al. (2020) and Davies et al. (2020). We report the final contact matrices considered for our analysis, in Appendix D, in terms of resultant contacts after implementation of the corresponding level of NPI. We quantify the reduction in contacts in different settings like Home, School, Work, other locations depending on the severity of the NPI level (1–5). For example, when an NPI includes school closure, we set the contacts made in the school setting to 0, and marginally reduce the other travel related contacts, whereas we marginally increase the contacts at home considering children spend more time at home. A reduction in social activity also implies a subsequent drop in economic activity. Hence, the implementation of an NPI, while ensuring a reduced transmission rate also reduces the labour productivity. Charpentier et al. (2020) consider that the economic loss due to lockdowns is a quadratic function of the fraction of population who are unable to work because of it. Other recent papers such as Richter and Throckmorton (2020), Djidjou-Demasse, Michalakis, Choisy, Sofonea and Alizon (2020), assume a quadratic cost impact of NPIs. Our characterisation of the drop in productivity and consequent economic loss because of NPIs is similar. After computing the impact of each NPI on the economic output, we fit it on a scale of for a simpler representation. Note that each of the NPI levels correspond to the degree of combination of NPIs, from 1 to 4, with 5 indicating a Lockdown. In Appendix C, we explain the details of NPI calibration. We also report the epidemic parameter values for France, based on French government archives and empirical papers in Appendix D.2

In order to reduce the computational complexity for executing the models, we have not considered age specific variables. Instead, we have factored in weighted averages for the age-centric epidemic parameters and contact data. For scenario we assume and for scenario we assume . For changeovers, we consider . For initiating the progression of the epidemic, we considered 0.004% of the population of each region to be infected already at the beginning of the planning horizon. We assume .

Gain of decentralization

The MILP and ILP models for decentralised and centralised strategies were executed on CPLEX 12.10 on a Mac OSX, 8 GB RAM system. In Table 5 , we report the results of the optimization model for minimization of infections, for a time horizon of 8 and 10 weeks, three different severity budgets (), and for each of the scenarios , and . We test the same for minimization of deaths in Appendix D. Not surprisingly given its higher flexibility, the decentralised strategy consistently exhibits a better objective function value in all scenarios and instances compared to the centralised strategy, with a reduction in infections up to 20% and reduction in deaths up to 15%. We can also observe that the objective function values follow the relationship in all instances: . This is logical because in scenarios and we add restrictions on lockdowns which results in more infections and deaths in comparison to .

Table 5

Gain of Decentralisation for Inf* (* = Optimization Criterion).

Scenario S0		Centralised			Decentralised			Variation
B	T	#Inf*	#Deaths	Time(s)	#Inf*	#Deaths	Time(s)	ΔInf	ΔDeaths
5	8	19,370	473	0.55	18,036	439	6.2	-7%	-7%
6.5	8	10,826	290	0.6	9819	276	6.4	-9%	-5%
8	8	7146	223	0.49	7008	222	6.4	-2%	-0.4%
5	10	25,746	621	1.12	22,251	553	66.5	-14%	-11%
6.5	10	10,906	308	0.84	10,221	293	65.3	-6%	-5%
8	10	7432	238	0.83	7038	230	65.3	-5%	-3%
Scenario S1		Centralised			Decentralised			Variation
5	8	19,370	473	0.78	18,036	439	6.4	-7%	-7%
6.5	8	10,826	290	0.53	10,028	284	6.2	-7%	-2%
8	8	8273	253	0.92	8273	253	6.0	0%	0%
5	10	25,746	621	2.49	22,251	553	66.5	-14%	-11%
6.5	10	10,906	308	1.09	10,537	303	68.2	-3%	-1.6%
8	10	8499	268	0.98	8499	268	67.2	0%	0%
Scenario S2		Centralised			Decentralised			Variation
5	8	19,370	473	0.73	18,244	444	10.9	-6%	-6%
6.5	8	11,591	321	0.8	10,821	310	11.2	-7%	-3%
8	8	10,219	301	0.88	8747	279	11.0	-14%	-7%
5	10	28,736	688	2.37	23,226	591	222.3	-19.2%	-14%
6.5	10	15,103	420	1.41	12,608	374	225.2	-17%	-11%
8	10	11,236	362	1.07	9780	321	223.1	-13%	-11.3%

Analysis of deaths due to ICU bed shortage

In Table 6 , we depict a particular instance, where we assume a skewed spread of initial infections, considering regions with higher population (5 millions) have a higher initial infections to population ratio in comparison to the rest. We consider the objective of minimization of deaths, for in this case. The table indicates the peak ICU demand to capacity ratio for the 8 weeks at a regional level. When the ratio is more than , it implies that there was a shortage in ICU beds, which resulted in as many deaths. In terms of NPI allocation, for the centralised strategy, we can see that every region gets the same allocation of severity index. As a result, there is a very high skew in demand to capacity ratio, as high as for Ile-de-France, whereas in the decentralised model, the ratios stay within . This leads to a much higher level of shortage in ICU beds in the centralised strategy, and hence more deaths. The total deaths due to shortage in the decentralised strategy is 1096, in comparison to 2704 for the centralised strategy, lesser. The weighted peak demand to capacity ratio is for Centralised and 119% for Decentralised strategy. In terms of allocation of NPIs, if we observe Bretagne and Normandie regions, both have the same population of 3.3 million, however, Bretagne’s severity index of NPI is 3.4, in comparison to 2.62 for Normandie. This can be explained by the fact that, Bretagne has a population density of 122.8 people per square kilometres and an ICU density of 6.9 beds per 100,000 people, in comparison to Normandie’s population density of 111.6 people per square kilometres and ICU density of 8.3 beds per 100,000 people. Population density directly impacts the infection spread, hence a higher population density coupled with lower ICU density indicates a greater risk, necessitating more stringent measures for safeguarding the sanitary situation. Hence, these kind of allocations at the regional level in the decentralised model result in a better management of the pandemic, with lesser infections and deaths, for the same economic loss, as compared to the centralised strategy.

Table 6

Heterogeneity of NPI severity & ICU shortage across regions (Scenario , B=4, T=8, Deaths*).

Regional Data					Peak Demand/Cap.		#Deaths (Shortage)		Severity
PD	ID	P	Region	Cap.	C	DC.	C	DC	C	DC
39.7	6.1	0.3	Corse	16	88%	138%	0	6	3.98	3.25
58.2	11.5	2.8	Bourgogne-Franche-Comté	240	59%	106%	0	14	3.98	2.95
71.4	7.9	6.0	Nouvelle-Aquitaine	400	141%	140%	162	190	3.98	3.85
81.5	10.0	5.9	Occitanie	480	116%	116%	79	75	3.98	3.85
65.4	6.8	2.6	Centre-Val-de-Loire	160	88%	119%	0	31	3.98	3.40
96.0	7.6	5.5	Grand Est	320	173%	120%	262	63	3.98	4.45
115.2	7.1	8.0	Auvergne-Rhône-Alpes	480	151%	150%	246	321	3.98	3.85
111.6	8.3	3.3	Normandie	240	57%	123%	0	55	3.98	2.62
122.8	6.9	3.3	Bretagne	160	101%	138%	2	60	3.98	3.40
118.5	7.5	3.8	Pays de la Loire	240	75%	114%	0	33	3.98	3.25
161.0	7.9	5.1	Provence-Alpes-Côte d’Azur	320	173%	118%	250	59	3.98	4.45
188.0	11.2	6.0	Hauts-de-France	560	121%	83%	119	0	3.98	4.45
1022.2	10.2	12.3	Île-de-France	1040	231%	118%	1584	189	3.98	4.79
				4656	146.4%	119%	2704	1096	3.98	4.0

*P=Population (in millions), PD=Population Density, ID=ICU Density, Cap=Capacity, C=Centralised, DC=Decentralised.

Gain of sequence‐based model

In Table 7 , we report the number of sequences generated for the ILP model compliant with Proposition 1 for scenarios , and with Proposition 2 for scenario . the table shows there is a reduction of from the total number of possible sequences. In Table 8 , we compare the computational efficiency between the time-based (MILP) and sequence-based (ILP) formulations. We can check that both formulations yield the exact same objective value for infections. However, there is a drastic reduction in CPU time, as high as in comparison to the time based MILP formulation.

Table 7

Number of sequences generated for scenarios .

Time	#Total Sequences	#Feasible Sequences in Seq
T	Total=5^T	S0	S1	S2
8	390,625	115	106	100
10	9,765,625	225	206	180

Table 8

Computational gain of the Sequence-based model (Inf*).

Scenario S0		Time based		Sequence based		CPU Time (s)
B	T	#Inf*	#Deaths	#Inf*	#Deaths	Time based	Sequence based	ΔTime
7	6	8109	228	8109	228	11.4	3.32	-70.9%
6.5	6	9076	248	9076	248	44.9	3.29	-92.7%
6	6	10,372	270	10,372	270	1081.3	3.3	-99.7%
7	7	8341	242	8341	242	11.18	6.23	-44.3%
8	8	7008	222	7008	222	3814.2	11	-99.7%
8.5	10	6630	222	6630	222	1820.1	295.1	-83.8%

Sensitivity analysis

In Fig. 4 , we highlight the optimal NPI sequences for scenarios for both objectives of infections and deaths minimization, for decentralised and centralised strategies respectively. We observe that the optimal NPI sequences for the two objectives are different. For the minimization of deaths, the total number of lockdowns across regions are greater with a more drastic shift to lower severity NPIs towards the end of the planning horizon, in comparison to the minimization of infections. In Fig. 5 , we demonstrate the evolution of infections for scenarios and for the centralised strategy, performing a sensitivity analysis on NPI budget for weeks. It can be clearly seen that the infections explode exponentially with a fall in budget . An unlimited budget would naturally allow continuous lockdowns over the weeks, which would not be accepted by the population for a large . Also, for , the pattern follows alternating rises and falls, owing to the consecutive lockdown condition, with higher peaks and troughs in comparison to . Finally, in Appendix D, we perform a sensitivity analysis based on , the reduction factor, by varying it in a range of the estimated values. We test it for scenario , three budgets and time horizon weeks. These three budgets and 8 correspond to different values of average severity per individual, between Level 4 corresponding to , and Level 5 corresponding to (the cost is non-linear with the NPI level; see Appendix D for more details on the scaling of severity costs and calibration of coefficients used in the calculus of ). We observe that the decentralised strategy still performs better for minimization of infections (up to 10%).

Fig. 4

Differentiated NPI decisions across regions, (, B=6.5, T=10).

Fig. 5

Evolution of the optimal number of infections with budget .

Managerial insights

Based on the output of the optimization model, we deduce some key takeaways for policymakers. Although there is a degree of ambiguity surrounding the epidemic parameter values and the exact economic and health impact of NPIs, we can identify some broad trends.

The first key finding is with respect to the optimal sequence of NPIs. For both the decentralised and centralised strategies, the NPIs follow a pattern of decreasing severity across the planning horizon. Although other studies like Ferguson et al. (2020), Walker et al. (2020) have established how a suppression strategy can control the explosion of infections, they compare the impact of different NPIs on infection spread across the planning horizon, considering each of them are applied throughout the duration in separate scenarios. However, our analysis further extends it by taking into account sequences of multiple NPIs and showing that a descending pattern of implementation is most effective. This indicates that the policymaker needs to impose heavier restrictions in the beginning, in order to check the transmission in the population and utilise that time to ramp up healthcare capacity to prevent a health crisis. Initially, a few countries like UK had contemplated implementing the herd immunity strategy (Titheradge & Kirkland, 2020) which would allow a sizeable proportion (around ) of the population to be infected and develop antibodies against it. However, the sharp increase in hospitalizations and overwhelming of healthcare capacity did not allow them to go ahead with it. Hence, the optimal strategy for controlling COVID-19 would be to put heavier NPIs in the beginning followed by less severe ones, for minimizing infections or deaths within a given budget constraint.

The second takeaway is the benefit of a decentralised strategy (without considering logistics aspects and for a given level of compliance of the population) for a given mean budget of severity. Based on our numerical analysis of the French regions, we observe that the decentralised strategy yields better results for all instances, with up to lesser infections and lesser deaths. It is natural that the objective value improves when giving more flexibility to the model, but the gain appears to be quite substantial here. This is an important consideration for policymaking in a pandemic, where vaccine development takes time and NPIs need to be sustained over a longer horizon to keep the infection transmission at a manageable level. Exhaustive confinement and lockdowns at a national level leads to social fatigue and is economically unsustainable. Given these limitations, it is pragmatic for the policymaker to follow a more customised strategy at a regional/ state/ city/ district level, considering the demographics, healthcare capacity and the transmission dynamics observed initially. We conclude that a decentralised strategy not only results in better performance in terms of epidemic control but also aids in mobilising the economy wherever feasible. Of course, with a severity constraint in terms of mean budget, the price to pay is the discrepancy of severity across individuals and regions, as hinted in the introduction. This geographical differentiation has already been implemented in France (in Ile-de-France region, the city of Nice, the ”département” of Moselle, to cite a few examples), and in other countries.

Finally, a very crucial aspect of pandemic management is to identify capacity risks and bottlenecks and take proactive measures. Through our numerical experiments we show that regions which have a relatively higher population density and lower healthcare resource availability per unit population are at greater risk than other regions (For e.g, Comparison between Normandie and Bretagne in Table 6). Since we know that a higher population density leads to higher infection transmission (Gerritse, 2020), the pressure on the healthcare system is evidently higher. This necessitates stricter measures in order to prevent shortages in healthcare resources and a consequent spike in infections and deaths. Hence, with the help of such scenario analyses the policymaker can foresee which regions need to be prioritised based on their available capacity and potential case load and take measures in advance. Our analysis of the decentralised and centralised strategies considering shortages in ICU beds, reveals how the former leads to a better allocation of NPIs, resulting in lesser shortages (up to 60%) and consequently lesser deaths. This policymaking tool allows the policy think tank to identify capacity risks based on various scenarios and increase capacities accordingly to avoid an unfortunate overwhelming of the healthcare system.

Conclusion

In this paper, we develop an epidemic optimization model for controlling COVID-19 infections in a given population, considering budget constraints on the set of NPIs decided by the policymaker, and healthcare resource shortages. The classical SIR epidemic compartmental model of Kermack and McKendrick (1932) is incorporated with certain extensions in order to capture the dynamics of COVID-19. We factor in healthcare resource shortages in terms of doctors, ICU beds, and regular beds. We linearize the initial MINLP time-based model, and then use the property (which we establish analytically) of the decreasing severity of NPIs over time to convert it into a ‘sequence-based’ ILP model, which can be solved in a much faster way (time reduction of up to 99%). Two strategies are analysed - centralised and decentralised, where NPIs are modelled for the entire country for the former, and at a regional level for the latter. We use social contact matrices from empirical studies and actual epidemic and demographic data from French archives for developing a case study on the spread of COVID-19 in France. Our results indicate that the decentralised model performs better in terms of lesser infections (up to 20%), lesser deaths (up to 15%), lesser shortages (up to 60%) for a given total or mean budget. Despite the evolving nature of the pandemic and the uncertainty around data, these findings offer some direction to the policymaker for implementing NPIs optimally, weighing both the sanitary and economic impact of their decisions. We believe that this analysis will be valuable not only from the perspective of COVID-19, but also for managing future occurrences of epidemics of any kind, by leveraging some of the insights as rule of thumb. In terms of limitations of our research, although our analysis establishes the decentralized strategy to be better, compliance of populations with too many differential restrictions and severity levels might pose challenges to the administration. Our motivation is to demonstrate the benefits of a decentralized strategy from a scientific standpoint, and the solutions obtained do not show too chaotic patterns across regions in any case. We discuss an extension to our model in Appendix E (Online Supplement) where we factor in population compliance and numerically show the conditions for which a decentralised containment strategy is more effective than a centralised one.

We suggest some further extensions to this research. First, application of machine learning methods in combination with mathematical programming can be used for a more data driven approach to policy making, where an adaptive rolling plan of NPIs can be modelled based on population compliance data. Data obtained from the first wave of the pandemic can be used to train algorithms to predict compliance of the population to NPIs in subsequent waves and build more robust, predictive models for policymaking. Second, capacity increase, facility location and layout, are some of the other dimensions of the epidemic control problem which can be explored. Third, for the sequence based model, column generation can be implemented when the number of various NPIs is large, in order to tackle the explosion of sequences.

Finally, this paper restricts to non-pharmaceutical interventions at the early stage of the spread for a new virus like COVID-19, before efficient vaccines can be found. An extension to the MILP model can be made for considering vaccine administration along with NPIs simultaneously. A considerable number of approved vaccines have emerged since December, 2020 for inoculating the global population against COVID-19, namely Pfizer/BioNtech, Moderna, AstraZeneca, Johnson and Johnson, and Sputnik V (Healthdata, 2021). One could modify the objective function applying to the rate of infections a discount factor based on the percentage of an age group vaccinated with dose 1 and dose 2 at each week and region, along with the effectiveness of each dose (% of people getting sufficient immunity by dose 1 and 2). We introduce this model extension in Appendix E (Online Supplement) where we account for the speed of vaccine roll-out and the compliance of population to NPIs along with it. However, another extension to the vaccine problem could be explored where the quantity of each dose made available at each time period in each region is modelled as a supply chain problem, considering lead time, capacity, logistics and budget constraints. A new decision model could then indicate the optimal vaccine production and distribution plan along with the optimal NPI sequence.

Table 9

Comparison of Infection levels between Linear and Non-Linear epidemic models.

S0	B=3			B=1.5
t	Linear	Non-Linear	ΔInf	Linear	Non-Linear	ΔInf
1	3265	3265	0%	3265	3265	0%
2	3602	3602	0%	4401	4401	0%
3	3987	3987	0%	5956	5956	0%
4	4428	4428	0%	8097	8096	0.01%
5	4935	4934	0.01%	11,059	11,057	0.02%
6	7375	7374	0.02%	16,573	16,565	0.04%
7	11,085	11,081	0.04%	24,978	24,958	0.08%
8	16,761	16,750	0.07%	44,269	44,205	0.14%
9	25,497	25,468	0.11%	79,017	78,808	0.26%
10	39,029	38,957	0.19%	142,056	141,357	0.49%

Table 10

NPI Calibration based on Social Contact Matrices.

Parameter	Description
g,g^′	Index for age group, g,g^′∈{1,2,3,4,5}
i	Index for NPI level, i∈N={1,2,3,4,5}
a	Index for location category,a∈A={1,2,3,4},1 = Home,2 = School, 3 = Work, 4 = Others
d_ia	Ratio of contacts for NPI level i to baseline for location a.
o	Baseline contact.
w_a	Weight of contact location a, a∈A.
c_igg^′.	Number of contacts between age group g and g^′ for NPI level i
cagg^′o	Number of contacts between age group g and g^′ in baseline scenario o, for location a.
p_g	Population of age group g