Robust control and data reconstruction for nonlinear epidemiological models using feedback linearization and state estimation

1.   Introduction

Epidemic control has been the target of intensive research in recent decades, and its significance became even more clear during the COVID-19 pandemic ^[1,2,3,4]. In most cases, the dynamical modeling of disease spread on the population level is done by using compartmental models on different levels of detail in terms of disease stages or population distribution in nonlinear ordinary differential equation (ODE) form ^[5,6,7]. It has to be noted, however, that the assumption of the homogeneity of the population and constant model parameters are often oversimplifying and may not allow the sufficiently accurate description of the epidemic process in complex environments. Therefore, it was already proposed in ^[8] that the probabilities of infection and death should depend on the different stages of the disease. To address the homogeneity problem, ^[9] gives a comprehensive overview on the modeling of infection mechanisms in heterogeneous networks using the notions of network science. In ^[10], the authors separate the infection dynamics for symptomatic and asymptomatic COVID-19-carriers, and also assume heterogeneous connectivity between people. With this modeling setup, it is shown that a more realistic picture of the pandemic in Italy can be reconstructed than by using classical models. Stochastic agent-based models can be very useful as well in handling time and/or space dependent probabilities, and also in tracking individuals and events and testing various conditions and restrictions ^[11,12].

Typical control goals are to keep the peaks and/or the sum of infected (or hospitalized) under predefined limits during a time interval to reduce the burden on the healthcare system. To handle the often conflicting goals and constraints related to economic and societal considerations, model predictive control can be an obvious and sufficiently realistic choice for control design ^[13,14,15]. However, general nonlinear model predictive control (recognizing the practical significance and applicability of the method) is a primarily numerical approach which sometimes does not give enough insight into certain fundamental properties and limitations of the solution. The sliding mode approach has also been successfully applied in the field of epidemic control ^[16,17].

Feedback linearization is typically applied for nonlinear systems where the models are known with sufficient precision, such as (electro) mechanical systems ^[18,19,20]. The robustness analysis of feedback linearized loops is known to be a difficult problem, and it has been the target of intensive research for several decades ^[21,22]. A critical point in such an analysis is checking the stability of the zero dynamics ^[23]. Therefore, feedback linearization is usually not the first choice for epidemiological or biological systems, where researchers are often faced with non-negligible model uncertainties and measurement data of questionable quantity or quality. However, there exist some successful applications in this area such as ^[24] and ^[25], where vaccination inputs were computed through the feedback linearization of nonlinear models. Moreover, the stability of the zero dynamics for any meaningful parameter values of a susceptible-exposed-infected-recovered (SEIR)-type model was also proved in ^[24].

For model-based feedback control, the estimation of non-measured states is a fundamentally important task ^[26]. Epidemic models (depending on the level of detail) typically contain several non-measured states such as the number of susceptible, exposed or asymptomatic people. The literature on the state estimation of epidemics is wide, where popular techniques include Kalman filtering, particle filtering, and maximum likelihood estimation ^[27,28,29]. We also highlight ^[30] and ^[31], where Kalman filter-based state estimation is successfully combined with feedback control. In ^[32], continuous time observers with guaranteed convergence were proposed for low dimensional epidemic models. A discrete time observer for SEIR models with efficiently computable stability conditions was presented in ^[33]. Feedback linearization combined with a state observer was proposed in ^[34] to compute a vaccination input which guarantees the asymptotic eradication of an epidemic described by a SEIR model. The convergence of the observation error to zero, the feasibility of the input, as well as the asymptotic stability of the desired state were proved.

The preliminaries of this study are ^[35,36], where promising initial results were obtained using feedback linearization assuming that full state information is available. However, this assumption is not realistic enough, since the majority of the state variables in epidemic models is not actually observed. Therefore, the main new contributions of this paper are the following. First, state estimation is involved in the control/reconstruction task in several different scenarios assuming significant modeling uncertainties as well. Second, the improved methodology is illustrated and tested on the data of another country (Sweden) besides Hungary.

The structure of the paper is as follows. Section 2 summarizes the compartmental models used for simulation and control design, respectively. The applied different control setup cases and the approach for state estimation are described in Sections 3 and 4, respectively. Section 5 contains the detailed computation results and the corresponding discussion, while the most important conclusions are given in Section 6.

2.   Process modeling

2.1.   Simulation model (SLPIAHRD)

To show the applicability of our control and estimation setup in realistic scenarios, we use different nonlinear dynamic models and parameters for the simulation of the epidemic process and for the estimation and control steps.

The epidemic simulation relies on a comprehensive eight-compartment model, abbreviated as SLPIAHRD (for susceptible, latent, presymptomatic, infected, asymptomatic, hospitalized, recovered, and deceased), which was designed with key characteristics of COVID-19 in mind. This model has been effectively used for data reconstruction and epidemic control in Hungary across multiple waves of the pandemic. Initially introduced in ^[15] and subsequently expanded in ^[37,38,39], the model accommodates various factors such as different latent periods (infectious and noninfectious), the impact of presymptomatic and asymptomatic infections on the spread, the waning of immunity, and countermeasures like symptom-based testing and quarantines. Additionally, it can be adjusted to account for new virus variants by gradually modifying the model parameters based on their prevalence in the population.

In this paper, we utilize the original model outlined in ^[15], which categorizes the population of $N$ individuals into 8 disjoint groups based on their roles in the epidemic spread. These groups include susceptibles (${\mathbf S}$), who lack protection against the virus; latently infected individuals (${\mathbf L}$), who carry the disease but are not yet spreading it; presymptomatic individuals (${\mathbf P}$), who can transmit the virus without showing symptoms and thus remain unaware of their infection ^[40]; and asymptomatic individuals (${\mathbf A}$), who eventually recover, possibly without ever realizing they were infected. Symptomatic individuals (${\mathbf I}$) may either recover (${\mathbf R}$) or require hospitalization. Hospitalized patients (${\mathbf H}$) can either recover or succumb to the disease (${\mathbf D}$). The compartmental structure and potential transmission pathways are illustrated in Figure 1, with a more extensive explanation available in ^[15].

Formally, the model is given by the following system of ODEs:

The state variables in the above model are the continuous numbers of individuals belonging to the compartments. For readability, we omit time arguments of the variables unless they are essential for understanding.

Our simulations rely on public epidemiological data from Hungary and Sweden during the second and third waves of the epidemic (15/08/2020 – 15/08/2021) to generate reference outputs for tracking. The data source is the Our World in Data (OWID) database ^[41], which keeps records of reported infections and hospitalizations provided by the responsible governmental institutions. The model parameters, shown in Table 1 and detailed in ^[38], are aligned with the dominant virus variants (original and Alpha) of this period, to ensure that the simulations are as realistic as possible. Parameter values for the specific variants and time periods in Hungary were obtained from relevant scientific literature, computed using standard parameter estimation or statistics-based techniques. Although they are slightly different, we decided to use the same parameters for Hungary and Sweden in this model to simplify the approach and to illustrate the robustness of the computations.

2.2.   Control design model (SEIR)

For feedback computation and part of the state observer design, we will use a largely simplified model compared to the previously described simulation model. This will be the well-known SEIR nonlinear compartmental model, which is widely used for describing epidemic processes, including COVID-19, despite its known limitations.

The model identifies four stages of infection for each individual, categorizing the population into four compartments: Susceptibles ($\mathbf {\overline{S}}$), who are unprotected, and may become Exposed ($\mathbf {\overline{E}}$) to the disease upon contact with Infected ($\mathbf {\overline{I}}$) individuals, who are already spreading the virus. Exposed individuals carry the disease, but take some time to begin spreading it, at which point they move to the $\mathbf {\overline{I}}$ compartment. Infected individuals eventually recover and move to the Recovered ($\mathbf {\overline{R}}$) compartment, gaining full immunity against future infections. The model operates with three parameters: the transmission rate $\beta$ of the virus (which is influenced by administrative countermeasures and thus considered a bounded manipulable input for the system), the latent period $k_2^{-1}$ (time spent in the Exposed compartment), and the infectious period $k_3^{-1}$ (time spent in $\mathbf {\overline{I}}$).

To consistently denote all variables, we use overline notation for the compartments in the control design model. With this notation, the differential equations of the continuous model can be expressed as:

Nominal values for the model parameters $k_2$ and $k_3$ were selected by approximately lumping the corresponding compartments of the simulation model based on their role in the epidemic spread ($\mathbf {\overline{E}} \approx {\mathbf L}$, $\mathbf {\overline{I}} \approx {\mathbf P} + {\mathbf I} + {\mathbf A}$), summing the time an individual spends in them. Additionally, we took the averages of the parameters for the original and Alpha virus variants, obtaining a time-invariant model. While this matching process could easily be refined, the applied model mapping is straightforward to use in real-world scenarios and underscores the robustness of our feedback approach to parameter uncertainty. Choosing a simpler model for control computation introduces model mismatch, a common issue in realistic setups where no single model can exactly describe the epidemic spread.

3.   Control design

We implemented the control of the simulation model in the following steps: estimation or computation (depending on control setup) of the control model's states based on output(s) of the simulation model, feedback linearization of the control model, computation of a control signal for the linearized control model, and, finally, application of this control input to the simulation model. Based on how we compute/estimate the states of the control model, we present three different control setups which are different in their degree of reality. This approach is also useful to compare how the results may be affected by the availability of data.

3.1.   Control setup 0: no estimation

We originally introduced control setup 0 in ^[36], although we had previously experimented with simple feedback-linearization of epidemic models in ^[35]. The structure of the feedback loop is shown in Figure 2, depicting a setup where no estimation step is carried out, and we assume full knowledge of the simulation model's state.

Our goal is to control the simulation model to track the reference signal $r(t)$ with the number of individuals capable of spreading the disease (i.e., minimizing $|r - ({\mathbf P} + {\mathbf I} + {\mathbf A})| = |r- \mathbf {\overline{I}}|$ by selecting the appropriate $\beta$). To compute the input, first we apply a modified feedback linearization based on [23, Chap 4.3], calculated for the control design model, as detailed below in Subsections 3.4 and 3.5. This technique requires the full state information of the SEIR model, which we compute directly from the states of the known SLPIAHRD states, and results in a transformed system (with state variables $z$, acting as a simple integrator cascade between the transformed input $v$ and output $y$. Next, we compute a simple full-state feedback gain for this linearized control model, which has a proportional-integral-derivative (PID) structure due to the obtained set of state variables.

3.2.   Control setup Ⅰ: partial state estimation

As knowledge of the full state information of the simulation model is problematic in a real-world application, we introduced control setup Ⅰ. in ^[42]. The feedback loop, as illustrated in Figure 3, contains one important modification compared to the previous approach, namely, we estimate the control model's state information $(\mathbf {\overline{S}}, \mathbf {\overline{E}}, \mathbf {\overline{I}})$ using an extended Kalman filter (EKF).

It must be noted that this setup assumes that the control model state $\mathbf {\overline{I}} = {\mathbf P}+ {\mathbf I}+ {\mathbf A}$ is a directly measurable output of the simulation model. Even though this might seem unrealistic at first as, e.g., the number of infected and virus-spreading asymptomatic individuals (${\mathbf A}$) is hardly measurable, during the COVID-19 pandemic, several methods were proposed and tested for the estimation of these variables ^[43]. Thus, we assume in this scenario that a reliable estimation for $\mathbf {\overline{I}}$ is already available.

Apart from this modification, we use the same extended feedback-linearization technique as in control setup 0, which again results in a PID-like control structure for the transformed, linearized model.

3.3.   Control setup Ⅱ: full state estimation

In order to overcome the limitations discussed above (i.e., the need for an extra estimation step) and to take into account the effect of that estimation's error on the whole system, we created control setup Ⅱ. While similar to setup Ⅰ., in this case we assume the only measurable output of the simulation model to be the daily number of hospitalizations (${\mathbf H}(t)$). We apply an EKF constructed using the same model structure as the simulation model, but (not to move away from reality) having some of its parameters mismatched; we estimate the states of the simulation model, and finally we 'compute' the states of the control model by merging and matching the compartments of the simulation model as described above in the modeling section. Setup Ⅱ can be seen in Figure 4.

3.4.   Feedback linearization and asymptotic output tracking

Based on ^[36], we introduce the technique for feedback linearization following [23, Chap 4.3].

We derive the equations for nonlinear single input single output systems, written as:

where $u\in\mathbb{R}$ is the system's input, $x\in\mathbb{R}_+^n$ is the state, and $y\in\mathbb{R}$ is the observed output. We assume that $f$, $g$, and $h$ are well-defined and sufficiently smooth functions. Furthermore, we define a reference signal $r(t)$, as our ultimate control objective is calculating an input $u$ such that $y$ (asymptotically) converges to $r$.

Using the well-known definition of the Lie derivative of a scalar-valued function $h(x)$ along the vector field $f(x)$ as

and the relative degree of a system being $\rho$ if

a state transformation $z$ can be defined as:

Applying the nonlinear state feedback:

for the transformed system results in the closed-loop (from input $v$ to output $y$) acting in a way equivalent to a cascade of $\rho$ integrators:

Using the fact that such an integrator cascade can be stabilized by a simple state feedback, we can construct a stabilizing input driving the output to $0$ in the form of:

where the constants $K = (K_0, K_1, ..., K_{\rho-1})$ can be arbitrary values that make $(A-BK)$ a stability matrix. In fact, the value of $K$ can be obtained by applying any linear stabilization method to the system (3.5) such as pole placement, linear–quadratic regulator (LQR), or any other more advanced feedback approach.

Now, given a reference signal $r(t)$ to track with the output of the system, we can define the error output $e(t) = y(t) - r(t)$ and the corresponding state transformation as:

Analogously to (3.6), we can stabilize the system using the simple state feedback. Thus, the application of:

will result in $\lim\limits_{t\rightarrow \infty}y(t) - r(t) = \lim\limits_{t\rightarrow \infty}e(t) = 0$, if $K$ is chosen such that $(A-BK)$ is a stability matrix [23, Chap 4.5]. For the sake of completeness, substituting $v$ into the input mapping Equation (3.4), we compute the asymptotic reference tracking input, which should be applied to the original system as:

3.5.   Robust feedback linearization and output tracking

As it was previously shown in ^[36], the robustness of the control loop may be significantly improved by the introduction of an extra integrator step into this cascade, given the same model structure as in Equation (3.1), and effectively increasing the relative degree $\rho$ by 1.

We define the new, transformed output $\overline{y}$ as $\dot{\overline{y}} = h(x) - r$, and modify the transformed system as:

Analogously to the original case, we achieve asymptotic output tracking with a full state feedback gain $K = (K_0, K_1, ..., K_{\rho})$, where $K$ is an arbitrary constant making the augmented tracking error dynamics asymptotically stable.

By introducing the extra integrator into Equation (3.4) and using the feedback $v = -Kz$ (Equation (3.8)), we can conclude that the asymptotic output tracking input can be calculated as:

It is clearly visible, that the application of the technique to our SEIR control design model requires the state information $x = (\mathbf {\overline{S}}, \mathbf {\overline{E}}, \mathbf {\overline{I}}, \mathbf {\overline{R}})$ at any given time.

4.   State estimation

In this section, we will describe the state estimation approach with the corresponding control input computation.

4.1.   Extended Kalman filter for state estimation

For the state estimation, we apply a discrete-time nonlinear system model formally described as:

where $w$ and $v$ are additive process and output noise terms, i.e., independent random variables with Gaussian distribution, having zero mean and covariance matrices $Q$ and $R$, respectively.

Around this model, we construct a state estimator using the discrete-time EKF algorithm. The EKF has two key attributes:

● $\hat{x}_{k|j}$ represents the state estimated for time $k$ at time $j$ ($j \leq k$), initialized as $\hat{x}_{0|-1} = \mathbb{E}_{x_0}$, which is the expected value based on prior knowledge of the system, without any measurements taken.

● $P_{k|j}$ denotes the state estimation error covariance matrix (i.e., the uncertainty of the estimation) for time $k$ calculated at time $j$, initialized as $P_{0|-1} = \mathbb{E}{(x_0 - \hat{x}_{0|-1})(x_0 - \hat{x}_{0|-1})^T}$ and also based on prior knowledge.

At each time step $k = 0, 1, 2, ...$ the EKF algorithm executes two main steps:

1. Update (correction) of the current state estimate $\hat{x}_{k|k}$ and current uncertainty $P{k|k}$, based on the previously predicted state $\hat{x}_{k|k-1}$, previously assumed uncertainty $P_{k|k-1}$, and the current measurement $y_k$.

2. Prediction of the state at the next time instance $\hat{x}_{k+1|k}$ using the currently available state estimate $\hat{x}_{k|k}$ and the state transition function of the estimated model, i.e., $f(x_k, u_k)$ from Equation (4.1).

Due to the fact, that functions ($f, h$) can be nonlinear, the EKF extends the operating principle of the linear Kalman filter algorithm by approximating the state transition and output matrices ($A$ and $C$) of a linear time-invariant (LTI) system using the Jacobians $A_k = \left.\frac{\partial f}{\partial x}\right|_{\hat{x}_{k|k}}$ and $C_k = \left.\frac{\partial h}{\partial x}\right|_{\hat{x}_{k|k-1}}$. Using this notation, the above steps can be computed as:

1. Correction:

2. Prediction:

4.2.   Control input computation with estimated SEIR states

First of all, to obtain a suitable model (corresponding to the structure (4.1)), we introduce the notations $x_k = (\mathbf {\overline{S}}_k \; \mathbf {\overline{E}}_k \; \mathbf {\overline{I}}_k \; \mathbf {\overline{R}}_k)$, $y_k = \mathbf {\overline{I}}_k$, $u = \overline{\beta}$, and discretize our control design model (SEIR) using a simple forward Euler method and time step $dt$:

As we need the Jacobians for substitution into 4.2 and 4.3, we compute them as:

Second, we derive the state transformation and the required output mapping as detailed in (3.10) - (3.11). We compute the required Lie-derivatives for the continuous system (3.1) as follows

and conclude from the definition (3.2) that $\rho = 2$. However, as our technique for robustifying the controller defines a transformed output (3.10) and effectively increases the relative degree by 1, we have to consider the relative degree to be $\rho = 3$, and thus construct the input mapping required for asymptotic output tracking as:

creating a system equivalent to:

Now, by introducing the notation $\hat{x}_k = (\hat{x}_{1, k} \; \hat{x}_{2, k} \; \hat{x}_{3, k}) = (\hat{ \mathbf {\overline{S}}}_k \; \hat{ \mathbf {\overline{E}}}_k \; \hat{ \mathbf {\overline{I}}}_k)$ for the state information computed by the EKF, we can substitute the estimated state values into the above formulae, resulting in:

Moreover, applying the input required for the asymptotic output tracking, we obtain:

Thus, the obtained feedback has the following PID structure:

Here, we must take two important notes. First, the formulae for feedback linearization requires the derivatives of signals $r$ and $y$, which, due to the EKF estimation, are available only at discrete time instants. Therefore, we use the numerical derivatives, computed as $\dot{y}_k = (y_{k+1} - y_{k})/dt$, $\dot{r}_k = (r_{k+1} - r_{k})/dt$, and provide a two-times continuously differentiable reference signal for the system. If the reference signal is not differentiable, we can ensure the required smoothness by basic preprocessing techniques (see, e.g., ^[37]). Second, we emphasize, that despite the PID-like structure of (4.11), the constants of the state feedback can be computed using a simple pole placement (computed for the system presented in Equation (4.8)). In our experiments, we chose the desired poles of the controlled system to be $-(0.339, 0.407, 0.475)$, values obtained empirically during the creation of our first modified feedback-linearization based controller in ^[35]. For the EKF, we used $dt = 0.5$ (days) and assumed $Q = \text{diag}(1, \; 1, \; 1, \; 1)$, $R = 0.01$, and initial estimation error covariance $P_{0|-1} = 0.01$. We chose initial state different from the true state value as follows: $\hat{x}_{0|-1} = {\mathbf N}\cdot((1-10^{-5}), \; 5\cdot10^{-6}, \; 5\cdot10^{-6}, \; 0)^T$, inducing a short transient period at the beginning of the simulation. The estimator used the same (possibly mismatched) $k_2$ and $k_3$ parameters as the SEIR-based controller.

4.3.   Control input computation with estimated SLPIAHRD states

For control setup Ⅱ, we apply the discrete state estimation to the detailed simulation model given in Equation (2.1), and compute the estimated control model states by merging the respective compartments. Compared to the previous case, instead of substituting Equation (4.5) into Equation (4.1), we derive the discrete transition and measurement functions (using the Euler method), as well as their Jacobi matrices for the SLPIAHRD model as follows:

For the controller, we used the same parameter combinations as above. For the EKF, we used $dt = 0.5$ (day) and assumed $Q = \text{diag}(10^{-6}, \; 0, \; 0, \; 10^{-6}, \; 10^{-6}, \; 10^{-6}, \; 0, \; 0)$, $R = 1$, and initial estimation error covariance $P_{0|-1} = 10^{-3}$. For the initial state, we used $\hat{x}_{0|-1} = {\mathbf N}\cdot(0.9999869, \; 10^{-5}, \; 10^{-6}, \; 10^{-6}, \; 10^{-6}, \; 10^{-7}, \; 0, \; 0)^T$. The estimator used the same parameters as the simulation model, apart from the mismatched parameters $\gamma$ and $\delta$, detailed in the next section.

4.4.   Synthesis of continuous and discrete-time subsystems

As we are presenting both discrete-time and continuous models for the control and estimation steps in this paper, we would like to clarify the connection of all the components in this subsection. We carry out all computations regarding the feedback-linearization and control on the continuous SEIR model. The output of the simulation model, which is continuous as well, is sampled for the discrete-time EKF, which in turn produces discrete state estimations. These estimations are used for computing the input of the system by substituting the obtained values into the state-feedback and input-mapping formulae. The obtained discrete-time input signal is fed back into the simulation model using zero-order hold. The derivatives of the discrete time signals (where needed) are approximated by the differences of the consecutive data samples.

5.   Results and discussion

To illustrate and test the operation of the previously presented control setups, we carried out five experiments for the control and reconstruction of historical epidemic data, altering some parameters of the control model and the EKF (compared in Table 2).

We chose Sweden's and Hungary's estimated infection data from the interval 15/08/2020 – 15/08/2021 as a base for our reference signals, an interval containing two epidemic peaks for both countries. For the sake of simplicity, in both cases, we assumed the same model parameters as in the case of Hungary (as shown in Table 1, and obtained through parameter estimation and statistical methods detailed in ^[15]). We set the parameters of the simulation model as slowly-varying with the weighted average of the parameter values being used according to the dominance of the respective virus variant at the given time (detailed in ^[35]).

The simulations were run in Matlab/Simulink using the ode45 solver with relative tolerance $10^{-10}$, on a computer with processor Intel Core i7-8565u (1.8–4.6 GHz) and 16GB RAM, requiring between 45–90s per simulation (appr. 1.5 hours per control setup).

To illustrate the robustness of the different scenarious, we altered the parameters of the model(s) used for control and estimation in the different controller setups. The comparison of the applied parameter errors in the different scenarios is shown in Table 2.

5.1.   Contol setup 0: robustness without state estimation on Sweden's data

In this experiment, we used control setup 0 (Section 3.1) for reference tracking control of the SLPIAHRD model using Sweden's data. This setup assumes that all state variables of the simulation model are known at all times, and we computed the states of the control model through the lumping of the respective compartments.

We altered the parameters $k_2$ and $k_3$ on a $10 \times 10$ grid (i.e., $(k_2, k_3) \in \{(0.2\cdot m_2k_2, 0.2\cdot m_3k_3) | m_2, m_3 = 1\dots10 \}$ applying $20-200\%$ relative error compared to the nominal value. In total, we ran 100 simulations, the results (nominal trajectories and the computed standard deviation) being illustrated in Figures 5–7. We can observe, that the reference signal is closely followed by the controlled system's output even for significant changes of the parameter values, showing the robustness of the proposed controller setup. Additionally, Figure 7 confirms that even for the simulations with high inaccuracies, there is no significant change in the computed input, which apart from the transient period at the beginning remains in a feasible range for the observed epidemic.

5.2.   Control setup Ⅰ: robustness with SEIR state estimation on Sweden's data

In this scenario, we examined the performance of the system in control setup Ⅰ. (Section 3.2), using an EKF constructed around the SEIR model and Sweden's infection data as a reference signal.

Similarly to the previous scenario, we altered the relative error of parameters $k_2, k_3$ for the control model (affecting both the state estimation and the feedback linearization) on a $9\times9$ size grid $\{0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 2.5, 3\}^2$ (the middle of the grid $(1, 1)$ being the nominal parameter values).

The simulation results are illustrated in Figures 8–10. Figure 8 compares the estimated and simulated values of susceptible and recovered individuals: compartments ${\mathbf S}$ and ${\mathbf R}$ of the simulation model plotted alongside the estimated compartments $\hat{ \mathbf {\overline{S}}}$ and $\hat{ \mathbf {\overline{R}}}$ of the control design model, as computed by the EKF algorithm.

Similarly, Figure 9 displays the same comparison for estimated and simulated $\mathbf {\overline{E}} = {\mathbf L}$ and $\mathbf {\overline{I}} = {\mathbf P} + {\mathbf I} + {\mathbf A}$ compartments. It should be noted, however, that deviations of the simulated compartments are scarcely visible in any of the figures (being several magnitudes lower than other plotted quantities). Finally, Figure 10 contains the computed input fed into the simulation model.

Qualitatively matching the results described in ^[42], where an experiment with similar setup but different aim and data (i.e., stricter reference curve, to evaluate performance in case of more intense epidemic control based on Hungary's data), we can observe that the controller performed in an acceptable range for all examined parameter combinations.

First, as we can see looking at the continuous and dashed lines in Figures 8 and 9, the EKF produces precise estimation of states in the case of nominal parameters, regardless of the fact that these nominal values were obtained by the simple intuitive matching of compartments based on their functionality, rather than any kind of formal parameter computation (e.g., by fitting through optimization). As for the mismatched parameters, expectedly they produce significant estimation errors, especially in compartments which do not affect the system's output (like $\mathbf {\overline{R}}$). On the other hand, simulation model trajectories show minimal deviance (approx. 2 magnitudes lower, considering the mean tracking error). The input visible in Figure 10 also shows small variation, emphasizing the error-prone operation of the proposed control law and EKF-based estimation setup.

5.3.   Control setup Ⅱ: robustness with SLPIAHRD state estimation on Sweden's data

In this experiment, we examined the performance of the system in control setup Ⅱ (Section 3.3). Differently from the previous case, this time we used an EKF constructed around the SLPIAHRD model with mismatched parameters. Inspired by our experience during the Hungarian COVID-19 waves, we chose the uncertain parameters to be $\gamma$ (probability of symptomatic infection) and $\delta$ (relative infectiousness of asymptomatic individuals), as these are the hardest to accurately determine using statistical methods. Similarly to the previous scenario, we also altered the parameters $k_2, k_3$ of the control model.

We ran 81 simulations, altering $k_2, k_3, \gamma, \delta$ on a $3\times3\times3\times3$ sized grid: $k_2, k_3 \in \{0.5, 1, 1.5\}$ of their nominal values. We used $\gamma \in \{0.25, 0.6, 0.99\}$ and $\delta \in \{0.3, 0.75 0.99\}$ (these being parameters showing probability, thus having a feasible range $[0, 1]$; the middle value for both parameters corresponds to their nominal value, as shown in the parameter table).

The results of the computations can be seen in Figures 11–17. From Figures 11 and 12 we can see that the applied parametric uncertainties cause significant deviations in the state variable values of the simulation model. Specifically, the presymptomatic and asymptomatic numbers are affected, while the number of hospitalized people shows a relatively smaller variation. Naturally, these uncertainties have a huge impact on the state estimations of the EKF which uses the nominal model, as it is shown in Figures 13 and 14. The deviation of the infected compartment is shown separately in Figure 15. We can observe that the uncertainty is the largest at the most critical point, i.e., around the peak of the first massive wave. However, as Figure 16 shows, this deviation can be greatly reduced using the proposed feedback scheme, even with a highly uncertain state estimation. Finally, Figure 17 shows the computed transmission rate. Although the deviation of $\beta$ is naturally larger than in the previous scenario, the estimation quality is still acceptable in this realistic situation.

From the results of control setup Ⅱ, two significant conclusions can be drawn. First, as it is visible in Figure 15, asymptotically tracking a reference for a non-measurable 'output' of the simulation model (${\mathbf P}+ {\mathbf I}+ {\mathbf A}$) estimated using mismatched parameters will result in the estimation error directly appearing in the reference tracking performance - hence the high deviation in the $\mathbf {\overline{I}}$ compartment. However, as shown in Figure 16, the controller – in spite of the serious state estimation errors, combined with the also mismatched control model parameters distorting the input mapping – can effectively minimize the difference between the estimated model output and the reference signal. Not surprisingly, as we can observe in Figures 11 and 12, this estimation bias (effectively causing the simulation model to follow different references in the simulations) can be seen in the deviation of the model states as well. However, as the much higher deviance values in Figures 13 and 14 show (especially in the case of compartments ${\mathbf I}$ and ${\mathbf A}$ highly affected by the uncertain parameters), the controller even in this case greatly reduces the spread caused by the incorrect state estimation, and maintains all important qualitative properties of the trajectory to be tracked. As a second observation, the small deviance values of the input $\beta$ (visible in Figure 17) must be emphasized, showing that the incorrect controller and estimation parameters do not influence the input significantly. Consequently, even though imperfectly estimated output used for reference tracking can naturally decrease the tracking performance, we haven't experienced problematic inputs (abrupt changes, out-of-bound values, etc).

5.4.   Testing: reconstruction of Hungarian data using SLPIAHRD state estimation

As it is summarized e.g., in ^[44], such a reference tracking scenario can also be viewed as an inversion of a dynamical model. By trying to track / reproduce the actual output of the epidemic process (as observed during the pandemic) with our model, the model's input (desired reproduction rate, computed by the controller) effectively becomes an estimation of the actual reproduction rate of the epidemic. In this way, the reproduction of model trajectories from the output can be transformed into a control problem, solved by the presented feedback-linearization technique.

In this scenario, we use control setup Ⅱ. to track the transmission rate and the number of active infections measured (or more precisely, estimated from multiple factors) in Hungary between 15/08/2020 – 15/08/2021. Similarly to the previous scenario, we alter the parameters $k_2, k_3, \gamma, \delta$ on the same grid, totaling 81 simulations.

Unfortunately, we cannot use $\beta$ directly to compare our results to those published in the relevant literature, thus we compute the effective reproduction number of the epidemic from it. As explained in ^[7], the time-varying effective reproduction number of an epidemic corresponds to the number of susceptible individuals an infectious person will infect on average, given the circumstances at the given time, and can be derived for the SLPIAHRD model as:

The effective reproduction number obtained as a result of this simulation can be seen in Figure 18, along with the effective reproduction number estimated for Hungary by the OWID ^[41] and by the ÁtlóTeam ^[45], both using a classical statistical approach.

For comparison, we also include our own previous estimation published in ^[36], created using the same robust controller setup, but without any state estimation (assuming all states of the simulation model to be directly measurable). In that case, we used the SEIR model to compute the effective reproduction number from $\beta$, derived as:

The results of those simulations are shown in Figure 19, compared again with the results of the OWID and the ÁtlóTeam.

5.5.   Testing: reconstruction using public data of Sweden

Similarly to the Hungarian reconstruction scenario above, for the sake of completeness and validation we also compute the effective reproduction number for Sweden, using control setup Ⅱ (EKF based on the SLPIAHRD model, with altered $k_2, k_3, \gamma, \delta$ parameters). In this case, we could only compare our results to the reproduction number computed by OWID. As it is visible in Figure 20, following a brief transient period (while the state estimator and the controller converge), we could track the reproduction number with high accuracy (and low deviation, in spite of the serious parameter mismatch).

5.6.   Comparison: Hungarian and Swedish data

During the COVID-19 outbreak, Hungary and Sweden opted for drastically different measures for dealing with the high number of infections and the burden it put on the country's healthcare system. Hungary initially opted for a strict suppression strategy, practically eliminating the first epidemic wave, and gradually releasing restrictions later, resulting in serious peaks in the second (around 10/20-11/20) and third (02/21) waves. Sweden on the other hand adopted a more permissive intervention strategy in the beginning, trying not only to mitigate the effects of the pandemic, but also minimizing the negative effects of the restrictions on the economy. As a result, more serious first and second waves were followed by smaller outbreaks. This difference of strategy appears in the reconstructions of the pathogen's reproduction numbers in the two countries, as illustrated in Figure 21.

6.   Conclusions

The goal of this paper was to study the possibilities of feedback linearization-based epidemic control combined with state estimation in the presence of parametric and modeling uncertainties. For this, we first proposed a reference tracking control scheme for a nonlinear 8-compartment epidemic model. The system input is the transmission rate, while the controlled output is the number of infected people within the population. The model used for control design is a simplified 4-compartment SEIR system where the state variables are mapped to that of the more detailed model. The control law is based on feedback linearization and asymptotic output tracking in a servo configuration. The state information for the feedback is computed by an extended Kalman filter. Two control setups containing state estimation were tested: in the first one, the state estimator was based on the simplified model, while in the second scenario, the EKF also used the detailed simulation model. The robustness of both scenarios were studied for several different parameter combinations in a range of 10% –300% of their nominal values. Although the state estimation errors were significant in the case of uncertain parameters and model mismatch, the tracking performance was good in all studied cases. The proposed approach was also illustrated as a possible tool for retrospective data (especially, transmission rate) reconstruction. From this respect, we can say that the applied model is 'well invertible' in the sense that the variance of the transmission rate was small in all computations. In summary, the most important message of the study is that the combination of a state estimator and an appropriately designed nonlinear state feedback can be successful even if the control model contains significant modeling and/or parametric uncertainties. This is promising from the point of view of possible applications.

The main limitations of the study are the following. First, the number of infected people was used as observed 'performance' output, while it is known that it is challenging to determine it in real time. However, there are several statistical and dynamical-model-based methods to refine this number from, e.g., the number of hospitalized people and/or from wastewater samples ^[46,47,48]. Second, the manipulated input is assumed to be a variable which can be changed in each sampling instant. This makes the proposed methodology more suitable for data reconstruction as it was illustrated in the case studies. However, the generated inputs and trajectories can be used as feasible starting points for more realistic approaches such as model predictive control (MPC). Moreover, MPC generally also requires state estimation in a similar setup, and an important message of this study is that an appropriate feedback may still be operational with high estimation uncertainties caused by model mismatch. Additionally, we have to note that the properties of the entire control loop were not mathematically analyzed. In this respect, the most relevant results from the literature mentioned in the Introduction are ^[32,33,34]. However, model uncertainties are not considered in these articles.

Further work will be focused on the application of more advanced robust controllers on the linearized system as well as on the translation of the computed input to specific epidemiological measures using our agent-based simulator ^[11].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was partially supported by the Hungarian National Laboratory of Health Security (RRF-2.3.1-21-2022-00006) and the Post-COVID-2022 grant of the Hungarian Academy of Sciences. B. Cs. acknowledges the support of the ÚNKP-23-3-Ⅱ-PPKE-10 New National Excellence Program and the EKOP-24-3-Ⅱ-PPKE-135 University Research Fellowships of the Ministry for Culture and Innovation from the source of the National Research, Development and Innovation Fund (NKFIH). G. Sz. acknowledges the support of the grant NKFIH-OTKA 145934. The authors thank the anonymous Reviewers for their constructive comments and suggestions. The authors acknowledge the assistance of Kiley Kesterson for their careful and meticulous proofreading of the manuscript.

Conflict of interest

All authors declare no conflicts of interest in this paper.