Gaussian Process Model Predictive Control of Unknown Nonlinear
Systems
Gang Cao1,* , Edmund M-K Lai2 , Fakhrul Alam1

arXiv:1612.01211v1 [] 5 Dec 2016

1

School of Engineering and Advanced Technology, Massey University, Auckland, New Zealand.
Department of Information Technology and Software Engineering, Auckland University of Technology, Auckland, New Zealand.
*
g.cao@massey.ac.nz
2

Abstract: MPC of an unknown system that is modelled by GP techniques is studied in this paper.
Using GP, the variances computed during the modelling and inference processes allow us to take
model uncertainty into account. The main issue in using MPC to control systems modelled by GP
is the propagation of such uncertainties within the control horizon. In this paper, two approaches
to solve this problem, called GPMPC1 and GPMPC2, are proposed. With GPMPC1, the original SMPC problem is relaxed to a deterministic nonlinear MPC based on a basic linearized GP
local model. The resulting optimization problem, though non-convex, can be solved by the SQP.
By incorporating the model variance into the state vector, an extended local model is derived. This
model allows us to relax the non-convex MPC problem to a convex one which can be solved by
an active-set method efficiently. The performance of both approaches is demonstrated by applying
them to two trajectory tracking problems. Results show that both GPMPC1 and GPMPC2 produce
effective controls but GPMPC2 is much more efficient computationally.

1.

Introduction

MPC, also known as receding horizon control, is a class of computer control algorithms that predicts future responses of a plant based on its system model, and computes optimized control inputs
by repeatedly solving a finite horizon optimization problem [1]. The advantages of MPC mainly
lie in its conceptual simplicity for multiple variable problems, and its ability to handle input and
output “hard-constraints” that are commonly encountered in practice but are not well addressed by
other control methods. It has been applied to many different types of control problems [2, 3].
The performance of MPC is highly dependent on the accuracy of the system model that describes its dynamics. Traditionally, these models are derived mathematically. More recently, datadriven modelling approaches based on computational intelligence and machine learning techniques
are becoming popular [4, 5]. This approach is especially suitable for complex and highly nonlinear systems where complete knowledge of the system dynamics is seldom available, giving rise
to unmodelled dynamics or model uncertainty. From the MPC perspective, attempts to address
the issue of model uncertainty has been made through RMPC schemes such as open-loop “minmax” MPC [6], closed-loop “min-max” MPC [7] and tube-based MPC [8]. “Min-max” MPC is
conceptually simple. However, its control laws are computed based on worst-case scenarios and
are therefore considered too conservative. Tube-based MPC overcomes this problem by combining a conventional MPC for the nominal system and a local feedback control law that steers the

states of the unknown system to the inside of a “tube” centered on the nominal trajectory [9]. This
“tube”, which relates to the uncertainty bounds, must be carefully defined. Otherwise, there may
not be a feasible solution. The major problem with RMPC is that model uncertainties are assumed
to be deterministic even though they are typically stochastic.
SMPC is an alternative where model uncertainties are assumed to be stochastic with an underlying probability distribution [10, 11, 12, 13]. Control laws are computed by solving a stochastic
optimization problem. Furthermore, since the state or output constraints are also probabilistic, they
can be satisfied with a predefined level of confidence. This effectively alleviates the conservatism
of “min-max” MPC. Furthermore, it is possible to trade-off control performance with robustness
against model uncertainties by adjusting these probabilistic constraints. A key problem with SMPC
is the propagation of uncertainties over a finite prediction horizon. The most common solution is
to use sampling-based MC simulation techniques. However, they are computationally demanding.
More recently, a technique known as polynomial chaos expansions has been proposed to lighten
the computation burden [14, 13].
A model known as GP has become very useful in statistical modelling [15]. The GP variances
which are computed as part of the modelling process provide a useful indication of the accuracy
of the model. These variances can also be propagated in multiple-step ahead predictions. The
hyperparameters of these models are learnt from data by maximizing the log-likelihood function.
This optimization problem is unconstrained, nonlinear and non-convex optimization. It is typically
solved by CG [15] or by PSO techniques [16, 17, 18].
A GP based MPC scheme was first introduced in [19]. Subsequently, an SMPC scheme using
GP was proposed in [20]. Even though GP is a probabilistic model, the cost functions used in these
papers are deterministic. Consequently, the variances could only be treated as slack variables of
the state constraints. This indirect way of handling GP variances leads to a nonlinear optimization
problem that is very computationally demanding to solve. More recently, in [21, 22, 23], variances
are included in the cost function and can be directly handled in the optimization process. However,
only unconstrained MPC have been considered.
In this paper, two new GP based MPC approaches, referred to as GPMPC1 and GPMPC2, are
proposed for the control of unknown nonlinear dynamical systems with input and state constraints.
The GPMPC1 approach is similar to those in [19, 20] in the sense that the GP variances are considered as a slack variable in the state constraints. The main difference is that the resulting non-convex
optimization problem is solved by using a SQP based method together with a linearized GP model
which is called the basic GP based local model in this paper. The constrained stochastic problem
is then relaxed to a deterministic one by specifying the confidence level. With GPMPC2, the nonlinear MPC problem is reformulated to a convex optimization problem. In contrast with earlier
methods, GP variances are directly included in the cost function of the optimization. The solution
method makes use of a modified version of the basic local model which includes the variance in
the modified state variable of the system. The resulting MPC problem is efficiently solved by using
an active-set method. The effectiveness of these approaches are demonstrated by applying them to
two trajectory tracking problems.
The rest of this paper is organized as follows. Section 2 introduces the modelling of the unknown nonlinear system by using GP models. The basic and extended GP based local dynamical
models are presented in the Section 3. In Section 4, the proposed GPMPC1 and GPMPC2 are presented for the general trajectory tracking problem of the unknown nonlinear system. In addition,
the feasibility and stability of the proposed algorithms are also discussed. The simulation results
are next reported to demonstrate the performance of the proposed algorithms in Section 5. Finally,

Section 6 draws the conclusions.
2.

Unknown system modelling using GP

Consider a discrete-time nonlinear dynamical system described by the following general form:
xk+1 = f (xk , uk ) + wk

(1)

where f : Rn × Rm → Rn is a nonlinear function, w ∈ Rn represents additive external disturbances, x ∈ Rn denotes the state vector, and u ∈ Rm are control signals. In this paper, we assume
that f is totally unknown but can be represented by a GP model. The uncertainty of a GP model
can be measured by the GP variances. Therefore, a disturbance observer will not be required. The
hyperparameters of a GP model is learnt from a set of training data consisting of inputs to the
system and the system’s response as target.
To model a system given by (1), a natural choice of the model inputs and their targets are
the state-control tuple x̃k = (xk , uk ) ∈ Rn+m and the next state xk+1 respectively. Let ∆xk =
xk+1 − xk ∈ Rn . In practice, the variation between ∆xk+1 and ∆xk is much less the variation
between xk+1 and xk , for all k. Therefore it is more advantageous to use ∆xk as the model target
instead [24]. This will be assumed in the rest of this paper.
2.1.

GP Modelling

A GP model is completely specified by its mean and covariance function [15]. Assuming that the
mean of the model input x̃k is zero, the squared exponential covariance is given by K(x̃i , x̃j ) =
1
σs2 exp(− (x̃i − x̃j )T Λ(x̃i − x̃j )) + σn2 . The parameters σs2 , σn2 and the entries of matrix Λ are
2
referred to as the hyperparameters θ of a GP model. Given D training inputs X̃ = [x̃1 , · · · , x̃D ]
and their corresponding training targets y = [∆x1 , · · · , ∆xD ]T , the joint distribution between y
and a test target ∆x∗k for training input x̃∗k is assumed to follow a Gaussian distribution. That is,
!
!
y
K(X̃, X̃) + σn I K(X̃, x̃∗k )
p
∼ N 0,
(2)
∆x∗k
K(x̃∗k , X̃)
K(x̃∗k , x̃∗k )
In addition, the posterior distribution over the observations can be obtained by restricting the joint
distribution to only contain those targets that agree with the observations [15]. This is achieved
by conditioning the joint distribution on the observations, and results in the predictive mean and
variance function as follows:
m(x̃∗k ) = Ef [∆x∗k ] = K(x̃∗k , X̃)K−1
σ y
σ

2

(x̃∗k )

=

Varf [∆x∗k ]

= K(x̃∗k , x̃∗k )
∗
− K(x̃∗k , X̃)K−1
σ K(X̃, x̃k )

(3a)
(3b)

where Kσ = K(X̃, X̃) + σn I. The state at the next sampling time also follows a Gaussian distribution. Thus,
p(xk+1) ∼ N (µk+1, Σk+1)
(4)

where
µk+1 = xk + m(x̃∗k )

(5a)

2

(5b)

Σk+1 = σ

(x̃∗k )

Typically, the hyperparameters of the GP model are learned by maximizing the log-likelihood
function given by
1
D
1
log K−1
− log(2π)
log p(y|X̃, θ) = − yT K−1
σ y−
σ
2
2
2

(6)

This results in a nonlinear non-convex optimization problem that is traditionally solved by using CG or BFGS algorithms. Recently, PSO based algorithms that minimizes the model error
instead of the log-likelihood function have been shown in [18] to be more efficient and effective.
2.2.

Uncertainty propagation

With the GP model obtained, one-step-ahead predictions can be made by using (3) and (5). When
multiple-step predictions are required, the conventional way is to iteratively perform multiple onestep-ahead predictions using the estimated mean values. However, this process does not take into
account the uncertainties introduced by each successive prediction. This issue has been shown to
be important in time-series predictions [25].
The uncertainty propagation problem can be dealt with by assuming that the joint distribution
of the training inputs is uncertain and follows a Gaussian distribution. That is,

p(x̃k ) = p xk , uk ∼ N (µ̃k , Σ̃k )
(7)

with mean and variance given by

µ̃k = [µk , E [uk ]]T


Σk
Cov [xk , uk ]
Σ̃k =
Cov [uk , xk ]
Var [uk ]

(8a)
(8b)

where Cov [xk , uk ] = E [xk uk ] − µk E [uk ]. Here, E [uk ] and Var [uk ] are the mean and variance of
the system controls.
The exact predictive distribution of the training target could then be obtained by integrating
over the training input distribution:
Z
∗
p(∆xk ) = p(f (x̃∗k )|x̃∗k )p(x̃∗k )dx̃∗k
(9)
However, this integral is analytically intractable. Numerical solutions can be obtained using
Monte-Carlo simulation techniques. In [26], a moment-matching based approach is proposed to
obtain an analytical Gaussian approximation. The mean and variance at an uncertain input can
be obtained through the laws of iterated expectations and conditional variances respectively [24].
They are given by
h 
i
(10a)
m(x̃∗k ) = Ex̃∗k Ef ∆x∗k
h
i
h
i




(10b)
σ 2 (x̃∗k ) = Ex̃∗k Varf ∆x∗k + Varx̃∗k Ef ∆x∗k

Equation (5) then becomes
µk+1 =µk + m(x̃∗k )
2

(x̃∗k )

Σk+1 =Σk + σ




+ Cov xk , ∆xk + Cov ∆xk , xk

(11a)
(11b)

The computational complexity of GP inference using (10) is O(D 2 n2 (n + m)) which is quite
high. Hence, GP is normally only suitable for problems with limited dimensions (under 12 as
suggested by most publications) and limited size of training data. For problems with higher dimensions, sparse GP approaches [27] are often used.
3.

GP Based Local Dynamical Models

When dealing with the control of nonlinear systems, it is common practice to obtain local linearized
models of the system around operating points. The main purpose is to reduce the computation
involved in the nonlinear control problem. The same technique is used here for the GP based MPC
optimization problem. The main difference here is that the model of the system is probabilistic
rather than deterministic. Thus there is more than one way by which the GP model could be
linearized.
In [28], a GP based local dynamical model allows standard robust control methods to be used
on the partially unknown system directly. Another GP based local dynamical model is proposed
in [29] to integrate GP model with dynamic programming. In these two cases, the nonlinear optimization problems considered are unconstrained.
In this section, we shall present two different GP based local models. They will be applied to
the constrained nonlinear problems presented in Section 4.
3.1.

Basic GP-based Local Model

Linearization can be done based on the mean values in the GP model. In this case we replace the
state vector xk by its mean µk . Then (1) becomes
µk+1 = F (µk , uk )

(12)

Let (µ∗k , u∗k ) be the operating point at which the linearized model is to be obtained. Given that
∆µk = µk − µ∗k and ∆uk = uk − u∗k are small, from (12), we have
∂F
∂F
∆µk +
∆uk
∂µk
∂uk
∂µk+1
∂µk+1
∆µk +
∆uk
=
∂µk
∂uk

∆µk+1 =

(13a)
(13b)

∂µk+1
∂µk+1
and
are the Jacobian state and input matrices respectively. Using the chain
∂µk
∂uk
rule, we get
where

∂µk+1 ∂ µ̃k ∂µk+1 ∂ Σ̃k
∂µk+1
=
+
∂µk
∂ µ̃k ∂µk
∂ Σ̃k ∂µk
∂µk+1
∂µk+1 ∂ µ̃k ∂µk+1 ∂ Σ̃k
=
+
∂uk
∂ µ̃k ∂uk
∂ Σ̃k ∂uk

(14a)
(14b)

where

∂ µ̃k ∂ Σ̃k ∂ µ̃k ∂ Σ̃k
∂µk+1
,
,
,
can be easily obtained based on (8). Elaborations of
and
∂µk ∂µk ∂uk ∂uk
∂ µ̃k

∂µk+1
can be found in [24].
∂ Σ̃k
3.2.

Extended GP-based Local Model

Model uncertainties are characterized by the variances. However, the basic local model derived
above only involves the mean values. The extended local model aims to take into account model
uncertainties. Similar to what
p we have done2to derive the basic model, we replace the state vector
xk in (1) by sk = [µk , vec( Σk )]T ∈ Rn+n which shall be known as the “extended state”. Here,
vec(·) denotes the vectorization of a matrix 1 . Hence (1) becomes
sk+1 = F ′ (sk , uk )
Linearizing at the operating point (s∗k , u∗k ) where s∗k = [µ∗k , vec(
∆sk+1

(15)
p

Σ∗k )]T , we have

∂F ′
∂F ′
∆sk +
∆uk
=
∂sk
∂uk

Here, ∆sk = sk − s∗k and ∆uk = uk − u∗k .

∂µk+1
∂F ′ 
√∂µk
=
 ∂ Σk+1
∂sk
∂µk
 ∂µ
k+1
∂F ′  √∂uk
=
∂ Σk+1
∂uk
∂uk

The Jacobian matrices are

∂µk+1
√
∂√ Σk 
 ∈ R(n+n2 )×(n+n2 )
∂ Σk+1 
√
∂ Σk


(n+n2 )×m
∈R

(16)

(17a)

(17b)

with the entries given by
∂µk+1
√
∂ Σk
√
∂ Σk+1
∂µ
√ k
∂ Σk+1
√
∂ Σk
√
∂ Σk+1
∂uk

∂µk+1 ∂Σk
√
∂Σk ∂ Σk
√
∂ Σk+1 ∂Σk+1
=
∂Σk+1 ∂µk
√
∂ Σk+1 ∂Σk+1 ∂Σk
√
=
∂Σk+1 ∂Σk ∂ Σk
√
∂ Σk+1 ∂Σk+1
=
∂Σk+1 ∂uk
=

(18a)
(18b)
(18c)
(18d)

1
Σk is a real symmetric matrix therefore can be diagonalized. The square root of a diagonal matrix can simply be obtained by computing the
square roots of diagonal entries.

√
√
1
∂ Σk+1
1
∂ Σk
, they can be expressed as
= √
and
= √
Since
∂Σk
∂Σk+1
2 Σk+1
2 Σk
∂µk+1
∂µk+1 ∂ µ̃k
∂µk+1 ∂ Σ̃k
=
+
∂Σk
∂ µ̃k ∂Σk
∂ Σ̃k ∂Σk
∂Σk+1 ∂ µ̃k ∂Σk+1 ∂ Σ̃k
∂Σk+1
=
+
∂µk
∂ µ̃k ∂µk
∂ Σ̃k ∂µk

(19a)
(19b)

∂Σk+1
∂Σk+1 ∂ µ̃k
∂Σk+1 ∂ Σ̃k
=
+
∂Σk
∂ µ̃k ∂Σk
∂ Σ̃k ∂Σk
∂Σk+1 ∂ µ̃k ∂Σk+1 ∂ Σ̃k
∂Σk+1
=
+
∂uk
∂ µ̃k ∂uk
∂ Σ̃k ∂uk

(19c)
(19d)

∂ µ̃k
∂ Σ̃k
∂Σk+1
∂Σk+1
and
can be easily obtained based on (8). Elaborations of
and
can be
∂Σk
∂Σk
∂ µ̃k
∂ Σ̃k
found in [24].
4.

Model Predictive Control based on GP

A discrete-time nonlinear dynamical system defined by (1) is required to track a trajectory given by
{rk } for k = 1, 2, · · · . Using MPC with a prediction horizon H ≥ 1, the optimal control sequence
can be obtained by solving the following problem:
Vk∗ = min J (xk , uk−1 , rk )

(20a)

s.t. xk+i|k = f (xk+i−1|k , uk+i−1 )
xmin ≤ xk+i|k ≤ xmax
umin ≤ uk+i−1 ≤ umax
i = 1, · · · , H

(20b)
(20c)
(20d)

u(·)

where only the first control action uk of the resulting control sequence u(·) = [uk , · · · , uk+H−1 ]T
is applied to the system at time k. xmin ≤ xmax and umin ≤ umax are the upper and lower bounds of
the system states and control inputs, respectively.
In the rest of this paper, the cost function J (xk , uk−1 , rk ) shall be rewritten as J (xk , uk−1) for
brevity. The quadratic cost function given by
J (xk , uk−1) =

H n
X
i=1

xk+i − rk+i

2
Q

+ uk+i−1

2
R

o

(21)

will be used. Here, · Q and · R denote the two 2-norms weighted by positive definite matrices
Q ∈ Rn×n and R ∈ Rm×m respectively. The control horizon will be assumed to be equal to the
prediction horizon.

4.1.

GPMPC1

4.1.1. Problem Reformulation: We assume that the system function f (·) is unknown and it is
replaced by a GP model. Consequently, problem (20) becomes a stochastic one [30]:


Vk∗ = min E J (xk , uk−1 )
(22a)
u(·)

s.t. p(xk+1 |xk ) ∼ N (µk+1, Σk+1 )
umin ≤ uk+i−1 ≤ umax

p xk+i|k ≥ xmin ≥ η

p xk+i|k ≤ xmax ≥ η

(22d)

µk+i − 2Σk+i ≥ xmin
µk+i + 2Σk+i ≤ xmax

(23a)
(23b)

(22b)
(22c)
(22e)

where η denotes a confidence level. For η = 0.95, the chance constraints (22d) and (22e) are
equivalent to

Using (21) as the cost function, we get


E J (xk , uk−1 )
H n
hX
=E
xk+i − rk+i

2
Q

+ uk+i−1

i=1

H
h
X
=
E xk+i − rk+i

=

i=1
H 
X

E

i=1

h

xk+i −

2
Q

+ uk+i−1

2
rk+i Q

i

+E

h

2
R

2
}
R

i

i

2
uk+i−1 R

(24)
i

 
In practice, the controls are deterministic. Hence, E u2k = u2k and (24) becomes


E J (xk , uk−1 )

H n 
o
X
2
2
E xk+i − rk+i Q + uk+i−1 R
=
i=1

=

H n
X
i=1

µk+i − rk+i

=h (µk , uk−1 )

2
Q

+ uk+i−1

2
R

+ trace QΣk+i

(25)

o

Now we have a deterministic cost function which involve the model variance Σ that allows
model uncertainties to be explicitly included in the computation of the optimized controls. Note
that for multiple-step predictions with uncertainty propagation, the computational complexity of
problem (22) will not increase even though the GP model becomes more complicated. This is
because the modelling and the control processes are independent of each other.

4.1.2. Nonlinear Optimization Solution: With the cost function (25) and the state constraints
(23), the original stochastic optimization problem (22) has been relaxed to a deterministic constrained nonlinear optimization problem. But it is typically non-convex. This is usually solved
by derivative-based approaches, such as Lagrange multipliers [31] based on first-order derivatives
(gradient), or by SQP and interior-point algorithms based on second-order derivatives (Hessians
matrix) [32]. When the derivative of the cost function is unavailable or is too difficult to compute,
it could be approximated iteratively by a sampling method [33, 34]. Alternatively, evolutionary
algorithms, such as PSO [35] and GA [36], could be used to solve the problem. This approach is
able to handle general constrained optimization problems. However, there is no guarantee that the
solutions obtained are near optimum. A review of nonlinear optimization techniques for the MPC
problem can be found in [32].
A suitable technique for solving our MPC problem is the FP-SQP algorithm proposed in [37]. It
can be explained using the following general form of a constrained nonlinear optimization problem:
min h(z)
z

z ∈ Rn+m

s.t. c(z) = 0
d(z) ≤ 0

(26a)
(26b)
(26c)

where h : Rn+m → R is the objective function, c : Rn+m → Rn and d : Rn+m → Rn+m represents
the corresponding equality and inequality constraints, respectively. FP-SQP generates a sequence
of feasible solutions {zj }j=0,1,2,··· by splitting the original problem into several QP sub-problems.
In particular, a step ∆zj from current iterate zj to the next one zj+1 can be obtained by solving the
following QP subproblem:
1
T
min▽h(zj )T ∆zj + ∆zj Hj ∆zj
∆zj
2
j
j T
s.t. c(z ) + ▽c(z ) ∆zj = 0
d(zj ) + ▽d(zj )T ∆zj ≤ 0

(27a)
(27b)
(27c)

under the trust-region constraint
∆zj ≤ γ j

(28)

where ▽h(·) denotes the first-order derivative of the objective function at zj , ▽c(·) and ▽d(·)
are two linearised Jacobian matrices at the zj . The matrix Hj ∈ R(n+m)×(n+m) is an exact or
approximated Lagrangian Hessian matrix and γ j represents the trust-region radius. To guarantee
the feasibility of ∆zj , its corresponding perturbation ∆z̃j which satisfies the following conditions
need to be computed:
zj + ∆z̃j ∈ Π
1
˜
k∆zk2 ≤ ∆z
2

(29a)
2

≤

3
k∆zk2
2

(29b)

where Π denotes the feasible points set of problem (26). A method to obtain such a perturbation
is proposed in [38]. An acceptability value of ∆zj defined by:
ρj =

h(zj+1) − h(zj )
−▽h(zj )T ∆zj − 12 ∆zj T Hj ∆zj

(30)

1

Initialization
feasible point z0 ∈ Π,
Hessian matrix H0 ,
trust region upper bound γmax > 0,
initial trust region radius γ 0 = k▽h(z0 )k
τ = 0, 0 < τ1 < τ2 < 1

2

for j = 0, 1, 2, · · · , J < ∞ do
Obtain step ∆zj by solving the problem (27);
1
T
if ▽h(zj )T ∆zj + ∆zj Hj ∆zj = 0 then
2
Stop;
else
Update ρj by using (30);
Update zj+1 , τ :
n j
z + ∆zj , ρj ≥ τ1
zj+1 =
zj ,
otherwise

3
4
5
6
7
8

k∆zj k
, ρj ≥ τ1
τ = { k▽h(zj+1) − h(zj )k
τ /4,
otherwise

9

Update trust region radius:
γ

j+1


min τ k▽h(zj+1 )k, γmax , ρj ≥ τ2
={
τ k▽h(zj+1 )k,
otherwise

Update Hessian matrix Hj+1 by using (35);
11
j = j + 1;
12
end
13 end
Algorithm 1: The Feasibility-Perturbed Sequential Quadratic Programming used in the
GPMPC1 algorithm
10

If this value is not acceptable, then the trust-region radius γ j will need to be adjusted. An adaptive method to adjust γ j can be found in [39]. The complete FP-SQP algorithm is described in
Algorithm 1.
4.1.3. Application to GPMPC1: Applying FP-SQP to the GPMPC1 problem (22), it should be
noted that the constraints (23) are linear. Therefore it is possible to simply use ∆z̃j = ∆zj . The
next iterate then can be obtained by
zj+1 = zj + ∆zj

(31)

Expressing the cost function (25) as (26a), define zk = [µTk , uTk−1 ]T ∈ Rn+m . Hence,
h(zk ) =

H n
X
i=1

zTk+i



Q 0
0 R



o

zk+i + trace QΣk+i

(32)

One key issue in using FP-SQP is the local linearisation at ∆zj . The basic GP based local
model derived Section 3.1 shall be used to derive the QP subproblem as:
H n
X
∂h
1
min
∆zk+i + ∆zTk+i Hk ∆zk+i
∆zk ,∆Σk
∂zk+i
2
i=1
o

+ trace Q(Σk+i + ∆Σk+i )

s.t. ∆µk+i+1 = Ak+i ∆µk+i + Bk+i ∆uk+i
umin ≤ uk+i + ∆uk+i ≤ umax
µk+i + ∆µk+i − 2(Σk+i + ∆Σk+i ) ≥ xmin
µk+i + ∆µk+i + 2(Σk+i + ∆Σk+i ) ≤ xmax
k∆zk+i k ≤ γ

(33a)
(33b)
(33c)
(33d)
(33e)
(33f)
(33g)

∂µk+i+1
∂µk+i+1
and Bk+i =
are the two Jacobian matrices of the basic GP
∂µk+i
∂uk+i
based local model (12).
The computation of the Hessian matrix Hk of the Lagrangian in (27) is another key issue when
using the FP-SQP algorithm. The exact Hessian matrix is usually obtained by
Note that Ak+i =

2

Hk = ▽ h(zk ) +

n
X

2

αi ▽ c(zk ) +

i=1

n+m
X

βi ▽2 d(zk )

(34)

i=1

where α and β are two Lagrange multipliers applied to the equality and the inequality constraints
respectively. This allows rapid local convergence but requires the second-order derivatives ▽2 c(zk )
which are generally not available. When the system is represented by a GP model, these derivatives
are mathematically computable 2 but are computationally expensive to obtain. In addition, the exact
Hessian matrix may be not positive definite. To address these issues, approximation approaches
have been proposed in [40]. In our work, Hk is approximately updated by using a Quasi-Newton
method based on the BFGS. The update equation is given by
Hk+1 = Hk −

Hk ∆zk ∆zTk Hk
yk ykT
+
∆zTk Hk ∆zk
ykT ∆zk

(35)

where ∆zk = zk+1 − zk and yk = µk+1 − µk .
4.2.

GPMPC2

With GPMPC1, model uncertainty was introduced through the variance term into the objective
function in (25). But this is an indirect way to handle model uncertainties. A more direct approach is to introduce the variance into that state variable. This can be done through the use of
the extended GP based local model (16). In this way, the variances are directly handled in the
optimization process.
Another disadvantage of GPMPC1 is that the MPC optimization problem (22) is non-convex.
Due to the recursive nature of SQP optimizations, the process could be time consuming. With
2

As shown in [24], the first-order derivatives are functions of µ̃ and Σ̃, the second-order derivatives therefore can be obtained by using the
chain-rule.

GPMPC2, the non-convex problem is relaxed to a convex one, making it much easier to solve.
Sensitivity to initial conditions is reduced and in most cases exact solutions can be obtained [32].
This convex optimization problem can be solved offline by using mp-QP [41] where the explicit
solutions are computed as a lookup table of nonlinear controllers. An example can be found in [20].
However, the size of the table grows exponentially with the number of states. Hence it is only
suitable for problems with less than 5 states [42]. Using the extended GP based local model (16),
the problem can be solved efficiently by an online active-set algorithm.
4.2.1. Problem Reformulation:
state variable as

Based on the extended local model in Section 3.2, define the


T
2
sk+1|k , · · · , sk+H|k ∈ RH(n+n )
p
p
= [µk+1 , Σk+1 , · · · , µk+H , Σk+H ]T

Zk+1 =
Also, let

Uk = [uk , · · · , uk+H−1]T ∈ RHm

r∗k+1

T

= [rk+1, 0, · · · , rk+H , 0] ∈ R

(36)
(37)

H(n+n2 )

(38)

Problem (22) then becomes
min
U

n

Zk+1 − r∗k+1

2
Q̃

+ kUk+1 k2R̃

s.t. IHn xmin ≤ Mz Zk+1 ≤ IHn xmax
IHm umin ≤ Uk+1 ≤ IHm umax

o

(39a)
(39b)
(39c)

where
Q̃ =diag{[Q, diag{vec(Q)}, · · · , Q, diag{vec(Q)}]} ∈ RH(n+n

a

R̃ =diag{[R, · · · , R]} ∈ R

Hm×Hm

2 )×H(n+n2 )

,

Ia ∈ R is the identity vector, and
 T
0 ··· 0
In 2ITn2 0
T
 0
0 In 2ITn2 · · · 0

Mz =  ..
..
..
..
..
..
 .
.
.
.
.
.
T
0
0
0 · · · In 2ITn2




2

 ∈ RH×H(n+n )


,

(40a)
(40b)

(41)

Let Tu ∈ RHm×Hm be a lower triangular matrices with unit entries. Then,
Uk = IHm uk−1 + Tu ∆Uk

(42)

∆Zk+1 can be expressed as
∆Zk+1 = Ã∆sk + B̃∆Uk
based on the extended local model, with the state and control matrices given by

T
2
Ã = A, A2, · · · , AH ∈ RH(n+n )


B
0
··· 0
B
··· 0 
 AB
2
∈ RH(n+n )×Hm
B̃ = 
.
.
..
.. 


..
..
.
.
AH−1 B AH−2 B · · · B

(43)
(44a)

(44b)

where A and B are the two Jacobian matrices (17) and (18) respectively. The corresponding state
variable Zk+1 is therefore given by


Zk+1 = sk + Tz Ã∆sk + B̃∆Uk
(45)
2

2

where Tz ∈ RH(n+n )×H(n+n ) denotes a lower triangular matrix with unit entries.
Based on (42) and (45), problem (39) can be expressed in a more compact form as
1
min k∆Uk k2Φ + ψ T ∆Uk + C
∆U 2


Tu
s.t. ∆Umin ≤
∆Uk ≤ ∆Umax
Tz B̃

(46a)
(46b)

where
Φ =B̃T TTz Q̃Tz B̃ + TTu R̃Tu ∈ RHm×Hm
T

ψ =2(sk Q̃Tz B̃ + ∆sk Ã Q̃B̃

C

− r∗k+1 Q̃Tz B̃ + uk−1 R̃Tu ) ∈ RHm
=(s2k + r∗k+1 )Q̃ + 2sk ∆sk Q̃Tz Ã
+ u2k−1 R̃ + ∆s2k ÃT Q̃Ã
− 2r∗k+1 (sk Q̃ − ∆sk Q̃Tz Ã)



IHm (umin − uk−1 )
IH(n+n2 ) (xmin − sk − Tz Ã∆sk )


IHm (umax − uk−1 )
=
IH(n+n2 ) (xmax − sk − Tz Ã∆sk )

∆Umin =
∆Umax



(47a)
(47b)
(47c)
(47d)
(47e)
(47f)
(47g)
(47h)

Since Q̃, R̃, Tz and Tu are positive definite, Φ is also positive definite. Hence (46) is a constrained QP problem and is strictly convex. The solution will therefore be unique and satisfies
the KKT conditions.
4.2.2. Quadratic Programming Solution: The optimization problem (46) can be solved by an
active-set method [43]. It iteratively seeks an active (or working) set of constraints and solve an
equality constrained QP problem until the optimal solution is found. The advantage of this method
is that accurate solutions can still be obtained even when they are ill-conditioned or degenerated.
In addition, it is conceptually simple and easy to implement. A warm-start technique could also be
used to accelerate the optimization process substantially.
Let G = [Tu , Tz B̃]T , the constraint (46b) can be written as




G
∆Umax
∆U ≤
(48)
−G
−∆Umin
Ignoring the constant term C, problem (46) becomes
1
min k∆Uk k2Φ + ψ T ∆Uk
∆U 2
˜U
s.t. G̃∆Uk ≤ ∆

(49a)
(49b)

˜ U = [∆Umax , −∆Umin ]T ∈ R2H(m+n+n ) .
where G̃ = [G, −G]T ∈ R2H(m+n+n )×Hm and ∆
Let Π∆U be the set of feasible points, and I = {1, · · · , 2H(m + n + n2 )} be the constraint
index set. For a feasible point ∆U∗k ∈ Π∆U , the index set for the active set of constraints is defined
as
˜ U,i }
A(∆U∗k ) = {i ⊆ I|G̃i ∆U∗k = ∆
(50)
˜ U,i is the ith row of the ∆
˜ U . The inactive set is therefore given
where G̃i is the ith row of G̃ and ∆
by
B(∆U∗k ) = I \ A(∆U∗k )
(51)
˜ U,i }
= {i ⊆ I|G̃i ∆U∗k < ∆
2

2

Given any iteration j, the working set Wkj contains all the equality constraints plus the inequality
constraints in the active set. The following QP problem subject to the equality constraints w.r.t.
Wkj is considered given the feasible points ∆Ujk ∈ Π∆U :
1
2
∆Ujk + δ j Φ + ψ T (∆Ujk + δ j )
δj 2
1 j 2
δ Φ + (ψ + Φ∆Ujk )T δ j
= min
j
2
δ
1
2
∆Ujk Φ + ψ T ∆Ujk
+
|2
{z
}
min

s.t.

constant
G̃i (∆Ujk + δ j )

(52a)
(52b)

˜ U,i , i ∈ W j
=∆
k

(52c)

This problem can be simplified by ignoring the constant term to:
1 j 2
δ Φ + (ψ + Φ∆Ujk )T δ j
2
1 jT
δ Φδ j + (ψ + Φ∆Ujk )T δ j
= min
δj 2
˜ U,i − G̃i ∆Uj , i ∈ W j
s.t. G̃i δ j = ∆
k
k
min

(53a)

δj

(53b)
(53c)

By applying the KKT conditions to problem (53), we can obtain the following linear equations:

 j  

Φ G̃TA
−ψ − Φ∆Ujk
δ
(54)
= ˜
λ∗k
G̃A 0
∆U,A − G̃A ∆Ujk
|
{z
}
Lagrangian Matrix

2
˜ U,A ⊂ ∆
˜U
where λ∗k ∈ R2H(m+n+n ) denotes the vector of Lagrangian multipliers, G̃A ⊆ G̃ and ∆
j
are the weighting matrix and the upper bounds of the constraints w.r.t. Wk . Let the inverse of
Lagrangian matrix be denoted by


−1 
Φ G̃TA
L1 LT2
(55)
=
L2 L3
G̃A 0

If this inverse exists, then the solution is given by
˜ U,A − G̃A ∆Uj )
δ j = −L1 (ψ + Φ∆Ujk ) + LT2 (∆
k
j
j
∗
˜ U,A − G̃A ∆U )
λ = −L2 (ψ + Φ∆U ) + L3 (∆
k

k

k

(56a)
(56b)

where
L1 = Φ−1 − Φ−1 G̃TA (G̃A Φ−1 G̃TA )−1 G̃A Φ−1

(57a)

L2 = Φ−1 G̃TA (G̃A Φ−1

(57b)

L3 = −(G̃A Φ−1 G̃TA )−1

(57c)

If δ j =
6 0, then the set of feasible points ∆Ujk fails to minimize problem (49). In this case, the
next set of feasible point is computed for the next iteration by ∆Uj+1
= ∆Ujk + κj δ j with step
k
size
(
)
˜ U,i − G̃i ∆Uj
∆
k
κj = min 1, min
(58)
G̃i δ j
i∈B(∆Ujk )
˜ U,i − G̃i ∆Uj
∆
k
should be “activated”,
j
j
G̃i δ
i∈B(∆Uk )

If κj < 1, the inequality constraint with index q = argmin

giving the working set Wkj+1 = Wkj ∪ q. Otherwise, we have Wkj+1 = Wkj .
Alternatively, if the solution gives δ j = 0, then the current feasible points ∆Ujk could be the
optimal solution. This can be verified by checking the Lagrangian multiplier λ∗k = min λ∗k,i . If
i∈Wkj ∩I

λ∗k

≥ 0, the optimal solution of the (49) at sampling time k is found. Otherwise, this inequality
constraint indexed by p = argmin λ∗k,i should be removed from the current working set, giving us
i∈Wkj ∩I

Wkj+1 = Wkj \ p. Algorithm 2 summarizes the active set algorithm used in the GPMPC2.
4.2.3. Implementation Issues: The key to solving the linear equations (54) is the inverse of
the Lagrangian matrix. However, G̃A is not always full ranked. Thus the Lagrangian matrix is
not always invertible. This problem can be solved by decomposing G̃A using QR factorization
technique, giving us GTA = Q [R 0]T where R ∈ Rm1 ×m1 is an upper triangular matrix with
m1 = rank(G̃A ). Q ∈ RHm×Hm is an orthogonal matrix that can be further decomposed to
Q = [Q1 Q2 ] where Q1 ∈ RHm×m1 and Q2 ∈ RHm×(Hm−m1 ) . Thus, GTA = Q = Q1 R and
L1 = Q2 (QT2 ΦQ2 )−1 QT2
T

L2 = Q1 R−1 − L1 ΦQ1 R−1
L3 = R−1 QT1 ΦL2

(59a)
T

(59b)
(59c)

The second issue relates to using the appropriate warm-start technique to improve the convergence rate of the active-set method. For GPMPC2, since the changes in the state between two
successive sampling instants are usually quite small, we simply use the previous ∆U∗k as the starting point ∆U0k+1 for the next sampling time k + 1. This warm-start technique is usually employed
in MPC optimizations because of its proven effectiveness [42].
4.3.

Stability

The stability of the closed-loop controller is not guaranteed because the MPC problem is openloop. This can be demonstrated by the stability analysis of the proposed algorithms.

1

2
3
4
5

Initialization
the feasible point ∆U0k ∈ Π∆U ;
the working set W 0 = A(∆U0k );

for j = 0, 1, 2, · · · do
Compute the δ j and λ∗k by solving the linear equations (54);
if δ j = 0 then
λ∗k = min λ∗k,i ,
i∈Wkj ∩I

6
7
8
9
10
11
12
13
14
15
16

17
18
19
20
21
22
23
24
25

p = argmin λ∗k,i
λ∗k

i∈Wkj ∩I

≥ 0 then
∆U∗k = ∆Ujk ;
Stop.
else
Wkj+1 = Wkj \ p;
∆Uj+1
= ∆Ujk ;
k
end
else
Compute the step length κj by (58),
˜ U,i − G̃i ∆Uj
∆
k
q = argmin
j
G̃i δ
i∈B(∆Ujk )
j
if κ < 1 then
∆Uj+1
= ∆Ujk + κj δ j ;
k
j
A(∆Uj+1
k ) = A(∆Uk ) ∪ q;
else
∆Uj+1
= ∆Ujk + δ j ;
k
j
A(∆Uj+1
k ) = A(∆Uk );
end
end
end
Algorithm 2: Active set method for solving the GPMPC2 problem
if

In particular, for the MPC problem (22) in the GPMPC1 algorithm, the objective (25) can be
directly used as the Lyapunov function. Therefore, it can be known that
∗

V (k) =

H n
X
i=1

∆µ∗k+i

2
Q

+ u∗k+i−1

2
R

+ trace(QΣ∗k+i )

o

(60)

where ∆µ∗k+i = µ∗k+i − rk+i , u∗ is the optimal control inputs, and µ∗k+i and Σ∗k+i represent the
corresponding optimal means and variances of the GP model at time k. The Lyapunov function at

time k + 1 is subsequently obtained by,
V(k + 1)
H n
o
X
=
k∆µk+1+i k2Q + kuk+i k2R + trace(QΣk+1+i )

(61a)
(61b)

i=1

= V ∗ (k) − ∆µ∗k+1

2
Q

− ku∗k k2R − trace(QΣ∗k+1 )

(61c)

+ k∆µk+1+H k2Q + kuk+H k2R + trace(QΣk+1+H )

It is easy to know that V ∗ (k + 1) ≤ V(k + 1) due to the nature of the optimization. Furthermore,
the following inequality can be obtained,
V ∗ (k + 1) ≤V(k + 1)

≤V (k) + k∆µk+1+H k2Q
+ kuk+H k2R + trace(QΣk+1+H )
∗

(62a)
(62b)

2

because of ∆µ∗k+1 Q ≥ 0, ku∗k k2R ≥ 0 and trace(QΣ∗k+1) ≥ 0. The stability result of the problem
(46) in the GPMPC2 algorithm can be obtained in the same way.
The result in (62) shows that, to guarantee the stability, additional terminal constraints on the
means and variances of the GP model, as well as the control inputs are required such that,
µk+H+1|k − rk+H+1 = 0
Σk+H+1|k = 0
uk+H|k = 0

(63a)
(63b)
(63c)

However, it should be noted that, these newly added constraints altered the optimization problem.
Hence its feasibility will need to be analysed. Another approach to provide the guaranteed stability
is by introducing a terminal cost into the objective function [1].
5.

Numerical Simulations

GPMPC1 and GPMPC2 are applied to two trajectory tracking problems of a MIMO nonlinear system with time-varying parameters. For each problem, 50 independent simulations are performed
on a computer with a 3.40GHz Intelr CoreTM 2 Duo CPU with 16 GB RAM, using Matlabr
version 8.1. The average simulation results of these 50 trials are presented here.

−3

8

x 10

0.15
0.1
1

Y (k)

4

1

Y (k)

6

2

0

0
−2

0.05

0

20

40

60

80

100

120

140

160

−0.05

180

0

20

40

60

80

100

120

140

160

180

0

20

40

60

80
100
120
sample time k

140

160

180

−3

5

x 10

0.1
0.05

2

2

Y (k)

Y (k)

0

−5

−10

0
−0.05

0

20

40

60

80
100
120
sample time k

140

160

180

−0.1

(a) “Step” trajectory

(b) “Lorenz” trajectory

Fig. 1. Training errors of the system outputs for the two trajectory tracking problems.
The MIMO nonlinear system in [44] is used for our simulations. It is described by:
x1 (k)2
+ 0.3x2 (k)
1 + x1 (k)2
x1 (k)2
x2 (k + 1) =
1 + x2 (k)2 + x3 (k)2 + x4 (k)2
+ a(k)u1 (k)
x3 (k)2
x3 (k + 1) =
+ 0.2x4 (k)
1 + x3 (k)2
x3 (k)2
x4 (k + 1) =
1 + x1 (k)2 + x2 (k)2 + x4 (k)2
+ b(k)u2 (k)
y1 (k + 1) = x1 (k + 1) + ω1
y2 (k + 1) = x3 (k + 1) + ω2

x1 (k + 1) =

(64)

where x1 , x2 , x3 and x4 are system states, u1 , u2 and y1 , y2 denote system inputs and outputs,
respectively. ω1 , ω2 ∼ N (0, 0.01) are independent Gaussian white noise. In addition, the timevarying parameters a(k) and b(k) are given by
a(k) = 10 + 0.5 sin(k)
10
b(k) =
1 + exp(−0.05k)
5.1.

(65a)
(65b)

“Step” Trajectory Tracking

The objective of the first experiment is to steer the nonlinear system to follow a step trajectory
shown as the reference in Figure 2a. The system inputs are subjected to the following constraints:
0 ≤ u1 (k) ≤ 5,

0 ≤ u2 (k) ≤ 5

0.4

Reference
GPMPC1
GPMPC2

0.2
0

0

20

40

60

80

100

120

140

160

1
GPMPC1
GPMPC2

1

1.5

IAE on Y

2

0.6
u1(k)

Y1(k)

0.8

1
0.5
0

180

1

20

40

60

80

100

120

140

160

0.2

180

3

1

0

20

40

60

80
100
120
sample time k

140

160

0

180

80

100

120

140

160

180

1
0.8

0

20

40

60

80
100
120
sample time k

140

160

0.4

180

GPMPC1
GPMPC2

40

60

80

100

120

140

160

−4

180

10

10

5
u (k)

20

2

0

−5

−20

−10

0

20

40

60

80
100
120
sample time k

140

160

180

(d) Controlled Outputs – “Lorenz”

80
100
120
sample time k

140

160

180

IAE on Y1

60
40

GPMPC1
GPMPC2

20
0

20

40

60

80

100

120

140

160

0

180

0

20

40

60

80

100

120

140

160

180

0

20

40

60

80
100
120
sample time k

140

160

180

150

0

−10

60

80

0
−2

20

40

100

IAE on Y2

−10

0

20

(c) IAE – “Step”

GPMPC1
GPMPC2

2
u1(k)

0

0

(b) Control Inputs – “Step”
4

Reference
GPMPC1
GPMPC2

10
1

60

0.6

20

Y (k)

40

2

IAE on Y

u2(k)

Y2(k)

0.4

(a) Controlled Outputs – “Step”

2

20

1.2

0.2

Y (k)

0

1.4

2

0.6

−20

0.6
0.4

0

0.8

0

0.8

0

20

40

60

80
100
120
sample time k

140

160

180

100

50

0

(e) Control Inputs – “Lorenz”

(f) IAE – “Lorenz”

Fig. 2. Simulation results of the two trajectory tracking problems.
To generate the observations for GP modelling, this problem is first solved by using the NMPC
strategy proposed in [45]. 189 observations are collected and are used to train the GP models. The
learning process took approximately 2.1 seconds, with a training MSE of 9.9114 × 10−5. Figure 1a
shows the training errors for the 189 samples. These results show that the system is accurately
learnt by using the GP models.
The MPC parameters in this simulation are: initial states x0 = [0, 0, 0, 0]T and initial control
inputs u0 = [0, 0]T , weighting matrix Q = I4×4 and R = I2×2 . In addition, the prediction horizon
H is 10. Theoretically, a long enough H is necessary to guarantee the stability of MPC controllers.
However, the complexity of MPC problem increases exponentially with increasing H. This value
of H is chosen as a trade-off between the control performance and computational complexity.
The resultant controlled outputs and control inputs by using GPMPC1 and GPMPC2 are shown
in Figures 2a and 2b, respectively. They show that both algorithms exhibit equally good control
performances in this task since they both produced outputs close to the target. The IAE values can
be found in Figure 2c.
GPMPC1 takes on average 34.1 seconds to compute the 189 optimized control inputs. However,
GPMPC2 only requires 4.51 seconds which is more than 8 times more efficient than GPMPC1.
This shows the advantage in our formulation of the problem as convex optimization.
5.2.

“Lorenz” Trajectory Tracking

The second problem is to track a “Lorenz” trajectory as shown in Figure 2d. In this case, the
constraints on the control inputs are:
− 4 ≤ u1 (k) ≤ 4,

−7 ≤ u2 (k) ≤ 7

Similar to the previous experiment, the NMPC method is used to generate 189 observations for
training the GP model. Training time is approximately 2.4 seconds with a training MSE of 0.0196.

15

20
Reference
Nonlinear GPMPC
GPMPC1
GPMPC2

10

Reference
Nonlinear GPMPC
GPMPC1
GPMPC2

15

10
5

Y (k)

2

1

Y (k)

5
0

0
−5
−5
−10

−15

−10

0

20

40

60

80
100
120
sample time k

140

160

180

−15

0

(a) Y1 – “Lorenz”

20

40

60

80
100
120
sample time k

140

160

180

(b) Y2 – “Lorenz”

Fig. 3. Comparison of tracking performance of GPMPC1, GPMPC2 and nonlinear GPMPC in
[46] for the Lorenz trajectory.
Table 1 MSE values for Lorenz trajectory tracking problem with GPMPC1 and GPMPC2 using GP models with
different amount of training data.

Training Data
Y1
Y2

Model for GPMPC1
60%
80%
100%
4.7831
1.36
0.0528
10.7518
1.0960
0.2995

Model for GPMPC2
60%
80%
100%
4.6493
0.6879
0.0539
6.6748
2.1522
0.3085

Figure 1b shows the training error.
The MPC parameters are: initial states x0 = [0, 0, 0, 0]T , initial control inputs u0 = [0, 0]T ,
prediction horizon H = 10, weighting matrix Q = I4×4 and R = diag{[21000, 27000]}.
The tracking results can be found in Figures 2d, 2e and 2f. They demonstrate again that the
control performance GPMPC1 and GPMPC2 are virtually the same. In this case, on average
GPMPC2 is about 5 times more efficient than GPMPC1 (5.38 seconds versus 24.72 seconds).
The performance of the two proposed algorithms is compared with the nonlinear GPMPC proposed in [46]. Even though problem (22) with cost function (25) is more complicated than the
one considered in [46], they are essentially similar. Tracking results for H = 1 are shown in Figure 5.2. show that the both two proposed algorithms outperform than the nonlinear GPMPC in
the “Lorenz” trajectory tracking problem. In addition, the GPMPC1 and GPMPC2 only require
approximately 5 and 7 seconds to compute all 189 control actions, compared to 150 seconds used
in nonlinear GPMPC.
5.3.

Sensitivity to Training Data

Since the closed-loop stability of proposed GPMPC1 and GPMPC2 are not guaranteed as discussed in Section 4.3, it is necessary to test them with different models. Here, both GPMPC1
and GPMPC2 are each tested with three separate GP models for the Lorenz trajectory tracking
problem. These models are trained by using 60%, 80% and 100% of all of 179 observations respectively. Figure 4 shows how well each model track the reference outputs. Table 1 shows the
tracking MSE values. These results indicate that while the models trained with 100% and 80%
observations perform quite well, the ones trained with 60% data are inadequate.

15

20
Reference
GPMPC1−60%
GPMPC1−80%
GPMPC1−100%
GPMPC2−60%
GPMPC2−80%
GPMPC2−100%

10

5

Reference
GPMPC1−60%
GPMPC1−80%
GPMPC1−100%
GPMPC2−60%
GPMPC2−80%
GPMPC2−100%

15

10

Y2(k)

1

Y (k)

5
0

0
−5
−5
−10

−15

−10

0

20

40

60

80
100
sample time k

(a) Y1 – “Lorenz”

120

140

160

180

−15

0

20

40

60

80
100
sample time k

120

140

160

180

(b) Y2 – “Lorenz”

Fig. 4. Comparisons between proposed approaches with different learnt GP models in the
“Lorenz” trajectory tracking problem.
6.

Conclusions

Two GP based MPC approaches (GPMPC1 and GPMPC2) have been presented for the trajectory
tracking problem of an unknown nonlinear dynamical system. The system is modelled using GP
techniques offline. These two approaches handle the model uncertainties in the form of GP variances in different ways. GPMPC1 formulated the MPC optimization problem in such a way that
model uncertainties are treated as the slack variables of GP mean constraints and are included
in the objective function as the penalty term. The resulting SMPC problem is relaxed to a deterministic non-convex nonlinear optimization problem. The solution of the resultant problem is
obtained using the FP-SQP method based on a linearized GP local model. With GPMPC2, the
variance forms part of the state vector. This allows model uncertainties to be directly included
in the computation of the optimized controls. By using the extended linearized GP local model,
the non-convex optimization problem is relaxed to a convex one which is solved using an activeset method. Simulation results on two different trajectories show that both approaches perform
equally well. However, GPMPC2 is several times more efficient computationally compared with
GPMPC1, especially for a longer horizon. A brief discussion on how closed-loop stability could
be guaranteed reveals that the resulting optimization problem will be different from the one considered in this paper. This issue will be addressed in future work.
7.

References

[1] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive
control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000.
[2] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,”
Control Engineering Practice, vol. 11, no. 7, pp. 733–764, 2003.
[3] D. Q. Mayne, “Model predictive control: Recent developments and future promise,” Automatica, vol. 50, no. 12, pp. 2967–2986, 2014.

[4] D. P. Solomatine and A. Ostfeld, “Data-driven modelling: some past experiences and new
approaches,” Journal of hydroinformatics, vol. 10, no. 1, pp. 3–22, 2008.
[5] O. Nelles, Nonlinear system identification: from classical approaches to neural networks and
fuzzy models. Springer Science & Business Media, 2013.
[6] T. Alamo, D. M. de La Peña, D. Limón, and E. F. Camacho, “Constrained min-max predictive control: Modifications of the objective function leading to polynomial complexity,” IEEE
Transactions on Automatic Control, vol. 50, no. 5, pp. 710–714, 2005.
[7] D. Limón, T. Alamo, F. Salas, and E. F. Camacho, “Input to state stability of min–max MPC
controllers for nonlinear systems with bounded uncertainties,” Automatica, vol. 42, no. 5, pp.
797–803, 2006.
[8] W. Langson, I. Chryssochoos, S. Raković, and D. Q. Mayne, “Robust model predictive control
using tubes,” Automatica, vol. 40, no. 1, pp. 125–133, 2004.
[9] L. Zhang, S. Zhuang, and R. D. Braatz, “Switched model predictive control of switched linear
systems: Feasibility, stability and robustness,” Automatica, vol. 67, pp. 8–21, 2016.
[10] A. T. Schwarm and M. Nikolaou, “Chance-constrained model predictive control,” American
Institute of Chemical Engineers, vol. 45, no. 8, pp. 1743–1752, 1999.
[11] D. Bernardini and A. Bemporad, “Scenario-based model predictive control of stochastic constrained linear systems,” in IEEE Proceedings of International Conference on Decision and
Control. IEEE, 2009, pp. 6333–6338.
[12] M. Cannon, B. Kouvaritakis, S. V. Raković, and Q. Cheng, “Stochastic tubes in model predictive control with probabilistic constraints,” IEEE Transactions on Automatic Control, vol. 56,
no. 1, pp. 194–200, 2011.
[13] A. Mesbah, S. Streif, R. Findeisen, and R. Braatz, “Stochastic nonlinear model predictive
control with probabilistic constraints,” in American Control Conference. IEEE, 2014, pp.
2413–2419.
[14] L. Fagiano and M. Khammash, “Nonlinear stochastic model predictive control via regularized
polynomial chaos expansions,” in IEEE Proceedings of International Conference on Decision
and Control. IEEE, 2012, pp. 142–147.
[15] C. Rasmussen and C. Williams, Gaussian Processes for Machine Learning.
MA, USA: MIT Press, 1 2006.

Cambridge,

[16] F. Zhu, C. Xu, and G. Dui, “Particle swarm hybridize with Gaussian process regression for
displacement prediction,” in IEEE Proceedings of International Conference on Bio-Inspired
Computing: Theories and Applications. IEEE, 2010, pp. 522–525.
[17] D. Petelin and J. Kocijan, “Control system with evolving Gaussian process models,” in IEEE
Workshop on Evolving and Adaptive Intelligent Systems (EAIS). IEEE, 2011, pp. 178–184.
[18] G. Cao, E. M.-K. Lai, and F. Alam, “Particle swarm optimization for convolved Gaussian
process models,” in International Joint Conference on Neural Networks (IJCNN). IEEE,
6-11 July 2014, pp. 1573–1578.

[19] J. Kocijan, R. Murray-Smith, C. E. Rasmussen, and A. Girard, “Gaussian process model
based predictive control,” in American Control Conference, vol. 3. IEEE, 2004, pp. 2214–
2219.
[20] A. Grancharova, J. Kocijan, and T. A. Johansen, “Explicit stochastic nonlinear predictive
control based on Gaussian process models,” in European Control Conference, 2007, pp. 2340–
2347.
[21] E. D. Klenske, M. N. Zeilinger, B. Scholkopf, and P. Hennig, “Gaussian process-based predictive control for periodic error correction,” IEEE Transactions on Control Systems Technology,
2015.
[22] G. Cao, E. M.-K. Lai, and F. Alam, “Gaussian process based model predictive control for
linear time varying systems,” in International Workshop on Advanced Motion Control (AMC
Workshop). IEEE, 22-24 April 2016.
[23] ——, “Gaussian process model predictive control of Unmanned Quadrotors,” in International
Conference on Control, Automation and Robotics (ICCAR). IEEE, 28-30 April 2016.
[24] M. P. Deisenroth, “Efficient reinforcement learning using Gaussian processes,” Ph.D. dissertation, Karlsruhe Institute of Technology, 2010.
[25] A. Girard, C. E. Rasmussen, J. Q. Candela, and R. Murray-Smith, “Gaussian process priors
with uncertain input – Application to multiple-step ahead time series forecasting,” in Advances
in Neural Information Processing Systems. MIT, 2003, pp. 545–552.
[26] J. Q. Candela, A. Girard, J. Larsen, and C. E. Rasmussen, “Propagation of uncertainty in
bayesian kernel models-application to multiple-step ahead forecasting,” in IEEE Proceedings
of International Conference on Acoustics, Speech, and Signal Processing (ICASSP), vol. 2.
IEEE, 2003, pp. II–701.
[27] J. Quiñonero-Candela and C. E. Rasmussen, “A unifying view of sparse approximate Gaussian process regression,” Journal of Machine Learning Research, vol. 6, pp. 1939–1959, 2005.
[28] F. Berkenkamp and A. P. Schoellig, “Learning-based robust control: Guaranteeing stability while improving performance,” in IEEE/RSJ Proceedings of International Conference on
Intelligent Robots and Systems (IROS), 2014.
[29] Y. Pan and E. Theodorou, “Probabilistic differential dynamic programming,” in Advances in
Neural Information Processing Systems, 2014, pp. 1907–1915.
[30] A. Grancharova, J. Kocijan, and T. A. Johansen, “Explicit stochastic predictive control of
combustion plants based on Gaussian process models,” Automatica, vol. 44, no. 6, pp. 1621–
1631, 2008.
[31] F. Tröltzsch, “Regular Lagrange multipliers for control problems with mixed pointwise
control-state constraints,” SIAM Journal on Optimization, vol. 15, no. 2, pp. 616–634, 2005.
[32] M. Diehl, H. J. Ferreau, and N. Haverbeke, “Efficient numerical methods for nonlinear MPC
and moving horizon estimation,” in International Workshop on assessment and future directions on Nonlinear Model Predictive Control. Pavia, Italy: Springer, 2008, pp. 391–417.

[33] S. Lucidi, M. Sciandrone, and P. Tseng, “Objective-derivative-free methods for constrained
optimization,” Mathematical Programming, vol. 92, no. 1, pp. 37–59, 2002.
[34] G. Liuzzi, S. Lucidi, and M. Sciandrone, “Sequential penalty derivative-free methods for
nonlinear constrained optimization,” SIAM Journal on Optimization, vol. 20, no. 5, pp. 2614–
2635, 2010.
[35] L. Yiqing, Y. Xigang, and L. Yongjian, “An improved PSO algorithm for solving non-convex
NLP/MINLP problems with equality constraints,” Computers & chemical engineering, vol. 31,
no. 3, pp. 153–162, 2007.
[36] O. Yeniay, “Penalty function methods for constrained optimization with genetic algorithms,”
Mathematical and Computational Applications, vol. 10, no. 1, pp. 45–56, 2005.
[37] S. J. Wright and M. J. Tenny, “A feasible trust-region sequential quadratic programming
algorithm,” SIAM journal on optimization, vol. 14, no. 4, pp. 1074–1105, 2004.
[38] Y.-h. Peng and S. Yao, “A feasible trust-region algorithm for inequality constrained optimization,” Applied mathematics and computation, vol. 173, no. 1, pp. 513–522, 2006.
[39] X. Zhang, J. Zhang, and L. Liao, “An adaptive trust region method and its convergence,”
Science in China Series A: Mathematics, vol. 45, no. 5, pp. 620–631, 2002.
[40] M. J. Tenny, S. J. Wright, and J. B. Rawlings, “Nonlinear model predictive control via
feasibility-perturbed sequential quadratic programming,” Computational Optimization and Applications, vol. 28, no. 1, pp. 87–121, 2004.
[41] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic
regulator for constrained systems,” Automatica, vol. 38, no. 1, pp. 3–20, 2002.
[42] Y. Wang and S. Boyd, “Fast model predictive control using online optimization,” IEEE Transactions on Control Systems Technology, vol. 18, no. 2, pp. 267–278, 2010.
[43] R. Fletcher, Practical methods of optimization, 2nd ed.
1987.

Wiley-Interscience Publication,

[44] Y. Pan and J. Wang, “Model predictive control of unknown nonlinear dynamical systems
based on recurrent neural networks,” IEEE Transactions on Industrial Electronics, vol. 59,
no. 8, pp. 3089–3101, 2012.
[45] L. Grüne and J. Pannek, Nonlinear model predictive control–Theory and Algorithms.
don, U.K: Springer-Verlag, 2011.

Lon-

[46] J. Kocijan and R. Murray-Smith, “Nonlinear predictive control with a Gaussian process
model,” in In R. Murray-Smith and R. Shorten (eds.), Switching and Learning in Feedback
Systems. Springer, 2005, pp. 185–200.

