arXiv:2009.10262v1 [eess.SY] 22 Sep 2020
Safety-Critical Control of Compartmental Epidemiologica
l Models
with Measurement Delays
Tam ́as G. Moln ́ar
1
, Andrew W. Singletary
2
, G ́abor Orosz
1
,
3
and Aaron D. Ames
2
Abstract
— We introduce a methodology to guarantee safety
against the spread of infectious diseases by viewing epidem
i-
ological models as control systems and by considering human
interventions (such as quarantining or social distancing)
as con-
trol input. We consider a generalized compartmental model t
hat
represents the form of the most popular epidemiological mod
els
and we design safety-critical controllers that formally gu
arantee
safe evolution with respect to keeping certain populations
of
interest under prescribed safe limits. Furthermore, we dis
cuss
how measurement delays originated from incubation period a
nd
testing delays affect safety and how delays can be compensat
ed
via predictor feedback. We demonstrate our results by syn-
thesizing active intervention policies that bound the numb
er
of infections, hospitalizations and deaths for epidemiolo
gical
models capturing the spread of COVID-19 in the USA.
I. INTRODUCTION
The rapid spreading of COVID-19 across the world forced
people to change their lives and practice mitigation effort
s
at a level never seen before, including social distancing,
mask-wearing, quarantining and stay-at-home orders. Thes
e
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. Ther
e-
fore, it is important to focus mitigation efforts and determ
ine
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 mitiga-
tion efforts minimal. We use compartmental epidemiologica
l
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 epidemiologica
l
models as control systems has been proposed in the literatur
e
recently [3], [4], [5], and various models with varying
transmission rate [6], [7], [8], [9] have appeared to quanti
fy
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 contr
ol
strategies that guide human interventions so that certain
safety criteria (such as keeping infection, hospitalizati
on and
death below given limits) are fulfilled with minimal miti-
gation efforts. The approach is based on the framework of
1
Department of Mechanical Engineering, University of Michi
gan,
Ann Arbor, MI 48109, USA
molnart@umich.edu,
orosz@umich.edu
2
Department of Mechanical and Civil Engineering, Californi
a Institute
of Technology, Pasadena, CA 91125, USA
ames@caltech.edu,
asinglet@caltech.edu
3
Department of Civil and Environmental Engineering, Univer
sity of
Michigan, Ann Arbor, MI 48109, USA
Fig. 1. Illustration of the SIR model as control system and it
s fit to US
COVID-19 data [10]. Model parameters were estimated from co
mpartmental
data (right) by accounting for a measurement delay
τ
. The transmission rate
and the corresponding control input (left) were fitted to mob
ility data.
control barrier functions [11], [12] that leverages the the
ory
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 in-
troduces a generalized compartmental model, which covers
the class of the most popular epidemiological models. Sec-
tion III introduces safety critical control without consid
ering
measurement delays, while Sec. IV is dedicated to delay
compensation. Conclusions are drawn in Sec. V.
II. GENERALIZED COMPARTMENTAL MODEL
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
:
̇
w
(
t
) =
f
(
w
(
t
)) +
g
(
w
(
t
))
u
(
t
)
,
̇
z
(
t
) =
q
(
w
(
t
)) +
r
(
z
(
t
))
,
(1)
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 epi-
demiological 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 popula-
tions. 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
w
=
S
I
, f
(
w
) =
−
β
0
N
SI
β
0
N
SI
−
γI
, g
(
w
) =
β
0
N
SI
−
β
0
N
SI
,
z
=
R,
q
(
w
) =
γI,
r
(
z
) = 0
.
(2)
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 ove
r 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
w
=
S
E
I
, f
(
w
) =
−
β
0
N
SI
β
0
N
SI
−
σE
σE
−
γI
, g
(
w
) =
β
0
N
SI
−
β
0
N
SI
0
,
z
=
R,
q
(
w
) =
γI,
r
(
z
) = 0
.
(3)
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 hospitalizat
ion
rate
λ >
0
, the recovery rate
ν >
0
in hospitals and the death
rate
μ >
0
. Equation (1) yields the SIHRD model for
w
=
S
I
, f
(
w
)=
−
β
0
N
SI
β
0
N
SI
−
(
γ
+
λ
+
μ
)
I
, g
(
w
)=
β
0
N
SI
−
β
0
N
SI
,
z
=
H
R
D
, q
(
w
) =
λI
γI
μI
,
r
(
z
) =
−
νH
νH
0
.
(4)
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 cap-
turing 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,
I
(
t
−
τ
) +
R
(
t
−
τ
)
, 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 o
f
infection. The effects of the delay
τ
on safety-critical control
and its compensation will be discussed in Sec. IV.
III. SAFETY-CRITICAL CONTROL
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
:
S
:=
{
x
∈
R
n
+
m
:
h
(
x
)
≥
0
}
,
(5)
where
x
= [
w
T
, z
T
]
T
. 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
u
(
t
) =
A
(
x
(
t
))
(6)
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 withi
n
prescribed limits. First, we consider safety for multiplic
ative
compartments, and then for outlet compartments.
A. Safety Guarantees for Multiplicative Compartments
Consider keeping the
i
-th multiplicative compartment
(
1
≤
i
≤
n
) below a safe limit given by
C
i
, i.e., we prescribe
h
(
x
) =
C
i
−
w
i
,
(7)
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 follow-
ing safety-critical active intervention controller guara
ntees
that
S
is forward invariant (safe) under dynamics (1) if
g
i
(
w
)
6
= 0
,
∀
w
∈
R
n
:
u
(
t
)=
A
i
(
x
(
t
))=
−
sign(
g
i
(
w
(
t
)))ReLU
φ
i
(
w
(
t
))
|
g
i
(
w
(
t
))
|
,
(8)
where
ReLU(
·
) = max
{
0
,
·}
is the rectified linear unit,
φ
i
(
w
) =
f
i
(
w
)
−
α
(
C
i
−
w
i
)
(9)
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 con-
dition of forward set invariance is given by
1
̇
h
(
x
(
t
))
≥ −
αh
(
x
(
t
))
,
(10)
∀
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 condit
ion
−
φ
i
(
w
(
t
))
−
g
i
(
w
(
t
))
u
(
t
)
≥
0
,
(11)
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:
u
(
t
) =
A
i
(
x
(
t
)) = argmin
u
∈U
u
2
s
.
t
.
(11)
.
(12)
Based on the KKT conditions [33], the explicit solution is
u
(
t
)=
A
i
(
x
(
t
))=
(
0
if
−
φ
i
(
w
(
t
))
≥
0
,
−
φ
i
(
w
(
t
))
g
i
(
w
(
t
))
if
−
φ
i
(
w
(
t
))
<
0
,
(13)
if
g
i
(
w
(
t
))
6
= 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
A
I
(
x
) = ReLU
1
−
α
I
(
I
max
−
I
) +
γI
β
0
SI/N
.
(14)
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
I
max
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 safety-
critical 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 ze
ro.
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).
B. Safety Guarantees for Outlet Compartments
Now consider the case where the
j
-th outlet compartment
(
1
≤
j
≤
m
) needs to be kept within the safe limit
C
j
:
h
(
x
) =
C
j
−
z
j
.
(15)
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 follow-
ing safety-critical active intervention controller guara
ntees
that
S
is forward invariant (safe) under dynamics (1) if
̇
h
(
x
(
0)) +
αh
(
x
(0))
≥
0
and if
L
g
q
j
(
w
)
6
= 0
,
∀
w
∈
R
n
:
u
(
t
)=
A
j
(
x
(
t
))=
−
sign(
L
g
q
j
(
w
(
t
)))ReLU
φ
e
j
(
x
(
t
))
|
L
g
q
j
(
w
(
t
))
|
,
(16)
where
L
g
q
j
(
w
) =
∂q
j
∂w
(
w
)
g
(
w
)
,
φ
e
j
(
x
) =
∂q
j
∂w
(
w
)
f
(
w
) +
∂r
j
∂z
(
z
)
q
(
w
) +
r
(
z
)
+ (
α
+
α
e
)
q
j
(
w
) +
r
j
(
z
)
−
α
e
α
(
C
j
−
z
j
)
,
(17)
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 appea
rs:
h
e
(
x
(
t
)) :=
̇
h
(
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 satisfi
es
̇
h
e
(
x
(
t
))
≥ −
α
e
h
e
(
x
(
t
))
,
(19)
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
−
φ
e
j
(
x
(
t
))
−
L
g
q
j
(
w
(
t
))
u
(
t
)
≥
0
,
(20)
where
φ
e
j
is defined by (17). This can be satisfied by a min-
norm controller obtained from the quadratic program:
u
(
t
) =
A
j
(
x
(
t
)) = argmin
u
∈U
u
2
s
.
t
.
(20)
.
(21)
The explicit solution of the quadratic program is
u
(
t
)=
A
j
(
x
(
t
))=
(
0
if
−
φ
e
j
(
x
(
t
))
≥
0
,
−
φ
e
j
(
x
(
t
))
L
g
q
j
(
w
(
t
))
if
−
φ
e
j
(
x
(
t
))
<
0
,
(22)
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:
A
H
(
x
) = ReLU
1
−
α
e
H
α
H
(
H
max
−
H
)
λβ
0
SI/N
−
(
ν
−
α
H
−
α
e
H
)(
λI
−
νH
) + (
γ
+
λ
+
μ
)
λI
λβ
0
SI/N
,
(23)
whereas prescribing
h
(
x
) =
D
max
−
D
ensures safety by
upper bounding deaths via:
A
D
(
x
) = ReLU
1
−
α
e
D
α
D
(
D
max
−
D
)
μβ
0
SI/N
−
(
γ
+
λ
+
μ
−
α
D
−
α
e
D
)
μI
μβ
0
SI/N
.
(24)
C. Safety Guarantees for a Combination of Compartments
Having synthesized controllers that keep selected com-
partments safe, let us now guarantee safety for multiple
compartments at the same time: a set of multiplicative com-
partments
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 ca
n
only be solved numerically and one may need relaxation
terms to satisfy multiple constraints [12]. However, analy
tical
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 com-
partments 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 As-
sumption 1 and the controllers (8) and (16) that keep indi-
vidual 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 guarante
es safety
for all compartments at the same time:
u
(
t
) = max
i
∈I
,j
∈J
{
A
i
(
x
(
t
))
, A
j
(
x
(
t
))
}
.
(25)
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:
min
i
∈I
,j
∈J
−
φ
i
(
w
(
t
))
|
g
i
(
w
(
t
))
|
,
−
φ
e
j
(
x
(
t
))
|
L
g
q
j
(
w
(
t
))
|
+
u
(
t
)
≥
0
.
(26)
Then, one can solve the quadratic program:
u
(
t
) =
A
IJ
(
x
(
t
)) = argmin
u
∈U
u
2
s
.
t
.
(26)
(27)
in the form:
u
(
t
)=ReLU
−
min
i
∈I
,j
∈J
−
φ
i
(
w
(
t
))
|
g
i
(
w
(
t
))
|
,
−
φ
e
j
(
x
(
t
))
|
L
g
q
j
(
w
(
t
))
|
.
(28)
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 ac-
count 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, safety-
critical active intervention control is applied to limit bo
th 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 parame-
ters 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
. Safety-
critical control is able to reduce mitigation efforts while
keeping the system below the prescribed hospitalization an
d
death bounds and preventing a second wave of the virus.
IV. SAFETY UNDER MEASUREMENT DELAYS
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,
Fig. 3. The dynamics of the SIHRD model fitted to US COVID-19 da
ta
under the safety-critical active intervention policy that
keeps hospitalization
and deaths under the prescribed limits
H
max
and
D
max
.
the instantaneous state
x
(
t
)
may not be available for feed-
back, 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 endan-
ger 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.
A. Safety of Predictor Feedback Control
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 apply-
ing
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 data-
based methods; see Example 4 for instance. At this point we
only assume that the prediction error defined by
e
(
t
) :=
x
p
(
t
)
−
x
(
t
)
(29)
is bounded in the sense that
k
e
(
t
)
k
∞
≤
ε
for some
ε
≥
0
.
The prediction error leads to an input disturbance
d
(
t
) :=
A
(
x
p
(
t
))
−
A
(
x
(
t
))
(30)
relative to the nominal control input
u
(
t
) =
A
(
x
(
t
))
, which
yields the closed control loop
̇
w
(
t
) =
f
(
w
(
t
)) +
g
(
w
(
t
))(
u
(
t
) +
d
(
t
))
,
̇
z
(
t
) =
q
(
w
(
t
)) +
r
(
z
(
t
))
.
(31)
Note that for a Lipschitz continuous controller
A
with
Lipschitz constant
c
, the input disturbance is upper bounded
by
k
d
(
t
)
k
∞
≤
c
k
e
(
t
)
k
∞
≤
cε
=:
δ
. The following theorem
summarizes how the disturbance affects safety via the notio
n
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 sys-
tem (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
k
d
(
t
)
k
∞
≤
δ
. Then, set
S
is input-to-state safe in the
sense that a larger set
S
d
⊇
S
given by
S
d
:=
{
x
∈
R
n
+
m
:
h
d
(
x
)
≥
0
}
(32)
is forward invariant (safe) under dynamics (31), where
h
d
:
R
n
+
m
→
R
is defined by
h
d
(
x
) :=
h
(
x
) +
δ
α
k
g
i
(
w
(
t
))
k
∞
.
(33)
Proof.
Similarly to (10) and (19), the necessary and suffi-
cient condition for the invariance of
S
d
is given by
̇
h
d
(
x
(
t
))
≥ −
αh
d
(
x
(
t
))
.
(34)
Substituting (33) and taking the derivative along the solut
ion
of (31) yields
−
φ
i
(
w
(
t
))
−
g
i
(
w
(
t
))(
u
(
t
)+
d
(
t
))+
δ
k
g
i
(
w
(
t
))
k
∞
≥
0
,
(35)
which indeed holds, since (11) and
k
d
(
t
)
k
∞
≤
δ
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 imple-
ments a delayed state feedback controller without predicti
on
(
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
(
θ
))
,
(36)
where
x
p
= [
w
T
p
z
T
p
]
T
and the initial condition for integra-
tion is
x
p
(
t
−
τ
) =
x
(
t
−
τ
)
. The results in Figs. 2 and 3
involve this kind of predictor feedback to compensate the de
-
lay
τ
. The simulations were carried out without considering
uncertainties in the delay or other parameters, therefore t
he
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 contr
ol
input. By the framework of control barrier functions, we syn
-
thesized optimal safety-critical active intervention pol
icies
that formally guarantee safety against the spread of infect
ion
while keeping mitigation efforts minimal. We highlighted
that time delays arising during state measurements can
significantly affect safety-critical control, and we propo
sed
predictor feedback to compensate the delays while preservi
ng
a certain level of input-to-state safety. We demonstrated o
ur
results on compartmental models fitted to US COVID-19
data, where we synthesized controllers to keep infection,
hospitalization and deaths within prescribed limits. Thes
e
controllers can help guide policy makers to decide when and
how much mitigation efforts shall be reduced or increased.
ACKNOWLEDGMENT
The authors would like to thank Franca Hoffmann for
her insights into compartmental epidemiological models an
d
G ́abor St ́ep ́an for discussions regarding non-pharmaceut
ical
interventions in Europe. This research is supported in part
by the National Science Foundation, CPS Award #1932091.
R
EFERENCES
[1] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Fi
lippo,
A. Di Matteo, and M. Colaneri, “Modelling the COVID-19 epide
mic
and implementation of population-wide interventions in It
aly,”
Nature
Medicine
, 2020.
[2] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmu
nds,
S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday
et al.
, “Early
dynamics of transmission and control of COVID-19: a mathema
tical
modelling study,”
The Lancet Infectious Diseases
, vol. 20, pp. 553–
558, 2020.
[3] M. Elhia, M. Rachik, and E. Benlahmar, “Optimal control o
f an SIR
model with delay in state and control variables,”
ISRN Biomathematics
,
vol. 2013, p. 403549, Aug 2013.
[4] L. Bolzoni, E. Bonacini, C. Soresina, and M. Groppi, “Tim
e-optimal
control strategies in SIR epidemic models,”
Mathematical Biosciences
,
vol. 292, pp. 86–96, Oct 2017.
[5] F. Casella, “Can the COVID-19 epidemic be controlled on t
he basis
of daily test reports?”
arXiv preprint
, p. 2003.06967, 2020.
[6] S. C. Anderson, A. M. Edwards, M. Yerlanov, N. Mulberry, J
. Stock-
dale, S. A. Iyaniwura, R. C. Falcao, M. C. Otterstatter, M. A.
Irvine,
N. Z. Janjua, D. Coombs, and C. Colijn, “Estimating the impac
t
of COVID-19 control measures using a Bayesian model of physi
cal
distancing,”
medRxiv
, p. 2020.04.17.20070086, Jan 2020.
[7] R. Dandekar and G. Barbastathis, “Quantifying the effec
t of quaran-
tine control in Covid-19 infectious spread using machine le
arning,”
medRxiv
, p. 2020.04.03.20052084, 2020.
[8] E. Franco, “A feedback SIR (fSIR) model highlights advan
tages and
limitations of infection-based social distancing,”
arXiv preprint
, p.
2004.13216, 2020.
[9] J. S. Weitz, S. J. Beckett, A. R. Coenen, D. Demory, M. Domi
nguez-
Mirazo, J. Dushoff, C.-Y. Leung, G. Li, A. M ̆ag ̆alie, S. W. Pa
rk,
R. Rodriguez-Gonzalez, S. Shivam, and C. Y. Zhao, “Modeling
shield
immunity to reduce COVID-19 epidemic spread,”
Nature Medicine
,
2020.
[10] A. D. Ames, T. G. Moln ́ar, A. W. Singletary, and G. Orosz,
“Safety-
critical control of active interventions for COVID-19 miti
gation,”
IEEE
Access
, 2020, under review.
[11] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barri
er function
based quadratic programs with application to adaptive crui
se control,”
in
53rd IEEE Conference on Decision and Control
. IEEE, 2014, pp.
6271–6278.
[12] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control
barrier
function based quadratic programs for safety critical syst
ems,”
IEEE
Transactions on Automatic Control
, vol. 62, no. 8, pp. 3861–3876,
2016.
[13] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sre
enath,
and P. Tabuada, “Control barrier functions: Theory and appl
ications,”
in
2019 18th European Control Conference (ECC)
. IEEE, 2019, pp.
3420–3431.
[14] H. Brezis, “On a characterization of flow-invariant set
s,”
Communica-
tions on Pure and Applied Mathematics
, vol. 23, no. 2, pp. 261–263,
1970.
[15] R. Abraham, J. E. Marsden, and T. Ratiu,
Manifolds, tensor analysis,
and applications
. Springer Science & Business Media, 2012, vol. 75.
[16] S. Prajna, “Barrier certificates for nonlinear model va
lidation,”
Auto-
matica
, vol. 42, no. 1, pp. 117–126, 2006.
[17] J.-P. Aubin,
Viability theory
. Springer Science & Business Media,
2009.
[18] F. Blanchini and S. Miani,
Set-theoretic methods in control
. Springer,
2008.
[19] Y. Chen, J. Cheng, Y. Jiang, and K. Liu, “A time delay dyna
mical
model for outbreak of 2019-nCoV and the parameter identicat
ion,”
arXiv preprint
, p. 2002.00418, Feb 2020.
[20] S. Pei and J. Shaman, “Initial simulation of SARS-CoV2 s
pread
and intervention effects in the continental US,”
medRxiv
, p.
2020.03.21.20040303, 2020.
[21] P. Boldog, T. Tekeli, Zs. Vizi, A. D ́enes, F. A. Bartha, a
nd G. R ̈ost,
“Risk assessment of novel coronavirus COVID-19 outbreaks o
utside
China,”
Journal of Clinical Medicine
, vol. 9, no. 2, p. 571, 2020.
[22] W. Michiels and D. Roose, “Time delay compensation in un
stable
plants using delayed state feedback,” in
Proceedings of the IEEE
Conference on Decision and Control
, Orlando, FL, USA, 2001.
[23] M. Krstic and N. Bekiaris-Liberis, “Compensation of in
finite-
dimensional input dynamics,”
Annual Reviews in Control
, vol. 34,
no. 2, pp. 233–244, 2010.
[24] I. Karafyllis and M. Krstic,
Predictor feedback for delay systems:
Implementations and approximations
. Basel: Birkh ̈auser, 2017.
[25] H. W. Hethcote, “The mathematics of infectious disease
s,”
SIAM
review
, vol. 42, no. 4, pp. 599–653, 2000.
[26] M. Batista, “Estimation of the final size of the COVID-19
epidemic,”
medRxiv
, p. 2020.02.16.20023606, 2020.
[27] K. B. Blyuss and Y. N. Kyrychko, “Eects of latency and age
struc-
ture on the dynamics and containment of COVID-19,”
medRxiv
, p.
2020.04.25.20079848, Apr 2020.
[28] L. R. Lopez and X. Rodo, “A modified SEIR model to predict t
he
COVID-19 outbreak in Spain and Italy: simulating control sc
enarios
and multi-scale epidemics,”
medRxiv
, p. 2020.03.27.20045005, Jan
2020.
[29] J. Fern ́andez-Villaverde and C. I. Jones, “Estimating
and simulating
a SIRD model of COVID-19 for many countries, states, and citi
es,”
National Bureau of Economic Research, Working Paper 27128,
2020.
[30] R. Humphries, M. Spillane, K. Mulchrone, S. Wieczorek,
M. O’Riordain, and P. Hoevel, “A metapopulation network mod
el for
the spreading of SARS-CoV-2: Case study for Ireland,”
medRxiv
, p.
2020.06.26.20140590, 2020.
[31] “The COVID Tracking Project,” https://covidtracking
.com/, 2020,
[Online].
[32] “SafeGraph,” https://www.safegraph.com, 2020, [Onl
ine].
[33] S. Boyd and L. Vandenberghe,
Convex Optimization
. Cambridge
University Press, 2004.
[34] Q. Nguyen and K. Sreenath, “Exponential control barrie
r functions
for enforcing high relative-degree safety-critical const
raints,” in
2016
American Control Conference (ACC)
. IEEE, 2016, pp. 322–328.
[35] S. Kolathaya and A. D. Ames, “Input-to-state safety wit
h control
barrier functions,”
IEEE Control Systems Letters
, vol. 3, no. 1, pp.
108–113, 2019.