key: cord-0489585-vgq5duh3
authors: Molnar, Tamas G.; Singletary, Andrew W.; Orosz, Gabor; Ames, Aaron D.
title: Safety-Critical Control of Compartmental Epidemiological Models with Measurement Delays
date: 2020-09-22
journal: nan
DOI: nan
sha: 86696ef8089cb93989043fd3aee484ca3a41954a
doc_id: 489585
cord_uid: vgq5duh3

We introduce a methodology to guarantee safety against the spread of infectious diseases by viewing epidemiological models as control systems and by considering human interventions (such as quarantining or social distancing) as control input. We consider a generalized compartmental model that represents the form of the most popular epidemiological models and we design safety-critical controllers that formally guarantee safe evolution with respect to keeping certain populations of interest under prescribed safe limits. Furthermore, we discuss how measurement delays originated from incubation period and testing delays affect safety and how delays can be compensated via predictor feedback. We demonstrate our results by synthesizing active intervention policies that bound the number of infections, hospitalizations and deaths for epidemiological models capturing the spread of COVID-19 in the USA.

The rapid spreading of COVID-19 across the world forced people to change their lives and practice mitigation efforts at a level never seen before, including social distancing, mask-wearing, quarantining and stay-at-home orders. These human actions played a key role in reducing the spreading of the virus, although such interventions often have economic consequences, lose of jobs and physiological effects. Therefore, it is important to focus mitigation efforts and determine when, where and what level of intervention needs to be taken.

This research provides a methodology to determine the level of active human intervention needed to provide safety against the spreading of infection while keeping mitigation efforts minimal. We use compartmental epidemiological models to describe the spreading of the infection [1] , [2] , and we view these models as control systems where human intervention is the control input. Viewing epidemiological models as control systems has been proposed in the literature recently [3] , [4] , [5] , and various models with varying transmission rate [6] , [7] , [8] , [9] have appeared to quantify the level of human interventions in the case of COVID- 19. In this paper, we build on our recent work [10] and use a safety-critical control approach to synthesize control strategies that guide human interventions so that certain safety criteria (such as keeping infection, hospitalization and death below given limits) are fulfilled with minimal mitigation efforts. The approach is based on the framework of Fig. 1 . Illustration of the SIR model as control system and its fit to US COVID-19 data [10] . Model parameters were estimated from compartmental data (right) by accounting for a measurement delay τ . The transmission rate and the corresponding control input (left) were fitted to mobility data. control barrier functions [11] , [12] that leverages the theory of set invariance [13] for dynamical [14] , [15] and control systems [16] , [17] , [18] . We take into account that data about the spreading of the infection may involve significant measurement delays [5] , [19] , [20] , [21] due to the fact that infected individuals may not show symptoms and get tested for quite a few days. We use predictor feedback control [22] , [23] , [24] to compensate these delays, and we provide safety guarantees against errors in delay compensation.

The outline of the paper is as follows. Section II introduces a generalized compartmental model, which covers the class of the most popular epidemiological models. Section III introduces safety critical control without considering measurement delays, while Sec. IV is dedicated to delay compensation. Conclusions are drawn in Sec. V.

Compartmental models describe how the size of certain populations of interest evolve over time. Consider n + m compartments, given by x ∈ R n+m , which are separated into two groups: n so-called multiplicative compartments, given by w ∈ R n , and m outlet compartments, given by z ∈ R m . The evolution of these compartments over time t can be given by the following generalized compartmental model:

where x = [w T z T ] T , initial conditions are x(0) = x 0 , and f, g : R n → R n , q : R n → R m and r : R m → R m depend on the choice of the model; see Examples 1, 2 and 3.

In model (1), the multiplicative compartments w are populations that essentially describe the transmission of the infection. The transmission can be reduced by active interventions, whose intensity is quantified by a control input u ∈ U ⊂ R. The outlet compartments z, on the other hand, do not actively govern the transmission, but rather indicate its effects, as they are driven by the evolution of the multiplicative compartments.

