Stochastic data-driven model predictive control using Gaussian processes

TL;DR

This paper introduces a Gaussian process-based stochastic model predictive control method that accounts for uncertainty, ensuring constraint satisfaction and improved performance in nonlinear control systems.

Contribution

It proposes a novel GP-based NMPC algorithm with offline Monte Carlo sampling for constraint tightening, enhancing online evaluation speed and handling uncertainty more effectively.

Findings

Successfully applied to a semi-batch bioprocess case study.

Guarantees chance constraint satisfaction through offline Monte Carlo sampling.

Reduces conservativeness by combining online learning with state-dependent uncertainty.

Abstract

Nonlinear model predictive control (NMPC) is one of the few control methods that can handle multivariable nonlinear controlsystems with constraints. Gaussian processes (GPs) present a powerful tool to identify the required plant model and quantifythe residual uncertainty of the plant-model mismatch. It is crucial to consider this uncertainty, since it may lead to worsecontrol performance and constraint violations. In this paper we propose a new method to design a GP-based NMPC algorithmfor finite horizon control problems. The method generates Monte Carlo samples of the GP offline for constraint tighteningusing back-offs. The tightened constraints then guarantee the satisfaction of chance constraints online. Advantages of our proposed approach over existing methods include fast online evaluation, consideration of closed-loop behaviour, and thepossibility to alleviate conservativeness by considering both online learning and state dependency of the uncertainty. The algorithm is verified on a challenging semi-batch bioprocess case study.

Tables

Table 1. Table 1: Parameter values for ordinary differential equation system in Equation 34 .

Parameter	Value	Units
$u_m$	0.0572	${\text{h}}^{-1}$
$u_d$	0.0	${\text{h}}^{-1}$
$K_N$	393.1	mg.L-1
Y${}_{\frac{N}{X}}$	504.5	mg.g-1
$k_m$	0.00016	mg.g-1.h-1
$k_d$	0.281	h-1
$k_s$	178.9	$\mu {\text{mol.m}}^{-2}{\text{.s}}^{-1}$
$k_i$	447.1	$\mu {\text{mol.m}}^{-2}{\text{.s}}^{-1}$
$k_{sq}$	23.51	$\mu {\text{mol.m}}^{-2}{\text{.s}}^{-1}$
$k_{iq}$	800.0	$\mu {\text{mol.m}}^{-2}{\text{.s}}^{-1}$
$K_{Np}$	16.89	mg.L-1

Table 2. Table 2: The mean of the back-off values for the nitrate concentration constraints g 1 / g 3 subscript 𝑔 1 subscript 𝑔 3 g_{1}/g_{3} and the ratio of bioproduct to biomass constraint g 2 subscript 𝑔 2 g_{2} from the final back-off iteration.

Algorithm variation	Mean back-off $g_1/g_3$ (mg/L)	Mean back-off $g_2$ (mg/L)
GP NMPC 50	205.8	0.022
GP NMPC 60	34.2	0.008
GP NMPC 100	38.8	0.003
GP NMPC learning 50	54.0	0.007
GP NMPC 50 SD	23.3	0.002
GP NMPC 50 NSD	36.2	0.004

Table 3. Table 3: Lower bound on the probability of satisfying the joint constraint β ^ l ​ b subscript ^ 𝛽 𝑙 𝑏 \hat{\beta}_{lb} , the approximate satisfaction probability from the final simulation β ^ ^ 𝛽 \hat{\beta} , average computational times to solve a single optimal control problem (OCP) for the GP NMPC, and the average computational time required to complete one back-off iteration.

Algorithm variation	Probability ${\hat{\beta }}_{lb}$	Probability $\hat{\beta }$	OCP time (ms)	Back-off iteration time (s)
GP NMPC 50	0.99	1.00	65	782
GP NMPC 60	0.89	0.91	54	753
GP NMPC 100	0.91	0.93	135	1626
GP NMPC learning 50	0.91	0.93	69	825
GP NMPC 50 SD	0.89	0.91	48	574
GP NMPC 50 NSD	0.91	0.93	49	584

Equations

$${\delta }_{ij}:=\{\begin{array}{ll}1,& \text{if}i=j\\ 0,& \text{otherwise}\end{array}$$

$$\mathrm{diag}(a_0,a_1,\dots ,a_n):=\left[\begin{array}{cccc}{a}_0& 0& \dots & 0\\ 0& a_1& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \dots & 0& a_n\end{array}\right]$$

$$\mathbb{E}[\text{upphi}]:={\int }_{\Omega }\text{upphi}d{p}_{\text{upphi}}$$

$$1\{C\le c\}:=\{\begin{array}{ll}1,& \text{if}C\le c\\ 0,& \text{otherwise}\end{array}$$

$$\mathbb{P}\{\text{upphi}\in \mathcal{A}\}:={\int }_{\text{upphi}\in \mathcal{A}}\text{upphi}d{p}_{\text{upphi}},\mathbb{P}\{\text{upphi}\le c\}:=\mathbb{E}[1\{\text{upphi}\le c\}]$$

$$\text{betainv}(P,A,B)\in F_{\varphi }^{-1}(P|A,B)=\{C:F_{\varphi }(C|A,B)=P)\}$$

$$F_{\varphi }(C|A,B):=\frac{1}{\mathcal{B}(A,B)}{\int }_0^C{t}^{A-1}(1-t{)}^{B-1}dt,\mathcal{B}(A,B)={\int }_0^1{t}^{A-1}(1-t{)}^{B-1}dt$$

$${\mathbf{x}}_{t+1}=\mathbf{f}({\mathbf{x}}_t,{\mathbf{u}}_t)+{\text{upomega}}_t,{\mathbf{x}}_0\sim \mathcal{N}({\text{upmu}}_{{\mathbf{x}}_0},{\mathbf{\Sigma }}_{{\mathbf{x}}_0})$$

$$\mathbf{y}=\mathbf{f}(\mathbf{x},\mathbf{u})+\text{upnu}$$

$$V_T({\mathbf{x}}_0,\mathbf{U})=\mathbb{E}\left[\sum _{t=0}^{T-1}\ell ({\mathbf{x}}_t,{\mathbf{u}}_t)+{\ell }_f({\mathbf{x}}_T)\right]$$

$${\mathbf{u}}_t\in {\mathbb{U}}_t$$

$$\mathbb{P}\left\{\bigcap _{t=0}^T\{{\mathbf{x}}_t\in {\mathbb{X}}_t\}\right\}\ge 1-ϵ$$

$${\mathbb{X}}_t=\{\mathbf{x}\in {\mathbb{R}}^{n_{\mathbf{x}}}\mid g_j^{(t)}(\mathbf{x})\le 0,j=1,\dots ,n_g\}$$

$$\mathbf{Z}=[{\mathbf{z}}^{(1)},\dots ,{\mathbf{z}}^{(N)}{]}^{\mathsf{T}}\in {\mathbb{R}}^{N\times n_{\mathbf{z}}}$$

$$\mathbf{Y}=[{\mathbf{y}}^{(1)},\dots ,{\mathbf{y}}^{(N)}{]}^{\mathsf{T}}\in {\mathbb{R}}^{N\times n_{\mathbf{x}}}$$

$$y=f(\mathbf{z})+\nu$$

$$f(\cdot )\sim \mathcal{G}\mathcal{P}(m(\cdot ),k(\cdot ,\cdot ))$$

$$m(\mathbf{z}):=0$$

$$k(\mathbf{z},{\mathbf{z}}^{\prime }):={\zeta }^2\exp \left(-\frac{1}{2}(\mathbf{z}-{\mathbf{z}}^{\prime }{)}^{\mathsf{T}}{\mathbf{\Lambda }}^{-2}(\mathbf{z}-{\mathbf{z}}^{\prime })\right)$$

$$y\sim \mathcal{G}\mathcal{P}(m(\mathbf{z}),k(\mathbf{z},{\mathbf{z}}^{\prime })+{\sigma }_{\nu }^2{\delta }_{\mathbf{z}{\mathbf{z}}^{\prime }})$$

$$\mathbf{Y}\sim \mathcal{N}(0,{\mathbf{\Sigma }}_{\mathbf{Y}})$$

$$\mathcal{L}(\mathbf{\Psi }):=-\frac{1}{2}\log (\det ({\mathbf{\Sigma }}_{\mathbf{Y}}))-\frac{1}{2}{\mathbf{Y}}^{\mathsf{T}}{\mathbf{\Sigma }}_{\mathbf{Y}}^{-1}\mathbf{Y}$$

$$\left[\begin{array}{c}\mathbf{Y}\\ f(\mathbf{z})\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\mathbf{\Sigma }}_{\mathbf{Y}}& {\mathbf{k}}^{\mathsf{T}}(\mathbf{z})\\ \mathbf{k}(\mathbf{z})& k(\mathbf{z},{\mathbf{z}}^{\prime })\end{array}\right]\right)$$

$$f(\mathbf{z})|\mathcal{D}\sim \mathcal{N}({\mu }_f(\mathbf{z};\mathcal{D}),{\sigma }_f(\mathbf{z};\mathcal{D}))$$

$${\mu }_f(\mathbf{z};\mathcal{D}):={\mathbf{k}}^{\mathsf{T}}(\mathbf{z}){\mathbf{\Sigma }}_{\mathbf{Y}}^{-1}\mathbf{Y}$$

$${\sigma }_f^2(\mathbf{z};\mathcal{D}):={\zeta }^2-{\mathbf{k}}^{\mathsf{T}}(\mathbf{z}){\mathbf{\Sigma }}_{\mathbf{Y}}^{-1}\mathbf{k}(\mathbf{z})$$

