State Observer for Time-Delay Systems Applied a Sir Compartmental Epidemiological Model for the COVID-19

1. Introduction

Since its origins, humanity has been immersed in the emergence of infectious diseases and/or its variants, bringing with it a high degree of confusion among the population due to the ignorance of its spread in a susceptible population. The use of epidemiological mathematical models has become a highly recurrent task today for researchers and/or health organizations, since they are a tool that allows describing the dynamics of an infectious disease and knowing some parameters involved in its spread. Among the best-known models for this purpose are the compartmental epidemiological mathematical models which differ according to the number of compartments, populations or subpopulations, considered in each model, as SIS (Susceptible-Infectious), SIR (Susceptible-Infectious-Recovered), and SEIR (Susceptible-Exposed-Infectious-Recovered), to mention a few, [1].

Historically, various infectious diseases have emerged that have seriously affected humanity, such as the so-called spanish flu; black death plague; HIV; among others. Table 1 shows some of the most relevant data on these infectious diseases. In addition to a small compilation of the mathematical models type SIS, SIR and/or SEIR used to study these pandemics.

Through a mathematical model, it is possible to improve the understanding of a disease and lower its prevalence and incidence, and based on the behavior of the model, it is possible to assist in the design of strategies to stabilize the transmission rate of the disease and, if possible, stop the spread of an infectious disease [11]. However, there is not always available information for all the population compartments considered in these models, so it is necessary to use tools of the control theory, such as state observers, to estimate the missing information. It is well known that the main objective of the synthesis of a state observer is the complete reconstruction of the state of the system (mathematical model) using the available information and initial conditions [12], so in the literature it is common to read about the use of this tool for a wide variety of systems.

Throughout history, various preconceptions have been proposed to estimate the state of a system, e.g., [13], but in [14] the term observer is proposed for the first time. This observer is known Luenberger observer. Here it is shown that the state vector of a linear system can be reconstructed by observing the inputs and outputs of the system.

While on observers for time-delay systems, in [15] an extension of the results given in [13] is presented for systems with time delays. A new approach for the construction of a state observer for nonlinear systems when the output measurements are available for computations after a non negligible time delay is proposed in [16]. There, the globally drift-observability is secured by a diffeomorphism using Lie derivatives. In [17] an extension of the observer given in [18] is proposed for a class of nonlinear systems, which present multiple noncommensurate time delays as well as distributed delay terms. In [19] a state observer for drift observable nonlinear systems when output measurements are affected by a known and bounded time varying delay is given. Such observer has a simple structure, it is a generalization of an observer for delay-free systems, the globally drift-observability is secured by a diffeomorphism using Lie derivatives and the convergence of the observation error is based on the Lyapunov–Razumikhin approach. In [20] a combined controller-observer design procedure for a class of time-delay uncertain systems is proposed using matrix’s Singular Value Decomposition (SVD) and Linear Matrix Inequality (LMI) techniques. The stability of the closed-loop systems presented is based on Lyapunov-Krasovskii functionals. In [21] a new method to design observers in a class of nonlinear time-delay systems is suggested. This design procedure of the observer involves two steps: an extension of the pseudo-linearization technique and solving a state-dependent Riccati equation. In [22] is presented a robust observer-based control of nonlinear one-sided Lipschitz systems in the presence of time delay, model uncertainties, and unknown energy-bounded exogenous disturbances. A observer structure requires a predetermined nominal/linear part. The observer-based finite-time robust control problem of a general class of nonlinear time-delay systems via Hamiltonian function method, that proposes some new finite-time robust stabilization results is studied in [23].

Derived from the current health crisis, the scientific community has focused efforts to understand, predict, and mitigate COVID-19. For this, various investigations have been proposed in all areas of science and technology. One of these is the construction of observers for epidemiological models, although this is not new, since in [24] the state estimation for a delay-free SIR model subject to seasonal variations and uncertainties in the transmission rate is considered. Here an asymptotic convergence of the observer is ensured using Lyapunov functions. While in [25] the same approach is presented but now for a nonlinear adaptive observer. The construction of backstepping observers for two linearized Kermack–McKendrick models: a first-order hyperbolic partial differential equation and a SIR model, both delay free models is proposed in [26]. In [27] a SIR model with two time-delays: $h_1$ in the transmission rate and $h_2$ in the number of infectious population is considered. Here an observer for the delay-free SIR model which is then uplifted into a predictor to compensate for time delays in the input and output is proposed. In [28] a high-gain observer design for a class of quasi-linear hyperbolic systems is presented. The exponential stability of the observation error is ensured using a Lyapunov-based analysis.

This manuscript presents the design of a Luenberger-type state observer for a class of nonlinear systems with multiple delays to estimate compartmental populations of epidemiological models. As well as an SIR compartmental epidemiological model with three time-delay. The model considers that a part of the recovered population becomes susceptible again after a period of time after its recovery plus three time-delays due to the exchange of individuals between the population compartments: ${\tau }_{1,2,3}$, for the dead-times of recovery, immunity loss and incubation, respectively. By making a change of variables, this model with three time-delays is transformed into a model with nominal part as a function of the time-delays. The latter allows to design a Luenberger-type state observer for estimation of the compartmental populations of $S\left(t\right)$Susceptible and $R\left(t\right)$Recovered not reported by the WHO. This only with the availability of real data from the Infectious compartmental population. In addition, using Lyapunov-Krasovskii functionals and the feasibility of Linear Matrix Inequalities, an exponential decay can be ensured on the observation error, thus guaranteeing a convergence between the real data and the estimated. To illustrate the efficiency and performance of the proposed results, without loss of generality real data reported by WHO is used to estimate $S\left(t\right)$ and $R\left(t\right)$ in a 55-day time window with the greatest impact on the population of Mexico. The data considered in the present contribution were collected prior to the implementation of the vaccination campaign. This is mainly to the fact that the SIR compartmental epidemiological model under consideration does not include this variable.

2. Problem Statement and Contributions

Next, the statement of the problem and specific contributions are given.

2.1. Problem Statement

Consider the following SIR fundamental model of Kermack and McKendrick (KMK) used in mathematical epidemiology [29]:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}S\left(t\right)=& -\beta S\left(t\right)I\left(t\right),\hfill \\ \hfill {\displaystyle \frac{d}{dt}}I\left(t\right)=& \beta S\left(t\right)I\left(t\right)-\gamma I\left(t\right),\hfill \\ \hfill {\displaystyle \frac{d}{dt}}R\left(t\right)=& \gamma I\left(t\right),\hfill \end{array}$$

(1)

where $S\left(t\right)$, $I\left(t\right)$ and $R\left(t\right)$ are the populations of Susceptible, Infectious and Recovered, respectively. The total population is denoted as $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$ and it is considered constant. Also, $\beta$ denotes the transmission coefficient or infection rate (recording infectiousness and the probability of contact between people or number of people infected per unit of time) and $\gamma$ the recovery/death rate in the infected population, this can be estimated as the reciprocal of the average number of days for transition from $I\left(t\right)$ to $R\left(t\right)$, see [30]. In addition, it is considered that if the disease does not induce mortality ($\gamma$ is the recovery rate), then the population $R\left(t\right)$ denotes those who have recovered and are immune to reinfection. On the other hand, if $\gamma$ records mortality ($\gamma$ is the dead rate), then the population $R\left(t\right)$ records the number of deaths induced by the epidemic, see [31].

