# Estimating parameters and states of an epidemic model using observers

Frédéric Hamelin, Abderrahman Iggidr, Alain Rapaport, Gauthier Sallet, Max

O. Souza

► To cite this version:

Frédéric Hamelin, Abderrahman Iggidr, Alain Rapaport, Gauthier Sallet, Max O. Souza. Estimating parameters and states of an epidemic model using observers. 2022. hal-04038203

HAL Id: hal-04038203

https://inria.hal.science/hal-04038203

Preprint submitted on 20 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Estimating parameters and states of an epidemic model using observers

F. Hamelin* A. Iggidr** A. Rapaport*** G. Sallet**
M. Souza****

* Institut Agro, Rennes, France (e-mail: frederic.hamelin@agrocampus-ouest.fr)

** Univ. de Lorraine, CNRS, Inria, IECL, F-57000 Metz, France
(e-mail: Abderrahman.Iggidr@inria.fr, gauthier.sallet@univ-lorraine.fr)

*** MISTEA, Univ. Montpellier, INRAE, Inst. Agro, France (e-mail: alain.rapaport@inrae.fr)

**** Instituto de Matemática e Estatística, Universidade Federal Fluminense, Niterói - RJ, 24210-201, Brasil (e-mail: maxsouza@id.uff.br)

**Abstract:** We consider the practical parameters and state estimation problem for the well-known SIR epidemic model with an additional compartment Q of the sub-population of infected individuals that are placed in quarantine (SIQR model). We assume that the flow of individuals placed in quarantine and the size of the quarantine population are known at any time. Using this information as the measurable output of the system, we build two kinds of observers that give dynamical estimates of the parameters and the size of infected individuals.

**Keywords:** Parameters estimation, identifiability, observability, parameters identification.

## 1. INTRODUCTION

Many papers in epidemiology proposing a mathematical model using dynamical systems, face the problem of parameter estimation. In general, some parameters are given, extracted from the literature, while the remaining unknown parameters are estimated by fitting the model to some observed data, usually by means of an optimization algorithm based on least squares or maximum likelihood methods. Here we are interested in the design of some observers to estimate some parameters as well as the unmeasured state variable for a classical epidemic model

The practical parameters identification with state reconstruction using observers for epidemic models has been the goal of, relatively, few papers: Diaby et al. (2015); Abdelhedi et al. (2014); Bichara et al. (2014); Degue et al. (2021); Bichara et al. (2019); Iggidr and Souza (2019); Bliman et al. (2018); Diouf et al. (2019); Hamelin et al. (2020) but these approaches are rarely considered in the literature compared to least-square methods, differently to other applications domains, such as bio-processing or manufacturing industries.

The global health crisis of COVID-19 outbreak has led to a spectacular resurgence of interest in this type of models, but with specificities related to the detection and isolation of infected individuals. This is why we revisit here the issues of parameters and state estimation for an extended 'SIQR' model Hethcote et al. (2002) for which such an analysis had not yet been performed (to the best of our knowledge). In particular, we tackle the task of proposing

a practical strategy for reconstructing the values of the parameters.

## 2. THE MODELS

Inspired by Hethcote et al. (2002); Nuño et al. (2008), we consider the classical SIR model (see for instance Li (2018)), where S, I, R denote the size of the populations of respectively susceptible, infected and recovered individuals, with an additive compartment where Q denotes the size of the population of identified and isolated infectious individuals that have been removed from the infected population and placed in quarantine:

$$ 
\left\{ 
\begin{aligned}
\dot{S} &= -\beta S \frac{I}{N} \\
\dot{I} &= \beta S \frac{I}{N} - (\rho + \alpha)I \\
\dot{Q} &= \alpha I - \rho Q \\
\dot{R} &= \rho I + \rho Q
\end{aligned}
\right. 
\qquad (1) $$

Note that one has

$$ S(t) + I(t) + Q(t) + R(t) = N, \quad \forall t \geq 0. $$

Since $S(t) + I(t) + Q(t) + R(t)$ is constant and R has no influence on the first three equations, it is sufficient to consider the following system

$$ 
\left\{ 
\begin{aligned}
\dot{S} &= -\beta S \frac{I}{N} \\
\dot{I} &= \beta S \frac{I}{N} - (\rho + \alpha)I \\
\dot{Q} &= \alpha I - \rho Q
\end{aligned}
\right. 
\qquad (2) $$
It can be proved that the open set

$$ \Omega = \{S > 0, I > 0, Q > 0, S + I + Q < N\} $$

is positively invariant.

These models have three parameters: the infectivity parameter $\beta$, the recovery rate $\rho$, that we assume to be identical for the infected populations placed in quarantine or not, and the rate of placement in quarantine $\alpha$. These parameters are unknown but we assume the following hypothesis.

**Assumption 1.** The basic reproduction number $R_0$ verifies

$$ R_0 := \frac{\beta}{\rho + \alpha} > 1. $$

This assumption implies that the epidemic can spread in the population i.e. at initial time with $S(0) = N - I(0)$ close to $N$ one has $\dot{I}(0) > 0$.

### 3. PARAMETERS ESTIMATION USING OBSERVERS

The application of classical least-square method to estimate the parameters of model (2) do not give satisfactory results. This is why we have opted for a dynamical estimation with the help of observers. The use of observers, although usually dedicated to state estimation possesses the advantage to tune the speed of error decay. Moreover, the choice of a speed-accuracy compromise can be balanced thru simulations with synthetic data corrupted by noise, before facing real data.

Note also that for large times, the solutions of (2) converge asymptotically to non-observable non-identifiable states when $I$ and $Q$ are both null. Consequently, we do not look precisely for results about asymptotic convergence of the error (as this is usually done in observers theory), but rather for an exponential decay of the error during initial transients.

In this section, we shall consider the additional hypothesis

**Assumption 2.** One has

$$ \alpha \leq \rho $$

that means that the rate of placement in quarantine is not larger than the recovery rate, which is often the case in epidemic regimes.

#### 3.1 Using High-Gain observers

Parameters $\alpha, \beta, \rho$ are unknown but are positive and bounded. Let us assume that some bounds, even loose, are available.

**Assumption 3.** There exist known positive numbers $\alpha_m < \alpha_M, \beta_m < \beta_M, \rho_m < \rho_M$ such that the unknown parameters $\alpha, \beta, \rho$ satisfy

$$ \alpha \in [\alpha_m, \alpha_M], \quad \beta \in [\beta_m, \beta_M], \quad \rho \in [\rho_m, \rho_M]. $$

**Theorem 4.** Assume that the parameters $\alpha, \beta, \rho$ fulfill Assumption 3. If the parameter $\rho$ is known, then for $\epsilon > 0$ sufficiently small, there exists $T_{\epsilon} > 0$ with $T_{\epsilon} \to +\infty$ when $\epsilon \to 0$ such that for $\theta > 0$ large enough the following quantities provide exponential estimators for the

parameters $\alpha, \beta$ and variables $S(t), I(t)$, respectively, for $t \in [0, T_{\epsilon}]$:

$$ \hat{\alpha}(t) = \frac{\hat{z}_3(t)^2}{\max\left(\hat{z}_4(t) - \hat{z}_2(t)\hat{z}_3(t), \epsilon \frac{\beta_m}{\alpha_M N} y_1(t)\right)} - \hat{z}_2(t) - \rho \quad (3) $$

$$ \hat{\beta}(t) = N e^{-\hat{z}_1(t)} \left( \frac{(\hat{z}_2(t) - \rho)(\hat{z}_4(t) - \hat{z}_2(t)\hat{z}_3(t))}{\min(\hat{z}_3(t), -\epsilon)} - \hat{z}_3(t) \right) \quad (4) $$

$$ \hat{S}(t) = N \frac{\hat{z}_2(t) + \rho + \hat{\alpha}(t)}{\max(\hat{\beta}(t), \beta_m)} \quad (5) $$

$$ \hat{I}(t) = \frac{e^{\hat{z}_1(t)}}{\max(\hat{\alpha}(t), \alpha_m)} \quad (6) $$

where the functions $\hat{z}_i(t)$ are solutions of the following observer dynamical system

$$ \dot{\hat{z}} = \begin{bmatrix} \dot{\hat{z}}_1 \\ \dot{\hat{z}}_2 \\ \dot{\hat{z}}_3 \\ \dot{\hat{z}}_4 \end{bmatrix} = \begin{bmatrix} \hat{z}_2 \\ \hat{z}_3 \\ \hat{z}_4 \\ \hat{z}_3^2 + \hat{z}_2 \hat{z}_4 + (\hat{z}_4 - \hat{z}_2 \hat{z}_3) \left( \hat{z}_2 + \frac{\hat{z}_4}{\hat{z}_3} \right) \end{bmatrix} - \begin{bmatrix} 4\theta \\ 6\theta^2 \\ 4\theta^3 \\ \theta^4 \end{bmatrix} (\hat{z}_1 - \log(y_1(t))), \quad (7) $$

where

$$ \begin{aligned} \hat{z}_2 &= sat_{-\beta_M - \rho_M - \alpha_M} (\hat{z}_2), \quad \hat{z}_3 = sat_{-\beta_M} (\hat{z}_3), \\ \hat{z}_4 &= sat_{-\beta_M (\beta_M + \rho_M + \alpha_M)} (\hat{z}_4). \end{aligned} \quad (8) $$

The proof of Theorem 4 will be made in several steps.

We first consider an observer with the single measurement $y_1(t) = \alpha I(t)$. This observer will deliver estimates of the following functions of parameters and state variable of the model (2): $\beta \frac{S(t)}{N} - \rho - \alpha, -\beta^2 \frac{S(t)I(t)}{N^2}$, and $-\beta^2 \frac{S(t)I(t)}{N^2} (\beta \frac{S(t)-I(t)}{N} - \rho - \alpha)$. More precisely we have the following result.

**Proposition 5.** Fix parameters $N, \alpha, \beta, \rho$ that fulfill Assumption 3 and initial condition $(S_0, I_0) \in (\mathbb{R}_+^*)^2$ with $S_0 + I_0 \le N$. Let $(S(\cdot), I(\cdot))$ be the solution of the system (2) and $y_1(\cdot) = \alpha I(\cdot)$. For $\epsilon > 0$ sufficiently small, there exists $T_{\epsilon} > 0$ with $T_{\epsilon} \to +\infty$ when $\epsilon \to 0$ such that for $\theta > 0$ large enough the solutions of the observer dynamical system given by (7-8) are such that the vector

$$ E(t) = \hat{z}(t) - \begin{bmatrix} \log(y_1(t)) \\ \beta \frac{S(t)}{N} - \rho - \alpha \\ -\beta^2 \frac{S(t)I(t)}{N^2} \\ -\beta^2 \frac{S(t)I(t)}{N^2} \left( \beta \frac{S(t) - I(t)}{N} - \rho - \alpha \right) \end{bmatrix} \quad (9) $$

converges exponentially to 0 on $[0, T_{\epsilon}]$, in the sense that there exists $K > 0$ such that

$$ \|E(t)\| \le K \|E(0)\| \exp \left( -\frac{\theta}{3} t \right), \quad t \in [0, T_{\epsilon}] \quad (10) $$

for any initial condition $\hat{z}_0 \in \mathbb{R}^4$.
**Proof.** Clearly, the output $y_1(t)$ remains positive for any $t \ge 0$. One can then consider the variable

$$z_1(t) = \log(y_1(t)), \quad t \ge 0 \quad (11)$$

and its three first time derivatives, denoted $z_2, z_3, z_4$. From equation (2), one has the following expressions

$$z_2 = \dot{z}_1 = \beta \frac{S}{N} - \rho - \alpha, \quad (12)$$

$$z_3 = \dot{z}_2 = -\beta^2 \frac{SI}{N^2}, \quad (13)$$

$$z_4 = \dot{z}_3 = -\beta^2 \frac{SI}{N^2} \left( \beta \frac{S}{N} - \rho - \alpha \right) + \beta^3 \frac{SI^2}{N^3} \\ = z_2 z_3 - \beta \frac{I}{N} z_3. \quad (14)$$

Note that $z_3(t)$ is negative for any $t \ge 0$ along the solutions of (2) with $S(0) > 0, I(0) > 0$, and thus cannot reaches 0 in finite time. Therefore, for any $\epsilon > 0$ such that $\epsilon < -z_3(0)$, there exists $T_\epsilon > 0$ such that $z_3(t) < -\epsilon$ for any $t \in [0, T_\epsilon]$, with $T_\epsilon \to +\infty$ when $\epsilon \to 0$. Moreover, from expressions (12), (13), (14) and Assumption 3, one obtains that the variables $z_2, z_3, z_4$ are bounded with

$$z_2(t) \in [-\rho_M - \alpha_M, \beta_M - \rho_M - \alpha_M], \quad (15)$$

$$z_3(t) \in [-\beta_M, -\epsilon], \quad (16)$$

$z_4(t) \in [-\beta_M(\beta_M - \rho_M - \alpha_M), \beta_M(\beta_M + \rho_M + \alpha_M)], (17)$ for any $t \in [0, T_\epsilon]$. Let us now write the time derivative of $z_4$ and show that it can be expressed as a function of $z_2, z_3, z_4$:

$$\dot{z}_4 = \dot{z}_3^2 + z_2 z_4 + \frac{z_4 - z_2 z_3}{z_3} (z_2 z_3 + z_4) \\ z_3^2 + z_2 z_4 + (z_4 - z_2 z_3) \left(z_2 + \frac{z_4}{z_3}\right) := \phi(z).$$

Knowing the bounds (15),(16),(17), one can write equivalently

$$\dot{z}_4(t) = \tilde{\phi}(z(t)) := \phi(\sigma(z(t))), \quad t \in [0, T_\epsilon]$$

with

$$\sigma(z) = \begin{bmatrix} z_1 \\ \text{sat}_{-\rho_M - \alpha_M}^{\beta_M - \rho_M - \alpha_M}(z_2) \\ \text{sat}_{-\beta_M}^{-\epsilon}(z_3) \\ \text{sat}_{-\beta_M(\beta_M + \rho_M + \alpha_M)}^{\beta_M(\beta_M + \rho_M + \alpha_M)}(z_4) \end{bmatrix}.$$

Note that the function $\tilde{\phi}$ is well-defined and globally Lipschitz on $\mathbb{R}^4$. Therefore, one can consider the dynamical system

$$\dot{z} = F(z) := \underbrace{\begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{bmatrix}}_{A_4} z + \underbrace{\begin{bmatrix} 0 \\ 0 \\ 0 \\ \tilde{\phi}(z) \end{bmatrix}}_{C_4}, \quad \tilde{y} = \underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix}}_{C_4} z$$

