key: cord-0773217-2nc61ros
authors: Mondi'e, Sabine; Castanos, Fernando
title: Observer-based predictor for a SIR model with delays
date: 2020-10-22
journal: International Journal of Robust and Nonlinear Control
DOI: 10.1002/rnc.5522
sha: c76c37c78f4d2e38e4124e55fc46d5f3ccb131bc
doc_id: 773217
cord_uid: 2nc61ros

We propose an observer for a SIR epidemic model. The observer is then uplifted into a predictor to compensate for time delays in the input and the output. Tuning criteria are given for tuning gains of the predictor, while the estimation-error stability is ensured using Lyapunov-Krasovskii functionals. The predictor's performance is evaluated in combination with a time-optimal control. We show that the predictor nearly recovers the performance level of the delay-free system.

Epidemics are a good example of how reality challenges researchers, offering the opportunity to show the strength of existing techniques and develop new ones in fields as varied as medicine, biology, computational sciences, and mathematical system theory.

Epidemiological models have been primarily used for prediction purposes, while mitigation policies are usually decided based on exhaustive simulations. From the perspective of control theory, an epidemic is viewed as a dynamical system with controlled variables. Its model is an instrument for designing a control action that will achieve the desired outcome. Depending on the context, different assumptions on the model and other control objectives can be formulated. Most works focus on vaccination or treatment policies, with the goals expressed in an optimal-control framework [22] . In the context of the current covid pandemic, with vaccines unavailable at present, attention is shifting towards intervention policies based on social distancing measures. Considering the dramatic effect that extended lockdowns have on people and countries economies, a minimumtime control using social distancing measures and considering hospital capacity restrictions was recently presented for the SIR model [1] .

The SIR model, which we also consider in this paper, is arguably the simplest epidemiological model. However, it already exhibits many of the nonlinear characteristics that are present in more elaborate models. We make the model more realistic by adding features such as inaccurate and partial state measurements, and input and measurement delays. In recent months, delays in measurements and policy implementation have proved to be critical in the success or failure of government strategies. The former correspond to the time taken for the tests to be carried out, processed, verified, and made available in centralized databases. The latter correspond to the time it takes the population to adopt restrictions such as quarantine, social distancing habits, and mask use.

Due to the prevalence of delays in the feedback loops of control systems and the associated detrimental effects on performance, input and output delays have received sustained attention in the past decades. Some early strategies for compensating the delays are the Smith controller [23] , the transformation-based reduction approach [2] , and the time-domain predictor-based designs [16] . Based on present state information, predictors based on Cauchy's formula provide the state ahead of time [3] . They were formally shown to ensure closed-loop stability [14] , but their practical implementation reveals instabilities due to neutral phenomena related to the integrals' discretization in Cauchy's formula. These issues inspired new proposals such as filtered predictors [18, 13] , and truncated predictors [24] . If the present state is not entirely measurable, it can be replaced by an estimation, provided that the system is observable [11] .

A recent approach to the compensation of delays consists of modifying an observer to predict future states. It was first introduced for the case of full-state information [19] and later extended to the partial information scenario [15] . This approach, called observer-predictor, suffers from some drawbacks: It loses the exact nature of predictions obtained with Cauchy's formula and requires the inclusion of extra sub-predictors designed via LMI techniques. However, it has significant advantages: The observer has the same structure as the system (modulo an output injection term), thus avoiding integrals in the prediction formulae. Also, it is readily applicable to the case of partial state information in observable systems, especially when observers are readily available. Systems with state delays [25] or nonlinear systems [7] can be modified easily to successfully compensate for input or output delays.

To tackle the complexity due to partial state availability, delay, and measurement errors, we resort to a wide array of tools available to specialists in the field of control of dynamical systems. For the control, we use a recent optimal law [1] . The objective is not to steer the epidemics towards a desired equilibrium. Instead, the aim is to track an optimal trajectory. As a result, the dynamics for the estimation error are timevarying and time-delayed. The stability of such dynamics is addressed from both the perspective of classical frequency-domain quasipolynomial analysis [20, 17] , and from the perspective of time-domain Lyapunov-Krasowskii analysis [9, 8] . In particular, the system's time-varying nature is taken into account by embedding the system into a model with polytopic uncertainty [4, 10] .

In Section 2, we introduce the SIR model and discuss the issues we want to overcome. An ad hoc change of variable allows designing an observer addressing incomplete state information for the delay-free system. In Section 3, this observer is developed into an observer-based predictor for the system with input and output delays. The next two sections are devoted to tuning the observer: A simple criterion to tune the observer gains is given in Section 4, and conditions for the stability of the predic-tion error dynamics are given in Section 5. The impact of measurement errors is discussed in Section 6. We conclude with some remarks. We show the validity of our approach by discussing a SIR case study along with the paper, which is of interest in its own right.

Allow us to recall some standard notation used in the literature of time-delay systems.

Notation. P C([−η, 0], R n ) is the set of piece-wise continuous functions defined on the interval [−η, 0]. Consider a time-delay differential equatioṅ

The time-varying delay, η(τ ), is bounded by 0 < η(τ ) ≤η. Given an initial function ϕ ∈ P C([−η(0), 0], R n ) the solution is denoted by ε(τ, ϕ). The restriction of ε(τ, ϕ) on the interval [τ − η(τ ), τ ] is denoted by

We will make use of the trivial function 0η : θ → 0, θ ∈ [−η, 0]. We use the Euclidean norm · for vectors and the corresponding induced norm for matrices. For ϕ ∈ P C([−η, 0], R n ) we use the norm

ϕ(θ) .

The notation Q > 0 means that the symmetric matrix Q is positive definite.

We consider a state-space SIR model

Here, S > 0, I > 0 denote the susceptible and the infected, respectively. The model is normalized, hence S + I ≤ 1. The transmission rate, β ∈ [βmin, βmax] with βmax > βmin > 0, can be controlled by applying social distancing measures, but such measures take effect h1 units of time later. The only information available at time t is the number of infected people at time t − h2. The recovery/death rate, γ > 0, and the time delays, h1, h2, are assumed to be known. The are of course more sophisticated models. It is possible to include exposed individuals (infected but not infectious), to distinguish between symptomatic and asymptomatic, dead and recovered, etc. However, for epidemics the parameters of which have large levels of uncertainty, such as covid-19, a simple model with fewer parameters is preferable, as long as it is able to reproduce the main features of the epidemics (hospital saturation, lock-down effects, herd immunity, and so forth). A simpler model is also preferable when the objective is to devise decision strategies, rather that simulating long-term behavior.

Using only the history of y, we wish to produce predictionsŜ,Î such that lim t→∞ (Ŝ(t)),Î(t)) = (S(t + h1), I(t + h1)) .

Our motivation is that, if we have a feedback β * (S, I) that is known to perform correctly on the system without delays, we can set β = β * (Ŝ,Î) and expect to recover a similar performance 1 . For concreteness, we consider the optimal-time control strategy described by Angulo et al [1] . For the SIR model, the basic (unmitigated) reproduction number is computed as R0 = βmax/γ, while the controlled (mitigated) reproduction number is Rc = βmin/γ. Suppose that the health system capacity of a given city is limited to Imax infected people. The strategy that ensures I(t) ≤ Imax and achieves herd immunity in a minimal time is

We consider a recovery rate γ = 1/7 with time units given in days [1] . For illustration purposes, we consider the case of Mexico City. The number of ICU beds is such that Imax = 12.63 × 10 −3 (see [1] ). We take the reproduction numbers as R0 = 1.7 and Rc = 1.1 2 . This gives βmax = 1.7/7 and βmax = 1.1/7 .

The optimal response, achieved with full state-feedback in the absence of delays, is shown in Fig. 1 . The optimal strategy is to allow the epidemic to run free until it reaches the sliding curve I = Φ(S). The state is then driven along this curve towards the region of herd immunity (yellow rectangle) where the intervention finally stops.

Suppose now that there is a delay of h1 = 3 days in the control action. Fig. 2 confirms the appearance of the so-called chattering effect, which should not be surprising given the discontinuous nature of (2). We can expect the performance of the closed-loop system to deteriorate even further when only the number of infected people is available for measurement, and when such measurements are also subject to important delays. The objective of the predictor, developed in the following section, is to mitigate these unfavorable effects. 1 The full analysis would have to be performed, of course. 

To attain our objective we follow the approach in which an observer for the delay-free system is constructed in a first step, and then developed into a predictor in a second one [19, 25, 6 ].

We begin by making the temporary assumption h1 = h2 = 0. Note that, by setting x1 = ln(y) and x2 = S, we obtain the model

with ρ(x) = −x2e x 1 . The main advantage over the original model is that the first equation is affine in the state. We can then write the simple observer

where α1, α2 will be defined below (Proposition 1). Consider the error x = x −x. By using the expansion

we can write the error dynamics as d dt

with

.

Note that, since we are linearizing the estimation error around a trajectory (rather than an equilibrium), the linearized system is time-varying. Fortunately, we can ensure its stability with a simple quadratic Lyapunov function.

Then (3) is quadratically stable 3 . Proof. Consider the candidate Lyapunov function V (x) =x Px with P = P > 0 the solution of the Lyapunov equation

The time-derivative of V along the trajectories of (3) iṡ

The restrictions (S, I) ∈ [0, 1] 2 , S + I ≤ 1 imply that SI ≤ 1/4, so the first leading principal minor of Q(S, I), 2(1 − SI)α1α2, is strictly positive.

Regarding the second leading principal minor, we have

Using again SI ≤ 1/4 we obtain the bound

Finally, the condition (4) ensures that |Q(S, I)| > 0, so that V is indeed a Lyapunov function.

In the original coordinates, the observer takes the form

A simulation of the observer's performance is shown in Fig. 3 . The observer gains, α = α1 α2 = 4 1 were chosen to satisfy (4). At first, the incidence of infection exceeds Imax by about 11 %, but then the epidemics behave as desired.

Consider again the input delay h1 = 3 days, and suppose now there is a measurement delay of h2 = 7 days [5] . The combined effect of both delays is disastrous. As illustrated in Fig. 4 , the hospital capacity is exceeded by 140 %. This motivates the design of the predictor presented in the sequel. 

Let us rewrite the observer in the original coordinates and remove the zero-delays assumption. This gives the predictor

with h = h1 + h2. Considering that y(t) = I(t − h2), the error variables

evolve according to the dynamics

with ψ(S, I,x) = ((S −x2)e −x 1 − S)I .

System (9) can be written in the general form

where the matrices are defined by

and

Since β multiplies all the right-hand side of (11), we can do away with it by rescaling time, much in the spirit of perturbation theory [12, Ch. 10] . Define the new time-scale

Since β is strictly positive, g is strictly increasing, invertible and τ is indeed a time-scale. Let us define the new state ε(τ ) =x(g −1 (τ )) and note that, by the Inverse Function Theorem, we have

Applying the chain rule and (14), we see that the new state evolves according to d dτ

with

and

The last equation follows from the condition

There are two main difficulties in establishing the stability of (15): The time-varying nature of B0, and the presence of the delayed state. The former difficulty will be addressed by formulating (15) in the framework of polytopic differential inclusions [4, Ch. 5] . In order to do so, we will focus on the state-space rectangle [0, 1]×[0,Ī], where Imax ≤Ī ≤ 1. When the state is restricted to such rectangle, B0 varies within a fixed polytope of matrices, i.e.,

where

and Co stands for convex closure, that is,

An obvious necessary condition for the stability of the polytopic model (15)-(16) is the stability of its linearized vertices, ε(τ ) = Ciε(τ ) + B1ε(τ −η) , i = 1, 2, 3.

The characteristic equations of the vertices are p1(s) = s 2 + (sα1 + α2) e −ηs p2(s) = s 2 + (s +Ī)α1 + α2 e −ηs + sĪ p3(s) = s 2 + (s +Ī)α1 + α2 e −ηs + (s + 1)Ī .

We will work out the stability/instability boundaries of these quasipolynomials in the space of parameters (α1, α2). According to the Dpartition method [20] , these boundaries correspond to roots crossing the imaginary axis of the complex plane. When the crossing root is real (s = 0), the boundaries are If a crossing root is imaginary (s = jω), the boundaries satisfy p1 :

ω sin(ωη) cos(ωη) ω cos(ωη) − sin(ωη) α = ω 2 0 p2 : ω sin(ωη) +Ī cos(ωη) cos(ωη) ω cos(ωη) −Ī sin(ωη) − sin(ωη) α = ω 2 −ωĪ p3 : ω sin(ωη) +Ī cos(ωη) cos(ωη)

The form and disposition of the stability regions depend onη and I. Figure 5 shows their boundaries for a relatively small delay,η = 0.5, and the largest incidence,Ī = 1. Figure 6 shows the boundaries of the stability regions for a larger delay,η = 2.5, but a smaller incidence,Ī = 0.03. Since stability of the three vertices is a necessary condition for the stability of (15)-(16), we require α to be placed at the intersection of the three regions (boundaries drawn in black), for example, at an approximate centroid of the intersection (marked with ).

Setting α as described in the previous section, i.e., at the intersection of the stability regions, only ensures that a necessary condition for stability is satisfied. We will now exploit the polytopic nature of (3) and the fact that stabilty is ensured by the existence of a Lyapunov-Krasowskii functional that is common to all the vertices of the polytope.

We will begin by summarizing a general assertion from the book by Fridman [8] .

.

where the terms are such that the overall matrix is symmetric.

The lemma is easily proved by differentiating V and applying Jensen's lemma [8, Ch. 3] . 

where C, P, R, S, P2 and P3 are 2 × 2 matrices. Suppose that there exists P > 0, R ≥ 0, S ≥ 0 and P2, P3 such that Q(Ci, P, R, S, P2, P3) < 0 , i = 1, 2, 3 .

Then, the trivial solution of the observer error-dynamics (15) is asymptotically stable.

Proof. The proof follows the descriptor approach [8] , but we pay special attention to the nonlinear terms in (15) and incorporate the time-varying nature of the system. Consider again (18) and note that, for any P2, P3 ∈ R 2×2 , 2 ε(τ ) P 2 +ε(τ ) P 3 · B0(τ )ε(τ ) + B1ε(τ − η(τ )) + H(τ, ε(τ ))ε(τ ) −ε(τ ) = 0 or, in matrix form,