Epidemiological models such as SIR can be used to better understand what has happened to the current pandemic, such as rates of infection, recovery or/and spread; estimates of population dynamics, among others, this during short periods/windows of time. However, according to WHO reports and statistics, it is necessary to adapt the fundamental model of KMK to the conditions given by the current pandemic [32], such as:

I. 

Take into account that recovered individuals may again be susceptible to the disease (there is not immunity).

II. 

Consider the dead-times/periods it takes an individual to move from one compartmental population to another.

In addition, the lack of information in the reported data limits the study of this disease and its variants, therefore it is necessary:

III. 

Estimate the dynamics/data of compartmental populations not reported by the authorities responsible for this matter.

Reason why in the present research work, the following is proposed:

2.2. Contributions of this Manuscript

• A Luenberger-type state observer for a class of nonlinear systems with multiple delays is proposed to estimate the dynamics of compartmental populations not reported by the WHO.
• A SIR epidemiological model with three time-delays is presented which considers that a part of the recovered population becomes susceptible again after a period of following recovery. Besides, the model consider the dead-times it takes an individual to move from one compartmental population to another, that is:

  -

  ${\tau }_1$ is the time it takes for an infectious individual to stop being so, in order to recover from the disease and acquire immunity (dead-time by recovery).

  -

  ${\tau }_2$ is the time it takes for a recovered individual to return to the susceptible population (dead-time by immunity loss).

  -

  ${\tau }_3$ is the time it takes for a susceptible individual to present symptoms and affect the spread of the disease (dead-time by incubation).

• To guarantee the estimation of the previous observer, a convergence analysis is presented using Lyapunov-Krasovskii functionals and the feasibility of Linear Matrix Inequalities to ensure that the observation error decreases with an exponential decay $\sigma >0$.
• The proposed results are validated through simulations (Simulink-Matlab), and data reported by the WHO for the COVID-19 pandemic in Mexico.
• The generality of the theoretical results proposed suggest that they can be applied to other types of epidemiological models that meet the conditions stated here; and to other countries (population data).

3. Preliminary Results and Observer Design

Next, some classic results found in the literature regarding time delay systems and state observers are stated; which will serve as the basis for the development of this proposal.

3.1. Systems with Several Delays

Let us consider a general retarded functional differential equations with multiple delays of the form