Example 1. SIR model. One of the most fundamental epidemiological models is the SIR model [25] , [26] that consists of susceptible, S, infected, I, and recovered, R, populations. The SIR model captures the spread of the infection based on the interplay between the susceptible and infected populations. Thus, S and I are multiplicative compartments, while R, that measures the number of recovered (or deceased) individuals, is an outlet compartment. The model uses three parameters: the transmission rate β 0 > 0, the recovery rate γ > 0 and the total population N . Active interventions given by the control input u ∈ [0, 1] allow the population to reduce the transmission to an effective rate β = β 0 (1 − u), where u = 0 means no intervention and u = 1 means total isolation of infected individuals. This puts the SIR model with active intervention to form (1) where

Example 2. SEIR model. The SEIR model [27] , [28] is an extension of the SIR model that incorporates an exposed population E apart from the S, I and R compartments. The exposed individuals are infected but not yet infectious over a latency period given by 1/σ > 0. Since the latency affects the transmission, E is a multiplicative compartment. The SEIR model can be described by (1) with

Example 3. SIHRD model. The SIHRD model [10] adds two more outlet compartments to the SIR model: hospitalized population H and deceased population D. Their evolution is captured by three additional parameters: the hospitalization rate λ > 0, the recovery rate ν > 0 in hospitals and the death rate µ > 0. Equation (1) yields the SIHRD model for

There exist several other compartmental models of form (1) which involve further compartments, such as the SIRD [29] , SIRT [7] , SIXRD [30] or SIDARTHE [1] models. More complex models can provide higher fidelity, although they involve more parameters that need to be identified. In what follows, we show applications of the SIR and SIHRD models and we discuss the occurrence of time delays related to incubation and testing. We omit further discussions on latency, the SEIR model or other more complex models. Fig. 1 shows the performance of the SIR model in capturing the spread of COVID-19 for the case of US national data. The parameters β 0 = 0.33 day −1 , γ = 0.2 day −1 and N = 33 × 10 6 of the SIR model were fitted following the algorithm in [10] to the recorded number of confirmed cases I + R [31] between March 25 and August 9, 2020, while the control input u(t), that represents the level of quarantining and social distancing, was identified from mobility data [32] based on the medium time people spent home. The fitted control input (blue) follows the trend of the mobility data (gray) well, especially in March where stay-at-home orders came into action. While the fit (blue) captures the data about confirmed cases (gray), the model also has predictive power (orange); see more details about forecasting in [10] .

Note that once an individual gets infected by COVID-19, it takes a few days of incubation period to show symptoms and an additional few days to get tested for the virus [5] , [19] , [20] , [21] . Therefore, the measured number of confirmed cases represents a delayed state of the system,

, and thus we involved a time delay τ in the model identification process, which was found to be τ = 11 days by fitting [10] . The delay-free counterpart of the fit (purple) shows that the measurement delay can lead to a significant error in identifying the true current level of infection. The effects of the delay τ on safety-critical control and its compensation will be discussed in Sec. IV.

Formally, safety can be translated into keeping system (1) within a safe set S ⊂ R n+m that is the 0-superlevel set of a continuously differentiable function h : R n+m → R:

Function h prescribes the condition for safety: for example, if one intends to keep the infected population I under a limit I max for the SIR, SEIR or SIHRD models, the safety condition is h(x) = I max − I ≥ 0.

To guarantee safety, we design a controller

that ensures that the set S in (5) is forward invariant under the dynamics (1), i.e., if x(0) ∈ S (h(x(0)) ≥ 0), then x(t) ∈ S (h(x(t)) ≥ 0) for all t > 0. Below we use the framework of control barrier functions [11] , [12] to synthesize controllers that are able to keep certain compartments of interest within prescribed limits. First, we consider safety for multiplicative compartments, and then for outlet compartments.

Consider keeping the i-th multiplicative compartment (1 ≤ i ≤ n) below a safe limit given by C i , i.e., we prescribe

where C i is an upper bound for w i . A lower bound could also be considered similarly, by taking h(x) = w i − C i .

Theorem 1: Consider dynamical system (1), function h in (7) and the corresponding set S given by (5) . The following safety-critical active intervention controller guarantees that S is forward invariant (safe) under dynamics (1) if g i (w) = 0, ∀w ∈ R n :

where ReLU(·) = max{0, ·} is the rectified linear unit,

and α > 0. Furthermore, the controller is optimal in the sense that it has minimum-norm control input.

Proof. According to [12] , the necessary and sufficient condition of forward set invariance is given by 1

∀t ≥ 0, where the derivative is taken along the solution of (1). If there exists a control input u(t) so that (10) is satisfied, then h is called a control barrier function. Substitution of (7) and (1) into (10) gives the safety condition

