Data-driven Optimized Control of the COVID-19 Epidemics
DData-driven Optimized Control of the COVID-19 Epidemics
Afroza Shirin , Yen Ting Lin , and Francesco Sorrentino Department of Mechanical Engineering, University of New Mexico, New Mexico, 87131 Electrical and Computer Engineering, University of New Mexico, New Mexico, 87131 Information Sciences Group, Computer, Computational and Statistical Sciences Division(CCS-3), Los Alamos National Laboratory, New Mexico, 87544September 9, 2020
Abstract
Optimizing the impact on the economy of control strategies aiming at containing the spread of COVID-19 is a critical challenge. We use daily new case counts of COVID-19 patients reported by local healthadministrations from different Metropolitan Statistical Areas (MSAs) within the US to parametrize a modelthat well describes the propagation of the disease in each area. We then introduce a time-varying controlinput that represents the level of social distancing imposed on the population of a given area and solvean optimal control problem with the goal of minimizing the impact of social distancing on the economy inthe presence of relevant constraints, such as a desired level of suppression for the epidemics at a terminaltime. We find that with the exception of the initial time and of the final time, the optimal control inputis well approximated by a constant, specific to each area, which is in contrast with the current system ofreopening ‘in phases’. For all the areas considered, this optimal level corresponds to stricter social distancingthan the level estimated from data. Proper selection of the time period for application of the control actionoptimally is important: depending on the particular MSA this period should be either short or long orintermediate. Finally, we see that small deviations from the optimal solution can lead to dramatic violationsof the constraints.
The fast propagation of the COVID-19 pandemic has attracted unprecedented attention from both the publicand the scientific community. Due to the high fatality rate of SARS-CoV-2 [1], governments throughout theworld have adopted measures such as lock-down, stay-at-home, and shelter-in-place, which in turn have ledto substantial economic losses, see e.g., [2]. In many countries control interventions have been articulated inphases, usually phase 1 to phase 3, with higher phases corresponding to progressively lower restrictions [3]. Afundamental challenge is to balance the need to suppress the spread of COVID-19 and the need to contain the1 a r X i v : . [ q - b i o . P E ] S e p conomic impact of measures aiming at limiting the spread of the disease. In this manuscript, we apply optimalcontrol theory on a mathematical model for the propagation of the epidemics, parametrized by real-world datadescribing different regions, and compute control strategies which are optimal for each region from an economicstandpoint.A number of papers and reports have focused on both modeling and controlling the pandemic. Flaxman etal. [4] looked at the effect of non-pharmaceutical interventions including school closures, banning of mass gather-ing, social distancing, etc. on the reproductive number R of COVID-19. Sanche et al. [5] used a mathematicalmodel with data on individual cases, real-time human travel, and infections, as well as estimated epidemiologyparameters to compute R and found that it is higher than initially estimated. Chang et al. [6] adopted anagent-based model to determine the efficacy of several intervention strategies on the spread of COVID-19 inAustralia. Anderson et al. [7] performed the Bayesian analysis to estimate the impact of social distancing onnumber of reported cases and hospitalizations in British Columbia and found that when 78% participation in so-cial distancing has been accomplished the cases would decrease; it also noted that if the participation were below45%, an exponential growth would restart. Morris et al. [8] explored COVID-19 intervention methods which arerobust to implementation errors and found that these methods in conjunction with optimal time-limited methodsderived from the standard SIR model can be used to mitigate the spread of the virus. Another study analyzedan open-loop optimal control policy updated weekly using real-world feedback, and found that this method iseffective in reducing fatalities even if some measurements are inaccurate [9]. A study published in March 2020estimated the ICU occupancy and ventilator use from a statistical model under the conditions of social distancingand found that the demand for hospital beds and ventilators will exceed the current supply [10]. Another studyexplored optimal policies for decreasing economic cost and mortality rates from a multi-risk SIR model and foundthat strict lock-down policies which specifically target the elderly population were most effective in minimizingdeaths and economic losses [11]. Previous work has not computed optimal controls for data driven models. Onthe one hand, Refs. [6], [8], [9], [11] explicitly compute optimal control strategies, but their models were stylizedand not parametrized/calibrated by data; On the other hand, Refs. [4] [5], [7], [10] used data to parametrize themodels, but the models do not have controllers which can be used to infer optimal control strategies.In what follows, we first construct a mathematical model, which is parametrized by historical and regionaldaily new case report. After parametrization, the data-driven model is capable to reproduce the regional progres-sion of the COVID-19 epidemics up to the present. Then, we apply optimal control theory to the parametrizedmodel to compute an optimal control strategy over the course of a pre-determined period into the future tosuppress the epidemics to a desired level, while minimizing economic costs. This type of approach is suitablefor long-term planning (i.e., over the course of several months) as opposed to short-term planning, which can bedifficult to implement by the government and by businesses.While it appears likely that an effective vaccine for SARS-CoV-2 is under way, most estimates indicate thatthe vaccine will be distributed to the population in summer or fall 2021 [12], which points out the need forplanning interventions to contain the epidemics for several months to come in the absence of a vaccine. Thus ourproposed workflow aims to bridge the gap until the time T vac when an effective vaccine is discovered, massivelymanufactured, and administered to the majority of the population. Another temporal consideration regards the2ime T herd at which a population achieves herd immunity in the absence of a vaccine and without overflowingthe medical facilities. While herd immunity from COVID-19 is the subject of much ongoing discussion [13],in the Methods we provide a rough estimate of T herd from available data. In this paper we proceed under theassumption that the inference period T inf and the control horizon T cont are such that T inf + T cont < T vac and T inf + T cont < T herd .We set out the analysis by first introducing a compartmental model which describes key features of theCOVID-19 epidemics. The whole population is divided into the following compartments. The susceptible popu-lation compartment (S) includes the people who can contract the pathogen SARS-CoV-2. The exposed population(E) are those who have been infected but have not progressed long enough into the disease to transmit it tosusceptible people. Those who can transmit the disease (‘carriers’) are divided into the asymptomatic group(A) who do not show symptoms and the infected symptomatic groups (I) who show symptoms. Both the symp-tomatic and asymptomatic groups can transmit the disease, but with different infectiousness—the asymptomaticpeople are less infectious. The infected symptomatic population I is divided into three sub-compartments.The first sub-compartment I sq includes those who just self quarantine and do not get tested. The second sub-compartment I tp includes those who get tested, result positive, and get quarantined. In contrast to the previoustwo sub-compartments, the third sub-compartment I s includes the rest of the symptomatic population, who donot change their behavior until they recover (or decease). Those who were tested positive and those who decidedto self-quarantine are moved into a quarantined compartment (Q), and stop interacting with other populations.The removed compartment (R) includes those who are completely recovered from the disease and have acquiredimmunity, and those who have died because of the disease. Both groups are not susceptible to reinfection andare removed from the system.Mathematically, the time evolution of the population density—defined as the compartmental populationnormalized by the total regional population—of each compartment is governed by the following set of coupledordinary differential equations,˙ S ( t ) = − βP ( t ) S ( t ) [ I s ( t ) + I sq ( t ) + I tp ( t ) + µA ( t )] (1a)˙ E ( t ) = βP ( t ) S ( t ) [ I s ( t ) + I sq ( t ) + I tp ( t ) + µA ( t )] − λE ( t ) (1b)˙ A ( t ) = λ (1 − σ ) E ( t ) − γ A A ( t ) (1c)˙ I sq ( t ) = p sq λσE ( t ) − [ γ I + γ sq ] I sq ( t ) (1d)˙ I tp ( t ) = p test λσE ( t ) − [ γ I + γ tp ] I tp ( t ) (1e)˙ I s ( t ) = (1 − p sq − p test ) λσE ( t ) − γ I I s ( t ) (1f)˙ Q ( t ) = γ tp I tp ( t ) + γ sq I sq ( t ) + p sq λσE ( t ) − γ I Q ( t ) (1g)˙ R ( t ) = γ A A ( t ) + γ I [ I s ( t ) + I tp ( t ) + Q ( t )] , (1h)where β is the transmissibility, λ is the transition rate from the exposed compartment to either the asymptomaticor symptomatic compartments, γ I is the transition rate from the infected compartment to the recovered com-partment, γ A is the transition rate from the asymptomatic compartment to the recovered compartment, µ is the3elative infectiousness of asymptomatic individuals (compared to symptomatic individuals), σ is the fraction ofexposed people who transition into the symptomatic compartments, p sq is the fraction of symptomatic peoplewho will self-quarantine, γ sq is the transition rate from from the I sq compartment into the quarantined compart-ment, γ tp is the transition rate from the I tp compartment into the quarantined compartment, and p test is thefraction of infected who get tested (but only a fraction of them will be able to get a test and be confirmed as apositive, I tp ). The model assumes that testing resources are not scarce, i.e., availability of a testing kit to everyperson in I tp ; we also assume positive people are identified as such with 100% accuracy. The case with limitedtesting resources is discussed in the SI. Realistically, the time scale associated with γ test can be several days, sowe assume γ test = 0 . ≤ P ( t ) ≤
1, which measuresthe reduction of contact probabilities between the susceptible and the infectious populations (which include bothA and I). The model is structurally similar to the models analyzed in Refs. [7, 14], but simplified to allow forefficient calculations of optimal control solutions. Figure 1 illustrates a schematic diagram of the model.In order to reduce the dimensionality of the dynamical system, we introduce the following simplification:we treat γ sq as a very large number; for γ sq → ∞ , we assume that those becoming symptomatic immediatelytransition into the Q compartment, yielding the following reduced-order model:˙ S ( t ) = − βP ( t ) S ( t ) [ I s ( t ) + I tp ( t ) + µA ( t )] (2a)˙ E ( t ) = βP ( t ) S ( t ) [ I s ( t ) + I sq ( t ) + I tp ( t ) + µA ( t )] − λE ( t ) (2b)˙ A ( t ) = λ (1 − σ ) E ( t ) − γ A A ( t ) (2c)˙ I tp ( t ) = p test λσE ( t ) − [ γ I + γ tp ] I tp ( t ) (2d)˙ I s ( t ) = (1 − p sq − p test ) λσE ( t ) − γ I I s ( t ) (2e)˙ Q ( t ) = γ tp I tp ( t ) + p sq λσE ( t ) − γ I Q ( t ) (2f)˙ R ( t ) = γ A A ( t ) + γ I [ I s ( t ) + I tp ( t ) + Q ( t )] . (2g)The schematic diagram which fits the above model in Eq. (2) is shown in the Supplementary Note 1 1.To bridge the model and the data, we also integrate an auxiliary variable C I which counts the cumulativeconfirmed cases and evolves according to ˙ C tp ( t ) = γ tp I tp ( t ) , (3)We will fit ∆ C tp ( t ) := C tp ( t + 1) − C tp ( t ) to the new case counts reported on each day, detailed in the section onparametrization below.From an economic point of view, the measures of social distancing P ( t ) and quarantining Q ( t ) have veryelevated costs for society. For example, the US real gross domestic product (GDP) dropped by roughly one thirdfrom year to year in the second quarter of 2020 [2], due mainly to the COVID-19 pandemic, see e.g., [15].We thus formulate the following optimal control problem,min P J = c p (cid:90) t f t i − P ( t ) P ( t ) dt + c q (cid:90) t f t i Q ( t )1 − Q ( t ) dt, (4)4igure 1: Compartmental model corresponding to Eq. 2.
The transition from the S (susceptible) com-partment to the E (exposed) compartment is affected by the population densities S ( t ), I tp ( t ), I s ( t ), A ( t ) and bythe time-varying control input P ( t ). c p , c q ≥
0, subject to Eq. (2), the initial conditions, the terminal constraint E ( t f ) + I s ( t f ) + I tp ( t f ) + A ( t f ) < (cid:15), (5)and the following path constraint I s ( t ) + I tp ( t ) ≤ I max , (6)where I max is an upper limit on the number of infected people that can receive proper treatment in the hospitals. t i is the time at which we perform the inference and start optimizing the control action, t f is the final time ofthe optimization, the previously defined control horizon T cont = ( t f − t i ) . The objective function (4) takes into account an economic cost for social distancing with an appropriatecoefficient c p > c q >
0. We set the costassociated with social distancing per unit time to be modeled by (1 − P ( t )) /P ( t ). By choosing this functionalform, the cost is linearly dependent on the scale of the reduction (1 − P ( t )) when 1 − P ( t ) (cid:28)
1, and nonlinearlydiverges near total shut-down ( P ( t ) (cid:28) − P ( t ) to indicate that strictermeasures of social distancing affect more and more ‘essential’ workers, and so increasingly larger parts of theeconomy. For example, the cost of limiting large gatherings of people like concerts or sport events is lower than5he cost of limiting customers’ access to stores and restaurants. Similarly, imposing quarantine requires resourcesthat are linearly proportional to the quarantined population when Q is small, and diverge nonlinearly when thequarantined population is close to the entire population Q = 1. We model such a cost by the functional form Q ( t ) / (1 − Q ( t )).Both social distancing and quarantining have a cost associated with the lack of economic return generatedby limiting person-to-person interactions. It is reasonable to assume c q ≥ c p , as strict quarantining requiressupervision costs as well as costs due to lowered productivity [16] while social distancing only incurs costs dueto lowered productivity [17].We will also consider the alternative objective function,min P J alt = c p (cid:90) t f t i (1 − P ( t )) dt + c q (cid:90) t f t i Q ( t ) dt, (7)for which the integrand functions are linear in P ( t ) and Q ( t ). While the formulation (7) is mathematically simpler,it does not take into account the fact that stricter measures of social distancing may result in progressively largereconomic losses. The alternative objective function (7) is here mainly introduced in order to compare the resultsto those obtained with (4).The tunable parameter (cid:15) represents the desired suppression level for the epidemics at the final time t f .In general, selected values for (cid:15) may depend on a number of factors, such as the time horizon over whichoptimization is performed and the particular stage of the epidemics (initial exponential growth, intermediategrowth, or plateau.) A possibility is to require complete eradication of the epidemics, which corresponds tosetting (cid:15) = 1 /N , where N is the number of individuals in the population. However, the high basic reproductivenumber of COVID-19 makes eradication unlikely; instead we set (cid:15) to a small number indicating suppression of thedisease. The other constraint given by I max represents the need to contain the infected population below a giventhreshold at all times (or ’flatten the curve’.) In what follows, unless differently stated, we set the terminalconstraint (cid:15) = 10 − , which corresponds to imposing that the number of infected people is below one personper 100000 population. This number is derived from the guidelines of European countries about re-opening,see for example [18], indicating that reopening occurred at about 2 × − . Also, European countries haveofficial guidelines for reimposing stricter lock-down measures, see: [19], which sets a critical population equal to50 / × − . The values of I max are set regionally based on the capacity of the ICU beds in differentmetropolitan areas and are summarized in Table 3. 6 arametrization of the model for different US metropolitan areas. We partitioned the model parameters into two sets: a set of fixed parameters and a set of inferred parameterswhich are estimated by the daily case counts reported by regional health administration and registered in therepository curated by the New York Times [20].The fixed parameter sets includes S , λ , γ I , γ A , µ , σ , p test , and p sq . S is the regional total population, and weused the US Census Bureau-estimated regional population of each of the Metropolitan Statistical Areas (MSAs)or ‘cities’, which are delineated by the US Office of Management and Budget [21]. Lauer et al. [22] estimated themedian of the incubation period to be about 5.2 days, however, there is evidence that patients become infectiousroughly two days before the onset of symptoms [23], which corresponds to approximately a three-day progressioninto contagiousness . We thus set a rough estimation for λ to be 1 / γ I and γ A , which are the transition rate to recovery of the symptomatic andasymptomatic populations, are estimated to be 0.12 (1/d) [24] and 0.26 (1/d) [25]. The relative infectiousness µ is set at 0 . γ I , γ A , µ , and σ are consistent with a more complex model which was used to perform daily forecasts of thedisease spread [14]. We assume p test = 0 .
25 and p sq = 0 .
4, noting that these parameters were able to reproducethe infected population at the time of the inference (we estimated that about 15 to 20% of the total populationinfected in the New York City MSA on 07-Jul-2020 when we parametrize the model.)We used the data from 21-Jan-2020 to 07-Jul-2020 to infer the inferred parameters, which corresponds to thepreviously defined inference period T inf . We assume that the social distance function P ( t ) before 07-Jul-2020 ispiecewise-linear to avoid over-fitting (due to observation noise). At the time when the analysis was performed,multiple MSAs had shown clear second-phase resurrection of the epidemics [14, 28]. We found that a two-phasepiecewise linear function is sufficient to reproduce the data of each of the MSAs: P ( t ) = , t ≤ t , , t < t ≤ t p − t − t ( t − t ) , t < t ≤ t p , t < t ≤ t p + p − p t − t ( t − t ) , t < t ≤ t p t < t < t i , where t is the time when the disease was introduced into a specific MSA, the monotonic t , t , t , and t ( t j ≤ t j +1 ) are the times when the social-distancing behavior changes, and p and p are the two social-distancingstrengths. We define ∆ t j := t j − t j − for j = 1 , . . .
4. The time at which we perform the inference is t , whichis also the time after which optimization of the control action begins. The two-phase piecewise linear modelis the minimal model that we found capable to reproduce the data, and can be validated by more rigorousmodel-selection procedure detailed in [14].To fit the model by the noisy daily report new counts, we adopt a negative-binomial noise model. We briefthe procedure here, noting that the procedure is similar to the inference method detailed in [14]. Given a set7f parameters θ (a stylized notation of the set of the inferred model parameters), the model Eq. (2) predictsa deterministic trajectory of the new positive tested case on day j , ∆ I tp ( j ; θ ). This deterministic prediction isinterpreted as the mean of a stochastic outcome, modeled by a negative binomial distribution NB( r, p j ) where r and p j are the parameter of the distribution, and p j is constrained by the deterministic model prediction p j = rr + ∆ C tp ( j ; θ ) . (8)Here, the r is the dispersion parameters which describes how disperse the noise distribution is; in the limit r → ∞ , the negative binomial statistics converges to Poissonian, and in the limit r → N daily reported new case counts { ∆ C tp ( j ) } N − j =0 can be formulated [14]: L (cid:16) θ ; { ∆ C tp ( j ) } N − j =0 (cid:17) ≡ N − (cid:89) j =0 (cid:18) ∆ C tp ( j ) + r − C tp ( j ) − (cid:19) p rj (1 − p j ) ∆ C tp ( j ) . (9)In summary, the inferred parameters θ include the regional-specific disease transmissibility β , onset of thedisease spread t , behavioral switching times t , t , t , t , the strengths of two social-distancing episodes p , p ,and a dispersion parameter r of the negative binomial noise model. These parameters were inferred by the dailycase reports from 21-Jan-2020 to 07-Jul-2020. With the formulated likelihood function (9), we used the standardMarkov chain Monte Carlo procedure (MCMC) with an adaptive sampler ([29], detailed in [14]) to estimate themaximum likelihood estimator of the parameters θ for each of the interested MSA’s.Table 1 summarizes the set of model parameters that were estimated for each MSA.8 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - , , , , , A N e w c a s ec o un t s NYCDataParametrized model - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - , , , , , B N e w c a s ec o un t s LADataParametrized model - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - , , , , C N e w c a s ec o un t s HoustonDataParametrized model - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - D N e w c a s ec o un t s SeattleDataParametrized model
Figure 2:
New case counts from January 21, 2020 to July 8, 2020 in four Metropolitan StatisticalAreas within the US.
Plus signs are daily new case counts reported by local health administrations. The solidline is the daily new cases obtained by integrating Eq. (2) after parametrization of the model.9able 1: Parameters by Metropolitan Statistical Area (MSA)Parameter New York Los Angeles Houston Seattle S λ γ I γ A µ σ p test p sq T D r β t t t t t P P esults: Optimal Control Solutions Following previous work by the authors [30, 31], we use pseudo-spectral optimal control (see Supplementary Note7) to compute solutions for the problem formulated in Eqs. (2),(4),(5),(6). We set t i = 169 (days) and t f = 259,corresponding to a ninety day control horizon, ( t f − t i ) = 90. We parametrize the solutions in I max , (cid:15) and c q ,after setting without loss of generality c p = 1.We focus on four different US metropolitan statistical areas: New York City (NYC), Los Angeles (LA),Houston, and Seattle. NYC is the largest US metropolitan area; it emerged as the main early hotspot of theepidemics in the US in March and April, but since early June has achieved stable control of the epidemics. LA isthe second largest metropolitan area in the US, it has seen a steady rise in the number of cases by the time whenwe performed the inference. Houston has seen a rapid increase of the cases in June and July. Seattle was also avery early hotspot, which has seen a decrease in the number of cases in April and May, followed by an increasein June and July. We chose these four MSA’s to cover a wide range of different dynamics before the time of theinference.The solution of the optimal control problem is the function of time P ∗ ( t ) that minimizes the objective function(4) subject to the constraints (2), (5), and (6). Different from the observation period t ≤ t i , for which we set P ( t )to be a piecewise linear function, in the optimization period t ∈ [ t i , t f ] we let P ( t ) be an arbitrary function of t , forthe purpose of computing the optimal control solution.The optimal control solutions that we obtain for differentMSAs are shown in Figs. 3. These solutions are robust to parameter variations (such as different values of thecoefficient c q , see the Supplementary Note 2) and also qualitatively consistent for different MSA’s. Robustnessis also found with respect to the choice of the constraint I max , see the Supplementary Note 2. Typically we firstobserve a quick drop in P ∗ ( t ) (initial tightening), followed by a long nearly constant trend at P ∗ s (steady socialdistancing), and by another drop near the final time (final stranglehold). We remark that P ∗ s corresponds tothe level of reduction needed to suppress the time-varying reproduction number R t ≈
1, and the cost per dayassociated with this level of reduction is (1 − P ∗ s ) /P ∗ s . It follows that for each MSA, there is an almost constantvalue of P ∗ s which well describes the optimal solution, except for the initial time and the final time. The valueof P ∗ s appears to be approximately the same for LA, Houston, and Seattle, while NYC has typically a slighterhigher value of P ∗ s . This also implies that the optimal cost function J ∗ is lower for NYC than for the other cities.In all the four metropolitan areas, P ∗ s is lower than the ‘current’ value of P ( t ) estimated from data, as can beseen from the initial dip in P ∗ ( t ). NYC has the smallest initial drop, indicating that the control action at thetime at which we performed the inference is the closest to optimal, followed by (in the order) Seattle, LA, andHouston. For each metropolitan area, I max was estimated from available data about ICU beds for different USstates, as shown in Table 3 in the Methods.The most important parameter of the optimal control problem is the terminal suppression constraint (cid:15) , whichdescribes the level at which one is trying to suppress the epidemic. The lower is (cid:15) , the closer the goal is toeradication of the disease. The optimal value attained by the objective function J ∗ versus (cid:15) is shown in Figs.4(A) and (B). It should be noted that for each value of (cid:15) , J ∗ shown in Fig. 4 is the lowest possible cost that can beattained. For all the US cities considered, this lowest possible cost increases dramatically as (cid:15) is reduced, which11 . . . A P ∗ ( t ) NYC( I max = 0 . . . . B P ∗ ( t ) LA( I max = 0 . . . . C P ∗ ( t ) Houston( I max = 0 . . . . D P ∗ ( t ) Seattle( I max = 0 . . . . . E S ( t ) R ( t ) 0 . . . . F . . . . G . . . . H - -
20 08 - -
20 09 - -
20 10 - - . . . . I max Q ( t ) Itp + Is - -
20 08 - -
20 09 - -
20 10 - - . . . . - -
20 08 - -
20 09 - -
20 10 - - . . . . - -
20 08 - -
20 09 - -
20 10 - - . . . . Figure 3:
Optimal control solutions for different Metropolitan Statistical Areas within the US. (A-D)Optimal control inputs for the metropolitan cities of NYC, LA, Houston, Seattle (from left to right). (E-H) Timeevolution of the states subject to the optimal control inputs. I max are chosen as the minimum of the range inTable 3 ( ρ = 2 /
3) and c p = c q = 1. The legends in (F-H) are same as the legend in (E).12xemplifies the dilemma between saving human lives and saving the economy. Again, here we are assumingthat the measures of social distancing are optimally implemented, while the cost would be higher in case ofnon-optimal implementation (discussed in the next section). For large values of the suppression constraint (cid:15) ,in all four cities considered, a different type of solution characterized by low cost emerges, which is discussedin more detail in what follows. Qualitatively similar results were obtained when the control input was chosenthat minimizes the alternative objective function (7). A complete study of the effects of varying the suppressionconstraint parameter (cid:15) can be found in the Supplementary Note 3. There are deeper implications behind Figs.4 (A) and (B), namely setting a larger value of (cid:15) corresponds to delaying the economic cost in time, rather thanremoving it.From Fig. 4 (C) and (D) we also see the effects of changing the control horizon ( t f − t i ) of the optimalcontrol problem, from a minimum of 60 days to a maximum of 120 days. We see here that different cities behavedifferently. For Houston and LA we see that a longer t f corresponds to a lower value of the optimal cost J ∗ ,while for NYC the cost increases with t f . A complete study of the effects of varying t f can be found in theSupplementary Note 4. In particular, we see that increasing t f has two contrasting effects on the objectivefunction. On the one hand the cost is integrated over a longer time period, on the other hand the longer timeperiod can be exploited to allow for less stringent measures of social distancing at any time (larger values of P ∗ ),which may result in a lower value for the objective function overall. Thus it appears that finely tuning the timehorizon of the objective function may be used to critically and selectively affect different cities. The implicationsare particularly significant for those cities, Houston and LA, that seem to need a longer time period to suppressthe epidemics. For Seattle we see that J ∗ is minimized for intermediate values of t f in the interval [255 , J or J alt .) The reason for this is thatwhen the integrand in (4) is nonlinear the optimal solution maintains P ∗ ( t ) as close as possible to the linearregime (high values of P ∗ ), which is why we do not see much difference with the linear case (7).In general, the main constraints of the problem are provided by I max and (cid:15) . However, for each case considered,typically either one of these two constraints is dominant. In all the simulations shown in Figs. 3, the dominantconstraint is provided by (cid:15) , with I s ( t ) + I tp ( t ) remaining well below I max over the entire period [ t i , t f ]. Thesesolutions, which we will refer to as the type-1 solution, are characterized by strong measures of social distancing(low P ∗ ) throughout the whole time interval considered, and the dominant constraint is to achieve suppressionof the epidemics at the prescribed terminal time ( t f ).We have also seen the emergence of different solutions, which we refer to as the type-2 solution, when thedominant constraint is given by I max . In these solutions there is at least one time t at which the constraint (6) issatisfied with the equal sign. Overall, type 2 solutions are less expensive economic-wise than type 1 solutions, i.e.,the value of J ∗ is lower. There are also cases when the optimal solution is actively affected by both constraints.In order to better understand the transition between type 1 and type 2 solutions, we have investigated theoptimal control problem (4), for the cases of LA, NYC, and Seattle, as both the suppression parameter (cid:15) andthe control horizon ( t f − t i ) are varied. The transition is characterized by a gradual change in J ∗ , high for type13 − − − A (cid:15) J ∗ NYCLAHoustonSeattle 10 − − − B (cid:15) J a l t ∗ NYCLAHoustonSeattle60 80 100 120100300500700 C t f − t i J ∗ NYCLAHoustonSeattle 60 80 100 120300500700 D t f − t i J a l t ∗ NYCLAHoustonSeattle
Figure 4:
Effects of varying the terminal suppression constraint (cid:15) and the control horizon t f − t i . (A) The optimal cost J ∗ obtained as the solution of the optimization problem (4) vs. the parameter (cid:15) . (B)The optimal cost J alt ∗ obtained as the solution of the optimization problem (7) vs. the parameter (cid:15) . (C) Theoptimal cost J ∗ obtained as the solution of the optimization problem vs. the control horizon ( t f − t i ). (D) Theoptimal cost J alt ∗ obtained as the solution of the optimization problem (7) vs. the control horizon ( t f − t i ). Theparameter c q and c p are both set to 1. 14 +++A − − − (cid:15) ( t f − t i ) J ∗ +++B − − − (cid:15) ( t f − t i ) J ∗ ++ +C − − − (cid:15) ( t f − t i ) J ∗ Figure 5:
The optimal cost J ∗ in the ( t f − t i ) , (cid:15) plane . (A) shows the case of the Los Angeles MetropolitanStatistical Area. (B) shows the case of the NYC Metropolitan Statistical Area. (C) shows the case of the SeattleMetropolitan Statistical Area. For each city, I max is chosen as the maximum of the range in Table 3. Type 1solutions (in green) are more expensive than type 2 solutions (in blue.) The regions in yellow/red correspond tothe transition between the two types of solutions. The parameter c p and c q are both set to 1.1 solutions (in green) and low for type 2 solutions (in blue), with a transition area in between shown in yellowand red. This is illustrated in Fig. 5 and in more detail in Fig. 19 (LA), Fig. 20 (Seattle) and Fig. 21 (NYC) ofthe Supplementary Note 5. Our results show that considerations about the timescale of the control action applydifferently to different cities. As can be seen, the particular emergence of type-1 or type-2 solutions is affectedby both the choice of (cid:15) and ( t f − t i ). From Fig. 5 (A) for LA we see that the most expensive control solutions areachieved when one is trying to suppress the epidemics to a low level in a short time (small (cid:15) and short t f − t i .)This is opposite to what seen for NYC in Fig. 5 (B), where longer time horizons are usually associated withmore expensive solutions. The most convenient solution arises when the suppression constant (cid:15) is large but thecontrol horizon is short (area shown in blue.) Finally, Fig. 5 (C) for Seattle shows that in the case of small (cid:15) ,expensive solutions are obtained when the control horizon is either short or long, while an intermediate controlhorizon is to be preferred. We also see that setting the suppression constraint to a larger value and the controlhorizon to be short can reduce the cost considerably (which is similar to NYC). In general, these results pointout the importance of carefully choosing the timescale over which one is seeking to suppress the epidemics, aswell as the suppression level. Implementation of non-optimal control solutions
We have provided an approach to robustly minimize the effects on the economy of social distancing measuresin the presence of relevant constraints. However, it is possible that a number of considerations may limit the15mplementation of such optimized interventions. We are interested in the effects of implementation of non-optimalcontrols solutions. We thus consider application of a variation of the optimal solution˜ P ( t ) = min { , (1 + α ) P ∗ ( t ) } , (10) α >
0, and analyze violations of the constraints I max and (cid:15) as α is varied. Increasing values of α indicatestronger deviations of the control action from the optimal one. Table 2 summarizes application of such non-optimal interventions in all four cities, with α varying from 0 . .
5. The letter ‘S’ stands for constraintsatisfied, otherwise we report the percentage by which the constraint is violated. We see that the constraint on I max remains always satisfied (this is expected as this constraint is not dominant) but strong violations of thesuppression constraint (cid:15) are otherwise recorded. From Table 2 we see that for α = 0 . (cid:15) inthe four areas of interest is violated by an amount that varies from roughly 4000% for Houston to 80000% forNYC (which corresponds to a fraction of infected people at the final time in the order of 10 − .) These resultsare opposite to those observed previously. The city that is less resilient to variation in the control input is NYC,which was previously reported to have received closer to optimal control interventions. This is due to the higher β for NYC (see Table 1) and to the fact that the optimal P ∗ ( t ) for NYC is higher than for the other cities(corresponding to less strict social distancing required); as a result, increasing P ( t ) further leads to the poorestoutcome. This also highlights the risk for resurgences of the epidemics after partial suppression has been achieved[32], which is currently seen in different parts of the world.16able 2: Violation in the constraints I max and (cid:15) for the non-optimal solutions. NYC: I max = 0 . , c q = 1 , , α I max (cid:15) J . . . . . I max = 0 . , c q = 1 , , α I max (cid:15) J . . . . . I max = 0 . , c q = 1 , , α I max (cid:15) J . . . . . I max = 0 . , c q = 1 , , α I max (cid:15) J . . . . . Discussion
In this paper, we have proposed an approach to optimize social distancing measures in time in order to contain thespread of the COVID-19 epidemics, while minimizing the impact on the economy. Our analysis has been applied tofour different metropolitan statistical areas within the US, but can be directly applied to other geographical areas.Each MSA was shown to be well described by a stylized mathematical model whose parameters were inferred bydaily new case counts reported by local health administrations. Other works have applied optimal control to theCOVID-19 epidemics but none has considered data-driven models. Our approach based on modeling, inferringmodel parameters from real data, and computing optimal control solutions for the inferred model is general andmay be applied to other complex systems of interest.Our choice of the objective function is quite simplified, namely we assume that increasing levels of socialdistancing and quarantining result in progressively higher economical costs. In addition, we did not accountfor the cost associated with medical treatments, which is important yet significantly smaller than the cost ofshutting down the economical activities of the entire society. Our approach can be easily generalized to othermore complex, more realistic types of objective functions, see e.g., [31, 30]. Different from this study, theseobjective functions may also be specific to given regions or try to capture a particular socio-economical model ofinterest.We have found that the optimal control solutions are quite robust to the specific choice of the objective17unction and of its parameters, such as c q . These control solutions tend to be qualitatively similar for differentcities. However, they are affected by the choice of the constraints, in particular the constraint that we haveassociated to suppression of the epidemic, and by the time horizon of the control action. When these are varied,different cities behave differently, which points out the importance of our data-driven approach. We have alsoseen that small deviations from the optimal solution can lead to dramatic violations of the constraints.It is possible to translate these optimal interventions in actual measures that can be imposed on the population,such as restricting the access to certain businesses or venues. While implementation of a time-varying controlmay be challenging in practice, we found that the optimal solution is typically characterized by an initial drop(due essentially to the non-optimality of current control interventions), followed by a nearly constant control(specific to each city) for a long time, and by a final drop, which is needed to achieve the desired suppression ofthe epidemics at the final time. Thus the optimal control solution is almost constant except for the initial timeand the final time, which substantially increases the applicability of this study. The constant part of the solutioncould be practically achieved by dynamical regulation to control the time-varying reproduction number R t ≈ (cid:15) grows or the control horizon ( t f − t i ) grows. A key observation is that the optimal solution P ∗ is generally lowerthan the value of P inferred by data. The initial drop in P ∗ may provide a measure of non-optimality of currentinterventions. This drop was seen to be smallest for NYC compared to other US cities, but was present in allthe cities we have analyzed.Our analysis reveals that the cost of eradicating the disease—i.e., suppressing the number of infected indi-viduals down to a certain critical threshold, such as (cid:28) I max constraint appear to be quitelong (see Methods) and finally comparable with the timescale over which we expect availability of a vaccine. Bysetting a very large control horizon in our simulations, we have seen the emergence of these ‘herd immunitysolutions’, shown in the Supplementary Note 6. However, the parameters of our realistic scenarios, estimatedfrom real data, are far away from the point at which such solutions are optimal. In general, our study allowsto quantify how ‘far’ one is in the parameter space from the point at which a herd immunity solution becomes18ptimal.The time-scale of control interventions, i.e., the control horizon, appears to play a fundamental role. Con-siderations about the control horizon may vary from area to area and may be affected by a number of factors,including expectations about the time at which a vaccine will become available. It appears that cities that haveseen an increase of cases during the inference period need a longer control horizon to suppress the epidemicsoptimally. In certain instances, the impact on the economy can be minimized by tuning the control horizon; forexample, for the city of Seattle we found that an optimal control horizon was equal to roughly 90 days whenthe suppression constraint (cid:15) = 10 − . We wish to emphasize that given the very large economical impact ofsocial distancing measures, even a small improvement in the control strategy can lead to considerable economicalbenefits.One relevant question is whether policy-makers can assess whether a currently employed control action isoptimal or not. The so-called HAMVET procedure, initially proposed in [39] and presented in the SupplementaryNote 8, can be used to validate a control strategy and evaluate its optimality. References [1] Safiya Richardson, Jamie S. Hirsch, Mangala Narasimhan, James M. Crawford, Thomas McGinn, Karina W.Davidson, , and the Northwell COVID-19 Research Consortium. Presenting Characteristics, Comorbidities,and Outcomes Among 5700 Patients Hospitalized With COVID-19 in the New York City Area.
JAMA ,323(20):2052–2059, 05 2020.[2] Lisa Mataloni, Dave Wasshausen, Erich Strassner, Jeannine Aversa. 2020. ”Gross Domestic Product, SecondQuarter 2020 (Advance Estimate) and Annual Update”. U.S. Bureau of Economic Analysis. . August 2020.[3] Scott Gottlieb, Caitlin Rivers, Mark B McClellan, Lauren Silvis, and Crystal Watson. National coronavirusresponse: a road map to reopening.
AEI Paper & Studies , 2020.[4] Seth Flaxman, Swapnil Mishra, Axel Gandy, H Unwin, Helen Coupland, T Mellan, Harisson Zhu, TresniaBerah, J Eaton, P Perez Guzman, et al. Report 13: Estimating the number of infections and the impact ofnon-pharmaceutical interventions on covid-19 in 11 european countries. 2020.[5] Steven Sanche, Yen Ting Lin, Chonggang Xu, Ethan Romero-Severson, Nicolas W Hengartner, and RuianKe. The novel coronavirus, 2019-ncov, is highly contagious and more infectious than initially estimated. arXiv preprint arXiv:2002.03268 , 2020.[6] Sheryl L Chang, Nathan Harding, Cameron Zachreson, Oliver M Cliff, and Mikhail Prokopenko. Modellingtransmission and control of the covid-19 pandemic in australia. arXiv preprint arXiv:2003.10218 , 2020.197] Sean C Anderson, Andrew M Edwards, Madi Yerlanov, Nicola Mulberry, Jessica Stockdale, Sarafa A Iyani-wura, Rebeca C Falcao, Michael C Otterstatter, Michael A Irvine, Naveed Z Janjua, et al. Estimating theimpact of covid-19 control measures using a bayesian model of physical distancing. medRxiv , 2020.[8] Dylan H Morris, Fernando W Rossine, Joshua B Plotkin, and Simon A Levin. Optimal, near-optimal, androbust epidemic control. arXiv preprint arXiv:2004.02209 , 2020.[9] Johannes K¨ohler, Lukas Schwenkel, Anna Koch, Julian Berberich, Patricia Pauli, and Frank Allg¨ower.Robust and optimal predictive control of the covid-19 outbreak. arXiv preprint arXiv:2005.03580 , 2020.[10] IHME COVID, Christopher JL Murray, et al. Forecasting covid-19 impact on hospital bed-days, icu-days,ventilator-days and deaths by us state in the next 4 months.
MedRxiv , 2020.[11] Daron Acemoglu, Victor Chernozhukov, Iv´an Werning, and Michael D Whinston. A multi-risk sir modelwith optimally targeted lockdown. Technical report, National Bureau of Economic Research, 2020.[12] Shawn Radcliffe. 2020. ”Here’s Exactly Where we Are with Vaccines and Treatments for COVID-19”. Healthline. . August 2020.[13] Haley E Randolph and Luis B Barreiro. Herd immunity: Understanding covid-19.
Immunity , 52(5):737–741,2020.[14] Yen Ting Lin, Jacob Neumann, Ely F Miller, Richard G Posner, Abhishek Mallela, Cosmin Stafa, JaideepRay, Gautam Thakur, Supriya Chinthavali, and William S Hlavacek. Daily Forecasting of New Cases forRegional Epidemics of Coronavirus Disease 2019 with Bayesian Uncertainty Quantification. medRxiv , 2020.[15] Maria Nicola, Zaid Alsafi, Catrin Sohrabi, Ahmed Kerwan, Ahmed Al-Jabir, Christos Iosifidis, MalihaAgha, and Riaz Agha. The socio-economic implications of the coronavirus pandemic (covid-19): A review.
International journal of surgery (London, England) , 78:185, 2020.[16] Anu G Gupta, Cheryl A Moyer, and David T Stern. The economic impact of quarantine: Sars in torontoas a case study.
Journal of Infection , 50(5):386–393, 2005.[17] Mikl´os Koren and Rita Pet˝o. Business disruptions from social distancing. arXiv preprint arXiv:2003.13983 ,2020.[18] Lother Wieler, Ute Rexroth, Rene Gottschalk. 2020. ”Emerging COVID-19 success story: Germany’s strongenabling environment”. Our World in data. https://ourworldindata.org/covid-exemplar-germany .August 2020. 2019] France 24. 2020. ”Germany reimposes local lockdowns after regional coronavirus outbreak”.France 24. . August 2020.[20] New York Times repository of Covid-19 data in the United States.[21] U.S. Office of Management and Budget Delineation Files.[22] Stephen A Lauer, Kyra H Grantz, Qifang Bi, Forrest K Jones, Qulu Zheng, Hannah R Meredith, Andrew SAzman, Nicholas G Reich, and Justin Lessler. The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application.
Annals of internal medicine ,2019, 2020.[23] Hiroshi Nishiura, Natalie M Linton, and Andrei R Akhmetzhanov. Serial interval of novel coronavirus(covid-19) infections.
International journal of infectious diseases , 2020.[24] Roman W¨olfel, Victor M Corman, Wolfgang Guggemos, Michael Seilmaier, Sabine Zange, Marcel A M¨uller,Daniela Niemeyer, Terry C Jones, Patrick Vollmar, Camilla Rothe, and Others. Virological assessment ofhospitalized patients with COVID-2019.
Nature , 581(7809):465–469, 2020.[25] Aki Sakurai, Toshiharu Sasaki, Shigeo Kato, Masamichi Hayashi, Sei-ichiro Tsuzuki, Takuma Ishihara,Mitsunaga Iwata, Zenichi Morise, and Yohei Doi. Natural History of Asymptomatic SARS-CoV-2 Infection.
New England Journal of Medicine , 0(0):null, 2020.[26] Nguyen Van Vinh Chau, Vo Thanh Lam, Nguyen Thanh Dung, Lam Minh Yen, Ngo Ngoc QuangMinh, Le Manh Hung, Nghiem My Ngoc, Nguyen Tri Dung, Dinh Nguyen Huy Man, Lam Anh Nguyet,Le Thanh Hoang Nhat, Le Nguyen Truc Nhu, Nguyen Thi Han Ny, Nguyen Thi Thu Hong, Evelyne Kestelyn,Nguyen Thi Phuong Dung, Tran Chanh Xuan, Tran Tinh Hien, Nguyen Thanh Phong, Tran Nguyen HoangTu, Ronald B Geskus, Tran Tan Thanh, Nguyen Thanh Truong, Nguyen Tan Binh, Tang Chi Thuong, GuyThwaites, Le Van Tan, and Oxford University Clinical Research Unit COVID-19 Research Group. The Nat-ural History and Transmission Potential of Asymptomatic Severe Acute Respiratory Syndrome Coronavirus2 Infection.
Clinical Infectious Diseases , 2020.[27] Ministry of Health, Official report on the cruise ship Diamond Princess, May 1, 2020.[28] LANL COVID-19 Prediction GitHub.[29] Christophe Andrieu and Johannes Thoms. A tutorial on adaptive mcmc.
Statistics and computing , 18(4):343–373, 2008.[30] Afroza Shirin, Feng Song, Yen-Ting Lin, William S Hlavacek, and Sorrentino S. Prediction of optimal drugschedules for controlling autophagy.
Scientific Reports , 9(1428), 2019.2131] Afroza Shirin, Fabio Della Rossa, Isaac Klickstein, John Russell, and Francesco Sorrentino. Optimal reg-ulation of blood glucose level in type i diabetes using insulin and glucagon.
PloS one , 14(3):e0213665,2019.[32] Nick Warren Ruktanonchai, JR Floyd, Shengjie Lai, Corrine Warren Ruktanonchai, Adam Sadilek, PedroRente-Lourenco, Xue Ben, Alessandra Carioli, Joshua Gwinn, JE Steele, et al. Assessing the impact ofcoordinated covid-19 exit strategies across europe.
Science , 2020.[33] Allan Saul, Nick Scott, Brendan S Crabb, Suman S Majundar, Benjamin Coghlan, and Margaret E Hellard.Victoria’s response to a resurgence of covid-19 has averted 9,000-37,000 cases in july 2020.
The MedicalJournal of Australia , page 1, 2020.[34] Lorenzo Tondo. Italy at crossroads as fears grow of covid-19 second wave.
The Guardian , 2020.[35] Ashifa Kassam. Spain warned of dire impact of second coronavirus lockdown.
The Guardian , 2020.[36] Jasmine C. Lee, Sarah Mervosh, Yuriria Avila, Barbara Harvey, and Alex Leeds Matthews. Spain warnedof dire impact of second coronavirus lockdown.
The New York Times , 2020.[37] Eric JW Orlowski and David JA Goldsmith. Four months into the covid-19 pandemic, sweden’s prized herdimmunity is nowhere in sight.
Journal of the Royal Society of Medicine , 113(8):292–298, 2020.[38] Kin On Kwok, Florence Lai, Wan In Wei, Samuel Yeung Shan Wong, and Julian WT Tang. Herd immunity–estimating the level required to halt the covid-19 epidemics in affected countries.
Journal of Infection ,80(6):e32–e33, 2020.[39] I Michael Ross.
A primer on Pontryagin’s principle in optimal control . Collegiate publishers, 2015.[40] Javier Perez-Saez, Stephen A Lauer, Laurent Kaiser, Simon Regard, Elisabeth Delaporte, Idris Guessous,Silvia Stringhini, Andrew S Azman, Serocov-POP Study Group, et al. Serology-informed estimates ofsars-cov-2 infection fatality risk in geneva, switzerland. medRxiv , 2020.[41] Pipetius Quah, Andrew Li, and Jason Phua. Mortality rates of patients with covid-19 in the intensive careunit: a systematic review of the emerging literature.
Critical Care , 24:1–4, 2020.[42] Janet Christenbury. 2020. ”Study in ICU Finds 30.9% Mortality Rate From COVID-19”. Futurity. . August 2020.[43] Xiaofan Liu, Hong Zhou, Yilu Zhou, Xiaojun Wu, Yang Zhao, Yang Lu, Weijun Tan, Mingli Yuan, XuhongDing, Jinjing Zou, et al. Risk factors associated with disease severity and length of hospital stay in covid-19patients.
Journal of Infection , 81(1):e95–e97, 2020.22 ethods
Estimation of I max The I max values can be approximated with simple assumptions as seen in Table 3 for a selection of major U.S.cities. These I max values are a ratio of ICU beds to the population of people which require them. Following [14],we estimate that the probability of death conditioned on symptomatic infection is equal to f H × (1 − f R ) = 0 . f H and f R were independently computed in [40] and [1], respectively. Data [41, 42]shows that the mortality rate for patients sent to ICU is between 30% and 40%, thus, it is reasonable to assumethat an overall fraction of infected people equal to 0 . / .
35 = 0 . .