1

Cautious Model Predictive Control using Gaussian

Process Regression
Lukas Hewing, Juraj Kabzan, Melanie N. Zeilinger

Abstract—Gaussian process (GP) regression has been widely learning the system dynamics as opposed to disturbances,
used in supervised machine learning due to its flexibility and however, imposes the challenge of propagating probability
inherent ability to describe uncertainty in function estimation. distributions of the state over the prediction horizon of the
In the context of control, it is seeing increasing use for modeling
of nonlinear dynamical systems from data, as it allows the direct controller. Efficient methods to evaluate GPs for Gaussian
inputs have been developed in [6], [7], [8]. By successively
arXiv:1705.10702v4 [[Link]] 29 Dec 2019

assessment of residual model uncertainty. We present a model

predictive control (MPC) approach that integrates a nominal employing these approximations at each prediction time step,
system with an additive nonlinear part of the dynamics modeled predictive control with GP dynamics has been presented in [9].
as a GP. Approximation techniques for propagating the state In [10] a piecewise linear approximate explicit solution for the
distribution are reviewed and we describe a principled way of
formulating the chance constrained MPC problem, which takes MPC problem of a combustion plant was presented. Applica-
into account residual uncertainties provided by the GP model tion of a one-step MPC with a GP model to a mechatronic
to enable cautious control. Using additional approximations for system was demonstrated in [11] and the use for fault-tolerant
efficient computation, we finally demonstrate the approach in a MPC was presented in [12]. A constrained tracking MPC
simulation example, as well as in a hardware implementation for for robotic applications, which uses a GP to improve the
autonomous racing of remote controlled race cars, highlighting
improvements with regard to both performance and safety over dynamics model from measurement data, was shown in [13]. A
a nominal controller. variant of this approach was realized in [14], where uncertainty
of the GP prediction is robustly taken into account using
Index Terms—Model Predictive Control, Gaussian Processes,
Learning-based Control confidence bounds. Solutions of GP-based MPC problems
using sequential quadratic programming schemes is discussed
in [15]. A simulation study for the application of GP-based
I. I NTRODUCTION predictive control to a drinking water network was presented
Many modern control approaches depend on accurate model in [16] and data efficiency of these formulations for learning
descriptions to enable safe and high performance control. Iden- controllers was recently demonstrated in [17].
tifying these models, especially for highly nonlinear systems, The goal of this paper is both to provide an overview of
is a time-consuming and complex endeavor. Often times, how- existing techniques and to propose extensions for deriving a
ever, it is possible to derive an approximate system model, e.g. systematic and efficiently solvable approximate formulation of
a linear model with adequate accuracy close to some operating the MPC problem with a GP model, which incorporates con-
point, or a simple model description from first principles. straints and takes into account the model uncertainty for cau-
In addition, measurement data from previous experiments or tious control. Unlike most previous approaches, we specifically
during operation is often available, which can be exploited to consider the combination of a nominal system description with
enhance the system model and controller performance. In this an additive GP part which can be of different dimensionality
paper, we present a model predictive control (MPC) approach, as the nominal model. This corresponds to a common scenario
which improves such a nominal model description from data in controller design, where an approximate nominal model
using Gaussian Processes (GPs) to safely enhance performance is known and can be improved on, offering a multitude of
of the system. advantages. A nominal model allows for rudimentary system
Learning methods for automatically identifying dynamical operation and the collection of measurement data, as well as
models from data have gained significant attention in the past the design of pre-stabilizing ancillary controllers which reduce
years [1] and more recently also for robotic applications [2], the state uncertainty in prediction. Additionally, it allows for
[3]. In particular, nonparametric methods, which have seen learning only specific effects from data, potentially reducing
wide success in machine learning, offer significant potential the dimensionality of the machine learning task. For example,
for control [4]. The appeal of using Gaussian Process regres- most dynamical models derived from physical laws include
sion for model learning stems from the fact that it requires an integrator chain, which can be directly represented in the
little prior process knowledge and directly provides a measure nominal dynamics and for which no nonlinear uncertainty
of residual model uncertainty. In predictive control, GPs were model has to be learned from data—it would rather add
successfully applied to improve control performance when unnecessary conservatism to the controller. The separation
learning periodic time-varying disturbances [5]. The task of between system and GP model dimension is therefore key for
the application of the GP-based predictive control scheme to
All authors are with the Institute for Dynamic Systems and Control, ETH higher order systems.
Zürich. [lhewing|kabzanj|mzeilinger]@[Link]
This work was supported by the Swiss National Science Foundation under The paper makes the following contributions. A review
grant no. PP00P2 157601 / 1. and compact summary of approximation techniques for prop-
 2

agating GP dynamics and uncertainties is provided and ex- Σw = diag([σ12 , . . . , σn2 d ]). We assume that both f and g are
tended to the additive combination with nominal dynamics. differentiable functions.
The approximate propagation enables a principled formulation The system is subject to state and input constraints X ⊆
of chance constraints on the resulting state distributions in Rnx , U ⊆ Rnu , respectively. The constraints are formulated
terms of probabilistic reachable sets [18]. In addition, the as chance constraints, i.e. by enforcing
nominal system description allows for reductin the GP model
learning to a subspace of states and inputs. We discuss sparse Pr(x(k) ∈ X ) ≥ px ,
GPs and a tailored selection of inducing points for MPC to Pr(u(k) ∈ U) ≥ pu ,
reduce the computational burden of the approach. The use
of some or all of these techniques is crucial to make the where px , pu are the associated satisfaction probabilities.
computationally expensive GP approach practically feasible
for control with sampling times in the millisecond range. C. Gaussian Process Regression
We finally present two application examples. The first is a
Gaussian process regression is a nonparametric framework
simulation of an autonomous underwater vehicle, illustrating
for nonlinear regression. A GP is a probability distribution
key concepts and advantages in a simplified setting. Second,
over functions, such that every finite sample of function values
we present a hardware implementation for autonomous racing
is jointly Gaussian distributed. We apply Gaussian processes
of remote controlled cars, showing the real-world feasibility
regression to infer the noisy vector-valued function g(x, u)
of the approach for complex high performance control tasks.
in (1) from previously collected measurement data of states
To the best of our knowledge this is the first hardware
and inputs {(xj , uj ), j = 0, . . . , M }. State-input pairs form
implementation of a Gaussian Process-based predictive control
the input data to the GP and the corresponding outputs are
scheme to a system of this complexity at sampling times of
obtained from the deviation to the nominal system model:
20 ms.
yj = g(xj , uj ) + wj = Bd† (xj+1 − f (xj , uj )) ,
II. P RELIMINARIES
A. Notation where Bd† is the Moore-Penrose pseudo-inverse. Note that the
measurement noise on the data points wj corresponds to the
The i-th element of a vector x is denoted [x]i . Similarly,
process noise in (1). With zj := [xTj , uTj ]T , the data set of the
[M ]ij denotes element ij of a matrix M , and [M ]·i , [M ]i·
GP is therefore given by
its i-th column or row, respectively. We use diag(x) to refer
to a diagonal matrix with entries given by the vector x. T
D = {y = [y0 , . . . , yM ] ∈ RM ×nd ,
The squared Euclidean norm weighted by M , i.e. xT M x T
z = [z0 , . . . , zM ] ∈ RM ×nz } .
is denoted kxk2M . We use boldface to emphasize stacked
quantities, e.g. a collection of vector-valued data in matrix Each output dimension is learned individually, meaning that
form. ∇f (z) is the gradient of f evaluated at z. The Pontryagin we assume the components of each yj to be independent, given
set difference is denoted A B = {a |a + b ∈ A ∀b ∈ B}. A the input data zj . Specifying a GP prior on g in each output
normally distributed vector x with mean µ and variance Σ dimension a ∈ {1, . . . , nd } with kernel k a (·, ·) and prior mean
is given by x ∼ N (µ, Σ). The expected value of x is E(x), function ma (·) results in normally distributed measurement
cov(x, y) is the covariance between vectors x and y, and the data with
variance var(x) = cov(x, x). We use p(x), (p(x|y)) to refer to [y]·,a ∼ N (m(z), Kzz a
+ Iσa2 ) , (2)
the (conditional) probability densities of x. Similary Pr(E),
(Pr(E | A)) denotes the probability of an event E (given A). where Kzz a
is the Gram matrix of the data points, i.e.
T
Realized quantities during closed loop control are time [Kzz ]ij = k a (zi , zj ) and ma (z) = [ma (z0 ), . . . , ma (zM )] .
a

indexed using parenthesis, i.e. x(k), while quantities in pre- The choice of kernel functions k a and its parameterization
diction use subscripts, e.g. xi is the predicted state i-steps is the determining factor for the inferred distribution of g
ahead. and is typically specified using prior process knowledge and
optimization [19], e.g. by optimizing the likelihood of the
B. Problem formulation observed data distribution (2). Throughout this paper we
consider the squared exponential kernel function
We consider the control of dynamical systems that can be
represented by the discrete-time model
 
2 1 T −1
k (zi , zj ) = σf,a exp − (zi − zj ) La (zi − zj ) , (3)
a

x(k+1) = f (x(k), u(k))+Bd (g(x(k), u(k))+w(k)) , (1) 2

where x(k) ∈ Rnx is the system state and u(k) ∈ Rnu the in which La is a positive diagonal length scale matrix and
control inputs at time k. The model is composed of a known
2
σf,a the signal variance. It is, however, straightforward to use
nominal part f and additive dynamics g that describe initially any other (differentiable) kernel function.
unknown dynamics of the system, which are to be learned The joint distribution of the training data and an arbitrary
from data and are assumed to lie in the subspace spanned by test point z in output dimension a is given by
Bd . We consider i.i.d. process noise wk ∼ N (0, Σw ), which Kzz + Iσa2 Kzz
 a   a a

m (z)
is spatially uncorrelated, i.e. has diagonal variance matrix p([y]a , [y]·,a ) = N , , (4)
ma (z) Kzza a
Kzz
 3

T
where [Kzz a
]j = k a (zj , z), Kzza
= (Kzza
) and similarly III. MPC C ONTROLLER D ESIGN
Kzz = k (z, z). The resulting distribution of [y]a condi-
a a
We consider the design of an MPC controller for system (1)
tioned on the observed data points
is again Gaussian with using a GP approximation d of the unknown function g:
p([y]a | [y]·,a ) = N µda (z), Σda (z) and
xi+1 = f (xi , ui ) + Bd (d(xi , ui ) + wi ) . (8)
−1
µda (z) = Kzz
a a
(Kzz + Iσa2 ) [y]·,a , (5a)
At each time step, the GP approximation evaluates to a
−1
Σda (z) = Kzz
a a
− Kzz a
(Kzz + Iσa2 ) a
Kzz . (5b) stochastic distribution according to the residual model uncer-
tainty and process noise, which is then propagated forward
The resulting multivariate GP approximation of the unknown in time. A stochastic MPC formulation allows the principled
function g(z) is then simply given by stacking the individual treatment of chance constraints, i.e. by imposing a prescribed
output dimensions, i.e. maximum probability of constraint violation. Denoting a state
and input reference trajectory by X r = {xr0 , . . . , xrN }, U r =
d(z) ∼ N µd (z), Σd (z) (6)


{ur0 , . . . , urN −1 }, respectively, the resulting stochastic finite
time optimal control problem can be formulated as
with µd = [µd1 , . . . , µdnd ]T and Σd = diag([Σd1 , . . . , Σdnd ]T ).
−1
N
!
Evaluating mean and variance in (5) has cost O(nd nz M ) X
min E lf (xN −xrN ) + l(xi −xri , ui −uri ) (9a)
and O(nd nz M 2 ), respectively and thus scales with the number Π(x) i=0
of data points. For large amounts of data points or fast
s.t. xi+1 = f (xi , ui ) + Bd (d(xi , ui ) + wi ), (9b)
real-time applications this can limit the use of GP models.
To overcome these issues, various approximation techniques ui = πi (xi ), (9c)
have been proposed, one class of which is sparse Gaussian Pr(xi+1 ∈ X ) ≥ px , (9d)
processes using inducing inputs [20], briefly outlined in the Pr(ui ∈ U) ≥ pu , (9e)
following.
x0 = x(k) , (9f)
for all i = 0, . . . , N − 1 with appropriate terminal cost
D. Sparse Gaussian Processes lf (xN ) and stage cost l(xi , ui ), where the optimization is
carried out over a sequence of input policies Π(x) =
Many sparse GP approximations make use of inducing
{π0 (x), . . . , πN −1 (x)}.
targets yind , inputs zind and conditional distributions to approx-
In the form (9), the optimization problem is computationally
imate the joint distribution (4) [21]. Many such approximations
intractable. In the following sections, we present techniques
exist, a popular of which is the Fully Independent Training
for deriving an efficiently solvable approximation. Specifically,
Conditional (FITC) [22], which we make use of in this
we use affine feedback policies with precomputed linear gains
paper. Given a selection of inducing inputs zind and using
−1 and show how this allows simple approximate propagation
the shorthand notation Qaζζ̃ := Kζz a
(Kzaind zind ) Kza ζ̃ the
ind ind of the system uncertainties in terms of mean and variance in
approximate posterior distribution is given by
prediction. Using these approximate distributions, we present
−1 a framework for reformulating chance constraints determin-
µ̃da (z) = Qazz (Qazz + Λ) [y]·,a , (7a)
istically and briefly discuss the evaluation of possible cost
−1
Σ̃da (z) = a
Kzz − Qazz (Qazz + Λ) Qzz (7b) functions.

with Λ = diag(Kzza
− Qazz + Iσa2 ). Concatenating the output
dimensions similar to (5) we arrive at the approximation A. Ancillary Linear State Feedback Controller
  Optimization over general feedback policies Π(x) is an in-
˜ ∼ N µ̃d (z), Σ̃d (z) .
d(z) finite dimensional optimization problem. A common approach
that integrates well in the presented framework is to restrict
Several of the matrices used in (7) can be precomputed and the policy class to linear state feedback controllers
the evaluation complexity becomes independent of the number πi (xi ) = Ki (µxi − xi ) + µui ,
of original data points. With M̃ being the number of inducing
points, this results in O(nd nz M̃ ) and O(nd nz M̃ 2 ) for the where µxi is the predicted mean of the state distribution. Using
predictive mean and variance, respectively. pre-selected gains Ki , we then optimize over µui , the mean of
There are numerous options for selecting the inducing the applied input. Note that due to the ancillary control law the
inputs, e.g. heuristically as a subset of the original data points, input ui applied in prediction is a random variable, resulting
by treating them as hyperparameters and optimizing their from the affine transformation of xi .
location [22], or letting them coincide with test points [23], In Fig. 1, the effect of state feedback on the propagation of
which is often referred to as transductive learning. In Sec- uncertainty is exemplified. It shows the evolution over time of
tion III-C we make use of such transductive ideas and propose a double integrator with a quadratic friction term inferred by a
a dynamic selection of inducing points, with a resulting local GP. The right plot displays mean and 2-σ variance of a number
approximation tailored to the predictive control task. of trajectories from repeated simulations of the system with
 4

