# Reconstruction of Epidemiological Data in Hungary Using Stochastic Model Predictive Control
## 4 Statistical Analysis for Normally Distributed Model Parameters
### Initial Solution for the \(\mu\Sigma\)-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 \(\mu\Sigma\)-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(\mu_{k}^{x})\). 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}\gets\mu_{0}^{x}\) and parameter vector \(\theta\leftarrow\mu^{\theta}\). 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}^{*}=\mu_{0}^{x}\), \(x_{k}^{*}\), \(k\in\mathbb{T}_{1}^{T-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:  
\[A_{k}=\tfrac{\partial f}{\partial x}(\mu_{k},v_{k}),\ B_{k}=\tfrac{\partial f} {\partial u}(\mu_{k},v_{k}). \tag{20}\]  
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:  
\[\sum_{t=k}^{\infty}\big{(}x_{t}^{\top}\,\mathbf{Q}_{k}^{\text{LQR}}\,x_{t}+u_ {t}^{\top}\,\mathbf{R}_{k}^{\text{LQR}}\,u_{t}\big{)},\ \text{with}\ u_{t}=-K_{k}\,x_{t}. \tag{21}\]  
For a DT-LTI state-space model \(x_{t+1}=A_{k}x_{t}+B_{k}u_{t}\) (with \(t=0,1,2,\ldots\), 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 d1 qr of the Control System Toolbox [70] for MATLAB. When selecting the weight matrices \(\mathbf{Q}_{k}^{\text{LQR}}\) and \(\mathbf{R}_{k}^{\text{LQR}}\) of the LQR problem at time \(k\in\mathbb{T}_{0}^{T-1}\), we needed to take into consideration that the value of \(K_{k}\) quantifies the trade-off between the uncertainty of \(\hat{x}_{k+1}\) and \(\hat{u}_{k}\). If the locally stabilizing gain has a higher value, the input \(\hat{u}_{k}\) is more adaptive (hence, more uncertain), but the uncertainty of \(\hat{x}_{k+1}\) is smaller. However, the chance constraint (16) does not allow the uncertainty of \(\hat{u}_{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 \(\mathbf{R}_{k}^{\text{LQR}}\). Finally, with the knowledge of \(\mu_{k}\), \(v_{k}\), \(K(\mu_{k}^{x})=K_{k}\), \(k\in\mathbb{I}_{0}^{T-1}\), \(\Sigma_{0}^{x}\), and \(\Sigma_{0}^{x\theta}=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 \(\mu\Sigma\)-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}\in\mathbb{R}\) and a simple LQR weight selection. ```
1: Fix \(\delta_{0}\sim\mathcal{N}(\mu_{0}^{x},\Sigma_{0}^{x})\), \(\Sigma_{0}^{x\theta}\gets 0\), \(\mu^{\theta}\), and \(\Sigma^{\theta}\). (Optionally, fix \(\mu_{0}^{u}\).)