key: cord-269363-drjj705k
authors: Nenchev, Vladislav
title: Optimal quarantine control of an infectious outbreak
date: 2020-07-28
journal: Chaos Solitons Fractals
DOI: 10.1016/j.chaos.2020.110139
sha: 
doc_id: 269363
cord_uid: drjj705k

This paper studies the optimal control of an infectious spread based on common epidemic models with permanent immunity and no vaccine availability. Assuming limited isolation control and capacity constraints on the number of infections, an optimal quarantine control strategy that balances between the total number of infections and the overall isolation effort is derived from necessary optimality conditions. The specific optimal policy is then obtained by optimizing the switching times of this generalized strategy. In the case of a newly emerged disease, these results can be used in a data-driven receding horizon manner to improve actions as more data becomes available. The proposed approach is applied to publicly available data from the outbreak of SARS-CoV-2 in Germany. In particular, for minimizing the total number of infections or the number of isolated individuals, the simulations indicate that a sufficiently delayed and controlled release of the lock-down are optimal for overcoming the outbreak. The approach can support public health authorities to plan quarantine control policies.

Since the first recorded COVID-19 infections in Wuhan, China, the disease has spread worldwide. There has been a plurality of attempts to interpret and predict the development of the disease based on public data using first-principles epidemiological [1] [2] [3] , data-driven [4] , or mixed models [5, 6] . Based on these predictions, a wide variety of temporal isolation policies have been established in most affected countries. However, evaluating the effectiveness of these responses is notoriously difficult due to the reliability of available data, and since their impact on the outbreak evolution can only be measured with a substantial delay.

Mathematical analysis of first-principle models of biological systems can provide valuable insight, e.g., on their general properties [25] , or how to control infectious disease outbreaks [7] . In the simplest model, the population is divided into the susceptible S , infected I and recovered R groups (therefore called the SIR model), and their evolution over time is represented as a set of coupled ordinary differential equations. Augmenting the state by exposed E individuals, leads to the so-called SEIR model. Applications of outbreak controls minimizing various different cost functions have been proposed for the SIR model, e.g., [8] [9] [10] [11] [12] to name a few, as well as heuristically tailored policies [13, 14] . Although S(E)IR models can be too simplistic for many real outbreaks, they often allow deriving complete analytical solutions, thus providing a foundation E-mail address: vladislav.nenchev@gmail.com for rigorous mathematical results upon which more complex models can be built.

An issue of practical concern for many disease outbreaks without an available vaccine, such as for SARS-CoV-2 as of June 2020, is minimizing the overall quarantine effort or the final outbreak size, while respecting control and capacity constraints on the current number of infections. To account for quarantine measure (similar to [5] and [15] ), the SIR model is extended by a quarantined subpopulation Q . Infected individuals are moved depending on the applied control rate from I to Q , where they cannot cause additional infections. Then, using the direct adjoining approach [16] and applying Pontryagin's Minimum Principle (PMP) to the corresponding state-and input-constrained control problem, in this paper it is shown that bang-bang or bang-boundary controls are admitted in a particular sequence to solve the corresponding optimal control problem. To obtain the specific optimal policy, the switching times between the bang-boundary subarcs are determined in an induced optimization problem. Analogous results can be obtained for the SEIR model. Upon an outbreak of a previously unknown disease, better model parameter estimates can be obtained as more data becomes available, and the induced optimization problem can be recomputed in a data-driven receding horizon manner to improve actions. The proposed methods are applied to publicly available data until June 2020 from the outbreak of SARS-CoV-2 in Germany. In particular, minimizing the total outbreak size or the total number of isolated individuals are studied. The findings of this work can be used to evaluate if the policies that have been established in practice resemble the optimal policy, and support public health authorities in implementing effective future quarantine strategies.

The remainder of this paper is structured as follows. In Section 2 , a mathematical statement of the problem is provided. In Section 3 , the optimal quarantine policy under isolation control and capacity constraints, and, in Section 4 , the data-driven receding horizon control scheme are described. Then, in Section 5 , the approach is applied based on COVID-19 data from Germany until June 2020. Finally, the main results and future work are discussed in Section 6 .

