1. INTRODUCTION

Governments around the world are using the reproduction numbers as criteria when deciding public health policies during COVID-19 pandemic. In principal, there are two types of reproduction numbers: basic reproduction number and effective reproduction number. The basic reproduction number (denoted by R0) shows the average expected number of cases generated by one case in a population where all individuals are susceptible. The effective reproduction number (denoted by Rt) shows the average expected number of cases generated in the current state of a population. In practice, R⁰ is used to determine how many population that needs to be immune to reach herd immunity. In this case, the herd-immunity threshold is 1 − 1 R⁰ . On the other hand, R^t is used to monitor transmissions of the disease in a population during the outbreak. Hence, this number is usually used as one of the main criteria to evaluate the public health policies. Theoretically, R^t needs to be below 1 to ensure the transmission is under control.

Many researchers have provided calculations to estimate R^t using different approaches, e.g., Bayesian estimation [1], serial interval [2], Extended Kalman filter [3], and parameter fitting [4]. If the estimation parameters are tuned properly, all of these approaches will provide similar pattern with small variation. Once the estimated R^t is obtained, we can create short-term forecasts to determine different reopening scenarios [5], [6].

The aim of this paper is to provide a novel approach to estimate the effective reproduction number R^t of infectious diseases. To this end, we design a Nonlinear Observer (NLO) from a discrete-time augmented Susceptible-Infectious-Removed (SIR) model. The method is efficient in the sense that the epidemiological data is injected directly into the model once we found a constant observer gain. The observer gain is obtained by solving a Linear Matrix Inequality (LMI). While the majority of estimation methods are based on stochastic process, this new approach is deterministic. The Confidence Interval (CI) provided in our estimation is inherited from the uncertainty in the infectious time and not from the method itself.

The paper is organized as follow. In Section II, we derive a discrete-time augmented SIR model. In this section, we consider the time-varying R^t as an augmented state. Furthermore, we assume its value is a piece-wise constant function. In Section III, we derive a sufficient condition for the observer gain in terms of LMI. Simulation results using epidemiological data from Jakarta is presented in Section IV. Extension to other compartmental model can be found in Section V. Finally, conclusions are given in Section VI.

40 Hasan, A.

2. MATHEMATICAL MODEL

We use a simple SIR model in this paper for several reasons. First, it can be used to described transmissions of many infectious diseases. Furthermore and most importantly, the actual data of the three compartments are available in most outbreak events. The removed compartment consists of individuals who are either recovered or deceased. Assuming constancy of population N, the SIR model can be written as follow [7]:

\[\dot{S}(t) = -\frac{\beta I(t)S(t)}{N},\tag{1}\]

\[\dot{I}(t) = \frac{\beta I(t)S(t)}{N} - \gamma I(t), \tag{2}\]

\[\dot{R}(t) = \gamma I(t), \tag{3}\] where S denotes the number of susceptible individuals, I denotes the number of infected individuals or active cases, and R denotes the number removed individuals. The model has two parameters: the transmission rate \(\beta\) and the removal rate \(\gamma\). By definition, \(\beta\) is the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious individual. Thus, in principal \(\beta\) is time-varying due to interventions. For this reason, in the remaining of this paper we consider \(\beta\) as a parameter that depends on time t and is unknown. On the other hand, the removal rate \(\gamma\) is an inverse of the average infectious time, i.e., \(\gamma = \frac{1}{T_i}\). The infectious time can be obtained through medical data and is usually known together with its Confidence Interval (CI). Taking into account reduction in the number of susceptible individuals, the effective reproduction number can be estimated as follow [8]:

\[\mathcal{R}_t(t) = \frac{S(t)}{N} \left( \frac{\beta(t)}{\gamma} \right) \approx \frac{\beta(t)}{\gamma}.\] (4)

Discretizing (1)-(3) using the Euler discretization method, substituting \(\beta(t) = \gamma \mathcal{R}_t(t)\) into the model, and augmenting \(\mathcal{R}_t(t)\) as a new state variable, we obtain the following discrete-time augmented SIR model:

\[S(k+1) = S(k) - \frac{\gamma \Delta t}{N} \mathcal{R}_t(k) I(k) S(k), \tag{5}\]

\[I(k+1) = (1 - \gamma \Delta t)I(k) + \frac{\gamma \Delta t}{N} \mathcal{R}_t(k)I(k)S(k), \tag{6}\]

\[R(k+1) = R(k) + \gamma \Delta t I(k), \tag{7}\]