$$\begin{array}{cc}\hfill \dot{z}\left(t\right)& =f\left(z\left(t\right),z(t-{\tau }_1),z(t-{\tau }_2),\dots ,z(t-{\tau }_m)\right),\hfill \\ \hfill y\left(t\right)& =h\left(z\left(t\right),z(t-{\tau }_1),z(t-{\tau }_2),\dots ,z(t-{\tau }_m\right),z\left(\theta \right)\hfill & \hfill =\varphi \left(\theta \right),\theta \in [-\tau ,0],\end{array}$$

(2)

where $0={\tau }_0<{\tau }_1<\cdots <{\tau }_m=\tau$ are delays and $f:{\mathbb{R}}^{n\times (m+1)}\to {\mathbb{R}}^n$ is a continuous nonlinear function and satisfies a local Lipschitz condition regarding the elements of the argument and $h:{\mathbb{R}}^{n\times (m+1)}\to {\mathbb{R}}^l$ is continuously differentiable function with respect to its argument. For $t\ge 0$, $y\left(t\right)\in {\mathbb{R}}^l$ is the only measured output of the system, denote by $z(t,\varphi )\in {\mathbb{R}}^n$ the system solution with initial function (condition) $\varphi \in \mathfrak{C}$, and by $\mathfrak{C}:=C([-\tau ,0],{\mathbb{R}}^n)$ the Banach space with norm ${\parallel \varphi \parallel }_{\tau }:=\underset{\theta \in [-\tau ,0]}{max}\parallel \varphi \left(\theta \right)\parallel$. Here, $\parallel ·\parallel$ denotes the euclidean norm. As a natural extension of the initial function, a solution segment of $z(t,\varphi )$ in a time interval $[t-\tau ,t]$ is denoted by

$$z_t=z(t+\theta )\in \mathfrak{C},\theta \in [-\tau ,0];$$

and called state of the system (2).

Now, let $z_t^{\ast }$ is a non zero constant equilibrium point of (2), namely, $f(z^{\ast }\left(t\right),z^{\ast }(t-{\tau }_1),\dots$$,z^{\ast }(t-{\tau }_m))=0$, for all $t>0$. Using the variable transformation/change

$$x(t-{\tau }_i)=z(t-{\tau }_i)-z^{\ast }(t-{\tau }_i),i=0,1,\dots ,m,$$

(3)

the nonlinear system (2) can be rewritten in the form

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =f\left(x_t+z_t^{\ast }\right)=L\left(x_t\right)+g\left(x_t\right),\hfill \\ \hfill y\left(t\right)& =C\left(x_t\right);x\left(\theta \right)=\phi \left(\theta \right)=\varphi \left(\theta \right)-z^{\ast }\left(\theta \right),\theta \in [-\tau ,0],\hfill \end{array}$$

(4)

where $f:\mathfrak{C}\to {\mathbb{R}}^n$, and $g:\mathfrak{C}\to {\mathbb{R}}^n$ are nonlinear operators, and $C:\mathfrak{C}\to {\mathbb{R}}^n$ and $L:\mathfrak{C}\to {\mathbb{R}}^n$ are linear operator, all properly defined. Here, $L(·)$ is called the nominal system of the nonlinear system (2) and this is of the form

$$L\left(x_t\right)=\sum _{i=0}^m{A}_{i}x(t-{\tau }_i),\mathrm{and}C\left(x_t\right)=\sum _{i=0}^m{C}_{i}x(t-{\tau }_i),$$

(5)

with $A_i\in {\mathbb{R}}^{n\times n}$ and $C_i\in {\mathbb{R}}^{l\times n}$.

Remark 1. 

It is well-known that just as (2) was taken to (4), (3) can be used to take (4) to (2).

Theorem 1. 

[33] Consider the system given in (2) and let $u,v,w:{\overline{\mathbb{R}}}^{+}\to {\overline{\mathbb{R}}}^{+}$ be continuous nondecreasing functions, where additionally $u\left(s\right)$ and $v\left(s\right)$ are positive for $s>0$, and $u\left(0\right)=v\left(0\right)=0$. If there exists a continuous differentiable functional $V:{\mathbb{R}}^{+}\times \mathfrak{C}\to \mathbb{R}$ such that

• $u\left(\right|\left|\phi \left(0\right)\right|\left|\right)\le V(t,\phi )\le {v\left(\right|\left|\phi \right||}_{\tau })$, and
• $\dot{V}(t,\phi )\le -w\left(\right|\left|\phi \left(0\right)\right|\left|\right)$,

then the trivial solution of (2) is uniformly stable. If $w\left(s\right)>0$ for $s>0$, then it is uniformly asymptotically stable. If, in addition, $\underset{s\to \infty }{lim}u\left(s\right)=\infty$, then it is globally uniformly asymptotically stable.

Theorem 2. 

[34] Consider system given in (2), the trivial solution is exponential stable if there exist positive constants ${\alpha }_1$, ${\alpha }_2$ and σ; and a functional $V(t,\phi )$ such that satisfies the following two conditions:

• ${\alpha }_1{\left|\right|\phi \left(0\right)\left|\right|}^2\le V(t,\phi )\le {\alpha }_2{\left|\right|\phi \left|\right|}_{\tau }^2,t\ge 0$;
• $\dot{V}(t,x_t)+2\sigma V(t,x_t)\le 0,t\ge 0$.

3.2. Luenberger-Type State Observer for Nonlinear Time-Delay Systems

Consider a system of the form (4)-(5):

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\sum _{i=0}^m{A}_{i}x(t-{\tau }_i)+g\left(x_t\right),\hfill \\ \hfill y\left(t\right)& =\sum _{i=0}^m{C}_{i}x(t-{\tau }_i);x\left(\theta \right)=\phi \left(\theta \right),\theta \in [-\tau ,0].\hfill \end{array}$$

(6)

Thus, the dynamic Luenberger-type state observer of the system (6) is proposed as

$$\begin{array}{cc}\hfill \dot{\widehat{x}}\left(t\right)& =\sum _{i=0}^m{A}_i\widehat{x}(t-{\tau }_i)+\sum _{i=0}^m{L}_i\left[y\left(t\right)-\widehat{y}\left(t\right)\right]+g\left({\widehat{x}}_t\right),\hfill \end{array}$$

(7)

where $L_i\in {\mathbb{R}}^{n\times l}$, $\widehat{y}\left(t\right)=\sum _{i=0}^m{C}_i\widehat{x}(t-{\tau }_i)$. Defining the observation error by $e\left(t\right)=x\left(t\right)-\widehat{x}\left(t\right)$ and using the equations (6) and (7), the error dynamics can be expressed as

$$\dot{e}\left(t\right)=\sum _{i=0}^m\left[A_i-L_i{C}_i\right]e(t-{\tau }_i)+\eta (x_t,{\widehat{x}}_t).$$

(8)

Here $\eta (x_t,{\widehat{x}}_t)=g\left(x_t\right)-g\left({\widehat{x}}_t\right)$ is such that

$$\left|\right|\eta (x_t,{\widehat{x}}_t)\left|\right|\le M\left|\right|e_t{\left|\right|}_{\tau },M>0.$$

(9)

Remark 2. 

Due to the previous term $\eta (·)$, the uniformly exponentially convergent of the system response (8) cannot be guaranteed. However, the practical stability is refereed as uniformly exponentially convergent of the system response (8) to a ball $B_{\mu }$ with radius $\mu >0$ [35].

Next, sufficient conditions are proposed to guarantee practical stability on (8).

Proposition 1. 

If there exits $0<\sigma \in \mathbb{R}$, $R_0$, $R_i\in {\mathbb{R}}^{n\times l}$, $0<P$, $0<Q_i\in {\mathbb{R}}^{n\times n}$, with $i=1,\dots ,m$; such that the following linear matrix inequality (LMI)

$$\mathcal{M}(P,Q_1,\dots ,Q_m,R_0,\dots ,R_m)+2\sigma \mathcal{N}\left(P\right)<0,$$

(10)

where

$$\begin{array}{cc}\hfill & \mathcal{N}\left(P\right)=\mathrm{diag}\left\{P,0_n,\dots ,0_n\right\}\in {\mathbb{R}}^{n(m+1)\times (m+1)n},\hfill \\ \hfill & \mathcal{M}(P,Q_1,\dots ,Q_m,R_0,\dots ,R_m)=\mathit{diag}\left\{P+\sum _{i=1}^m{Q}_i,-{\mathsf{e}}^{-2\sigma {\tau }_1}{Q}_1,\dots ,-{\mathsf{e}}^{-2\sigma {\tau }_m}{Q}_m\right\}\hfill \\ \hfill & +{\left[P\mathrm{\Lambda }-\mathrm{\Gamma }\right]}^T{\rm Y}+{{\rm Y}}^T\left[P\mathrm{\Lambda }-\mathrm{\Gamma }\right]\in {\mathbb{R}}^{n(m+1)\times (m+1)n}.\hfill \end{array}$$

with $\mathrm{\Lambda }=\left[A_0{A}_1\dots A_m\right]\in {\mathbb{R}}^{n\times (m+1)n}$, $\mathrm{\Gamma }=\left[R_0{C}_0{R}_1{C}_1\dots R_m{C}_m\right]\in {\mathbb{R}}^{n\times (m+1)n}$, ${\rm Y}=\left[I_n{0}_n\dots 0_n\right]\in {\mathbb{R}}^{n\times (m+1)n}$, $I_n\in {\mathbb{R}}^{n\times n}$ identity matrix and $0_n\in {\mathbb{R}}^{n\times n}$ zero matrix. Then, the solution of the observation error (8) satisfies the following exponential decay

$$\left|\right|e(t,\phi )\left|\right|\le \left(\sqrt{{\displaystyle \frac{{\alpha }_2}{{\alpha }_1}}}{\left|\right|\phi \left|\right|}_{\tau }-\mu \right){\mathrm{e}}^{-\sigma t}+\mu ,\forall t\ge 0.$$

(11)

Furthermore

$$\left|\right|e(t,\phi )\left|\right|\le \mu ,\forall t\ge T(\mu ,\phi ),$$

where $\mu ={\displaystyle \frac{\gamma {\alpha }_{3}M}{\sigma {\alpha }_1}}$, with ${\alpha }_3=\left|\right|P\left|\right|:={\lambda }_{max}\left(P\right)$ is the maximum eigenvalue of P, and $T(\mu ,\varphi )={\displaystyle \frac{1}{\sigma }}ln\left(\sqrt{{\displaystyle \frac{{\alpha }_2}{{\alpha }_1}}}{\left|\right|\phi \left|\right|}_{\tau }-\mu \right)$. The observer gains in (21) are $L_i=P^{-1}{R}_i\in {\mathbb{R}}^{n\times l}$, $i=0,\dots ,m$.

Proof. 

Consider the following candidate Lyapunov-Krasovskii functional [36]:

$$v\left(e_t\right)=e^T\left(t\right)Pe\left(t\right)+\sum _{i=1}^m{\int }_{-{\tau }_i}^0{e}^T(t+\theta ){\mathsf{e}}^{2\sigma \theta }{Q}_{i}e(t+\theta )d\theta .$$

(12)

Note that

$${\alpha }_1{\left|\right|e\left(t\right)\left|\right|}^2\le v\left(e_t\right)\le {\alpha }_2\left|\right|e_t{\left|\right|}_{\tau }^2,$$

(13)

where ${\alpha }_1={\lambda }_{min}\left(P\right)$ is the minimum eigenvalue of P, and ${\alpha }_2={\lambda }_{max}\left(P\right)+{\displaystyle \sum _{i=1}^m}{\tau }_i{\lambda }_{max}\left(Q_i\right)$. Thus, the time derivative of the functional (12) along the trajectories of system (8) is

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}v\left(e_t\right)=& 2e^T\left(t\right)P\dot{e}\left(t\right)+\sum _{i=1}^m{e}^T\left(t\right)Q_{i}e\left(t\right)-\sum _{i=1}^m{e}^T(t-{\tau }_i){\mathsf{e}}^{-2\sigma {\tau }_i}{Q}_{i}e(t-{\tau }_i)\hfill \\ \hfill & -2\sigma \sum _{i=1}^m{\int }_{-{\tau }_i}^0{e}^T(t+\theta ){\mathsf{e}}^{2\sigma \theta }{Q}_{i}e(t+\theta )d\theta .\hfill \end{array}$$

(14)

Using (8), notice that