in whole $\mathbb{R}^4$ with output $\tilde{y}$, which is in the "canonical observable form" Gauthier and Kupka (2001). We recall that for such system, the dynamical system

$$\dot{z} = F(\hat{z}) - S_{\theta}^{-1}(C\hat{z} - \tilde{y}), \quad \hat{z}(0) = \hat{z}_0$$

where the symmetric matrix $S_\theta$ is solution of the Lyapunov equation

$$\theta S_{\theta} + A_{4}^{\top} S_{\theta} + S_{\theta} A_{4} - C_{4}^{\top} C_{4} = 0, \quad \theta > 0 \quad (18)$$

is a globally exponential observer, provided that $\theta$ is large enough. The proof is based on the Lyapunov function $V(e) = e^\top S_\theta e$ which verifies $\frac{d}{dt}V(e) \le -(\theta/3)V(e)$ for large enough $\theta$, where $e = \hat{z} - z$ is the error vector (see Gauthier et al. (1992); Gauthier and Kupka (2001) for more details), whatever is $\hat{z}_0 \in \mathbb{R}^4$. A lengthy but straightforward computation gives the following expression of the inverse of the solution $S_\theta$ of equation (18)

$$S_{\theta}^{-1} = \begin{bmatrix} 4\theta & 6\theta^2 & 4\theta^3 & \theta^4 \\ 6\theta^2 & 14\theta^3 & 11\theta^4 & 3\theta^5 \\ 4\theta^3 & 11\theta^4 & 10\theta^5 & 3\theta^6 \\ \theta^4 & 3\theta^5 & 3\theta^6 & \theta^7 \end{bmatrix}$$

from which one deduces the equations (7), (8) and the estimation error (10) of the vector (9).

Note that any solution of the system (2) converge asymptotically to a disease-free equilibrium i.e. with $I = 0$, for which the system is not identifiable nor observable. Therefore one cannot expect to have an exact convergence of the observer (7) when $t$ tends to $+\infty$. However, in practice, its speed of convergence is tuned tanks to the choice of the parameter $\theta$ to obtain a small estimation error much before the state reaches a neighborhood of the asymptotic steady-state. Moreover, the innovation term $\hat{z}_1(t) - \log(y_1(t))$ informs about the effective convergence: when it remains small, we know that the estimation error $E(t)$ is small.

From the results of Proposition 5, one gets the following estimators.

**Corollary 6.** Under Assumptions of Proposition 5, the numbers