\[\mathcal{R}_t(k+1) = \mathcal{R}_t(k). \tag{8}\]

Remark that, in the last equation we assume \(\mathcal{R}_t\) as a piece-wise constant function with jumps every time new data come in. To simplify the model, let us define:

\[\boldsymbol{x}(k) = (S(k) \quad I(k) \quad R(k) \quad \mathcal{R}_t(k))^{\mathsf{T}}. \tag{9}\]

The discrete-time augmented SIR model (5)-(8), can be written as the following nonlinear state-space representation:

\[x(k+1) = Ax(k) + f(x(k)), \tag{10}\]

where

\[\mathbf{A} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 - \gamma \Delta t & 0 & 0 \\ 0 & \gamma \Delta t & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \text{ and } \mathbf{f}(\mathbf{x}(k)) = \begin{pmatrix} -\frac{\gamma \Delta t}{N} \mathcal{R}_t(k) I(k) S(k) \\ \frac{\gamma \Delta t}{N} \mathcal{R}_t(k) I(k) S(k) \\ 0 \\ 0 \end{pmatrix}. \tag{11}\]

Since data for S(k), I(k), and R(k) are available, the measurement vector \(y(k) \in \mathbb{R}^3\) is given by:

\[y(k) = Cx(k), \tag{12}\] where

\[C = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}. \tag{13}\]

3. METHOD

We design the NLO as follow:

\[\bar{\boldsymbol{x}}(k+1) = \boldsymbol{A}\bar{\boldsymbol{x}}(k) + \boldsymbol{f}(\bar{\boldsymbol{x}}(k)) + \boldsymbol{K}(\boldsymbol{y}(k) - \bar{\boldsymbol{y}}(k)), \qquad (14)\]

where \(\bar{\boldsymbol{x}}(k) \in \mathbb{R}^4\) is the estimated state from the NLO and \(\boldsymbol{K} \in \mathbb{R}^{4 \times 3}\) is the observer gain to be determined later. Let \(\tilde{\boldsymbol{x}}(k) = \boldsymbol{x}(k) - \bar{\boldsymbol{x}}(k)\), \(\tilde{\boldsymbol{y}}(k) = \boldsymbol{y}(k) - \bar{\boldsymbol{y}}(k)\), and \(\Delta \boldsymbol{f}(k) = \boldsymbol{f}(\boldsymbol{x}(k)) - \boldsymbol{f}(\bar{\boldsymbol{x}}(k))\). Subtracting (10) with (14), we obtain:

\[\tilde{\boldsymbol{x}}(k+1) = (\boldsymbol{A} - \boldsymbol{K}\boldsymbol{C})\,\tilde{\boldsymbol{x}}(k) + \Delta\boldsymbol{f}(k). \tag{15}\]

The problem is to find K such that the error \(\tilde{x}(k)\) asymptotically converges toward zero, which guarantees the estimated states converge to the actual states. Designing an NLO for a nonlinear system is not trivial and sometime impossible. Researchers usually simplify the problem using assumptions, for example by assuming the non-linearity is locally Lipschitz and bounded. These assumptions are common when designing NLO for nonlinear systems, e.g., see [9], [10].

Assumption 1. The nonlinear function f is a one-sided Lipschitz, i.e., it satisfies

\[\epsilon_1 \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix}^{\mathsf{T}} \begin{pmatrix} -\boldsymbol{I}_4 & -\frac{1}{2}\boldsymbol{I}_4 \\ -\frac{1}{2}\boldsymbol{I}_4 & \boldsymbol{0}_n \end{pmatrix} \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix} \ge 0, \tag{16}\]

for \(\epsilon_1 > 0\), where \(I_4\) denotes the identity matrix of size 4

Assumption 2. The nonlinear function f satisfies the quadratic inner-boundedness condition, i.e.,

\[\epsilon_{2} \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix}^{\mathsf{T}} \begin{pmatrix} \boldsymbol{I}_{4} & \frac{\nu}{2} \boldsymbol{I}_{4} \\ \frac{\nu}{2} \boldsymbol{I}_{4} & -\boldsymbol{I}_{4} \end{pmatrix} \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix} \ge 0. \tag{17}\]

for \(\epsilon_2 > 0\) and \(\nu \in \mathbb{R}\).

Our augmented model satisfies both assumptions. From here, we can derive a sufficient condition for the observer gain K, as stated in Theorem III.2.

Lemma 3.1. The discrete-time augmented SIR model (5)-(8) is locally Lipschitz and bounded.

Proof: Since \(\mathcal{R}_t(k+1) = \mathcal{R}_t(k)\), the Jacobian of (5)-(6) is linear. Thus, (5)-(6) are locally Lipschitz. Furthermore from (5), for S(0), I(0) > 0 the solutions S(k) and I(k) will remain bounded. See detailed explanation in [11].

Theorem 3.2. For \(\epsilon_1, \epsilon_2 > 0\) and \(\nu \in \mathbb{R}\), the error dynamics (15) is asymptotically stable if there exist matrices \(G = G^{\mathsf{T}} > 0\) and \(R \in \mathbb{R}^{3 \times 4}\) such that the following LMI holds

\[\begin{pmatrix} -G + (\epsilon_2 - \epsilon_1)I_4 & A^{\mathsf{T}}G - C^{\mathsf{T}}R + \frac{\nu\epsilon_2 - \epsilon_1}{2}I_4 & A^{\mathsf{T}}G - C^{\mathsf{T}}R \\ GA - R^{\mathsf{T}}C + \frac{\nu\epsilon_2 - \epsilon_1}{2}I_4 & G - \epsilon_2I_4 & 0 \\ GA - R^{\mathsf{T}}C & 0 & -G \end{pmatrix} < 0.\] (18)

Furthermore, the observer gain is given by:

\[K = G^{-1}R^{\mathsf{T}}. (19)\]

Proof: Let us define a Lyapunov function:

\[V(k) = \tilde{\boldsymbol{x}}(k)^{\mathsf{T}} \boldsymbol{G} \tilde{\boldsymbol{x}}(k). \tag{20}\]

Thus, if \(\Delta V(k+1) = V(k+1) - V(k)\), then we have:

\[\Delta V(k+1) = \tilde{\boldsymbol{x}}(k)^{\mathsf{T}} \left( (A - \boldsymbol{K}\boldsymbol{C})^{\mathsf{T}} \boldsymbol{G} (A - \boldsymbol{K}\boldsymbol{C}) - \boldsymbol{G} \right) \tilde{\boldsymbol{x}}(k) + 2\tilde{\boldsymbol{x}}(k)^{\mathsf{T}} \left( A - \boldsymbol{K}\boldsymbol{C} \right)^{\mathsf{T}} \boldsymbol{G} \Delta \boldsymbol{f}(k) + \Delta \boldsymbol{f}(k)^{\mathsf{T}} \boldsymbol{G} \Delta \boldsymbol{f}(k).\] (21)

Expressing the right hand side of (21) as a matrix multiplication, we have:

\[\Delta V(k+1) = \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix}^{\mathsf{T}} \begin{pmatrix} (A - \boldsymbol{K}\boldsymbol{C})^{\mathsf{T}} \boldsymbol{G} (A - \boldsymbol{K}\boldsymbol{C}) - \boldsymbol{G} & (A - \boldsymbol{K}\boldsymbol{C})^{\mathsf{T}} \boldsymbol{G} \\ \boldsymbol{G} (A - \boldsymbol{K}\boldsymbol{C}) & \boldsymbol{G} \end{pmatrix} \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix}. \quad (22)\]

42 Hasan, A.

Adding the left hand side of (16) and (17) into (22), we have:

\[\Delta V(k+1) \leq \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix}^{\mathsf{T}} \boldsymbol{\Phi} \begin{pmatrix} \tilde{\boldsymbol{x}}(k) \\ \Delta \boldsymbol{f}(k) \end{pmatrix}, \tag{23}\] where

\[\mathbf{\Phi} = \begin{pmatrix} (A - \mathbf{K}\mathbf{C})^{\mathsf{T}} \mathbf{G} (A - \mathbf{K}\mathbf{C}) - \mathbf{G} + (\epsilon_2 - \epsilon_1) \mathbf{I}_4 & (A - \mathbf{K}\mathbf{C})^{\mathsf{T}} \mathbf{G} + \frac{\nu \epsilon_2 - \epsilon_1}{2} \mathbf{I}_4 \\ \mathbf{G} (A - \mathbf{K}\mathbf{C}) + \frac{\nu \epsilon_2 - \epsilon_1}{2} \mathbf{I}_4 & \mathbf{G} - \epsilon_2 \mathbf{I}_4 \end{pmatrix}.\](24)

Substituting (19) and applying Schur complement to (24), then \(\Phi < 0\) is equivalent to (18). This completes the proof.

Remark 1. Even though the observer guarantees the stability of the estimate, finding matrix R can be difficult. Unfortunately, there is not systematic methods to obtain it. As for matrix G, it can be chosen as long as the matrix is symmetric and positive definite.

4. ESTIMATION OF \(\mathcal{R}_t\) IN JAKARTA

Daily epidemiological data of COVID-19, such as the number of active case and the number of removed case between April 2020 until August 2020, are used in our estimation. Code and data are available in: https://github.com/agusisma/COVIDNLO. In this simulation, we use the following parameters: \(\epsilon_1=10\), \(\epsilon_2=1\), and \(\nu=9\). Furthermore, the symmetric matrix \({\bf G}\) is chosen as \({\bf G}=0.1{\bf I}_4\), while the matrix \({\bf R}\) is chosen as:

\[\mathbf{R} = \begin{pmatrix} 0.0363 & -0.0032 & 0 & -0.0013 \\ -0.0320 & 0.0745 & 0 & 0.0057 \\ 0 & 0 & 0.0916 & 0 \end{pmatrix}. \tag{25}\]

We assume the average infectious time \(\frac{1}{\gamma} = 12\) days with standard deviation of 3 days. Using these parameters, the LMI (18) is negative definite with the largest and smallest Eigenvalues are -0.0989 and -9.1294, respectively. The observer gain is then given by:

\[\boldsymbol{K} = \begin{pmatrix} 0.3630 & -0.3200 & 0\\ -0.0320 & 0.7453 & 0\\ 0 & 0 & 0.9160\\ -0.0134 & 0.0573 & 0 \end{pmatrix}. \tag{26}\]

Figure 1: Real-time data fitting from the NLO for active and removed case with their estimation errors.

Figure 1 shows estimation results from the NLO for the daily number of active and removed case. It can be observed from the estimation errors that the NLO estimates these numbers reasonably accurate. The NLO is compared with the EKF method presented in [3] and the results can be seen from Figure 2. It can be observed that the estimation results are virtually almost identical. We should note, however, that these results are obtained after a lengthy process of trial and error when determining the matrix R. Unfortunately, there is no method that can be used to determine matrix R systematically. Having said that, the main advantage of using NLO is its stability and efficiency compare to EKF, since the NLO does not require calculation of inverse matrices.

Figure 2: Comparison of the effective reproduction numbers \(\mathcal{R}_t\) from NLO and EKF. The Confidence Interval (CI), indicated by the band, is a result from uncertainty in the infectious time \(\frac{1}{2}\).

5. ESTIMATION OF \(\mathcal{R}_t\) FROM OTHER COMPARTMENTAL MODEL

A natural question is then what happen to the estimation when we have unmeasured compartment? For example, let us consider the following Susceptible-Exposed-Infectious-Removed (SEIR) model:

\[\dot{S}(t) = -\frac{\beta I(t)S(t)}{N}, \qquad (27)\] \[\dot{E}(t) = \frac{\beta I(t)S(t)}{N} - \mu E(t), \qquad (28)\]

\[\dot{E}(t) = \frac{\beta I(t)S(t)}{N} - \mu E(t), \tag{28}\]

\[\dot{I}(t) = \mu E(t) - \gamma I(t), \tag{29}\]

\[\dot{R}(t) = \gamma I(t). \tag{30}\]

In this model, suppose that we do not have data on the number of exposed population. Applying the NLO into the SEIR model, we obtain the estimated number of exposed case in Figure 3. It can be seen that the number is quite reasonable and it follows the pattern of the infections compartment. However, the estimated \(\mathcal{R}_t\) in Figure 4 is lower than in Figure 2. When it comes to a question which one is a better estimate? we argue that the estimated \(\mathcal{R}_t\) from Figure 2 is a better estimate since we cannot confirm that the estimated number of exposed population is close to reality.

44 Hasan, A.

Figure 3: Real-time data fitting using NLO from the SEIR model.

Figure 4: Estimation of \(\mathcal{R}_t\) using NLO from the SEIR model.

6. CONCLUSION

In this paper, we have presented a new approach to estimate the effective reproduction number \(\mathcal{R}_t\) of infectious diseases. The idea is to inject the discrete-time augmented SIR model with epidemiological data, such as active and removed case. The main challenge is to find a suitable observer gain, since there are no methods to solve the LMI systematically. However, once the observer gain is found, the method is comparable with EKF. The main advantage of using NLO is its stability compared to EKF. Furthermore, it does not require computation of inverse covariance matrices, which makes NLO more efficient. Extension to other compartmental models can be done, without guarantee on the quality of the estimated parameters. For general parameter estimation problems, we advise to use compartmental models where all data are available.