$$\begin{array}{cc}\hfill 2e^T\left(t\right)P\dot{e}\left(t\right)=& 2e^T\left(t\right)P\left[\sum _{i=0}^m\left[A_i-L_i{C}_i\right]e(t-{\tau }_i)+\eta (x\left(t\right),{\widehat{x}}_t)\right]\hfill \\ \hfill =& {\epsilon }_m\left(t\right)\left[{\left[P\mathrm{\Lambda }-\mathrm{\Gamma }\right]}^T{\rm Y}+{{\rm Y}}^T\left[P\mathrm{\Lambda }-\mathrm{\Gamma }\right]\right]{\epsilon }_m^T\left(t\right)+2e^T\left(t\right)P\eta (x_t,{\widehat{x}}_t),\hfill \end{array}$$

(15)

where ${\epsilon }_m\left(t\right)=[e\left(t\right)e(t-{\tau }_1)\dots e(t-{\tau }_m)]$ and $R_i=P{L}_i$, $i=1,\dots ,m$. Therefore, the observer gain are $L_i=P^{-1}{R}_i$, $i=1,\dots ,m$. While, if $\left|\right|e_t{\left|\right|}_{\tau }<\gamma$ by (9) we get

$$\begin{array}{cc}\hfill 2e^T\left(t\right)P\eta (x_t,{\widehat{x}}_t)& \le 2\left|\right|e\left(t\right)\left|\right|\left|\right|P\left|\right|\left|\right|\eta (x_t,{\widehat{x}}_t)\left|\right|\le 2M\left|\right|e\left(t\right)\left|\right|\left|\right|P\left|\right|\left|\right|e_t{\left|\right|}_{\tau }\hfill \\ \hfill & \le 2\gamma {\alpha }_{3}M\left|\right|e\left(t\right)\left|\right|={\mathcal{O}\left(\right|\left|e\left(t\right)\right||}_{\tau }).\hfill \end{array}$$

Note that, ${\displaystyle \frac{\mathcal{O}\left(\right|\left|e_t\right|{|}_{\tau })}{\left|\right|e_t{\left|\right|}_{\tau }}}$ approaches zero as $\left|\right|e_t{\left|\right|}_{\tau }$ approaches zero. However, using (13) the above equation also satisfies that

$$\begin{array}{cc}\hfill \mathcal{O}\left(\right|\left|e\right(t\left)\right|\left|\right)& =2\gamma {\alpha }_{3}M\left|\right|e\left(t\right)\left|\right|\le {\displaystyle \frac{2\gamma {\alpha }_{3}M}{\sqrt{{\alpha }_1}}}\sqrt{v\left(e_t\right)}.\hfill \end{array}$$

(16)

Thus, substituting (12), (15) and (16) in (14) we have that

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d}{dt}}v\left(e_t\right)+2\sigma v\left(e_t\right)-{\displaystyle \frac{2\gamma {\alpha }_{3}M}{\sqrt{{\alpha }_1}}}\sqrt{v\left(e_t\right)}\le {\epsilon }_m\left(t\right)\left[\mathcal{M}(P,Q_1,\dots ,Q_m,R_0,\dots ,R_m)+2\sigma \mathcal{N}\left(P\right)\right]{\epsilon }_m^T\left(t\right).\hfill \end{array}$$

If (10) holds, then

$${\displaystyle \frac{d}{dt}}v\left(e_t\right)+2\sigma v\left(e_t\right)-{\displaystyle \frac{2\gamma {\alpha }_{3}M}{\sqrt{{\alpha }_1}}}\sqrt{v\left(e_t\right)}<0.$$

Now, if both sides of the above equation are multiplied by ${\mathrm{e}}^{2\sigma \theta }$, then it can be rewritten as

$${\displaystyle \frac{\partial }{\partial \theta }}\left({\mathrm{e}}^{2\sigma \theta }v\left(e_{\theta }\right)\right)\le {\displaystyle \frac{2\gamma {\alpha }_{3}M{\mathrm{e}}^{\sigma \theta }}{\sqrt{{\alpha }_1}}}\sqrt{{\mathrm{e}}^{2\sigma \theta }v\left(e_{\theta }\right)}.$$

Thus, integrating from 0 to t gives

$$\begin{array}{cc}\hfill & {\int }_0^t{\displaystyle \frac{\frac{\partial }{\partial \theta }\left({\mathrm{e}}^{2\sigma \theta }v\left(e_{\theta }\right)\right)}{\sqrt{{\mathrm{e}}^{2\sigma \theta }v\left(e_{\theta }\right)}}}d\theta \le {\int }_0^t{\displaystyle \frac{2\gamma {\alpha }_{3}M{\mathrm{e}}^{\sigma \theta }}{\sqrt{{\alpha }_1}}}d\theta \hfill \\ \hfill & \Rightarrow \sqrt{{\mathrm{e}}^{2\sigma t}v\left(e_t\right)}-\sqrt{v\left(\phi \right)}\le {\displaystyle \frac{2\gamma {\alpha }_{3}M}{2\sigma \sqrt{{\alpha }_1}}}\left({\mathrm{e}}^{\sigma t}-1\right).\hfill \end{array}$$

Therefore

$$\sqrt{v\left(e_t\right)}\le {\displaystyle \frac{\gamma {\alpha }_{3}M}{\sigma \sqrt{{\alpha }_1}}}\left(1-{\mathrm{e}}^{-\sigma t}\right)+{\mathrm{e}}^{-\sigma t}\sqrt{v\left(\phi \right)}.$$

Applying the property (13) to the previous inequality we have that

$$\begin{array}{cc}\hfill & \left|\right|e(t,\phi )\left|\right|\le \sqrt{v\left(e_t\right)}\le {\displaystyle \frac{\gamma {\alpha }_{3}M}{\sigma {\alpha }_1}}\left(1-{\mathrm{e}}^{-\sigma t}\right)+{\mathrm{e}}^{-\sigma t}\sqrt{{\displaystyle \frac{{\alpha }_2}{{\alpha }_1}}}{\left|\right|\phi \left|\right|}_{\tau }\hfill \\ \hfill & \le \left(\sqrt{{\displaystyle \frac{{\alpha }_2}{{\alpha }_1}}}{\left|\right|\phi \left|\right|}_{\tau }-{\displaystyle \frac{\gamma {\alpha }_{3}M}{\sigma {\alpha }_1}}\right){\mathrm{e}}^{-\sigma t}+{\displaystyle \frac{\gamma {\alpha }_{3}M}{\sigma {\alpha }_1}},\forall t>0.\hfill \end{array}$$

Further, $\left|\right|e(t,\phi )\left|\right|$ converges to the limit $\mu ={\displaystyle \frac{\gamma {\alpha }_{3}M}{\sigma {\alpha }_1}}$ as t tends to infinity. Note that

$$\begin{array}{cc}\hfill & \left|\right|e(t,\phi )\left|\right|\le \left(\sqrt{{\displaystyle \frac{{\alpha }_2}{{\alpha }_1}}}{\left|\right|\phi \left|\right|}_{\tau }-\mu \right){\mathrm{e}}^{-\sigma t}+\mu \le \mu \hfill \\ \hfill & ⟺t:=T(\mu ,\phi )\ge {\displaystyle \frac{1}{\sigma }}ln\left(\sqrt{{\displaystyle \frac{{\alpha }_2}{{\alpha }_1}}}{\left|\right|\phi \left|\right|}_{\tau }-\mu \right).\hfill \end{array}$$

□

4. Compartmental Model and Its State Observer