0.5 open-loop where Σui = Ki Σxi KiT , Σxu i = Σxi KiT and ? denotes terms
1
closed-loop given by symmetry. Considering the covariances between
d([x]2 )

states, inputs and GP, i.e. Σxd and Σud , is of great importance

[x]1
0 0
for an accurate propagation of uncertainty when the GP
model d is paired with a nominal system model f . Using a
−1
−0.5 linearization of the nominal system dynamics around the mean
−2 0 2 0 10 20    x 
x u x u x µ
[x]2 k f (x, u) ≈ f (µ , µ ) + ∇f (µ , µ ) − u ,
u µ
Fig. 1: Propagation of uncertainty for double integrator with similar to extended Kalman filtering, this permits simple
GP trained on a nonlinear friction term. The left plot displays update equations for the state mean and variance based on
the GP with 2-σ confidence bound, while the right plot shows affine transformations of the Gaussian distribution in (10)
the mean and 2-σ variance of repeated simulation runs of the
system under an open-loop control sequence and closed-loop µxi+1 = f (µxi , µui ) + Bd µdi ,
state feedback control. T
Σxi+1 = [∇f (µxi , µui ) Bd ] Σi [∇f (µxi , µui ) Bd ] .

2) Gaussian Process Prediction from Uncertain Inputs: In

different noise realizations. Mean and variance under an open- order to define µdi , Σdi , Σxd
i and Σi different approximations
ud
loop control sequence derived from a linear quadratic infinite of the posterior of a GP from a Gaussian input have been
time optimal control problem are shown in red. The results proposed in the literature. In the following, we will give a brief
with a corresponding LQR feedback law applied in closed- overview of the most commonly used techniques, for details
loop are shown in blue. As evident, open-loop input sequences please refer to [6], [7], [8]. For comparison, we furthermore
can lead to rapid growth of uncertainty in the prediction, and state the computational complexity of these methods. For
thereby cause conservative control actions in the presence of notational convenience, we will use
chance constraints. Linear ancillary state feedback controllers  x  x
Σi Σxu
  xd 
are therefore commonly employed in stochastic and robust z µi z i zd Σi
µi = u , Σi = , Σi = .
MPC [24]. µi ? Σui Σudi
The adequate choice of ancillary feedback gains Ki is
a) Mean Equivalent Approximation: A straightforward
generally a hard problem, especially if the system dynamics
and computationally cheap approach is to evaluate (6) at the
are highly nonlinear. A useful heuristic is to consider a
mean, such that
linearization of the system dynamics around an approximate
prediction trajectory, which in MPC applications is typically µd = µd (µzi ) , (12a)
available using the solution trajectory of the previous time  zd i  
Σi 0
step. Feedback gains for the linearized system can be derived = . (12b)
Σdi Σd (µzi )
e.g. by solving a finite horizon LQR problem [25]. For mild
or stabilizing nonlinearities, a fixed controller gain Ki = K In [27] it was demonstrated that this can lead to poor ap-
can be chosen to reduce computational burden, as e.g. done proximations for increasing prediction horizons as it neglects
for the example in Fig. 1. the accumulation of uncertainty in the GP. The fact that this
approach neglects covariance between d and (x, u) can in
B. Uncertainty Propagation addition severely deteriorate the prediction quality when paired
with a nominal system f (x, u). The computationally most
1) Approximation as Normal Distributions: Because of
expensive operation is the matrix multiplication in (5) for each
stochastic process noise and the representation by a GP
output dimension of the GP, such that the complexity of one
model, future predicted states result in stochastic distributions.
prediction step is O(nd nz M 2 ).
Evaluating the posterior of a GP from an input distribution
is generally intractable and the resulting distribution is not b) Taylor Approximation: Using a first-order Taylor ap-
Gaussian [6]. While under certain assumptions on g, some proximation of (5), the expected value, variance and covari-
strict over-approximations of the resulting distributions ex- ance of the resulting distribution results in
ist [26], they are typically very conservative and computation- µdi = µd (µzi ) , (13a)
ally demanding. We focus instead on computationally cheap  zd  " T
#
and practical approximations at the cost of strict guarantees. Σi Σzi (∇µd (µzi ))
= T . (13b)
State, control input and nonlinear disturbance are approxi- Σdi Σd (µzi ) + ∇µd (µzi )Σzi (∇µd (µzi ))
mated as jointly Gaussian distributed at every time step.
h iT Compared to (12) this leads to correction terms for the
xTi uTi (di + wi )
T
∼ N (µi , Σi ) posterior variance and covariance, taking into account the
 x   x xu
gradient of the posterior mean and the variance of the input
Σxd

µi Σi Σi i of the GP. The complexity of one prediction step amounts
= N µui  ,  ? Σui Σud
i
 , (10) to O(nd n2z M 2 ). Higher order approximations are similarly
d
µi ? ? Σi + Σ w
d
possible at the expense of increased computational effort.
 5

0 depend on only a subset of states and inputs, such that the

computational burden can be significantly reduced. This is due
[x]2 + d([x]2 )

−0.5 to a reduction in the effective input dimension nz , and more

importantly due to a reduction the in necessary training points
Mean Eq.
M for learning in a lower dimensional space.
−1 Taylor Appr.
Exact MM Remark 1. With only slight notational changes, the presented
M. Carlo approximation methods similarly apply to prior mean and
−1.5
−3 −2 −1 0 1 kernel functions that are functions of only a subset of states
[x]2
and inputs.
Another significant reduction in the computational complex-
Fig. 2: Comparison of prediction methods for Gaussian input
ity can be achieved by employing sparse GP approaches, for
[x]2 . The posterior distribution [x]2 + d([x]2 ) is evaluated with
which in the following we present a modification tailored to
the different approximation methods of Section III-B2. For
predictive control.
reference, the true distribution is approximated by Monte Carlo
simulation.
C. Dynamic Sparse GPs for MPC
To reduce the computational burden one can make use of
c) Exact Moment Matching: It has been shown that sparse GP approximations as outlined in Section II-D, which
for certain mean and kernel functions, the first and second often comes with little deterioration of the prediction quality. A
moments of the posterior distribution can be analytically property of the task of predictive control is that it lends itself to
T
computed [6]. Under the assumption that zi = [xTi , uTi ] is transductive approximation schemes, meaning that we adjust
normally distributed, we can define the approximation based on the test points, i.e. the (future)
µd = E (d(zi )) , (14a) evaluation locations in order to achieve a high fidelity local
 zd i approximation. In MPC, some knowledge of the test points
cov (zi , d(zi ))
 
Σi
= (14b) is typically available in terms of an approximate trajectory
Σdi var (d(zi )) .
through the state-action space. This trajectory can, e.g., be
In particular, this is possible for a zero prior mean function given by the reference signal or a previous solution trajectory
and the squared exponential kernel. As the procedure exactly which will typically lie close to the current one. We therefore
matches first and second moments of the posterior distribution, propose to select inducing inputs locally at each sampling
it can be interpreted as an optimal fit of a Gaussian to the true time according to the prediction quality on these approximate
distribution. Note that the model formulation directly allows trajectories. Ideally, the local inducing inputs would be opti-
for the inclusion of linear prior mean functions, as they can be mized to maximize predictive quality, which is, however, not
expressed as a nominal linear system in (8), while preserving
exact computations of mean and variance. The computational
complexity of this approach is O(n2d nz M 2 ) [8]. In Fig. 2 the
different approximation methods are compared for a one-step
prediction of the GP shown in Fig. 1 with normally distributed
x2

x2

input. As evident, the predictions with Taylor Approximation

and Moment Matching are qualitatively similar, which will
typically be the case if the second derivative of the posterior
mean and variance of the GP is small over the input distribu- x1 x1
tion. It is furthermore apparent that the posterior distribution 0.4
full GP
Σd (xi )

is not Gaussian and that the approximation as a Gaussian

sparse GP
distribution leads to prediction error, even if the first two 0.2
moments are matched exactly. This can lead to the effect that
p

locally the Taylor Approximation provides a closer fit to the 0

2 4 6 8 10
underlying distribution. The Mean Equivalent Approximation i
is very conservative by neglecting the covariance between [x]2
and d([x]2 ). Note that depending on the sign of the covariance Fig. 3: Illustration of dynamic sparse approximation [28].
it can also be overly confident. Countor plot of the posterior variance of the full GP (top
All presented approximations scale directly with the input left) and dynamic sparse approximation (top right) with corre-
and output dimensions, as well as the number of data points sponding data points as black crosses. Trajectories planned by
and thus become expensive to evaluate in high dimensional an MPC are shown as solid red lines, while the dashed lines
spaces. This presents a challenge for predictive control ap- show the prediction of the previous time step and used in the
proaches with GP models, which in the past have mainly approximation, with the chosen inducing points indicated by
focused on relatively small and slow systems [9], [12]. Using black circles. The bottom plot shows the respective variances
a nominal model, however, it is possible to consider GPs that along the planned trajectory.
 6

computationally feasible in the targeted millisecond sampling can be performed online. For brevity, we concentrate on state
times. Instead, inducing inputs are selected heuristically along constraints, as input constraints can be treated analogously.
the approximate state input trajectory, specifically we focus 1) Half-space constraints: Consider the constraint set X
here on the MPC solution computed at a previous time step. given by a single half-space constraint X hs := x hT x ≤ b ,
 
To illustrate the procedure, we consider a simple double h ∈ Rn , b ∈ R+ . Considering the marginal distribution of the
integrator system controlled by an MPC. Fig. 3 shows the error in the direction of the half-space hT exi ∼ N (0, hT Σxi h)
variance Σd (x) of a GP trained on the systems’ states. Addi- enables us to use the quantile function of a standard Gaus-
tionally, two successive trajectories from an MPC algorithm sian random variable φ−1 (px ) at the needed probability of
are displayed, in which the solid red line is the current constraint satisfaction px to see that
prediction, while the dashed line is the prediction trajectory   q 
from the previous iteration. The plot on the left displays the x x −1
 T
R (Σi ) := e h e ≤ φ (px ) h Σi h
 T x

original GP, with data points marked as crosses, whereas on

the right we have the sparse approximation resulting from is an i-step PRS of probability level px . In this case, evaluating
placing inducing points along the previous solution trajectory, the Pontryagin difference in (15) is straightforward and we can
indicated by circles. The figure illustrates how full GP and directly define the tightened constraint on the state mean
sparse approximation match closely along the predicted tra-   q 
jectory of the system, while approximation quality far away hs x −1
 T
Z (Σi ) := z h z ≤ b − φ (px ) h Σi h .
 T x
from the trajectory deteriorates. Since current and previous
trajectory are similar, however, local information is sufficient Remark   2. A tightening for slab constraints
for computation of the MPC controller. X sl = x |hT x| ≤ b can be similarly derived as
   q 
px +1
Z sl (Σxi ) := z |hT z| ≤ b − φ−1

D. Chance Constraint Formulation hT Σxi h .
2
The tractable Gaussian approximation of the state and input
distribution over the prediction horizon in (10) can be used to 2) Polytopic Constraints: Consider polytopic constraints,
approximate the chance constraints in (9d) and (9e). given as the intersection of nj half-spaces
We reformulate the chance constraints on state and input X p = x hTj x ≤ bj ∀j = 1, . . . , nj .
 
w.r.t. their means µx , µu using constraint tightening based on
the respective errors exi = µxi − xi and eui = Ki (µxi − xi ). Making use of ellipsoidal PRS for Gaussian random vari-
For this we make use of probabilistic reachable sets [18], ables one can formulate a semidefinite program to tighten
an extension of the concept of reachable sets to stochastic constraints [31], [32], which is computationally demanding
systems, which is related to probabilistic set invariance [29], to solve online. Computationally cheaper approaches typically
[30]. rely on more conservative bounds. One such bound is given
by Boole’s inequality
Definition 1 (Probabilistic i-step Reachable Set). A set R is
said to be a probabilistic i-step reachable set (i-step PRS) of Pr(E1 ∧ E2 ) ≤ Pr(E1 ) + Pr(E2 ) , (16)
probability level p if
which allows for the definition of polytopic PRS based on
Pr(ei ∈ R | e0 = 0) ≥ p . individual half-spaces. We present two possibilities, the first
of which is based on bounding the violation probability of the
Given i-step PRS Rx of probability level px for the state
individual polytope faces, and a second which considers the
error exi and similarly Ru for the input eui , we can define
marginal distributions of exi .
tightened constraints on µxi and µui as
a) PRS based on polytope faces: Similar to the treatment
µxi ∈ Z = X Rx , (15a) of half space constraints, we can define a PRS on the state error
µui ∈V =U R , u
(15b) ex which is aligned with the considered constraints. Using
Boole’s inequality (16), we have that
where denotes the Pontryagin set difference. Satisfaction
of the tightened constraints (15) for the mean thereby implies Rx (Σxi ) =
satisfaction of the original constraints (9d) and (9e), i.e. when
    
 T −1 1−px q T x
e hj e ≤ φ 1− hj Σi hj ∀j = 1, . . . , nj
µxi ∈ Z we have Pr(xi = µxi +exi ∈ X ) ≥ Pr(exi ∈ Rx ) ≥ px .

nj
Under the approximation of a normal distribution of xi , the
is a PRS of probability level px . The tightening results in
uncertainty in each time step is fully specified by the variance
matrix Σxi and Σui . The sets can then be computed as functions Z p (Σxi ) =
of these variances, i.e. Rx (Σxi ) and Ru (Σxi ), for instance as    
1−px q T x

−1
 T
ellipsoidal confidence regions. The online computation and z hj z ≤ bj − φ
 1− hj Σi hj ∀j = 1, . . . , nj
nj
respective tightening in (15), however, is often computationally
prohibitive for general convex constraint sets. which scales with the number of faces of the polytope and
In the following, we present important cases for which can therefore be quite conservative if the polytope has many
a computationally cheap tightening is possible, such that it (similar) faces.
 7

b) PRS based on marginal distributions: Alternatively, Further examples include a saturating cost [8], risk sensitive
one can define a PRS based on the marginal distribution of costs [34] or radial basis function networks [17].
the error ex in each dimension, which therefore scales with the We refer to the evaluation of the expected value (9a) in
state dimension. We use Boole’s inequality (16) with Remark 2 terms of mean and variance as
directly on marginal distributions in each dimension of the
state to define a box-shaped PRS of probability level px E (l(xi − xri , ui − uri )) = c(µxi − xri , µui − uri , Σxi )
  q  and similarly cf (µxN − xrN , ΣxN ) for the terminal cost.
Rx (Σxi ) = e |[e]j | ≤ φ−1 (p̄) [Σxi ]j,j ∀j = 1, . . . , nx .

Remark 4. While many performance indices can be expressed

x+1
 as expected values, there exist various stochastic cost mea-
with p̄ = 1 − n1x − p2n . To compute the Pontryagin
x sures, such as conditional value-at-risk. Given the approximate
difference (15) we make use of the following Lemma. state distributions, many can be similarly considered.
Lemma  1. Let A = {x |Hx ≤ b} be a polytope and
B = e −r ≤ e ≤ r, r ∈ Rn+ , a box in Rn . F. Tractable MPC Formulation with GP Model
Then A B = {x |Hx ≤ b − |H|r}, where |H| is the
By bringing together the approximations of the previous
element-wise absolute value.
sections, the following tractable approximation of the MPC
Proof. From definition of the Pontryagin difference we have problem (9) can be derived:
N −1
A B = {x |x + e ∈ A ∀e ∈ B} X
min cf (µxN −xrN , ΣxN ) + c(µxi −xri , µui −uri , Σxi )
= {x |H(x + e) ≤ b ∀e ∈ B} {µui } i=0
  
s.t. µxi+1 = f (µxi , µui ) + Bd µdi ,

= x [H]j,· x ≤ [b]j − max[H]j,· e ∀j = 1, . . . , n
e∈B T
Σxi+1 = [∇f (µxi , µui ) Bd ] Σi [∇f (µxi , µui ) Bd ] ,
= {x |[H]j,· x ≤ [b]j − |[H]j,· |r ∀j = 1, . . . , n}
µxi+1 ∈ Z(Σxi+1 ),
= {x |Hx ≤ b − |H|r} ,
µui ∈ V(Σxi ),
proving the result. µdi , Σi according to (10) and (12), (13) or (14),
The resulting tightened set is therefore µx0 = x(k), Σx0 = 0
n  o (18)
Z p (Σxi ) = z Hz ≤ b̃(Σxi ) for i = 0, . . . , N − 1. The resulting optimal control law is


obtained in a receding horizon fashion as κ(x(k)) = µu∗ 0 ,

with b̃(Σxi ) = b − |H|φ−1 (p̄) diag(Σxi ), where diag(·) is the
p
where µu∗ 0 is the first element of the optimal control sequence
vector of diagonal elements and the square root of the vector
0 , . . . , µN } obtained from solving (18) at state x(k).
{µu∗ u∗
is taken element wise. The presented formulation of the MPC problem with a GP-
Remark 3. Treating polytopic constraints through (16) can based model results in a non-convex optimization problem.
lead to undesired and conservative individual constraints [33]. Assuming twice differentiability of kernel and prior mean
Practically, it can therefore be beneficial to constrain the function, second-order derivative information of all quantities
probability of violating each individual half-space constraint. is available. Problems of this form can typically be solved to
local optima using Sequential Quadratic Programming (SQP)
or nonlinear interior-point methods [35]. There exist a number
E. Cost Function
of general-purpose solvers that can be applied to solve this
Given the approximate joint normal distribution of state and class of optimization problems, e.g. the openly available
input, the cost function (9a) can be evaluated using standard IPOPT [36]. In addition, there are specialized solvers that
stochastic formulations. For simplicity, we focus here on costs exploit the structure of predictive control problems, such as
in form of expected values, which encompass most stochastic ACADO [37] or FORCES Pro [38], [39], which implements a
objectives typically considered. While the expected value for nonlinear interior point method. Derivative information can be
general cost functions needs to be computed numerically, there automatically generated using automated differentiation tools,
exist a number of functions for which evaluation based on such as CASADI [40].
mean and variance information can be analytically expressed In the following, we demonstrate the proposed algorithm
and is computationally cheap. The most prominent example and its properties using one simulation and one experimental
for tracking tasks is a quadratic cost on states and inputs example. The first is an illustrative example of an autonomous
l(xi − xri , ui − uri ) = kxi − xri k2Q + kui − uri k2R (17) underwater vehicle (AUV) which, around a trim point, is well
described by nominal linear dynamics but is subject to nonlin-
with appropriate weight matrices Q and R, typically satisfying ear friction effects for larger deviations. The second example
Q  0 and R  0, resulting in is a hardware implementation demonstrating the approach for
E (l(xi − xri , ui − uri )) = autonomous racing of miniature race cars which are described
by a nonlinear nominal model. While model complexity and
kµxi − xri k2Q + tr(QΣxi ) + kµui − uri k2R + tr(RΣui ) . small sampling times require additional approximations in this
 8

second example, we show that we can leverage key benefits

of the proposed approach in a highly demanding real-world 0.5
application.
[xr ]1
0 [x]1
IV. O NLINE L EARNING FOR AUTONOMOUS U NDERWATER u
V EHICLE
2 3 4 5 6
We consider the depth control of an autonomous under-
time [s]
water vehicle (AUV). The vehicle is described by a nonlinear
continuous-time model of the vehicle at constant surge velocity Fig. 4: Predicted pitch angle [x]1 and rudder deflection u.
as ground truth and for simulation purposes [41]. Around a Dotted lines are the mean prediction and shaded the 2-σ
trim point of purely horizontal movement, the states and input confidence region. The dashed lines show the corresponding
of the system are state and input constraints, while the black line is the pitch
[x]1 the pitch angle relative to trim point [rad], angle reference.
[x]2 the heave velocity relative to trim point [m/s],
[x]3 the pitch velocity [rad/s],
u the stern rudder deflection around trim point [rad]. We use a quadratic stage cost as in (17), with weight
We assume that an approximate linear system model is matrices Q = diag([1, 0, 10, 0.5T ]), R = 20 and a pre-
given, which in practice can be established using methods diction horizon N = 35 to track a pitch angle reference.
of linear system identification around the trim point of the The terminal cost P is chosen according to the solution of
system. Using a zero-order hold discretization of the linear the associated discrete-time algebraic Riccati equation, i.e.
part of the model with Ts = 100 ms we obtain the LQR cost. As ancillary linear controller Ki , an infinite
horizon LQR control law based on the linear nominal model
x(k+1) = Ax(k) + Bu(k) + Bd (g(x(k)) + w(k)) , is designed using the same weights as in the MPC and used
in all prediction steps i. This stabilizes the linear system and
which is in the form of (1). The nonlinearity results from
reduces uncertainty growth over the prediction horizon. To
friction effects when moving in the water, which is therefore
propagate the uncertainties associated with the GP we make
modeled as only affecting the (continuous time) velocity states,
use of the Taylor Approximation outlined in Section III-B2.
Ts2 Constraints are introduced based on Remark 3 considering
 
0 2
Bd = Ts 0  , g(x(k)) = g([x(k)]2 , [x(k)]3 ) : R2 → R2 . a maximum probability of individual constraint violation of
0 Ts 2.28%, corresponding to a 2-σ confidence bound. The MPC
optimization problem (18) is solved using FORCES Pro [38],
Note that the eigenvalues of A are given by λ = [39].
{1.1, 1.03, 1, 0.727}, i.e. the linear system has one integrating
and two unstable modes.
B. Results
In the considered scenario, the goal is to track reference
changes from the original zero set point to a pitch angle of Fig. 4 shows the prediction of the GP-based MPC for the
30◦ , back to zero and then to 45◦ . We furthermore consider a first reference change from 0◦ to 30◦ . Since no data was
safety state constraint on the pitch angle of at least 10◦ below collected on the state trajectory necessary for this change, the
the reference, as well as input constraints corresponding to predicted state and input trajectories are uncertain and a safety
±20◦ rudder deflection. In this example, we treat the case of margin to the state constraint at the end of the horizon is
online learning, that is we start without any data about the enforced. The resulting closed-loop trajectory, during which
nonlinearity g, collect measurement data during operation and the system learns from online data, is displayed in the top
enhance performance online. plot of Fig. 5a. For comparison, we run the same simulation
with a soft constrained linear MPC formulation with the
same cost function, which does not consider nonlinearities
A. GP-based Reference Tracking Controller or uncertainties in the dynamics, shown in the bottom plot.
GP data D = {y, z} is generated by calculating the The results demonstrate improved performance of the GP-
deviation of the linear model from the measured states, as based MPC over the linear MPC, especially with regard to
described in Section II-C, where the input data is reduced to constraint satisfaction. The constraint on minimum pitch angle
the velocity states zj = [[xj ]2 , [xj ]3 ]T . Data is continuously is violated under the linear MPC control law during both
updated during the closed-loop run. Specifically, we consider reference changes, even though the soft constraint is chosen
a new data point every 5 time steps, and keep track of 30 as an exact penalty function such that constraints are always
data points by discarding the oldest when adding a new point. satisfied, if possible. Fig. 5b exemplifies the residual model
The training data is initialized with 30 data points of zero error during the closed-loop simulation applying the GP-based
input and zero output. We employ a squared exponential ker- controller, as well as the predicted 2-σ residual error bound of
nel (3) for both output dimensions with fixed hyperparameters the GP. We observe the largest error during reference changes,
L1 = L2 = diag([0.35, 0.15]T ), and variances σf,1 2
= 0.04 which are anticipated by the GP-uncertainty, overall matching
and σf,2 = 0.25.
2
the resulting residual errors well.
 9

1 are given purely by kinematic relationships, we assume that

0.5
[x]1 the model uncertainty, as well as the process noise w, only
u
affect the velocity states vx , vy and ω of the system. From
0 physical considerations, we can also infer that the unmodeled
−0.5 dynamics do not depend on the position states, i.e. we assume
1
0.5 Bd = [0 I]T , g(x, u) = g(vx , vy , ω, p, δ) : R5 → R3 ,
0 resulting in
−0.5
0 10 20 30 40 x(k+1) = f (x(k), u(k))+Bd (g(x(k), u(k))+w(k)) .
time [s]
(a) GP-based MPC (top) and linear MPC (bottom). The solid gray line The system is subject to input constraints U, i.e. the steering
shows the reference value of the pitch angle, dashed in the respective angle is limited to lie in ±δmax and the duty cycle has to lie
color the state and input constraints. [−0.1, 1], where the negative values correspond to negative
·10−2 applied torques, i.e. breaking. Additionally, operation requires
4 the vehicle to stay on track, which is expressed as a state
2 constraint X on the car’s position.
e(k)

0
−2
B. GP-based Racing Controller
−4
0 10 20 30 40 We consider a race track given by its centerline and a fixed
time [s] track width. The centerline is described by a piecewise cubic
(b) Residual model error with GP-based MPC. Measured model error spline polynomial, which is parameterized by the path length
as black dots, 2-σ residual model uncertainty from GP shaded blue. Θ. Progress along the track is characterized by Θi , which
is introduced as an additional state and enters the considered
Fig. 5: Simulation results of GP-based MPC for autonomous
cost function linearly, encouraging a maximization of progress
underwater vehicle in an online learning scenario.
along the track. Given a Θi , we can evaluate the corresponding
T
centerline position C(Θ) = [Xc (Θ), Yc (Θ)] , such that the
V. AUTONOMOUS R ACING constraint for the car to stay within the track boundaries can
be expressed as
As a second example we consider an autonomous racing
scenario, in which the goal is to drive a car around a X (Θi ) := xXY  xXY − C(Θi ) ≤ r ⊂ R2 ,
 
track as quickly as possible, while keeping the vehicle safe,
i.e. while avoiding collision with the track boundaries. The where r is half the track width.
controller is based on a model predictive contouring control Based on approximate Gaussian distributions of the state
formulation [42], [43] which has been applied to the problem in prediction we have for the position error eXY
i = xXY
i −
of autonomous racing in [44]. Preliminary simulations results µi ∼ N (0, Σi ) and define an ellipsoidal PRS as
XY XY
were published in [28]. 
−1
n o
Rell (ΣXY ) = eXY (eXY )T (ΣXY ) eXY ≤ χ22 (px ) ,

i i
A. Car Dynamics
The race cars are modeled by continuous-time nominal where χ22 (px ) is the quantile function of the chi-squared
dynamics ẋ = f c (x, u) obtained from a bicycle model distribution with two degrees of freedom. In order to simplify
with nonlinear tire forces given by a simplified Pacejka tire constraint tightening of X (Θi ) we find an outer approximation
model [45], which results in the following states and inputs by a ball Rb ⊆ Rell as

x = [X, Y, Φ, vx , vy , ω]T , u = [p, δ]T ,
 
b XY XY  XY
 q
2
R (Σi ) = e XY
 e ≤ λmax Σi χ2 (px ) ,
T
with position xXY = [X, Y ] , orientation Φ, longitudinal and
lateral velocities vx and vy , and yaw rate ω. The inputs to the
system are the motor duty cycle p and the steering angle δ.
For details on the system modeling please refer to [44], [28].
For use in the MPC formulation, we discretize the sys-
tem using a Runge-Kutta method with a sampling time of
Ts = 20 ms. In order to account for model mismatch due to
inaccurate parameter choices and limited fidelity of this simple
model, we add g(x, u) capturing unmodeled dynamics, as well Fig. 6: Illustration of the constraint tightening procedure.
as additive Gaussian white noise w. Due to the structure of the The effective track radius is adjusted based on the predicted
nominal model, i.e. since the dynamics of the first three states position uncertainty.
 10

3 3

Velocity [m/s]
Velocity [m/s]
2 2

1 1

0 0
(a) Nominal controller (b) GP-based controller
Fig. 7: Comparison of racelines with nominal and GP-based controller. The color indicates the 2-norm of the velocity, the red
cross is the starting point of the race car.

where λmax (·) is the maximum eigenvalue which can be erratic and there are a number of small collisions with the track
readily computed since ΣXYi is a 2 by 2 matrix. The necessary boundaries. This can be observed in the racelines plotted in
constraint tightening can therefore be expressed as Fig. 7a showing 20 laps run with the nominal controller. Using
the collected data, we train the GP error model d using 325
Z(Θi , ΣXY
i ) = X (Θi ) Rb (ΣXY ) (19) data points. We infer the hyperparameters in (3) as well as the
 XY  XY i XY


= z  z − C(Θi ) ≤ r̃ Σi
, noise level Σw using maximum likelihood optimization. The
q resulting racelines of 20 laps with the GP-based controller
where r̃ ΣXY = r − χ22 (px )λmax ΣXY . Fig. 6 exem-



i i are displayed in Figure 7b, generally showing a much more
plifies the predicted evolution of the cars position and the consistent and safe racing behavior. In particular, it can be
resulting constraint tightening. seen that almost all of the systematic and persistent problems
The state is extended by previously applied inputs and large in the raceline of the nominal controller can be alleviated.
input changes are additionally penalized. We make use of the Fig. 8 shows the encountered dynamics error in the yaw-rate
Taylor approximation (13) to propagate uncertainties without and the predicted error during the first lap with the sparse GP-
the use of an ancillary state feedback controller, i.e. K = 0. based controller. Mean and residual uncertainty predicted by
The prediction horizon is chosen as N = 30 and we formulate the GP matches the encountered errors well. It is important
the chance constraints (19) with χ22 (px ) = 1. To reduce to note that the apparent volatility in the plot is not due to
conservatism of the controller, constraints are instead only overfitting, but is instead due to fast changes in the input
tightened for the first 20 prediction steps and are applied to the and matches the validation data, i.e. the measured errors. To
mean for the remainder of the prediction horizon, similar to the quantify performance of the proposed controllers we compare
method used in [46]. We reduce computation times by making in Table I average lap time T l , minimum lap time Tl,min as
use of the dynamic sparse approximations with 10 inducing well as the average 2-norm error in the system dynamics kek,
points as outlined in Section III-C, placing the inducing inputs i.e. the difference between the mean state after one prediction
regularly along the previous solution trajectory.
To ensure real-time feasibility of the approach for the
sampling time of 20ms two additional approximations are 4
applied. The variance dynamics are pre-evaluated based on the
previous MPC solution, which enables the pre-computation 2
dω [rad/s]

of state constraints (19) such that they remain fixed during

optimization. Additionally, we neglect the mean prediction 0
of the lateral velocity error vy since the state is difficult to
−2
estimate reliably and the error generally small, such that we
observed no improvement when including it in the control −4
formulation. 6 8 10 12
time [s]

C. Results Fig. 8: Dynamic sparse GP compensation of the yaw-rate error

We start out racing the car using the nominal controller, i.e. with 10 inducing inputs during the first race lap. The black
the controller without an added GP term, which therefore does dots show the measured error on the yaw rate at each time
not consider uncertainties for constraint tightening. Since the step, while the blue line shows the error predicted by the GP.
nominal model is not well tuned, driving behavior is somewhat The shaded region is the 2-σ confidence interval.
 11

TABLE I: Experimental results R EFERENCES

Controller T l [s] Tl,min [s] kek [−] Tc [ms] Tc < 20 ms [1] K. J. Hunt, D. Sbarbaro, R. Zbikowski, and P. J. Gawthrop, “Neural
networks for control systems-A survey,” Automatica, vol. 28, no. 6, pp.
Nominal 10.32 9.65 0.73 18.2 76.3% 1083–1112, 1992.
GP-based 9.61 9.27 0.33 17.2 99.8% [2] D. Nguyen-Tuong and J. Peters, “Model learning for robot control: a
survey,” Cognitive Processing, vol. 12, no. 4, pp. 319–340, Nov. 2011.
[3] O. Sigaud, C. Salan, and V. Padois, “On-line regression algorithms
for learning mechanical models of robots: A survey,” Robotics and
step and the realized state, e(k+1) = µx1 − x(k+1). We see Autonomous Systems, vol. 59, no. 12, pp. 1115–1129, 2011.
[4] G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, and L. Ljung,
that the GP-based controller is able to improve significantly on “Kernel methods in system identification, machine learning and function
all these quantities, with an average lap time improvement of estimation: A survey,” Automatica, vol. 50, no. 3, pp. 657–682, 2014.
0.71 s, or almost 7%, which constitutes a large improvement [5] E. D. Klenske, M. N. Zeilinger, B. Schölkopf, and P. Hennig, “Gaussian
in the considered racing task. This is in part due to the process-based predictive control for periodic error correction,” IEEE
Transactions on Control Systems Technology, vol. 24, no. 1, pp. 110–
improved system model, as evident in the average dynamics 121, 2016.
error kek, but also due to the cautious nature of the controller, [6] J. Quionero-Candela, A. Girard, and C. Rasmussen, “Prediction at an
which helps to further reduce collisions and large problems uncertain input for Gaussian processes and relevance vector machines -
application to multiple-step ahead time-series forecasting,” Max Planck
in the raceline. Due to the cautious nature, the minimum lap Institute for Biological Cybernetics, Tübingen, Germany, Tech. Rep.
time gains are slightly less pronounced. In fact, the nominal IMM-2003-18, 2003.
controller consistently displays higher top speeds, which often [7] M. Kuß, “Gaussian process models for robust regression, classification,
and reinforcement learning,” Ph.D. dissertation, Technische Universität
times, however, lead to significant problems at the breakpoint Darmstadt, Darmstadt, 2006.
to a slow corner. Computation times are reported as average [8] M. Deisenroth, “Efficient reinforcement learning using Gaussian pro-
solve times T c and the percentage of solutions in under 20 ms, cesses,” Ph.D. dissertation, Karlsruhe Institute for Technology, Karl-
sruhe, 2010.
Tc < 20 ms. Average solution times of nominal and GP-based [9] J. Kocijan, R. Murray-Smith, C. E. Rasmussen, and A. Girard, “Gaussian
controller are similar. The percentiles of solutions in under process model based predictive control,” in American Control Confer-
20 ms, however, differ significantly. This is mainly due to ence, Jun. 2004, pp. 2214–2219.
frequent large re-planning actions occurring with the nominal [10] A. Grancharova, J. Kocijan, and T. A. Johansen, “Explicit stochastic
predictive control of combustion plants based on Gaussian process
controller. models,” Automatica, vol. 44, no. 6, pp. 1621–1631, 2008.
The results therefore demonstrate that the presented GP- [11] G. Cao, E. M.-k. Lai, and F. Alam, “Gaussian process model predic-
based controller can significantly improve performance while tive control of unmanned quadrotors,” in International Conference on
Control, Automation and Robotics, 2016.
maintaining safety in a hardware implementation of a complex [12] X. Yang and J. M. Maciejowski, “Fault tolerant control using Gaussian
system with small sampling times. processes and model predictive control,” International Journal of Ap-
plied Mathematics and Computer Science, vol. 25, no. 1, pp. 133–148,
2015.
VI. C ONCLUSION [13] C. J. Ostafew, A. P. Schoellig, T. D. Barfoot, and J. Collier, “Learning-
based nonlinear model predictive control to improve vision-based mobile
robot path tracking,” Journal of Field Robotics, vol. 33, no. 1, pp. 133–
The paper discussed the use of Gaussian process regression 152, 2016.
to learn nonlinearities for improved performance in model [14] C. J. Ostafew, A. P. Schoellig, and T. D. Barfoot, “Robust constrained
predictive control. Combining GP dynamics with a nominal learning-based NMPC enabling reliable mobile robot path tracking,” The
International Journal of Robotics Research, vol. 35, no. 13, pp. 1547–
system allows for learning only parts of the dynamics, which 1563, 2016.
is key for keeping the required number of data points and [15] G. Cao, E. M.-K. Lai, and F. Alam, “Gaussian process model predictive
computational complexity of the GP regression feasible for control of unknown non-linear systems,” IET Control Theory & Appli-
cations, vol. 11, pp. 703–713(10), 2017.
online control. Approximation methods for the propagation
[16] Y. Wang, C. Ocampo-Martinez, and V. Puig, “Stochastic model pre-
of the state distributions over the prediction horizon were dictive control based on Gaussian processes applied to drinking water
reviewed and it was shown how this enables a principled networks,” IET Control Theory Applications, vol. 10, no. 8, pp. 947–955,
treatment of chance constraints on both states and inputs. 2016.
[17] S. Kamthe and M. P. Deisenroth, “Data-efficient reinforcement learning
A simulation example as well as a hardware implementation with probabilistic model predictive control,” in International Conference
have shown how the proposed formulations provide cautious on Artificial Intelligence and Statistics, 2018.
control with improved performance for medium-sized systems [18] L. Hewing and M. N. Zeilinger, “Stochastic model predictive control
for linear systems using probabilistic reachable sets,” Conference on
with low sampling times. In particular, we have demonstrated Decision and Control, 2018, accepted.
in experiment that both performance and safety in an au- [19] C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine
tonomous racing setting can be significantly improved by using learning. The MIT Press, 2006.
[20] J. Quiñonero-Candela, C. E. Rasmussen, and C. K. Williams, “Approx-
cautious data-driven techniques. imation methods for Gaussian process regression,” Microsoft Research,
Tech. Rep. MSR-TR-2007-124, 2007.
[21] J. Quiñonero-Candela, C. E. Rasmussen, and R. Herbrich, “A unifying
ACKNOWLEDGEMENT view of sparse approximate Gaussian process regression,” Journal of
Machine Learning Research, vol. 6, pp. 1935–1959, 2005.
We would like to thank the Automatic Control Laboratory [22] E. Snelson and Z. Ghahramani, “Sparse gaussian processes using
(IfA) at ETH Zurich and in particular Alexander Liniger for his pseudo-inputs,” in Advances in Neural Information Processing Systems,
Y. Weiss, B. Schölkopf, and J. C. Platt, Eds., 2006, pp. 1257–1264.
valuable input and support with the hardware implementation [23] V. Tresp, “A Bayesian committee machine,” Neural Computation,
on the miniature race cars. vol. 12, no. 11, pp. 2719–2741, Nov. 2000.
 12

[24] A. Bemporad and M. Morari, “Robust model predictive control: A

survey,” Robustness in identification and control, vol. 245, pp. 207–226,
1999.
[25] J. Rawlings and D. Mayne, Model Predictive Control: Theory and
Design. Nob Hill Pub., 2009.
[26] T. Koller, F. Berkenkamp, M. Turchetta, and A. Krause, “Learning-
based model predictive control for safe exploration and reinforcement
learning,” arXiv:0707.3168, 2018.
[27] A. Girard, C. E. Rasmussen, and R. Murray-Smith, “Gaussian process
priors with uncertain inputs : Multiple-step-ahead prediction,” Depart-
ment of Computing Science, University of Glasgow, Glasgow, Tech.
Rep., 2002.
[28] L. Hewing, A. Liniger, and M. N. Zeilinger, “Cautious NMPC with
Gaussian process dynamics for miniature race cars,” European Control
Conference, 2018.
[29] E. Kofman, J. A. De Don, and M. M. Seron, “Probabilistic set invariance
and ultimate boundedness,” Automatica, vol. 48, no. 10, pp. 2670–2676,
2012.
[30] L. Hewing, A. Carron, K. Wabersich, and M. N. Zeilinger, “On a
correspondence between probabilistic and robust invariant sets for linear
systems,” European Control Conference, 2018.
[31] Z. Zhou and R. Cogill, “Reliable approximations of probability-
constrained stochastic linear-quadratic control,” Automatica, vol. 49,
no. 8, pp. 2435 – 2439, 2013.
[32] G. Schildbach, P. Goulart, and M. Morari, “Linear controller design for
chance constrained systems,” Automatica, vol. 51, pp. 278 – 284, 2015.
[33] M. Lorenzen, F. Dabbene, R. Tempo, and F. Allgwer, “Constraint-
tightening and stability in stochastic model predictive control,” IEEE
Transactions on Automatic Control, vol. 62, no. 7, pp. 3165–3177, 2017.
[34] P. Whittle, “Risk-sensitive linear/quadratic/gaussian control,” Advances
in Applied Probability, vol. 13, no. 4, p. 764777, 1981.
[35] M. Diehl, H. J. Ferreau, and N. Haverbeke, “Efficient numerical methods
for nonlinear MPC and moving horizon estimation problem formula-
tion,” Nonlinear model predictive control, pp. 391–417, 2009.
[36] A. Wächter and L. T. Biegler, “On the implementation of an interior-
point filter line-search algorithm for large-scale nonlinear programming,”
Mathematical programming, vol. 106, no. 1, pp. 25–57, 2006.
[37] B. Houska, H. J. Ferreau, and M. Diehl, “ACADO Toolkit – An Open
Source Framework for Automatic Control and Dynamic Optimization,”
Optimal Control Applications and Methods, vol. 32, no. 3, pp. 298–312,
2011.
[38] A. Domahidi and J. Jerez, “FORCES Professional,” embotech GmbH
([Link] Jul. 2014.
[39] A. Zanelli, A. Domahidi, J. Jerez, and M. Morari, “FORCES NLP:
an efficient implementation of interior-point methods for multistage
nonlinear nonconvex programs,” International Journal of Control, 2017.
[40] J. Andersson, “A general-purpose software framework for dynamic
optimization,” Ph.D. dissertation, KU Leuven, 2013.
[41] M. S. Naik and S. N. Singh, “State-dependent Riccati equation-based
robust dive plane control of AUV with control constraints,” Ocean
Engineering, vol. 34, no. 11-12, pp. 1711–1723, 2007.
[42] T. Faulwasser, B. Kern, and R. Findeisen, “Model predictive path-
following for constrained nonlinear systems,” Conference on Decision
and Control and Chinese Control Conference, pp. 8642–8647, 2009.
[43] D. Lam, C. Manzie, and M. Good, “Model predictive contouring
control,” in Conference on Decision and Control, 2010, pp. 6137–6142.
[44] A. Liniger, A. Domahidi, and M. Morari, “Optimization-based au-
tonomous racing of 1:43 scale RC cars,” Optimal Control Applications
and Methods, vol. 36, no. 5, pp. 628–647, 2015.
[45] H. B. Pacejka and E. Bakker, “The magic formula tyre model,” Vehicle
System Dynamics, vol. 21, no. sup001, pp. 1–18, 1992.
[46] J. V. Carrau, A. Liniger, X. Zhang, and J. Lygeros, “Efficient implemen-
tation of randomized MPC for miniature race cars,” European Control
Conference, pp. 957–962, 2016.