$${\mu }_f^{+}(\mathbf{z};{\mathcal{D}}^{+}):={\mathbf{k}}^{+\mathsf{T}}(\mathbf{z}){\mathbf{\Sigma }}_{\mathbf{Y}}^{+-1}{\mathbf{Y}}^{+}$$

$${\sigma }_f^{2+}(\mathbf{z};{\mathcal{D}}^{+}):={\zeta }^2-{\mathbf{k}}^{+\mathsf{T}}(\mathbf{z}){\mathbf{\Sigma }}_{\mathbf{Y}}^{+-1}{\mathbf{k}}^{+}(\mathbf{z})$$

$${\mathbf{k}}^{+}(\mathbf{z})=[{\mathbf{k}}^{\mathsf{T}}(\mathbf{z}),k(\mathbf{z},{\mathbf{z}}^{+}){]}^{\mathsf{T}}$$

$${\mathbf{Y}}^{+}=[{\mathbf{Y}}^{\mathsf{T}},y^{+}{]}^{\mathsf{T}},{\mathbf{Z}}^{+}=[{\mathbf{Z}}^{\mathsf{T}},{\mathbf{z}}^{+}{]}^{\mathsf{T}}$$

$${\mathbf{\Sigma }}_{\mathbf{Y}}^{+-1}={\left[\begin{array}{cc}{\mathbf{\Sigma }}_{\mathbf{Y}}& {\mathbf{k}}^{\mathsf{T}}({\mathbf{z}}^{+})\\ \mathbf{k}({\mathbf{z}}^{+})& k({\mathbf{z}}^{+},{\mathbf{z}}^{+})(+{\sigma }_{\nu }^2)\end{array}\right]}^{-1}$$

$$\mathbf{f}(\mathbf{z})|\mathcal{D}\sim \mathcal{N}({\text{upmu}}_f(\mathbf{z};\mathcal{D}),{\mathbf{\Sigma }}_f(\mathbf{z};\mathcal{D}))$$

$${\text{upmu}}_f(\mathbf{z};\mathcal{D})=[{\mu }_f(\mathbf{z};{\mathcal{D}}_1),\dots ,{\mu }_f(\mathbf{z};{\mathcal{D}}_{n_{\mathbf{x}}}){]}^{\mathsf{T}}$$

$${\mathbf{\Sigma }}_f(\mathbf{z};\mathcal{D})=\mathrm{diag}({\sigma }_f^2(\mathbf{z};{\mathcal{D}}_1),\dots ,{\sigma }_f^2(\mathbf{z};{\mathcal{D}}_{n_{\mathbf{x}}}))+{\mathbf{\Sigma }}_{\text{upomega}}$$

$${\text{upmu}}_f^{+}(\mathbf{z};{\mathcal{D}}^{+})=[{\mu }_f^{+}(\mathbf{z};{\mathcal{D}}_1^{+}),\dots ,{\mu }_f^{+}(\mathbf{z};{\mathcal{D}}_{n_{\mathbf{x}}}^{+}){]}^{\mathsf{T}}$$

$${\mathbf{\Sigma }}_f^{+}(\mathbf{z};{\mathcal{D}}^{+})=\mathrm{diag}({\sigma }_f^{2+}(\mathbf{z};{\mathcal{D}}_1^{+}),\dots ,{\sigma }_f^{2+}(\mathbf{z};{\mathcal{D}}_{n_{\mathbf{x}}}^{+}))+{\mathbf{\Sigma }}_{\text{upomega}}$$

$$\text{mathpzc}{u}_i^{(s)}=\kappa \left({\text{upchi}}_i^{(s)},i\right)$$

Full text

††journal: Computers $\&$ Chemical Engineering

1 Introduction

Model predictive control (MPC) describes an advanced control method that has found a wide range of applications in industry. MPC employs an explicit dynamic model of the plant to determine a finite sequence of control actions to take at each sampling time. The main advantage of MPC is its ability to deal with multivariate plants and process constraints explicitly [39]. Linear MPC is relatively mature and well-established in practice, however many systems display strong nonlinear behaviour motivating the use of nonlinear MPC (NMPC) [3]. NMPC is becoming progressively more popular due to the advancement of improved non-convex optimization algorithms [8], in particular in chemical engineering [9].

The performance of MPC is however greatly influenced by the accuracy of the plant model, which is one of the main reasons why MPC is not more widely used in industry [37]. The development of an accurate plant model has been cited to take up to 80% of the MPC commissioning effort [59]. NMPC algorithms exploit various types of models, commonly developed by first principles or based on process mechanisms [47]. Many mechanistic and empirical models are however often too complex to be used online and in addition have often high development costs. Alternatively, black-box identification models can be exploited instead, such as support vector machines [68], fuzzy models [28], neural networks (NNs) [50], or Gaussian processes (GPs) [32]. For example, recently in Wu et al. [67, 66] recurrent NNs are utilised for an extensive NMPC approach with proofs on closed-loop state boundedness and convergence applied to a chemical reactor. In addition, in Wu et al. [65] the approach was extended updating the recurrent NNs online to further improve the effectiveness.

GPs are an interpolation technique developed by Krige [35] that were popularized by the machine learning community [52]. While GPs have been predominantly used to model static nonlinearities, there are several works that apply GPs to model dynamic systems [22, 30, 15]. GP predictions are given by a Gaussian distribution. The mean of this distribution can be viewed as a deterministic prediction, while the variance can be interpreted as a measure of uncertainty for this deterministic prediction. This uncertainty measure is generally difficult to obtain by nonlinear parametric models [32] and may in part explain the relative popularity of GPs. For control applications this uncertainty measure can be utilized to efficiently learn a dynamic model by exploring unknown regions or avoiding regions with too high uncertainty to improve robustness [7]. So far, GPs have been exploited in a multitude of ways in the control community, including reinforcement learning [21], designing robust linear controllers [61], or adaptive control [17]. In particular, GPs have been shown to be an efficient approach to attain approximate plant models for NMPC.

The use of GPs for NMPC was first proposed in Murray-Smith et al. [44], which updates a GP model online for reference tracking without constraints. In Kocijan et al. [32], Kocijan and Murray-Smith [31] the GP is instead identified offline and utilized online, in which the variance is constrained to prevent the NMPC from steering the dynamic system into regions with high uncertainty. A GP plant model is updated online in Klenske et al. [29] and in Maciejowski and Yang [40] to overcome unmodeled periodic errors or changes to the dynamic system after a fault has occurred respectively. GPs have been shown in Wang et al. [64] to be an efficient means for disturbance forecasting for a linear stochastic MPC approach applied to a drinking water network. GPs have further been applied to approximate the mean and variance required in stochastic NMPC [11]. Grancharova et al. [23] derived an explicit solution for GP-based NMPC. In Cao et al. [16] a GP dynamic model is employed for the control of an unmanned quadrotor, while in Likar and Kocijan [36] a GP dynamic model is exploited to control a gas-liquid separation process. While these and other works show the feasibility of GP-based MPC, there is a lack of efficient approaches to account for the uncertainty measure provided. Model uncertainty can lead to constraint violations and worse performance. To mitigate the effect of uncertainty on MPC, robust MPC [6] and stochastic MPC [42, 24] methods have been developed.

The majority of works for GP-based MPC indeed consider the uncertainty measure provided, however most proposed algorithms employ stochastic uncertainty propagation to achieve this, for example [32, 31, 26, 16, 23, 64]. Hewing et al. [25] give an overview of the various stochastic propagation techniques available. These approaches have some considerable disadvantages, which are:

No known methods to exactly propagate stochastic uncertainties through GP models. Instead, only approximations are available relying on linearization or statistical moment-matching.

• 2.

Increased computational time of GP-based MPC due to the propagation approach itself.

• 3.

Most works consider only open-loop propagation of uncertainties, which is often prohibitively conservative due to open-loop growth of uncertainties.

Recently some papers have proposed different robust GP-based MPC algorithms. In Koller et al. [34] a NMPC algorithm is introduced based on propagating ellipsoidal sets using linearization, that provides closed-loop stability guarantees. This approach may however suffer from increased computational times, since the ellipsoidal sets are propagated online. Furthermore, the method may be relatively conservative due to the use of Lipschitz constants. Maiworm et al. [41] propose the use of a robust MPC approach by bounding the one-step ahead error, while the determination of the required parameters seems to be relatively difficult. Soloperto et al. [56] suggest a robust control approach for linear systems, in which the GP is used to represent unmodeled nonlinearities. The approach is shown to stabilize the linear system despite these uncertainties, which however may have no solution if the difference between the linearized system and the actual nonlinear system is too large.

In this paper we extend an algorithm first introduced in Bradford et al. [13], for which the following extensions were made:

Inclusion of uncertainty for the initial state.

• 2.

Adding additive disturbance noise to the problem definition.

• 3.

Accounting for state dependency on the GP noise.

• 4.

Improved algorithm to obtain the required back-offs using root-finding as opposed to the inverse CDF, which has superior convergence and leads to improved satisfaction of the required probability bounds.

The aim of this approach is to take into account the uncertainty given by a GP state space model for a NMPC finite-horizon control problem. Due to the issues using stochastic uncertainty propagation for the NMPC formulation as highlight previously, we base the NMPC only on cheap evaluations of the GP. This leads to considerably faster evaluation times with little effect on the performance. The proposed method utilizes explicit back-offs, which were recently proposed in [33, 48] to account for stochastic uncertainties in NMPC. These methods generally rely on generating closed-loop Monte Carlo (MC) samples offline from the plant to attain the required back-off values. To obtain exact MC samples of the GP dynamic models we exploit results from Conti et al. [19], Umlauft et al. [60]. There are several important advantages of this new method:

Back-offs are attained using closed-loop simulations, therefore the issue of open-loop growth of uncertainties is avoided.

• 2.

Required computations are carried out offline, such that the online computational times are nearly unaffected.

• 3.

Independence of samples allows some probabilistic guarantees to be given.

• 4.

Explicit consideration of online learning and state dependency of the uncertainty to alleviate conservativeness.

The proposed method is a data-driven NMPC approach relying on black-box identification of a GP model from input/output data pairs. The required data may be obtained either by simulations from a high-fidelity model or by experiments using for example step tests. In the algorithm the state dependency of the uncertainty is accounted for by introducing a penalty term on the variance in the objective. This variance for GPs is a function of the states and hence leads to a trade-off between exploiting the GP model to optimize the objective and avoiding uncertain regions to reduce the spread of the trajectories. The algorithm proposed is aimed at finite horizon control problems, for which batch processes are a particularly important example. They are utilized in many different chemical engineering sectors due to their inherent flexibility to deal with variations in feedstock, product specifications, and market demand. Frequent highly nonlinear behaviour and unsteady-state operation of batch processes have led to the increased acceptance for advanced control solutions, such as NMPC [46]. Works on batch process NMPC accounting for uncertainties include an extended and Unscented Kalman filter based algorithms for uncertainty propagation [45, 10], a NMPC algorithm using min-max successive linearization [62], NMPC algorithms that employ PCEs to account for possible parametric uncertainties [43, 12], and multi-stage NMPC [38].

The paper is comprised of the following sections. In Section 2 the problem definition is given. In Section 3 a general outline of GPs is given including the sampling procedure used. Section 4 shows how the GPs can be exploited to solve the defined problem. In Section 5 the semi-batch bioprocess case study is described, while in Section 6 the results and discussions for this case study are given. Section 7 concludes the paper.

Notation

$\mathbb{N}$ and $\mathbb{R}$ represent the sets of natural numbers and real numbers respectively. The variable $\delta_{ij}$ denotes the Kronecker delta function, such that:

[TABLE]

The notation $\mathop{\mathrm{diag}}(a_{0},a_{1},\ldots,a_{n})$ is used to represent the following diagonal matrix:

[TABLE]

We represent the Gaussian distribution with mean $\bm{\upmu}$ and covariance $\bm{\Sigma}$ as $\mathcal{N}(\bm{\upmu},\bm{\Sigma})$. Further, $\bm{\upphi}\sim\mathcal{N}(\bm{\upmu},\bm{\Sigma})$ denotes that the random variable $\bm{\upphi}$ follows a Gaussian distribution with mean $\bm{\upmu}$ and covariance $\bm{\Sigma}$.

The expected value of a random variable $\bm{\upphi}$ is denoted as:

[TABLE]

where $p_{\bm{\upphi}}$ the probability density function of $\bm{\upphi}$ over the sample space $\Omega$.

Further, we define the indicator function and the probability measure of random variable $\bm{\upphi}$ as follows:

[TABLE]

where $\mathcal{A}$ is a set defining an event on $\bm{\upphi}$ and $\mathbb{P}\{\bm{\upphi}\leq c\}$ is the probability that $\bm{\upphi}$ is less than or equal to $c$.

Lastly, we require the definition of the beta inverse cumulative distribution function (cdf) for a random variable $\phi$. This function $\text{betainv}(P,A,B)$ returns a value $C$ of $\upphi$ following a beta distribution with parameters $A,B$ that has a probability of $P$ to be less than or equal to $C$. The definition is as follows:

[TABLE]

2 Problem definition

The dynamic system in this paper is given by a discrete-time nonlinear equation system with additive disturbance noise:

[TABLE]

where $t$ is the discrete time, $\mathbf{x}\in\mathbb{R}^{n_{\mathbf{x}}}$ is the state, $\mathbf{u}\in\mathbb{R}^{n_{\mathbf{u}}}$ are the control inputs, $\mathbf{f}:\mathbb{R}^{n_{\mathbf{x}}}\times\mathbb{R}^{n_{\mathbf{u}}}\rightarrow\mathbb{R}^{n_{\mathbf{x}}}$ are nonlinear equations, and $\bm{\upomega}$ represents Gaussian distributed additive disturbance noise with zero mean and diagonal covariance matrix $\bm{\Sigma}_{\bm{\upomega}}$. The initial condition $\mathbf{x}_{0}$ is assumed to be Gaussian distributed with mean $\bm{\upmu}_{\mathbf{x}_{0}}$ and covariance matrix $\bm{\Sigma}_{\mathbf{x}_{0}}$.

We assume measurements of the states to be available with additive Gaussian noise expressed as:

[TABLE]

where $\mathbf{y}\in\mathbb{R}^{n_{\mathbf{x}}}$ is the measurement of $\mathbf{f}(\mathbf{x},\mathbf{u})$ perturbed by additive Gaussian noise $\bm{\upnu}\sim\mathcal{N}(\mathbf{0},\bm{\Sigma}_{\bm{\upnu}})$ with zero mean and a diagonal covariance matrix $\bm{\Sigma}_{\bm{\upnu}}=\mathop{\mathrm{diag}}(\sigma^{2}_{\nu_{1}},\ldots,\sigma^{2}_{\nu_{n_{\mathbf{x}}}})$.

The aim of the control problem is to minimize a finite-horizon cost function:

[TABLE]

where $T\in\mathbb{N}$ is the time horizon, $\mathbf{U}=[\mathbf{u}_{0},\ldots,\mathbf{u}_{T-1}]^{\sf T}\in\mathbb{R}^{T\times n_{\mathbf{u}}}$ is a joint vector of $T$ control inputs, $\ell:\mathbb{R}^{n_{\mathbf{x}}}\times\mathbb{R}^{n_{\mathbf{u}}}\rightarrow\mathbb{R}$ is the stage cost, and $\ell_{f}:\mathbb{R}^{n_{\mathbf{x}}}\rightarrow\mathbb{R}$ denotes the terminal cost.

The control problem is subject to hard constraints on the inputs:

[TABLE]

The states are subject to a joint chance constraint that requires the satisfaction of a nonlinear constraint set up to a certain probability, which can be stated as:

[TABLE]

The joint chance constraints are formulated such that the joint event over all $t\in\{0,\ldots,T\}$ of all $\mathbf{x}_{t}$ fulfilling the nonlinear constraint sets $\mathbb{X}_{t}$ has a probability greater than $1-\epsilon$.

For convenience we define the tuple $\mathbf{z}=(\mathbf{x},\mathbf{u})\in\mathbb{R}^{n_{\mathbf{z}}}$ with joint dimension $n_{\mathbf{z}}=n_{\mathbf{x}}+n_{\mathbf{u}}$. The dynamic system in Equation 1 is assumed to be unknown. Instead, we are only given a finite number of noisy measurements according to Equation 2. The available data can then be denoted by the following two matrices:

[TABLE]

where $\mathbf{z}^{(i)}$ represents the input of the $i$-th data point with corresponding noisy observation $\mathbf{y}^{(i)}$, $N$ denotes the overall number of training data points, $\mathbf{Z}$ is a collection of input data, and the corresponding noisy observations are collected in $\mathbf{Y}$.

It should be noted that the uncertainty in this problem arises partially from the uncertain initial condition $\mathbf{x}_{0}$ and the additive disturbance noise $\bm{\upomega}$. Most of the uncertainty however comes from the fact that we do not know $\mathbf{f}(\mathbf{x},\mathbf{u})$ and are only given noisy observations of $\mathbf{f}(\mathbf{x},\mathbf{u})$ instead. To solve this problem we train a GP to approximate $\mathbf{f}(\cdot)$ using the available data in Equation 6. The GP methodology is introduced for this purpose in the next section. This GP then represents a distribution over possible functions $\mathbf{f}(\cdot)$ given the available data, which can be exploited to attain stochastic constraint satisfaction of the closed-loop system.

3 Gaussian processes

3.1 Regression

In this section we introduce the use of GPs to infer a latent function $f:\mathbb{R}^{n_{\mathbf{z}}}\rightarrow\mathbb{R}$ from noisy data. For a more complete overview refer to Rasmussen and Williams [52]. Let the noisy observations $y$ of $f(\cdot)$ be given by:

[TABLE]

where $\mathbf{z}\in\mathbb{R}^{n_{\mathbf{z}}}$ is the argument of $f(\cdot)$ and $y$ is a perturbed observation of $f(\mathbf{z})$ with additive Gaussian noise $\nu\sim\mathcal{N}(0,\sigma^{2}_{\nu})$ with zero mean and variance $\sigma^{2}_{\nu}$.

GPs can be considered a generalization of multivariate Gaussian distributions to describe a distribution over functions. A GP is fully specified by a mean function and a covariance function. The mean function represents the ”average” shape of the function, while the covariance function specifies the covariance between any two function values. We write that a function $f(\cdot)$ is distributed as a GP with mean function $m(\cdot)$ and covariance function $k(\cdot,\cdot)$ as:

[TABLE]

The prior GP distribution is defined by the choice of the mean function and covariance function. In this study we apply a zero mean function and the squared-exponential (SE) covariance function defined as:

[TABLE]

where $\mathbf{z},\mathbf{z}^{\prime}\in\mathbb{R}^{n_{\mathbf{z}}}$ are arbitrary inputs, $\zeta^{2}$ denotes the covariance magnitude, and $\bm{\Lambda}^{-2}\coloneqq\mathop{\mathrm{diag}}(\lambda_{1}^{-2},\ldots,\lambda_{n_{\mathbf{z}}}^{-2})$ is a scaling matrix.

Remark** (Prior assumptions).**