In this section, a SIR epidemiological model with time-delays is proposed for the study of infectious diseases. In addition, because sometimes some population compartments are not available, a state observer is presented for their estimate.

Since March 12, 2020, the World Health Organization (WHO) has declared the infectious disease COVID-19 a pandemic [37]. Because this disease, caused by the SARS-CoV-2 virus, spreads rapidly, it caused severe collapses in various sectors of the population, especially in the hospital, public, and private health sectors. There were effects on the supply of imported inputs as a consequence of the global measures to combat the pandemic [38,39]. This pandemic continues today, in addition to the fact that the SARS-CoV-2 virus is constantly mutating, changing the disease patterns due to new variant strains. The latter present different characteristics, in aspects such as: the prevalence time of the disease in people, the amount of viral load, the effectiveness of the virus to infect cells, the ease of propagation, the efficacy of vaccines, among others, [40]. Therefore, it is important to understand the dynamics of COVID-19 as it evolves.

4.1. SIR Type Compartmental Model with Time-Delays

Model (1) has set the standard for the study of numerous infectious diseases. In order to be used for the COVID-19 pandemic, this section begins with the proposal of a model with time-delays, where it is considered that a part of the recovered population becomes susceptible again after a period of time after their recovery. In addition, three dead-times are incorporated due to the exchange of individuals between the population compartments:

$$\begin{array}{cc}\hfill \dot{S}\left(t\right)=& -\beta S\left(t\right)I\left(t-{\tau }_3\right)+\alpha R(t-{\tau }_2),\hfill \\ \hfill \dot{I}\left(t\right)=& \beta S\left(t\right)I\left(t-{\tau }_3\right)-\gamma I\left(t-{\tau }_1\right),\hfill \\ \hfill \dot{R}\left(t\right)=& \gamma I\left(t-{\tau }_1\right)-\alpha R(t-{\tau }_2),\hfill \end{array}$$

(17)

where $\alpha$ is the rate of loss of immunity and becoming susceptible.

This model is proposed because the current COVID-19 pandemic has shown that individuals who pass from the infected to recovered compartments do not acquire total or permanent immunity, that is, they can be reinfected by another variant of the disease after a time after having recovered. Thus, $\alpha R(t-{\tau }_2)$ represents the individuals that in a dead-time ${\tau }_2$ lose immunity, so they become susceptible again. While $\gamma I\left(t-{\tau }_1\right)$ is the fraction of the infectious that in a dead-time ${\tau }_1$ cease to be infectious and therefore they will recover. Finally, $\beta S\left(t\right)I\left(t-{\tau }_3\right)$ are all those susceptible individuals who became infected and after a dead-time ${\tau }_3$ start to become infectious.

On the other hand, the effective reproduction number ${\mathcal{R}}_e$ can define the threshold behavior of epidemiological models of the form (17): if ${\mathcal{R}}_e<1$ the epidemic event will tend to disappear quickly, while if ${\mathcal{R}}_e>1$ the epidemic event is expected to continue. Here, the higher the value of ${\mathcal{R}}_e$, the more difficult it is to control an epidemic [41]. From model (17), it can be seen that ${\displaystyle \frac{d}{dt}}I\left(t\right)>0$ if and only if $\beta S\left(t\right)I\left(t-{\tau }_3\right)-\gamma I\left(t-{\tau }_1\right)>0$, or equivalently if ${\displaystyle \frac{\beta S\left(t\right)}{\gamma }}>{\displaystyle \frac{I\left(t-{\tau }_1\right)}{I\left(t-{\tau }_3\right)}}$, or ${\displaystyle \frac{\beta S\left(t\right)}{\gamma }}{\displaystyle \frac{I\left(t-{\tau }_3\right)}{I\left(t-{\tau }_1\right)}}>1$. Therefore, an epidemic/pandemic continues if ${\mathcal{R}}_e={\displaystyle \frac{\beta S_0}{\gamma }}{\displaystyle \frac{I\left(t-{\tau }_3\right)}{I\left(t-{\tau }_1\right)}}>1$, where $S_0$ is the initial condition of the susceptible population. Note that ${\displaystyle \frac{\beta S_0}{\gamma }}={\mathcal{R}}_0$ is the basic reproductive number of the classical SIR model (1), which is obtained as follows:

$$\begin{array}{c}\hfill \dot{I}\left(t\right)=(\beta S_0-\gamma )I\left(t\right)\Rightarrow I\left(t\right)=I_0{\mathrm{e}}^{(\beta S_0-\gamma )t}=I_0{\mathrm{e}}^{\gamma \left({\displaystyle \frac{\beta S_0}{\gamma }}-1\right)t}=I_0{\mathrm{e}}^{\gamma \left({\mathcal{R}}_0-1\right)t}.\end{array}$$

Thus, an infectious disease can be classified as an epidemic if ${\mathcal{R}}_e>1$, this will be extinguished when ${\mathcal{R}}_e<1$ and this will remain constant if ${\mathcal{R}}_e=1$, see [42].

4.2. Compartmental Model Transformation

Next, a change of variable is proposed to obtain a new model from (17) with a nominal part and a non-linear part which will serve to a large extent for stability analysis.

Let us consider the following change of variable

$$x_1\left(t\right)=S\left(t\right)-S^{\ast },x_2\left(t\right)=I\left(t\right)-I^{\ast }\mathrm{and}{x}_3\left(t\right)=R\left(t\right)-R^{\ast },$$

(18)

where $S^{\ast }={\displaystyle \frac{\gamma }{\beta }}$, $I^{\ast }=\alpha$, and $R^{\ast }=\gamma$ constitute the equilibrium point of (17). Then, by defining $x^T\left(t\right)=\left[\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)\end{array}\right]$ model (17) can be rewritten in form (6) as follows

$$\begin{array}{c}\hfill \dot{x}\left(t\right)=\sum _{i=0}^3{A}_{i}x(t-{\tau }_i)+g(x_1\left(t\right),x_2(t-{\tau }_3)),\end{array}$$

(19)

where

$$\begin{array}{c}{A}_0=\left[\begin{array}{ccc}-\alpha \beta & 0& 0\\ \alpha \beta & 0& 0\\ 0& 0& 0\end{array}\right],A_1=\left[\begin{array}{ccc}0& 0& 0\\ 0& -\gamma & 0\\ 0& \gamma & 0\end{array}\right],A_2=\left[\begin{array}{ccc}0& 0& \alpha \\ 0& 0& 0\\ 0& 0& -\alpha \end{array}\right],A_3=\left[\begin{array}{ccc}0& -\gamma & 0\\ 0& \gamma & 0\\ 0& 0& 0\end{array}\right],\\ \mathrm{and}g(x_1\left(t\right),x_2(t-{\tau }_3))=\left[\begin{array}{c}-\beta x_1\left(t\right)x_2(t-{\tau }_3)\\ \beta x_1\left(t\right)x_2(t-{\tau }_3)\\ 0\end{array}\right].\end{array}$$

Remark 3. 

Since only the infectious population $I\left(t\right)$ is available as an output from the SIR system, thus $y\left(t\right)=I\left(t\right)-I^{\ast }$. Note also that the entire history of the population is available, then $y\left(t\right)=I(t-{\tau }_i)$. Considering $C_i=\left[010\right]$, $i=0,1,2,3$, the output vector can be expressed as

$$\begin{array}{c}\hfill y\left(t\right)=C_{0}x\left(t\right)+C_{1}x(t-{\tau }_1)+C_{2}x(t-{\tau }_2)+C_{3}x(t-{\tau }_3).\end{array}$$