where ϕ i is given by (9) . The control input u(t) must satisfy (11) for all t ≥ 0. To keep control efforts minimal, one can achieve this by solving the quadratic program:

Based on the KKT conditions [33] , the explicit solution is

if g i (w(t)) = 0, which can be simplified to (8) . We remark that if g i (w) = 0, safety can be ensured by the help of extended control barrier functions as discussed for the safety guarantees of outlet compartments in Sec. III-B.

For example, to keep the infected population I below the limit I max for the SIR model given by (2) , one shall prescribe h(x) = I max − I, and (8) leads to the controller

Fig . 2 shows the dynamics of the closed control loop for the COVID-19 model fitted in Fig. 1 by prescribing 1 More precisely, α must be chosen as an extended class K function [12] , but we use a constant for simpler discussion and without loss of generality. Fig. 2 . Safety-critical active intervention control of the SIR model fitted in Fig. 1 to US COVID-19 data. The controller keeps the infected population under the prescribed limit Imax as opposed to the second wave of infection that was experienced during the summer of 2020.

I max = 200, 000 and using α = γ/10. Indeed, the safetycritical controller (red) applied from June 1, 2020 is able to keep the level of infection under the safe limit (red dashed), while gradually reducing mitigation efforts to zero. Meanwhile, the US experienced a second wave of infections (gray) in the summer of 2020, which was caused by the drop in mitigation efforts in June (see the blue control input).

Now consider the case where the j-th outlet compartment (1 ≤ j ≤ m) needs to be kept within the safe limit C j :

In the following theorem, we use a dynamic extension of control barrier functions to guarantee safety.

Theorem 2: Consider dynamical system (1), function h in (15) and the corresponding set S given by (5) . The following safety-critical active intervention controller guarantees that S is forward invariant (safe) under dynamics (1) iḟ h(x ( 0)) + αh(x(0)) ≥ 0 and if L g q j (w) = 0, ∀w ∈ R n :

and α > 0, α e > 0. Furthermore, the controller is optimal in the sense that it has minimum-norm control input.

Proof. We again use (10) as the necessary and sufficient condition for safety, where the following expression appears:

h e (x(t)) :=ḣ(x(t)) + αh(x(t)) = −(q j (w(t)) + r j (z(t))) + α(C j − z j (t)), (18) which puts the safety condition into the form h e (x(t)) ≥ 0, ∀t ≥ 0. However, the control input does not explicitly show up in (18) . Still, if there exists a control input that satisfieṡ h e (x(t)) ≥ −α e h e (x(t)),

then h e is an extended control barrier function [34] , [13] , whose 0-superlevel is forward invariant, that is, h e (x 0 ) ≥ 0 implies h e (x(t)) ≥ 0, ∀t > 0. Substitution of (18), (15) and (1) into (19) gives the extended safety condition

where ϕ e j is defined by (17) . This can be satisfied by a minnorm controller obtained from the quadratic program: The explicit solution of the quadratic program is

which is equivalent to (16) .

As an example of keeping outlet compartments safe, consider limiting the number of hospitalizations below H max and deaths below D max for the SIHRD model given by (4) . By choosing h(x) = H max − H, one can guarantee safety in terms of hospitalization based on (16) by the controller:

whereas prescribing h(x) = D max − D ensures safety by upper bounding deaths via:

Having synthesized controllers that keep selected compartments safe, let us now guarantee safety for multiple compartments at the same time: a set of multiplicative compartments I ⊂ {1, . . . , n} and a set of outlet compartments J ⊂ {1, . . . , m}. To formulate the safety condition, one can utilize (11) for any multiplicative compartment i ∈ I and (20) for any outlet compartment j ∈ J . Then, one needs to solve the corresponding quadratic program subject to all these constraints. In general, the quadratic program can only be solved numerically and one may need relaxation terms to satisfy multiple constraints [12] . However, analytical solutions can be found in some special cases, such as the one given by the following assumption. Assumption 1. Assume that the following terms have the same sign: sign(g i (w(t))) = sign(L g q j (w(t))) = −1, ∀i ∈ I, ∀j ∈ J , ∀t ≥ 0.

This assumption often holds for models where compartments need to be upper bounded for safety, e.g., the assumption holds for keeping E, I, R, H or D below a safe limit in the SIR, SEIR or SIHRD models. Under this assumption, one can state the following proposition.

Proposition 3: Consider dynamical system (1) with Assumption 1 and the controllers (8) and (16) that keep individual multiplicative compartments w i , i ∈ I ⊂ {1, . . . , n} and outlet compartments z j , j ∈ J ⊂ {1, . . . , m} safe using the control barrier functions in (7) and (15) . The following safety-critical active intervention controller guarantees safety for all compartments at the same time:

That is, one needs to take the maximum of the individual control inputs that keep each individual compartment safe.

Proof. If Assumption 1 holds, the safety conditions in (11) and (20) can be combined into one inequality:

Then, one can solve the quadratic program:

in the form:

This can be simplified to (25) based on (8) and (16) . Fig. 3 shows the closed loop response of the SIHRD model given by (4) that was fitted to US COVID-19 data [10] . The data about confirmed cases were scaled by the cube root of the positivity rate (positive per total tests) to account for the significant under-reporting of cases during the first wave of the virus (and cube root was applied to scale less aggressively). Starting from June 1, safetycritical active intervention control is applied to limit both the hospitalizations below H max = 40, 000 and the deaths below D max = 400, 000. Based on (25), we utilize the controller A HD (x) = max{A H (x), A D (x)} where A H and A D are given by (23) and (24) . The model and controller parameters are β 0 = 0.53 day −1 , γ = 0.14 day −1 , λ = 0.03 day −1 , ν = 0.14 day −1 , µ = 0.01 day −1 , N = 15×10 6 , τ = 9 days, α D = α e D = α H = (γ + λ + µ)/10 and α e H = ν/10. Safetycritical control is able to reduce mitigation efforts while keeping the system below the prescribed hospitalization and death bounds and preventing a second wave of the virus.

Controller (6) in Sec. III is designed based on feeding back the instantaneous state x(t) of the compartmental model. However, data about certain compartments is measured with delay due to the incubation period and testing delays. Thus, the instantaneous state x(t) may not be available for feedback, but the delayed state x(t − τ ) with measurement delay τ shall be used. If one implements A(x(t − τ )) instead of A(x(t)) for active intervention, a significant discrepancy between the delayed and instantaneous states can endanger safety. For example, the delay was identified to be τ = 11 days for the US COVID-19 data in Fig. 1 , while the infected population grew from a few thousands to more than a hundred thousand within 11 days in mid March. This difference significantly impacts safety-critical control. Thus, we propose a method to compensate delays via predicting the instantaneous state from the delayed one and we analyze how the prediction error affects safety.

We use the idea of predictor feedback control [22] , [23] , [24] to overcome the effect of delays. Namely, at each time moment t we use the data that are available up to time t − τ and we calculate a predicted state x p (t) that approximates the instantaneous state: x p (t) ≈ x(t). Then, we use the predicted state in the feedback law by applying A(x p (t)) ≈ A(x(t)). If the prediction is perfect (i.e., x p (t) = x(t)), safety is guaranteed even in the presence of delay according to Sec. III. Below we analyze how errors in the prediction affect safety.

The prediction can be done by any model-based or databased methods; see Example 4 for instance. At this point we only assume that the prediction error defined by

is bounded in the sense that e(t) ∞ ≤ ε for some ε ≥ 0. The prediction error leads to an input disturbance

relative to the nominal control input u(t) = A(x(t)), which yields the closed control looṗ w(t) = f (w(t)) + g(w(t))(u(t) + d(t)), z(t) = q(w(t)) + r(z(t)).

Note that for a Lipschitz continuous controller A with Lipschitz constant c, the input disturbance is upper bounded by d(t) ∞ ≤ c e(t) ∞ ≤ cε =: δ. The following theorem summarizes how the disturbance affects safety via the notion of input-to-state safety [35] . For simplicity, we state this theorem only for the safety of multiplicative compartments.

Theorem 4: Consider the closed loop dynamical system (31) , function h in (7) and the corresponding set S given by (5) . Assume that the nominal controller u(t) guarantees safety without the input disturbance d(t) by satisfying (11), while the input disturbance d(t) defined by (30) is bounded by d(t) ∞ ≤ δ. Then, set S is input-to-state safe in the sense that a larger set S d ⊇ S given by

is forward invariant (safe) under dynamics (31) , where