We address epidemics with no vaccine availability, in which the only available control is based on isolation. Let the entire population p be divided into three sub-populations: susceptible S ; infected I ; and recovered R . Births' and deaths' effects are neglected. The recovered population is to be interpreted as those who can no longer infect others, and, who cannot become reinfected by the disease. This is known as the classic epidemiological SIR model [7] . To include quarantine or isolation control, the model is augmented by introducing a time varying quarantine strength control u ( t ) and a corresponding quarantined (and also recovered from the disease) population Q . Then, with x = [ S, I, R, Q] T ⊆ R 4 + the sub-populations' evolution is governed by the following system of coupled nonlinear ordinary differential equations:

where the population p = 4

i =1 x i (t) , and the infection b and recovery m rates are assumed to be constant in time. Note that the model assumes homogeneous mixing among the sub-populations, as well as uniform susceptibility and disease progress for every individual in the population. The SIR model assumes a direct transition from susceptible to infected. It can be easily augmented by including exposed individuals E to account for a significant incubation period during which individuals have been infected but are not yet infectious themselves. This so-called SEIR epidemiological model [7] has been employed in several prior studies, such as the SARS outbreak, e.g., [17, 18] , as well as the COVID-19 outbreak [1] [2] [3] . Since the SEIR model did not provide any additional qualitative insights in the studied problem, for simplicity, the SIR model is analyzed, and details for the SEIR model can be found in the appendix. An important characteristic of epidemics is the reproduction number

Thus, applying the control u ( t ) results in reducing the effective reproduction number R t . Determining the maximal quarantine control u max is highly dependent on local quarantine policies and population specifics, directly corresponding to the minimally achievable reproduction number (2) . Another important parameter of epidemics is the fraction α s of severe conditions among infected individuals requiring hospitalization and admission to an intensive care unit. Highly infectious diseases (with a high R t ) causing many severe conditions in the infected population (corresponding to a high α s ), may quickly lead to an overload of available treatment resources and, thus, cause a large number of additional potentially avoidable deaths. Assuming that the overall available capacity is described by α c , this capacity limitation can be described as the

In this work, the goal is to obtain an optimal quarantine control policy u ( t ), t ∈ [0, t f ] for a fixed final time t f , that minimizes a weighted combination of the total number of infections and the overall number of quarantined individuals at time t f . Note that the overall number of quarantined individuals corresponds to the applied level of containment and can, thus, be linked to the economical impact of the outbreak. The total number of infections is denoted by the sum of the final recovered and final quarantined populations, and is directly correlated to the overall number of deaths. This can be captured by the cost function

where w ∈ [0, 1]. The control is limited to u ( t ) ∈ [0, u max ] at all times, such that u = 0 corresponds to an uncontrolled epidemics evolution, and u = u max to the case, when maximal possible isola-

t f ] exists, this leads to the following (Mayer type) optimal control problem. (4) is minimized subject to the dynamics (1) , while respecting the constraint (3) at all times.

In the following, using the direct adjoining approach and applying Pontryagin's Minimum Principle, necessary optimality conditions are derived for the optimal control problem above. Further, it is shown that the obtained locally optimal solutions constitute an optimal generalized strategy under certain assumptions.

As (3) is a pure state constraint, (higher order) time derivatives are used to obtain the control u ( t ) that yields a system evolution along the active state constraint. Since u ( t ) appears linearly in the state equation of x 2 (1) , the first total time derivative of the function h ( x ( t )) (3) contains the control explicitly:

. If τ 1 and τ 2 are maximal, let τ 1 be the entry-time and τ 2 the exit-time of the boundary arc. On the boundary arc,

Thus, the control on a boundary arc is obtained from

The boundary control (6) depends on x 1 and is, therefore, not guaranteed to be in the feasible control set [0 , u max ] . Thus, assume that

To derive first order necessary optimality conditions by applying PMP, a Lagrange multiplier η(t) ∈ R associated with the state constraint (3) is introduced. For that, based on the aforementioned regularity and boundedness conditions on a boundary arc, and supposing that the state space constraint is not active at the initial and final time, i.e., h ( x (0)) < 0, h ( x ( t f )) < 0, the direct adjoining approach [16] is applied. The corresponding augmented Hamiltonian is given by

with the adjoint variable λ ∈ R 4 . As shown in Maurer et al. [19] , there exist an absolutely continuous λ( t ) and a piecewise absolutely continuous multiplier function η( t ), such that the following

For the SIR model, the following holds:

As the Hamiltonian (7) is linear in u ( t ), the optimal control will depend on the sign of the switching function, given by

As suggested in Ledzewicz and Schättler [20] , with the total time derivative of (11) and the costates (8) ,

holds on a boundary arc with t ∈ [ τ 1 , τ 2 ]. This yields the explicit expression for the multiplier

Assuming that a solution of Problem 1 exists and the control u ( t ) has finitely many switching times, the following proposition holds.

The generalized optimal control policy for Problem 1 is

If t 3 > t 2 , the tangency constraints (16) must also hold.

Proof. On interior arcs, where the state constraint is inactive h ( x ( t )) < 0, from the minimum condition (10) , the optimal control is

The switching function (11) is zero ⇐⇒ λ 2 (t) = 1 or x 2 (t) = 0 .

The latter case is trivial and will be neglected in the following. To fulfill λ 2 (t f ) = 0 , λ 2 (0) > 1 and ˙ λ 2 (0+) ≤ 0 , such that the optimal control policy may start with a zero segment for [ t 0 , t 1 ). If

) > 0 and σ (t 1 −) < 0 assuming u (t 1 −) = 0 , and as −˙ σ (t 1 )(0 − u max ) > 0 , a second segment for t ∈ [ t 1 , t 2 ) may have the control u max . On a boundary arc, given the regularity and boundedness conditions for the boundary control, the minimum condition

) ≤ u max (from the boundedness condition of the boundary control) and ˙ x 2 (t 2 +)( ū b ) = 0 , and thus, the continuity conditions on the boundary σ (

Substituting relevant equations into (13) yields

, resulting from the boundedness condition of the boundary control, completes the additional tangency constraints (16) .

If h ( x ( t )) < 0 and σ (t) = 0 for t ∈ [ ˜ τ 1 , ˜ τ 2 ) , a singular arc may occur. Since ∂ σ /∂ u = 0 holds trivially, the singular control u s can be obtained from d σ /d t = 0 with η(t) = 0 and λ 2 (t) = 1 , which yields u s = b p x 1 − m = u b . The optimality of the corresponding singular arc can be checked with the generalized Legendre-Clebsch condition of first order, which yields − ∂ ∂u

Transitioning from an interior arc to a singular arc and vice versa at time τ ∈ [0, t f ] violates the continuity conditions 

If the condition above holds for any junction time τ , then μ( t 2 , τ ) > 0 for the boundary arc and μ( τ , τ ) ≤ 0 for the singular arc would have to hold simultaneously. Thus, a singular arc will not be part of the optimal solution.

Finally, since h (x (t 3 −)) = 0 and h 1 (x (t 3 −)) = 0 if t 2 = t 3 , ˙ x 2 (t 3 )(u = 0) < 0 and the optimal control (14) ends with a zero segment for t ∈ [ t 3 , t f ], which also satisfies the transversality condition (9) .

If u b (x 1 (t 2 )) ≤ u max does not hold, by analogous derivations to the proof above, it can be shown that adding a finite number n of sequences of a zero arc followed by an u max -arc, until a switch to the boundary arc at time t 2 n +2 becomes possible with u b (x 1 (t 2 n +2 )) = u max , fulfills the necessary optimality conditions.

For simplicity, this case is excluded in the following elaborations.

As a consequence of the above proposition, solutions of Problem 1 will exhibit the structure of the generalized bangboundary control policy (14) . Note that the specific optimal policy strongly depends on the choice of w as well as the model parameters. What remains to be determined are the optimal switching times. Instead of directly optimizing the switching times t 1 , t 2 , t 3 , t 4 , introduce the arc durations ξ i = t i − t i −1 , assuming t 4 = t f . Thus, the induced optimization problem reads

(1) , (3) , (14) , (15) , (16) 

This optimization problem can be simplified further. If w = 1 , i.e., the goal is to minimize only the total number of infections, the op- (3) is violated. For other values of w , the optimal control is given by (14) , which denotes two possible policies -one containing a boundary arc, and one without a boundary arc. Note that (1) is explicitly solved on the boundary arc, as already used in the proof of Proposition 1 . Thus, solving (17) can be replaced by solving the three aforementioned boundary value problems and choosing the policy that yields the minimal cost.

Corresponding derivations for the SEIR model can be found in the appendix. The induced optimization problem (17) can be recomputed in a data-driven manner to obtain an improved quarantine action, as shown in the following.

While a fast response upon an epidemic outbreak is essential, for a newly emerged disease, like in the case of COVID-19 at the beginning of 2020, disease characteristics are highly uncertain. As more data and insights become available over time, it can be expected that the parameter estimates of the models can be improved continuously. This can be achieved, e.g., by minimizing the mean square error loss function between measured state data ˜ 

(1) .

Analogously, this can be done for the SEIR model with P = { b, m, c} . At any time τ , using the parameter estimates P τ obtained by (18) and the measured state ˜

x (τ ) , the induced optimization problem (17) can be re-solved with x 0 = ˜ x (τ ) . Assume that the parameter estimates are updated whenever new data becomes available. Based on that, a data-driven receding horizon approach is adopted, where the control is obtained by re-computing (17) , only when the current parameter estimates deviate from the parameter estimates used in the previous computation step evaluated by a function ζ by a pre-defined threshold δ. The similarity function ζ was assumed to be a simple absolute deviation in this work, but it can contain more complex considerations like, e.g., the covariance of the estimates. This leads to Algorithm 1 . Note that, in general, receding horizon control approaches with uncertain parameters cannot guarantee that hard constraints will not be violated. Thus, to increase satisfaction probability, the capacity constraint (3) can be tightened conservatively based on the confidence of the current parameter estimate P t .

The proposed method is applied to data from the SARS-CoV-2 outbreak in Germany until June 2020. Two cases are considered -w = 0 and w = 0 . 5 , which correspond to minimizing only the overall number of isolated individuals and a trade-off between the overall number of isolated individuals and the total size of the outbreak, respectively. The optimization problems are solved with SciPy's bvp and minimize routines.

Most parameters are taken from a report of the German Robert Koch Institute [21] -a susceptible population of x 1 (0) = 80 . 10 6 and an initial number of infections x 2 (0) = 10 0 0 , corresponding to the number of registered infections in Germany on March 8, 2020. The initial numbers of recovered and quarantined individuals are both 0. For the SEIR model, an initial exposed individuals number of 100 is assumed. Further, a basic reproduction number R t = 2 . 0 is assumed in the uncontrolled models, as well as an inverse mean infectious period in days of m = 0 . 091 . The percentage of infected individuals admitted to an intensive care unit is 1.2% and the mortality is 0.5%. The overall number of intensive care units in Germany was α c = 40 . 10 3 by the end of March 2020 [22] . The assumed maximum control value u max is defined such that it leads to an effective reproduction number ˜ R t = 1 . 0 , i.e., with Eq. (2) given by u max = 0 . 091 .

For the receding horizon control approach, publicly available data from RKI from March 1 until June 1, 2020 to re-estimate the parameters P 1 = { b, c, m } for this time period, and assume the actual parameters R t = 2 . 5 , c = 0 . 25 and u max = 0 . 068 thereafter. The control is re-computed according to Algorithm 1 .

Figs. 1 and 2 show the model state trajectories, the capacity constraint and the control over time for the uncontrolled (left) and the optimally controlled (right) with the aforementioned parameters for the SIR and SEIR model, respectively. It is notable that for both cases, i.e., for w = 0 and w = 0 . 5 , respectively, the optimal solution is given by a bang-boundary control, when using the SIR or the SEIR model.

Applying the optimal control to the SIR model under the given assumptions, the pandemic is predicted to lead to an expected total number of cases of 54 084 981, and a total number of deaths of 270 424. The simulations also indicate that the turnaround (final peak of active cases) should occur 26 weeks after the start of the outbreak. The overall outbreak should be completely overcome after 44 weeks, without a second wave of infections. This closely resembles predictions made in an der Heiden and Buchholz [21] and [23] . Fig. 3 shows the result with the data-driven receding horizon control. Due to the high estimate of the reproductive number in the initial phase, the maximal isolation control starts earlier, and, overall, applies more isolation due to the parameter uncertainty. Note that the control maintains the capacity constraint. In practice, a toggling between bang-bang controls could be prevented by introducing a hysteresis, when a new control policy is computed in Algorithm 1 .

Applying the optimal control under the given assumptions, the pandemic is predicted to lead to an expected total number of cases of 49 004 201, and a total number of deaths of 245 021. Note that the lower number of infections is a consequence of the longer period of applied quarantine control compared to the preceding simulation example. Another consequence of the longer quarantine is a slightly delayed turnaround that should occur 30 weeks after the start of the outbreak. The overall outbreak should be completely overcome after 50 weeks, without a second wave of infections, and while control and capacity constraints are maintained at all times.

The optimal quarantine control of an infectious spread was studied assuming an SIR epidemic model with permanent immunity and no vaccine availability. The control strategy, which min- Parameters : initial parameters P 0 , initial times t 0 , 1 , t 0 , 2 , t 0 , 3 , threshold δ, final time t f Output : optimal switching times t 1 , t 2 , t 3 1: if P t−1 = ∅ then 2: (18) with data x | [0 ,τ ] and initial parameter value P t−1 . 5: if ζ (P t−1 , P t ) > δ then 6: t 1 , t 2 , t 3 ← solve (17) initialized with previous t 1 , t 2 , t 3 7: P t−1 ← P t . 8: end if 9: return t 1 , t 2 , t 3 imizes the final outbreak state assuming that only limited quarantine control is available and capacity constraints have to be respected, was derived from necessary optimality conditions. The solution is to possibly delay the control action first, and then, if needed, apply maximal isolation. If the capacity boundary is reached, a boundary control maintaining the infection numbers may be applied afterwards, until the infection numbers begin to drop. This result was used in a data-driven receding horizon control approach to improve actions as more data becomes available for a previously unknown disease. The application of the proposed schemes to publicly available data from the outbreak of SARS-CoV-2 in Germany indicate that a sufficiently delayed and controlled release of the lockdown are optimal, and would result in overcoming the outbreak within one year.

Note that a constant overall population was assumed throughout this work. Effects resulting from population exchange could be addressed, e.g., as suggested in Xue et al. [24] , leading to a stochas-tic induced optimization problem. An interesting remaining question is how to determine the maximal isolation control value, as well as how to map the boundary control values to practical isolation actions. A promising approach for addressing both topics could be using machine learning techniques, e.g., building upon ideas from [5] . Another research direction could be to embed the parameter estimation uncertainty into the data-driven optimization to achieve an improved exploration-exploitation trade-off.

The author declares that he has no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. and c denote the exposure and infection rates, respectively. The reproduction number remains as defined in (2) . Analogously to the elaborations provided for the SIR model, the most important equations are derived as:

Following an analogous approach, Proposition 1 can be proven for the SEIR model, yielding an identically structured optimal policy.

Novel corona virus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions. medRxiv

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

Prediction of criticality in patients with severe COVID-19 infection using three clinical features: a machine learning-based prognostic model with clinical data in Wuhan. medRxiv

Quantifying the effect of quarantine control in COVID-19 infectious spread using machine learning. medRxiv

Early dynamics of transmission and control of COVID-19: a mathematical modelling study

Infectious diseases of humans: dynamics and control

An optimal isolation policy for an epidemic

Optimal isolation policies for deterministic and stochastic epidemics

An explicit optimal isolation policy for a deterministic epidemic model

Optimal control of epidemics with limited resources

Time-optimal control strategies in sir epidemic models

Strategic release of lockdowns in a covid infection model. medRxiv

Optimal policies for control of the novel coronavirus disease (COVID-19) outbreak

Effective containment explains subexponential growth in recent confirmed COVID-19 cases in China

A survey of the maximum principles for optimal control problems with state constraints

Modelling the SARS epidemic by a lattice-based montecarlo simulation

Extension and verification of the SEIR model on the 2009 influenza A(H1N1) pandemic in Japan

Optimization methods for solving bang-bang control problems with state constraints and the verification of sufficient conditions

A local field of extremals for single input systems with state space constraints

Modellierung von Beispielszenarien der SARS-CoV-2-Epidemie 2020 in Deutschland

Projecting the spread of COVID19 for Germany. medRxiv2020

A data-driven network model for the emerging COVID-19 epidemics in Wuhan, Toronto and Italy

Robustification as a tool in modeling biochemical reaction networks

The author thanks M.C. Lüffe for enriching discussions.

Including quarantine control analogously to (1) , the subpopulationsâ evolution in the SEIR model is governed by the following system of coupled nonlinear ordinary differential equations: