Linearized Gaussian Processes for Fast Data-driven Model Predictive Control

TL;DR

This paper introduces a linearized Gaussian Process approach combined with Sequential Convex Programming to enhance the speed and scalability of Gaussian Process-based Model Predictive Control, making it more suitable for real-time applications.

Contribution

The paper proposes a novel linearized Gaussian Process method with Sequential Convex Programming to significantly improve GP-based MPC's computational efficiency and scalability.

Findings

The new method is faster than existing GP-MPC approaches.

It scales better with larger models and maintains predictable performance.

Demonstrated effectiveness in a numerical example.

Abstract

Data-driven Model Predictive Control (MPC), where the system model is learned from data with machine learning, has recently gained increasing interests in the control community. Gaussian Processes (GP), as a type of statistical models, are particularly attractive due to their modeling flexibility and their ability to provide probabilistic estimates of prediction uncertainty. GP-based MPC has been developed and applied, however the optimization problem is typically non-convex and highly demanding, and scales poorly with model size. This causes unsatisfactory solving performance, even with state-of-the-art solvers, and makes the approach less suitable for real-time control. We develop a method based on a new concept, called linearized Gaussian Process, and Sequential Convex Programming, that can significantly improve the solving performance of GP-based MPC. Our method is not only faster but also much more scalable and predictable than other commonly used methods, as it is much less influenced by the model size. The efficiency and advantages of the algorithm are demonstrated clearly in a numerical example.

Tables

Table 1. TABLE I: Solving time statistics for several GP model sizes.

	$N=500$		$N=700$		$N=900$		$N=1100$		$N=1300$		$N=1500$
	$\mu$	$\sigma$	$\mu$	$\sigma$	$\mu$	$\sigma$	$\mu$	$\sigma$	$\mu$	$\sigma$	$\mu$	$\sigma$
linGP	0.096	0.0011	0.112	0.0015	0.134	0.0013	0.154	0.0020	0.178	0.00291	0.209	0.0042
linGP-inner	0.024	0.0002	0.024	0.0003	0.024	0.0002	0.024	0.0003	0.024	0.0002	0.024	0.0003
Knitro	0.246	0.0081	0.538	0.0172	0.896	0.0208	1.390	0.0175	2.040	0.0365	2.740	0.0390
Ipopt	0.331	0.0385	0.751	0.0909	1.29	0.1700	2.010	0.2580	2.950	0.3680	4.010	0.4680

Equations

$$m(x;\theta )$$

$$k(x,x^{\prime };\theta )$$

$${\mu }_{\star }$$

$${\sigma }_{\star }^2$$

$$f(x+{\Delta }_x)\approx {\tilde{f}}_x({\Delta }_x)=f(x)+{\Delta }_x^T{\nabla }_{x}f(x),{\Delta }_x\in {\mathbb{R}}^n\text{.}$$

$${\tilde{f}}_x({\Delta }_x)={\left[\begin{array}{c}1\\ {\Delta }_x\end{array}\right]}^T\left[\begin{array}{c}f(x)\\ g(x)\end{array}\right]={\hat{x}}^T\hat{f}(x)\text{.}$$

$$\left[\begin{array}{c}Y\\ f\\ g\end{array}\right]\sim N\left(0,\left[\begin{array}{ccc}K(X,X)+{\sigma }_n^2\mathrm{I}& k(X,x)& K^{(0,1)}(X,x)\\ k(x,X)& k(x,x)& K^{(0,1)}(x,x)\\ K^{(1,0)}(x,X)& K^{(1,0)}(x,x)& K^{(1,1)}(x,x)\end{array}\right]\right)\text{.}$$

$$\hat{f}|X,Y,x\sim N\left({\hat{m}}_x,{\hat{V}}_x\right)$$

$${\hat{m}}_x=\left[\begin{array}{c}k(x,X)\\ K^{(1,0)}(x,X)\end{array}\right]{\left(K(X,X)+{\sigma }_n^2\mathrm{I}\right)}^{-1}Y$$

$$\begin{array}{rl}{\hat{V}}_x=& \left[\begin{array}{cc}k(x,x)& K^{(0,1)}(x,x)\\ K^{(1,0)}(x,x)& K^{(1,1)}(x,x)\end{array}\right]-\\ & \left[\begin{array}{c}k(x,X)\\ K^{(1,0)}(x,X)\end{array}\right]{\left(K(X,X)+{\sigma }_n^2\mathrm{I}\right)}^{-1}\left[\begin{array}{cc}k(X,x)& K^{(0,1)}(X,x)\end{array}\right]\end{array}$$

$${\tilde{f}}_x({\Delta }_x)\sim \mathcal{N}\left({\hat{m}}_x^T\hat{x},{\hat{x}}^T{\hat{V}}_x\hat{x}\right)\text{.}$$

$${\tilde{\mu }}_{f|x}({\Delta }_x)={\hat{m}}_x^T\hat{x},{\tilde{\sigma }}_{f|x}^2({\Delta }_x)={\hat{x}}^T{\hat{V}}_x\hat{x}\text{,}$$

$$y_k=f(y_{k-l},\dots ,y_{k-1},u_{k-m},\dots ,u_k)\text{.}$$

$$\begin{gathered} x_k^T=[\overline{y}{}_{k-l},\overline{y}{}_{k-l+1},\dots ,\overline{y}{}_{k-1}, \\ {u}_{k-m}^T,u_{k-m+1}^T,\dots ,u_{k-1}^T] \end{gathered}$$

$$y_k\sim N\left(\overline{y}{}_k,{\sigma }_{y,k}^2\right)$$

$$\overline{y}{}_k={\mu }_f(x_k,u_k),{\sigma }_{y,k}^2={\sigma }_f^2(x_k,u_k)\text{.}$$

$$\underset{{\mathcal{U}}_t,{\mathcal{Z}}_t}{\mathrm{minimize}}J\left(\overline{\mathcal{Y}}{}_t,{\mathcal{S}}_{y,t},{\mathcal{U}}_t,{\mathcal{Z}}_t\right)$$

$$\overline{y}{}_k={\mu }_f(x_k,u_k)$$

$${\sigma }_{y,k}^2={\sigma }_f^2(x_k,u_k)$$

$$z_k\in Z,u_k\in U$$

$$g_i\left(\overline{\mathcal{Y}}{}_t,{\mathcal{S}}_{y,t},{\mathcal{U}}_t,{\mathcal{Z}}_t\right)\le 0,i=1,\dots ,n_{\mathrm{ieq}}$$

$$h_i\left(\overline{\mathcal{Y}}{}_t,{\mathcal{S}}_{y,t},{\mathcal{U}}_t,{\mathcal{Z}}_t\right)=0,i=1,\dots ,n_{\mathrm{eq}}$$

$$\underset{\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t}{\mathrm{minimize}}J\left(\overline{\mathcal{Y}}{}_t,{\mathcal{S}}_{y,t},{\mathcal{U}}_t^{\star }+\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t\right)$$

$$\overline{y}{}_k={\mu }_f(x_k^{\star }+\Delta x_k,u_k^{\star }+\Delta u_k)$$

$${\sigma }_{y,k}^2={\sigma }_f^2(x_k^{\star }+\Delta x_k,u_k^{\star }+\Delta u_k)$$

$$z_k\in Z,u_k^{\star }+\Delta u_k\in U$$

$$g_i\left(\overline{\mathcal{Y}}{}_t,{\mathcal{S}}_{y,t},{\mathcal{U}}_t^{\star }+\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t\right)\le 0,i=1,\dots ,n_{\mathrm{ieq}}$$

$$h_i\left(\overline{\mathcal{Y}}{}_t,{\mathcal{S}}_{y,t},{\mathcal{U}}_t^{\star }+\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t\right)=0,i=1,\dots ,n_{\mathrm{eq}}$$

$$\begin{gathered} \Delta x_k^T=[\Delta {\tilde{y}}_{k-l},\Delta {\tilde{y}}_{k-l+1},\dots ,\Delta {\tilde{y}}_{k-1}, \\ \Delta u_{k-m}^T,\Delta u_{k-m+1}^T,\dots ,\Delta u_{k-1}^T]\text{.} \end{gathered}$$

$$\Vert \Delta u_k\Vert \le \rho ,\Vert \Delta {\tilde{y}}_k\Vert \le \rho$$

$$\underset{\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t}{\mathrm{minimize}}J\left({\tilde{\mathcal{Y}}}_t,{\tilde{\mathcal{S}}}_{y,t},{\mathcal{U}}_t^{\star }+\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t\right)$$

$${\tilde{y}}_k={\hat{m}}_{x_k^{\star },u_k^{\star }}^T{\xi }_k,{\tilde{\sigma }}_{y,k}^2={\xi }_k^T{\hat{V}}_{x_k^{\star },u_k^{\star }}{\xi }_k$$

$${\xi }_k^T=\left[1,\Delta x_k^T,\Delta u_k^T\right]$$

$$\Delta {\tilde{y}}_k={\tilde{y}}_k-\overline{y}{}_k^{\star },\text{Δxk defined in eqrefeq:diff-state-vector}$$

$$\Vert \Delta u_k\Vert \le \rho ,\Vert \Delta {\tilde{y}}_k\Vert \le \rho$$

$$z_k\in Z,u_k^{\star }+\Delta u_k\in U$$

$$g_i\left({\tilde{\mathcal{Y}}}_t,{\tilde{\mathcal{S}}}_{y,t},{\mathcal{U}}_t^{\star }+\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t\right)\le 0,i=1,\dots ,n_{\mathrm{ieq}}$$

$$h_i\left({\tilde{\mathcal{Y}}}_t,{\tilde{\mathcal{S}}}_{y,t},{\mathcal{U}}_t^{\star }+\Delta {\mathcal{U}}_t,{\mathcal{Z}}_t\right)=0,i=1,\dots ,n_{\mathrm{eq}}$$

$$\phi\left(\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t},\mathcal{S}_{y,t},\mathcal{U}_{t},\mathcal{Z}_{t}\right)=J\left(\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t},\mathcal{S}_{y,t},\mathcal{U}_{t},\mathcal{Z}_{t}\right)+\\ \lambda\big{(}\textstyle\sum_{i=1}^{n_{\mathrm{ieq}}}\max\left(0,g_{i}\left(\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t},\mathcal{S}_{y,t},\mathcal{U}_{t},\mathcal{Z}_{t}\right)\right)+\\ \textstyle\sum_{i=1}^{n_{\mathrm{eq}}}\left\lvert h_{i}\left(\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t},\mathcal{S}_{y,t},\mathcal{U}_{t},\mathcal{Z}_{t}\right)\right\rvert\big{)}$$

$$\displaystyle\underset{\Delta\mathcal{U}_{t},\mathcal{Z}_{t}}{\text{minimize}}\quad\phi\big{(}\tilde{\mathcal{Y}}_{t},\tilde{\mathcal{S}}_{y,t},\mathcal{U}^{\star}_{t}+\Delta\mathcal{U}_{t},\mathcal{Z}_{t}\big{)}$$

$$\delta^{(j)}=\phi^{(j)}-\phi\big{(}\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t},\mathcal{S}_{y,t},\tilde{\mathcal{U}}_{t},\tilde{\mathcal{Z}}_{t}\big{)}$$

Full text

I Introduction

Model predictive control (MPC) is an advanced control approach that utilizes a mathematical model of the controlled process to predict its future responses over a finite time horizon, then, based on these predictions, minimizes a cost function to obtain optimized control inputs subject to input and state constraints. This procedure is repeated at every time step, resulting in an optimal control strategy. The MPC approach is attractive due to its ability to effectively and intuitively handle complex dynamics and system constraints [1]. Stochastic MPC (SMPC) can efficiently handle stochastic uncertainties with underlying probability distributions in optimal control problems [2]. For these reasons, MPC and its variants have been applied widely in many practical control problems [3].

A major caveat of MPC is that it requires a reasonably accurate mathematical model of the system because its performance is highly dependent on the accuracy of the model-based predictions. Traditionally, these models are developed using first principles based on physics. However, for complex systems where knowledge of the system dynamics is often incomplete or hard to obtain fully, the required effort for such model development and engineering would pose significant difficulty and limit the use of physics-based models for MPC [4]. To overcome this challenge, data-driven modeling for MPC has recently been studied. This approach uses machine learning techniques to learn a system model from data, possibly with certain prior knowledge about the system, and uses that model in place of the traditional mathematical model in MPC. Such data-driven MPC approach has become increasingly appealing, especially for complex and large-scale systems, as a result of recent advancements in machine learning, optimization, and computation.

Gaussian Processes (GPs) – a type of statistical machine learning models – have been used for modeling dynamical systems and for MPC [5]. GPs are highly flexible and able to capture complex behaviors with fewer parameters than other machine learning techniques, hence they generally work well with small data sets [6]. More importantly, a GP provides an estimate of uncertainty or doubt in the predictions through the predictive variance, which can be used to assess or guarantee the performance of a learning-based system. Because of these advantages, GPs have been employed successfully in MPC frameworks and applications [5, 7, 4, 8]. In this paper, we will consider the GP-based MPC approach for data-driven predictive control, which we will call GP-MPC henceforth.

Computational complexity is a major challenge of GPs, thus of GP-MPC. GP computations generally scale cubically with the size of the training data [6]. When the training data size increases, GP-MPC quickly becomes very expensive to solve. This is particularly true when the GP variance, which represents model uncertainty, is used in the cost function or the constraints, leading to a non-convex and highly demanding optimization problem. As a result, a general-purpose nonlinear program (NLP) solver is usually used to solve GP-MPC, for instance by Interior-Point methods (IPM), Sequential Quadratic Programming (SQP) methods, and active-set methods [5]. However, the performance of these methods in GP-MPC problems is usually unsatisfactory. For instance, it was reported in [7] that a moderate-sized GP-MPC problem took up to a minute to solve with Ipopt. Section V reports the solving time of three general NLP methods for a numerical GP-MPC example, illustrating the computational performance issues of these methods.

Our goal is to make GP-MPC fast, scalable, and more predictable compared to other NLP methods. In particular, we aim to make solving GP-MPC less influenced by the training data size. We propose a new concept called linearized Gaussian Process (linGP), which is a (valid) Gaussian Process of the linearization of the latent function about a given input. This is not the same as linearizing the GP mean and variance because a linGP is still a GP and therefore it has mean, variance, and moments of any order. We then develop a Sequential Convex Programming (SCP) algorithm based on linGP to solve a general class of GP-MPC problems very efficiently and, most importantly, much more independently of the training data size compared to other methods. The proposed SCP algorithm is inspired by the successive convexification approach in [9] to solve non-convex optimal control problems and similar in concept to successive linearization methods to solve nonlinear MPC (such as in [10]). However, our approach is different from other attempts to accelerate GP-MPC and machine-learning-based MPC in general, for example [8, 10], in that we approximate not the output equation of the machine learning model but the underlying latent process of the model.

The concept of linGP is introduced in Section II, followed by the considered class of GP-MPC problems in Section III. We present the linGP-based SCP algorithm in Section IV. Finally, the efficacy of our method is demonstrated in a numerical example in Section V.

II Linearized Gaussian Processes (linGP)

This section introduces the concept of a linearized Gaussian Process (linGP) and develops the mathematical derivation of a linGP from a given Gaussian Process (GP). We will begin with a brief introduction to GPs and the GP-related notations used in this paper.

II-A Gaussian Processes: Basics and Notations

Given $N$ noisy observations $y^{(i)}$ of an underlying function $f:\mathbb{R}^{n}\mapsto\mathbb{R}$ through a Gaussian noise model: $y^{(i)}=f\left(x^{(i)}\right)+\epsilon^{(i)}$, with inputs $x^{(i)}\in\mathbb{R}^{n}$, $\epsilon^{(i)}\sim\operatorname*{\mathcal{N}}\left(0,\sigma_{n}^{2}\right)$, and $i=1,\dots,N$ . A GP of $f$, considered as a distribution over functions [6], is essentially a probability distribution on the observations of all possible realizations of $f$. The prior $\mathbb{P}(f)$ represents the initial belief about this distribution. Let $\mathcal{D}=(X,Y)$ be the finite set of observation data of $f$, where $X=[x^{(1)},\dots,x^{(N)}]\in\mathbb{R}^{n\times N}$ collects the regression vectors and $Y=[y^{(1)},\dots,y^{(N)}]\in\mathbb{R}^{N}$ contains the corresponding observed outputs. By conditioning the GP on $\mathcal{D}$, the posterior $\mathbb{P}(f|\mathcal{D})$ is the updated distribution after seeing the observation data, which allows probabilistic inference at a new input.

The GP assigns a joint Gaussian distribution to any finite subset of random variables $\{f^{(1)},\dots,f^{(M)}\}$ corresponding to inputs $\{x^{(1)},\dots,x^{(M)}\}$. It is fully specified by its mean function $m(x)$ and covariance function $k(x,x^{\prime})$, parameterized by the hyperparameters $\theta$,

[TABLE]

The mean function is often set to zero without loss of generality.

The output $f_{\star}$ of the GP corresponding to an input $x_{\star}$ is a random variable $f_{\star}\sim\operatorname*{\mathcal{N}}\left(\mu_{\star},\sigma_{\star}^{2}\right)$, where

[TABLE]

in which $K_{\star}=[k(x_{\star},x^{(1)}),\dots,k(x_{\star},x^{(N)})]\in\mathbb{R}^{N}$, $K_{\star\star}=k(x_{\star},x_{\star})\in\mathbb{R}$, and $K\in\mathbb{R}^{N\times N}$ is the covariance matrix with elements $K_{ij}=k(x^{(i)},x^{(j)})$. Note that the mean and variance are nonlinear in $x_{\star}$ and their computations scale cubically with the size of $X$ and $Y$.

In theory, the hyperparameters $\theta$ are also random variables, whose posterior distributions are obtained by conditioning them on $\mathcal{D}$ using the Bayes’ theorem. In practice, however, $\theta$ are often obtained by maximizing the likelihood: $\operatorname*{arg\,max}_{\theta}\Pr(Y|X,\theta)$. For an in-depth treatment of GPs, interested readers are referred to [6].

II-B Linearized Gaussian Processes

Consider a GP of a nonlinear function $f(x):\mathbb{R}^{n}\to\mathbb{R}$, with covariance function $k(x,x^{\prime}):\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$. We will denote the GP of $f$ by $\mathcal{G}_{f}$, whose predictive mean and variance at an input $x\in\mathbb{R}^{n}$ are given by (1) and are denoted by $\mu_{f}(x)$ and $\sigma_{f}^{2}(x)$, respectively. For notational simplicity we will take the mean function $m(x)$ to be zero. If the GP has a non-zero mean function, it is equivalent to the sum of $m(x)$ and a zero-mean GP, hence our results can be extended to this case.

The GP inference at any input $x$, given by (1), is non-convex due to the covariance function $k(\cdot,\cdot)$. This makes any optimization involving GP regressions non-convex and generally difficult to solve. However, in many of these problems, non-convexity only comes from the GP regressions involved. One typical approach to overcome this issue is to approximate the mean and variance by linear or quadratic equations using Taylor’s series expansions, which often results in an approximated linear or quadratic program. For instance, this approach is employed by SQP solvers. However, as shown in Section I, this approach often does not work well. In this paper, we take a different approach, described below.

Given an input vector $x$, we assume that $f$ is differentiable at $x$. We also assume that the covariance function $k(\cdot,\cdot)$ is differentiable, which is typically the case [11]. Recall that the GP of $f$ is essentially a (posterior) distribition over the space of realizations of $f$. We can then linearize $f$ around $x$ as

[TABLE]

Define $g\equiv\nabla_{x}f$ as the gradient of $f$. The above linearization can be rewritten as, with $\hat{x}=[1,\Delta_{x}^{T}]^{T}$,

[TABLE]

It is important to note here that $\hat{f}(x)=[f(x),g^{T}(x)]^{T}$ is a random vector of length $(n+1)$ which, through the inner product with the vector $\hat{x}$, results in the random variable $\tilde{f}_{x}(\Delta_{x})$ of the linearized process. In other words, we approximate the original GP of ${f}$ around ${x}$ by a process of linearized function ${\tilde{f}_{x}}$ of ${f}$, which we will call a linearized Gaussian Process or linGP. The linGP of the GP $\mathcal{G}_{f}$ at input $x$ is denoted by $\tilde{\mathcal{G}}_{f|x}$.

As differentiation is a linear operator, $\hat{f}$ is a (multivariate) Gaussian Process derived from the original GP. Therefore, $\hat{f}(x)$ is a Gaussian random vector that defines a posterior distribution, derived from the posterior $\mathbb{P}(f|X,Y)$, over the values of $f(x)$ and its gradiant at $x$. Consequently, the linGP $\tilde{\mathcal{G}}_{f|x}$ has a distribution derived from $\mathbb{P}(f|X,Y)$ over the linear approximations of the function $f$ near $x$ . To characterize this linGP, we need to characterize the GP of $\hat{f}$.

Let $K^{(1,0)}=(\nabla_{x}k)$ – the gradient of $k$ with respect to the first argument – and $K^{(0,1)}=(D_{x^{\prime}}k)$ – the Jacobian of $k$ with respect to the second argument. Let $K^{(1,1)}(x,x^{\prime})$ be the $n\times n$ matrix such that $K^{(1,1)}_{i,j}(x,x^{\prime})=\frac{\partial^{2}}{\partial x_{i}\partial x^{\prime}_{j}}k(x,x^{\prime})$. Following [11], we can write the joint distribution of the observation data, the process $f$, and the gradient process $g$ as

[TABLE]

Here, we use the convention that when a function is applied with $X$, it is broadcast along the corresponding dimensions. For example, $K(X,X)$ is a $N\times N$ matrix with $K_{i,j}(X,X)=k(x^{(i)},x^{(j)})$; similarly, the $i^{\mathrm{th}}$ row of $K^{(0,1)}(X,x)$ is $K^{(0,1)}(x^{(i)},x)$. The posterior distribition of $\hat{f}$ conditioned on the observation data $(X,Y)$ is then

[TABLE]

where $\mathrm{I}$ denotes an identity matrix of appropriate dimensions. The linearized process in (2) can then be written as follows.

[TABLE]

The mean of the original GP at $(x+\Delta_{x})$ is approximated by a linear function of $\hat{x}$, while its variance is approximated by a quadratic function of $\hat{x}$. Because $\hat{f}$ is a valid GP, $\hat{V}_{x}$ is positive semi-definite (PSD) and, therefore, it is guaranteed that the variance of $\tilde{f}_{x}(\Delta_{x})$ is non-negative. We emphasize again that the linGP is in fact a Gaussian Process of ${\tilde{f}_{x}}$, which approximates the original GP around ${x}$. A key advantage of the linGP is that its mean and variance are respectively linear and convex quadratic functions. We will use $\tilde{\mu}_{f|x}(\Delta_{x})$ and $\tilde{\sigma}_{f|x}^{2}(\Delta_{x})$ as the notations for the mean and variance of the linGP $\tilde{\mathcal{G}}_{f|x}$ at input $(x+\Delta_{x})$ in a neighborhood of $x$. They are calculated as

[TABLE]

in which $\hat{m}_{x}$ and $\hat{V}_{x}$ are given in (4), and $\hat{x}$ is defined in (2).

Remark 1

The above derivations can easily be extended to the case when a subset of elements in $x$ are fixed. In that case, the elements of $\hat{x}$, $K$, $K^{(0,1)}$, $K^{(1,0)}$, and $K^{(1,1)}$ corresponding to the fixed elements of $x$ are removed.

Remark 2

The mean and variance of the predictive output of the GP can be approximated directly based on Taylor series, as commonly employed by SQP-based NLP methods. A key distinction between such an approximation approach and linGP is that while the former linearizes the mean and variance of the probabilistic predictive output (i.e., linearization happens after prediction), linGP linearizes the latent function directly and derives the probabilistic distribution of such linearized functions (i.e., linearization happens before prediction). The mean $\tilde{\mu}_{f|x}$ obtained above is exactly the same as if the exact mean equation (1a) is linearized directly. However, the variance $\tilde{\sigma}_{f|x}^{2}$ is not the same as if the exact variance equation (1b) is approximated to the second order by a Taylor series. Indeed, while $\tilde{\sigma}_{f|x}^{2}$ is a true variance and guaranteed to be non-negative, it is not necessarily the case with a second-order Taylor series approximation. As we will explain later in Section IV-C, this key difference is an important advantage of the linGP approach, that could lead to a substantial performance improvement of our algorithm compared to other methods.

III GP-based Model Predictive Control

In this section, we will consider a general MPC formulation using GP models.

III-A Gaussian Processes for Dynamical Systems

GPs can be used for modeling nonlinear dynamical systems [5]. This can be achieved by feeding autoregressive, or time-delayed, input and output signals back to the model as regressors. Specifically, it is common to model a discrete-time dynamical system by an autoregressive nonlinear function

[TABLE]

Here, $k$ denotes the time step, $u\in\mathbb{R}^{n_{u}}$ the control input, $y\in\mathbb{R}$ the output, and $l$ and $m$ are respectively the lags for autoregressive outputs and control inputs. The regressor vector includes not only the current control input $u_{k}$ but also the past (i.e., autoregressive) values of $y$ and $u$. The vector of all autoregressive inputs can be thought of as the current state of the model. The nonlinear function $f$ of the dynamical system can be learned by a GP $\mathcal{G}_{f}$, trained from the system’s data in the same way as any other GPs. MPC formulations that use GP models of the involved dynamical systems will be called GP-MPCs in this paper. We assume that a single GP is used in a GP-MPC, however our results readily extend to the case when multiple GPs are involved.

Due to the predictive nature of a GP-MPC, it is necessary to simulate the GP $\mathcal{G}_{f}$ over $H$ future steps, where $H$ is the MPC horizon, and predict its multistep-ahead behavior. Because the output of $\mathcal{G}_{f}$ is a distribution rather than a point estimate, the autoregressive outputs fed to the model beyond the first step are random variables, resulting in more and more complex output distributions as we go further. Therefore, a multistep simulation of $\mathcal{G}_{f}$ involves the propagation of uncertainty through the model. There exist several methods for approximating uncertainty propagation in GPs [5].

It was shown in [7] that the zero-variance method, which replaces the autoregressive outputs with their corresponding expected values and therefore does not propagate uncertainty, could achieve sufficient prediction accuracy. Its computational simplicity is attractive, especially in GP-MPC where the GP must be simulated for many time steps. In this paper, we will assume the zero-variance method for predicting future GP outputs in MPC. Let us define the state vector $x$ such that, at time step $k$,

[TABLE]

where $\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{i}$ is the expected value of $y_{i}$. Note that all $\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{k}$ for $k<t$, where $t$ is the current time, can be observed.

The output of $\mathcal{G}_{f}$ at time $k\geq t$ then reads

[TABLE]

III-B Model Predictive Control with Gaussian Process Models

In this paper, we consider a general GP-MPC formulation described below. Let $t$ be the current time step and $H>0$ the MPC horizon length. Define the collections $\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t}=\{\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{k}\,|\,k\in\mathcal{I}_{t}\}$, $\mathcal{S}_{y,t}=\{\sigma_{y,k}\,|\,k\in\mathcal{I}_{t}\}$, and $\mathcal{U}_{t}=\{u_{k}\,|\,k\in\mathcal{I}_{t}\}$, where $\mathcal{I}_{t}=\{t,\dots,t+H-1\}$ contains the indices of all time steps in the current MPC horizon. We also introduce additional variables $z\in\mathbb{R}^{n_{z}}$ that are used in the MPC formulation but not in the GP model, for example to model non-GP dynamics in the system. The collection of these variables in the MPC horizon is denoted by $\mathcal{Z}_{t}=\{z_{k}\,|\,k\in\mathcal{I}_{t}\}$.

The GP-MPC reads

[TABLE]

in which $k$ ranges from $t$ to $t+H-1$ and

• •

$J(\cdot)$ is the objective function and typically has the form $J=\textstyle\sum_{k=t}^{t+H-1}c_{k}\left(\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{k},\sigma_{y,k},u_{k},z_{k}\right)$ where $c_{k}(\cdot)$ are the stage cost functions;

• •

Constraints (10b) and (10c) represent the GP dynamics;

• •

Sets $Z$ and $U$ in (10d) are convex constraint sets of the variables $z$ and $u$, respectively;

• •

There are $n_{\mathrm{ieq}}$ inequality constraints of the form (10e);

• •

There are $n_{\mathrm{eq}}$ equality constraints of the form (10f);

We make the following assumptions about the GP-MPC formulation (10).

Assumption 1

Suppose that, for every $k$, $\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{k}$ is affine and $\sigma_{y,k}^{2}$ is convex quadratic in the control variables $\mathcal{U}_{t}$. Then $J$ is convex, each $g_{i}$ is convex, and each $h_{i}$ is affine in the optimization variables $\mathcal{U}_{t}$ and $\mathcal{Z}_{t}$.

This assumption holds in many applications of GP-MPC because the non-convexity of (10) usually comes solely from the GP dynamics [5, 7, 4]. It also usually holds in stochastic GP-MPC with chance constraints approximated by convex deterministic constraints (see, e.g., [7, 4]).

Because of the non-convexity of the general GP-MPC (10), a non-convex NLP solver is usually employed. However, as discussed in Section I and illustrated in Section V, the performance of these solvers on GP-MPC problems are usually unsatisfactory, especially in time-critical control applications.

IV linGP-based Model Predictive Control

In this section, an algorithm based on linGP is developed to solve the GP-MPC problem (10) with significantly improved performance compared to general NLP solvers. More importantly, its complexity is much less influenced by the size of the GP model (i.e., the size of the training data set).

IV-A Local GP-MPC problem with linGP

We first rewrite the GP-MPC problem in an equivalent form. Suppose that nominal control inputs $u^{\star}_{k}$, for $k\in\mathcal{I}_{t}=\{t,\dots,t+H-1\}$, are given. One can then simulate the GP model $\mathcal{G}_{f}$ to obtain the nominal outputs $\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu^{\star}_{k}$. As the vectors $x_{k}$ defined in (8) are constructed from the control inputs and GP outputs, we also obtain their nominal values $x^{\star}_{k}$. We note that the exact control inputs and process outputs before the current time $t$ are available (calculated or measured). With these nominal values, define $u_{k}=u^{\star}_{k}+\Delta u_{k}$ and $x_{k}=x^{\star}_{k}+\Delta x_{k}$, for all $k\in\mathcal{I}_{t}$. All $\Delta u_{k}$ are collected in $\Delta\mathcal{U}_{t}=\{\Delta u_{k}\,|\,k\in\mathcal{I}_{t}\}$. With a slight abuse of notations, we will write $\mathcal{U}_{t}=\mathcal{U}^{\star}_{t}+\Delta\mathcal{U}_{t}$, where $\mathcal{U}^{\star}_{t}$ collects all nominal control inputs. The GP-MPC formulation (10) is then equivalent to

[TABLE]

With this equivalent formulation, the GP-MPC problem can be approximated locally around the nominal values by replacing the GP model $\mathcal{G}_{f}$ with its linGP $\tilde{\mathcal{G}}_{f|(x^{\star}_{k},u^{\star}_{k})}$ at each $k$. Specifically, $\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{k}$ and $\sigma^{2}_{y,k}$ are replaced by the mean $\tilde{y}_{k}$ and variance $\tilde{\sigma}^{2}_{y,k}$ of the linGP at nominal inputs $x^{\star}_{k}$ and $u^{\star}_{k}$, as $\tilde{y}_{k}=\tilde{\mu}_{f|x^{\star}_{k},u^{\star}_{k}}(\Delta x_{k},\Delta u_{k})$ and $\tilde{\sigma}_{y,k}^{2}=\tilde{\sigma}_{f|x^{\star}_{k},u^{\star}_{k}}^{2}(\Delta x_{k},\Delta u_{k})$ respectively. It follows from the derivation of the linGP, particularly Eqs. (6), that $\tilde{y}_{k}$ is linear and $\tilde{\sigma}^{2}_{y,k}$ is convex quadratic in $\Delta x_{k}$ and $\Delta u_{k}$: $\tilde{y}_{k}=\hat{m}^{T}_{x^{\star}_{k},u^{\star}_{k}}\xi_{k}$ and $\tilde{\sigma}_{y,k}^{2}=\xi_{k}^{T}\hat{V}_{x^{\star}_{k},u^{\star}_{k}}\xi_{k}$, where $\xi_{k}^{T}=\left[1,\Delta x_{k}^{T},\Delta u_{k}^{T}\right]$. The vector $\hat{m}$ and the matrix $\hat{V}$ are defined in (4). Let $\Delta\tilde{y}_{k}=\tilde{y}_{k}-\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu^{\star}_{k}$. Then $\Delta x_{k}$ is derived from $\Delta\tilde{y}_{i}$ and $\Delta u_{i}$ for $i<k$, see (8), as

[TABLE]

As the approximate problem, in particular the linGP model, is only accurate in a neighborhood of the nominal GP inputs $x^{\star}_{k}$ and $u^{\star}_{k}$, we introduce a trust region of size $\rho\geq 0$. Specifically, we impose the bounds, for all $k\in\mathcal{I}_{t}$

[TABLE]

where $\left\lVert\cdot\right\rVert$ is an appropriate norm. Typically, box constraints, i.e., $\left\lVert\cdot\right\rVert_{\infty}$, are used. Note that since $x_{k}$ is derived from $u_{i}$ and $y_{i}$ with $i<k$, the trust region is defined on $u_{k}$ and $\tilde{y}_{k}$. The trust region size $\rho$ is adapted by the sequential convex programming algorithm presented in Section IV-B.

Putting everything together, the complete local problem, which we will call linGP-MPC, is given in (14) below.

[TABLE]

where $k\in\mathcal{I}_{t}$, $\tilde{\mathcal{Y}}_{t}$ and $\tilde{\mathcal{S}}_{y,t}$ respectively collect $\tilde{y}_{k}$ and $\tilde{\sigma}^{2}_{y,k}$.

By Assumption 1, the linGP-MPC problem (14) is convex and therefore can be solved efficiently with a suitable solver. By sequentially linearizing the GP with linGP, solving the convex linGP-MPC problem (14), and updating the nominal solution and trust region appropriately, the non-convex GP-MPC problem can be solved. The next subsection will present such a sequential convex programming (SCP) algorithm.

IV-B Sequential convex programming for linGP-MPC

An issue that can arise in the SCP approach is that the local linGP-MPC problem (14) can be infeasible even though the non-convex GP-MPC problem (10) is feasible. The nominal control inputs $u^{\star}_{k}$ are assumed to satisfy the control constraints (14f): $u^{\star}_{k}\in U$; otherwise we can project them onto $U$. However, the nominal values may not be feasible for the constraints (14g) and (14h) and, with a trust region bound $\rho$ small enough, problem (14) is infeasible. This issue often occurs during the early iterations of the SCP algorithm.

To address this issue, we will use exact penalty formulations of the GP-MPC and linGP-MPC problems. A penalized cost function is defined as

[TABLE]

where $\lambda>0$ is the penalty weight. If $\lambda$ is sufficiently large, the original GP-MPC problem can be solved by minimizing $\phi$ subject to the constraints (10b) to (10d). The exact penalty formulation of the linGP-MPC problem is

[TABLE]

Note that while both exact penalty formulations use the penalized cost function $\phi$ in (15), $\phi$ is evaluated at the exact GP outputs ($\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t}$ and $\mathcal{S}_{y,t}$) in GP-MPC, whereas in linGP-MPC, it is evaluated at the approximate linGP outputs ($\tilde{\mathcal{Y}}_{t}$ and $\tilde{\mathcal{S}}_{y,t}$).

We now present the linGP-SCP algorithm, inspired by the successive convexification algorithm in [9]. At iteration $j$, for $k\in\mathcal{I}_{t}$, let $u_{k}^{(j)}$ and $z_{k}^{(j)}$ be the current solution, $\mkern 2.0mu\overline{\mkern-2.0muy\mkern-1.0mu}\mkern 1.0mu_{k}^{(j)}$ and $x_{k}^{(j)}$ the corresponding outputs and states, $\rho^{(j)}$ the current trust region size, and $\phi^{(j)}$ the current exact penalized cost. We obtain the linGP model $\tilde{\mathcal{G}}_{f|x_{k}^{(j)},u_{k}^{(j)}}$, specifically the vectors $\hat{m}_{{x_{k}^{(j)},u_{k}^{(j)}}}$ and the SDP matrices $\hat{V}_{{x_{k}^{(j)},u_{k}^{(j)}}}$, for each $k$. The convex problem (16) is then solved, resulting in the approximate solution $(\tilde{\mathcal{Y}}_{t},\tilde{\mathcal{S}}_{y,t},\tilde{\mathcal{U}}_{t}=\mathcal{U}_{t}^{(j)}+\Delta\tilde{\mathcal{U}}_{t},\tilde{\mathcal{Z}}_{k})$. By simulating the GP model $\mathcal{G}_{f}$ with the inputs $\tilde{\mathcal{U}}_{t}$, we obtain the exact output mean $\mkern 2.0mu\overline{\mkern-2.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu_{t}$ and variance $\mathcal{S}_{y,t}$. To judge the algorithm progress, or the quality of the linGP approximation used in linGP-MPC, we compare the actual cost reduction

[TABLE]

to the predicted cost reduction

[TABLE]

If $|\tilde{\delta}^{(j)}|\leq\epsilon$, where $\epsilon>0$ is a predefined tolerance, the solution is considered converged and the algorithm is terminated. Otherwise, the ratio $r^{(j)}=\delta^{(j)}/\tilde{\delta}^{(j)}$ is inspected. Given three predefined thresholds $0<r_{0}<r_{1}<r_{2}<1$, there are four possibilities:

$r^{(j)}<r_{0}$: the approximation is considered too inaccurate, hence the approximate solution is rejected and the trust region size $\rho^{(j)}$ is contracted by a predefined factor $\beta^{\mathrm{fail}}<1$; 2. 2.

$r_{0}\leq r^{(j)}<r_{1}$: the approximation is deemed inaccurate but acceptable, hence the approximate solution is accepted but the trust region size $\rho^{(j)}$ is still contracted by $\beta^{\mathrm{fail}}$; 3. 3.

$r_{1}\leq r^{(j)}<r_{2}$: the approximation is considered sufficiently accurate, hence the approximate solution is accepted and the trust region size $\rho^{(j)}$ is retained; 4. 4.

$r^{(j)}\geq r_{2}$: the approximation is deemed highly accurate or even conservative, hence the approximate solution is accepted and the trust region size $\rho^{(j)}$ is enlarged by a predefined factor $\beta^{\mathrm{succ}}>1$.

The iteration is repeated until convergence or until a maximum number of iterations $j_{\mathrm{max}}>0$ is reached. The above linGP-SCP algorithm is detailed in Algorithm 1.

IV-C Complexity of linGP-SCP algorithm

The linGP-SCP algorithm 1 consists of a main loop and a convex local subproblem (i.e., the linGP-MPC problem). Because a general MPC formulation is assumed, resulting in a general convex subproblem, it is not possible to determine the specific complexity of the subproblem. However, we note that the convex subproblem can usually be solved efficiently. In this subsection, we focus instead on the complexity of the algorithm with respect to the size of the GP model (i.e., the size $N$ of its training data).

An important feature of the linGP-SCP algorithm is that the subproblem’s complexity is independent of the size of the GP model. This is obvious from inspecting the linGP-MPC subproblem in (14), where the parameters $\hat{m}$ and $\hat{V}$ do not depend on $N$. Indeed, the numerical example in Section V will show that the time spent on solving the subproblem is steady regardless of the size $N$ of the GP model. As a consequence, the inner iterations of the linGP-SCP algorithm, which solve the convex subproblem in each outer iteration of the main loop, do not depend on $N$.

The main loop involves, among other computations, the simulation and linearization of the GP model, which scale cubically with $N$. In fact, for a sufficiently large GP model, $N$ is the largest factor determining the complexity and the solving time of the GP-MPC problem. Therefore, we need to consider the number of iterations typically required by the linGP-SCP algorithm. Unlike standard trust-region algorithms, which usually perform a line search to achieve cost reduction, the linGP-SCP algorithm optimally solves the convex subproblem to reduce the cost in each main succession. Hence, the number of inner iterations in each outer iteration is increased, while the number of outer iterations is decreased. In other words, by achieving a greater cost reduction in each outer iteration, at the expense of increased subproblem solving time, the linGP-SCP algorithm can reduce the number of outer iterations. Because the inner iterations are independent of $N$ while the outer iterations scale cubically with $N$, this feature of the linGP-SCP algorithm helps reduce the overall complexity and solving time significantly, particularly for large GP models (i.e., large $N$).

Compared to SQP-based methods, the linGP-SCP algorithm has a crucial advantage, stemming from the use of linGP models to approximate the original GP model. An SQP-based method typically uses a computationally expensive technique such as the BFGS update to estimate the Hessian of the cost function and probably other functions. Moreover, because the exact GP mean and variance, as computed by (1), are used directly, their Hessians need to be computed, whose computational complexity grow exponentially with $N$. By contrast, the linGP approximation, as expressed in (4) and (6), relies only on the Jacobians of the covariance function and therfore is much cheaper to compute. Another distinction between the linGP-SCP algorithm and SQP-based methods is that the linGP approximation, as a valid GP in itself, always guarantees the positive definiteness of its variance (cf. (5) and the paragraph following it). Whereas there is no such guarantee in the direct approximation of the GP variance used by SQP-based methods. Therefore, SQP-based methods must employ additional measures to ensure the positive definiteness of the estimated Hessian so that the SQP subproblem is convex, incurring more computational cost.

In summary, the use of the linGP approximation, with first-order computation and guaranteed convexity, and its resulting convex linGP-MPC problem in the linGP-SCP algorithm leads to a significant reduction in computational complexity compared to other methods. This advantage is demonstrated in the numerical example in the next section.

V Numerical Example

To demonstrate the efficacy of the linGP-SCP algorithm, we considered a GP-MPC problem of a nonlinear dynamic system, adapted from [5, Example 4.3, p 183]. The process to be controlled is a discrete-time scalar system with noisy output of the form:

[TABLE]

where $w_{k}$ is a white noise with standard deviation $0.025$ and zero mean. The following input constraints are imposed

[TABLE]

This nonlinear process was learned by GP models using the autoregressive input $[y_{k-1},u_{k}]$ and the standard squared exponential covariance function (see [5, 6]). The training data were obtained by feeding a random input signal that satisfies the constraints (20) to the process [5]. To illustrate the performance of solving the GP-MPC problem with respect to the training data size $N$, we trained 15 GP models with $N$ incremented by 100 from 100 to 1,500. It is worth pointing out that a training data set larger than $N=400$ did not improve the GP model accuracy significantly; however, we considered large values of $N$ to illustrate the computational benefits of our algorithm.

The control goal is to track a given reference signal $r_{k}$. The GP-MPC problem is formulated as follows [5]

[TABLE]

The state constraints (21c) and terminal constraints (21d) and (21e) are deterministic approximations of stochastic constraints that, with high probability, bound the process state. The MPC parameters are $H=12$, $Q=10$, $R=0.1$, and $\delta=0.075$. The closed-loop control system was simulated for 100 steps with the following reference signal:

[TABLE]

All simulations were run in Matlab R2017b on an iMac with Intel Core i7 4.2 GHz processor and 16 Gb RAM.

V-A Solving the GP-MPC

Figure 1 shows a box plot of the solving time of the GP-MPC problem (21) for increasing training data size and with three NLP methods: IPM, SQP, and active-set. In all cases, the solving time grows quickly with the training data size. While the IPM solving time does not vary much between time steps (i.e., low variance), the SQP and active-set methods have noticeably varying solving time (i.e., high variance with far-away outliers). In other words, the performance of IPM is steady though slow, while SQP and active-set methods are unpredictable in terms of performance. For these reasons, general NLP methods are not well suited for real-time GP-MPC applications, where the optimization solving performance must be both fast and consistent.

To compare the performance of the linGP-SCP algorithm with state-of-the-art solvers, we solved the GP-MPC problem (21), for all 15 GP models, with different solvers:

• •

Ipopt: the GP-MPC was formulated in CasADi [12] for Matlab and solved by the open-source solver Ipopt [13] with the optimized linear solver MA57. All Jacobians and Hessians were calculated automatically by CasADi in C code, well-known to be very fast.

• •

Knitro: similar to the previous, however the commercial nonlinear solver Knitro was used. Knitro is commonly regarded as one of the best nonlinear solvers.

• •

linGP-SCP: the linGP-SCP algorithm was prototyped purely in Matlab, where the subproblem was modeled by Yalmip [14] and solved by Gurobi.

The solvers were configured to similar accuracy levels; indeed, given the same parameters, their solutions were always within $1\text{\times}{10}^{-6}$ from each other. We make two remarks.

Remark 3

We also used the commercial solver fmincon, included in Matlab, with its SQP and active-set methods. However, as discussed earlier and illustrated in Figure 1, its performance was unpredictable and often worse than other solvers. Therefore we will not report its results here.

Remark 4

Unlike the highly optimized solvers Ipopt and Knitro, which are usually implemented in C/C++ and/or Fortran, the linGP-SCP prototype in Matlab is certainly not an optimized implementation. In addition, the Yalmip toolbox is known to be slower than the C-based CasADi library. Therefore, there is significant room to improve the performance of the linGP-SCP implementation.

V-B Simulation results and discussions

In all cases, the GP-MPC controller was able to track the given reference. Because our main interest is in the solving performance, we omit the tracking control results and only report the performance results of the used solvers.

Figure 2 is a box plot of the solving time of each MPC step by Knitro and linGP-SCP, for all 15 GP models. We do not report the timing results of Ipopt in Figure 2 because they were significantly worse than the other two. Observe that linGP-SCP’s solving time was highly consistent with low variance and few small outliers.

Most importantly, while the solving time of Knitro grew quickly with the GP model size $N$, the solving time of linGP-SCP grew much slower. This observation was more obvious in Figure 3, which plots the median solving time. Clearly, the solving time of Ipopt and Knitro grew exponentially with $N$ while the solving time of linGP-SCP did not increase significantly. Figure 4 shows the relative ratio of the median solving time of Ipopt and Knitro to that of linGP-SCP. It can be seen that Ipopt and Knitro became increasingly slower than linGP-SCP when $N$ increased. These results confirm the complexity analysis in Section IV-C that, compared to other methods, the complexity of linGP-SCP is less affected by the GP model size and therefore grows slowly with $N$.

In Table I, the mean and standard deviation of the execution time (in seconds) per MPC step for several model sizes are listed. Again, these numbers support the previous claim about the performance advantage of the linGP-SCP algorithm. We can also observe that the solving time of the linGP-MPC subproblem (the “linGP-inner” row in the table) was constant regardless of the GP model size, as stated in the complexity analysis in Section IV-C.

VI Conclusion

We developed a fast and scalable approach to solve a general class of GP-based MPC problems. The approach is based on the concept of linearized Gaussian Process (linGP), proposed in this paper, and Sequential Convex Programming. Our approach not only solves GP-MPC faster than other NLP methods but also is much less influenced by the GP training data size – a key factor affecting the computational complexity of GPs and GP-MPC. Therefore, our algorithm is more scalable and predictable, which makes it more suitable for real-time control.

We are developing a more efficient implementation of the method and will conduct real-time control experiments. We are also studying the convergence properties of the linGP-SCP algorithm, and investigating different applications of the linGP concept other than MPC.