Proof. Similarly to (10) and (19), the necessary and sufficient condition for the invariance of S d is given bẏ

Substituting (33) and taking the derivative along the solution of (31) yields −ϕ i (w(t))−g i (w(t))(u(t)+d(t))+δ g i (w(t)) ∞ ≥ 0, (35) which indeed holds, since (11) and d(t) ∞ ≤ δ hold. How much larger set S d is compared to set S depends on the size δ of the disturbance that is related to the prediction error ε. If the prediction is perfect (x p (t) = x(t)), then ε = 0, δ = 0 and S d recovers S. However, if one implements a delayed state feedback controller without prediction (x p (t) = x(t − τ )), then ε and δ can be large, while S d can be significantly larger than the desired set S. Example 4. A possible model-based prediction can be done as follows. At each time moment t, we take the most recent available measurement x(t − τ ) and calculate the predicted state x p (t) by numerically integrating the ideal delay-free closed loop over the delay interval θ ∈ [t − τ, t]: w p (θ) = f (w p (θ)) + g(w p (θ))A(x p (θ)), z p (θ) = q(w p (θ)) + r(z p (θ)),

where x p = [w T p z T p ] T and the initial condition for integration is x p (t − τ ) = x(t − τ ). The results in Figs. 2 and 3 involve this kind of predictor feedback to compensate the delay τ . The simulations were carried out without considering uncertainties in the delay or other parameters, therefore the predicted state x p (t) matched the instantaneous state x(t) up to the numerical accuracy of integration. This allowed us to guarantee safety even in the presence of a significant delay.

V. CONCLUSIONS In this paper, we viewed compartmental epidemiological models as control systems where human actions (such as quarantining or social distancing) are considered as control input. By the framework of control barrier functions, we synthesized optimal safety-critical active intervention policies that formally guarantee safety against the spread of infection while keeping mitigation efforts minimal. We highlighted that time delays arising during state measurements can significantly affect safety-critical control, and we proposed predictor feedback to compensate the delays while preserving a certain level of input-to-state safety. We demonstrated our results on compartmental models fitted to US COVID-19 data, where we synthesized controllers to keep infection, hospitalization and deaths within prescribed limits. These controllers can help guide policy makers to decide when and how much mitigation efforts shall be reduced or increased.

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

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

Optimal control of an SIR model with delay in state and control variables

Time-optimal control strategies in SIR epidemic models

Can the COVID-19 epidemic be controlled on the basis of daily test reports

Estimating the impact of COVID-19 control measures using a Bayesian model of physical distancing

Quantifying the effect of quarantine control in Covid-19 infectious spread using machine learning

A feedback SIR (fSIR) model highlights advantages and limitations of infection-based social distancing

Modeling shield immunity to reduce COVID-19 epidemic spread

Safetycritical control of active interventions for COVID-19 mitigation

Control barrier function based quadratic programs with application to adaptive cruise control

Control barrier function based quadratic programs for safety critical systems

Control barrier functions: Theory and applications

On a characterization of flow-invariant sets

Barrier certificates for nonlinear model validation

Viability theory

Set-theoretic methods in control

A time delay dynamical model for outbreak of 2019-nCoV and the parameter identication

Initial simulation of SARS-CoV2 spread and intervention effects in the continental US

Risk assessment of novel coronavirus COVID-19 outbreaks outside China

Time delay compensation in unstable plants using delayed state feedback

Compensation of infinitedimensional input dynamics

Predictor feedback for delay systems: Implementations and approximations

The mathematics of infectious diseases

Estimation of the final size of the COVID-19 epidemic

Eects of latency and age structure on the dynamics and containment of COVID-19

A modified SEIR model to predict the COVID-19 outbreak in Spain and Italy: simulating control scenarios and multi-scale epidemics

Estimating and simulating a SIRD model of COVID-19 for many countries, states, and cities

A metapopulation network model for the spreading of SARS-CoV-2: Case study for Ireland

SafeGraph

Convex Optimization

Exponential control barrier functions for enforcing high relative-degree safety-critical constraints

American Control Conference (ACC)

Input-to-state safety with control barrier functions

The authors would like to thank Franca Hoffmann for her insights into compartmental epidemiological models and Gábor Stépán for discussions regarding non-pharmaceutical interventions in Europe. This research is supported in part by the National Science Foundation, CPS Award #1932091.