$$\beta^{\dagger}(t) = \frac{\hat{z}_3(t)^2}{\max\left(\hat{z}_4(t) - \hat{z}_2(t)\hat{z}_3(t), \epsilon \frac{\beta_m}{\alpha_M N} y_1(t)\right)} \quad (19)$$

$$\mathcal{R}_0^{\dagger}(t) = \frac{\hat{z}_3(t)^2}{\max\left(\hat{z}_3(t)^2 - \hat{z}_2(t)[\hat{z}_4(t) - \hat{z}_2(t)\hat{z}_3(t)], \chi(t)\right)} \quad (20)$$
with $\chi(t) = \epsilon(\rho_m + \alpha_m) \frac{\beta_m}{\alpha_M N} y_1(t)$,

are estimators of the quantities $\beta \frac{S(t)}{N}$ and $\mathcal{R}_0 \frac{S(t)}{N}$, respectively, for $t \in [0, T_\epsilon]$.

**Proof.** Recall first that $\epsilon > 0$ is a lower bound of the quantity

$$-z_3(t) = \beta^2 \frac{S(t)I(t)}{N^2}, \quad t \in [0, T_\epsilon].$$

From the convergence of the vector $E(t)$ defined in (9), the quantity $\hat{z}_4(t) - \hat{z}_2(t)\hat{z}_3(t)$ is an estimation of

$$\chi_1(t) = \beta^3 \frac{S(t)I(t)^2}{N^3} > \epsilon \frac{\beta_m}{\alpha_M N} y_1(t), \quad t \in [0, T_\epsilon]$$
and $\hat{z}_3(t)^2$ is an estimation of

$$ \nu(t) = \beta^4 \frac{S(t)^2 I(t)^2}{N^4}. $$

Then, the expression $\beta^\dagger(t)$ given in (19) is an estimation of $\nu(t)/\chi_1(t) = \beta \frac{S(t)}{N}$ for $t \in [0, T_\epsilon]$.

Consider now the quantity $\hat{z}_3(t)^2 - \hat{z}_2(t)(\hat{z}_4(t) - \hat{z}_2(t)\hat{z}_3(t))$ which is an estimation, for $t \in [0, T_\epsilon]$, of

$$ \chi_2(t) = (\rho + \alpha)\beta^3 \frac{S(t)I(t)^2}{N^3} > \epsilon(\rho_m + \alpha_m) \frac{\beta_m}{\alpha_M N} y_1(t). $$

Then, the expression $\mathcal{R}_0^\dagger(t)$ given by (20) is an estimation of $\nu(t)/\chi_2(t) = \mathcal{R}_0 \frac{S(t)}{N}$ for $t \in [0, T_\epsilon]$.

**Remark 7.** Consider the dynamics of the ratio of susceptible $s(t) = S(t)/N$:

$$ \dot{s} = -\beta s \frac{I(t)}{N} = -\frac{y_1(t)}{N} \frac{\beta}{\alpha} s $$

and assume that $N$ is large. When $I(0) \ll N$ and the ratio $s(0)$ is close to 1, expressions (19) and (20) are estimators of the parameters $\beta$ and $\mathcal{R}_0$, as long as the observation $y_1(t)$ is small compared to $N$.

**Corollary 8.** Assume that the parameter $\rho$ is known. Under Assumptions of Proposition 5, the formulas (3), (4), (5), and (6) provide estimators of parameters $\alpha, \beta$ and variables $S(t), I(t)$, respectively, for $t \in [0, T_\epsilon]$.

**Proof.** From the convergence of the vector $E$ given by expression (9) in Proposition 5, $\hat{z}_2(t)$ is an estimation of the quantity

$$ \beta \frac{S(t)}{N} - \rho - \alpha $$

and from the estimation (19) of the quantity $\beta S(t)/N$ provided by Corollary 6, one obtains the estimation (3) for the parameter $\alpha$. Again, from the expression (9) of vector $E(t)$, one obtains that $-\hat{z}_3(t)$ is an estimation of the quantity

$$ \beta \left( \frac{\beta S(t)}{N} \right) \frac{e^{\hat{z}_1(t)}}{\alpha N} $$

or equivalently

$$ Ne^{-\hat{z}_1(t)} \alpha \left( \frac{\beta S(t)}{N} \right)^{-1} $$

is an estimation of $\beta$. From the estimations (19) for $\beta S(t)/N$ and (3) for $\alpha$, one obtains the estimation (4) of the parameter $\beta$. Finally, the estimation of $S(t)$ is obtained from $\hat{z}_2(t)$ and the former estimation of $\beta$, and the estimation of $I(t)$ from $\hat{z}_1(t)$ and the estimation of $\alpha$.

When $\rho$ is known (or estimated with the observation $y_2$) and $N$ is large with $S(0)/N$ close to 1, we show now that one can consider a reduced order observer as long as $S(t)/N$ remains close to 1. For this purpose, we introduce the new parameter

$$ k = \frac{\beta^2}{\alpha}, \qquad (21) $$

and we consider the variable

$$ v(t) = \beta \frac{S(t)}{N} - \rho - \alpha \qquad (22) $$

One can write $\dot{v} = -\frac{\beta^2}{\alpha} \frac{S(t)}{N^2} y_1(t) \simeq -\frac{\beta^2}{\alpha N} y_1(t) = -\frac{k}{N} y_1(t)$ (this last approximation being valid as long as the quantity $1 - \frac{S(t)}{N}$ is small compared to 1).

Then, we have

$$ \dot{y}_1 = vy_1 \qquad (23) $$

$$ \dot{v} = -\frac{k}{N} y_1 \qquad (24) $$

$$ \dot{k} = 0 \qquad (25) $$

We look now for an observer for the system (23)-(24)-(25) with the observation $y_1$.

