Stochastic Model Predictive Control Utilizing Bayesian Neural Networks
Integrating measurements and historical data can enhance control systems
through learning-based techniques, but ensuring performance and safety is
challenging. Robust model predictive control strategies, like stochastic model
predictive control, can address this by accounting for uncertainty. Gaussian
processes are often used but have limitations with larger models and data sets.
We explore Bayesian neural networks for stochastic learning-assisted control,
comparing their performance to Gaussian processes on a wastewater treatment
plant model. Results show Bayesian neural networks achieve similar performance,
highlighting their potential as an alternative for control designs,
particularly when handling extensive data sets.
J. Pohlodek, H. Alsmeier, B. Morabito, C. Schlauch, A. Savchenko, R. Findeisen
# Stochastic Model Predictive Control Utilizing Bayesian Neural Networks

###### Abstract

Integrating measurements and historical data can enhance control systems through learning-based techniques, but ensuring performance and safety is challenging. Robust model predictive control strategies, like stochastic model predictive control, can address this by accounting for uncertainty. Gaussian processes are often used but have limitations with larger models and data sets. We explore Bayesian neural networks for stochastic learning-assisted control, comparing their performance to Gaussian processes on a wastewater treatment plant model. Results show Bayesian neural networks achieve similar performance, highlighting their potential as an alternative for control designs, particularly when handling extensive data sets.

## I Introduction

Optimal, flexible, safe, and reliable control close to the feasibility boundary is becoming increasingly important to achieve many processes' economic, energy-efficient, and sustainable operation. Optimization-based control, such as model predictive control (MPC) [1, 2], is, in principle, well suited to achieve these objectives, as it allows the formulation of elaborate control goals while incorporating state and input constraints. However, the MPC performance heavily depends on the quality of the underlying process model [1, 3], as it is employed to predict the system's behavior and determine the optimal control action. This brings forward the question of counteracting inevitable model-plant mismatch and measurement uncertainties in real systems. Though nominal MPC exhibits inherent robustness properties [1, 4] through repeated optimization in the closed loop, uncertainties degrade its performance and can potentially lead to constraint violations or stability loss [1].

Learning-supported MPC approaches aim to reduce the uncertainty in model dynamics by augmenting it with data-driven parts. However, process uncertainty and stochasticity are often not considered in the MPC formulation if combined with learning approaches; one assumes that the learned part adjusts and compensates for these uncertainties.

In contrast, robust and stochastic MPC approaches [1, 5, 6] are explicitly designed for such a scenario, incorporating the information on model uncertainties directly into the MPC formulation. These formulations typically trade off some performance for guaranteeing constraint satisfaction under assumptions on the uncertainty, which is especially important in safety-critical scenarios.

Combining the two strategies -- robust model predictive control strategies and learning -- has been a subject of recent work [7, 8, 9, 10, 11]. The learned model is often expected to provide some measure of the uncertainty, which is then delivered to the robust or stochastic MPC. Gaussian processes (GPs) have been widely employed in this role, as they explicitly yield a posterior variance along with the regression mean [12]. Though GPs are generally easy to tune, their computational complexity grows cubic with the dataset size, restricting their application to relatively small sets [12]

An alternative is the use of Bayesian neural networks (BNNs) designed to provide a measure of uncertainty besides the regression [13]. Belonging to the family of deep learning models, it is computationally efficient, especially as most of the computational effort is spent offline during training. However, compared to other deep learning methods, BNNs are not yet as widely researched -- it is generally challenging to infer the distributions analytically. Instead, one relies on posterior approximations [13].

In this work, we aim to determine the suitability of BNNs for (stochastic) predictive control applications, with a particular focus on comparing them to state-of-the-art GPs. BNNs have been used in combination with model predictive control, e.g., [14] employs BNNs in combination with a hierarchical MPC approach for the control of a surgical robot to reduce the uncertainty. Here, the kinematics and dynamics of a highly nonlinear robotic system were modeled with the help of BNNs. The uncertainty information provided by the BNN is used in a hierarchical MPC scheme, achieving superior performance. In [15], a learning-based adaptive-scenario-tree model predictive control (MPC) approach is used to achieve probabilistic safety guarantees using BNNs to learn the model uncertainty.

We produced all results in this work using the open-source Python toolbox HILO-MPC1[16]. In a simple-to-use way, the toolbox allows combining (robust) predictive, and optimization-based control and estimation approaches with methods from machine learning, such as Gaussian processes and Bayesian neural networks. 2

Footnote 1: [https://www.ccps.tu-darmstadt.de/research_ccps/hilo_mpc/](https://www.ccps.tu-darmstadt.de/research_ccps/hilo_mpc/)

Footnote 2: The case study used as an example, including all code, will be made freely available in HILO-MPC.

## II Problem Formulation and Methods

We consider dynamic systems which can be described by

\[x_{k+1}=f(x_{k},u_{k})+B\left(d\left(x_{k},u_{k}\right)+w_{k}\right). \tag{1}\]

Here, \(x\in\mathbb{R}^{n_{x}}\) are the dynamical states, \(u\in\mathbb{R}^{n_{u}}\) are the inputs to the system and \(B\in\mathbb{R}^{n_{x}\times n_{d}}\) is a known
matrix. The function \(f\colon\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}\) describes the known part of the system dynamics, while the function \(d\colon\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}\) describes an unknown or difficult to model effect. This unknown effect will be learned using specific machine learning algorithms. The variable \(w\in\mathbb{R}^{n_{d}}\) is assumed to be zero-mean normally distributed process noise \(w\sim\mathcal{N}\left(0,\Sigma^{w}\right)\) with the variance matrix \(\Sigma^{w}=\mathrm{diag}\left(\left[\sigma_{1}^{2}\quad\ldots\quad\sigma_{n_{d} }^{2}\right]\right)\). We focus on Gaussian processes and Bayesian neural networks to learn the unknown effect, which we briefly introduce in the following.

Gaussian ProcessWe briefly review the main concepts of Gaussian processes. For a more detailed introduction to GPs, we refer to [12]. A GP is a stochastic supervised machine learning algorithm used in regression and classification tasks. GPs are less prone to overfitting and naturally provide uncertainty measures on predictions. We formulate a regression task as a mapping \(\psi\colon\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}:\chi\rightarrow\varphi \left(\chi\right)+\nu\) with the input vector \(\chi\in\mathbb{R}^{n_{\chi}}\), the output \(\psi\in\mathbb{R}\) and the zero-mean normally distributed noise \(\nu\sim\mathcal{N}\left(0,\sigma^{2}\right)\) affecting the output. The unknown function \(\varphi\) is assumed to be normally distributed \(\varphi\left(\chi\right)\sim\mathcal{N}\left(m\left(\chi\right),k\left(\chi, \chi\right)\right)\), with the mean function \(m\colon\mathbb{R}^{n_{\chi}}\rightarrow\mathbb{R}\) and covariance function \(k\colon\mathbb{R}^{n_{\chi}}\times\mathbb{R}^{n_{\chi}}\rightarrow\mathbb{R}\). These functions are assumed to depend on hyperparameters, that are determined by GP training procedure, resulting in scalar-valued predictions. For simplicity, multiple outputs are handled as separate one-dimensional GPs independently.

Bayesian Neural NetworkAs a second learning method, we consider Bayesian neural networks [17], which can also provide a measure of uncertainty similar to GPs. In contrast to GPs, computational complexity of BNNs does not grow with the number of data points. Unlike traditional feed-forward NNs, BNNs learn distributions over the output instead of single value predictions. To accomplish this task, the weights of conventional neural networks are replaced by distributions over the weight in each layer. Thus, BNNs extend the class of NNs to estimate posterior distributions instead of most likely values, cf. Fig. 1.

A common way of implementing the weights \(\mathcal{W}\) as distributions is to assume that they are independently Gaussian distributed with zero mean and some arbitrary variance \(\lambda\)[13] as follows

\[p(\mathcal{W})=\prod_{l=i}^{L}\prod_{i=0}^{V_{l}}\prod_{j=1}^{V_{l-1}+1}\mathcal{ N}(w_{i,j}|0,\lambda^{-1}). \tag{2}\]

Here \(L\) denotes here the number of layers in the network and \(V_{l}\) the nodes in each layer \(l\). The likelihoods with regard to the network weights are given by

\[p(y|\mathcal{W},\mathcal{Z})=\prod_{i=1}^{M}\mathcal{N}(y_{i}|d(z_{i}|\mathcal{ W}),\gamma^{-1}), \tag{3}\]

where \(\gamma\) is the precision and \(\mathcal{Z}\) is the collection of nodes \(z_{i}\). Now we can estimate the posterior distribution over the weights given the data via Bayes' rule

\[p(\mathcal{W}|\mathcal{D})=\frac{p(y|\mathcal{W},\mathcal{Z})p(\mathcal{W})}{p( y|\mathcal{Z})}, \tag{4}\]

where we can calculate the marginal likelihood via

\[p(y|\mathcal{Z})=\int p(y|\mathcal{W},\mathcal{Z})\mathcal{W}d\mathcal{W}. \tag{5}\]

In reality, the marginal likelihood and equation (5) cannot be calculated, it becomes intractable because of the nonlinearities in the network caused by the activation \(d(\cdot,\mathcal{W})\). This results in the need to approximate the posterior distribution. This is a common challenge in the Bayesian inference domain, and different methods exist to approximate the needed posterior. The most popular are Markov chain Monte Carlo (MCMC), Laplace approximation (LA), probabilistic backpropagation (PBP), and variational inference (VI), for a survey see [18] and [19].

Nominal Model Predictive ControlAs a comparison and baseline, we consider a nominal model predictive controller, which solves a finite-horizon optimal control problem at the sampling times [1, 2] based on the nominal model3:

Footnote 3: In principle the nominal MPC could also use the learned “nominal” model part, which would be the mean of the predicted state for case of the GP model. This is avoided here, due to space limitations.

\[\min_{\begin{subarray}{c}\{u_{k}\}\end{subarray}} J(\{x_{k}\},\{u_{k}\}) \tag{6a}\] \[\mathrm{s.t.} x_{k+1}=f\left(x_{k},u_{k}\right),\quad x_{0}=\tilde{x}_{j},\] (6b) \[x_{k}\in\mathcal{X},\quad u_{k}\in\mathcal{U}, \tag{6c}\]

where \(J(\{x_{k}\},\{u_{k}\})=\sum_{k=0}^{N}L\left(x_{k},u_{k}\right)+E\left(x_{k}\right)\) is the cost function, \(N\) is the control horizon, \(L\colon\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}\) is the stage cost, \(E\colon\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}\) is the terminal cost, \(\mathcal{X}\) and \(\mathcal{U}\) are the feasible input and state sets. The first input \(u_{k}^{*}\) of the obtained optimal input sequence is applied. The optimization problem is repeated at every sampling time updating \(x_{k}\) with the measured state \(\tilde{x}_{j}\).

Stochastic Model Predictive ControlTo exploit the stochastic uncertainty information of the model, we use stochastic MPC. Compared to the nominal case, we consider chance constraints, allowing for a violation of the constraints

Fig. 1: a) Bayesian neural network with two features (inputs), one hidden layer, and one label (output). b) Single neuron.
with a certain probability [5]. Assuming a cost that contains the expected value of the now stochastic state and control variables minimizing a sequence of optimal policies \(\Pi(x)=\{\pi_{0}(x),...,\pi_{N}(x)\}\) instead of the optimal open-loop inputs, one obtains the following stochastic optimal control problem

\[\min_{\Pi\left(x\right)} \mathbb{E}\left[J(\{x_{k}\},\{u_{k}\})\right] \tag{7a}\] \[\mathrm{s.t.} x_{k+1}=f_{\text{prob}}\left(x_{k},u_{k}\right),\quad x_{0}=\tilde{x }_{j},\] (7b) \[p\left(x_{k}\in\mathcal{X}\right)\geq p_{x},\quad\forall k\in 0,\, \ldots,\,N,\] (7c) \[p\left(u_{k}\in\mathcal{U}\right)\geq p_{u},\quad\forall k\in 0, \,\ldots,\,N,\] (7d) \[u_{k}=\pi\left(x_{k}\right). \tag{7e}\]

with \(f_{\text{prob}}\) defining the probabilistic model of the system, \(p_{x}\) and \(p_{u}\) are the probabilities with which the chance constraints are allowed to be violated and \(k\) describes the time steps of the discrete model.

Since the stochastic optimal control problem is generally intractable, we reformulate and approximate it to make it tractable and include the learning-based hybrid model. The problem's intractability results from the fact that we need to consider infinitely many state trajectories if we solve for an optimal sequence of policies \(\Pi\).

The reformulations and approximations are based on the approximation presented in [9] and the theory of reachable sets presented in [8]. We omit the necessary assumptions for the sake of brevity and state the derived outcome: the system state, control input, and uncertainty dynamics, that is to be learned by a GP or a BNN, are jointly Gaussian distributed:

\[\mathcal{N}(\mu_{k},\Sigma_{k})=\mathcal{N}\left(\begin{bmatrix}\mu_{k}^{x}\\ \mu_{k}^{u}\end{bmatrix},\begin{bmatrix}\Sigma_{k}^{x}&\Sigma_{k}^{xu}&\Sigma_{ k}^{xd}\\ \Sigma_{k}^{x}&\Sigma_{k}^{x}&\Sigma_{k}^{xd}\\ \Sigma_{k}^{dx}&\Sigma_{k}^{du}&\Sigma_{k}^{d}+\Sigma_{k}^{u}\end{bmatrix} \right). \tag{8}\]

Furthermore, we assume that the chance constraints are tightening around the mean of our input and state variables \(\mu^{x}\) and \(\mu^{u}\), which are defined by the assumed Gaussian distributions. With these assumptions, we can formulate of tractable stochastic optimal control problem of the following structure (as we only want to show the basic structure and due to limited space, we omit introducing all symbols)

\[\min_{\mu_{k}^{u}} \mathbb{E}\left[J(\{x_{k}\},\{u_{k}\})\right]\] (9a) \[\mathrm{s.t.} \mu_{k+1}^{x}=\hat{f}\left(x_{k},u_{k},k\right)+B_{d}\mu_{k}^{d},\] (9b) \[\Sigma_{k+1}^{x}\!=\!\Big{[}\nabla\hat{f}\left(\mu_{k}^{x},\mu_{k }^{u}\right)B_{d}\Big{]}\Sigma\!\Big{[}\nabla\hat{f}\left(\mu_{k}^{x},\mu_{k} ^{u}\right)B_{d}\Big{]}^{\!\top}\!,\] (9c) \[\mu_{0}^{x}=\tilde{x}_{j},\mu_{k}^{d}=\mu^{d}\left(\mu_{k}^{x}, \mu_{k}^{u}\right),\Sigma_{0}^{x}=0,\] (9d) \[\Sigma_{k}=\text{according to eq.~{}\eqref{eq:eq
cost \(E\) are defined as \(L\left(x_{k},u_{k}\right)=\left\|x_{k}-x_{\text{ref}}\right\|_{Q_{q}}^{2}+ \left\|u_{k}-u_{\text{ref}}\right\|_{R_{x}}^{2}+\left\|u_{k}-u_{k-1}\right\|_{R _{x}}^{2},E\left(x_{k}\right)=\left\|x_{N}-x_{\text{ref}}\right\|_{Q_{t}}^{2},\) with the reference values \(x_{\text{ref}}=\left[X_{\text{ref}}\text{ }S_{\text{ref}}\right]=\left[1046.28\,\,\,101.615 \right],u_{\text{ref}}=F_{\text{ref}}=0.714286\), and the weights \(Q_{s}=Q_{t}=\left[\begin{smallmatrix}10&0\\ 0&10\end{smallmatrix}\right],\quad R_{s}=1,\quad R_{c}=5\cdot 10^{3}\). The reference of the flow rate \(F_{\text{ref}}\) translates to a hydraulic retention time of \(\tau=7\,\text{d}\), i.e., the average time a volume of wastewater will remain in a particular part of the plant, and lies within the range of hydraulic retention times referenced in [20]. The reference values of the biomass concentration \(X_{\text{ref}}\) and the substrate concentration \(S_{\text{ref}}\) are the steady states of the open loop simulation using the chosen reference flow rate \(F_{\text{ref}}\). The initial conditions of the states were \(X_{0}=0.2\) and \(S_{0}=0\), and both states were constrained to be at least zero over the whole process time. Furthermore, the input was constrained to lie within the range \(0\leq F(t)\leq 2\). The sampling interval was \(\Delta t=0.125\,\text{d}\) and the length of the control horizon was \(N=80\), which translates to a time of \(10\,\text{d}\). The influent substrate concentration was calculated using (III) and was kept constant over the control horizon.

We ran six closed loop simulation using this nominal MPC for a simulation time of \(T_{\text{sim}}=70\,\text{d}\) for each simulation, which ensured that the steady state was reached in each simulation. Additionally, we sampled the reference states for each simulation from a uniform distribution around the actual reference states to gain a sufficient distribution of data in the state space region of interest \(\tilde{X}_{\text{ref}}\sim\mathcal{U}\left(0.9X_{\text{ref}},1.1X_{\text{ref} }\right),\tilde{S}_{\text{ref}}\sim\mathcal{U}\left(0.9S_{\text{ref}},1.1S_{ \text{ref}}\right).\) Both states are assumed to be measured at all times and serve as the features of the machine learning models. The discrepancy between the true plant model and the nominal control model, the \(d\) part of (1)), is used as the label for the machine learning models. Overall, \(2800\) data points were assembled and split between training and test sets with the ratio \(0.8:0.2\). The data was scaled to zero-mean and standard deviation to improve training performance.

Gaussian ProcessSince the training of a GP is computationally expensive, a sparser subset of training data was generated. This was done by iterating through all observations in the training data and disregarding all the observations that have an Euclidean distance to the current observation below a certain threshold. This way we were able to obtain more observations from regions that have been observed very little and disregard observations from well-observed regions in the state space. Since the data was normalized, the threshold was set to \(0.2\). The new training data set obtained using this threshold amounts to \(92\) data points.

As the model of the wastewater treatment plant has two states, we trained two GPs. The noise variance of each GP was fixed to \(\sigma_{X}^{2}=0.1\) and \(\sigma_{S}^{2}=0.01\), respectively. We chose a squared exponential kernel with automatic relevance detection as the covariance function for both GPs \(k\left(\chi_{i},\chi_{j}\right)=\sigma_{f}^{2}\exp\left(-\frac{\left(\chi_{i}- \chi_{j}\right)^{2}}{2\ell^{2}}\right),\) where \(\sigma_{f}\) is the signal variance and \(\ell\) are the length scales. The mean function was assumed to be zero for both GPs, as is typically done [23]. The hyperparameters of the trained GPs can be found in Table II. The length scales of both GPs are relatively close to each other, indicating that every input to the GPs is equally important. Fig. 2 shows the resulting predictive standard deviation for GP\({}_{X}\) in the top row. Its model uncertainty is small in the region around the training data points. In regions where there are no observations available, it predicts a very high uncertainty. Fig. 3 shows the calibration curves of both GPs [24]. The calibration curve for GP\({}_{X}\) is very close to the ideal calibration line, indicating reliable predictions of the model-plant mismatch \(d_{X}\) for the biomass concentration. This can also be shown by calculating the miscalibration area, i.e., the area between the curve and the ideal calibration line. The closer the calibration curve is to the perfect calibration line, the lower the miscalibration area, resulting in higher prediction reliability. The resulting miscalibration area

\begin{table}
\begin{tabular}{|c|c|c|} \hline  & Length scales \(\ell\) & Signal variance \(\ell\) \\ \hline GP\({}_{X}\) & \(2.35425\) & \(2.91528\) & \(5.86936\) \\ GP\({}_{S}\) & \(2.82233\) & \(2.44933\) & \(10.493\) \\ \hline \end{tabular}
\end{table} TABLE II: Hyperparameters of the trained GPs.

Fig. 3: Calibration curves, the ideal calibration is shown red.

Fig. 2: Heat map of the predictive standard deviation for the biomass \(X\). The black dots are the training data points.
is listed in Table III. The calibration curve for the GP trained on the model-plant mismatch in the substrate concentration (GP\({}_{S}\)) is a further away from the ideal calibration line, leading to a higher miscalibration area (see Table III). This indicates a lower prediction reliability of GP\({}_{S}\) compared to GP\({}_{X}\). Two additional metrics, the root-mean-square error (RMSE) as well as the negative log-likelihood (NLL), are listed in Table III. The RMSE reflects the correctness of the predictions, and the NLL is another indicator of the predictions' reliability. Both metrics respectively show that GP\({}_{X}\) has higher accuracy and reliability in the forecasts compared to GP\({}_{S}\), as both values are lower for GP\({}_{X}\).

#### Iv-B1 Bayesian Neural Network

In this work, we used probabilistic backpropagation [25] to approximate the posterior distribution (4) of the BNN, which is a closed-form approximation based on the assumed density filtering approach [26]. Its full derivation is beyond the scope of this work, so we refer the reader to [25] for details. In contrast to our work, [14] and [15] employ variational inference to approximate the posterior distribution.

In probabilistic backpropagation, the prior precision \(\lambda\) and the noise precision \(\gamma\) are assumed to be Gamma distributed hyperpriors \(p\left(i|\alpha_{i},\beta_{i}\right)\sim\Gamma\left(\alpha_{i},\beta_{i}\right)\) where \(\alpha_{i}\) are the shape parameters and \(\beta_{i}\) are the inverse scale parameters. Here, we specified the weight hyperpriors with shape \(\alpha_{\lambda}=6\) and inverse scale \(\beta_{\lambda}=6\) for both outputs. The noise hyperpriors were \(\alpha_{\gamma,X}=\beta_{\gamma,X}=40\) and \(\alpha_{\gamma,S}=\beta_{\gamma,S}=6\), respectively. As with the GPs, the BNNs also need to be trained for each output individually. Both BNNs had one hidden layer with \(50\) nodes and the ReLU activation function [25]. The BNNs were trained for \(10\) epochs each.

Fig. 2 shows the resulting predictive distribution for BNN\({}_{X}\) in the bottom row. Similar to GP\({}_{X}\), it predicts a very low uncertainty in the regions where the training data are located. Unlike GP\({}_{X}\), BNN\({}_{X}\) centers its posterior on the best observed region, shows higher uncertainty in regions with fewer observations, but lower uncertainties in regions without observations. The calibration curves of both BNNs are also very close to the ideal calibration line, indicating high reliability in the predictions. Overall, BNNs and GPs show competitive performance, which is also reflected when comparing the metrics (RMSE, NLL, miscalibration area -- see Table III).

#### Iv-B Open Loop Simulation of Hybrid Models

Given the trained GPs and BNNs, we generated the respective hybrid models. To illustrate their recovery of the model-plant mismatch, we compared open loop step response versus the true plant model and the nominal control model without the process noise. The result for a the step response is shown in Fig. 4. Both hybrid models significantly improve the nominal control model, especially substrate concentration. The performance improvement compared to the nominal control model is substantial, while the mismatch in the steady states of the biomass concentration of a true model is negligible if process noise were considered.

#### Iv-B Closed Loop Control

Finally, learning-supported stochastic MPC for both learning approaches is evaluated. Box constraints of the nominal MPC augmented with a time-varying constraint that ensures survival of the biomass despite disturbed substrate concentration \(S_{f}\). This constraint is realistic for practical applications that assume model-plant mismatches, since drastic drops in \(S_{f}\) could wash out the biomass and drastic increases of \(S_{f}\) could lead to unmodeled behavior, e.g., cell aggregation. Hence, we demand the biomass concentration to be within once the steady state is reached \(X\left(t\right)\in[(X_{\text{ref}}-20)b_{l},(X_{\text{ref}}+20)b_{u}]\). For \(t=30\,\mathrm{d}\) after the steady state is reached, the parameters \(b_{l}\) and \(b_{u}\) are set to \(b_{l}\left(t\right)=\tanh\left(0.1t+0.01\right)\), \(b_{u}\left(t\right)=\tanh\left(0.1t+1\right).\) For the stochastic MPC we reformulated the constraint into a chance constraint with an acceptance probability of \(99\%\). In order to ease the computational burden, the control horizon was shortened to \(N=32\) time steps.

Fig. 5 shows the closed loop behavior of the control designs. The difference between the approaches is mainly

Fig. 4: Open loop step response of the true plant model.

Fig. 5: Closed loop simulations.

\begin{table}
\begin{tabular}{||c|c|c|c||} \hline  & \multicolumn{1}{c||}{RMSE} & \multicolumn{1}{c||}{NLL} & \multicolumn{1}{c||}{Mechanical parameters} \\ \hline GP\({}_{X}\) & \(0.173\) & \(\mathbf{-0.277}\) & \(0.117\) \\ \hline BNN\({}_{X}\) & \(\mathbf{0.165}\) & \(0.478\) & \(\mathbf{0.107}\) \\ \hline GP\({}_{S}\) & \(\mathbf{0.205}\) & \(0.148\) & \(0.135\) \\ \hline BNN\({}_{S}\) & \(0.264\) & \(\mathbf{0.043}\) & \(\mathbf{0.101}\) \\ \hline \end{tabular}
\end{table} TABLE III: Root-mean-square error (RMSE), negative log likelihood (NLL) and miscalibration area of the trained GPs and BNNs. Bold indicates the lowest value for each metric.
present in the vicinity of the constraints. Disturbances in \(s_{f}\) lead to learning-supported MPC violating lower constraint, showing its inability to guarantee constraint satisfaction. In contrast, the learning-supported stochastic MPC designs satisfy the constraints.

Fig. 6 compares the computation time for both stochastic MPC controllers. Multiple GPs were trained and run in the closed loop control setup, varying the number of training points. One can see the increase in total computation time of the corresponding stochastic MPC, and around 130 training data points it surpasses the computation time of the stochastic MPC using the BNN.

## IV Conclusion and Outlook

State-of-the-art learning-supported stochastic MPC requires uncertainty prediction, often using Gaussian processes. However, their computational complexity increases with data set size. We explored Bayesian neural networks (BNNs) for hybrid models in learning-supported stochastic MPC, comparing their performance to Gaussian processes. BNNs achieved similar performance, minimized model-plant mismatch in open-loop simulations, and proved effective in closed-loop simulations. BNNs offer a valuable alternative, efficiently handling large data sets. All results used the open-source Python toolbox HILO-MPC [16].

Future research will investigate BNNs for more complex models and provide strict stability and performance guarantees for specific model classes.

## Acknowledgment

The authors acknowledge funding of the DIGIPOL project (Magdeburg, Saxony-Anhalt) funded in the EU-ERDF scheme and the KI-Embedded project funded by the BMWK.