so adding (22) to (19) giveṡ

Because of (16), there exists continuous real-valued functions ci such that

Since C appears affinely in (20) , we have [8, Remark 3.6]

By (21) , the derivative of V is negative definite in a neighborhood of the origin and asymptotic stability follows.

For our example, we take a conservative approach and set η = 5 > (h1 + h2)βmax = 2.4 andĪ = 30 × 10 −3 (recall that Imax = 12.63 × 10 −3 ). The gain α = 0.115 0.005 lays within the intersection of the stability regions of p1, p2 and p3.

The LMI (21) was solved using SCS [21] . A solution is so the prediction error converges to zero asymptotically. This is illustrated in Fig. 7 . The hospital capacity is now exceeded by only 7.8 %.

6 On the effect of measurement errors A frequent situation in epidemics is poor output variable measurement, mainly due to underregistration. A sound assumption is that the output is proportional to the measured variable, I(t − h2). The proportion may be time-varying but always less than 1. It is described as

The logarithmic term in the prediction error dynamics (8) is now ln

Define d(t) = ln(m(t)) and observe that, by (23) , is such that |d(t)| ≤d withd = ln 1 1 − δ .

The prediction error now has the dynamicṡ ε(τ ) = B0(τ )ε(τ ) + B1ε(τ − η(τ )) + H(τ, ε(τ ))ε(τ ) − αd(τ − η(τ )) . (24) To analyze the effect of the measurement error on the system response, we compute the time derivative of the functional (18), now along the trajectories of system (24) . Following the same steps as in the proof of Theorem 1, we obtaiṅ

ci(τ )E (τ )Q(Ci; P, R, S, P2, P3)E(τ ) − 2 ε(τ ) P 2 +ε(τ ) P 3 αd(τ − η(τ )) + O( E(τ ) 3 ) .

From the order relation 2 ε(τ ) P 2 +ε(τ ) P 3 αd(τ − η(τ )) = O(d · E(τ ) ) we conclude that, ford small enough, the solutions are ultimately bounded with ultimate bound proportional tod. We now simulate the effects of the measurement noise described above. We take δ = 0.5 and define m(t) a random variable with a beta distribution. We perform simulations for the predictor and the observer (see Fig. 8 ). The measurement errors result in a higher number of infected people. The poor performance of the observer had already been established. Noise simply makes it worse. In the case of the predictor, the increment is relatively low if we take into account the high amplitude of the error.

We have presented a predictor for systems with large input delays. The predictor is designed using the observer-predictor methodology introduced in the past years. In contrast with existing proposals, we present simple tuning criteria for ensuring the simultaneous stability of the polytopic model's vertices. The stability of the complete polytopic model is then established formally with the help of a Lyapunov-Krasovskii functional. The functional can also be used to assess the sensitivity of the predictor with respect to error measurements. It is worthy of mentioning that it can also be used to find estimates of the domain of attraction, robustness bounds for the delay, and parameter uncertainties, among other problems of interest.

Our current research concerns the extension of the present results, obtained for two-dimensional systems, to n-dimensional ones, and to apply them to the study of more elaborated epidemic models involving states such as vaccinated individuals and asymptomatic ones.

A simple criterion to design optimal nonpharmaceutical interventions for epidemic outbreaks. medRxiv

Linear systems with delayed controls: a reduction

Differential-difference equations

Linear Matrix Inequalities in System and Control Theory

Forecasting hospital demand during covid-19 pandemic outbreaks

Prediction-based control for nonlinear systems with input delay

prediction based control for nonlinear systems with input delay. Mathematical Problems in Engineering

Introduction to Time-Delay Systems, Analysis and Control

Stability of timedelay systems

Parameterdependent lyapunov functional for stability of time-delay systems with polytopic-type uncertainties

Stabilization of nonlinear delay systems using approximate predictors and high-gain observers

Nonlinear Systems

Predictor-based controls: the implementation problem

Boundary control of PDEs: A course on backstepping designs

Prediction-based control for lti systems with uncertain time-varying delays and partial state knowledge

Finite spectrum assignment problem for systems with delays

Stability and Stabilization of Time-Delay Systems. Society for Industrial and Applied Mathematics

Finite spectrum assignment of unstable time-delay systems with a safe implementation

Closed-loop control of dead time systems via sequential sub-predictors

D-subdivisions and spaces of quasipolynomials

Conic optimization via operator splitting and homogeneous self-dual embedding

Optimal control in epidemiology

A controller to overcome dead time

Truncated predictor feedback for linear systems with long time-varying input delays

Stabilization of linear systems with both input and state delays by observer-predictors