# A Data-Driven Automatic Tuning Method for MPC under Uncertainty using
Constrained Bayesian Optimization

Farshud Sorourifar  Georgios Makrygirgos  Ali Mesbah  Joel A. Paulson
Department of Chemical and Biomolecular Engineering, The Ohio State
University, Columbus, OH 43210, USA  Department of Chemical and Biomolecular
Engineering, University of California, Berkeley, CA 94720, USA

###### Abstract

The closed-loop performance of model predictive controllers (MPCs) is
sensitive to the choice of prediction models, controller formulation, and
tuning parameters. However, prediction models are typically optimized for
prediction accuracy instead of performance, and MPC tuning is typically done
manually to satisfy (probabilistic) constraints. In this work, we demonstrate
a general approach for automating the tuning of MPC under uncertainty. In
particular, we formulate the automated tuning problem as a constrained black-
box optimization problem that can be tackled with derivative-free
optimization. We rely on a constrained variant of Bayesian optimization (BO)
to solve the MPC tuning problem that can directly handle noisy and expensive-
to-evaluate functions. The benefits of the proposed automated tuning approach
are demonstrated on a benchmark continuously stirred tank reactor example.

###### keywords:

Model predictive control; Constrained Bayesian optimization; Automated tuning

##  1 Introduction

Model predictive control (MPC) is one of the most widely used methods for the
control of constrained multivariable systems (Rawlings and Mayne, 2009). The
closed-loop performance of MPC strongly depends on (i) the quality of its
underlying process model used to make dynamic predictions; (ii) the
formulation of the objective and constraints; and (iii) the choice of several
tuning parameters (e.g., prediction horizon, weights in cost function, and
constraint backoff terms) (Garriga and Soroush, 2010; Paulson and Mesbah,
2018). A major challenge in MPC tuning arises from the non-trivial
relationships between the tuning parameters and the closed-loop control
performance and constraint satisfaction (Lu et al., 2020). As such, MPC tuning
via trial-and-error or other heuristic strategies may require a significant
number of closed-loop simulations, which can quickly become prohibitive
especially when system uncertainties are considered.

Recently, there has been a renewed interest in automated strategies for
controller tuning using Bayesian optimization (BO); e.g., see Berkenkamp et
al. (2016); Bansal et al. (2017); Neumann-Brosig et al. (2019); Forgione et
al. (2019); Khosravi et al. (2020); Lu et al. (2020). BO has emerged as a
powerful derivative-free method for optimizing black-box functions in various
applications (Shahriari et al., 2015), most notably for hyperparameter
selection of machine learning algorithms (Snoek et al., 2012). BO is deemed
particularly useful for solving MPC tuning problems since it can accommodate a
mixture of continuous and discrete decision variables and also overcome the
limitations of alternative derivative-free optimization methods such as
genetic algorithms and particle swarm optimization (Garriga and Soroush,
2010).

In any model-based control approach, the chosen model plays a pivotal role in
the design of the controller; however, identification of the model has been
traditionally separated from controller design (even in the context of
automated tuning strategies). An alternative idea that has been gaining
popularity in recent years is to treat the identification process as a “tuning
parameter” that should take into account the intended control application. The
identification for control (I4C) rationale has been heavily studied in the
context of fixed-order controllers for linear time-invariant systems (Gevers,
2005). Recently, the I4C methodology was extended to MPC in (Piga et al.,
2019) wherein BO is used to search for the best (parametrized) prediction
model for MPC by directly optimizing closed-loop performance determined from
experimental (or high-fidelity simulation) data. We take a similar perspective
in this work and show how the framework is applicable to more general MPC
formulations (including economic, nonlinear, robust, and/or output feedback).

The main focus of this paper is on integrated performance-driven model
learning and MPC tuning under uncertainty. The explicit incorporation of
uncertainty into the formulation of the automated tuning optimization problem
(through both the objective and constraints) is one of the main contributions
of this work. In particular, we show how independent uncertainty samples can
be used to obtain effectively noisy measurements of objective and constraint,
which can be accounted for in the BO method through proper estimation of the
noise variance. The second contribution is to leverage a constrained variant
of BO (Gardner et al., 2014; Hernández-Lobato et al., 2016) to directly handle
output constraints. Constrained BO has been shown to overcome challenges with
barrier methods (that penalize constraint violation in the objective) since
the constraints are separately modeled and thus provide an independent
representation of the feasibility region.

##  2 Problem Formulation

Consider an uncertain, time-invariant nonlinear system of the following form

$\displaystyle x_{k+1}$ $\displaystyle=f_{\text{plant}}(x_{k},u_{k},w_{k}),$
(1a) $\displaystyle y_{k}$ $\displaystyle=h_{\text{plant}}(x_{k},v_{k}),$ (1b)

where $k\in\mathbb{N}_{0}$ denotes the discrete time index,
$x_{k}\in\mathbb{R}^{n_{x}}$ denotes the system state,
$u_{k}\in\mathbb{R}^{n_{u}}$ denotes the control input,
$y_{k}\in\mathbb{R}^{n_{y}}$ denotes the measured output,
$w_{k}\in\mathbb{R}^{n_{w}}$ denotes the process noise, and
$v_{k}\in\mathbb{R}^{n_{v}}$ denotes the measurement noise. The initial state
$x_{0}$, the process noise sequence $w_{0:T-1}$, and the measurement noise
sequence $v_{0:T}$ are assumed to be random variables for some finite horizon
$T\in\mathbb{N}$. Control performance $J(y_{1:T},u_{1:T})$ is assumed to be
some known function of the output and input sequences. Since the output
depends on the realization of the uncertainty
$\delta=(x_{0},w_{0:T-1},v_{0:T})$, our actual objective is to minimize
expected performance cost

$\displaystyle\mathbb{E}_{\delta}\\{J(y_{1:T},u_{1:T})\\}=\int_{\Delta}J(y_{1:T},u_{1:T})p(\delta)d\delta,$
(2)

where $p(\delta)$ denotes the probability density function and $\Delta$ is the
support of the random vector $\delta$. The controlled system should also
satisfy hard input constraints

$\displaystyle
u_{0:T}\in\mathbb{U}_{0:T}=\mathbb{U}_{0}\times\cdots\times\mathbb{U}_{T},$
(3)

and joint chance constraints on the output

$\displaystyle\mathbb{P}_{\delta}\\{y_{1:T}\in\mathbb{Y}_{1:T}=\mathbb{Y}_{1}\times\cdots\times\mathbb{Y}_{T}\\}\geq
1-\epsilon,$ (4)

where $\epsilon\in[0,1]$ denotes the allowed violation probability. The ideal
controller design problem can be formulated in terms of the following
stochastic optimization problem:

$\displaystyle\min_{\pi\in\Pi}$
$\displaystyle~{}~{}\mathbb{E}_{\delta}\\{J(y_{1:T},u_{1:T})\\},$ (5a) s.t.
$\displaystyle~{}~{}x_{k+1}=f_{\text{plant}}(x_{k},u_{k},w_{k}),$
$\displaystyle\forall k\in\mathbb{N}_{0}^{T-1},$ (5b)
$\displaystyle~{}~{}y_{k}=h_{\text{plant}}(x_{k},v_{k}),$
$\displaystyle\forall k\in\mathbb{N}_{0}^{T},$ (5c)
$\displaystyle~{}~{}u_{k}=\pi_{k}(y_{0:k})\in\mathbb{U}_{k},$
$\displaystyle\forall k\in\mathbb{N}_{0}^{T},$ (5d)
$\displaystyle~{}~{}\mathbb{P}_{\delta}\\{y_{1:T}\in\mathbb{Y}_{1:T}\\}\geq
1-\epsilon,$ (5e)

where $\pi=\\{\pi_{0},\pi_{1},\ldots,\pi_{T}\\}$ denotes the control policy
that is composed of a sequence of control laws
$\pi_{k}:\mathbb{R}^{n_{y}(k+1)}\to\mathbb{R}^{n_{u}}$, which are arbitrary
functions of the history of measured outputs and $\Pi$ denotes the set of
possible control policy structures. This problem is extremely challenging for
several reasons including the set of functions $\pi\in\Pi$ is infinite
dimensional and the probabilistic operators $\mathbb{E}_{\delta}\\{\cdot\\}$
and $\mathbb{P}_{\delta}\\{\cdot\\}$ cannot be computed exactly for general
nonlinear systems. Another key complication considered in this work is that
functions $f_{\text{plant}}(\cdot)$ and $h_{\text{plant}}(\cdot)$ and
uncertainty distribution $\delta\sim p(\delta)$ may not be explicitly known.
Thus, (5) belongs to the challenging class of constrained black-box
optimization under uncertainty problems.

Since the optimization problem (5) cannot be solved directly, one must resort
to a heuristic controller design strategy such as model predictive control
(MPC) (Rawlings and Mayne, 2009) The standard MPC design strategy involves:
(i) identification of a control-relevant model of the system (1) using first-
principles or system identification techniques; (ii) specification of the MPC
controller using this control-relevant model to make internal predictions
about the future system behavior to select the optimal control inputs; and
(iii) trial-and-error experimentation to select the remaining MPC tuning
parameters (e.g., prediction horizon, input and output weight matrices, and
constraint backoffs) that can have a strong influence on closed-loop
performance. Not only are the tuning parameters selected in an ad hoc way in
this approach, it is difficult to know how the quality of the predictions of
the identified model will translate onto control performance. For example, the
I4C methodology mentioned previously is based on the fact that the best model
for model-based controller design may not be the one that provides the
smallest output prediction errors (Gevers, 2005).

To address these challenges, we take a performance-driven perspective in this
work in which the model of (1) is treated as a degree of freedom that can be
used to minimize the closed-loop performance cost. The proposed MPC
parametrization is summarized in the next section, which is followed by a
description of an automated procedure for simultaneously selecting the
internal model and other relevant MPC tuning parameters.

##  3 Policy Approximation: Parametrized Model Predictive Control

One of the most critical components of MPC is the model used to predict the
outputs given a sequence of inputs. Let the control-relevant model be defined
as

$\displaystyle s_{k+1}$
$\displaystyle=f_{\text{mod}}(s_{k},u_{k};\theta_{m}),$ (6a) $\displaystyle
y_{k}$ $\displaystyle=h_{\text{mod}}(s_{k};\theta_{m}),$ (6b)

where $s_{k}\in\mathbb{R}^{n_{s}}$ is the state of the control-relevant model
at time step $k\in\mathbb{N}_{0}$ and $\theta_{m}\in\mathbb{R}^{n_{\theta,m}}$
is the set of model parameters. Although written in state-space form, (6)
could be a realization of an input-output black-box model such as a nonlinear
auto-regressive model with exogenous variables (NARX) in which case
$\theta_{m}$ would represent the coefficients of the selected basis functions.

Given the approximate input-output system model (6), the MPC controller solves
the following finite-horizon optimal control problem at each sampling time
$k\in\mathbb{N}_{0}^{T}$:

$\displaystyle\min_{\mathbf{U}_{k}}$
$\displaystyle~{}~{}\textstyle\sum_{i=1}^{N_{p}}\tilde{\ell}(y_{i|k},u_{i|k};\theta_{o}),$
(7a) s.t.
$\displaystyle~{}~{}s_{i+1|k}=f_{\text{mod}}(s_{i|k},u_{i|k};\theta_{m}),$
$\displaystyle\forall i\in\mathbb{N}_{0}^{N_{p}-1},$ (7b)
$\displaystyle~{}~{}y_{i|k}=h_{\text{mod}}(s_{i|k};\theta_{m}),$
$\displaystyle\forall i\in\mathbb{N}_{1}^{N_{p}},$ (7c)
$\displaystyle~{}~{}(y_{i|k}+\theta_{b},u_{i|k})\in\mathbb{Y}_{k+i}\times\mathbb{U}_{k+i},$
$\displaystyle\forall i\in\mathbb{N}_{1}^{N_{p}},$ (7d)
$\displaystyle~{}~{}s_{0|k}=\hat{s}_{k},$ (7e)

where $N_{p}$ denotes the prediction horizon; $s_{i|k}$ and $u_{i|k}$ are,
respectively, the predicted state and control input at time $k+i$ given
information at time $k$; $\mathbf{U}_{k}=\\{u_{0|k},\ldots,u_{N-1|k}\\}$ is
the control input sequence; $\hat{s}_{k}\in\mathbb{R}^{n_{s}}$ is the current
state estimate;
$\tilde{\ell}:\mathbb{R}^{n_{y}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\theta,o}}\to\mathbb{R}$
is the stage cost function parametrized by
$\theta_{o}\in\mathbb{R}^{n_{\theta,o}}$; and
$\theta_{b}\in\mathbb{R}^{n_{\theta,b}}$ are the constraint backoff parameters
that can be used to improve the robustness properties of the controller, as
discussed in, e.g., (Paulson and Mesbah, 2018). The closed-loop performance
index $J(y_{1:T},u_{1:T})$ is not necessarily the same as the MPC cost
function (7a). In most cases $J(y_{1:T},u_{1:T})$ will be some function that
reflects engineering and economic goals, while $\tilde{\ell}(\cdot)$ can be
substantially simplified to facilitate online optimization.

Let
$\\{u^{\star}_{i|k}(\hat{s}_{k};\theta_{m},\theta_{o},\theta_{b},N_{p})\\}_{i=0}^{N_{p}}$
denote the optimal solution to (7) for a given state estimate and set of
tuning parameters. The receding-horizon control law, implicitly defined in
terms of the solution to (7), is given by

$\displaystyle\kappa_{c}(\hat{s}_{k};\theta_{m},\theta_{o},\theta_{b},N_{p})=u^{\star}_{0|k}(\hat{s}_{k};\theta_{m},\theta_{o},\theta_{b},N_{p}).$
(8)

Since the control-relevant state may not be directly measured, we assume a
parametrized state estimator of the following form exists:

$\displaystyle\hat{s}_{k+1}=\kappa_{e}(\hat{s}_{k},u_{k},y_{k+1};\theta_{e}),$
(9)

where $\theta_{e}\in\mathbb{R}^{n_{\theta,e}}$ denotes any free parameters in
the estimator. For example, $\theta_{e}$ could represent the process and
measurement noise matrices in an extended Kalman filter (EKF) (Hoshiya and
Saito, 1984). By combining the specified controller (8) and estimator (9)
structures with the system dynamics (1), we can represent the closed-loop
system in terms of the augmented state $z_{k}=(x_{k},\hat{s}_{k})$

$\displaystyle z_{k+1}$ $\displaystyle=f_{\text{cl}}(z_{k},w_{k};\theta),$
(10a) $\displaystyle u_{k}$ $\displaystyle=\kappa_{\text{cl}}(z_{k};\theta),$
(10b) $\displaystyle y_{k}$ $\displaystyle=h_{\text{cl}}(z_{k},v_{k}),$ (10c)

where
$f_{\text{cl}}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{w}}\times\mathbb{R}^{n_{\theta}}\to\mathbb{R}^{n_{z}}$
defines the autonomous closed-loop dynamics,
$\kappa_{\text{cl}}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{\theta}}\to\mathbb{R}^{n_{u}}$
is the control function in the augmented space,
$h_{\text{cl}}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{v}}\to\mathbb{R}^{n_{y}}$
is the measurement function in the augmented space,
$\theta=(\theta_{m},\theta_{o},\theta_{b},N_{p},\theta_{e})$ denotes the
complete set of tuning parameters appearing in the control structure, and
$n_{z}=n_{x}+n_{s}$ and
$n_{\theta}=n_{\theta,m}+n_{\theta,o}+n_{\theta,b}+n_{\theta,e}+1$ are the
augmented state and tuning parameter dimensions, respectively. Under this
restricted class of control policies, we can construct an approximation to (5)
as follows:

$\displaystyle\min_{\theta\in\Theta}$
$\displaystyle~{}~{}\mathbb{E}_{\delta}\\{J(y_{1:T},u_{1:T})\\},$ (11a) s.t.
$\displaystyle~{}~{}\text{\eqref{eq:closed-loop-
system}},~{}\mathbb{P}_{\delta}\\{y_{1:T}\in\mathbb{Y}_{1:T}\\}\geq
1-\epsilon.$ (11b)

Note that the control input constraints (3) are not included in this
formulation as they are directly enforced by the MPC law (7). The key
difference between (11) and (5) is that the proposed approximation (11)
optimizes over a finite dimensional space
$\Theta\subset\mathbb{R}^{n_{\theta}}$. The problem is also a single stage one
in which a full simulation can be carried out once the tuning parameters
$\theta$ are fixed. As such, we do not have to worry about adapting the
control policy over the time steps $k=1,\ldots,T$, as this is implicitly done
through the MPC law defined above. However, this problem is still not
immediately solvable due to the presence of the probabilistic operators. We
attempt to address this challenge next using a simulation optimization (SO)
paradigm (Amaran et al., 2016).

###### Remark 1

It is important to note that the number of model parameters $n_{\theta,m}$ in
(6) can be quite large in generic black-box models, which inherently increases
the complexity of the problem. As such, it is advised to encode as much prior
knowledge as possible in a particular problem at hand, which will typically
result in a gray-box model defined in a reduced parametric space.

###### Remark 2

Although we focused on a nominal MPC formulation above, this could
straightforwardly be replaced with robust or stochastic MPC methods that
directly account for uncertainty within the predictions. However, this will
come at the cost of more expensive closed-loop simulations as well as the need
to develop an uncertainty description. Again, one can imagine that some
parameters of the uncertainty distribution $\theta_{u}$ could also be treated
as controller tuning parameters.

###### Remark 3

Whenever the state of the model (6) is measurable then we can ignore the state
estimator (9). It is worthwhile noting, however, that the estimator could
still be useful whenever it is desired to adapt some of the model parameters
online due to unexpected drifts or faults.

##  4 Proposed Closed-loop Simulation Optimization Method

Let $f(\theta)=\mathbb{E}_{\delta}\\{J(y_{1:T},u_{1:T})\\}$ denote the
expected performance and
$c(\theta)=\mathbb{P}_{\delta}\\{y_{1:T}\in\mathbb{Y}_{1:T}\\}-1+\epsilon$ be
a shorthand for the chance constraint. We can now restate (11) in the
following compact manner:

$\displaystyle\min_{\theta\in\Theta}~{}f(\theta)~{}~{}\text{s.t.}~{}~{}c(\theta)\geq
0.$ (12)

As discussed in (Paulson and Mesbah, 2018), the uncertainty propagation steps
needed to evaluate $f(\theta)$ and $c(\theta)$ can be performed using a
variety of different techniques. Since the number of uncertainties
$n_{\delta}=n_{x}+n_{w}T+n_{v}(T+1)$ grows with the number of time steps $T$,
Monte Carlo (MC) sampling is a likely a good choice since its convergence rate
is known to be independent of the $n_{\delta}$ and instead only depends on the
number of samples. Using MC sampling, the objective and constraints in (12)
can then be approximated as (Kleywegt et al., 2002)

$\displaystyle f(\theta)$ $\displaystyle\approx
M^{-1}\textstyle\sum_{j=1}^{M}J(y_{1:T}^{j},u_{1:T}^{j}),$ (13a)
$\displaystyle c(\theta)$ $\displaystyle\approx
M^{-1}\textstyle\sum_{j=1}^{M}\big{(}\mathbf{1}_{\mathbb{Y}_{1:T}}(y_{1:T}^{j})-1+\epsilon\big{)},$
(13b)

where $M$ denotes the number of samples, $\mathbf{1}_{A}(x)$ denotes the
indicator function over the set $A$ (1 when $x\in A$ and 0 otherwise),
$\\{\delta^{1},\ldots,\delta^{M}\\}$ are independent and identically
distributed (i.i.d.) samples of $\delta\sim p(\delta)$, and
$(y_{1:T}^{j},u^{j}_{1:T})$ are the simulated closed-loop output and input
sequences given the $j$th uncertainty sample $\delta^{j}$. It is worth noting
that our only assumption about the system thus far is that we have the ability
to execute (10) for any choice of tuning parameters $\theta\in\Theta$, and
every simulation is performed under an i.i.d. sample of $\delta$.

Since a new set of uncertainty samples are drawn every time we run the closed-
loop “simulator” (10), the sample average approximations in (13) produce
stochastic/noisy observations of the objective and the constraints

$\displaystyle y^{f}$
$\displaystyle=f(\theta)+\varepsilon^{f},~{}~{}y^{c}=c(\theta)+\varepsilon^{c},$
(14)

where $\varepsilon^{f}$ and $\varepsilon^{c}$ represent the observation noise
for $f(\theta)$ and $c(\theta)$, respectively. For any $M$, these estimators
are unbiased, i.e.,
$\mathbb{E}_{\delta^{1},\ldots,\delta^{M}}\\{y^{f}\mid\theta\\}=f(\theta)$ and
$\mathbb{E}_{\delta^{1},\ldots,\delta^{M}}\\{y^{c}\mid\theta\\}=c(\theta)$;
however, these observations may have relatively high variance unless $M$ is
large. We assume it is not possible to select a large $M$ due to the fact that
the closed-loop simulation or experiment is prohibitively costly (either from
a computational or monetary cost point-of-view). Instead these functions can
only be evaluated on the order of 100 times or less, meaning $M$ cannot be
large enough to ignore noise in the observations. As such, we want to utilize
an algorithm that can systematically explore the tuning parameter space
(relative to random or grid-based search methods) that can also accommodate
noisy objective and constraint observations.

Bayesian optimization (BO) is a family of algorithms that can solve black-box
optimization problems in the presence of noisy observations. The basic idea in
BO is to construct a surrogate model of the objective function $f(\theta)$
using a set of initial $n$ observations denoted by
$\mathcal{D}_{n}=\\{(\theta_{1:n},y^{f}_{1:n})\\}$. The statistical surrogate
model thus provides a posterior distribution of the objective function that
can be combined with an acquisition function to decide where to sample next.
The acquisition function is defined in a way to tradeoff between uncertainty
(related to variance) and performance (related to mean) at unexplored points
$\theta\in\Theta$. The new observation is added to the data set, i.e.,
$\mathcal{D}_{n+1}=\\{\mathcal{D}_{n},(\theta_{n+1},y^{f}_{n+1})\\}$, which is
then used to update the surrogate model. This process continues until the
surrogate model converges to the global solution or the maximum number of
iterations is reached. See, e.g., Shahriari et al. (2015) for a recent review
of BO and Lu et al. (2020) for an application to MPC tuning in central
heating, ventilation, and air conditioning (HVAC) plants.

First, we discuss the traditional BO strategy for (12) that neglects
constraints $c(\theta)\geq 0$ and then describe an extension that modifies the
acquisition function to account for the probability that constraints are
satisfied. In this work, we exclusively use Gaussian process (GP) surrogate
models (Rasmussen, 2003). In particular, we assume that the objective function
has a GP prior of the form $f(\theta)\sim\mathcal{GP}(\mu_{0},k)$ where
$\mu_{0}:\Theta\to\mathbb{R}$ is the prior mean and
$k:\Theta\times\Theta\to\mathbb{R}$ is the prior covariance function. There
are many possible choices for the covariance function; in this work we
selected the squared exponential function

$\displaystyle
k(\theta,\theta^{\prime})=\exp\left(-\frac{1}{2l^{2}}\|\theta-\theta^{\prime}\|^{2}\right),$
(15)

where $l$ is a length-scale parameter. Under the GP prior and assumed i.i.d.
Gaussian noise $\varepsilon^{f}\sim\mathcal{N}(0,\sigma^{2})$, the function
evaluations $y^{f}_{1:n}$ are jointly Gaussian:

$\displaystyle
y^{f}_{1:n}\sim\mathcal{N}(\mathbf{m},\mathbf{K}+\sigma^{2}I_{n}),$ (16)

where $[\mathbf{m}]_{i}=\mu_{0}(\theta_{i})$ are the elements of the mean
vector and $\mathbf{K}\in\mathbb{R}^{n\times n}$ is a symmetric covariance
matrix with elements $[\mathbf{K}]_{i,j}=k(\theta_{i},\theta_{j})$. This
implies the corresponding function value $f(\theta)$ at any test point
$\theta$ must be jointly Gaussian with $y^{f}_{1:n}$, i.e.,

$\displaystyle\begin{bmatrix}y^{f}_{1:n}\\\
f(\theta)\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{m}\\\
\mu_{0}(\theta)\end{bmatrix},\begin{bmatrix}\mathbf{K}+\sigma^{2}I_{n}&\mathbf{k}(\theta)\\\
\mathbf{k}(\theta)^{\top}&k(\theta,\theta)\end{bmatrix}\right).$ (17)

where $\mathbf{k}(\theta)\in\mathbb{R}^{n}$ is a vector of covariance terms
between $\theta_{1:n}$ and $\theta$, i.e.,
$[\mathbf{k}(\theta)]_{i}=k(\theta_{i},\theta)$. Due to the properties of
joint Gaussian random variables, we find that the posterior distribution
$p(f(\theta)|y^{f}_{1:n},\theta_{1:n},\theta)$ of the objective given all
available noisy observations is Gaussian with the following mean and
covariance (Rasmussen, 2003):

$\displaystyle\mu_{n}(\theta)$
$\displaystyle=\mu_{0}(\theta)+\mathbf{k}(\theta)^{\top}(\mathbf{K}+\sigma^{2}I_{n})(y^{f}_{1:n}-\mathbf{m}),$
(18a) $\displaystyle\sigma^{2}_{n}(\theta)$
$\displaystyle=k(\theta,\theta)-\mathbf{k}(\theta)^{\top}(\mathbf{K}+\sigma^{2}I_{n})^{-1}\mathbf{k}(\theta).$
(18b)

Given this posterior distribution, we still need to optimize over an
acquisition function $\alpha(\theta;\mathcal{D}_{n})$ to search for the next
optimal sampling point $\theta_{n+1}$. Although several are available, we
focus on the expected improvement (EI) criteria that measures the expected
amount by which the current objective will be improved over some incumbent (or
current best) solution. The EI acquisition function has been shown to be
analytically computable for a GP model

$\displaystyle\alpha_{\text{EI}}(\theta;\mathcal{D}_{n})=\sigma_{n}(\theta)(z_{n}(\theta)\Phi(z_{n}(\theta))+\phi(z_{n}(\theta)),$
(19)

where $z_{n}(\theta)=(\eta-\mu_{n}(\theta))/\sigma_{n}(\theta)$, $\Phi$ and
$\phi$ are, respectively, the standard Gaussian cumulative density and
probability density functions, and $\eta$ is the incumbent solution. The
“best” choice of $\eta$ depends on the context; for deterministic objective
evaluations, it is often set to the best observed value
$\eta=\max_{i\in\\{1,\ldots,n\\}}y_{i}^{f}$. However, as discussed in (Wang
and de Freitas, 2014), this choice can be quite fragile in the stochastic
setting and thus we use the best mean value
$\eta=\min_{\theta\in\Theta}\mu_{n}(\theta)$ as a more reasonable alternative.

To handle the black-box constraints, we need to modify the acquisition
function to show improvement only when $c(\theta)\geq 0$ holds. Similarly to
the objective, we model the constraint function $c(\theta)$ with a GP prior
whose evaluations are corrupted with Gaussian noise. We must then weight the
original EI in (19) by the probability of the constraints being satisfied.
This results in the expected improvement with constraints (EIC) that can be
analytically computed as follows (Gardner et al., 2014):

$\displaystyle\alpha_{\text{EIC}}(\theta;\mathcal{D}_{n})$
$\displaystyle=\alpha_{\text{EI}}(\theta;\mathcal{D}_{n})\mathbb{P}\\{c(\theta)\geq
0|\mathcal{D}_{n}^{c}\\},$ (20)
$\displaystyle=\alpha_{\text{EI}}(\theta;\mathcal{D}_{n})\Phi\left(\frac{\mu_{n}^{c}(\theta)}{\sigma_{n}^{c}(\theta)}\right),$

where $\mathcal{D}_{n}^{c}=\\{(\theta_{1:n},y^{c}_{1:n})\\}$ is the set of $n$
noisy constraint observations and $\mu^{c}_{n}(\theta)$ and
$\sigma_{n}^{c}(\theta)$ are, respectively, the posterior predictive mean and
standard deviation of the GP surrogate model for $c(\theta)$ – similarly
defined to (18) for the objective. Due to the presence of constraints, we must
modify the definition of the incumbent value $\eta$ to be the best mean value
such that the constraints are satisfied. Since we cannot guarantee the
constraints will be satisfied given only noisy observations, it has been
suggested to require constraint satisfaction with some probability:

$\displaystyle\eta=\min_{\theta\in\Theta}~{}\mu_{n}(\theta)~{}~{}\text{s.t.}~{}~{}\Phi\left(\frac{\mu_{n}^{c}(\theta)}{\sigma_{n}^{c}(\theta)}\right)\geq
1-\beta,$ (21)

for a relatively small $\beta\in(0,1)$. However, we still cannot guarantee a
feasible solution to (21) exists due to, e.g., limited data available for the
constraints. To overcome this challenge, one can ignore the factor
$\alpha_{\text{EI}}(\theta;\mathcal{D}_{n})$ in (20) to search only for a
feasible point by maximizing the probability of constraint satisfaction
(Gelbart et al., 2014).

In our proposed strategy, the next sampling point $\theta_{n+1}$ is obtained
by solving the following optimization problem:

$\displaystyle\theta_{n+1}=\operatorname*{argmax}_{\theta\in\Theta}~{}\alpha_{\text{EIC}}(\theta;\mathcal{D}_{n}).$
(22)

Since the acquisition function is cheap to evaluate (relative to $f$ and $c$),
its maximization can be carried out more efficiently. An overview of the
various methods used to solve this optimization that can avoid local solutions
is provided in (Gardner et al., 2014). An illustration of the proposed
automated MPC tuning strategy under uncertainty using constrained BO is shown
in Figure 1.

###### Remark 4

There are several hyperparameters in the GP regression models that may appear
in the mean function, covariance kernel (e.g., $l$), or the observation model
(e.g., noise level $\sigma^{2}$). There are two main ways for handling these
hyperparameters. The most involved way is to marginalize out the uncertainty
in these hyperparameters by determining their posterior distribution given the
data $\mathcal{D}_{n}$ using Bayes’ rule. The acquisition function must then
be integrated over this posterior, which can be approximated using some
sampling strategy. A simpler approach is to approximate this integral in terms
of a Dirac delta measure for the posterior at the maximum likelihood or
maximum a posteriori point estimate. These point estimates can be found by
either maximizing the marginal likelihood or the unnormalized posterior for
which standard optimization methods are readily available.

###### Remark 5

The assumed Gaussian noise model with constant variance for the objective and
constraint evaluations in (14) is an approximation made for simplicity. In
reality, the noise level may depend on the specific value of $\theta$, which
should be accounted for to obtain a more accurate GP surrogate model. One way
to do this is to model the noise with a second GP model, as discussed in
(Goldberg et al., 1998). Future work is needed to best incorporate input-
dependent noise models into constrained BO.

Figure 1: Illustration of how constrained BO method can be used for data-
driven automatic MPC tuning under constraints. The blue box denotes the
closed-loop experiment or simulation that must be performed 1 or more times
for fixed tuning parameters. The red box shows the constrained BO steps, which
must be executed to decide the next best tuning parameter values to use for
closed-loop data collection.

##  5 Numerical Example

###  5.1 Black-box plant model description

We demonstrate the proposed approach on a nonlinear continuously stirred tank
reactor (CSTR) problem adapted from (Paulson and Mesbah, 2018). The following
set of reactions are assumed to take place in the CSTR

A $\displaystyle\to\text{B}\to\text{C},~{}~{}2\text{A}\to\text{D}.$

The plant model (1) is assumed to be a discretized version of the following
set of ordinary differential equations

$\displaystyle\dot{c}_{\mathrm{A}}$
$\displaystyle=F(c_{\mathrm{A}0}-c_{\mathrm{A}})-k_{1}c_{\mathrm{A}}-k_{3}c_{\mathrm{A}}^{2},$
(23a) $\displaystyle\dot{c}_{\mathrm{B}}$
$\displaystyle=-Fc_{\mathrm{B}}+k_{1}c_{\mathrm{A}}-k_{2}c_{\mathrm{B}},$
(23b) $\displaystyle\dot{T}_{\mathrm{R}}$
$\displaystyle=F(T_{\text{in}}-T_{\mathrm{R}})+\frac{k_{\mathrm{W}}A}{\rho
c_{\mathrm{p}}V_{\mathrm{R}}}(T_{\mathrm{K}}-T_{\mathrm{R}})$ (23c)
$\displaystyle~{}~{}~{}-\frac{k_{1}c_{\mathrm{A}}\Delta
H_{\mathrm{AB}}+k_{2}c_{\mathrm{B}}\Delta
H_{\mathrm{BC}}+k_{3}c_{\mathrm{A}}^{2}\Delta H_{\mathrm{AD}}}{\rho
c_{\mathrm{p}}},$ $\displaystyle\dot{T}_{\mathrm{K}}$
$\displaystyle=\frac{1}{m_{\mathrm{K}}c_{\mathrm{pK}}}(\dot{Q}_{\mathrm{K}}+k_{\mathrm{W}}A(T_{\mathrm{R}}-T_{\mathrm{K}})).$
(23d)

where $c_{\mathrm{A}}$ and $c_{\mathrm{B}}$ denote the concentrations of
species A and B, respectively, $T_{\mathrm{R}}$ is the reactor temperature,
and $T_{\mathrm{K}}$ is the coolant temperature. The control input is the feed
flowrate $F=\dot{V}_{in}/V_{R}$ and the reaction kinetics exhibit the
following Arrhenius dependence:

$\displaystyle
k_{i}=k_{0,i}e^{\frac{-E_{A,i}}{R\left(T_{R}+273.15\right)}},~{}~{}\forall
i\in\\{1,2,3\\}.$ (24)

The sampling time was chosen to be $\delta t=0.005$ h. We assume that only
species B and the reactor temperature are measurable online such that
$y_{k}=(c_{\mathrm{B},k},T_{\mathrm{R},k})+v_{k}$ where
$v_{k}\sim\mathcal{N}(0,\Sigma_{v})$ is noise, with
$\Sigma_{v}=\text{diag}(\sigma_{B}^{2},\sigma_{R}^{2})$. For simplicity, we
did not consider process noise. The model parameters can be found in Table 1
and the initial conditions and relevant constraints are summarized in Table 2.
Since only noisy temperature measurements are available, the bounds on
$T_{\mathrm{R}}$ are enforced as a chance constraint. Instead of the single
constraint (4), we enforce constraints separately at each time point, i.e.,
$\mathbb{P}_{\delta}\\{y_{k}\in\mathbb{Y}_{k}\\}\geq 1-\epsilon_{k}$ where
$\epsilon_{k}=0.05$ for all $k=1,\ldots,T$. This does not change the
derivations shown in Section 4 because we can always define
$c(\theta)=\min_{k\in\\{1,\ldots,T\\}}\\{\mathbb{P}_{\delta}\\{y_{k}\in\mathbb{Y}_{k}\\}-1+\epsilon_{k}\\}\geq
0$ in terms of the smallest probability over time. We assume that our closed-
loop objective is to maximize the production of component B

$\displaystyle
n_{B,Total}=\int_{0}^{t_{f}}\dot{n}_{B}(t)dt\approx\sum_{k=0}^{T-1}V_{\text{in}}c_{\mathrm{B},k}F_{k}\delta
t,$ (25)

such that the performance cost is $J(y_{1:T},u_{1:T})=-n_{B,Total}$ where
$t_{f}=0.2$ h is the final time of the experiment and $T=40$ sampling time
intervals.

Table 1: Parameter values for CSTR case study. Parameter | Value | Unit
---|---|---
$k_{0,1}$ | $1.287~{}10^{12}$ | h-1
$k_{0,2}$ | $1.287~{}10^{12}$ | h-1
$k_{0,3}$ | $9.043~{}10^{9}$ | L mol-1 h-1
$E_{A,1}/R$ | $9758.3$ | K
$E_{A,2}/R$ | $9758.3$ | K
$E_{A,3}/R$ | $7704.0$ | K
$\Delta H_{\mathrm{AB}}$ | $4.2$ | kJ mol-1
$\Delta H_{\mathrm{BC}}$ | $-11.0$ | kJ mol-1
$\Delta H_{\mathrm{AD}}$ | $-41.85$ | kJ mol-1
$\rho$ | $0.9342$ | kg L-1
$c_{\mathrm{p}}$ | $3.01$ | kJ kg-1 K-1
$c_{\mathrm{pK}}$ | $2.0$ | kJ kg-1 K-1
$A$ | $0.215$ | m2
$V_{\mathrm{R}}$ | $10.01$ | L
$m_{\mathrm{k}}$ | $5.0$ | kg
$T_{\mathrm{in}}$ | $130.0$ | ${}^{\circ}\text{C}$
$k_{\mathrm{W}}$ | $4032$ | kJ h-1 m-2 K-1
$\dot{Q}_{\mathrm{K}}$ | $-4500$ | kJ h-1
Table 2: Initial conditions and constraints for CSTR. Variable | Init. cond. | Min. | Max. | Unit
---|---|---|---|---
$c_{\text{A}}$ | 1.0 | – | – | mol/L
$c_{\text{B}}$ | 1.0 | – | – | mol/L
$T_{\text{R}}$ | 100.0 | 100.0 | 150.0 | ${}^{\circ}\text{C}$
$T_{\text{K}}$ | 100.0 | – | – | ${}^{\circ}\text{C}$
$F$ | – | 5.0 | 35.0 | h-1

###  5.2 Parametrized formulation of MPC policy

We emphasize the fact that the plant model (23)–(24) is a black-box in this
work and so cannot be used in model-based control design. Instead, we look to
learn a control-relevant model of the form (6). In particular, we focus on an
NARX model structure that can be stated as follows

$\displaystyle
y_{k}=f_{\text{NARX}}(y_{k-1},,\ldots,y_{k-L_{y}},u_{k-1},\ldots,u_{k-L_{u}}),$
(26)

where $L_{y}$ and $L_{u}$ denotes the number of output and input lags,
respectively. The function $f_{\text{NARX}}$ can be any nonlinear function
such as a polynomial or a neural network. In this work, we selected a second-
order polynomial that ignores the interaction terms. Although in principle
$L_{y}$ and $L_{u}$ could be treated as model parameters $\theta_{m}$, we
decided to fix them to $L_{y}=L_{u}=1$ to ensure the tuning parameter space
was not too large. This implies that $\theta_{m}$ is composed of 14
coefficients – 7 for each $n_{y}=2$ output that multiplies the basis set
$\\{1,[y_{k}]_{1},[y_{k}]_{2},u_{k},[y_{k}]_{1}^{2},[y_{k}]_{2}^{2},u_{k}^{2}\\}$.
We fixed the MPC prediction horizon to be $N_{p}=10$, as performance was found
to be relatively insensitive to this choice. We selected an economic MPC
objective function related to the moles of B produced, i.e.,
$\tilde{\ell}(y_{i|k},u_{i|k})=V_{\text{in}}c_{\mathrm{B},i|k}F_{i|k}\delta
t$. The reactor temperature constraint is active in this problem such that it
can easily be violated whenever the prediction model is not very accurate. To
address this, we include a backoff parameter $\theta_{b}$ that we wish to
select using the proposed constrained BO method, which results in a total of
$n_{\theta}=14+1=15$ tuning parameters to be optimized.

###  5.3 Automated MPC tuning using constrained BO and comparison to
sequential model identification

Given the black-box plant description and the parametrized MPC policy, we can
now solve (12) using the constrained BO method summarized in Section 4; the
bayesopt function in Matlab can be used to handle the GP construction of the
objective and constraints as well as specify the EIC criteria using the
“coupled constraint” option. We explicitly set the objective and constraint
evaluations to being stochastic such that the variance is included as a
hyperparameter in the corresponding GP models. We budgeted a total of
$N_{\text{iter}}=40$ iterations, with the first 5 being randomly chosen within
the assumed parameter space $\Theta=[-2,2]^{14}\times[0,0.2]$, i.e., all NARX
parameters are bounded by $[-2.2]$ and the backoff
$\theta_{b}\in[0,0.2]$111The NARX model is constructed using scaled data, so
that all inputs and outputs are bounded within $[0,1]$..

Figure 2: The maximum performance (moles of B) versus number of BO iterations
for 100 replicate runs. The dark blue line represents the mean value, the
shaded region is +/- one standard deviation, and the error bars show the
minimum and maximum values.

To focus on constraint handling first, we initially set the measurement noise
variances to zero. The performance evolution, measured in terms of total moles
of B produced, over the number of iterations is shown in Figure 2 for 100
independent runs of the constrained BO algorithm. We see a consistent
improvement in performance up until around 30 iterations – we also see that
the variance in the final result reduces as number of iterations increases.
The resulting closed-loop temperature and feed flowrate profiles are shown in
Figure 3, which tightly satisfy constraints. To better contextualize these
results, we compare this automated model learning procedure to the traditional
open-loop NARX identification paradigm (see red dash-dotted lines in Figure
3). In particular, we performed pseudo random binary step tests on the plant
every 100 time steps to collect a total of 3000 input-output data points. This
open-loop data was used to train an NARX model (26) using Matlab’s System
Identification Toolbox, which achieved high accuracy on a holdout set of data,
i.e., prediction accuracy of 91.5%, which is equal to $1-(\text{normalized
root mean squared error})$.

(a) temperature profile

(b) feed flowrate profile

Figure 3: Closed-loop output and input profiles for the best solution obtained
with BO compared to an open-loop identified NARX model for the noise-free
case.

Even though the prediction accuracy is high, the resulting closed-loop
performance is poor in the sense that we observe significant variability in
the control input profile and minor temperature constraint violations. In
fact, the moles of B produced was 5.62, which is less than half of the 12.17
obtained with the constrained BO method. This is likely due to the fact that
the input-output data was not sufficiently informative in the regime of high-
performance control. As such, we can interpret the proposed strategy as a way
to select control-oriented prediction models that are suited to the task at
hand – this is important because the desired set of closed-loop trajectories
usually represents a much smaller slice of the input-output space than that
used in open-loop identification techniques.

Lastly, to demonstrate that the approach is capable of operating in the
presence of uncertainty, we repeated the analysis above with measurement noise
variances $\sigma_{B}=0.2$ mol and $\sigma_{R}=10^{\circ}\text{C}$. To keep a
small experimental budget, we selected $M=1$. Due to space limitations, we
only show the resulting closed-loop temperature profiles for the open-loop
NARX and constrained BO identified models in Figure 4. It is interesting to
note that the open-loop identified NARX model appears to be quite susceptible
to overfitting, which leads to much higher variability and constraint
violations than in the noise-free case. On the other hand, the constrained BO
method is able to select NARX coefficients that produce temperature profiles
with significantly lower variance and, through proper tuning of the backoff,
completely avoids constraint violations. We also emphasize that these results
were obtained using the same number of runs as in the deterministic case (only
1 experiment performed for a random noise sequence). This highlights the value
of using a GP model that can explicitly account for noisy objective and
constraint violations, as it helps guide the search in a way that is not
overly optimistic about any one observation.

Figure 4: Closed-loop temperature profile obtained using the proposed
constrained BO method compared to an open-loop identified NARX model in the
presence of measurement noise for 100 Monte Carlo realizations.

##  6 Conclusions and Future Work

We present an automated performance-driven MPC tuning strategy for black-box
systems in the presence of uncertainty. Since the MPC tuning parameters (which
includes the prediction model itself in our approach) affect the closed-loop
performance and constraints in complex (generally non-convex and non-smooth)
ways, we rely on derivative-free simulation optimization methods to automate
the overall tuning process. In particular, we rely on a constrained variant of
Bayesian optimization (BO) that utilizes state-of-the-art surrogate models for
the objective and constraint functions constructed from noisy observations. By
combining these surrogates with a so-called acquisition function, we can
iteratively determine the best tuning parameter combination to use for our
next closed-loop experiment that directly tradeoffs exploitation of the
current best tuning parameters with exploitation of the feasible tuning
parameter space (while also accounting for the probability of constraint
satisfaction). We demonstrated the effectiveness of the approach on a
nonlinear continuously stirred tank reactor case study wherein we observed a
more than two-fold improvement in performance compared to the traditional
open-loop model identification methods. There are several interesting
directions for future work including development of novel constrained BO
methods that can more readily scale to high-dimensional problems and that can
account for noise levels in the function evaluations that depend on the
specific tuning parameter values.

## References

  * Amaran et al. (2016) Amaran, S., Sahinidis, N.V., Sharda, B., and Bury, S.J. (2016).  Simulation optimization: A review of algorithms and applications.  _Annals of Operations Research_ , 240, 351–380.
  * Bansal et al. (2017) Bansal, S., Calandra, R., Xiao, T., Levine, S., and Tomiin, C.J. (2017).  Goal-driven dynamics learning via Bayesian optimization.  In _Proceedings of the IEEE 56th Annual Conference on Decision and Control_ , 5168–5173. Melbourne.
  * Berkenkamp et al. (2016) Berkenkamp, F., Schoellig, A.P., and Krause, A. (2016).  Safe controller optimization for quadrotors with Gaussian processes.  In _Proceedings of the IEEE International Conference on Robotics and Automation_ , 491–496. Stockholm.
  * Forgione et al. (2019) Forgione, M., Piga, D., and Bemporad, A. (2019).  Efficient calibration of embedded MPC.  _arXiv preprint arXiv:1911.13021_.
  * Gardner et al. (2014) Gardner, J.R., Kusner, M.J., Xu, Z.E., Weinberger, K.Q., and Cunningham, J.P. (2014).  Bayesian optimization with inequality constraints.  In _Proceedings of the International Conference on Machine Learning_ , 937–945. Beijing.
  * Garriga and Soroush (2010) Garriga, J.L. and Soroush, M. (2010).  Model predictive control tuning methods: A review.  _Industrial & Engineering Chemistry Research_, 49(8), 3505–3515.
  * Gelbart et al. (2014) Gelbart, M.A., Snoek, J., and Adams, R.P. (2014).  Bayesian optimization with unknown constraints.  In _Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence_ , 250–259.
  * Gevers (2005) Gevers, M. (2005).  Identification for control: From the early achievements to the revival of experiment design.  _European Journal of Control_ , 11, 335–352.
  * Goldberg et al. (1998) Goldberg, P.W., Williams, C.K., and Bishop, C.M. (1998).  Regression with input-dependent noise: A Gaussian process treatment.  In _Advances in neural information processing systems_ , 493–499.
  * Hernández-Lobato et al. (2016) Hernández-Lobato, J.M., Gelbart, M.A., Adams, R.P., Hoffman, M.W., and Ghahramani, Z. (2016).  A general framework for constrained Bayesian optimization using information-based search.  _Journal of Machine Learning Research_ , 17(160), 1–53.
  * Hoshiya and Saito (1984) Hoshiya, M. and Saito, E. (1984).  Structural identification by extended Kalman filter.  _Journal of Engineering Mechanics_ , 110, 1757–1770.
  * Khosravi et al. (2020) Khosravi, M., Behrunani, V., Myszkorowski, P., Smith, R.S., Rupenyan, A., and Lygeros, J. (2020).  Performance-driven cascade controller tuning with Bayesian optimization.  _arXiv preprint arXiv:2007.12536_.
  * Kleywegt et al. (2002) Kleywegt, A.J., Shapiro, A., and Homem-de Mello, T. (2002).  The sample average approximation method for stochastic discrete optimization.  _SIAM Journal on Optimization_ , 12, 479–502.
  * Lu et al. (2020) Lu, Q., Kumar, R., and Zavala, V.M. (2020).  MPC controller tuning using Bayesian optimization techniques.  _arXiv preprint arXiv:2009.14175_.
  * Neumann-Brosig et al. (2019) Neumann-Brosig, M., Marco, A., Schwarzmann, D., and Trimpe, S. (2019).  Data-efficient autotuning with Bayesian optimization: An industrial control study.  _IEEE Transactions on Control Systems Technology_ , 28(3), 730–740.
  * Paulson and Mesbah (2018) Paulson, J.A. and Mesbah, A. (2018).  Nonlinear model predictive control with explicit backoffs for stochastic systems under arbitrary uncertainty.  _IFAC-PapersOnLine_ , 51(20), 523–534.
  * Piga et al. (2019) Piga, D., Forgione, M., Formentin, S., and Bemporad, A. (2019).  Performance-oriented model learning for data-driven MPC design.  _IEEE Control Systems Letters_ , 3(3), 577–582.
  * Rasmussen (2003) Rasmussen, C.E. (2003).  Gaussian processes in machine learning.  In _Summer School on Machine Learning_ , 63–71. Springer.
  * Rawlings and Mayne (2009) Rawlings, J.B. and Mayne, D.Q. (2009).  _Model predictive control: Theory and design_.  Nob Hill Publishing, Madison, Wisconsin.
  * Shahriari et al. (2015) Shahriari, B., Swersky, K., Wang, Z., Adams, R.P., and De Freitas, N. (2015).  Taking the human out of the loop: A review of Bayesian optimization.  _Proceedings of the IEEE_ , 104(1), 148–175.
  * Snoek et al. (2012) Snoek, J., Larochelle, H., and Adams, R.P. (2012).  Practical Bayesian optimization of machine learning algorithms.  In _Advances in Neural Information Processing Systems_ , 2951–2959.
  * Wang and de Freitas (2014) Wang, Z. and de Freitas, N. (2014).  Theoretical analysis of Bayesian optimisation with unknown Gaussian process hyperparameters.  _arXiv preprint arXiv:1406.7758_.