(20)

Observer Design and Convergence Analysis

Following (7), the state observer for (19)–(20) is of the form

$$\begin{array}{c}\hfill \dot{\widehat{x}}\left(t\right)=\sum _{i=0}^3{A}_i\widehat{x}(t-{\tau }_i)+\sum _{i=0}^3{L}_i\left[C_{i}x(t-{\tau }_i)-C_i\widehat{x}(t-{\tau }_i)\right]+g({\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3)),\end{array}$$

(21)

where ${\widehat{x}}^T\left(t\right)=\left[{\widehat{x}}_1\left(t\right){\widehat{x}}_2\left(t\right){\widehat{x}}_3\left(t\right)\right]$, $L_i^T=\left[L_{i_1}{L}_{i_2}{L}_{i_3}\right]$ and $C_i=\left[C_{i_1}{C}_{i_2}{C}_{i_3}\right]$. Thus, by (8) the dynamics of the observation error of (19) and (21) is

$$\begin{array}{c}\hfill \dot{e}\left(t\right)=\sum _{i=0}^3\left[A_i-L_i{C}_i\right]e(t-{\tau }_i)+\eta (x_1\left(t\right),x_2(t-{\tau }_3),{\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3)),\end{array}$$

(22)

with $\eta (x_1\left(t\right),x_2(t-{\tau }_3),{\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3))=g(x_1\left(t\right),x_2(t-{\tau }_3))-g({\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3))$. Now, by Proposition 1, the observation error (22) satisfies

$$\begin{array}{cc}\hfill \dot{v}\left(e_t\right)+2\sigma & v\left(e_t\right)\le {\epsilon }_3\left(t\right)\left[\mathcal{M}(P,Q_1,Q_2,Q_3,R_0,\dots ,R_3)+2\sigma \mathcal{N}\left(P\right)\right]{\epsilon }_3^T\left(t\right)+\mathcal{O}\left(\right||e_t\left|\right|),\hfill \end{array}$$

(23)

where

$$\begin{array}{c}\hfill \left[\mathcal{M}(P,Q_1,Q_2,Q_3,R_0,\dots ,R_3)+2\sigma \mathcal{N}\left(P\right)\right]=\left[\begin{array}{cccc}{\mathrm{\Theta }}_{1,1}+2\sigma P& (\ast )& (\ast )& (\ast )\\ {\mathrm{\Theta }}_{2,1}& {\mathrm{\Theta }}_{2,2}& (\ast )& (\ast )\\ {\mathrm{\Theta }}_{3,1}& 0& {\mathrm{\Theta }}_{3,3}& (\ast )\\ {\mathrm{\Theta }}_{4,1}& 0& 0& {\mathrm{\Theta }}_{4,4}\end{array}\right].\end{array}$$

(24)

Here, ${\mathrm{\Theta }}_{1,1}=P{A}_0+A_0^{T}P-R_0{C}_0-C_0^T{R}_0^T+Q_1+Q_2+Q_3$, ${\mathrm{\Theta }}_{2,1}=A_1^{T}P-C_1^T{R}_1^T$, ${\mathrm{\Theta }}_{3,1}=A_2^{T}P-C_2^T{R}_2^T$, ${\mathrm{\Theta }}_{4,1}=A_3^{T}P-C_3^T{R}_3^T$, ${\mathrm{\Theta }}_{2,2}=-{\mathsf{e}}^{-2\sigma {\tau }_1}{Q}_1$, ${\mathrm{\Theta }}_{3,3}=-{\mathsf{e}}^{-2\sigma {\tau }_2}{Q}_2$, ${\mathrm{\Theta }}_{4,4}=-{\mathsf{e}}^{-2\sigma {\tau }_3}{Q}_3$ and $(\ast )$ denote the transpose of the symmetric elements. On the other hand, note that

$$\begin{array}{cc}\hfill \left|\right|\eta (x_1\left(t\right),x_2(t-{\tau }_3),{\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3))\left|\right|& =\left|\right|g(x_1\left(t\right),x_2(t-{\tau }_3))-g({\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3))\left|\right|\hfill \\ \hfill & =4{\beta }^2|x_1\left(t\right)e_2(t-{\tau }_3)+{\widehat{x}}_2(t-{\tau }_3)e_1\left(t\right)|,\hfill \end{array}$$

where $e_1\left(t\right)=x_1\left(t\right)-{\widehat{x}}_1\left(t\right)$ and $e_2(t-{\tau }_3)=x_2(t-{\tau }_3)-{\widehat{x}}_2(t-{\tau }_3)$. Now assuming that the populations are bounded, we observe that

$$\begin{array}{cc}\hfill \left|\right|\eta (x_1& \left(t\right),x_2(t-{\tau }_3),{\widehat{x}}_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3)\left)\right||\le 4{\beta }^{2}m|e_2(t-{\tau }_3)+e_1\left(t\right)|\le 4{\beta }^{2}m\left|\right|e_t\left|\right|=M\left|\right|e_t\left|\right|,\hfill \end{array}$$

with $m=max\{x_1\left(t\right),{\widehat{x}}_2(t-{\tau }_3)\}$, and $M=4{\beta }^{2}m$. As mentioned above, the total population $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$ is considered constant and this is normalized, thus $|x_i\left(t\right)|\le 1=m$, $i=1,2,3$, and $\left|\right|e\left(t\right)\left|\right|<\gamma =1$, therefore by (16)

$$\begin{array}{c}\hfill \mathcal{O}\left(\right|\left|e\left(t\right)\right|\left|\right)={\displaystyle \frac{2\gamma {\alpha }_{3}M}{\sqrt{{\alpha }_1}}}\sqrt{v\left(e_t\right)}={\displaystyle \frac{8{\alpha }_3{\beta }^2}{\sqrt{{\alpha }_1}}}\sqrt{v\left(e_t\right)}.\end{array}$$

(25)

Proposition 2. 

If there exist $\sigma \in {\mathbb{R}}^{+}$, vectors $R_0$, $R_i\in {\mathbb{R}}^{3\times 1}$ and positive definite matrices $P,Q_i>0$ for $i=1,2,3$; such that (24) is a negative definite matrix, then the observation error (22) is exponentially convergent to a ball $B_{\mu }$ with radius $\mu >0$.

Proof. 

The result follows from the Proposition 1. □

Now, we express the equivalent observer in original coordinates, to this end from (21) we have that

$$\begin{array}{cc}\hfill {\dot{\widehat{x}}}_1\left(t\right)=& -\alpha \beta {\widehat{x}}_1\left(t\right)+\alpha {\widehat{x}}_3(t-{\tau }_2)-\gamma {\widehat{x}}_2(t-{\tau }_3)-\beta {\widehat{x}}_1\left(t\right){\widehat{x}}_2(t-{\tau }_3)\hfill \\ \hfill & +\sum _{i=0}^3\sum _{j=1}^3{L}_{i_1}{C}_{i_j}\left[x_j(t-{\tau }_i)-{\widehat{x}}_j(t-{\tau }_i)\right],\hfill \\ \hfill {\dot{\widehat{x}}}_2\left(t\right)=& \alpha \beta {\widehat{x}}_1\left(t\right)-\gamma {\widehat{x}}_2(t-{\tau }_1)+\gamma {\widehat{x}}_2(t-{\tau }_3)+\beta {\widehat{x}}_1\left(t\right){\widehat{x}}_2(t-{\tau }_3)\hfill \\ \hfill & +\sum _{i=0}^3\sum _{j=1}^3{L}_{i_2}{C}_{i_j}\left[x_j(t-{\tau }_i)-{\widehat{x}}_j(t-{\tau }_i)\right],\hfill \\ \hfill {\dot{\widehat{x}}}_3\left(t\right)=& \gamma {\widehat{x}}_2(t-{\tau }_1)-\alpha {\widehat{x}}_3(t-{\tau }_2)+\sum _{i=0}^3\sum _{j=1}^3{L}_{i_3}{C}_{i_j}\left[x_j(t-{\tau }_i)-{\widehat{x}}_j(t-{\tau }_i)\right],\hfill \end{array}$$

and following (18) the observer of the model SIR given in (17) takes the form

$$\begin{array}{cc}\hfill \dot{\widehat{S}}\left(t\right)=& \sum _{i=0}^3{L}_{i_1}[C_{i_1}\left[S(t-{\tau }_i)-\widehat{S}(t-{\tau }_i)\right]+C_{i_2}\left[I(t-{\tau }_i)-\widehat{I}(t-{\tau }_i)\right]\hfill \\ \hfill & +C_{i_3}\left[R(t-{\tau }_i)-\widehat{R}(t-{\tau }_i)\right]]-\beta \widehat{S}\left(t\right)\widehat{I}(t-{\tau }_3)+\alpha \widehat{R}(t-{\tau }_2),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \dot{\widehat{I}}\left(t\right)=& \sum _{i=0}^3{L}_{i_2}[C_{i_1}\left[S(t-{\tau }_i)-\widehat{S}(t-{\tau }_i)\right]+C_{i_2}\left[I(t-{\tau }_i)-\widehat{I}(t-{\tau }_i)\right]\hfill \\ \hfill & +C_{i_3}\left[R(t-{\tau }_i)-\widehat{R}(t-{\tau }_i)\right]]+\beta \widehat{S}\left(t\right)\widehat{I}(t-{\tau }_3)-\gamma \widehat{I}(t-{\tau }_1),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \dot{\widehat{R}}\left(t\right)=& \sum _{i=0}^3{L}_{i_3}[C_{i_1}\left[S(t-{\tau }_i)-\widehat{S}(t-{\tau }_i)\right]+C_{i_2}\left[I(t-{\tau }_i)-\widehat{I}(t-{\tau }_i)\right]\hfill \\ \hfill & +C_{i_3}\left[R(t-{\tau }_i)-\widehat{R}(t-{\tau }_i)\right]]+\gamma \widehat{I}\left(t-{\tau }_1\right)-\alpha \widehat{R}(t-{\tau }_2).\hfill \end{array}$$

(28)

5. Implementation of Results to the Population of Mexico

In this section, the SIR type compartmental model with time-delays proposed in (17) and the Luenberger-type state observer given in (26)–(28) are used to estimate the compartmental populations of Susceptible $S\left(t\right)$ and Recovered$R\left(t\right)$ not reported by the WHO. Here, the observer only considers data of cases confirmed of Infectious$I\left(t\right)$ in Mexico, reported by [43]. Matlab^® software is used for these simulations.

5.1. Implementation Description

To estimate $S\left(t\right)$ and $R\left(t\right)$ the following assumptions are made:

A1. 

Compartmental populations are considered closed and normalized, that is, $S\left(t\right)+R\left(t\right)+I\left(t\right)=1$.

A2. 

Only the population of infectious $I\left(t\right)$ is available. The real data is denoted by $I_r\left(t\right)$.

By A2, the state observer (26)–(28) can be rewritten as:

$$\begin{array}{cc}\hfill \dot{\widehat{S}}\left(t\right)=& -\beta \widehat{S}\left(t\right)\widehat{I}(t-{\tau }_3)+\alpha \widehat{R}(t-{\tau }_2)+\sum _{i=0}^3{L}_{i_1}{C}_{i_2}\left[I(t-{\tau }_i)-\widehat{I}(t-{\tau }_i)\right],\hfill \\ \hfill \dot{\widehat{I}}\left(t\right)=& \beta \widehat{S}\left(t\right)\widehat{I}(t-{\tau }_3)-\gamma \widehat{I}(t-{\tau }_1)+\sum _{i=0}^3{L}_{i_2}{C}_{i_2}\left[I(t-{\tau }_i)-\widehat{I}(t-{\tau }_i)\right],\hfill \\ \hfill \dot{\widehat{R}}\left(t\right)=& \gamma \widehat{I}\left(t-{\tau }_1\right)-\alpha \widehat{R}(t-{\tau }_2)+\sum _{i=0}^3{L}_{i_3}{C}_{i_2}\left[I(t-{\tau }_i)-\widehat{I}(t-{\tau }_i)\right].\hfill \end{array}$$

(29)

By A1, the third equation of (29) can be replaced by $R\left(t\right)$=1-$\widehat{S}\left(t\right)$-$\widehat{I}\left(t\right)$. Next, the model (17) and the observer (29) are implemented in two processes: tuning of model parameters and of estimation of compartmental populations, which are described below.

5.1.1. First Process: Tuning of Model Parameters

Assume A1. For given initial conditions ${\alpha }_0,{\beta }_0,{\gamma }_0,\in {\mathbb{R}}^{+}$, and $S_0\left(\theta \right),I_0\left(\theta \right),R_0\left(\theta \right)\in \mathfrak{C}$, the epidemiological model (17) is used to tune the parameters $\alpha ,\beta ,\gamma ,{\tau }_1,{\tau }_2,{\tau }_3\in {\mathbb{R}}^{+}$ that best fit the curve of $I_r\left(t\right)$ reported by the WHO in a 55-day time window which is indeed the real permissible range for each parameter. In this case and without loss of generality, the 55-day time window was selected according to its variability/fruition in the curve/data and the greatest impact it had on the population of Mexico. This is only an illustrative example and can be changed according to the user’s preferences. This process is shown Figure 1 and detailed below:

• This process begins by proposing initial conditions (ic) for the populations $S_0\left(\theta \right)$, $R_0\left(\theta \right)$, $I_0\left(\theta \right)\in \mathfrak{C}$, with $I_0\left(\theta \right)$ close to $I_r\left(\theta \right)$. Also, knowing a permissible range of variation of each parameter, an initial tuning of the model parameters ${\alpha }_0,{\beta }_0,{\gamma }_0,{\tau }_{1_o},{\tau }_{2_o},{\tau }_{3_o}\in {\mathbb{R}}^{+}$ is proposed.
• The response of the model (17) is obtained by using the initial conditions and the model parameters.
• Given the real data $I_r\left(t\right)$, the Mean Squared Error (MSE) is calculated between $I\left(t\right)$ and $I_r\left(t\right)$, i.e., MSE$\left(e_2\right)=\mathrm{MSE}(I\left(t\right)-I_r\left(t\right))$.
• Given a $\delta >0$, if MSE$\left(e_2\right)<\delta$, then the parameters $\alpha ={\alpha }_k,\beta ={\beta }_k,\gamma ={\gamma }_k,{\tau }_1={\tau }_{1_k},{\tau }_2={\tau }_{2_k},{\tau }_3={\tau }_{3_k}$, for some $k=0,1,\dots$, are used for the second process and the first process ends; otherwise, new $kth$ model parameters are proposed within the allowed ranges and return to step (2).

