Linearized Gaussian Processes for Fast Data-driven Model Predictive Control
Truong X. Nghiem


[
    Hyunsuk Hong
    Received: date / Accepted: date
===================================

empty
empty




  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.























§ 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 <cit.>.
Stochastic MPC (SMPC) 
can efficiently handle stochastic uncertainties with underlying probability distributions in optimal control problems <cit.>.
For these reasons, MPC and its variants have been applied widely in many practical control problems <cit.>.


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 <cit.>.
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 <cit.>.
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 <cit.>.
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 
<cit.>.
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 <cit.>.
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 <cit.>. 
However, the performance of these 
methods in GP-MPC problems is usually unsatisfactory.
For instance, it was reported in <cit.> that a moderate-sized GP-MPC problem took up to a minute to solve with Ipopt.
Section <ref> 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 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 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 <cit.> to solve non-convex optimal control problems and similar in concept to successive linearization methods to solve nonlinear MPC (such as in <cit.>).
However, our approach is different from other attempts to accelerate GP-MPC and machine-learning-based MPC in general, for example <cit.>, in that we approximate not the output equation of the machine learning model but the underlying latent process of the model.   



The concept of is introduced in Section <ref>, followed by the considered class of GP-MPC problems in Section <ref>.
We present the -based SCP algorithm in Section <ref>.
Finally, the efficacy of our method is demonstrated in a numerical example in Section <ref>.











§ LINEARIZED GAUSSIAN PROCESSES (LINGP)


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



 §.§ Gaussian Processes: Basics and Notations









Given N noisy observations y^(i) of an underlying function f: ^n ↦ through a Gaussian noise model: y^(i) = f(x^(i)) + ϵ^(i), with inputs x^(i)∈^n, ϵ^(i)∼0σ_n^2, and i = 1,…,N .
A GP of f, considered as a distribution over functions <cit.>, is essentially a probability distribution on the observations of all possible realizations of f.
The prior (f) represents the initial belief about this distribution.
Let = (X, Y) be the finite set of observation data of f, where X = [x^(1), …, x^(N)] ∈^n × N collects the regression vectors and Y = [y^(1), …, y^(N)] ∈^N contains the corresponding observed outputs.
By conditioning the GP on , the posterior (f | ) 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), …, f^(M)} corresponding to inputs {x^(1), …, x^(M)}.
It is fully specified by its mean function m(x) and covariance function k(x,x'), parameterized by the hyperparameters θ,

    m(x; θ)    =  [f(x)] 
    
    k(x,x'; θ)    =  [(f(x)-μ(x)) (f(x') -μ(x'))] .


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


The output f_⋆ of the GP corresponding to an input x_⋆ is a random variable f_⋆∼μ_⋆σ_⋆^2, where



    μ_⋆   = m(x_⋆) + K_⋆ K^-1 (Y - μ(X)) 
    σ_⋆^2    = K_⋆⋆ - K_⋆ K^-1 K_⋆^T ,


in which K_⋆ = [k(x_⋆, x^(1)), …, k(x_⋆, x^(N))] ∈^N, K_⋆⋆ = k(x_⋆, x_⋆) ∈, and K ∈^N × N is the covariance matrix with elements K_ij = k(x^(i), x^(j)).
Note that the mean and variance are nonlinear in x_⋆ and their computations scale cubically with the size of X and Y.

In theory, 
the hyperparameters θ are also random variables, 
whose posterior distributions are obtained by conditioning them on  using the Bayes' theorem.
In practice, however, θ are often obtained by maximizing the likelihood: _θ(Y | X, θ).



For an in-depth treatment of GPs, interested readers are referred to <cit.>.




 §.§ Linearized Gaussian Processes



Consider a GP of a nonlinear function f(x): ^n→, with covariance function k(x, x') : ^n×^n→.
We will denote the GP of f by [f], whose predictive mean and variance at an input x ∈^n are given by (<ref>) and are denoted by f(x) and f(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 (<ref>), is non-convex due to the covariance function k(·, ·).
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 <ref>, 
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(·, ·) is differentiable, which is typically the case <cit.>.
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 

    f(x + Δ_x) ≈f̃_x(Δ_x) = f(x) + Δ_x^T∇_x f(x) ,   Δ_x∈^n .




Define g ≡∇_x f as the gradient of f.
The above linearization can be rewritten as, with x̂ = [1, Δ_x^T]^T,

    f̃_x(Δ_x) =
      [   1; Δ_x ]^T[ f(x); g(x) ] =
      x̂^Tf̂(x) .

It is important to note here that
f̂(x) = [f(x), g^T(x)]^T
is a random vector of length (n+1) which, through the inner product with the vector x̂, results in the random variable f̃_x(Δ_x) of the linearized process.
In other words, we approximate the original GP of f around x by a process of linearized function f̃_x of f, which we will call a linearized Gaussian Process or linGP.
The linGP of the GP [f] at input x is denoted by fx.

As differentiation is a linear operator, f̂ is a (multivariate) Gaussian Process derived from the original GP.
Therefore, f̂(x) is a Gaussian random vector that defines a posterior distribution, derived from the posterior (f | X,Y), over the values of f(x) and its gradiant at x.
Consequently, the linGP fx 
has a distribution derived from (f | X,Y) over the linear approximations of the function f near x .
To characterize this linGP, 
we need to characterize the GP of f̂.

Let K^(1,0) = (∇_x k) – the gradient of k with respect to the first argument – and K^(0,1) = (D_x' k) – the Jacobian of k with respect to the second argument.
Let K^(1,1)(x,x') be the n × n matrix such that K^(1,1)_i,j(x,x') = ∂^2/∂ x_i∂ x'_j k(x,x').
Following <cit.>, we can write the joint distribution of the observation data, the process f, and the gradient process g as

    [
        Y 
     f 
     g
      ] ∼0[[ K(X,X) + σ_n^2         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) ]].

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 × N matrix with K_i,j(X,X) = k(x^(i), x^(j)); similarly, the i^th row of K^(0,1)(X,x) is K^(0,1)(x^(i), x).
The posterior distribition of f̂ conditioned on the observation data (X,Y) is then

    f̂ | X,Y,x ∼m̂_xV̂_x
    m̂_x =
      [       k(x,X); K^(1,0)(x,X) ]( K(X,X) + σ_n^2)^-1 Y 
    V̂_x =    [[       k(x,x) K^(0,1)(x,x); K^(1,0)(x,x) K^(1,1)(x,x) ]] - 
       [[       k(x,X); K^(1,0)(x,X) ]]
        ( K(X,X) + σ_n^2)^-1[[       k(X,x) K^(0,1)(X,x) ]]

where  denotes an identity matrix of appropriate dimensions.
The linearized process in (<ref>) can then be written as follows.

    f̃_x(Δ_x) ∼(
        m̂_x^Tx̂,
        x̂^TV̂_x x̂) .

The mean of the original GP at (x + Δ_x) is approximated by a linear function of x̂, while its variance is approximated by a quadratic function of x̂.
Because f̂ is a valid GP, V̂_x is positive semi-definite (PSD) and, therefore, it is guaranteed that the 
variance of f̃_x(Δ_x) is non-negative.
We emphasize again that the linGP is in fact a Gaussian Process of 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 fx(Δ_x) and fx(Δ_x) as the notations for the mean and variance of the linGP fx at input (x+Δ_x) in a neighborhood of x.
They are calculated as

    fx(Δ_x) = m̂_x^Tx̂,   fx(Δ_x) = x̂^TV̂_x x̂,

in which m̂_x and V̂_x are given in (<ref>), and x̂ is defined in (<ref>).











  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 x̂, K, K^(0,1), K^(1,0), and K^(1,1) corresponding to the fixed elements of x are removed.



  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 is that while the former linearizes the mean and variance of the probabilistic predictive output (linearization happens after prediction), linearizes the latent function directly and derives the probabilistic distribution of such linearized functions (linearization happens before prediction).
  The mean fx obtained above is exactly the same as if the exact mean equation (<ref>) is linearized directly.
  However, the variance fx is not the same as if the exact variance equation (<ref>) is approximated to the second order by a Taylor series.
  Indeed, while fx 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 <ref>, this key difference is an important advantage of the approach, that could lead to a substantial performance improvement of our algorithm compared to other methods.  
  



























§ GP-BASED MODEL PREDICTIVE CONTROL


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




 §.§ Gaussian Processes for Dynamical Systems



GPs can be used for modeling nonlinear dynamical systems <cit.>.
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

    y_k = f(y_k-l, …, y_k-1, u_k-m, …, u_k) .

Here, k denotes the time step, u ∈^n_u the control input, 
y ∈ 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 (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 [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 [f] over H future steps, where H is the MPC horizon, and predict its multistep-ahead behavior.
Because the output of [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 [f] involves the propagation of uncertainty through the model.
There exist several methods for approximating uncertainty propagation in GPs <cit.>.

It was shown in <cit.> 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,

    x_k^T = [y_k-l, y_k-l+1,…, y_k-1,     
    
      u_k-m^T, u_k-m+1^T, …, u_k-1^T]

where y_i is the expected value of y_i.
Note that all y_k for k < t, where t is the current time, can be observed. 

The output of [f] at time k ≥ t then reads


    y_k∼y_kσ_y,k^2
    y_k = f(x_k, u_k),   σ_y,k^2 = f(x_k, u_k) .











 §.§ 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 _t = {y_k |  k ∈_t}, _y,t = {σ_y,k |  k ∈_t}, and _t = {u_k |  k ∈_t}, 
where _t = {t,…,t+H-1} contains the indices of all time steps in the current MPC horizon.
We also introduce additional variables z ∈^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 _t = {z_k |  k ∈_t}.

The GP-MPC reads



    (GP-MPC)
       _t, _tminimize   J(_t, _y,t, _t, _t) 
       subject to
           y_k = f(x_k, u_k) 
           σ_y,k^2 = f(x_k, u_k) 
            z_k∈ Z,    u_k∈ U 
            g_i(_t, _y,t, _t, _t) ≤ 0,   i = 1,…,n_ieq
            h_i(_t, _y,t, _t, _t) = 0,   i = 1,…,n_eq


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


  * J(·) is the objective function and typically has the form
  J = ∑_k=t^t+H-1 c_k( y_k, σ_y,k, u_k, z_k)
  where c_k(·) are the stage cost functions;

  * Constraints (<ref>) and (<ref>) represent the GP dynamics;

  * Sets Z and U in (<ref>) are convex constraint sets of the variables z and u, respectively;

  * There are n_ieq inequality constraints of the form (<ref>);

  * There are n_eq equality constraints of the form (<ref>);



We make the following assumptions about the GP-MPC formulation (<ref>).

  
  Suppose that, for every k, y_k is affine and σ_y,k^2 is convex quadratic in the control variables _t.
  Then J is convex, each g_i is convex, and each h_i is affine in the optimization variables _t and _t.

This assumption holds in many applications of GP-MPC because the non-convexity of (<ref>) usually comes solely from the GP dynamics <cit.>.
It also usually holds in stochastic GP-MPC 
with chance constraints approximated by convex deterministic constraints (see, <cit.>).

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











§ LINGP-BASED MODEL PREDICTIVE CONTROL


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




 §.§ Local GP-MPC problem with 



We first rewrite the GP-MPC problem in an equivalent form.
Suppose that nominal control inputs u^⋆_k, for k ∈_t = {t,…,t+H-1}, are given.
One can then simulate the GP model [f] to obtain the nominal outputs y^⋆_k.
As the vectors x_k defined in (<ref>) are constructed from the control inputs and GP outputs, we also obtain their nominal values x^⋆_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^⋆_k + Δ u_k and x_k = x^⋆_k + Δ x_k, for all k ∈_t.
All Δ u_k are collected in Δ_t = {Δ u_k |  k ∈_t}.
With a slight abuse of notations, we will write _t = ^⋆_t + Δ_t, where ^⋆_t collects all nominal control inputs.
The GP-MPC formulation (<ref>) is then equivalent to



    Δ_t, _tminimize   J(_t, _y,t, ^⋆_t + Δ_t, _t) 
       subject to
         y_k = f(x^⋆_k + Δ x_k, u^⋆_k+Δ u_k) 
         σ_y,k^2 = f(x^⋆_k + Δ x_k, u^⋆_k+Δ u_k) 
          z_k∈ Z,    u^⋆_k+Δ u_k∈ U 
          g_i(_t, _y,t, ^⋆_t + Δ_t, _t) ≤ 0,   i = 1,…,n_ieq
          h_i(_t, _y,t, ^⋆_t + Δ_t, _t) = 0,   i = 1,…,n_eq



With this equivalent formulation, the GP-MPC problem can be approximated locally around the nominal values by replacing the GP model [f] with its f(x^⋆_k, u^⋆_k) at each k.
Specifically, y_k and σ^2_y,k are replaced by the mean ỹ_k and variance σ̃^2_y,k of the at nominal inputs x^⋆_k and u^⋆_k, as
ỹ_k = fx^⋆_k, u^⋆_k(Δ x_k, Δ u_k)
and
σ̃_y,k^2 = fx^⋆_k, u^⋆_k(Δ x_k, Δ u_k)
respectively.
It follows from the derivation of the linGP, 
particularly Eqs. (<ref>), that ỹ_k is linear and σ̃^2_y,k is convex quadratic in Δ x_k and Δ u_k:

ỹ_k = m̂^T_x^⋆_k, u^⋆_kξ_k 
and

σ̃_y,k^2 = ξ_k^TV̂_x^⋆_k, u^⋆_kξ_k,

where
ξ_k^T = [ 1, Δ x_k^T, Δ u_k^T].
The vector m̂ and the 
matrix V̂ are defined in (<ref>).

Let Δỹ_k=ỹ_k - y^⋆_k.
Then Δ x_k is derived from Δỹ_i and Δ u_i for i < k, see (<ref>), as

    Δ x_k^T = [Δỹ_k-l, Δỹ_k-l+1,…, Δỹ_k-1,     
    Δ u_k-m^T, Δ u_k-m+1^T, …, Δ u_k-1^T].



As the approximate problem, in particular the model, is only 
accurate in a neighborhood of the nominal GP inputs x^⋆_k and u^⋆_k, we introduce a trust region of size ρ≥ 0.
Specifically, we impose the bounds, for all k ∈_t

    Δ u_k≤ρ,   Δỹ_k≤ρ

where · is an appropriate norm.
Typically, box constraints, ·_∞, 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 ỹ_k.
The trust region size ρ is adapted by the sequential convex programming algorithm presented in Section <ref>. 

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



    (linGP-MPC)
       Δ_t, _tminimize   J(_t, _y,t, ^⋆_t + Δ_t, _t) 
       subject to
         ỹ_k = m̂_x^⋆_k, u^⋆_k^Tξ_k,  
          σ̃_y,k^2 = ξ_k^TV̂_x^⋆_k, u^⋆_kξ_k
         ξ_k^T = [ 1, Δ x_k^T, Δ u_k^T] 
         Δỹ_k = ỹ_k - y_k^⋆,   Δ x_k defined in (<ref>)
         Δ u_k≤ρ,   Δỹ_k≤ρ
          z_k∈ Z,    u^⋆_k+Δ u_k∈ U 
          g_i(_t, _y,t, ^⋆_t + Δ_t, _t) ≤ 0,   i = 1,…,n_ieq
          h_i(_t, _y,t, ^⋆_t + Δ_t, _t) = 0,   i = 1,…,n_eq


where
k∈_t, 
_t and _y,t respectively collect ỹ_k and σ̃^2_y,k.




By Assumption <ref>, the linGP-MPC problem (<ref>) is convex 
and therefore can be solved efficiently with a suitable solver.
By sequentially linearizing the GP with , solving the convex linGP-MPC problem (<ref>), 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.




 §.§ Sequential convex programming for linGP-MPC



An issue that can arise in the SCP approach is that the local linGP-MPC problem (<ref>) can be infeasible even though the non-convex GP-MPC problem (<ref>) is feasible.
The nominal control inputs u^⋆_k are assumed to satisfy the control constraints (<ref>): u^⋆_k∈ U; otherwise we can project them onto U.
However, the nominal values may not be feasible for the constraints (<ref>) and (<ref>) and, with a trust region bound ρ small enough, 
problem (<ref>) 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

    ϕ(_t, _y,t, _t, _t) = 
      J(_t, _y,t, _t, _t) + 
    λ(
        ∑_i=1^n_ieqmax( 0, g_i(_t, _y,t, _t, _t)) + 
    ∑_i=1^n_eqh_i(_t, _y,t, _t, _t))

where λ > 0 is the penalty weight.
If λ is sufficiently large, the original GP-MPC problem can be solved by minimizing ϕ subject to the constraints (<ref>) to (<ref>).










The exact penalty formulation of the linGP-MPC problem is

    (linGP-MPC exact penalty)
       Δ_t, _tminimize  ϕ(_t, _y,t, ^⋆_t + Δ_t, _t) 
       subject to constraints (<ref>) to (<ref>).


Note that while both exact penalty formulations use the penalized cost function ϕ in (<ref>), ϕ is evaluated at the exact GP outputs (_t and _y,t) in GP-MPC, whereas in linGP-MPC, it is evaluated at the approximate outputs (_t and _y,t).

We now present the linGP-SCP algorithm, inspired by the successive convexification algorithm in <cit.>.
At iteration j, for k ∈_t, let u_k^(j) and z_k^(j) be the current solution, y_k^(j) and x_k^(j) the corresponding outputs and states, ρ^(j) the current trust region size, and ϕ^(j) the current exact penalized cost.
We obtain the linGP model fx_k^(j), u_k^(j), specifically the vectors m̂_x_k^(j), u_k^(j) and the SDP matrices V̂_x_k^(j), u_k^(j), for each k.
The convex problem (<ref>) is then solved, resulting in the approximate solution (_t, _y,t, _t = _t^(j) + Δ_t, _k).
By simulating the GP model [f] with the inputs _t, we obtain the exact output mean _t and variance _y,t.
To judge the algorithm progress, or the quality of the approximation used in linGP-MPC, we compare the actual cost reduction

    δ^(j) = ϕ^(j) - ϕ(_t, _y,t, _t, _t)

to the predicted cost reduction

    δ̃^(j) = ϕ^(j) - ϕ(_t, _y,t, _t, _t) .

If |δ̃^(j)| ≤ϵ, where ϵ > 0 is a predefined tolerance, the solution is considered converged and the algorithm is terminated.
Otherwise, the ratio r^(j) = δ^(j) / δ̃^(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 ρ^(j) is contracted by a predefined factor β^fail < 1;

  * r_0≤ r^(j) < r_1: the approximation is deemed inaccurate but acceptable, hence the approximate solution is accepted but the trust region size ρ^(j) is still contracted by β^fail;

  * r_1≤ r^(j) < r_2: the approximation is considered sufficiently accurate, hence the approximate solution is accepted and the trust region size ρ^(j) is retained;

  * r^(j)≥ r_2: the approximation is deemed highly accurate or even conservative, hence the approximate solution is accepted and the trust region size ρ^(j) is enlarged by a predefined factor β^succ > 1.

The iteration is repeated until convergence or until a maximum number of iterations j_max>0 is reached.
The above linGP-SCP algorithm is detailed in Algorithm <ref>.






 §.§ Complexity of linGP-SCP algorithm



The linGP-SCP algorithm <ref> consists of a main loop and a convex local subproblem (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 (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 (<ref>), where the parameters m̂ and V̂ do not depend on N.
Indeed, the numerical example in Section <ref> 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 (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 (<ref>), are used directly, their Hessians need to be computed, whose computational complexity grow exponentially with N.
By contrast, the linGP approximation, as expressed in (<ref>) and (<ref>), 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 ((<ref>) 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.










§ NUMERICAL EXAMPLE
  

To demonstrate the efficacy of the linGP-SCP algorithm, we considered a GP-MPC problem of a nonlinear dynamic system, adapted from <cit.>.

The process to be controlled is a discrete-time scalar system with noisy output of the form:

    x_k+1 = x_k - 0.5 tanh(x_k + u_k^3),    y_k = x_k + w_k

where w_k is a white noise with standard deviation 0.025 and zero mean.
The following input 
constraints are imposed

    -1 ≤ u_k≤ 1,    -0.5 ≤ u_k+1 - u_k≤ 0.5,   ∀ k .


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 <cit.>).
The training data were obtained by feeding a random input signal that satisfies the constraints (<ref>) to the process <cit.>.
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 <cit.>



    u_t, …, u_t+H-1minimize ∑_k=t^t+H-1( Q (y_k - r_t)^2 + R (u_k-u_k-1)^2) 
    
         subject to
         y_k = f(y_k-1, u_k),   σ_y,k^2 = f(y_k-1, u_k) 
          -1 ≤ u_k≤ 1,    -0.5 ≤ u_k+1 - u_k≤ 0.5 
          -1.2 ≤y_k - 2 σ_y,k,   y_k + 2 σ_y,k≤ 1.2 
          -δ≤y_t+H-1 - r_t - 2 σ_y,t+H-1
         y_t+H-1 - r_t + 2 σ_y,t+H-1≤δ


The state constraints (<ref>) and terminal constraints (<ref>) and (<ref>) 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 δ = 0.075.

The closed-loop control system was simulated for 100 steps with the following reference signal:

    r_k = -0.5,   k ∈ [0, 50];    r_k = -0.2,   k ∈ [51, 100].



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





 §.§ Solving the GP-MPC




Figure <ref> shows 
a box plot of the solving time of the GP-MPC problem (<ref>) 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 (low variance), the SQP and active-set methods have noticeably varying solving time (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 -SCP algorithm with state-of-the-art solvers, we solved the GP-MPC problem (<ref>), for all 15 GP models, with different solvers:


  * Ipopt: the GP-MPC was formulated in CasADi <cit.> for Matlab and solved by the open-source solver Ipopt <cit.> 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 <cit.> and solved by Gurobi. 

The solvers were configured to similar accuracy levels; indeed, given the same parameters, their solutions 
were always within 1e-6 from each other.
We make two remarks.

  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 <ref>, its performance was unpredictable and often worse than other solvers.
  Therefore we will not report its results here.


  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.
  





 §.§ 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 <ref> 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 <ref> 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 <ref>, 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 <ref> 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 <ref> 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 <ref>, 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 <ref>.

  










§ 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 (), 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 concept other than MPC.











IEEEtran






