**Proposition 9.** Consider positive numbers $N, y_{10}, \bar{v}_0, \bar{k}$. For $\theta > 0$ large enough, the dynamical system

$$ \dot{z} = \begin{bmatrix} \hat{z}_2 \\ \hat{z}_3 \\ \tilde{z}_2 \tilde{z}_3 \end{bmatrix} - \begin{bmatrix} 3\theta \\ 3\theta^2 \\ \theta^3 \end{bmatrix} (\hat{z}_1 - \log(y_1(t))), \quad \hat{z}(0) = \hat{z}_0 \qquad (26) $$

where

$$ \tilde{z}_2 = sat \frac{\bar{v}_0}{\sqrt{\bar{v}_0^2 + \frac{2\bar{k}}{N} y_{10}}} (\hat{z}_2), \quad \tilde{z}_3 = sat \frac{0}{\frac{\bar{v}_0^2}{2} - \frac{\bar{k}}{N} y_{10}} (\hat{z}_3) \qquad (27) $$

is an exponential observer for the system (23)-(24)-(25) with observation $y_1(\cdot)$ and initial condition $(y_{10}, v_0, k)$ such that $v_0 \le \bar{v}_0, k \ge \bar{k}$, in the sense that there exists $K > 0$ such that the vector

$$ E(t) = \hat{z}(t) - \begin{bmatrix} \log(y_1(t)) \\ v(t) \\ -\frac{k}{N} y_1(t) \end{bmatrix} $$

verifies

$$ \|E(t)\| \le K \|E(0)\| \exp \left( -\frac{\theta}{3} t \right), \quad t \ge 0 $$

**Proof.** Clearly, the solution $y_1(t)$ remains positive for any $t \in \mathbb{R}_+$. Let $z_1(t) = \log(y_1(t))$ and consider the variables $z_2(t) = \dot{z}_1(t)$ and $z_3(t) = \dot{z}_2(t)$. From (23)-(24)-(25), one gets $z_2(t) = v(t), z_3(t) = -\frac{k}{N} y_1(t)$ and obtains $\dot{z}_3 = -\frac{K}{N} y_1 v = z_2 z_3 := \phi(z)$. Note that the function

$$ V(y_1, v) = y_1 + \frac{N}{2k} v^2 $$

is constant along the solutions of (23)-(24)-(25) i.e.

$$ \frac{d}{dt} V(y_1(t), v(t)) = v(t)y_1(t) + \frac{N}{2k} 2v(t) \frac{k}{N} y_1(t) = 0 $$

from which one deduces that its solutions are bounded with

$$ z_2(t) \in \left[ -\sqrt{\frac{\bar{v}_0^2 + \frac{2\bar{k}}{N} y_{10}}{\bar{v}_0}}, z_3(t) \in \left[ -\frac{\bar{v}_0^2}{2} - \frac{\bar{k}}{N} y_{10}, 0 \right] $$
So one can write

$$ \dot{z}_3 = \dot{\phi}(z) := \phi(\sigma(z)) \text{ with } \sigma(z) = \begin{bmatrix} z_1 \\ \frac{\mathrm{sat}\,\bar{v}_0}{\sqrt{\bar{v}_0^2 + \frac{2\bar{k}}{N}y_{10}}} \\ \frac{\mathrm{sat}\,\bar{v}_0^0}{-\frac{\bar{v}_0^2}{2} - \frac{\bar{k}}{N}y_{10}} \\ z_3 \end{bmatrix} $$

Therefore, the solutions of system (23)-(24)-(25) with any initial condition $(y_{10}, v_0, k)$ such that $v_0 \le \bar{v}_0$, $k \le \bar{k}$ written in in coordinates $z = (z_1, z_2, z_3)$ are solutions of the system

$$ \dot{z} = F(z) := \underbrace{\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}}_{A_3} z + \underbrace{\begin{bmatrix} 0 \\ 0 \\ \tilde{\phi}(z) \end{bmatrix}}_{\tilde{y} = \underbrace{\begin{bmatrix} 1 & 0 & 0 \end{bmatrix}}_{C_3} z} , \quad \tilde{y} = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} z $$

which is in the "canonical observable form" Gauthier and Kupka (2001) for the output $\tilde{y}$, and where $\tilde{\phi}$ is globally Lipschitz on $\mathbb{R}^3$. We are then in position to state that the dynamical system

$$ \dot{\tilde{z}} = F(\tilde{z}) - \Sigma_{\theta}^{-1}(C\tilde{z} - \tilde{y}), \quad \tilde{z}(0) = \tilde{z}_0 $$

is a globally exponential observer on $\mathbb{R}^3$ for $\theta > 0$ large enough, where the symmetric matrix $S_\theta$ is solution of the Lyapunov equation

$$ \theta\Sigma_{\theta} + A_{3}^{\top}\Sigma_{\theta} + \Sigma_{\theta}A_{3} - C_{3}^{\top}C_{3} = 0 \quad (28) $$

and such that

$$ \|\hat{z}(t) - z(t)\| \le K \|\hat{z}_0 - z(0)\| \exp\left(-\frac{\theta t}{3}\right), \quad t \ge 0 $$

for some $K > 0$ (see Gauthier et al. (1992); Gauthier and Kupka (2001)). Finally, one can easily check that the inverse of the solution $\Sigma_\theta$ of (28) is given by the following expression

$$ \Sigma_{\theta}^{-1} = \begin{bmatrix} 3\theta & 3\theta^2 & \theta^3 \\ 3\theta^2 & 5\theta^3 & 2\theta^4 \\ \theta^3 & 2\theta^4 & \theta^5 \end{bmatrix} $$

which gives the equations of the observer (26).

**Corollary 10.** Let Assumptions 1 and 3 be fulfilled and assume that the parameter $\rho$ is known. When $N$ is large and $S(0)/N$ is close to 1 with $I(0) > 0$, then the observer (26) defined in Proposition 9 with $\bar{v}_0 = \beta_M - \rho - \alpha_m$ and $\bar{k} = \beta_M^2/\alpha_m$ provides the following estimators, as long as $S(t)/N$ remains close to 1.

$$ \hat{\beta}(t) = \frac{-\hat{z}_3(t)N}{2e^{\hat{z}_1(t)}} \sqrt{\frac{\max((\hat{z}_3(t)N)^2 + 4\hat{z}_3(t)N(\rho + \hat{z}_2(t)))e^{\hat{z}_1(t)}, 0}} \quad (29) $$

$$ \hat{\alpha}(t) = \hat{\beta}(t) - \rho - \hat{z}_2(t) \quad (30) $$

**Proof.** From the expressions (22), (21) of $v(t)$ and $k$, one obtains that the parameter $\beta$ is solution of the equation

$$ \beta^2 - k \frac{S(t)}{N} \beta + k(v(t) + \rho) = 0, \quad t \ge 0 $$

which has necessary real roots $\beta^-$, $\beta^+$ with $\beta^- \le k \frac{S(t)}{N} \le \beta^+$. However, having $\beta = \beta^+$ with $S(t)/N$ close to one implies $\beta - \rho - \alpha < 0$ which contradicts the hypothesis $\mathcal{R}_0 > 1$. So, $\beta$ has to be the smallest root. With $z_2(t) = v(t)$ and $z_3(t) = -\frac{k}{N}y_1(t)$, one can then write

$$ \beta(t) = \frac{-z_3(t)S(t)}{2y_1(t)} \sqrt{\frac{\max((z_3(t)S(t))^2 + 4z_3(t)(\rho + z_2(t)))y_1(t), 0}} \quad (29) $$ $$ := B(y_1(t), z_2(t), z_3(t), S(t)) $$

Note that the function $B$ is well defined on $\mathbb{R}_+^* \times \mathbb{R}^2 \times \mathbb{R}^+$. Then, from the convergence of the error vector $Es$ provided by Proposition 9, and as long as $S(t)/N$ remains close to 1, one obtains the estimator of $\beta$

$$ \hat{\beta}(t) = B(e^{\hat{z}_1(t)}, \hat{z}_2(t), \hat{z}_3(t)), N) $$

which gives the expression (29). Finally, an estimator of $\alpha$ is straightforwardly given by equation (30).

**Numerical simulations of the High-Gain estimators:** We illustrate the convergence of the observer (Proposition 9) and the estimations (Corollary 8) with the following values

<table><thead><tr><th>α</th><th>β</th><th>ρ</th><th>N</th><th>I(0)</th><th>R(0)</th><th>Q(0)</th></tr></thead><tbody><tr><td>0.1</td><td>0.4</td><td>0.1</td><td>10<sup>5</sup></td><td>1000</td><td>0</td><td>0</td></tr></tbody></table>

under the hypothesis that $\rho$ is known. On Figure 1, one can see on the first picture the trajectories of the system in dash-line and the estimation of the state variables provided by the observer in plain line. The second picture shows the estimations of the parameters $\alpha$ and $\beta$. The last four pictures give the trajectories of the four internal states $\hat{z}_i$ of the observer along with the trajectories of the system in $z$ coordinates.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Simulations with observer gain $\theta = 2.5$</figcaption></figure>

Simulations of the reduced order observer given in Proposition 9 with the estimations of $\alpha$ and $\beta$ provided by Corollary 10 are presented on Figure 2 for the same parameters but with the initial condition $I(0) = 10$.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Simulations of the reduced order observer with gain $\theta = 2.5$</figcaption></figure>

## 3.2 Using "Simpler" observers

We present a simpler estimator that is, moreover, more robust with respect to noise. We consider again the parameter (defined in (21)) $k = \frac{\beta^2}{\alpha}$.

**Proposition 11.** Let $\lambda$ and $\mu$ be two positive vectors (i.e., with positive components) in $\mathbb{R}^4$ and $\mathbb{R}^3$ respectively. For $t > 0$, consider the dynamical system

$$ 
\left\{ 
\begin{aligned}
\dot{\hat{z}}_1 &= \hat{\delta} - \hat{\rho} - K_1(\hat{z}_1 - \log(y_1(t))), \\
\dot{\hat{z}}_2 &= \frac{y_1(t)}{y_2(t)} - \hat{\rho} - K_2(\hat{z}_1 - \log(y_1(t))), \\
\dot{\hat{\delta}} &= -K_3(\hat{z}_1 - \log(y_1(t))), \\
\dot{\hat{\rho}} &= -K_4(\hat{z}_1 - \log(y_1(t))) - (\hat{z}_2 - \log(y_2(t))), \\
\dot{\hat{y}}_1 &= \hat{v}y_1(t) - K_5y_1(t)(\hat{y}_1 - y_1(t)), \\
\dot{\hat{v}} &= -\hat{k}\frac{y_1(t)}{N} - K_6y_1(t)(\hat{y}_1 - y_1(t)), \\
\dot{\hat{k}} &= -K_7Ny_1(t)(\hat{y}_1 - y_1(t)),
\end{aligned}
\right.
\quad (31) $$

with the gain vector

$$ K = \left[ \begin{array}{c} \displaystyle \sum_{i=1}^{4} \lambda_i \\ \displaystyle \sum_{i=1}^{4} \lambda_i + \sum_{1 \le i < j < k \le 4} \lambda_i \lambda_j \lambda_k \\ \displaystyle - \lambda_1 \lambda_2 \lambda_3 \lambda_4 \\ \displaystyle - \sum_{1 \le i < j \le 4} \lambda_i \lambda_j - \lambda_1 \lambda_2 \lambda_3 \lambda_4 - 1 \\ \displaystyle \mu_1 + \mu_2 + \mu_3 \\ \displaystyle \mu_1 \mu_2 + \mu_1 \mu_3 + \mu_2 \mu_3 \\ \displaystyle - \mu_1 \mu_2 \mu_3 \end{array} \right]. \quad (32) $$

Then, the output vector

$$ \hat{p}(t) = \left[ \begin{array}{c} \hat{\rho}(t) \\ \hat{\alpha}(t) = \hat{\beta}(t) - \hat{\delta}(t) \end{array} \right] \quad (33) $$

is an estimator of $p = [\rho, \beta, \alpha]^\top$, that is, $\hat{p}(t)$ converges exponentially towards $p$ as $t$ tends to infinity. Moreover, the exponential decay of the error can be made as fast as desired keeping the number

$$ l = \min_{i,j} (\lambda_i, \mu_j) $$

large, as long as $S(t)/N$ remains close to 1. Moreover, the state $I$ is estimated with

$$ \hat{I}(t) = \frac{y_1(t)}{\hat{\alpha}(t)}. \quad (34) $$

**Proof.** We define a new parameter $\delta$ by

$$ \delta := \beta - \alpha. \quad (35) $$

As far as $S/N$ remains close to 1, the size of the population $I$ is small compared to $S$, and the dynamics of the outputs $y_1 = \alpha I$ and $y_2 = Q$ can be approximated by the linear dynamics

$$ \left\{ \begin{aligned} \dot{y}_1 &= \delta y_1 - \rho y_1 \\ \dot{y}_2 &= y_1 - \rho y_2 \end{aligned} \right. \quad (36) $$

For $t > 0$ and $i = 1, 2$, consider the new outputs $z_i(t) = \log(y_i(t))$, whose dynamics is given by the system

$$ \left\{ \begin{aligned} \dot{z}_1 &= \delta - \rho \\ \dot{z}_2 &= \exp(z_1 - z_2) - \rho \end{aligned} \right. \quad (37) $$

For this sub-system with unknown parameters $\delta$ and $\rho$, we consider the following candidate observer in $\mathbb{R}^4$:

$$ \left\{ \begin{aligned} \dot{\hat{z}}_1 &= \hat{\delta} - \hat{\rho} - K_1(\hat{z}_1 - z_1), \\ \dot{\hat{z}}_2 &= \exp(z_1 - z_2) - \hat{\rho} - K_2(\hat{z}_1 - z_1), \\ \dot{\hat{\delta}} &= -K_3(\hat{z}_1 - z_1), \\ \dot{\hat{\rho}} &= -K_4(\hat{z}_1 - z_1) - (\hat{z}_2 - z_2), \end{aligned} \right. \quad (38) $$

The dynamics of the error $e_1 = (\hat{z}_1, \hat{z}_2, \hat{\delta}, \hat{\rho}) - (z_1, z_2, \delta, \rho)$ is given by the linear time invariant system $\dot{e}_1 = M_1 e_1$ with

$$ M_1 = \begin{pmatrix} -K_1 & 0 & 1 & -1 \\ -K_2 & 0 & 0 & -1 \\ -K_3 & 0 & 0 & 0 \\ -K_4 & -1 & 0 & 0 \end{pmatrix} $$

A calculation shows that with the choice of $K_i$ ($i = 1, ..4$) given by formula (32), the spectrum of $M_1$ is $\{-\lambda_i, i = 1 \cdots 4\}$. This shows that the first four equations of system (31) give the reconstruction of the parameters $\delta$ and $\rho$ with an exponential decay of the error larger than $\min_i \lambda_i$.
Consider the variable

$$v(t) := \beta \frac{S(t)}{N} - \rho - \alpha.$$

As far as $S(t)/N$ remains close to 1, one can write the approximation

$$\dot{v} = -\frac{\beta^2}{\alpha} \frac{S(t)}{N^2} y_1(t) \approx -\frac{\beta^2}{\alpha N} y_1(t)$$

Then, this amounts to approximate the dynamics of (2) or (1) in the $(y_1, v, k)$ coordinates by the following dynamical system

$$\left\{ \begin{aligned} \dot{y}_1 &= vy_1 \\ \dot{v} &= -\frac{k}{N}y_1, \quad \dot{k} = 0, \end{aligned} \right.$$

where $k$ is an unknown parameter, for which we consider the following candidate observer in $\mathbb{R}^3$

$$\left\{ \begin{aligned} \dot{\hat{y}}_1 &= \hat{v}y_1 - K_5y_1(\hat{y}_1 - y_1) \\ \dot{\hat{v}} &= -\hat{k}\frac{y_1}{N} - K_6y_1(\hat{y}_1 - y_1) \\ \dot{\hat{k}} &= -K_7Ny_1(\hat{y}_1 - y_1) \end{aligned} \right. \quad (39)$$

whose dynamics of the error $e_2 = (\hat{y}_1, \hat{v}, \hat{k})^\top - (y_1, v, k)^\top$ is given by the non-autonomous linear system

$$\dot{e}_2 = y_1(t)M_2 e_2 \quad (40)$$

with

$$M_2 = \begin{pmatrix} -K_5 & 1 & 0 \\ -K_6 & 0 & -\frac{1}{N} \\ -K_7N & 0 & 0 \end{pmatrix}$$

One can easily check that for the choice of $K_5, K_6$ and $K_7$ given by formula (32), the spectrum of $M_2$ is $\{-\mu_i, i = 1, \dots, 3\}$. Then, from (40), we obtain the upper bound on the error decrease

$$|\hat{k}(t) - k| \le \exp \left( -\left( \min_j \mu_j \right) \int_0^T y_1(\tau) d\tau \right) ||e_2(0)||$$

whose exponential decay can be made as large as desired with large $l$.

Finally, from the reconstruction of parameters $\delta, k$ by observers (38), (39) and expressions (35) and (21), the original parameters $\alpha, \beta$ are recovered as roots of

$$\beta^2 - k\beta + k\delta = 0 \Rightarrow \beta = \frac{k \pm \sqrt{k^2 - 4k\delta}}{2} \quad (41)$$

Note first that Assumption 1 implies $\beta > \alpha$ and thus $k^2 - 4k\delta = k^2 - 4k(\beta - \alpha) > 0$. Moreover, one has $k > \beta(1 + \frac{\rho}{\alpha})$, and by Assumption 2 one has $k > 2\beta$, which implies that only the smaller root of (41) is valid, leading to the expression (33) of the estimator. Note that this expression preserves the exponential decay of the error obtained for $\delta, k$ and $\rho$.

Let us make some comments about this observer. It consists in reconstructing functions of the parameters $\delta$ and $k$ and not directly the parameters $\alpha, \beta$. There is an apparent redundancy of variables $\hat{z}_1$ and $\hat{y}_1$ in dynamics (31), which reconstruct $\log y_1$ and $y_1$. Indeed, this allows to decouple the observer into two sub-systems of dimensions 4 and 3, which avoids the use of two large correction gains compared to a full order observer. Finally, outputs of these two sub-systems are coupled in expression (33) to reconstruct the original parameters.

**Numerical illustrations:** The proposed observer has been tested with synthetic data (simulated data) for a population size $N = 10^5$ with parameter values $\alpha = 0.07, \beta = 0.4, \rho = 0.1$ and initial condition $I(0) = 10, Q(0) = 5, R(0) = 0$ over a time horizon of 10 days (see Figure 3). The gains have been computed for the choice of vectors $\lambda = [1; 1.5; 2; 2.5]$ and $\mu = [1/(13.10^3); 1/(15.10^3); 1/(19.10^3)]$. Note that vector $\mu$ has been chosen quite small to avoid too large gains when multiplied by $N$ in the observer equations. Then, we have simulated a measurement noise

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Observer simulation without measurement noise.</figcaption></figure>
with a centered Gaussian law of variance equal to 5% of the signal (see Figure 4). Because the expression (34) of the

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 4. Observer simulation with measurement noise.</figcaption></figure>

estimation of the state *I* uses the unfiltered measurement *y*₁, we have applied a moving average smoothing to the estimation *Î* (see Figure 5). Finally, these simulations show that the method allows to reconstruct the parameter values in few days in a quite accurate manner.

## 4. CONCLUSION

This work shows that although the identifiability of the SIR-Q models has singularity points where measured variables are null, it is possible to design an observer with exponential decay of the estimation error during the first stage of the epidemics, and recover parameters in few days. Further investigations will concern real data of COVID epidemics provided by various territories.

## REFERENCES

* Abdelhedi, A., Boutat, D., Sbita, L., and Tami, R. (2014). Extended observer to estimate the spreading of conta-

<figure><img src="image_12.png" /><figcaption>Fig. 5. Smoothed estimation of the infected population in presence of measurement noise.</figcaption></figure>

* gious disease. In *2014 European Control Conference (ECC)*, 1843–1848. doi:10.1109/ECC.2014.6862172.

* Bichara, D., Cozic, N., and Iggidr, A. (2014). On the estimation of sequestered infected erythrocytes in *plasmodium falciparum* malaria patients. *Math. Biosci. Eng.*, 11(4), 741–759.

* Bichara, D., Guiro, A., Iggidr, A., and Ngom, D. (2019). State and parameter estimation for a class of schistosomiasis models. *Mathematical Biosciences*, 315.

* Bliman, P.A., Efimov, D., and Ushirobira, R. (2018). A class of nonlinear adaptive observers for sir epidemic model. In *ECC 2018-European Control Conference*, 6.

* Degue, K.H., Efimov, D., and Iggidr, A. (2021). Interval Observer Design for Sequestered Erythrocytes Concentration Estimation in Severe Malaria Patients. *European Journal of Control*, 58, 399–407.

* Diaby, M., Iggidr, A., and Sy, M. (2015). Observer design for a schistosomiasis model. *Math. Biosci.*, 269, 17–29.

* Diouf, M.L., Iggidr, A., and Souza, M.O. (2019). Stability and estimation problems related to a stage-structured epidemic model. *Mathematical Biosciences and Engineering*, 16(5), 4415–4432.

* Gauthier, J.P., Hammouri, H., and Othman, S. (1992). A simple observer for nonlinear systems applications to bioreactors. *IEEE Trans. Autom. Control*, 37(6), 875–880.

* Gauthier, J.P. and Kupka, I. (2001). *Deterministic Observation Theory and Applications*. Cambridge University Press. doi:10.1017/cbo9780511546648.

* Hamelin, F., Iggidr, A., Rapaport, A., and Sallet, G. (2020). Observability, identifiability and epidemiology of a survey. *arXiv preprint arXiv:2011.12202*.

* Hethcote, H., Zhien, M., and Shengbing, L. (2002). Effects of quarantine in six endemic models for infectious diseases. *Math. Biosci.*, 180, 141–160. John A. Jacquez memorial volume.

* Iggidr, A. and Souza, M.O. (2019). State estimators for some epidemiological systems. *Journal of Mathematical Biology*, 78, 225–256.

* Li, M.Y. (2018). *An introduction to mathematical modeling of infectious diseases*, volume 2 of *Mathematics of Planet Earth*. Springer, Cham.

* Nuño, M., Castillo-Chavez, C., Feng, Z., and Martcheva, M. (2008). Mathematical models of influenza: the role of cross-immunity, quarantine and age-structure. In *Mathematical epidemiology*, volume 1945 of *Lecture Notes in Math.*, 349–364. Springer, Berlin.