5.1.2. Second Process: Estimation of Compartmental Populations

Next, the parameters $\alpha ,\beta ,\gamma ,{\tau }_1,{\tau }_2,{\tau }_3\in {\mathbb{R}}^{+}$ obtained in the first process are used on the observer given in (29) to estimate the compartmental populations $S\left(t\right)$ and $R\left(t\right)$ not reported by the WHO. Here, it is assumed that only $I\left(t\right)$ is available, that is, $\widehat{y}\left(t\right)={\sum }_{i=0}^3{C}_{i_2}I(t-{\tau }_i)$. Also the same real data $I_r\left(t\right)$ from the first process is utilized. In Figure 2 this second process is schematized.

5.2. Application of Results via Simulation

In this subsection, what was previously explained on the first and second processes is applied, for which we use the solver Sedumi [44] in Matlab^® to determine feasibility of the LMI (24). By this and Proposition 2, it is possible to guarantee an exponential decay in the observation error $e\left(t\right)$, in addition to obtaining the observer gains $L_i$, $i=0,\dots ,3$. Finally, to corroborate the efficiency of the procedure and the obtained results, a simulation is performed using Simulink in Matlab^®.

5.2.1. LMI Feasibility

In the first process, we consider: the rate of immunity loss $\alpha =1$, the infection rate $\beta =0.0109$, the recovery rate $\gamma =0.085$, the dead-time by recovery ${\tau }_1=3$, the dead-time by immunity loss ${\tau }_2=5$ and the dead-time by incubation ${\tau }_3=6$.

Let $C_i=\left[\begin{array}{cc}0& 1\end{array}\right]$, $i=0,1,2,3$. Then, the LMI conditions (24), with the appropriate dimensions, is feasible for $\sigma =0.003$,

$$\begin{array}{cc}\hfill P=& \left[\begin{array}{cc}57.7477& -20.8676\\ -20.8676& 36.0528\end{array}\right],Q_1=\left[\begin{array}{cc}0.1890& -0.0028\\ -0.0028& 0.2022\end{array}\right],\hfill \\ \hfill Q_2=& \left[\begin{array}{cc}0.1206& -0.0028\\ -0.0028& 0.1170\end{array}\right]\mathrm{and}{Q}_3=\left[\begin{array}{cc}0.1207& -0.0029\\ -0.0029& 0.1173\end{array}\right],\hfill \end{array}$$

as well as

$$R_0=\left[\begin{array}{c}5.8150\\ 34.8704\end{array}\right],R_1=\left[\begin{array}{c}1.7490\\ -3.2965\end{array}\right],R_2=\left[\begin{array}{c}-0.0198\\ -0.2073\end{array}\right]\mathrm{and}{R}_3=\left[\begin{array}{c}-6.7004\\ 4.6352\end{array}\right].$$

Thus, by Proposition 2 the variables of the observation error $e_2\left(t\right)=I_r\left(t\right)-\widehat{I}\left(t\right)$ and $e_3\left(t\right)=R\left(t\right)-\widehat{R}\left(t\right)$ satisfy the exponential decay given in (11), with $\mu =0.4771$ and $T=770.1474$. Even more, the observer gains obtained are

$$L_0=\left[\begin{array}{c}0.5693\\ 1.2967\end{array}\right],L_1=\left[\begin{array}{c}-0.0035\\ -0.0935\end{array}\right],L_2=\left[\begin{array}{c}-0.0031\\ -0.0075\end{array}\right]\mathrm{and}{L}_3=\left[\begin{array}{c}-0.0880\\ 0.0776\end{array}\right].$$

(30)

5.2.2. Observer Simulation

For simulation purposes, consider that $N=350000$ is the total population, of which $S_0=1-I_0-R_0$, where $I_0=69025/N$ and $R_0=0$ are the normalized initial conditions of the model, corresponding to the Susceptible, Infectious and Recovered population, respectively. The initial conditions of the observer are ${\widehat{S}}_0=1-{\widehat{I}}_0-{\widehat{R}}_0$, ${\widehat{I}}_0=0.225$ and ${\widehat{R}}_0=0.05$. Besides the model parameters and the observer gains are given in (30). Next, simulations carried out in Simulink in Matlab^® are presented.

The simulation between the compartmental model (17) and the state observer (29) using a solution step size of $t=0.01$Day is shown in Figure 3. Here, given $\delta =0.5\times {10}^{-2}$, we have that the MSE$\left(e_1\left(t\right)\right)=$MSE$(S\left(t\right)-\widehat{S}\left(t\right))=9.524\times {10}^{-2}$, MSE$\left(e_2\left(t\right)\right)=$MSE$(I\left(t\right)-\widehat{I}\left(t\right))=0.443\times {10}^{-2}$ and MSE$\left(e_3\left(t\right)\right)=$MSE$(R\left(t\right)-\widehat{R}\left(t\right))=9.601\times {10}^{-2}$ in the observed dynamics of Susceptible, Infectious and Recovered, respectively.

Now, the observer’s behavior is analyzed using real data of the Infectious population $I_r\left(t\right)$ recorded by the WHO of the 13 January to March 8 of 2022 for the COVID-19 pandemic. This period was a time window with the greatest impact on the population of Mexico.

Estimates $\widehat{S}\left(t\right)$ and $\widehat{R}\left(t\right)$ of the compartmental populations of Susceptible and Recovered, respectively, not reported by the WHO are presented in Figure 4. Here, the solution step size used is $t=1$Day and we obtain a MSE$\left(e_2\left(t\right)\right)=$MSE$(I_r\left(t\right)-\widehat{I}\left(t\right))=1.0203\times {10}^{-2}$ in the dynamic of the Infectious error.

Remark 4. 

Note that the nonlinearity $\eta (x_t,{\widehat{x}}_t)$ of the observation error (22) implies that an asymptotic convergence of $e\left(t\right)$ to zero cannot be guaranteed, but a convergence to a neighborhood $\mu >0$ at time T can be ensured. Unfortunately this neighborhood $\mu \in [0,1]$ usually is very conservative and for normalized data this does not seem to make sense. However, although $\mu =0.4771$ and $T=770.1474$ days it seems that these do not affect convergence of the observation error $e\left(t\right)$. Since in Figure 4 it can be seen that the variables of $\widehat{I}\left(t\right)$ converge to the data $I_r\left(t\right)$, with a MSE$\left(e_2\right)=1.0203\times {10}^{-2}$.

6. Conclusions

In this paper, an epidemiological model of the SIR type and a Luenberger-type observer to estimate compartmental populations not reported by the WHO are presented. The model proposed is a nonlinear system with three time-delays, which considers the dead-times it takes for an individual to move from one compartmental population to another and that recovered individuals may again be susceptible to the disease. This model is used to determine parameters such as the $\beta$ infection rate, the $\gamma$ recovery rate and the $\alpha$ rate of immunity loss; and for the design of a Luenberger-type state observer with time-delays. To guarantee the estimation of the observer, a convergence analysis is presented using Lyapunov-Krasovskii functionals to ensure that the observation error decreases with an exponential decay $\sigma >0$. This observer is employed to estimate the compartmental populations of Susceptible$S\left(t\right)$ and Recovered$R\left(t\right)$ from only the availability of the real data $I_r\left(t\right)$, Infectious compartmental population reported by WHO in a 55-day time window with the highest impact on the population of Mexico. The results shown suggest that the proposed mathematical model and observer are a suggestive alternative to use to estimate missing compartmental populations.