Reconstruction of Epidemiological Data in Hungary Using Stochastic Model Predictive Control

In this paper, we propose a model-based method for the reconstruction of not directly measured epidemiological data. To solve this task, we developed a generic optimization-based approach to compute unknown time-dependent quantities (such as states, inputs, and parameters) of discrete-time stochastic nonlinear models using a sequence of output measurements. The problem was reformulated as a stochastic nonlinear model predictive control computation, where the unknown inputs and parameters were searched as functions of the uncertain states, such that the model output followed the observations. The unknown data were approximated by Gaussian distributions. The predictive control problem was solved over a relatively long time window in three steps. First, we approximated the expected trajectories of the unknown quantities through a nonlinear deterministic problem. In the next step, we fixed the expected trajectories and computed the corresponding variances using closed-form expressions. Finally, the obtained mean and variance values were used as an initial guess to solve the stochastic problem. To reduce the estimated uncertainty of the computed states, a closed-loop input policy was considered during the optimization, where the state-dependent gain values were determined heuristically. The applicability of the approach is illustrated through the estimation of the epidemiological data of the COVID-19 pandemic in Hungary. To describe the epidemic spread, we used a slightly modified version of a previously published and validated compartmental model, in which the vaccination process was taken into account. The mean and the variance of the unknown data (e.g., the number of susceptible, infected, or recovered people) were estimated using only the daily number of hospitalized patients. The problem was reformulated as a finite-horizon predictive control problem, where the unknown time-dependent parameter, the daily transmission rate of the disease, was computed such that the expected value of the computed number of hospitalized patients fit the truly observed data as much as possible.


Introduction
The recent and still ongoing COVID-19 pandemic has brought unprecedented challenges for most countries to protect human lives and operate the economy and society at an acceptable level at the same time [1,2]. In supporting the related difficult decisions, dynamical modeling of the epidemic process has had a key importance in all developed societies [3]. Depending on the modeling goal, a wide range of computational techniques is available to describe, forecast [4], and even control the epidemic process [5]. Here, we can only mention a few selected references from the related extensive literature. The majority of the modeling solutions are based on deterministic compartmental description derived from susceptible-exposed-infected-recovered (SEIR)-type models [6][7][8]. Due to the rapidly increasing computing power, agent-based models have also become popular for modeling epidemics and studying control possibilities [9,10]. Logistic wavelets were used in [11] to compute the possible cumulative number of people infected with COVID. The application of artificial intelligence and machine learning has also been successful in COVID-related modeling and prediction [12][13][14].
The online tracking of informative epidemic parameters and variables proved to be essential for the continuous monitoring and evaluation of the situation, making predictions, or planning interventions. A key parameter in such analyses is the time-varying reproduction number (R t ) which is the population-level transmission potential of the disease at time t, which is known to be non-trivial to estimate [15]. Therefore, several different computational approaches have been proposed to track this quantity for epidemic outbreaks [16,17]. In [18], the authors fit a stochastic compartmental model defined by a random walk to estimate the time-dependent reproduction number of the COVID-19 epidemic in Wuhan from publicly available datasets. A discrete-time Hawkes process was used for the estimation of R t in [19], which allowed the detection of events such as restrictions and their relaxation. The estimation of other non-measured variables such as the number of people in latent or asymptomatic stages is also a relevant problem in data reconstruction [20]. In [21], a state estimator with proven convergence was proposed to implement model predictive control (MPC) satisfying complex constraints for an eight-compartment model of the COVID-19 pandemic. Essentially the same model was used in [22] for an inversion-based estimation of R t from Hungarian data. Similar to [21], an MPC approach was presented in [23] to determine optimal social distancing rules (and hence, R t ) to mitigate the epidemic.
The majority of the data analysis approaches use primarily the daily infected data possibly together with the recovery statistics. However, it is well known that only a fraction of the true cases are detected, which also depends on the testing intensity [24][25][26]. Moreover, the recording of recoveries is also often not immediate and sometimes not precise enough as well. Therefore, similar to [22,27], we used the official data on the daily number of hospitalized people in Hungary, assuming that current testing rules and protocols in hospitals give sufficiently reliable and timely information. The basic system theoretic idea behind our proposed solution is that the transmission parameter β, which is closely related to R t , can be considered as an input of a nonlinear system describing epidemic spread, and its estimation can be traced back to a trajectory-tracking control problem, where the output to be tracked is the true reported number of hospitalized people. A straightforward choice for the solution is MPC, where we can computationally handle model nonlinearities and parametric uncertainties [28].
The theory of optimization-and prediction-based simultaneous state and parameter estimation is well founded in the literature [29][30][31][32], and it is applied in a wide spectrum of sciences, e.g., in geosciences [30], in medicine [33], in agriculture [34,35], in aerospace [36], or in meteorology [37]. The classical approaches [30,37] use variational principles and consider continuous-time models. Due to its simplicity and transparency, the MPC approaches with discrete-time model descriptions are widely used to solve optimal filtering problems, e.g., [32][33][34][35] addressed model predictive data assimilation, i.e., optimal state/parameter reconstruction to minimize the deviation between the measurement and model output. In [36], an optimal design parameter was computed through a constrained MPC for a small satellite system, which provides constraint satisfaction during operation time. Furthermore, optimal dosing of cancer chemotherapy was addressed in [33] by solving a predictive control problem with joint state and parameter estimation.
In general, the nonlinear MPC (NMPC) approaches result in a cumbersome optimization problem, especially when the model equations are stochastic. However, the available sequential convex programing approaches [38,39], the algorithmic differentiation techniques [40,41], and the numerical software tools [42,43] exploit the special structure of a typical MPC problem and provide an efficient toolkit to solve the nonlinear problems precisely in a reasonable time. Among the several approaches to cope with stochastic dynamics [44], we mention two groups of techniques, which are popular in control theory. First, the particle-based approaches [45][46][47][48][49] with scenario trees allow coping with general (not necessarily Gaussian) models. Secondly, the tube-based approaches [50][51][52][53] approxi-mate each predicted state and input by a Gaussian distribution. Therefore, the dynamic equations in these references are recast as a deterministic mean-variance recursion.
In this paper, we propose a generic optimization-based approach to reconstruct epidemiological data through the approximation of unknown time-dependent quantities (such as the state, unknown input, or parameter) of a class of discrete-time nonlinear stochastic dynamical models by a sequence of Gaussian distributions. The problem was formulated as a single stochastic NMPC (SNMPC) computation. However, the solution of an SNMPC over a relatively large prediction horizon is challenging. Therefore, the problem was solved in multiple steps. First, we approximated the expected input and state trajectories through a nonlinear deterministic problem. If the expected values of the unknown quantities are fixed, the variance matrix of the joint distribution is well defined. An optional state feedback gain computation is proposed to reduce the solution's conservatism, i.e., the estimated standard deviation of the unknown quantities. Finally, the computed mean and variance values serve as an initial solution for the SNMPC problem, which ensures a fast convergence for the nonlinear optimization.
The paper is organized as follows. First, we describe the applied compartmental model for the COVID-19 epidemic spread in Section 2. Then, we introduce our computational approach in two steps in Sections 3 and 4. The numerical results with a discussion are presented in Section 5.