Note the zero mean assumption can be easily achieved by normalizing the data beforehand. The SE covariance function is smooth and stationary, such that choosing the SE covariance function assumes the latent function $f(\cdot)$ to be smooth and stationary as well. The algorithm presented in this work can be utilised using any covariance function. In the case of highly non-stationary functions it may be necessary to use non-stationary covariance functions [54].

From the additive property of Gaussian distributions the measurements of $f(\cdot)$ also follow a GP accounting for measurement noise:

[TABLE]

We denote the hyperparameters defining the prior jointly by $\bm{\Psi}\coloneqq[\zeta,\lambda_{1},\ldots,\lambda_{n_{\mathbf{z}}},\sigma_{\nu}]^{\sf T}$, in which the variance $\sigma_{\nu}$ of the measurement noise is included in case it is unknown. Commonly the hyperparameters are unknown, such that they need to be inferred from the available data using for example maximum likelihood estimation (MLE).

Assume we are given $N$ noisy function evaluations according to Equation 7 denoted by $\mathbf{Y}\coloneqq[y^{(1)},\ldots,y^{(N)}]^{\sf T}\in\mathbb{R}^{N}$ as the result of the inputs given in $\mathbf{Z}=[\mathbf{z}^{(1)},\ldots,\mathbf{z}^{(N)}]^{\sf T}\in\mathbb{R}^{N\times n_{\mathbf{z}}}$. According to the prior GP assumption, the data follows a multivariate Gaussian distribution:

[TABLE]

where $[\bm{\Sigma}_{\mathbf{Y}}]_{ij}=k(\mathbf{z}^{(i)},\mathbf{z}^{(j)})+\sigma^{2}_{\nu}\delta_{ij}$ for each pair $(i,j)\in\{1,\ldots,N\}^{2}$.

The log-likelihood of the observations is consequently given by (ignoring constant terms):

[TABLE]

The MLE estimate of the hyperparameters $\bm{\Psi}$ is determined by maximizing Equation 13. Once the hyperparameters are known, we need to determine the posterior GP distribution of the latent function $f(\cdot)$. From the prior GP assumption we know that the training data and the value of $f(\cdot)$ at an arbitrary input $\mathbf{z}$ follow a joint multivariate normal distribution:

[TABLE]

where $\mathbf{k}(\mathbf{z})\coloneqq[k(\mathbf{z},\mathbf{z}^{(1)}),\ldots,k(\mathbf{z},\mathbf{z}^{(N)})]^{\sf T}$.

The posterior Gaussian distribution of $f(\mathbf{z})$ given the data $(\mathbf{Z},\mathbf{Y})$ can then be found by using the conditional distribution rule for multivariate normal distributions based on the joint normal distribution in Equation 14, which leads to:

[TABLE]

where $\mathcal{D}=(\mathbf{Z},\mathbf{Y})$ denotes the training data available to obtain the posterior Gaussian distribution. The mean function $\mu_{f}(\mathbf{z};\mathcal{D})$ in this context is the prediction of the GP at $\mathbf{z}$, while the variance function $\sigma_{f}^{2}(\mathbf{z};\mathcal{D})$ is a measure of uncertainty.

In Figure 1 we illustrate a prior GP in the top graph and the posterior GP in the bottom graph.

3.2 Recursive update

Often we are given a training dataset $\mathcal{D}$ initially to build a posterior GP and then obtain data points individually afterwards. For example in GP-based MPC we may build an initial GP offline following the procedure shown in Section 3.1, and then also update this model online as new data becomes available. In this work we keep the hyperparameters constant but update the mean function and variance function in Equations 15b-15c recursively given the new data points. Furthermore, the approach used for MC sampling of the GPs to be introduced in Section 3.4 requires recursive updates of this form as well.

Let the new dataset be given by $\mathcal{D}^{+}=(\mathcal{D},(\mathbf{z}^{+},y^{+}))$, where $\mathcal{D}$ is the training data for the initial GP, while $\mathbf{z}^{+}$ is the new input and $y^{+}$ the new corresponding output measurement. Then we will refer to the updated mean function and variance function as:

[TABLE]

The updated terms in Equations 16a-16b can be expressed as:

[TABLE]

where $\mathbf{k}(\mathbf{z})$, $\mathbf{Z}$ and $\mathbf{Y}$ refer to quantities of the initial GP. The noise term $\sigma_{\nu}^{2}$ in the lower diagonal is shown in brackets, since the new ”measurement” $y^{+}$ may be noiseless as is the case for GP MC samples. In this case the noise term should not be added to the new diagonal element.

Note the updates $\mathbf{k}^{+}(\mathbf{z})$, $\mathbf{Z}^{+}$, and $\mathbf{Y}^{+}$ are trivial, however the update of the inverse covariance matrix $\bm{\Sigma}^{+-1}_{\mathbf{Y}}$ is more involved. In essence we require the inverse of the previous covariance matrix after adding a horizontal row and a vertical row to the covariance matrix of the initial GP, see Equation 17c. For this process there are efficient formula available, one of which is introduced in A. These take advantage of the fact that we already know the inverse covariance matrix $\bm{\Sigma}_{\mathbf{Y}}^{-1}$ of the initial GP. Once the update has been carried out, these terms then define the new initial GP. This update procedure is then repeated for the next measurement.

3.3 State space model

In this section we briefly show how the previously introduced GP methodology can be utilized to identify unknown state space models in the form of Equation 1 based on the measurements (data) according to Equation 2. GPs are commonly applied to model scalar functions with vector inputs as shown in Section 3.1. To extend this to the multi-input, multi-output case as required it is common to build a separate independent GP for each output dimension, see for example Deisenroth and Rasmussen [21]. Let the function in Equation 1 be given by $\mathbf{f}(\mathbf{x},\mathbf{u})\coloneqq\mathbf{f}(\mathbf{z})\coloneqq[f_{1}(\mathbf{z}),\ldots,f_{n_{\mathbf{x}}}(\mathbf{z})]^{\sf T}$. We aim to build a separate GP for each function $f_{i}(\mathbf{z})\,\forall i\in\{1,\ldots,n_{\mathbf{x}}\}$ according to Section 3.1. For this purpose we are given observations $\mathbf{Y}_{i}=[y_{i}^{(1)},\ldots,y_{i}^{(N)}]^{\sf T}\,\forall i\in\{1,\ldots,n_{\mathbf{x}}\}$ and corresponding inputs $\mathbf{Z}=[\mathbf{z}^{(1)},\ldots,\mathbf{z}^{(N)}]^{\sf T}$, where $y_{i}$ refers to the $i$-th dimension of measurements obtained according to Equation 2. Let $\mathbf{Y}=[\mathbf{Y}_{1},\ldots,\mathbf{Y}_{n_{\mathbf{x}}}]$ correspond to the overall measurements available. The posterior Gaussian distribution of $\mathbf{f}(\cdot)$ at an arbitrary input $\mathbf{z}=(\mathbf{x},\mathbf{u})$ is:

[TABLE]

where $\mu_{f}(\mathbf{z};\mathcal{D}_{i})$ and $\sigma_{f}^{2}(\mathbf{z};\mathcal{D}_{i})$ are the mean function and variance function built according to Section 3.1 with datasets $\mathcal{D}_{i}=(\mathbf{Z},\mathbf{Y}_{i})\,\forall i\in\{1,\ldots,n_{\mathbf{x}}\}$ with $\mathcal{D}=(\mathbf{Z},\mathbf{Y})$.

Remark** (Additive disturbance noise).**

Note the additive disturbance noise defined in Equation 1 is simply added to the posterior covariance matrix due to the additive property of multivariate Gaussian distributions.

In addition, given an initial GP state space model built with a dataset $\mathcal{D}$ and a new data point $(\mathbf{z}^{+},\mathbf{y}^{+})$, we can update it recursively utilizing the method introduced in Section 3.2:

[TABLE]

where $\mathcal{D}_{i}^{+}=(\mathcal{D}_{i},(\mathbf{z}^{+},y_{i}^{+}))\,\forall i\in\{1,\ldots,n_{\mathbf{x}}\}$ and $\mathcal{D}^{+}=(\mathcal{D},(\mathbf{z}^{+},\mathbf{y}^{+}))$

3.4 Monte Carlo sampling

GPs are distribution over functions and hence their realizations yield deterministic functions, see for example the GP samples shown in Figure 1. In this section we show how to attain independent samples of GP state space models over a finite time horizon. Generating a MC sample of a GP would require sampling an infinite dimensional stochastic process, while there is no known method to achieve this. Instead, approximate approaches have been applied such as spectral sampling [14, 51]. Exact samples of GPs are however possible if the GP MC sample needs to be known at only a finite number of points, which is exactly the situation for state space models over a finite time horizon. This technique was first outlined in Conti et al. [19] and has been employed in Umlauft et al. [60] for the optimal design of linear controllers. We next outline how to obtain an exact sample of a GP state space model over a finite time horizon for an arbitrary feedback control policy.

Assume we are given a GP state space model as shown in Section 3.3 from the input-output dataset $\mathcal{D}=(\mathbf{Z},\mathbf{Y})$. The initial condition $\mathbf{x}_{0}$ is assumed to follow a known Gaussian distribution as defined in Equation 1. A GP state space model represents a distribution over possible plant models, for which each realization will lead to a different state sequence. The aim of this section is therefore to show how to obtain a single independent sample of such a state sequence, which can then be repeated to obtain multiple independent MC samples of the GP. Let $\bm{\mathcal{X}}^{(s)}=[\bm{\upchi}_{0}^{(s)},\bm{\upchi}_{1}^{(s)},\ldots,\bm{\upchi}_{T}^{(s)}]^{\sf T}$ denote such a state sequence, where $s$ denotes a particular GP realization and $\bm{\upchi}_{i}^{(s)}$ the realization of the state at discrete time $i$ in the sequence of MC sample $s$. The control inputs at discrete time $i$ for MC sample $s$ are denoted by $\mathpzc{u}_{i}^{(s)}$. The corresponding joint input of $\bm{\upchi}_{i}^{(s)}$ is represented by $\mathpzc{Z}_{i}^{(s)}=(\bm{\upchi}_{i}^{(s)},\mathpzc{u}_{i}^{(s)})$. We assume the control inputs to be the result of a feedback control policy, which we represent as $\kappa:\mathbb{R}^{n_{\mathbf{x}}}\times\mathbb{R}\rightarrow\mathbb{R}^{n_{\mathbf{u}}}$. The control actions $\mathpzc{u}_{i}^{(s)}$ are then given as follows:

[TABLE]

where $i$ is the current discrete time. Note the control policy depends on the discrete time directly, since it is a finite horizon control policy.

Consequently, the control actions over the finite time horizon $T$ are a function of $\mathcal{X}^{(s)}$ and denoted jointly as $\mathcal{U}^{(s)}=[\mathpzc{u}_{0}^{(s)},\ldots,\bm{\mathpzc{u}}_{T-1}^{(s)}]^{\sf T}=[\kappa\left(\bm{\upchi}_{0}^{(s)},0\right),\ldots,\kappa\left(\bm{\upchi}_{T-1}^{(s)},T-1\right)]^{\sf T}$. Note these control inputs are different for each MC sample $s$ due to feedback.

We start by sampling the Gaussian distribution of the initial state $\mathbf{x}_{0}\sim\mathcal{N}(\bm{\upmu}_{\mathbf{x}_{0}},\bm{\Sigma}_{\mathbf{x}_{0}})$ to obtain the realization $\bm{\upchi}_{0}^{(s)}$. The posterior Gaussian distribution of the next state in the sequence $\mathbf{x}_{1}$ is subsequently given by the GP of $\mathbf{f}(\cdot)$ as defined in Equation 18 dependent on $\bm{\upchi}_{0}^{(s)}$:

[TABLE]

The realization of $\mathbf{x}_{1}$ is obtained by sampling the above normal distribution, which we will denote as $\bm{\upchi}_{1}^{(s)}$. To obtain the next state in the sequence $\bm{\upchi}_{2}^{(s)}$ we need to first condition on $\bm{\upchi}_{1}^{(s)}$, since this is part of this sampled function path. This requires to treat $\bm{\upchi}_{1}^{(s)}$ similarly to a new training point, however without observation noise (i.e. no $\sigma^{2}_{\nu}$ is added to the kernel evaluation $k(\mathpzc{Z}_{0},\mathpzc{Z}_{0})$) and without changing the hyperparameters. Note that if the sampled function were to return to the same input it would lead to the exact same output, since it is conditioned on a noiseless output. This shows that the sampled function is deterministic, since it is the result of sampling. Next we draw $\bm{\upchi}_{2}^{(s)}$ according to the posterior Gaussian distribution obtained from adding the previously sampled data point to the training dataset $\mathcal{D}$ as a noiseless observation. This sample is then again added to the training dataset as a noiseless observation, from which the GP is updated and the next state is drawn. This process is repeated until the required time horizon $T$ has been reached. In this paper we consider a finite time horizon control problem and hence the GP state space model should for, moderate time horizons, not become too computationally expensive, since at most $T$ new data-points are added. Nonetheless, for large time horizons this could become a problem and approximate sampling approaches, such as spectral sampling should then be considered instead [14, 51].

This sampling approach is summarized in Algorithm 1 below and is illustrated in Figure 2. Each GP MC sample is defined by a state and corresponding control action sequence. Note that this gives us a single MC sample and subsequently needs to be repeated multiple times to obtain multiple realizations.

Lastly, we define a nominal trajectory by setting all samples in Algorithm 1 to their mean values. Let $\overline{\bm{\mathcal{X}}}=[\overline{\bm{\upchi}}_{0},\overline{\bm{\upchi}}_{1},\ldots,\overline{\bm{\upchi}}_{T}]^{\sf T}$ refer to this nominal state sequence and $\overline{\bm{\mathcal{U}}}=[\overline{\mathpzc{u}}_{0},\ldots,\overline{\mathpzc{u}}_{T-1}]^{\sf T}$ to the corresponding nominal control sequence, where $\overline{\bm{\upchi}}_{i}$ and $\overline{\mathpzc{u}}_{i}$ are the values of the states and control inputs of the nominal trajectory at discrete time $i$ respectively. Therefore by definition, $\overline{\bm{\upchi}}_{0}=\bm{\upmu}_{\mathbf{x}_{0}}$, $\overline{\mathpzc{u}}_{t-1}=\kappa(\overline{\bm{\upchi}}_{t-1},t-1)$, and $\overline{\bm{\upchi}}_{t}=\bm{\upmu}_{f}(\overline{\mathpzc{Z}}_{t-1};\mathcal{D})$, where $\overline{\mathpzc{Z}}_{t-1}=(\overline{\bm{\upchi}}_{t-1},\overline{\mathpzc{u}}_{t-1})$. Note updating $\bm{\upmu}_{f}(\cdot;\mathcal{D})$ with mean values has no effect.

In Section 4 a NMPC formulation is introduced, which uses the mean function of the GP as the prediction model. The control actions from this NMPC algorithm then define the control policy in Equation 20, while the MC samples represent possible plant responses. This is then exploited to tune the GP NMPC algorithm to attain the desired behaviour.

4 Solution approach

From the input-output dataset $\mathcal{D}=(\mathbf{Z},\mathbf{Y})$ defined in Equation 6, we fit a GP state space model as outlined in Section 3.3. The aim is to solve the problem defined in Section 2 based on this GP state space model. In this context the GP represents a distribution over possible plant models for the process given the available dataset. In Section 3.4 we have shown how to create a sample of this plant model over a finite time horizon $T$, which each lead to different state sequences and corresponding control sequences based on a control policy. In this paper we aim to design a NMPC algorithm based on the GP that acts as this control policy. The MC samples are utilized to tune the NMPC formulation by adjusting so-called back-offs to tighten the constraints and attain the closed-loop probabilistic constraint satisfaction. We next state the GP NMPC formulation, which is based on the tightened constraint set and predictions from the mean function $\bm{\upmu}_{f}(\cdot;\mathcal{D})$ and covariance function $\bm{\Sigma}_{f}(\cdot;\mathcal{D})$.

4.1 Finite-horizon Gaussian process model predictive control formulation

In this section we define the NMPC OCP based on the GP nominal model given the dataset $\mathcal{D}=(\mathbf{Z},\mathbf{Y})$. For the GP NMPC formulation the initial state $\mathbf{x}$ at each sampling time is assumed to be measured or estimated and propagated forward in time exploiting the GP mean function. The predicted states are then used to optimize the objective subject to the tightened constraints. Let the corresponding optimization problem be denoted as $P_{T}\left(\bm{\upmu}_{f}(\cdot;\mathcal{D}),\bm{\Sigma}_{f}(\cdot;\mathcal{D});\mathbf{x},t\right)$ for the current known state $\mathbf{x}$ at discrete time $t$ based on the mean function $\bm{\upmu}_{f}(\cdot;\mathcal{D})$ and covariance function $\bm{\Sigma}_{f}(\cdot;\mathcal{D})$:

[TABLE]

where $\hat{\mathbf{x}}$, $\hat{\mathbf{u}}$, and $\hat{V}_{T}(\cdot)$ refers to the states, control inputs, and control objective of the MPC formulation, $\hat{\mathbf{U}}_{t:T-1}=[\hat{\mathbf{u}}_{t},\ldots,\hat{\mathbf{u}}_{T-1}]^{\sf T}$, $\eta_{k}$ are weighting factors to penalize uncertainty, and $\overline{\mathbb{X}}_{k}$ is a tightened constraint set denoted by: $\overline{\mathbb{X}}_{k}=\{\mathbf{x}\in\mathbb{R}^{n_{x}}\;|\;g_{j}^{(k)}(\mathbf{x})+b_{j}^{(k)}\leq 0,\,j=1,\ldots,n_{g}\}$. The variables $b_{j}^{(k)}$ represent so-called back-offs, which tighten the original constraints $\mathbb{X}_{t}$ defined in Equation 5.

Remark** (Objective in expectation).**

It should be noted that the above objective in Equation LABEL:eq:nominalMPC does not exactly optimize the objective in Equation 3, since it is difficult to obtain the expectation of a nonlinear function. Approximations of this can be found in Hewing et al. [25], however these generally are considerably more expensive and often only lead to marginally improved performance.

Remark** (Scaling for state dependency factors).**

Note we have opted for scalar scaling factors $\eta_{k}$ to account for state dependency. For this to work reliably it is therefore necessary to normalize the data to ensure that all data has approximately the same magnitude, for example normalizing the data to have zero mean and unit variance.

The NMPC algorithm solves $P_{T}\left(\bm{\upmu}_{f}(\cdot;\mathcal{D}),\bm{\Sigma}_{f}(\cdot;\mathcal{D});\mathbf{x}_{t},t\right)$ at each sampling time $t$ given the current state $\mathbf{x}_{t}$ to obtain an optimal control sequence:

[TABLE]

Only the first optimal control action is applied to the plant at time $t$ before the same optimization problem is solved at time $t+1$ with a new state measurement $\mathbf{x}_{t+1}$. This procedure implicitly defines the following feedback control law, which needs to be repeatedly solved for each new measurement $\mathbf{x}_{t}$:

[TABLE]

It is explicitly denoted that the control actions depend on the GP model used. There are several important variations of the GP NMPC control policy. Firstly, it may seem reasonable to update the mean and covariance function using the previous state measurement and corresponding input by applying the recursive update rules introduced in Section 3.2. We will refer to this as learning. This may however lead to a more expensive and less reliable NMPC algorithm.

Secondly, the algorithm may want to avoid regions in which there is great uncertainty due to sparsity of data. This can be achieved by assigning some of the $\eta_{k}$ with non-zero values to penalize the algorithm moving into these regions with high variance. This explicitly takes advantage of the state dependency of the noise covariance function and is hence referred as state dependent. It should be noted that evaluation of the covariance function is computationally expensive. It was determined that setting only $\eta_{t+1}$ to a non-zero value is often sufficient due to the continued feedback update, i.e. penalizing only the variance for the one-step ahead prediction at time $t$. These variations can be summarized as follows:

Learning: Update the mean and variance function of the GP using the previous state measurement and the known corresponding input.

• 2.

State dependent: Set some $\eta_{k}$ not equal to zero, which will lead to the NMPC algorithm trying to find a path that has less variance and hence exploiting the state dependent nature of the uncertainty.

• 3.

Non learning: Keep the mean and variance function the same throughout the run.

• 4.

Non state dependent: Set all $\eta_{k}$ to zero and hence ignoring the possible state dependency of the uncertainty.

Remark** (Full state feedback).**

Note in the control algorithm we have assumed full state feedback, i.e. it is assumed that the full state can be measured without noise. This assumption can be dropped if required by introducing a suitable observer and introduced in the closed-loop simulations to account for this additional uncertainty.

4.2 Probabilistic guarantees

In this section we illustrate how to obtain probabilistic guarantees for the joint chance constraint introduced in Section 2 in Equation 5 based on independent samples of the GP plant model. For convenience we define a single-variate random variable $C(\cdot)$ representing the satisfaction of the joint chance constraint [20]:

[TABLE]

where $\mathbf{X}=[\mathbf{x}_{0},\ldots,\mathbf{x}_{T}]^{\sf T}$ defines a state sequence, and $\mathbb{X}_{t}=\{\mathbf{x}\in\mathbb{R}^{n_{\mathbf{x}}}\mid g_{j}^{(t)}(\mathbf{x})\leq 0,j=1,\ldots,n_{g}\}$.

The probability in Equation 25 is intractable, however a good nonparametric approximation is often achieved utilizing the so-called empirical cumulative distribution function (ecdf). We define the cdf to be approximated as follows:

[TABLE]

Assuming we are given $S$ independent and identically distributed MC samples of $\mathbf{X}$ and hence of $C(\mathbf{X})$, the ecdf estimate of the true cdf in Equation 26 is given by:

[TABLE]

where $\bm{\mathcal{X}}^{(s)}$ is the $s$-th MC sample and $\hat{F}_{C(\mathbf{X})}(c)$ is the ecdf approximation of the true cdf $F_{C(\mathbf{X})}(c)$.

The quality of the approximation in Equation 27 strongly depends on the number of samples used and it is therefore desirable to quantify the residual uncertainty of the sample approximation. This problem has been studied to a great extent in the statistics literature [18]. In addition, there are several works applying these results for chance constrained optimization, see for example Alamo et al. [2]. The main result applied in this study is given below in Theorem 1.

Theorem 1** (Confidence interval for empirical cumulative distribution function).**

Assume we are given a value of the ecdf, $\hat{\beta}=\hat{F}_{C(\mathbf{X})}(c)$, as defined in Equation 27 based on $S$ independent samples of $C(\mathbf{X})$, then the true value of the cdf, $\beta=F_{C(\mathbf{X})}(c)$, as defined in Equation 26 has the following lower and upper confidence bounds:

[TABLE]

Proof.

The proof uses standard results in statistics and can be found in Clopper and Pearson [18], Streif et al. [58]. The proof relies on the following observations. Firstly, $\mathbf{1}\{C(\bm{\mathcal{X}}^{(s)})\leq c\}$ for a fixed value of $c$ describes a Bernoulli random variable, in which either $C(\bm{\mathcal{X}}^{(s)})$ exceeds $c$ and takes the value [math] or otherwise takes the value $1$ with probability $F_{C(\mathbf{X})}(c)$. Secondly, the ecdf describes the number of successes of $S$ realizations of these Bernoulli random variables divided by the total number of samples and hence follows a binomial distribution, i.e. $\hat{F}_{C(\mathbf{X})}(c)\sim\frac{1}{N_{S}}\text{Bin}(S,F_{C(\mathbf{X})}(c))$. The confidence bound for the ecdf can consequently be determined from the Binomial cdf. This method was first introduced by Clopper and Pearson [18] as ”exact confidence intervals”. Due to the close relationship between beta distributions and binomial distributions, a simplified expression can be obtained using beta distributions instead, which leads to the theorem shown [58]. ∎

In other words the probability of $\beta$ exceeding the value $\hat{\beta}_{ub}$ has a probability of $\alpha$ and the probability of $\beta$ being less than or equal to $\hat{\beta}_{lb}$ has also a probability of $\alpha$. In particular, $\hat{\beta}_{lb}$ for small $\alpha$ represents a conservative lower bound on the true probability $\beta$. An illustration of the confidence bound for the ecdf is shown in Figure 3. In general more samples will lead to a tighter confidence bound as expected. The theorem provides a lower bound $\hat{\beta}_{lb}$ that accounts for the statistical error due to the finite sample estimate made, i.e. it gives us a conservative value that is less than or equal to the true probability of feasibility with a confidence level of $1-\alpha$.

We assume we are given $S$ independent samples of the trajectory $\mathbf{X}$ generated according to Section 3.4. Let the approximate ecdf be given by $\hat{\beta}=\hat{F}_{C(\mathbf{X})}(0)$ according to Equation 27, and $\hat{\beta}_{lb}$ the corresponding lower bound according to Theorem 1 with confidence level $1-\alpha$. From this the following Corollary follows:

Corollary 1** (Feasibility probability).**

Assuming the GP representation of the plant model to be a correct description of the uncertainty of the system and given a value of the ecdf $\hat{\beta}=\hat{F}_{C(\mathbf{X})}(0)$, as defined in Equation 27 based on $S$ independent samples and a corresponding lower bound $\hat{\beta}_{lb}\geq 1-\epsilon$ with a confidence level of $1-\alpha$, then the original chance constraint in Equation 5 holds true with a probability of at least $1-\alpha$.

Proof.

The GP MC sample described in Section 3.4 is exact and therefore each sample of the GP plant state space model leads to independent state trajectories $\mathbf{X}$ according to the GP distribution. From $S$ such samples a valid lower bound $\hat{\beta}_{lb}$ to the true cdf value $\beta$ can be determined from Theorem 1 with a confidence level of $1-\alpha$. If $\hat{\beta}_{lb}$ is greater than or equal to $1-\epsilon$, then the following probabilistic bound holds on the true cdf value $\beta$ according to Theorem 1: $\mathbb{P}\left\{\beta\geq\hat{\beta}_{lb}\geq 1-\epsilon\right\}\geq 1-\alpha$, which in other words means that $\beta=\mathbb{P}\left\{C(\mathbf{X})\leq 0\right\}\geq 1-\epsilon$ with a probability of at least $1-\alpha$. ∎

4.3 Determining back-off constraints

In this section we describe how to determine the required back-off values to tighten the constraints of the GP NMPC algorithm in Equation LABEL:eq:nominalMPC. The aim is to choose these values to obtain probabilistic guarantees for the state constraints defined in Equation 5 despite not knowing the exact dynamics. The GP provides a nominal model for the GP NMPC formulation in Equation LABEL:eq:nominalMPC using the mean function and a distribution of possible plant models given the initial dataset in Equation 6. This distribution is exploited to attain different possible plant model realizations to simulate the closed-loop response of the GP NMPC algorithm and then uses the response values to tighten the constraints, such that the original constraint set is satisfied with a high probability according to Equation 5.

We have shown how to obtain the closed-loop trajectory of a control policy according to the realization of a plant model GP distribution in Section 3.4. To this we apply the GP NMPC control policy defined in Equation 24. We propose to utilize $S$ independent MC samples of the GP distribution generated according to Section 3.4, which then in turn describe $S$ different possible plant models with corresponding state and control trajectories. The goal now is to adjust the back-offs, such that the $S$ different state trajectories adhere the original constraint set for all but a few samples to attain the required probability of constraint satisfaction. Let $\bm{\mathcal{X}}^{(s)}=[\bm{\upchi}_{0}^{(s)},\ldots,\bm{\upchi}_{T}^{(s)}]^{\sf T}$ refer to the state trajectory of sample $s$.

The update rule for adjusting the back-offs is based on two steps: First, we define an approximate constraint set, which is then adjusted by a constant factor to obtain the required constraint satisfaction using the ecdf of the joint chance constraint in Equation 30. The approximate constraint set in essence needs to reflect the difference in constraint values for the nominal model of the MPC and the realizations of the GP plant model. We first set all the back-off values to zero and run $S$ MC samples of the GP plant model. As defined in Section 3.4, let $\overline{\bm{\upchi}}_{t}$ refer to the states according to the nominal trajectory with the back-offs set to zero as well. Now assume we aim to obtain back-off values that imply the satisfaction of the following individual chance constraints:

[TABLE]

where $\delta$ is a tuning parameter and should be set to a reasonably low value.

The rule in other words aims to find approximate back-offs for the nominal predictions $\overline{\bm{\upchi}}_{t}$ utilized in the MPC formulation in Equation LABEL:eq:nominalMPC, such that the chance constraints holds for any possible GP plant model MC sample with a probability of $1-\delta$. The parameter $\delta$ in this case is a tuning parameter to obtain the initial back-off values. Note the considered individual chance constraint in Equation 29 is only there to obtain the initial constraint set and is unrelated to the joint chance constraint in Equation 5, which we fulfil by adjusting this approximate constraint set using a root-finding algorithm.

To accomplish this we define the following ecdf based on the $S$ MC samples available:

[TABLE]

where $\hat{F}_{g_{j}^{(t)}}(0)$ is a sample approximation of the chance constraint given in Equation (29) on the RHS.

In Paulson and Mesbah [48] it is proposed to employ the inverse ecdf to approximately fulfill the requirement given in Equation (29) using the $S$ MC samples available. The back-offs can then be stated as follows:

[TABLE]

where $\hat{F}_{g_{j}^{(t)}}^{-1}$ is the inverse of the ecdf defined in Equation 30 and $\tilde{b}_{j}^{(t)}$ refers to these initial back-off values. Note the inverse of an ecdf is given by the quantile function of the discrete $S$ values with probability $\delta$. This gives us the initial back-off values as required for the first step. Note that both the nominal trajectory $\overline{\bm{\upchi}}_{t}$ and the GP realizations depend on the back-off values, which is however ignored since we are only interested in obtaining some reasonable initial values. In the next step these back-off values are further adjusted using a constant back-off factor $\gamma$. The new back-offs are then defined as:

[TABLE]

We aim to change $\gamma$ until the lower bound of the ecdf $\hat{\beta}_{lb}$ as defined in the previous section for the joint chance constraint is equal to $1-\epsilon$ in Equation (5). This is a root finding problem, in which $\gamma$ is adjusted until $\hat{\beta}_{lb}$ reaches the required value:

[TABLE]

where the aim is to determine a value of $\gamma$, such that $h(\gamma)$ is approximately zero.

To attain the required $\gamma$ we use the so-called bisection method [5]. This method determines the root of a function in an interval $a_{\gamma}$ and $b_{\gamma}$, where $h(a_{\gamma})$ and $h(b_{\gamma})$ have opposite signs. In our case this is relatively easy. Setting the value of $\gamma$ too low returns generally a negative value of $h(\gamma)$ due to the constraint violations using low back-offs, while setting it too high leads to positive values leading to a too conservative solution. Note we generally set the initial $a_{\gamma}$ to zero since this corresponds to the $S$ MC samples used to determine $\tilde{b}_{j}^{(t)}$. The bisection method consists of repeatably bisecting the interval, in which the root is contained. The overall algorithm to determine the required back-offs in $n_{b}$ back-off iterations is summarized below as Algorithm 2. The output of the algorithm are the required back-offs with the corresponding lower bound on the probability of satisfying the state chance constraint. Note for $\textit{learning}=true$ the mean and covariance function of the GP are recursively updated utilizing the same procedure as for the update of the GP plant model MC sample.

Remark** (Conservativeness of chance constraint).**

Note to adjust the back-offs we use the ecdf, which does account for the true shape of the underlying probability distribution. This avoids the problem that is often faced in stochastic optimization utilizing Chebyshev’s inequality to robustly approximate chance constraints, which is often excessively conservative [49].

4.4 Algorithm

The overall algorithm proposed in this paper is summarized in this section. Firstly, the problem to be solved needs to be defined as outlined in Section 2. Thereafter, it needs to be decided if the GP NMPC should learn online or exploit the state dependency of the uncertainty as shown in Section 4. Once the GP NMPC has been formulated the back-offs are determined by running closed-loop simulations of the defined problem as shown in Section 4.3. Lastly, these back-offs then give us the tightened constraint set required for the GP NMPC online. This GP NMPC is then run online solving the problem initially outlined. An overall summary can be found in Algorithm 3.

Note the back-offs could be also updated online making use of the new initial conditions and updated prediction model in the case of ”learning”, which would lead to overall less conservativeness [48]. This would however require to carry-out the offline calculations online, which is expensive, and not computationally desirable.

5 Case study

The case study utilized in this paper deals with the photo-production of phycocyanin synthesized by cyanobacterium Arthrospira platensis. Phycocyanin is a high-value bioproduct and its biological function is to enhance the photosynthetic efficiency of cyanobacteria and red algae. It has been considered as a valuable compound because of its applications as a natural colorant to replace other toxic synthetic pigments in both food and cosmetic production. Furthermore, it has shown great promise for the pharmaceutical industry because of its unique antioxidant, neuroprotective, and anti-inflammatory properties. Using a simplified dynamic model we verify our GP NMPC algorithm by operating this process using a limited dataset. The GP NMPC problem is formulated with an economic objective aiming to directly maximize the bioproduct concentration of the final batch subject to two path constraints and one terminal constraint.

5.1 Semi-batch bioreactor model

The simplified dynamic system consists of three ODEs describing the evolution of the concentration of biomass, nitrate, and bioproduct. The dynamic model is based on the Monod kinetics, which describes microorganism growth in nutrient sufficient cultures, where intracellular nutrient concentration is kept constant because of the rapid replenishment. We assume a fixed volume fed-batch. Control inputs are given by the light intensity ($I$) in $\upmu\text{mol.m}^{-2}$.s*-1* and nitrate inflow rate ($F_{N}$) in mg.L*-1*.h*-1* [53]. To capture the effects of light intensity on microalgae growth and bioproduction (photolimitation, photosaturation, and photoinhibition phenomena) the Aiba model is used [1]. The balance equations are given as follows:

[TABLE]

where $C_{X}$ is the biomass concentration in g/L, $C_{N}$ is the nitrate concentration in mg/L, and $C_{q_{c}}$ is the phycocyanin (bioproduct) concentration in the culture in mg/L. The corresponding state vector and control vector are given by $\mathbf{x}=[C_{X},C_{N},C_{q_{c}}]^{\sf T}$ and $\mathbf{u}=[I,F_{N}]^{\sf T}$ respectively. The initial condition is denoted as $\mathbf{x}_{0}=[{C_{X}}_{0},{C_{N}}_{0},{C_{q_{c}}}_{0}]^{\sf T}$. The missing parameter values can be found in Table 1.

5.2 Problem set-up

The time horizon $T$ was set to $12$ with an overall batch time of $240\text{h}$, and consequently the sampling time is $20\text{h}$. Based on the dynamic system in Equation 34 we define the objective and the constraints according to the general problem definition in Section 2. The measurement noise matrix $\bm{\Sigma}_{\bm{\upnu}}$ and disturbance noise matrix $\bm{\Sigma}_{\bm{\upomega}}$ were set to:

[TABLE]

The mean $\bm{\upmu}_{\mathbf{x}_{0}}$ and covariance $\bm{\Sigma}_{\mathbf{x}_{0}}$ of the initial condition are given by:

[TABLE]

The control algorithm aims to maximize the amount of bioproduct produced $C_{q_{c}}$ with a penalty on the change of control actions. The corresponding stage and terminal cost can be stated as follows:

[TABLE]

where ${\bm{\Delta}}_{\mathbf{u}_{t}}=\mathbf{u}_{t}-\mathbf{u}_{t-1}$ and $\mathbf{R}=\mathop{\mathrm{diag}}(3.125\times 10^{-8},3.125\times 10^{-6})$. The overall objective is then defined by Equation 3.

For the case of state dependency we set all $\eta_{i}$ to zero except $\eta_{0}$, see Equation LABEL:eq:nominalMPC. The value of $\eta_{0}$ was set to 15. Note that for these factors to work properly it is important to normalize the data as we did in this case study.

There are two path constraints in the problem. The amount of nitrate is constrained to remain below $800$ mg/L, while the ratio of bioproduct to biomass may not exceed $11.0$ mg/g for high density biomass cultivation. These constraints can be stated as:

[TABLE]

Lastly, there is a terminal constraint on nitrate to reach a final concentration of below $150$ mg/L:

[TABLE]

The maximum probability for violating the joint chance constraint was set to $\epsilon=0.1$. The control inputs light intensity and nitrate inflow rate are constrained as follows:

[TABLE]

For the back-off iterations we employed $S=1000$ MC iterations with the initial back-offs computed according to Equation 31 with $\delta=0.1$ and $\alpha=0.01$. The maximum number of back-off iterations for the bisection algorithm was set to $n_{b}=16$.

5.3 Implementation details and initial dataset generation

The optimization problem for the GP NMPC in Equation LABEL:eq:nominalMPC is solved using Casadi [4] to obtain the gradients of the problem using automatic differentiation in conjunction with IPOPT [63]. The ”real” plant model was simulated using IDAS [27]. In the next section, different variations of the proposed algorithm are presented, for which two different type of datasets were collected. For the first type of dataset we designed the entire input data matrix $\mathbf{Z}$ according to a Sobol sequence [55] in the range $\mathbf{z}_{i}\in[0,20]\times[50,800]\times[0,0.18]\times[120,400]\times[0,40]$. The ranges were chosen for the data to cover the expected operating region. The corresponding outputs $\mathbf{Y}$ were then obtained from the IDAS simulation of the system perturbed by Gaussian noise as defined in the problem setup. In the second approach only the control inputs were set according to the Sobol sequence in the range $\mathbf{u}_{i}\in[120,400]\times[0,40]$ and the corresponding states $\mathbf{Y}$ were obtained from the trajectories of the ”real” system perturbed by noise using samples of the initial condition and the time horizon as defined in the problem setup based on these control inputs. The system was simulated in ”open-loop” using these control actions, i.e. without any feedback controller present. For both datasets the input data $\mathbf{Z}$ and output data $\mathbf{Y}$ are normalized to zero mean and a standard deviation of one. The reason we use two different types of datasets is to highlight the advantages of accounting for state dependency in two of the algorithm variations. In the first dataset the data is relatively evenly distributed and hence considering the state dependency of the uncertainty to avoid regions with high data sparsity has essentially no effect, while in the second approach there are clearly defined trajectories that can be followed by accounting for state dependency.

6 Results and discussions

In this section we present and discuss the results from the case study described in the previous section. For comparison purposes we compare six different variations of the proposed GP NMPC approach, which are as follows:

GP NMPC 50, 60, 100: GP NMPC approach without learning and without taking into account state dependency for dataset sizes of 50, 60, and 100 points using the first type of dataset.

• 2.

GP NMPC learning 50: GP NMPC approach with learning and without state dependency for a dataset size of 50 points, which will be compared to the above case of 50 data points without learning. The first type dataset is utilized.

• 3.

GP NMPC SD/NSD 50: GP NMPC approach with and without accounting for the state dependency for a dataset size of 50 points employing the second type of dataset.

In addition, we compare the approaches to a nominal NMPC algorithm based on the GP model to show the importance of employing back-offs to prevent constraint violations, i.e. we run the GP NMPC on the ”real” plant model, while setting the back-offs to zero. The results of the outlined runs are summarized in Figures 6-12, and in Tables 2-3. In Figures 6-7 we show the evolution of the back-off factor and the probability of constraint satisfaction $\hat{\beta}_{lb}$ over the 16 back-off iterations from Algorithm 2. The next two Figures 4-5 show the 1000 MC trajectories of the constraints with a line to highlight the nominal prediction of the GP NMPC. Next the GP NMPC was applied to the ”real” plant with back-offs from the final iteration and without back-offs referred to as nominal as shown in Figures 8-9. Figure 10 shows the probability density function of the objective values obtained from the ”real” plant, where in the figure larger objective values correspond to better objective values. Figure 11 shows representative control trajectories for GP NMPC 50, 60, 100 compared with the optimal trajectory obtained from solving the OCP of the ”real” plant ignoring uncertainties. Figure 12 shows the back-off values for the nitrate constraints $g_{1}$ and $g_{2}$ for GP NMPC 50 and GP NMPC 50 learning. Lastly, Table 2 shows the mean values for the back-offs averaged over time for the final back-off iteration, while Table 3 shows the attained probability of satisfaction $\hat{\beta}_{lb}$ together with the average computational times for solving a single GP NMPC optimization problem. We can draw the following conclusions from these results:

Figures 6-7 and Table 2 show that apart from GP NMPC 50 the other variations reach the required $\hat{\beta}_{lb}$ and hence successfully converge to a reasonable back-off factor. For these, as expected, a low back-off value leads to too low $\hat{\beta}_{lb}$ values near zero, while too high back-off values lead to too high $\hat{\beta}_{lb}$ values. The value of $\hat{\beta}_{lb}$ does vary by $\pm 0.01$ even on convergence, which is due to the randomness of the MC samples. Nonetheless, since $\hat{\beta}_{lb}$ is a sample robust value, it is high enough if at least $0.9$ is reached once over the 16 iterations, which is the case for all of them. Note that GP NMPC 50 does not converge, since even without back-offs the NMPC remains feasible for all MC trajectories and hence the bisection procedure fails. The performance of GP NMPC 50 is however also by far the worst. This is due to insufficient amounts of data in crucial areas for the control problem.

• 2.

From GP NMPC 50, 60 to 100 the objective values steadily increase and hence improve with increased number of data points as shown in Figure 10. This is as expected, since more data points should lead to a more accurate GP plant model and hence more optimal control actions. Lastly, a more accurate GP plant model should also require less conservative back-offs. For the constraint $g_{2}$ the mean of the back-off values steadily decreases from $0.022$ mg/L for GP NMPC 50 to $0.003$ for GP NMPC 100 as shown in Table 2. For the constraints $g_{1}/g_{3}$ on the other hand the mean of the back-offs decreases dramatically from GP NMPC 50 with a value of $250$mg/L to a value of $34.2$mg/L for GP NMPC 60, while slightly increasing again for GP NMPC 100 to $38.8$ mg/L. This is further illustrated in Figure 4, for which the spread of the trajectories decreases steadily from GP NMPC 50, 60 to 100. All in all, larger datasets lead to improved solutions.

• 3.

The learning approach GP NMPC 50 learning leads to a reasonable solution with an objective value that is on average higher and therefore an improvement over GP NMPC 60 as can be seen in Figure 10. Further, the sharper peak of the objective value suggests a more reliable performance. In contrast, GP NMPC 50 without learning is unable to determine a good solution and therefore has an objective value that is considerably worse than the remaining scenarios with an objective value that is on average over 30$\%$ lower. GP NMPC 50 is also unable to reach a $\hat{\beta}_{lb}$ value of $0.9$ and has instead a much higher value as shown in Table 3, but at the expense of performance. This is highlighted in Figure 11 in which the control trajectory of GP NMPC 50 decreases the nitrate flowrate early to satisfy the terminal constraint, which leads to a sub optimal solution. This is believed to be due to the high uncertainty of the $g_{1}/g_{2}/g_{3}$ constraint trajectories of GP NMPC 50, which can be seen by the large spread of the trajectories in Figure 4. GP NMPC 50 learning on the other hand has a spread of the constraints $g_{1}/g_{2}/g_{3}$ that is significantly less than GP NMPC 50. This is further highlighted by the considerably higher back-off values of GP NMPC 50 compared to GP NMPC 50 learning as can be seen in Table 2, where the $g_{1}/g_{3}$ back-off values are nearly 400$\%$ larger, while the $g_{2}$ back-off values are over $300\%$ larger. In conclusion, accounting for online learning has lead to a significantly better solution, although it should be noted that at larger datasets the effect is nearly negligible.

• 4.

GP NMPC 50 can be seen to be less erratic than GP NMPC 50 learning in Figure 4, which is however expected. GP NMPC 50 uses the same prediction model throughout and hence the control inputs are only influenced by changes in the initial condition. For GP NMPC 50 learning on the other hand the prediction model changes at each sampling point and hence this leads to more irregular behaviour, which can be seen by the increased oscillations. This then leads to overall improved control actions and performance due to exploiting the new data available, as can be seen in Figure 10.

• 5.

It can be seen in Figure 5 that GP NMPC NSD 50 has a larger spread of trajectories than GP NMPC SD 50. This is as expected, since GP NMPC SD 50 aims directly in its objective to minimize uncertainty. Consequently, the mean of the back-off values of constraints $g_{1}/g_{3}$ and $g_{2}$ in Table 2 are over $50\%$ larger and over $100\%$ greater than those for GP NMPC SD 50, respectively. Nonetheless, the attained average objective of GP NMPC SD 50 is marginally lower and hence worse as can be seen in Figure 10. This is expected, since accounting for state dependency reduces conservativeness, but may also lead to a sub optimal solution due to conflicting objectives.

• 6.

In Table 3 the average computational times of a single GP NMPC evaluation are shown, which range from $135$ms to $48$ms. Overall, it can be seen that by far the largest computational times are attributed to GP NMPC 100, which is reasonable since the complexity of GP plant models grows exponentially with the number of data points. GP NMPC 50 and GP NMPC 50 learning can be seen however to have higher computational times, which is due to a more complex optimization problem from the reduced amount of data. This is further highlighted by GP NMPC 50 SD and GP NMPC 50 NSD attaining the lowest average computational times due to the second type of dataset leading to easier optimization solutions. In Table 3 we can further see that the computational times required for a single back-off iteration is for all variations nearly solely determined by the GP NMPC evaluation time. In addition, in Table 3 we show the ecdf value $\hat{\beta}$ for the final back-off iteration. It can be seen that these probabilities are substantially higher than the required probability of $0.9$ ranging from $0.91$ to $0.93$ for the converged solutions, which is due to the conservativeness of the probabilistic lower bound used. This leads to higher back-off values than required and therefore worse objective values. This conservativeness can be reduced by setting $\alpha$ to a higher value.

• 7.

Lastly, in Figures 8-9 trajectories of the constraints are shown by applying the GP NMPC variants to the ”real” plant. It can be seen that the nominal variations of GP NMPC 50, 60, 100, and GP NMPC 50 learning with back-offs set to zero violate the nitrate constraint $g_{1}$ to remain below $800$mg/L by a substantial amount of up to $50$mg/L for GP NMPC 50 learning and GP NMPC 60. With back-offs on the other hand the approaches remain feasible throughout the run, which illustrates the importance of employing back-offs. GP NMPC NSD 50 nominal can be also seen to violate the nitrate constraint $g_{1}$ by $50$mg/L, while GP NMPC SD 50 nominal does not violate this constraint. This is likely due to GP NMPC SD 50 nominal following a feasible trajectory in the dataset. GP NMPC NSD 50 with back-offs remains feasible. Overall, it can be seen that back-offs are important to achieve feasibility given the presence of plant-model mismatch.

7 Conclusions

In conclusion, a new approach is proposed for finite-horizon control problems using NMPC in conjunction with GP state space models. The method utilizes the probabilistic nature of GPs to sample deterministic functions of possible plant models. Tightened constraints using explicit back-offs are then determined, such that the closed-loop simulations of these possible plant models are feasible to a high probability. In addition, it is shown how probabilistic guarantees can be derived based on the number of constraint violations from the simulations. Furthermore, it is shown that online learning and state dependency of the uncertainty can be taken into account explicitly in this method, which leads to overall less conservativeness. Moreover, the computational times are shown to be relatively low, since constraint tightening is performed offline. Finally, through the comprehensive semi-batch bioprocess case study, the efficiency and potential of this method for the optimisation of complex stochastic systems (e.g. biological processes) is well demonstrated.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 675215.