Notations and Abbreviations
All random variables are distinguished from the deterministic variables in notation by the accent•. Namely,x is a random variable, whereas, x is deterministic. Whenx is normally distributed with expected value µ and variance Σ, we write thatx ∼ N (µ, Σ). Whenx ∼ N (µ, σ 2 ) is a scalar-valued Gaussian variable, the confidence intervals µ ± 1σ = [µ − 1σ, µ + 1σ] and µ ± 2σ = [µ − 2σ, µ + 2σ] of probability levels 68.2% and 95.4% are called the 1σ and 2σ confidence intervals, respectively. The value of a time series x : {0, 1, . . . } → R n at time instant k is denoted by x k . Each constant or variable, which denotes a given number of people, is denoted by a boldface letter, e.g., N constitutes the number or people in a community or the population of a country. When A ∈ R n×m is a matrix, He{A} stands for A + A, where A is the transposition of A. Let x = (x 1 , x 2 , . . . , x n ) denote x = x 1 x 2 . . . x n . The matrix-valued functions ∂ f ∂x : R n+m → R p×n and ∂ f ∂(x,y) : R n+m → R p×(n+m) (with arguments (x, y) ∈ R n+m ) denote the Jacobian of function f : R n+m → R p with respect to (w.r.t.) x ∈ R n and (x, y) ∈ R n+m , respectively. Furthermore, the value of ∂ f ∂x at x 0 ∈ R n and y 0 ∈ R m is referred to as ∂ f ∂x (x 0 , y 0 ) ∈ R p×n . The Euclidean norm of a vector x ∈ R n is x , whereas the weighted norm of x w.r.t. the symmetric and positive definite matrix Q is referred to as x Q , namely x 2 Q = x Q x. Finally, let I b a = {a, a + 1, . . . , b} denote the set of integers between a and b.

Compartmental Model of the Spread of the COVID-19 Epidemic in Hungary
In this section, we present the compartmental model describing our knowledge on the dynamics of disease spreading.

Transitions between the Phases of the Disease
To capture the spread and the evolution of the COVID-19 epidemic, we considered a modified version of the compartmental model introduced in [21]. This model divides the population of N individuals into eight classes, representing the different stages of the illness. The compartments of the model correspond to the following subsets/groups of the population: susceptible individuals (S), infected people in the latent (L) and the presymptomatic (P) phases, infected people in the main sequence of the disease (A, I), infected people who need hospital treatment (H), and finally, the recovered (R) and deceased (D) people. The main phase of the disease is further divided into those who remain asymptomatic (A) and those who produce symptoms (I).
In this work, the model of [21] was complemented with a new compartment (U) comprising all individuals who became immune through vaccination. The members of this compartment are governed by the daily number of vaccinated people (V). A discrete-time (DT) version of the augmented continuous-time model was obtained through the explicit Euler method with a 1 d-long sampling period. The dynamic equations are given as follows: The transitions between the compartments are illustrated in Figure 1. It is worth mentioning that the recovered's compartment R contains the recovered people who are not yet vaccinated. To represent all the recovered people including the vaccinated, we can consider the following additional equation: which clearly does not affect the dynamics of other states. The detailed worldwide [54] and local serological [55] and dynamical [7,20,56,57] analysis results provide estimates for the average lengths of the phases of the illness and the probabilities of transitions between the compartments. These parameters are aligned with the Hungary-specific data and were presented in detail in [21]. After infection, the latent period of the disease (L) usually lasts approximately α −1 = 2.5 d. This period is followed by a pre-symptomatic phase (P) of ζ −1 = 3 d. A person in the main sequence of the disease (I or A) remains infectious for about ρ −1 The empirical probability of producing symptoms in the main sequence is γ = 0.6; furthermore, an η = 0.076 portion of the symptomatic cases require hospitalization. The average length of a hospital treatment is λ −1 = 10 d. A hospitalized patient dies with a probability of µ = 0.205, or recovers. The recovered people are assumed to be immune to reinfection. We assumed that the disease is transmitted by the members of compartments P, A, and I, such that the relative infectiousness of the asymptomatic individuals (A) is δ = 0.75, compared to those who produce symptoms (I). The transmission rate β of the disease is the most prominent parameter of the epidemic spread, which is typically time-dependent. The nominal values of the above-mentioned model parameters and their assumed uncertainty are collected in Table 1. For simplicity, each uncertain parameter was assumed to have a normal distribution, such that its nominal value (µ) is the expectation and its uncertainty (±a%) gives the 2σ confidence interval with the standard deviations σ = aµ/200. Here, we examined the evolution of the epidemic in a fixed time window between 1 March 2020 (k = 0) and 30 June 2021 (k = T). This interval contains the first three waves of the epidemic in Hungary. Subscript k ∈ {0, 1, . . . , T} denotes the number of days elapsed in the given time window of length T + 1.

Vaccination Model
For simplicity, our vaccination model assumed that only susceptible and recovered people are eligible for vaccination, and we neglected those rare cases when the shot is given during an unidentified infection. Based on the serological test data [58], our model assumed that a subject becomes immune T V = 21 d after the first dose with an average probability of ν = 0.75. Correspondingly, variable V k in (1) denotes the number of individuals who received the first dose of vaccine on day k − T V . In our model, an individual is said to be immune to the disease if he/she will not be infected within the modeled time horizon.
With this simplification, the people in the R and U compartments do not transmit the disease any more, as well as those who are still in the hospital. It is worth remarking that a positive IgG test does not necessarily ensure immunity in this sense. Serological tests suggest that even a relatively high IgG level does not exclude the possibility of reinfection [59].
On the other hand, we assumed that the willingness of the susceptible and recovered patients to vaccinate is roughly the same. Namely, the proportion of susceptible and recovered people vaccinated coincides with the proportion of all susceptibles and recovereds on each day. Therefore, the model counts with ν S k V k /(S k + R k ) susceptible and ν R k V k /(S k + R k ) recovered individuals who achieved immunity at time k. It is worth remarking that the value of V k is known, but cannot be manipulated as the computations were performed on past data of the epidemic spread. Correspondingly, V : k → V k can be considered as a preliminarily known input, or a scheduling variable [60], or a measured disturbance [36,61]. The official European vaccination data including Hungary are available at [62]. Official Hungarian COVID-related data with additional analyses are also available at [63].

Computing the Reproduction Number
To give meaningful estimates from an epidemiological perspective, we computed the basic and the time-dependent effective reproduction numbers. The basic reproduction number, namely the average number of new infections generated by a single infected individual in a fully susceptible population, can be given by the following closed-form expression [21]: It must be noted that using R 0 = 2.2, an early estimate for Hungary commonly used in the literature, and expressing the transmission rate β from the equation, a nominal value of β = 1/3 can be given. However, as this parameter is highly influenced by the circumstances varying in time (e.g., the stringency of restrictions, new variants of the virus, etc.), it is considered to be a time-dependent parameter and estimated as such (the corresponding R 0 being calculated afterwards).
The time-dependent effective reproduction number R t shows the average number of infections caused by a single individual, given the state of the model at time t (thus, taking into account the time-varying nature of beta and the decrease of the susceptible population). This is calculated as follows, It can be a base for comparison between different epidemic-handling strategies, incorporating the strictness of the restrictions as well. In accordance with the traditional notation "R t ", we used t ∈ I T 0 (instead of k) as the time parameter of the time-dependent reproduction number.
Inspired by [21], we considered the daily transmission rate β k of the disease an unknown time-dependent parameter. The past values of β k were computed such that the evolution of the epidemic spread matched the available observations.

Available Measurements
It is commonly stated in the literature [24,25] that the daily number of infected people is not well observable, as the measurement relies on aggressive and exhaustive contact tracing and testing strategies [64,65]. Though it is reasonable to assume that testing is wide-spread and quick enough in the hospitals, the registered numbers of hospitalized patients [66] are still influenced by practical considerations. The limited healthcare capacity on weekends and holidays usually results in a lower number of performed and documented tests, as well as in a delayed hospital discharge. Therefore, following the common engineering practice, we applied a 7 d-long moving average filter to the published data (H Off,raw k ) [66] to avoid biased estimates caused by these administrative inaccuracies. The smoothed time series H Off k is formally calculated as: Obviously, the 7 d-long sliding window must be truncated at both ends of the series. Finally, the filtered hospitalization data were considered as the single available processed measurement, which reveal relevant information about the time-evolution of the process.

State-Space Model Representation and Problem Statement
In Section 2.1, we presented the dynamic Equation (1) of the epidemic spread. We introduced a possible vaccination model in Section 2.2, where V k acts as a measured disturbance input of the dynamical model. In Section 2.3, we explained why parameter β k can be considered as an unknown input of the system. Finally, in Section 2.4, we proposed to consider the hospitalization data (H) as the model output. These "ingredients" allowed us to embed the epidemic spread model into the following discrete-time state-space representation: where is a measured disturbance, and y k = H k ∈ R s is the output with C = (0 0 0 0 0 1 0 0 0). In [22], we presented two possible linear time-invariant (LTI) methods to reconstruct the past states x k+1 and the unknown inputs u k using the measured output y k , k = I T−1 0 . Both techniques in [22] rely on the dynamic inversion of the LTI subsystem of Model (1). In this paper, we revisited the unknown input filtering problem and reformulated it as an optimal predictive tracking control problem. Namely, we computed an optimal input sequence u k , k = I T−1 0 such that the output y k of the system follows the reference signal r k = H Off k , which contains the past output measurements, i.e., the daily number of hospitalized patients (H Off k ). The simultaneous unknown input and state reconstruction can be formulated as the following optimization task.
Problem 1 (NMPC for epidemiological data reconstruction with fixed model parameters). Given the dynamical model (6) with initial condition x 0 , a vector of constants θ, a measured disturbance v k , and a reference output trajectory r k+1 to track (k ∈ I T−1 0 ), we looked for a sequence of inputs u k and states x k+1 that solve the state recursion (6), satisfy the constraint u k ∈ U , and minimize the following weighted cost function: where X = x 1 . . . x T and U = u 0 . . . u T−1 collect the free decision variables of the optimization, Q, R are positive definite weight matrices, and U is a closed subset of the input space R m .
In Problem 1, we formulated a data assimilation problem in the form of a nonlinear model predictive controller (NMPC) computation. The available numerical optimization tools [40][41][42][43] make it possible to solve Problem 1 precisely in a reasonable time. From an epidemiological point of view, the first term of the cost function (7) minimizes the deviation of the computed number of hospitalized patients from the official data, whereas the second term minimizes the slope of the transmission rate of the pathogen. In this way, the NMPC design provides an optimal smooth solution for the unknown transmission rate function β : {0, . . . , T − 1} → U = [0.06, 1], which does not have sudden changes.

Remark 1.
The daily transmission rate of the disease β k is an unknown time-dependent (but, supposedly not abruptly varying) parameter. During an outbreak of the epidemic, the number of infected people is not negligible, namely the sum P k + I k + δA k is significant. In this case, the transmission rate function β : k → β k influences the overall dynamics significantly and determines the shape of the epidemic wave. Therefore, parameter β k is generally well identifiable from the measurements during an outbreak of the epidemic, and it becomes uncertain when the number of infectious people is small. The unknown input filtering task becomes complicated when the model parameters and the initial state (from which the state reconstruction was performed) are probabilistic variables. In the next section, we address the stochastic extension of Problem 1.

Gaussian Assumptions
In this section, we allowed the model parameters to vary in time (θ : k →θ k ), but we assumed that the parameter processθ is a collection of independent identically distributed (i.i.d.) Gaussian random variables: where µ θ is the expected value ofθ k corresponding to the values presented in Section 2.1 and the diagonal Σ θ is its variance. The expected value of the parameter vector contains the nominal values from Table 1, whereas the variances are determined such that the uncertainty intervals from Table 1 resemble the 2σ confidence intervals. Moreover, we assumed that the initial state, from which the prediction was performed, is itself a random variable, namely:x Consequently, every further state and output are random variables, which obey the following stochastic recursion and output equation: Due to the nonlinear terms in the state transition function f , the distribution of the predicted statesx k becomes more and more complicated as we look forward in time (k = 1, 2, . . . ). Therefore, it is very inefficient to compute or at least approximate the non-Gaussian probability density functions of the predicted states for the nonlinear stochastic model (10). As is commonly done in the literature (see, e.g., [50][51][52][53]), we performed a tube-like trajectory estimation. With this technique, each predicted statex k is described by the first two moments, the expected value µ x k , and the variance Σ x k , namely the states are approximated by normal distributions:

Closed-Loop Control Policy
In the literature [44], the values of the optimal control input u are often searched as functions of the states as follows: where µ u k are free decision variables and K is (not necessarily a closed-form) function of the expected state. Thus, the control input is inherently a random variable and is normally distributed as the state (11) itself is approximated by a Gaussian. If Σ xθ k = (Σ θx k ) denotes the covariance betweenx k andθ k , the joint distribution ofx k ,û k , andθ k is: Remark 2 (Nonlinear state-dependent input policy). When K = 0, the optimal tracking problem is said to be an open-loop MPC problem [67], whereas K = 0 results in a so-called closedloop MPC problem [51], where the optimal input policy is parameterized by the state. A stabilizing state feedback (12) is typically useful when the predicted states are random variables, and their actual realizations may deviate from the predicted expectations. When the prediction model is stochastic, a sequence of deterministic input values (K = 0) may result in a diverging sequence of state variances, and hence in a conservative (overly cautious) prediction. When the input is parameterized by the state (K = 0), the adaptability of the input may reduce the uncertainty of the predicted states significantly if the feedback function (12) is determined appropriately. In this sense, the gain function quantifies the trade-off between the uncertainty of the state and the input. Unfortunately, it is not straightforward to compute a stabilizing gain function K for the nominal model (6). Later, in Section 4.5, we demonstrate that a reference state trajectory (if available) makes it possible to compute the values of K separately in each operating reference state through a classical LTI state feedback approach, e.g., a pole placement or a linear quadratic regulator (LQR) design ([68] Section 6.4.2).

Probabilistic Cost and Input Constraint
Problem 1 with the stochastic state Equation (10) and the joint distribution (13) results in a stochastic optimal control problem, where both the cost function (7) and the input constraint are probabilistic. Therefore, the inputs and the states are meant to be found such that they minimize the expected cost, namely: The expanded Formula (14) of the expectation of cost (7) is derived in Appendix A. (14) is typically a nonquadratic function of the mean and the variance of the joint distribution (13). This term introduces a potential difficulty to the optimization, which is addressed later in Section 4.4.
The conditions on the input can be formulated as chance constraints of the form Pr(û k ∈ U ) ≥ p u , where p u denotes the probability level of the confidence set U . Whenû k comprise a single input and the input domain U is an interval, the chance constraint Pr(û k ∈ [u, u]) ≥ p u is equivalent to the following deterministic interval constraint [51]: where Φ : R → (0, 1) denotes the (cumulative) distribution function of the standard normal distribution N (0, 1). This technique for the reformulation of a probabilistic condition is referred to as constraint tightening [69].

Linear Approximation of the State Dynamics around the Expectation
In the literature, there exist different stochastic sample-based optimization approaches for a predictive optimal controller design; see, e.g., [45][46][47][48][49][50]. However, these approaches are computationally tractable only for a shorter prediction horizon. Alternatively, we have the possibility to formulate deterministic recursions for the first two moments of the state vector, e.g., [51] proposed the state transition function f to be approximated by its firstorder Taylor polynomial around the expected values µ k = (µ x k , µ u k , µ θ ) of the probabilistic variables (x k ,û k ,θ k ), namely: This approach leads to a deterministic mean-variance ("µΣ") dynamics, which is typically nonlinear in the free variables µ x k and µ u k .
Note that the linear Taylor approximation ofx k+1 allowed us to express the nonquadratic term in the cost function (14) as follows: The dynamic equations in (18) constitute a possible deterministic prediction model for System (10) and result in the following nonlinear optimal predictive control problem.

Problem 2 (µΣ-NMPC for unknown-input state reconstruction).
Given the dynamical model (18) with an initial state distribution (9), an i.i.d. parameter process (8), a measured disturbance v k , an input policy (12) with a fixed gain function K, and a reference output trajectory r k+1 to track (k ∈ I T−1 0 ), we looked for a sequence of deterministic values µ u k , state moments µ x k+1 , Σ x k+1 , and covariance matrices Σ xθ k+1 with Σ xθ 0 = 0, which solve (18), satisfy the input constraint (16), and minimize the cost (14). The free variables of the optimization are collected in (15). Problem 2 is a stochastic data assimilation problem, reformulated as an optimal predictive tracking problem with a deterministic nonlinear µΣ-prediction model (10). Henceforth, we refer to Problem 2 as a Gaussian or mean-variance NMPC problem (µΣ-NMPC). In general, the variance dynamics (18b,c) significantly increase the complexity of the control problem. If n and p denote the dimension of the statex k and the parameterθ k , respectively, the equations in (18) comprise np + n(n + 1)/2 separate scalar equations, whereas the deterministic model (6) constitutes a system of n scalar equations. Therefore, the µΣ-NMPC in Problem 2 is typically (at least) an order of magnitude more demanding than the ordinary NMPC in Problem 1. However, an appropriate initial guess for the solution of µΣ-NMPC may reduce the computational complexity of the optimization substantially by providing a fast convergence of the solution.

Initial Solution for the µΣ-NMPC Problem
In this section, we compute a pseudo-optimal (i.e., feasible, but not necessarily optimal) solution of Problem 2, which satisfies the dynamic equations (10) and the input constraint (16), but it does not necessarily minimize the cost (14). The computed solution can be considered an initial value for the µΣ-NMPC problem. The solution relies on three observations. First, observe that the mean equation in (18a) resembles the deterministic state recursion in (6) as the mean dynamics is not affected by the variances nor the state-dependent feedback gain K(µ x k ). Therefore, Problem 2 simplifies to Problem 1 if we neglect the variances (S = 0) and their dynamics (18b,c) from the optimization. Accordingly, a possible guess for the expectation (M, V), which solves the mean Equation (18a), can be given by the optimal solution (X * , U * ) of Problem 1 with initial condition x 0 ← µ x 0 and parameter vector θ ← µ θ .
Secondly, we note that the gain function K depends (by design) on the expected states only. This allows computing an appropriate gain K k at each operating point x * 0 = µ x 0 , x * k , k ∈ I T−1 1 along the computed mean solution. We determined K k through the DT version of the LQR design applied to the controllable modes of the pair (A k , B k ), where: Through a sequence of DT-LQR computations, we selected a static feedback gain matrix K k at each time instant k, which minimizes the quadratic cost: For a DT-LTI state-space model x t+1 = A k x t + B k u t (with t = 0, 1, 2, . . . , but a fixed k), the constant gain matrix K k can be computed through simple linear algebra operations ([68] Section 6.4.2), which were implemented in function dlqr of the Control System Toolbox [70] for MATLAB. When selecting the weight matrices Q LQR k and R LQR k of the LQR problem at time k ∈ I T−1 0 , we needed to take into consideration that the value of K k quantifies the trade-off between the uncertainty ofx k+1 andû k . If the locally stabilizing gain has a higher value, the inputû k is more adaptive (hence, more uncertain), but the uncertainty ofx k+1 is smaller. However, the chance constraint (16) does not allow the uncertainty ofû k to increase beyond any bounds. Therefore, the gain should be selected carefully, such that it generates an input distribution satisfying (16). If the first computed value for K k does not result in an admissible input distribution, we are allowed to compute K k multiple times with a gradually decreasing value for R LQR k . Finally, with the knowledge of µ k , v k , K(µ x k ) = K k , k ∈ I T−1 0 , Σ x 0 , and Σ xθ 0 = 0, we computed the variances according to (18b,c), which give the value of S in (15). If the expected values and K are fixed, the variances are well-defined by the variance Equation (18b,c). By construction, the tuple (M, S, V) is a feasible solution for Problem 2 as it satisfies both the µΣ-Equation (18) and the input constraint (16). The computed solution is a good initial guess for the optimal solution of Problem 2. In Algorithm 1, we summarize the proposed operations with a single input u k ∈ R and a simple LQR weight selection.
Algorithm 1 Computing a pseudo-optimal solution for Problem 2. i ← 1.

5:
repeat 6: Compute K k for the pair (A k , B k ) given in (20) through a DT-LQR design with weight matrices Q LQR ← I n and R LQR k ← 2 i−1 ; i ← i + 1.

Results and Discussion
In this section, we present the numerical results we obtained through the MPCbased reconstruction of the unknown epidemiological data. The results were computed in the MATLAB environment with the Control System Toolbox [70]. For algorithmic differentiation, we used CasADi [40,41]. To solve nonlinear MPC problems, we used IPOPT [42], an interior point line search algorithm, with the MUltifrontal Massively Parallel sparse direct Solver (MUMPS) [71,72]. The MPC implementations are available online in the public repository [73].
To compute the unknown epidemiological data, we followed the operations of Algorithm 1 to find a pseudo-optimal solution for the µΣ-NMPC in Problem 2. On 1 March 2020 (k = 0), we assumed a susceptible and almost healthy population, with a small uncertainty as follows:  1, 1, 1, 1, 1, 1, 1, 0).
We note that the effect of the initial number of infected people on the reconstructed state vanishes after a transient period due to the stability properties of the compartmental model (1). Furthermore, we considered pairwise independent random parameters with µ θ and a diagonal Σ θ as presented in Table 1. We fixed the initial value for β to µ u 0 = 1/3 [21]. We solved the ordinary NMPC in Problem 1 to find the candidate mean functions forx andû. Then, we computed the gain matrices and the variances of the joint distribution (13). Using the obtained feasible solution as an initial guess, we solved the µΣ-NMPC in Problem 2. We learned that the local optimum (M * , S * , V * ) found for Problem 2 is qualitatively the same as the initial guess (M, S, V). The relative difference in the cost obtained by the optimal and the pseudo-optimal solution is negligible, namely: If X, X * ∈ R (n+m+np+n(n+1)/2)·T−1 denote the vectors of independent variables of (M, S, V) and (M * , S * , V * ), respectively, the relative difference between X and X * is: In Table 2, we present the dimensional differences between the ordinary NMPC and the µΣ-NMPC if the length of time window is T = 487 d. In Figure 2, we illustrate the computed marginal distributions of the transmission rate of the pathogen and the daily numbers of people in the different stages of the disease. The expectation for the states are presented in Plot 1 of Figure 2, which were computed through the ordinary NMPC design in Problem 1. In each of Plots 2-12 of Figures 2 and 3, the time evolution of the marginal distributions of scalar quantities are visualized, such that the shaded dark and light areas highlight the 1σ and 1σ confidence intervals, respectively. The shape of the epidemic curves clearly show the three waves of the epidemic until summer 2021. Table 2. Computational complexity of Problems 1 and 2 illustrated through the epidemiological data assimilation case study. In this comparison, the cumulative number of recovered people R (all) as an additional state variable is not considered. Accordingly, the number of state variables is n = 9, the number of uncertain parameters is p = 10, and the number of inputs and the measured disturbances are m = 1 and q = 1, respectively. The processing time was measured on a laptop PC with Intel Core i7-4710MQ CPU at 2.50 GHz and 16 GB of RAM. As was noted in Remark 1, the daily transmission rate becomes uncertain when the number of infectious people reduces significantly. In this case, the input constraint with the computed gain K k may be violated; therefore, we increased the input weight R LQR k to obtain a lower gain K k . These heuristic operations to compute an admissible gain were relevant only when the third wave of the epidemic suddenly dropped after the end of April 2021. The scaled logarithm of the input weights on each day k are illustrated in Plot 2 of Figure 2. In Plot 3 of Figure 2, we present the reconstructed number of hospitalized patients in comparison with the official (i.e., reference) data. Plots 4, 5, and 6 of Figure 2 illustrate three derived probabilistic quantitiesẑ k = h(x k ,û k ,θ), namely the time-dependent effective reproduction number (4) ([ẑ k ] 1 = R k ), the number of all infected people ([ẑ k ] 2 =L k +P k +Î k +Â k +Ĥ k ), and the number of daily new infections ([ẑ k ] 3 =β k (P k +Î k +δÂ k )Ŝ k N ). The first and third coordinates ofẑ k are nonlinear functions of random variables. Therefore, the mean and the variance ofẑ k were approximated by the first-order Taylor polynomial of function h, namely,

Quantitative Properties of the Optimization
The yellow curve in Plot 4 illustrates the estimated reproduction number published online by Atlo Team in [74]. In Plot 12 of Figure 3, we present three uncertain time series, namely the number of recovered, but not yet vaccinated people (blue), the cumulative number of recovered people (red), and the cumulative number of immune people (green).
In Figures 2 and 3, we can observe the successfully suppressed first wave with a dramatic effect of the strict lockdown introduced in March 2020. As a result, the disease was mainly confined to closed institutions such as certain hospital wards and elderly homes. This policy could not be maintained in the autumn of 2020, and therefore, a substantially larger second wave occurred, causing a huge burden on the healthcare system. Therefore, further restrictions (online education in secondary schools, closure of certain public spaces, banning of most gatherings, and curfew from 8 p.m. to 5 a.m.) had to be introduced in the first half of November 2020. These measures had the planned effect in terms of the significant reduction of the transmission rate, as is visible in Plot 2 of Figure 3. Then, from January 2021, R t began to increase again due to the appearance of the more contagious alpha (B.  [75]. We note that we can compare these data since they were estimated under the same restriction level. Plot 12 in Figure 3 shows the estimated number of people gaining immunity by infection and/or vaccination. According to this estimation, more than 30% of the population might have gone through the COVID infection until the end of June 2021. This suggests an approximately 26% detection rate. This is significantly lower than the value of certain European countries such as Germany, Italy, or Spain, reaching or sometimes exceeding 50%, but the number of performed tests per population has also been much lower in Hungary than in the mentioned countries [76].

Conclusions
In this paper, we proposed an optimization-based data assimilation approach to compute the unknown inputs and states of discrete-time compartmental epidemic models with uncertain normally distributed parameters. We started from the assumption that the joint state input parameter distribution is Gaussian. Then, a deterministic meanvariance recursion was developed, which made it possible to formulate the stochastic data assimilation problem as a single model predictive control design with a nonlinear mean-variance prediction model. We noted and demonstrated that the resulting µΣ-NMPC is computationally intensive, but its local optimum can be well approximated by a more efficient ordinary NMPC and further closed-form variance computations.
We proposed simple heuristics to predict appropriate feedback gains, which realize state-dependent control actions along the prediction horizon. In this way, a trade-off can be made between the uncertainty of the computed states and inputs as the predicted control action is scheduled by the deviation between the actual realization of the state and its predicted expectation. As the approach does not make a difference between the unknown parameters and inputs, the joint state observation, change detection, and parameter estimation are also possible.
Through the finite horizon predictive control computation, we estimated the unknown data of the past evolution of the COVID-19 epidemic spread within a fixed time window in Hungary. Among the unknown quantities, we considered the daily number of people in the different phases of the disease and the transmission rate of the pathogen, which highly depends on the actual social distancing rules, mobility restrictions, and virus mutations. The unknown time series were computed such that the expected value of the computed number of hospitalized patients fit the truly observed data as much as possible. Compared to our previous results [21,22], we considered an augmented and uncertain compartmental epidemiological model with normally distributed random model parameters and a simple vaccination model as well.
The main limitations of this study were the following. The length of hospital treatment was considered to be constant in the model, since no data have been published on the daily new hospital admissions with COVID-19, from which this parameter could be tracked efficiently. Moreover, no representative nationwide serological testing in Hungary has been organized since the summer of 2020. Such a result would definitely be helpful in making the estimate on the number of immune people more precise. Finally, the waning of immunity after infection or vaccination was also not taken into consideration in the model. However, such an extension does not affect the applicability of the proposed MPC-based estimation.
With a few modifications, the approach can be applied to compute multiple uncertain possibly time-dependent parameters, but also for the prediction of the future behavior of the epidemic spread. The proposed methodology is able to extract and reconstruct detailed information from the whole time horizon of the epidemic process beyond giving estimates for the cumulative number of infected and recovered people. Future work will be focused on the analysis of other European countries.
Funding: This work was supported by the grants "OTKA" SNN125739 and K131545 and the Thematic Excellence Program (TKP2020-NKA-11) of the National Research, Development, and Innovation (NRDI) Office. We also acknowledge the partial support of the Ministry for Innovation and Technology and the NRDI Office within the framework of the Autonomous Systems National Laboratory and the Epidemic Dynamics, Invasion Ecology, and Data-driven Health National Laboratory Programs.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The MPC implementation together with the collected hospitalization [66] and vaccination data [63] are all available online in the public repository [73]. Acknowledgments: P.P. is grateful to Roland Tóth and Tamás Péni for the useful discussions on the theory and computational aspects of stochastic and/or nonlinear MPC problems. We are also very thankful to Krisztina Latkóczy for the inspiring conversations about the microbial background of the ongoing COVID-19 pandemic.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: SEIR Susceptible-exposed-infected-recovered (compartmental model) IgG Immunoglobulin  (borrowed from (14)) The proof is given in multiple steps, but it is essentially based on a simple observation, which allows expressing the expectation of the squared weighted norm of a random variablê x as follows: To prove (A1), first, consider the following chain of identities: Then, we take the trace of the quantities in (A2) to obtain: Equality (A1) is a direct consequence of (A3). Accordingly, the squared weighted norm of the output error at time k + 1 can be expressed as follows: Secondly, the input variation cost is expressed as follows: The variance term in (A5) is further developed as follows: