# Coupling and selecting constraints in Bayesian optimization under
uncertainties

Julien Pelamatti  EDF R&D, France (formerly Ecole des Mines de Saint-Etienne)
Rodolphe Le Riche  CNRS LIMOS (Mines Saint-Etienne and UCA), France  Céline
Helbert  Université de Lyon - CNRS, UMR 2508, Institut Camille Jordan - Ecole
Centrale de Lyon, France  Christophette Blanchet-Scalliet  Université de Lyon
- CNRS, UMR 2508, Institut Camille Jordan - Ecole Centrale de Lyon, France

###### Abstract

We consider chance constrained optimization where it is sought to optimize a
function while complying with constraints, both of which are affected by
uncertainties. The high computational cost of realistic simulations strongly
limits the number of evaluations and makes this type of problems particularly
challenging. In such a context, it is common to rely on Bayesian optimization
algorithms. Assuming, without loss of generality, that the uncertainty comes
from some of the inputs, it becomes possible to build a Gaussian process model
in the joint space of design and uncertain variables. A two-step acquisition
function is then used to provide, both, promising optimization variables
associated to relevant uncertain samples. Our overall contribution is to
correlate the constraints in the GP model and exploit this to optimally
decide, at each iteration, which constraint should be evaluated and at which
point. The coupled Gaussian model of the constraints relies on an output-as-
input encoding. The constraint selection idea is developed by enabling that
each constraint can be evaluated for a different uncertain input, thus
improving the refinement efficiency. Constraints coupling and selection are
gradually implement in 3 algorithm variants which are compared to a reference
Bayesian approach. The results are promising in terms of convergence speed,
accuracy and stability as observed on a 2, a 4 and a 27-dimensional problems.

##  1 Introduction

Models of realistic physical or engineering systems are frequently impacted by
uncertain random parameters that represent model error, uncontrolled
environmental conditions, or numerical noise such as imperfectly converged
internal routines. When optimizing these systems, it is necessary to take into
account the presence of the random parameters to determine a sufficiently
robust solution in terms of both the objective and the constraint functions.
Depending on the application, different approaches to handle the uncertainties
within an optimization problem exist.

A distinction can be made between the stochastic and the robust approaches. In
the stochastic methods, the uncertainty in the optimization criteria
(objective and constraint functions) is modeled as a random noise corrupting
the observed values. Stochastic algorithms can, under conditions about their
own stochasticity, accomodate the noise in the observations and converge to
the problem solutions: this has given rise in the early 50’s to stochastic
descent algorithms [1, 2] that have since then experienced great developments
[3, 4] often in relation to machine learning [5]. Well-performing versions of
(stochastic) evolutionary algorithms for noisy functions have been identified
[6, 7] thanks to competitions on reference benchmarks [8].

Other approaches to optimization under uncertainty have emerged in many
communities (operations research [9], finance [10], engineering [11], machine
learning [12] and many others). They quantify with deterministic functions the
trade-off between the performance and the risk associated to the
uncertainties. In robust optimization, the set of uncertain parameters is
typically handled through worst-case scenarii or convex approximations of it
[9, 13]. In a variety of situations where the set of uncertainties is large
(possibly infinite) and the problem structure non-convex, over-conservatism of
the solutions and the computational tractability of the methods are inherent
difficulties of this family of problems. When the uncertainties are described
by probability distributions, it is now customary to optimize a risk measure
which takes into account the whole probability distribution.

In the presence of optimization constraints, the uncertainties affect the
feasibility of the solutions. Probabilistic models of the constraints are
called chance constraints [14] or reliability constraints [15]. The handling
of chance constraints, mostly in engineering, has spawned the field of
Reliability-Based Design Optimization (RBDO). In RBDO, a popular family of
approaches appproximate the probability of constraints satisfaction by
reliability indices which are computationally easier to calculate and suitable
to a coupling within an optimization problem (First/Second Order Reliability
Methods, [15]). However, when working with reliability indices, the failure
probability is often poorly estimated. Other RBDO techniques approximate the
failure probability with subset and importance sampling in conjunction with
metamodels which allow the inclusion of the chance constraint within an
optimization loop [16, 17].

In this article, the problem of minimizing the mean of a stochastic function
(also known as Bayes risk) under chance constraints (i.e., reliability
constraints) is addressed. It is assumed that the problem functions are costly
(i.e., cannot be calculated more than $\mathcal{O}(100)$ times) so that one
seeks to find the neighborhood of the optimal solution with as few function
evaluations as possible. Examples of expensive functions are nonlinear Finite
Element Models used for structural analyses, Computational Fluid-Dynamics
simulations and most real life experiments. The cost limitation rules out the
possibility of relying on most stochastic algorithms (like stochastic
gradients, evolution strategies) which require a large number of evaluations.
It is also assumed that there is no information regarding the convexity or the
gradient of the problem functions. These assumptions limit the type of
algorithms which are relevant to the class of surrogate-based methods, among
which Bayesian Optimization (BO) [18, 19] will be our focus.

The principle underlying BO is to learn the problem functions by Gaussian
processes (GPs) [20] which are computationally less expensive than the
functions. An acquisition criterion, based on the GPs and adapted to the
problem formulation, defines the most promising points of the design space.
Its repeated optimization, followed by an update of the GPs, constitutes the
BO algorithm.

Only recently has the handling of chance constraints been treated with BO [16,
17, 21, 22] (cf. Section 2 for a detailed discussion of past related works).
This article presents three contributions to BO in the situation where several
inequalities must be satisfied in probability. Similarly to [21, 23], the
starting point is a procedure which first maximizes the feasible improvement
and then reduces its variance in order to define new values of both the
optimization and the random variables.

The first contribution of this article is to account for the couplings between
the constraints through an “output as input” formulation. This is important as
the constraints are usually linked to each other, e.g., when they express
different criteria extracted from the same physical phenomenon, or when they
are competing anticorrelated criteria (such as price vs. availability of raw
materials).

The second contribution concerns the common situation where the objective and
the constraints are calculated with different codes. We will show how the
objective and the constraints can be refined independently in the random
search space. Indeed, the random variables that provide the largest amount of
information regarding i) the expected objective function and ii) the expected
feasibility of a given candidate solution are different.

The third contribution is to allow the optimization algorithm to select, at
each iteration, the constraint that provides the most information about the
location of the robust optimum. The motivation is that, when the constraints
can be computed independently, it is useless to refine a constraint for which
the (in)feasibility of the candidate solution is known with a sufficiently
high certainty.

The article continues with a more detailed review of past related works
followed, in Section 3, by a formulation of the problem. Section 4 describes
the output as input model and the BO acquisition criteria. The independent
calculation of the random variables and the selection of the constraint to be
refined are presented in Section 5. The last parts of the article are
dedictated to the implementation of the algorithms (Section 6.2) and numerical
tests (Section 7) on 2 analytical and an industrial problems. To assist the
reader, the main symbols and acronyms used in the article are summarized in
Appendix A.

##  2 Past related works

Our study rests on preliminary works concerning multi-output Gaussian
processes, Bayesian optimization with constraints and Bayesian optimization
with uncertainties. We now review this bibliography in turn.

###  2.1 Multi-output Gaussian Processes

A contribution of this article is the modeling of chance constraints with a
single multi-output Gaussian process rather than with independent models. The
idea is to take into account the correlation between outputs in order to
improve the modeling accuracy thus reducing the number of computer runs
required to converge to the problem optimum.

The model we introduce relies on the “output as input” principle where a
nominal input is added to the variables list to designate the selected
constraint. Other approaches found in the litterature can also model
correlations between outputs. A first parallel can be drawn with surrogates
for multi-fidelity modeling [24], [25]. These methods represent with a single
vector valued function the responses of a physical system associated to
different precisions (_e.g._ , varying mesh densities in finite elements).

Most of the so-called co-kriging approaches ([26], [27], [28]) define a
hierarchical relation between the various fidelity levels: the surrogate model
for the fidelity level $t$ is recursively defined as the Gaussian process
model of the previous fidelity level ($t-1$) summed with a second Gaussian
process, modeling the additional information provided by the next fidelity
level. However, such a hierarchical dependence hypothesis doesn’t always
stand. For example, when the components of the vectorial output represent
different physical quantities, it doesn’t necessarily make sense to model them
recursively. In this case, the order of the various outputs would be an
arbitrary modeling decision.

Alternative, non hierarchical GPs for vector valued functions include models
based on co-regionalization matrices between the outputs [29], [30], [31].
Recently, deep Gaussian processes have also been used to model multi-fidelity
phenomena, allowing to better capture non-stationary behaviors ([32], [33]).

The multi-output modeling methods described above have been applied to
computationally intensive optimization problems. These include multi-objective
optimization ([34], [35]) and multi-fidelity design optimization ([36], [37],
[38]). The ability to select which constraint should be evaluated at each
iteration in order to maximise the information gained is a salient feature of
the current article. A somewhat comparable idea was proposed in [39] to deal
with multi-model design. A specific sub-model (and its intervening variables)
was selected for evaluation based on a ratio of the “knowledge gradient” (the
predicted one-step-ahead gain in performance) to the submodel cost.

###  2.2 Constrained Bayesian Optimization

The current article handles the chance constraints by aggregating them with
the objective function through the expected feasible improvement acquisition
criterion that was first proposed in [40]. Other techniques to account for
constraints in Bayesian optimization exist.

In [41], the constraints are metamodeled but they are not aggregated with the
objective function. The acquisition of a new point is a constrained
optimization problem involving only the metamodels of the objective function
and the constraints. This problem, which is not as costly as the original one,
is treated with a nonlinear constrained optimization algorithm.

Augmented Lagrangians have been adapted to Bayesian optimization in the
presence of equality and inequality constraints in [42]. Slack variables are
introduced to transform the inequalities into equalities. An expected
improvement formulated in terms of augmented Lagrangian values is given.

In [43], a constrained Bayesian optimization is carried out inside a filter-
based generalized pattern search algorithm. Candidate points are created by
the sequential application of two acquisition criteria, an expected violation
of the constraints that must be below a threshold followed by an expected
improvement.

###  2.3 Bayesian optimization under uncertainty

It is no wonder that Bayesian optimization has often been adapted to account
for uncertainties. When the uncertainties are given as probability densities
and the performance is quantified by risk measures, Gaussian processes, which
are probabilistic mathematical objects, come as a natural approach. The
adaptation of BO to uncertainties affects either the Gaussian process (the
parameters over which it is defined, the function that it models) and/or the
acquisition criterion.

##### Robust optimization.

Many previous works have focused on the robust optimization of a function in
the absence of constraints. When dealing with noisy objective functions,
nugget effects have been dimensionned in proportion to the noise and added to
the Gaussian process [44]. The expected improvement acquisition criterion has
been modified in various ways to account for this noise [45]. New acquisition
criteria have been proposed to accomodate this noise. Risk adverse criteria
are based on quantiles [46, 47, 48] while a risk neutral criterion is the
correlated knowledge gradient [49]. A risk adverse approach is described in
[50] that encompasses both quantiles and expectiles: specific Gaussian
processes are given for the quantiles and expectiles along with two BO
strategies, a generalization of the upper-confidence bound and a batch
sampling strategy. The compromises between risk adverse and risk neutral
criteria have been explored through a multi-objective approach in [51].

A key idea in the treatment of uncertainties is the joined space, which is the
tensorization of the deterministic space of the optimization variables with
the space of the random parameters. It opens the possibility to define the
Gaussian processes and the acquisition in the joined space. At each iteration
then, a BO can provide not only a new value for the optimization variables,
but also a new sample for the random parameters. Such a procedure is a
simultaneous optimization and optimal sampling, where the optimality of the
sampling is in the sense of the acquisition criterion. An early example is
given in [52] where the mean objective function is minimized and the
acquisition criterion is the one-step-ahead variance of the objective function
at desirable locations. The same mean objective function is considered, still
in a joined space, in [53]. The variables and random parameters are samples of
an adapted importance sampling density which is calculated with the GP and
which accounts sequentially for the contribution to the objective function and
to the uncertainty in the objective function.
In [47], the acquisition criterion is a knowledge gradient applied to a
quantile (a conditional-value-at-risk) of the objective function. Because of
the quantile, this work bears a resemblance to [46, 50] but in the joined
space, therefore also providing random parameters values.

##### Optimization with reliability constraints.

The estimation of reliability constraints, outside of any optimization, with
the help of Gaussian processes has received a lot of attention [54, 55, 56,
57]. Such reliability analyses based on metamodels have then been included
within optimization searches. A recent survey and a global framework for RBDO
can be found in [58]. Several of these contributions involve kriging
metamodels built over the space of design variables [59, 60, 61].

More recently, like in the case of noisy objective functions, optimization
with reliability constraints has been tackled through Gaussian processes built
in a space that allows to select both the variables and the noise parameters.
In [16], [17] and [22] because the noise directly affects some of the
variables, a formulation in the _augmented space_ (not to be mistaken with the
joined space) is developped: the new optimization variables are parameters of
the noise density functions and the Gaussian processes are built in the space
of the original design vriables. The acquisition criteria, which command the
computationally intensive simulations, primarily aim at refining the
estimation of the reliability constraints. In [16], the acquisition criterion
is the probability of a point to be in the margin of uncertainty of the
constraint. In [17], it is the deviation number (a ratio of the average and
standard deviation of the kriging model of the constraint). In both works, the
distribution variables are updated independently of calls to the simulation by
a direct iteration on the optimization problem.

[21] provides a complete Bayesian problem formulation close to ours: the mean
objective function is minimized with reliability constraints; the Gaussian
processes are defined over the joined space; the acquisition criterion stems
from a similar two steps logic, the first step being a maximization of the
feasible improvement to define the optimization variable and the second step
being a maximization of the posterior (integrated) variance to define the
noise parameters. However the constraints are independent.

To the authors knowledge, this article is the first fully Bayesian formulation
of the problem of minimizing the expected value of a function under chance
constraints where the constraints are not independent, which allows to select
which constraint needs to be evaluated.

##  3 Problem formulation

We work in the joined space of the optimization variables,
$\mathbf{x}\in\mathcal{S_{X}}\in\mathbb{R}^{d}$, over which the user has
control (e.g., targeted dimensions, masses, pressures), and the random
parameters $\mathbf{U}\in\mathcal{S_{U}}\in\mathbb{R}^{m}$ (e.g., atmospheric
parameters, manufacturing errors, numerical noise) for which a probability
density function $\rho_{U}$ is known. Although the $\mathbf{U}$’s are not
controlled in the real world, they can be chosen during the numerical
simulations through their samples $\mathbf{u}$. We assume indeed that the
objective and the constraint functions depend solely on $\mathbf{x}$ and
$\mathbf{u}$ and can be calculated with simulation codes that are
deterministic. The postulate of $\rho_{U}$ is a strong hypothesis that
considerably simplifies the calculations. Less restrictive approaches to the
uncertainty description would still require expert knowledge in the form of a
parameterized family of probability density functions that must be selected.

We consider the following optimization problem where $f(\cdot)$ is the
objective function and where there are $l$ constraints
$\\{g_{1}(\cdot),\dots,g_{l}(\cdot)\\}$:

$\begin{split}&\text{find
}\quad\mathbf{x}^{\star}=\arg\min_{\mathbf{x}\in\mathcal{S_{X}}\subset\mathbb{R}^{d}}\mathbb{E}_{\mathbf{U}}[f(x,U)]\\\
&\text{such that }\quad\mathbb{P}_{\mathbf{U}}(g_{p}(\mathbf{x}^{\star},U)\leq
0,p=1,\dots,l)\geq 1-\alpha\end{split}$ (1)

The objective is to minimize the mean of $f(\cdot)$ with respect to
$\mathbf{U}$ (also called Bayes risk) while simultaneously satisfying the
constraints with a probability above a given threshold $\alpha$ (called
reliability or chance constraint). Other problem formulations that include
uncertainties exist [62]. Of course, the formulation conditions the
characteristics of the optimal design. In our work, the emphasis is placed on
constraint satisfaction with a guaranteed reliability. The effect of the
uncertain variables on the objective function is taken in expectation because
of its analytical and numerical tractability. The average performance is the
appropriate formulation in the frequent situations where the objective $f$
will be measured many times. Examples include the life span of electric bulbs,
the energy consumption of a machine, the speed of convergence of an algorithm.

The objective function is considered to be independent from the constraints,
while the constraints are taken as mutually dependent. Other formulations with
quantiles conditioned by feasibility have been described in [63], yet they
remain an open challenge for costly problems. Problem (1) can be rewritten by
expressing the probability of feasibility as the expectation of a Bernoulli
variable representing whether or not all of the constraints are satisfied:

$\displaystyle\mathbb{P}_{\mathbf{U}}(g_{p}(\mathbf{x},\mathbf{U})\leq
0,~{}p=1,\dots,l)\geq 1-\alpha$ $\displaystyle\Leftrightarrow$
$\displaystyle\mathbb{E}_{\mathbf{U}}[\mathbbm{1}_{\\{g_{p}(\mathbf{x},\mathbf{U})\leq
0,~{}p=1,\dots,l\\}}]\geq 1-\alpha,$ $\displaystyle\Leftrightarrow$
$\displaystyle
1-\alpha-\mathbb{E}_{\mathbf{U}}[\mathbbm{1}_{\\{g_{p}(\mathbf{x},\mathbf{U})\leq
0,~{}p=1,\dots,l\\}}]\leq 0,$ $\displaystyle\Leftrightarrow$ $\displaystyle
c(\mathbf{x})\leq 0.$

The optimization problem becomes

$\mathbf{x}^{\star}=\arg\min\limits_{\mathbf{x}\in\mathcal{S_{X}}}~{}z(\mathbf{x})~{}\text{s.t.}~{}c(\mathbf{x})\leq
0~{},$ (2)

where $z(.)=\mathbb{E}_{\mathbf{U}}[f(.,\mathbf{U})]$.

##  4 Robust Bayesian optimization with a multi-output model of the
constraints

Problems (1) or (2) have been treated with a robust Bayesian approach in [23]
where it was assumed that the constraints were uncorrelated. Such an
assumption has two main consequences: first of all, the constraints can be
modeled by independent GPs which are trained with separate data sets. Second,
the probability that all of the constraints are satisfied at a given location
in the design space can be computed as the product of the probabilities of
feasibility of each constraint.

However, the independence assumption may be overly simplistic and may
potentially limit the algorithm performance. Real constraints are often
positively or negatively correlated to each other because they are computed
from the same underlying physical phenomenon. An example is the correlation
between the maximum displacement and the Von Mises stress on a structure in
bending. Another example is the negative correlation between manufacturing
cost and the quality (e.g., the dimensional accuracy) of a product. Improving
the model accuracy in the context of the optimisation of the return of a
portofolio under a constraint on its variance is also a relevant example. The
expected return is maximized under a constraint on the return variance, but
the variance of the portfolio depends on the correlation between its financial
assets.

For the abovementioned reasons, in this article, we consider correlated
constraints. Not only will our model be more realistic, but it will have the
ability to share information between the constraints, thus improving the
accuracy of the model.

###  4.1 Gaussian regression for a multi-output model of the constraints

Let’s assume that the objective function $f$ and the constraints
$\mathbf{g}(\cdot)=\\{g_{1}(\cdot),\dots,g_{l}(\cdot)\\}$ are realizations of
a Gaussian process $F$ and a multi-output GP $\mathbf{G}$, respectively.
$\mathbf{G}$ is furthermore assumed to be independent of $F$, and both $F$ and
$\mathbf{G}$ are built in the joint space
$\mathcal{S_{X}}\times\mathcal{S_{U}}$. As GPs, $F$ and $\mathbf{G}$ are
characterized by a prior mean and a covariance function:

$\displaystyle F(\mathbf{x},\mathbf{u})$ $\displaystyle\sim$
$\displaystyle\mathcal{GP}(m_{F}(\mathbf{x},\mathbf{u}),k_{F}(\mathbf{x},\mathbf{u},\mathbf{x}^{\prime},\mathbf{u}^{\prime}))$
(3) $\displaystyle\mathbf{G}(\mathbf{x},\mathbf{u})=\left(\begin{array}[]{c
}G_{1}(\mathbf{x},\mathbf{u})\\\ \vdots\\\
G_{l}(\mathbf{x},\mathbf{u})\end{array}\right)$ $\displaystyle\sim$
$\displaystyle\mathcal{GP}(\mathbf{m}_{\mathbf{G}}(\mathbf{x},\mathbf{u}),\mathbf{k}_{\mathbf{G}}(\mathbf{x},\mathbf{u},\mathbf{x}^{\prime},\mathbf{u}^{\prime}))$
(7)

where $m_{F}$ and $\mathbf{m_{G}}$ are the mean functions and $k_{F}$ and
$\mathbf{k_{G}}$ are the covariance functions. While $k_{F}$ is scalar valued,
$\mathbf{k_{G}}$ is a matrix indexed in space,

$\mathbf{k_{G}}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime}):(\mathcal{S_{X}}\times\mathcal{S_{U}})\times(\mathcal{S_{X}}\times\mathcal{S_{U}})\rightarrow
l\times l$ (8)

The GPs in the joint space are indexed by the design variables $\mathbf{x}$
and the sampled values of the random variables $\mathbf{U}$, written
$\mathbf{u}$.

The notations $F^{(t)}$ and $\mathbf{G}^{(t)}$ denote Gaussian processes
conditioned by a set of $t$ observations:
$\mathbf{f}^{(t)}=(f(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,f(\mathbf{x}^{t},\mathbf{u}^{t}))$
and
$\mathbf{g}^{(t)}=(\mathbf{g}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,\mathbf{g}(\mathbf{x}^{t},\mathbf{u}^{t}))$
evaluated at
$\mathcal{D}^{(t)}=\\{(\mathbf{x}^{i},\mathbf{u}^{i})~{},~{}i=1,..,t\\}$.

The conditional GP mean and covariance have the following expressions:

$\displaystyle m_{F}^{(t)}(\mathbf{x},\mathbf{u}))$ $\displaystyle=$
$\displaystyle
m_{F}+\mathbf{k}_{F}(\mathbf{x},\mathbf{u})\mathbf{K_{F}}^{-1}(\mathbf{f}^{(t)}-\mathbf{m}_{F})$
(9) $\displaystyle
k_{F}^{(t)}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime})$
$\displaystyle=$ $\displaystyle
k_{F}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime})-\mathbf{k}_{F}(\mathbf{x},\mathbf{u})\mathbf{K_{F}}^{-1}\mathbf{k}_{F}^{\top}(\mathbf{x}^{\prime},\mathbf{u}^{\prime})$
(10)

where
$\mathbf{k}_{F}(\mathbf{x},\mathbf{u})=\\{k_{F}(\mathbf{x},\mathbf{u},;\mathbf{x}^{1},\mathbf{u}^{1}),\dots,k_{F}(\mathbf{x},\mathbf{u},;\mathbf{x}^{t},\mathbf{u}^{t})\\}$,
$(\mathbf{K}_{F})_{i,j}=k_{F}(\mathbf{x}^{i},\mathbf{u}^{i};\mathbf{x}^{j},\mathbf{u}^{j})$
and $\mathbf{m}_{F}$ is a vector containing $t$ times the value $m_{F}$.

The same equations apply to the posterior mean and covariance of the GP
$\mathbf{G}$:

$\mathbf{m}_{\mathbf{G}}^{(t)}(\mathbf{x},\mathbf{u})=\mathbf{m}_{\mathbf{G}}(\mathbf{x},\mathbf{u})+\mathbf{K_{G}}(\mathbf{x},\mathbf{u};\mathcal{D}^{(t)})(\mathbf{K_{G}}(\mathcal{D}^{(t)},\mathcal{D}^{(t)}))^{-1}(\mathbf{g}^{(t)}-\mathbf{M}_{\mathbf{G}})$
(11)
$\mathbf{K_{G}}^{(t)}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime})=\mathbf{K_{G}}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime})\\\
-\mathbf{K_{G}}(\mathbf{x},\mathbf{u};\mathcal{D}^{(t)})(\mathbf{K_{G}}(\mathcal{D}^{(t)},\mathcal{D}^{(t)}))^{-1}\mathbf{K_{G}}(\mathbf{x}^{\prime},\mathbf{u}^{\prime};\mathcal{D}^{(t)})^{\top}$
(12)

where $\mathbf{K_{G}}(\mathbf{x},\mathbf{u};\mathcal{D}^{(t)})$ is the
$l\times tl$ matrix containing the covariance values over the $l$ outputs
between the sample $\\{\mathbf{x},\mathbf{u}\\}$ and the training data set
$\mathcal{D}^{(t)}$. Note that the dimensions of the vectors and matrices in
the equations above are augmented with respect to the standard GP framework.
For a prediction at a single point,
$\mathbf{K_{G}}^{(t)}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime})\in\mathbb{R}^{l\times
l}$,
$\mathbf{K_{G}}(\mathcal{D}^{(t)},\mathcal{D}^{(t)})\in\mathbb{R}^{tl\times
tl}$,
$\mathbf{K_{G}}(\mathbf{x},\mathbf{u};\mathcal{D}^{(t)})\in\mathbb{R}^{l\times
tl}$.

Several approaches exist to model the structure of correlation between the $l$
components of $\mathbf{G}$. Without loss of generality, the “output as input”
(oai) approach is considered here. The multi-output GP is seen as a scalar GP
defined on an increased input space, with an additional nominal input variable
whose levels characterize one of the $l$ outputs. The additional nominal input
is considered as an integer variable ranging from $1$ to $l$ and referred to
with the letter $p$. Let $G_{p}$ and $G_{p^{\prime}}$ be two components of
$\mathbf{G}$,
$Cov(G_{p}(\mathbf{x},\mathbf{u}),G_{p^{\prime}}(\mathbf{x}^{\prime},\mathbf{u}^{\prime}))$
is assumed to be equal to
$k^{oai}_{\mathbf{G}}(\\{\mathbf{x},\mathbf{u},p\\},\\{\mathbf{x}^{\prime},\mathbf{u}^{\prime},p^{\prime}\\})$
where

$k^{oai}_{\mathbf{G}}(\\{\mathbf{x},\mathbf{u},p\\},\\{\mathbf{x}^{\prime},\mathbf{u}^{\prime},p^{\prime}\\}):(\mathcal{S_{X}}\times\mathcal{S_{U}}\times\mathbb{N})\times(\mathcal{S_{X}}\times\mathcal{S_{U}}\times\mathbb{N})\rightarrow\mathbb{R}$

Note that the “output as input” covariance kernel
$k^{oai}_{\mathbf{G}}(\cdot,\cdot)$ is a scalar function. Different ways of
characterizing the aforementioned kernel exist. For instance,
$k^{oai}_{\mathbf{G}}(\cdot,\cdot)$ can be expressed as the product between
the input dependent term and the output dependent term:

$k^{oai}_{\mathbf{G}}(\\{\mathbf{x},\mathbf{u},p\\},\\{\mathbf{x}^{\prime},\mathbf{u}^{\prime},p^{\prime}\\})=k_{x}(\mathbf{x},\mathbf{u};\mathbf{x}^{\prime},\mathbf{u}^{\prime})*k_{p}(p,p^{\prime})$
(13)

Although less common than the standard continuous kernels, a few kernels
allowing to characterize the covariance function between the values of a
discrete unordered variable such as $p$ exist in the literature. Different
examples can be found in [64, 65, 66]. Given that the main focus of this work
is not to optimize the modeling performance of a multi-output GP model, but
rather to study its effect within the framework of robust optimization, a
single discrete kernel parameterization known as the hypersphere decomposition
[67] is considered for the remainder of the article. The underlying idea is
detailed in Appendix B.

Returning to the optimization Problem (2), models for the functions $z$ and
$c$ in the design space $\mathcal{S_{X}}$ are deduced from the models of the
objective function and the constraints in the joint space. Since the
expectation is a linear operator, the averaged process over the uncertain
domain remains Gaussian. A GP modeling the average objective function has
therefore been introduced:

$Z^{(t)}(\mathbf{x})=\mathbb{E}_{\mathbf{U}}[F^{(t)}{(\mathbf{x},\mathbf{U})}]$

The mean and covariance function of $Z$ are:

$\displaystyle m^{(t)}_{Z}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{m}}m^{(t)}_{F}(\mathbf{x},\mathbf{u})\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u},$
(14) $\displaystyle k^{(t)}_{Z}(\mathbf{x},\mathbf{x^{\prime}})$
$\displaystyle=$
$\displaystyle\iint\limits_{\mathbb{R}^{m}}k^{(t)}_{F}(\mathbf{x},\mathbf{u},\mathbf{x^{\prime}},\mathbf{u^{\prime}})\rho_{\mathbf{U}}(\mathbf{u^{\prime}})\rho_{\mathbf{U}}(\mathbf{u^{\prime}})d\mathbf{u}d\mathbf{u^{\prime}}$
(15)

The integrals in (15) can be evaluated analytically under some assumptions
(see [52]). When this is not feasible a quadrature scheme can be used to
approximate these integrals.

The probability of feasibility with a given safety margin $\alpha$ is
estimated with the GPs $\mathbf{G}$ through the process

$C^{(t)}(\mathbf{x})=1-\alpha-\mathbb{E}_{\mathbf{U}}[\mathbbm{1}_{\mathbf{G}(\mathbf{x},\mathbf{U})\leq\mathbf{0}}]$

where $\mathbf{0}$ is a $l\times 1$ vector of zeros. It is important to note
that, by construction, $C$ is not Gaussian. This characteristic makes the
analytical tractability of the proposed method more complicated, as it will be
discussed later.

###  4.2 Robust Bayesian optimization

Let us now turn to the adaptation of Bayesian optimization to Problem (2). We
generalize the approach described in [23] to multi-output models of the
constraints. The algorithm performs both optimization and sampling in that the
design variables and the uncertain variables at which the exact problem
functions are to be evaluated are determined in the same iteration. The
working principle of this algorithm is to maximize the feasible improvement,
which stands as the acquisition criterion, but in a one-step-ahead uncertainty
reduction strategy. The method is detailed in the following paragraphs.
The most promising location of the design space, denoted $x_{targ}$ for
targeted $x$, is the set of design variables which maximizes the Expected
Feasible Improvement (EFI) :

$\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}\mathbb{E}[FI^{(t)}(\mathbf{x})]$
(16)

where
$FI^{(t)}(\mathbf{x})=\left(z_{min}^{feas}-Z^{(t)}(\mathbf{x})\right)^{+}\mathbbm{1}_{\\{C^{(t)}(\mathbf{x})\leq
0\\}}$.
By exploiting the independence between the objective function and the
constraints, Equation (16) can be written as:

$\displaystyle\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}EI^{(t)}(\mathbf{x})\mathbb{P}(C^{(t)}(\mathbf{x})\leq
0).$ (17)

The first term of (17), the Expected Improvement, is a common acquisition
criterion and it is written

$EI^{(t)}(\mathbf{x})=(z_{\min}^{\text{feas}}-m_{Z}^{(t)}(\mathbf{x}))\Phi\bigg{(}\frac{z_{\min}^{\text{feas}}-m_{Z}^{(t)}(\mathbf{x})}{\sigma_{Z}^{(t)}(\mathbf{x})}\bigg{)}+\sigma_{Z}^{(t)}(\mathbf{x})\phi\bigg{(}\frac{z_{\min}^{\text{feas}}-m_{Z}^{(t)}(\mathbf{x})}{\sigma_{Z}^{(t)}(\mathbf{x})}\bigg{)},$
(18)

where $\sigma_{Z}^{(t)}(\mathbf{x})$ is the standard deviation of the
conditioned GP $Z$ at $\mathbf{x}$, while $\Phi$ and $\phi$ are the normal
cumulative distribution and density functions, respectively.
$z_{\min}^{\text{feas}}$ designates the best feasible target value for the
mean function.

An expression of Equation (17) is given in [23] but it must be adapted to the
multi-output context. The second term of (17),
$\mathbb{P}(C^{(t)}(\mathbf{x})\leq 0)$, was estimated numerically with the
help of a Monte Carlo routine. This computation is extended to the multi-
output model: the simulation of independant GPs is replaced by the simulation
of multi-output GP trajectories. In Equation (18), the expression of
$z_{\min}^{\text{feas}}$ involves the cumulative density function (cdf) of a
multivariate Gaussian distriubution instead of the product of $l$ Gaussian
cdf’s. The details of these calculations are provided in Appendix C.

Once $x_{targ}$ is determined, the design of experiments $\mathcal{D}^{(t)}$
is enriched with the next iterate $\\{(\mathbf{x}^{t+1},\mathbf{u}^{t+1})\\}$.
This iterate is such that a proxy $S$ of the one-step-ahead variance of the
acquisition criterion (the feasible improvement) is minimized at $x_{targ}$.
$S$ will serve as a Sampling criterion for the random variables, hence the
associated symbol. In [52], it is shown that the optimal values of
$\mathbf{x}^{t+1}$ tend to be very close to $x_{targ}$. For this reason, in
order to reduce the computational overhead, it is assumed that
$\mathbf{x}^{t+1}=x_{targ}$. Then, only $\mathbf{u}^{t+1}$ needs to be
computed, which is done as follows:

$\mathbf{u}^{t+1}=\arg\min_{\mathbf{\tilde{u}}}S(\mathbf{x}_{targ},\mathbf{\tilde{u}})~{}.$
(19)

A global view of the entanglement of the above steps within a constrained,
robust, optimization method can be gained by looking at Algorithm 1 for which
other options will be introduced later (see Section 6.1).

The proxy to the one-step-ahead variance of the feasible improvement is
written

$S(\tilde{\mathbf{x}},\mathbf{\tilde{u}})=~{}Var(\big{(}z_{\min}^{\mbox{feas}}-Z^{(t+1)}(\mathbf{x}_{targ})\big{)}^{+})\\\
\int_{\mathbb{R}^{m}}Var\big{(}\mathbbm{1}_{\\{\mathbf{G}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u})\leq
0\\}}\big{)}\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}$ (20)

where $Z^{(t+1)}$ and $\mathbf{G}^{(t+1)}$ are conditioned on
$\mathcal{D}^{(t)}\cup\\{(\tilde{\mathbf{x}},\mathbf{\tilde{u}})\\}$.

The first term of $S(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$, i.e.
$Var\left(\big{(}z_{\min}^{\mbox{feas}}-Z^{(t+1)}(\mathbf{x}_{targ})\big{)}^{+}\right)$,
corresponds to the improvement variance. Its expression has some resemblance
to the one of the expected improvement and is given in [23]. It is expressed
with the probability and cumulated density functions of the standard Gaussian
distribution (see Appendix D).

In order to compute the second term of
$S(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$ with a multi-output constraint
model, one needs to determine $\mathbf{G}^{(t+1)}$ by conditioning
$\mathbf{G}^{(t)}$ on the new sample
$\\{\tilde{\mathbf{x}},\tilde{\mathbf{u}}\\}$ while taking into account the
correlation between the components of $\mathbf{G}^{(t)}$. By drawing on the
Kriging Believer assumption111 The Kriging Believer assumption states that the
value of a function at an unobserved point is equal to the kriging prediction
at that point. For our constraints, it means
$g_{p}(\tilde{\mathbf{x}},\tilde{\mathbf{u}})={{\mathbf{m}}_{\mathbf{G}}^{(t)}}_{p}(\tilde{\mathbf{x}},\tilde{\mathbf{u}})~{},~{}p=1,\ldots,l$.,
the updated prediction is

$\displaystyle\mathbf{m_{G}}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u})$
$\displaystyle=$
$\displaystyle\mathbf{m_{G}}^{(t)}(\mathbf{x}_{targ},\mathbf{u})$
$\displaystyle\mathbf{K_{G}}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u};\mathbf{x}_{targ},\mathbf{u})$
$\displaystyle=$
$\displaystyle\mathbf{K_{G}}^{(t)}(\mathbf{x}_{targ},\mathbf{u};\mathbf{x}_{targ},\mathbf{u})$
$\displaystyle-\mathbf{K_{G}}^{(t)}(\mathbf{x}_{targ},\mathbf{u};\tilde{\mathbf{x}},\mathbf{\tilde{u}})\mathbf{K}_{new}^{-1}\mathbf{K_{G}}^{(t)}(\mathbf{x}_{targ},\mathbf{u};\tilde{\mathbf{x}},\mathbf{\tilde{u}})^{\top}$

with
$\mathbf{K}_{new}=\mathbf{K_{G}}^{(t)}(\tilde{\mathbf{x}},\mathbf{\tilde{u}};\tilde{\mathbf{x}},\mathbf{\tilde{u}})$
the $l\times l$ matrix containing the covariance of the vector
$\mathbf{G}(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$, i.e., the covariance
values between the different outputs of the GP evaluated at
$(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$.

Now,
$\mathbbm{1}_{\mathbf{G}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}}$
is a Bernoulli variable with parameter $p(\mathbf{u})$ equal to

$p(\mathbf{u})=\Phi\left(\mathbf{0}-\mathbf{m_{G}}^{(t)}(\mathbf{x}_{targ},\mathbf{u}),\mathbf{K}^{(t+1)}_{\mathbf{G}}(\mathbf{x}_{targ},\mathbf{u};\mathbf{x}_{targ},\mathbf{u})\right)~{},$
(22)

where $\Phi(\cdot)$ is the cumulative distribution function of a multivariate
Gaussian distribution. As a consequence,

$\displaystyle\int_{\mathbb{R}^{m}}Var(\mathbbm{1}_{\mathbf{G}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}})\rho(\mathbf{u})d\mathbf{u}=\int_{\mathbb{R}^{m}}p(\mathbf{u})(1-p(\mathbf{u}))\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}$
(23)

As it will be seen in Section 4.3, modeling the correlation between the
constraints increases the accuracy of the surrogate. However, the difference
in performance between the two variants (with or without constraint
correlation) will largely depend on the actual correlation value between the
modeled functions. If the constraints are independent for example, the multi-
output model will be more complex than the standard approach while providing
no additional benefit. More generally, the potential gain in performance of
the coupled model depends on the ability of the correlation structure to
describe the constraints relationship. In any case, accounting for constraints
couplings results in a substantial increase of the covariance matrix size,
which may become a practical bottleneck with large data sets. When this
happens, a sparse GP covariance structure might be indicated.

###  4.3 Illustration of feasibility estimation with a multi-output GP

In order to illustrate the behaviour of the multi-output model of the
constraints, an analytical example where the constraints are designed to be
linked is considered:

$\begin{split}&\min_{\mathbf{x}}\mathbb{E}_{U}[f(\mathbf{x},\mathbf{U})]~{}\text{
such that }~{}\mathbb{P}(g_{p}(\mathbf{x},\mathbf{U})\leq 0,p=1,2)\geq 0.95\\\
&\text{where }\\\
&f(x_{1},x_{2},u_{1},u_{2})=5(x_{1}^{2}+x_{2}^{2})-(u_{1}^{2}+u_{2}^{2})+x_{1}(u_{2}-u_{1}+5)+x_{2}(u_{1}-u_{2}+3)\\\
&g_{1}(x_{1},x_{2},u_{1},u_{2})=-x_{1}^{2}+5x_{2}-u_{1}+u_{2}^{2}-1\\\
&g_{2}(x_{1},x_{2},u_{1},u_{2})=g_{1}(x_{1},x_{2},u_{1},u_{2})(x_{1}+5)/5-u_{1}-1\end{split}$
(24)

with $\mathbf{x}\in[-5,5]^{2}$ and $\mathbf{U}\sim\mathcal{U}([-5,5]^{2})$. A
correlation is induced between the two constraints by making the expression of
the second constraint a transformation of the expression of the first
constraint.

Depending on whether an independent or a dependent model of the constraints is
considered, the following expression of the Probability of Feasibility (PoF)
stands:

$\displaystyle\mbox{ Independent constraints: }\ \ \text{PoF}(\mathbf{x})$
$\displaystyle=$ $\displaystyle
E_{\mathbf{U}}[\mathbbm{1}_{\cap_{p=1}^{2}\\{G_{p}(\mathbf{x},\mathbf{U})\leq
0\\}}]$ $\displaystyle\mbox{ Dependent constraints: }\ \
\text{PoF}(\mathbf{x})$ $\displaystyle=$ $\displaystyle
E_{\mathbf{U}}[\mathbbm{1}_{\mathbf{G}(\mathbf{x},\mathbf{U})\leq\mathbf{0}}]$

The modeling performance is calculated as the absolute difference between the
PoF computed with the GP models, and the one computed on the true function at
$N=1000$ $\mathbf{u}^{j}$’s sampled with respect to $\mathcal{U}([-5,5]^{2})$,

$\text{PoF
error}(\mathbf{x})~{}=~{}|\text{PoF}(\mathbf{x})-\frac{1}{N}\sum_{j=1}^{N}1_{\cap_{p=1}^{2}\\{g_{p}(\mathbf{x},\mathbf{u}^{j})\leq
0\\}}|~{}.$

The results are plotted in Figure 1 where, for the sake of clarity, the input
variables are normalized between 0 and 1 in the plots. A data set of 30
samples per constraint function is used to train the scalar and the multi-
output GPs.

Figure 1: Contour plots of the error in the probability of feasibility on the
Problem (24). Left: independent modeling of each constraint. Right: multi-
output GP model of dependent constraints. The training points in the
$\mathcal{S_{X}}$ space are indicated with the orange dots.

The Figure shows that the areas of $\mathcal{S_{X}}$ characterized by the
largest accuracy (lowest error) are the same in both models. Overall however,
the multi-output GP provides a lower PoF estimation error than the independent
model. The gain in accuracy is particularly clear in the parabolic region
where the constraints are near 0, which is the considered feasibility
threshold. The advantage of the multi-output constraint model is further
highlighted in Figure 2 where the mean PoF error computed on 400 randomly
sampled values in $S_{X}$ over 10 repetitions of the DoE and GPs training is
shown for both approaches.

Figure 2: Average error in probability of feasibility on Problem (24). The
left boxplot is obtained with a an independent GP for each constraint, whereas
the right boxplot corresponds to a correlated multi-output GP.

The difference in modeling performance can be explained by the fact that the
multi-output GP benefits from a larger data set since it processes information
from the training data of both constraints. In the context of optimization, it
is reasonable to assume that a more accurate feasibility prediction will allow
a better search for robust optimal solutions, which highlights the potential
benefit of a multi-output model of the constraints.

##  5 Extensions to robust Bayesian optimization enabled by constraints
coupling

Two extensions to robust Bayesian optimization are now proposed. The first
extension is the possibility of refining the objective and constraint
functions with different sets of random variables (but identical design
variables point). The second extension allows to select a single constraint to
evaluate at $x_{targ}$ in order to save constraints evaluations. The chosen
constraint provides the largest amount of information regarding the
feasibility at the candidate optimum location $x_{targ}$. These extensions are
relevant when each output, at least the objectif function and the set of
constraints, can be obtained through separate simulators in which case a goal
oriented sampling strategy that sparingly selects which simulator should be
called will save computations.

###  5.1 Independent selection of uncertain variables

The robust Bayesian optimization algorithm of Section 4 associates to each
target design variable, $\mathbf{x}_{target}$, a single vector of random
variables $\mathbf{u}^{t+1}$ that will be used to refine the GPs of both the
objective and the constraints functions. In most cases, however, the most
informative $\mathbf{u}$’s from the objective function improvement and from
the constraint satisfaction perspectives are different. Selecting a single
$\mathbf{u}\in\mathbb{R}^{m}$ for a better accuracy in, both, the improvement
and the failure probability, is necessarily a compromise. For this reason, we
introduce the possibility of selecting different vectors, $\mathbf{u}_{f}$ and
$\mathbf{u}_{g}$, for the objective and the constraint functions. Again, the
underlying assumption in this duplication of the random samples is that the
computationally intensive codes used to evaluate the objective and the
constraints can be called independently.

Remember that the acquisition function of Equation (20),

$S(\mathbf{x}_{targ},\mathbf{u})=Var((z_{min}^{feas}-Z^{(t+1)}(\mathbf{x}_{targ}))\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}})\rho(\mathbf{u})d\mathbf{u}~{},$

is composed of two terms, the first related to the objective function and the
second to the constraints. It can be decomposed into two independent terms,
$S_{f}$ and $S_{g}$ involving different samples of the random variables :

$\displaystyle S_{f}(\mathbf{x}_{targ},\mathbf{u}_{f})$ $\displaystyle=$
$\displaystyle Var((z_{min}^{feas}-Z^{(t+1)}(\mathbf{x}_{targ})^{+})$
$\displaystyle=$ $\displaystyle
Var((z_{min}^{feas}-Z(\mathbf{x}_{targ}))^{+}|F(\mathcal{D}^{(t)}_{f})=f^{(t)},F(\mathbf{x}_{targ},\mathbf{u}_{f}))$
$\displaystyle S_{g}(\mathbf{x}_{targ},\mathbf{u}_{g})$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}})\rho(\mathbf{u})d\mathbf{u}$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}|\mathbf{G}(\mathcal{D}^{(t)}_{g})=\mathbf{g}^{(t)},\mathbf{G}(\mathbf{x}_{targ},\mathbf{u}_{g})})\rho(\mathbf{u})d\mathbf{u}$

where $\mathcal{D}^{(t)}_{f}$ and $\mathcal{D}^{(t)}_{g}$ are design of
experiments after $t$ evaluations that contain different samples of $u$. More
precisely if $t_{init}$ denotes the size of the initial DOE, we have
$\mathcal{D}^{(t_{init})}_{f}=\mathcal{D}^{(t_{init})}_{g}=\mathcal{D}^{(t_{init})}$.
For $t>t_{init}$ , the DOE for $f$ (for $g$, respectively) evolves as
$\mathcal{D}^{(t)}_{f}=\mathcal{D}^{(t-1)}_{f}\cup\\{(\mathbf{x}^{t},\mathbf{u}^{t}_{f})\\}$,
($\mathcal{D}^{(t)}_{g}=\mathcal{D}^{(t-1)}_{g}\cup\\{(\mathbf{x}^{t},\mathbf{u}^{t}_{g})\\}$,
respectively).

As a consequence, it is possible to separately optimize $S_{f}$ and $S_{g}$,
thus obtaining two separate vectors of uncertain variables, one for the next
objective function and one for the next constraints evaluation:

$\displaystyle\mathbf{u}_{f}^{t+1}$ $\displaystyle=$
$\displaystyle\arg\min_{\mathbf{u}_{f}}S_{f}(\mathbf{x}_{targ},\mathbf{u}_{f})$
(27) $\displaystyle\mathbf{u}_{g}^{t+1}$ $\displaystyle=$
$\displaystyle\arg\min_{\mathbf{u}_{g}}S_{g}(\mathbf{x}_{targ},\mathbf{u}_{g})$
(28)

This enables the algorithm to more efficiently refine the surrogate models of
the problem functions : the set of variables $\mathbf{u}_{f}$ minimizes the
one-step-ahead variance of the improvement, and $\mathbf{u}_{g}$ minimizes the
one-step-ahead variance of the probability of feasibility. The generic
Algorithm 1 may help the reader in positioning the constraint selection
operation within the global optimization procedure (see Section 6.1).

To highlight the effect of this modification, we return to the analytical
Problem (24). A data set of size $t=20$ is created, GPs of the objective
function and the two constraints are learned, and the most promising location
in the design space $\mathbf{x}_{target}$ determined by maximization of the
Expected Feasible Improvement (Equation (17)). The GPs and
$\mathbf{x}_{target}$ can then be used in the formulas for either $S$, the
proxy to the one-step ahead feasible improvement variance, or its terms
$S_{f}$ and $S_{g}$. At that point, these sampling criteria only depend on the
random parameters $u_{1}$ and $u_{2}$. A contour plot of $S$ as a function of
$u_{1}$ and $u_{2}$ is drawn in Figure 3. It is meant to be compared to the
contours of $S_{f}$ and $S_{g}$, the one-step-ahead variance of the expected
improvement and of the feasibility probability, given in Figure 4.

Figure 3: Contour plot of $S$, the proxy of the one-step-ahead feasible
improvement variance computed at a given $\mathbf{x}_{target}$.

Figure 4: Contour plots of the one-step-ahead variance of the improvement
(left) and the feasibility probability (right), computed at a given
$\mathbf{x}_{target}$.

Samples of the random parameters are generated by minimization of the one-
step-ahead variances, $S$, $S_{f}$ or $S_{g}$. When comparing the feasible
improvement variance proxy, $S$, with its two factors, $S_{f}$ and $S_{g}$, it
is observed that the improvement variance ($S_{f}$) is the dominating factor
in $S$ and drives the selection of the next samples $\mathbf{u}^{t+1}$ in the
generic algorithm of Equation (19). However, such a $\mathbf{u}^{t+1}$ will
poorly reduce the variance of the probability of feasibility ($S_{g}$ map). In
the above example, $\mathbf{u}^{t+1}=\\{1,0\\}$ when the best sample for the
constraints is $\mathbf{u}^{t+1}_{g}=\\{0.1,0.05\\}$.

Therefore a unique $\mathbf{u}$ may not be the best strategy to efficiently
learn where is the feasible domain.

###  5.2 Iterative optimal constraint selection

When the constraints relate to different physics or disciplines, it is common
to have separate simulation codes for each constraints function (or group of
constraints).

A direct extension of the previous section would then be to consider different
instances of the random variables for each of the constraint. The acquisition
function could be detailed for each of the constraint random sample
$\mathbf{u}_{g_{i}}$

$\displaystyle
S_{g}^{all}(\mathbf{x}_{targ},\mathbf{u}_{g_{1}},...,\mathbf{u}_{g_{l}})=$
(29)
$\displaystyle\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}}|\mathbf{G}(\mathcal{D}^{(t)}_{g})=g^{(t)},G_{1}(\mathbf{x}_{targ},\mathbf{u}_{g_{1}}),\ldots,G_{l}(\mathbf{x}_{targ},\mathbf{u}_{g_{l}}))\rho(\mathbf{u})d\mathbf{u}$

where
$\mathcal{D}^{(t_{init})}_{g}=\\{(\mathbf{x}^{i},\mathbf{u}^{i},1),\ldots,(\mathbf{x}^{i},\mathbf{u}^{i},l),\forall
1\leq i\leq t_{init}\\}$ and for $t>t_{init}$,
$\mathcal{D}^{(t)}_{g}=\mathcal{D}^{(t-1)}_{g}\cup\\{(\mathbf{x}^{t},\mathbf{u}^{t}_{g_{1}},1),\ldots,(\mathbf{x}^{t},\mathbf{u}^{t}_{g_{l}},l)\\}$.
Remark that the covariance matrix to be inverted in the update Equation (4.2),
$\mathbf{K}_{new}$, is slightly more complex as it is computed with respect to
data samples characterized by different values of both $\mathbf{u}$ and $p$.

When the codes of each constraint are autonomous, a more appealing extension
for numerical efficiency is to select only one constraint to update at each
iteration. Indeed, the refinement of some constraints may not provide any
additional information when the surrogate models are already accurate enough
to determine feasibility at a given location. Furthermore, the step-by-step
selection of the most relevant constraint to evaluate prevents an unnecessary
increase of the covariance matrix size. This is important in Bayesian
optimization because the dominating numerical cost comes from the repeated
$\mathcal{O}(t^{3})$ inversion of the covariance matrix during the likelihood
maximizations. Adding constraint evaluations one at a time is particularly
tempting when the constraints are modeled with a multi-output GP since all the
GPs improve with that one evaluation.

In order to select the constraint, the one-step-ahead variance of the
probability of feasibitily, $S_{g}$ (see Equation (5.1)), is modified by
including the index $p$ of the constraint function considered, $g_{p}$, and
the associated sampling location $\mathbf{u}_{g}$ :

$S_{g}(\mathbf{x}_{targ},\mathbf{u}_{g},p)=\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}|\mathbf{G}(\mathcal{D}^{(t)}_{g})=g^{(t)},G_{p}(\mathbf{x}_{targ},\mathbf{u}_{g})})\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}$
(30)

where for $t>t_{init}$,
$\mathcal{D}^{(t)}_{g}=\mathcal{D}^{(t-1)}_{g}\cup\\{(\mathbf{x}^{t},\mathbf{u}^{t}_{g},p^{t})\\}$.
In the above Equation, the one-step-ahead update of the GP includes a single
additional data sample for the constraint designated by $p$. In this case, the
covariance matrix of $\mathbf{G}(\mathcal{D}^{(t)})$ has a size smaller than
$tl\times tl$ and its update coming from Equation (4.2) is

$\mathbf{K_{G}}^{(t+1)}(\mathbf{x}_{targ},\mathbf{u};\mathbf{x}_{targ},\mathbf{u})=\mathbf{K_{G}}^{(t)}(\mathbf{x}_{targ},\mathbf{u};\mathbf{x}_{targ},\mathbf{u})-\mathbf{w}k_{new}^{-1}\mathbf{w}^{\top}$

with $k_{new}$ a scalar and $\mathbf{w}$ a $l\times 1$ vector,

$\displaystyle k_{new}$ $\displaystyle=$ $\displaystyle
Cov(G_{p}^{(t)}(\mathbf{x}_{targ},\mathbf{u}_{g}),G_{p}^{(t)}(\mathbf{x}_{targ},\mathbf{u}_{g}))$
$\displaystyle\mathbf{w}$ $\displaystyle=$
$\displaystyle\left(Cov(G_{1}^{(t)}(\mathbf{x}_{targ},\mathbf{u}),G_{p}^{(t)}(\mathbf{x}_{targ},\mathbf{u}_{g})),\ldots,Cov(G_{l}^{(t)}(\mathbf{x}_{targ},\mathbf{u}),G_{p}^{(t)}(\mathbf{x}_{targ},\mathbf{u}_{g}))\right)^{\top}$

For selecting the constraint at update time, the one-step-ahead feasibility
variance is the acquisition function which is simultaneously optimized with
respect to $\mathbf{u}_{g}$ and $p$. The solution is a pair characterizing
both the constraint and the parameter in the $\mathcal{U}$ space which
provides the largest amount of information regarding the feasibility of
$\mathbf{x}_{targ}$,

$p^{t+1},\mathbf{u}_{g}^{t+1}=\arg\min_{p,\mathbf{u}_{g}}S_{g}(\mathbf{x}_{targ},\mathbf{u}_{g},p)\quad.$
(31)

At each iteration, the random sample of the objective function,
$\mathbf{u}_{f}^{t+1}$ comes from the one-step-ahead improvement variance
minimization as expressed in Equation (27). An evaluation of the objective
function is performed together with each evaluation of the selected constraint
at $(\mathbf{x}_{targ},\mathbf{u}_{f}^{t+1})$ and
$(\mathbf{x}_{targ},\mathbf{u}_{g}^{t+1})$, respectively. Hence, for each
objective function evaluation there is a single constraint function
evaluation, rather than $l$ distinct constraint function evaluations. A reason
for this is that once it is known that a constraint is satisfied or violated,
there is no use in fine tuning it. On the contrary, further adjusting the
objective function model in high-performance regions is usually useful.

In order to highlight the potential benefit of selecting the most useful
constraint at every iteration, the analytical test case described in Equations
(24) is considered anew. For illustrative purposes, the $\mathbf{x}_{target}$
value ($\\{0,0.3\\}$) is manually selected so that the constraint $g_{1}$ is
almost certainly satisfied, whereas the probability of feasibility for $g_{2}$
is approximately $0.6$. Each constraint is modeled by an independent GP,
trained by relying on a data set of 30 randomly generated samples. The
quantity of interest is the reduction in one-step-ahead variance of the
feasibility which would be obtained by refining the model of either $g_{1}$
(when $p=1$) or $g_{2}$ ($p=2$),

$\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}|\mathbf{G}(\mathcal{D}^{(t)}_{g})=g^{(t)}})\rho(\mathbf{u})d\mathbf{u}~{}-~{}\\\
\int_{\mathbb{R}^{m}}Var(1_{\mathbf{G}(\mathbf{x}_{targ},\mathbf{u})\leq\mathbf{0}|\mathbf{G}(\mathcal{D}^{(t)}_{g})=g^{(t)},\mathbf{G}(\mathbf{x}_{targ},\mathbf{u}_{g},p)})\rho(\mathbf{u})d\mathbf{u}~{}.$

This quantity is plotted for both values of $p$ (i.e., for either a $g_{1}$ or
a $g_{2}$ enrichment) in Figure 5.

Figure 5: One-step-ahead feasibility variance reduction which would be
obtained by refining the model of either $g_{1}$ (left) or $g_{2}$ (right) at
the value $(\mathbf{x}_{targ},u_{1},u_{2})$, Problem (24)

The left part of the Figure shows that refining the $g_{1}$ constraint which
is surely feasible at $\mathbf{x}_{targ}$ (i.e., the PoF is almost equal to 1)
provides no additional information. On the contrary, the uncertainty reduction
is significant on the right plot when updating $g_{2}$. It is reasonable to
conclude that, for this example, only evaluating $g_{2}$ and not $g_{1}$ would
reduce the optimization cost by avoiding a computationally intensive
simulation ($g_{1}$) and by reducing the increase in data set size, while not
impacting the algorithm performance.

It should be stated that the constraint selection idea that we just introduced
in the context of robust optimization would also stand for deterministic
problems where the constraints are jointly modeled. In fact, it would be
simpler: the one-step-ahead variance at $\mathbf{x}_{targ}$ of Equation (30)
would be directly calculated without the integral on the $\mathbf{u}$’s. It
would be minimized on $p$ only.

##  6 Numerical implementation

###  6.1 Algorithms compared

Four algorithms that correspond to the main combinations of the ideas
presented above, i.e., coupled modeling of the constraints or not and
constraint selection or not, will be compared:

  * •

REF: separate model of each constraint, selection of a single joint value
$(\mathbf{x},\mathbf{u})$ at each iteration where the objective and all the
constraints are evaluated. This is the reference method which was first
proposed in [23], hence the acronym.

  * •

SMCS: Separate Model of the constraints and Constraint Selection. A single
constraint and the associated $\mathbf{u}$ are identified at each iteration to
provide the largest reduction in the probability of feasibility variance.

  * •

MMCU: Multi-output Model of the constraints and Common value of $\mathbf{u}$.
A single value of $\mathbf{u}$ is chosen at each iteration for which, after
composition with $\mathbf{x}^{t+1}$, all the constraints are evaluated.

  * •

MMCS: Multi-output Model and Constraint Selection.

The REF algorithm could also have been called, with the keywords underlying
the other acronyms, SMCU (Separate Model Common $\mathbf{u}$). These four
algorithms are described within a unified template in Algorithm 1 and detailed
in Appendix E. At the end of the four algorithms, the solution returned is the
predicted solution to Problem (1) as seen from the mean of the last GPs,

$\arg\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0$

The last term, $\mathbb{E}[C^{(t)}(\mathbf{x})]\leq 0$, is the expectation of
satisfying the constraints with confidence $1-\alpha$.

Algorithm 1 Template of the robust constrained BO algorithms
(the colors are associated to the variants: coupled or uncoupled constraints,
common $\mathbf{u}$ or constraint selection)

Define initial data sets
$\mathcal{D}_{f}^{(t_{init})},~{}\mathcal{D}_{g_{p}}^{(t_{init})},p=1,\ldots,l$
for the GPs of the objective function and constraints comprised of pairs
$(\mathbf{x},\mathbf{u})\in\mathcal{S_{X}}\times\mathcal{S_{X}}$

Evaluate the objective and the constraint functions through
$f(\mathbf{x},\mathbf{u})~{},~{}(\mathbf{x},\mathbf{u})\in\mathcal{D}_{f}^{(t_{init})}$
and
$g_{p}(\mathbf{x},\mathbf{u})~{},~{}(\mathbf{x},\mathbf{u})\in\mathcal{D}_{g_{p}}^{(t_{init})},~{}p=1,\ldots,l$
$\Rightarrow\mathbf{f}^{(t_{init})},\mathbf{g_{1}}^{(t_{init})},\ldots,\mathbf{g_{l}}^{(t_{init})}$.

Let $t\leftarrow t_{init}$

while $t\leq t_{max}$ do

Train the GPs, $F$ w.r.t. $(\mathcal{D}_{f}^{(t)},\mathbf{f}^{(t)})$

and $\mathbf{G}$ w.r.t.
$\left((\mathcal{D}_{g_{1}}^{(t)},\mathbf{g_{1}}^{(t)}),\ldots,(\mathcal{D}_{g_{l}}^{(t)},\mathbf{g_{l}}^{(t)})\right)$

in a coupled (MMCU, MMCS) or independent (REF, SMCS) way

Compute the targeted variable by maximizing the EFI (Equation (17)),

$\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}\mathbb{E}[FI^{(t)}(\mathbf{x})]$

Set $\mathbf{x}^{t+1}=\mathbf{x}_{targ}$

Compute the random parameters :

  * •

either (REF, MMCU) in a unified fashion for the objective and constraints by
minimizing the proxy to the one-step-ahead EFI variance (Equation (20)),
${\color[rgb]{0,.5,.5}\mathbf{u}^{t+1}}=\arg\min_{\mathbf{u}}S(\mathbf{x}_{targ},\mathbf{u})$

  * •

or (MMCS, SMCS) separately for the objective function and the selected
constraint:

    * –

$\mathbf{u}_{f}^{t+1}$ minimizes the one-step-ahead improvement variance
(Equation (5.1)),
$\mathbf{u}_{f}^{t+1}=\mbox{arg}\min_{\mathbf{u}}S_{f}(\mathbf{x}_{targ},\mathbf{u})$

    * –

$\mathbf{u}_{g}^{t+1}$ and the constraint index minimize the one-step-ahead
variance of the feasibility probability (Equation (30)) ,
${\color[rgb]{0,1,0}\\{\mathbf{u}_{g}^{t+1},p\\}}=\mbox{arg}\min_{\mathbf{u},p^{\prime}}S_{g}(\mathbf{x}_{targ},\mathbf{u},p^{\prime})$

Update $\mathcal{D}_{f}^{(t)},~{}\mathcal{D}_{g_{p}}^{(t)},p=1,\ldots,l$

Calculate $f(,)$ and (one – MMCS, SMCS – or many – REF, MMCU –)

constraints at the new joined point(s), update $\mathbf{f}^{(t)}$ and
$\mathbf{g_{1}}^{(t)},\ldots,\mathbf{g_{l}}^{(t)}$

$t\leftarrow t+1$

end while

return the best feasible point according to the means of the processes,

$\arg\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0$

###  6.2 Remarks about the implementation and code availability

The algorithms presented and compared in this article are implemented in R.
The scalar, independent, GPs are programmed thanks to the DiceKriging [68]
toolbox. The multi-ouput GPs rely on the kergp [69] toolbox which facilitates
the implementation of the kernel seen in Equation (13). Please note that the
source code allowing to reproduce the results presented in this paper on a
simple 2-dimensional test case can be found on the following repository :
https://github.com/JPelamatti/Robust_Bayesian_Optimization.
We now comment on their numerical implementations.

####  6.2.1 Multi-output GP likelihood optimization

When dealing with the multi-output model of the constraints, the conditioning
data set is often made of the various constraints evaluated at the same
location. That is, the data set may contain several points with the same
$\mathbf{x}$ and $\mathbf{u}$ which only differ in $p$, the constraint
indicator. This setting affects the algorithms which are eligible to maximize
the GP likelihood. Indeed, the numerical computation of the likelihood
gradient will return null gradient terms with respect to several GP hyper-
parameters, namely the parameterization of the hypersphere kernel (Appendix B)
for the discrete indices $p$ contains $\theta_{i,j}$ hyper-parameters whose
gradients vanish. For this reason, the training of the multi-output GP models
was performed with the COBYLA [70] optimization algorithm, which does not
require gradients while still being quite efficient. COBYLA searches were
restarted 20 times from random initial points to make the likelihood
maximization more robust.

####  6.2.2 Other internal optimizations

On top of the Gaussian process hyperparameter learning discussed above, the
algorithms of this article include several internal optimizations such as the
maximization of the feasible improvement (Equation (17)), the minimization of
the one-step-ahead variance (Equation (19)) and the computation of a feasible
mean performance (Equation (35)).

Unlike standard Bayesian optimization, these auxiliary optimizations have a
non-negligible computational cost because they require Monte Carlo simulations
to estimate expected values and probabilities of the candidate solutions (this
is for instance illustrated in Appendix D).

In order to afford a fair comparison between the algorithms studied, the
results of Section 7 are obtained by optimizing the acquisition functions
(Equations (17) and (19)) with respect to a predetermined set of candidate
solutions randomly generated through a Latin Hypercube Sampling (LHS)
procedure. The size of the considered candidate solutions set varies as
$\simeq 500\times$dimension. This implementation is not viable in larger
dimensional problems (typically in dimension 5 and above), as the number of
candidate solutions required to sufficiently fill the search space would be
too large. In these cases, proper optimization algorithms (_e.g.,_ gradient-
based, COBYLA, heuristic algorithms) are required to keep the number of
candidate solutions evaluations at a computable level.

Notice also that the evaluations of the acquisition functions at different
locations of the joint space are independent from each other. This property
could be exploited to evaluate several candidate solutions in parallel, thus
allowing toreduce the computation time required by the optimization of these
acquisition functions.

##  7 Experiments and results

This Section reports the experiments done with the four algorithms (REF, SMCS,
MMCU, MMCS) on two analytical test cases as well as on a more realistic
engineering design problem.

###  7.1 Experiments set-ups and performance metrics

In order to take into account the variability coming from the initial data
set, the runs are repeated 10 times with different initial data sets obtained
through an optimized LHS, while the random variables are independently sampled
from their specific distributions. The initial data set of the multi-
constraints model is the union of the data sets for the scalar outputs. It is
therefore $l$ times larger than each training set of the scalar GPs.

The same computational budget is allocated to every method. For comparing
methods with and without constraint selection, the budget is expressed in
terms of the number of constraints evaluations and an iteration is defined as
$l$ evaluations of one constraint, regardless of repeated evaluations.

At a given iteration, the performance of the algorithms is assessed in terms
of the best current objective function value that is feasible with a
probability of $1-\alpha$ where $\alpha=0.05$ or $\alpha=0.1$ depending on the
considered problem. This value is computed with respect to the true problem
functions.

On top of the previously discussed likelihood optimization, the BO algorithms
studied embed auxiliary optimizations: the identification of the current
feasible minimum, the optimization of the expected feasible improvement and
the minimization of the one-step-ahead variance. To limit the impact of the
randomness in these auxiliary calculations on the tests, the auxiliary
optimization are performed by determining the optimum among a fixed set of
points. A different set of points is used for each optimization repetition.

###  7.2 2-dimensional analytical test case

The first test case depends on 1 design variable and 1 random variable, and is
subject to two constraints:

$\begin{split}&\min_{\mathbf{x}}\mathbb{E}_{U}[f(\mathbf{x},\mathbf{U})]\ \
s.t.\ \ \mathbb{P}(g_{p}(\mathbf{x},\mathbf{U})\leq 0,p=1,2)\geq 0.95~{},\\\
&\text{where }\quad f(x,u)=(x-10)^{3}+(u-20)^{3}~{},\\\ &\phantom{\text{where
}\quad}g_{1}(x,u)=-(x-5)^{2}-(u-5)^{2}+500~{},\\\ &\phantom{\text{where
}\quad}g_{2}(x,u)=(x-6)^{2}+(u-5)^{2}-9000~{},\end{split}$ (32)

$\mathbf{x}\in[13,100]$ and $\mathbf{U}\sim\mathcal{U}([0,100])$. This test
case is intended to highlight the potential advantage of multi-output GP
modeling thanks to the large correlation between the constraints.
Additionally, given the expression of the two constraints, the covariance
between $p=1$ and $p=2$ within the multi-output GP is expected to be negative.
The discrete kernel parameterization should therefore be able to return
negative values, which happens with the hypersphere decomposition.

The GP models of the constraints and the objective function are initialized
with a data set of 6 samples per function. A budget of 40 constraint
evaluations is allocated to the optimization iterations. The convergence rates
of the 4 algorithms compared over the 10 repetitions are shown in Figure 6. An
increase in the median best feasible value is seen around the third iteration
for the MMCU algorithm. This can happen when some of the runs do not find a
feasible solution at the beginning and are thus not taken into account. When
an additional feasible point arises, the median best feasible objective may
increase.

Figure 6: Convergence of $\mathbb{E}_{\mathbf{U}}f$ for the best feasible
point of the 4 optimization algorithms on the 2-dimensional Problem (32).
Median and 25% and 75% quartiles estimated from 10 repetitions.

A complementary view on the convergence expressed in terms of the distance
between the current best location and the true optimum is provided in Figure
7. The reference (true) optimum was found by exhaustive enumeration on the
true functions and constraints.

Figure 7: Distance between the current best point and the optimum as a
function of the iteration, for the 4 compared BO algorithms on the
2-dimensional Problem (32). Median and 25% and 75% quartiles estimated from 10
repetitions.

The best mean objectives of feasible solutions found by the algorithms over
the 10 repetitions are summarized in Figure 8.

Figure 8: Summary statistics for the best mean objective,
$\mathbb{E}_{\mathbf{U}}f$, of feasible points obtained by the 4 optimizers on
the 2-dimensional Problem (32), 10 repetitions of the runs.

For this test case, the 3 methods introduced here, SMCS, MMCU and MMCS,
converge faster during the first iterations than the reference method (REF).
It is also observed that constraint selection, as exercised by SMCS and MMCS,
accelerates the initial convergence (MMCS is better than MMCU). The effect of
constraint selection is detailed in Figure 9 which shows an histogram of the
frequency with which each constraint is called. By construction, both the REF
and MMCU approaches evaluate $g_{1}(\cdot)$ and $g_{2}(\cdot)$ the same number
of times, i.e., they both represent 50% of the overall constraints calls. On
the contrary, both the SMCS and MMCS approaches evaluate $g_{1}(\cdot)$ more
often than $g_{2}(\cdot)$, 70% against about 30% of the time.

Figure 9: Proportions of $g_{1}(\cdot)$ and $g_{2}(\cdot)$ evaluations made by
the 4 algorithms on the 2-dimensional Problem (32), averaged over 10
repetitions.

An explanation for why the infill procedure prioritizes $g_{1}(\cdot)$ over
$g_{2}(\cdot)$ is found in Figure 10 which shows the profiles of
$\mathbb{P}_{U}(g_{1}(\mathbf{x},\mathbf{u})\leq 0)$,
$\mathbb{P}_{U}(g_{2}(\mathbf{x},\mathbf{u})\leq 0)$ and
$\mathbb{E}_{U}(f(\mathbf{x},\mathbf{u})\leq 0)$. It is seen that at the
constrained optimum $\mathbf{x}^{\star}$ (marked by the vertical grey dashed
line), $g_{1}(\cdot)$ is saturated while $g_{2}(\cdot)$ is above the threshold
$1-\alpha$ (marked by the horizontal grey dashed line). As a consequence, the
surrogate of $g_{1}(\cdot)$ is required to be more precise than that of
$g_{2}(\cdot)$ to accurately assess feasibility as measured, in our
approaches, by the variance of the feasibility probability through $S_{g}$.

Figure 10: Profiles of $P_{U}(g_{1}(\mathbf{x},\mathbf{u})\leq 0)$,
$P_{U}(g_{2}(\mathbf{x},\mathbf{u})\leq 0)$ and
$E_{U}(f(\mathbf{x},\mathbf{u}))$. The optimum $\mathbf{x}^{\star}$ is
indicated by the vertical grey dashed line. $E_{U}(f(\mathbf{x},\mathbf{u}))$
is rescaled for readability.

Later in the run, as seen in Figures 6 and 7, the multi-output model enables
to better refine the best solution, thus MMCS and MMCU consistently converge
to feasible points with lower mean objective than the algorithms with
independent constraint models, REF and SMCS.

The benefit of selecting independent $\mathbf{u}$’s for the evaluation of the
objective and the constraint functions is illustrated in Figure 11 where the
DoEs generated by one execution of the different algorithms is plotted. It is
observed that the algorithms that sample independently $u$ for the objective
and the constraint functions (SMCS and MMCS) have the objective function
samples $u_{f}$ that explores the entire search space while the constraints
samples $u_{g}$ cluster on the boundaries. Instead, by construction,
$u_{f}=u_{g}$ for REF and MMCU. For all algorithms, the distribution of the
variables $x$ has 2 clusters, one in the optimal region, $x^{\star}\approx
0.17$, and one where the second, inactive constraint,
$0.95-\left(\mathbb{P}_{U}(g_{2}(x,U)<0\right)$, crosses 0 which happens for
$x\approx 0.28$ (cf. Figure 10).

Figure 11:  DoE’s generated by one execution of the different algorithms on
the 2-dimensional Problem (32). $(x^{i},u^{i})~{},~{}i=1,\ldots,40$ for the
objective function (left column) and the constraints (right column), for REF,
SMCS, MMCU and MMCS from row 1 to 4, respectively.

###  7.3 4-dimensional analytical test case

Our second test case has 2 design variables, 2 random variables, i.e., the
models are built in 4 dimensions, and 2 constraints. It was already introduced
as Problem (24). It is also a modified version of a benchmark found in [23],
with the addition of a second constraint.

The GP models of either the constraints or the objective functions are
initialized with a data set of 30 points. Subsequently, a budget of 160
constraint evaluations is allocated to the optimization which, divided by 2
since there are 2 constraints, make 80 iterations. The convergence rate of the
4 algorithms averaged over the 10 repetitions are shown in Figure 12.

Figure 12: Convergence rate (median and interquartile interval estimated from
10 runs) of the 4 BO algorithms on the 4-dimensional analytical test case
(24).

Summary statistics of the objective,
$\mathbb{E}_{\mathbf{U}}f(\widehat{\mathbf{x}^{\star}},\mathbf{U})$ at the
best feasible point $\widehat{\mathbf{x}^{\star}}$ found by the algorithms
over the 10 repetitions are provided in Figure 13.

Figure 13: Summary statistics over 10 repetitions of the objective,
$\mathbb{E}_{\mathbf{U}}f(\widehat{\mathbf{x}^{\star}},\mathbf{U})$ at the
best feasible point $\widehat{\mathbf{x}^{\star}}$ obtained by the 4 BO
algorithms on the 4-dimensional analytical test case (24).

These results show that, similarly to the 2-dimensional test case, the three
proposed variants (SMCS, MMCU, MMCS) outperform the reference algorithm (REF),
both in terms of convergence speed and final performance of the solution while
the interquartile distances remain inferior to that of REF.

Figure 14: Distance (median and interquartile interval) between the best
solution and the reference solution obtained by the 4 BO algorithms on the 4-d
analytical test case (24).

The convergence rate of the algorithms is expressed in terms of Euclidean
distance between the current best point and the reference optimum in Figure
14. This reference optimum was found by exhaustive enumeration on the true
functions and constraints. With this metric, the difference in performance
between the various methods is less noticeable, but the proposed variants
still show a faster convergence rate compared to REF. The MMCS variant
provides the best results according to both metrics of comparison.

In order to further analyze the difference between the methods, Figure 15
shows how the evaluations are distributed between the two constraints. Whereas
the SMCS variant evaluates more often, on the average of the 10 runs, $g_{1}$
than $g_{2}$, the MMCS variant seems to evenly calculate the two constraints.
The number of evaluations is detailed for each run in Figure 15(b). It can be
seen that the apparent equal evaluations of MMCS are due to an averaging
effect: the MMCS focuses its computational effort on either one of the two
constraints depending on the run. This is because the two constraints are
strongly correlated and both are active. Therefore, either one of the two
constraints provides a useful information from an optimization perspective
which can be modeled by the multi-output GPs and exploited by the constraint
selection mechanism.

(a) Proportions of $g_{1}(\cdot)$ and $g_{2}(\cdot)$ evaluations averaged over
10 repetitions.

(b) Details of the number of calls to $g_{1}(\cdot)$ and $g_{2}(\cdot)$ for
each of the 10 repetitions.

Figure 15: Number of calls to $g_{1}(\cdot)$ vs. $g_{2}(\cdot)$ for the REF,
SMCS, MMCU and MMCS algorithms, 4-dimensional Problem (24).

###  7.4 Industrial test case: compressor rotor design

Our last test case is the design of an airplane engine compressor rotor. More
specifically, the NASA rotor 37, an axial flow compressor originally designed
and experimentally tested by Reid and Moore in 1978 [71, 72] is considered.
The objective of the design is to maximize the compressor polytropic
efficiency. The design concerns the shape of the rotor, as described by 5
quantities sketched in Figure 16: the chord, the maximum thickness, the
location of the maximum thickness, the pitch and the sweep angles. Each shape
feature evolves along the rotor axis and is parameterized with the help of a
spline function with 4 control points, thus resulting in 20 actual design
variables.

Figure 16: Schematic representation of the rotor test case and its design
variables, courtesy of Safran Tech.

Additionally, the rotor performance depends on 7 random parameters
$\mathbf{U}$: the tip gap and the blade rugosity (manufacturing
uncertainties), the inflow pressure, the temperature and the azimuthal
momentum (inflow uncertainties), the flow rate and the rotation speed
(operational uncertainties). Finally, the optimization problem is subject to 5
constraints related to the inlet and outlet relative flow angles
($g_{1},g_{2}$), the deceleration of the flow to prevent separation ($g_{3}$),
the Lieblin coefficient to avoid excessive loading ($g_{4}$) and the Mach
number of the rotor tip ($g_{5}$).

The rotor design problem is formulated as follows:

$\min_{\mathbf{x}}\mathbb{E}_{U}[f(\mathbf{x},\mathbf{U})]\ \ s.t.\ \
\mathbb{P}(g_{p}(\mathbf{x},\mathbf{U})\leq 0,p=1,\ldots,5)\geq 0.9$ (33)

with $\mathbf{x}\in\mathcal{S_{X}}\subset\mathbb{R}^{20}$ and
$\mathbf{U}\sim\mathcal{U}([\text{LB},\text{UB}]^{7})$, where LB and UB
represent the lower and upper bounds of the uniformly distributed random
parameters.

With 27 variables and parameters, and 5 constraints, the joint space of the
rotor design problem is considerably larger than those of the previous
analytical test cases. Robust optimizations are computationally much more
expensive. As a result, various adjustments in the experiments have been
deemed necessary. First, the reliability in constraint satisfaction has been
set to 90% ($\alpha=0.1$). Second, only a single compressor optimization is
performed, i.e., the tests are not repeated. The results provided hereafter
should be taken with some caution as the variability caused by the initial DoE
and the auxiliary optimizations (GP learning, acquisition criteria
optimizations) is not studied. Third, only two of the four methods are
considered: REF and MMCS.

The GP models are trained with initial DoEs of 100 data samples per function.
With 5 constraints, this results in the initial multi-output GP of MMCS being
conditioned on 500 data samples. Subsequently, 500 constraint functions
evaluations are performed for the actual optimization (i.e., 100 iterations).
The convergences of the REF and MMCS algorithms are plotted in Figure 17.

Figure 17: Convergence of the best feasible $\mathbb{E}_{\mathbf{U}}f$ for the
REF and MMCS algorithms, compressor rotor design test case (1 execution).

This result shows that the MMCS variant, in green, with its multi-output
constraints model and constraint selection, outperforms the reference (REF)
algorithm, in red, as it provides a faster initial convergence speed as well
as a better final mean objective. In order to further investigate the effects
of selecting the single most informative constraint at every iteration, the
number of evaluations of each constraint is recorded in Table 1.

Constraint | $g_{1}$ | $g_{2}$ | $g_{3}$ | $g_{4}$ | $g_{5}$
---|---|---|---|---|---
N∘ evaluations REF | 100 | 100 | 100 | 100 | 100
N∘ evaluations MMCS | 103 | 202 | 17 | 98 | 80
Table 1: Number of evaluations of each constraint in the compressor rotor
design test case carried out by the REF and MMCS methods.

It is seen in the Table that the MMCS focuses a larger part of the
computational budget on some constraints, in particular $g_{2}$ and to a
smaller extent $g_{1}$, and tends to neglect other constraints (e.g.,
$g_{3}$).

##  8 Conclusion and perspectives

In this paper, three compatible extensions to the robust Bayesian optimization
algorithm were proposed in order to solve computationally intensive chance
constrained problems. These extensions are based i) on the creation of a
coupled model of the constraints, ii) on the independent selection of the
constraint and the associated random parameter where it will be evaluated.
This was made possible thanks to a decomposition of the acquisition criterion
and to an output-as-input scheme. The results obtained on both academic and an
engineering design test cases show that the two changes contribute to a
convergence which is both more accurate and more stable with respect to the
initial data set when compared to the original method. We therefore think that
coupled constraint models and constraint selection are key features of BO
methods when dealing with costly design problems affected by the presence of
random variables.

Additional testing might be necessary to better assess the effects and
possible limitations of the proposed approaches. For example, it could be
interesting to assess whether relying on multi-output models can induce a loss
of performance when the constraints are not correlated.

From a methodological perspective the MM** algorithms should be further
developed to allow the identification and the processing of groups of
correlated constraints. Finally, another continuation to this work would be to
model the correlation between the objective and the constraint functions, with
the anticipated challenge of having to define an acquisition criterion fitted
to this assumption as well as computationally tractable.

#### Acknowledgements:

This work was partly supported by the OQUAIDO research Chair in Applied
Mathematics.

## References

  * [1] Jack Kiefer, Jacob Wolfowitz, et al.  Stochastic estimation of the maximum of a regression function.  The Annals of Mathematical Statistics, 23(3):462–466, 1952.
  * [2] Kenneth Joseph Arrow, Hirofumi Azawa, Leonid Hurwicz, and Hirofumi Uzawa.  Studies in linear and non-linear programming, volume 2.  Stanford University Press, 1958.
  * [3] James C Spall.  Introduction to stochastic search and optimization: estimation, simulation, and control, volume 65.  John Wiley & Sons, 2005.
  * [4] Laetitia Andrieu, Guy Cohen, and Felisa J Vázquez-Abad.  Gradient-based simulation optimization under probability constraints.  European Journal of Operational Research, 212(2):345–351, 2011\.
  * [5] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.
  * [6] Dimo Brockhoff, Anne Auger, Nikolaus Hansen, Dirk V. Arnold, and Tim Hohm.  Mirrored sampling and sequential selection for evolution strategies.  In PPSN 2010, pages 11–21. Springer, 2010.
  * [7] Ilya Loshchilov, Marc Schoenauer, and Michèle Sebag.  Self-adaptive surrogate-assisted covariance matrix adaptation evolution strategy.  In Proc. of the 14th GECCO, pages 321–328. ACM, 2012.
  * [8] Anne Auger, Nikolaus Hansen, Verena Heidrich-Meisner, Olaf Mersmann, Petr Posik, and Mike Preuss.  3rd GECCO workshop for real-parameter optimization.  https://coco.gforge.inria.fr/doku.php?id=bbob-2012-results, 2012\.
  * [9] Virginie Gabrel, Cécile Murat, and Aurélie Thiele.  Recent advances in robust optimization: An overview.  European journal of operational research, 235(3):471–483, 2014\.
  * [10] Frank J Fabozzi, Petter N Kolm, Dessislava A Pachamanova, and Sergio M Focardi.  Robust portfolio optimization and management.  John Wiley & Sons, 2007.
  * [11] Byeng D Youn, Christian Gogu, and Raphael T Haftka.  Structural and multidisciplinary optimization–special issue.  Structural and Multidisciplinary Optimization, 54(6):1365–1366, 2016.
  * [12] Constantine Caramanis, Shie Mannor, and Huan Xu.  14 robust optimization in machine learning.  Optimization for machine learning, page 369, 2012.
  * [13] Aharon Ben-Tal, Laurent El Ghaoui, and Arkadi Nemirovski.  Robust optimization, volume 28.  Princeton University Press, 2009.
  * [14] Arkadi Nemirovski.  On safe tractable approximations of chance constraints.  European Journal of Operational Research, 219(3):707–718, 2012\.
  * [15] Marcos A Valdebenito and Gerhart I Schuëller.  A survey on approaches for reliability-based optimization.  Structural and Multidisciplinary Optimization, 42(5):645–663, 2010\.
  * [16] Vincent Dubourg, Bruno Sudret, and Jean-Marc Bourinet.  Reliability-based design optimization using kriging surrogates and subset simulation.  Structural and Multidisciplinary Optimization, 44(5):673–690, 2011\.
  * [17] Maliki Moustapha, Bruno Sudret, Jean-Marc Bourinet, and Benoît Guillaume.  Quantile-based optimization under uncertainties using adaptive kriging surrogate models.  Structural and multidisciplinary optimization, 54(6):1403–1421, 2016.
  * [18] Donald R Jones, Matthias Schonlau, and William J Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global optimization, 13(4):455–492, 1998.
  * [19] Jonas Mockus.  Bayesian approach to global optimization: theory and applications, volume 37.  Springer Science & Business Media, 2012.
  * [20] Carl Edward Rasmussen.  Gaussian processes in machine learning.  In Summer school on machine learning, pages 63–71. Springer, 2003\.
  * [21] Justin J Beland and Prasanth B Nair.  Bayesian optimization under uncertainty.  In NIPS BayesOpt 2017 workshop, 2017.
  * [22] Alexis Cousin, Josselin Garnier, Martin Guiton, and Miguel Zuniga.  Chance constraint optimization of a complex system: Application to the fatigue design of a floating offshore wind turbine mooring system.  In WCCM-ECCOMAS2020, 2020.
  * [23] Reda El Amri, Rodolphe Le Riche, Céline Helbert, Christophette Blanchet-Scalliet, and Sébastien Da Veiga.  A sampling criterion for constrained bayesian optimization with uncertainties.  https://arxiv.org/abs/2103.05706, 2021.  arXiv preprint no. 2103.05706.
  * [24] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger.  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  Siam Review, 60(3):550–591, 2018.
  * [25] M Giselle Fernández-Godino, Chanyoung Park, Nam-Ho Kim, and Raphael T Haftka.  Review of multi-fidelity models.  arXiv preprint arXiv:1609.07196, 2016.
  * [26] Marc C Kennedy and Anthony O’Hagan.  Predicting the output from a complex computer code when fast approximations are available.  Biometrika, 87(1):1–13, 2000.
  * [27] Paris Perdikaris, Maziar Raissi, Andreas Damianou, Neil D Lawrence, and George Em Karniadakis.  Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2198):20160751, 2017.
  * [28] Loic Le Gratiet and Josselin Garnier.  Recursive co-kriging model for design of computer experiments with multiple levels of fidelity.  International Journal for Uncertainty Quantification, 4(5), 2014\.
  * [29] Chris Williams, Edwin V Bonilla, and Kian M Chai.  Multi-task gaussian process prediction.  Advances in neural information processing systems, pages 153–160, 2007.
  * [30] Michel Goulard and Marc Voltz.  Linear coregionalization model: tools for estimation and choice of cross-variogram matrix.  Mathematical Geology, 24(3):269–286, 1992.
  * [31] Mauricio A Alvarez, Lorenzo Rosasco, and Neil D Lawrence.  Kernels for vector-valued functions: A review.  arXiv preprint arXiv:1106.6251, 2011.
  * [32] Loïc Brevault, Mathieu Balesdent, and Ali Hebbal.  Overview of gaussian process based multi-fidelity techniques with variable relationship between fidelities, application to aerospace systems.  Aerospace Science and Technology, 107:106339, 2020.
  * [33] Kurt Cutajar, Mark Pullin, Andreas Damianou, Neil Lawrence, and Javier González.  Deep gaussian processes for multi-fidelity modeling.  arXiv preprint arXiv:1903.07320, 2019.
  * [34] Amar Shah and Zoubin Ghahramani.  Pareto frontier learning with expensive correlated objectives.  In International Conference on Machine Learning, pages 1919–1927. PMLR, 2016.
  * [35] Seyede Fatemeh Ghoreishi and Douglas Allaire.  Multi-information source constrained bayesian optimization.  Structural and Multidisciplinary Optimization, 59(3):977–991, 2019\.
  * [36] Mostafa Meliani, Nathalie Bartoli, Thierry Lefebvre, Mohamed-Amine Bouhlel, Joaquim RRA Martins, and Joseph Morlier.  Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design.  In AIAA aviation 2019 forum, page 3236, 2019.
  * [37] Christopher James Brooks, AIJ Forrester, AJ Keane, and S Shahpar.  Multi-fidelity design optimisation of a transonic compressor rotor.  2011\.
  * [38] Nicolas Garland, Rodolphe Le Riche, Yann Richet, and Nicolas Durrande.  Multi-fidelity for MDO using gaussian processes.  In Aerospace system analysis and optimization in uncertainty, pages 295–320. Springer, 2020.
  * [39] Siyu Tao, Anton Van Beek, Daniel W Apley, and Wei Chen.  Multi-model bayesian optimization for simulation-based design.  Journal of Mechanical Design, 143(11):111701, 2021.
  * [40] Matthias Schonlau, William J Welch, and Donald R Jones.  Global versus local search in constrained optimization of computer models.  Lecture Notes-Monograph Series, pages 11–25, 1998.
  * [41] Michael J Sasena, Panos Papalambros, and Pierre Goovaerts.  Exploration of metamodeling sampling criteria for constrained global optimization.  Engineering optimization, 34(3):263–278, 2002.
  * [42] Victor Picheny, Robert B Gramacy, Stefan Wild, and Sebastien Le Digabel.  Bayesian optimization under mixed constraints with a slack-variable augmented Lagrangian.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
  * [43] Charles Audet, J Dennis, Douglas Moore, Andrew Booker, and Paul Frank.  A surrogate-model-based method for constrained optimization.  In 8th symposium on multidisciplinary analysis and optimization, page 4891, 2000.
  * [44] Rodolphe Le Riche, Victor Picheny, Andre Meyer, Nam-Ho Kim, and David Ginsbourger.  Gears design with shape uncertainties using controlled monte carlo simulations and kriging.  In 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference 17th AIAA/ASME/AHS Adaptive Structures Conference 11th AIAA No, page 2257, 2009.
  * [45] Victor Picheny, Tobias Wagner, and David Ginsbourger.  A benchmark of kriging-based infill criteria for noisy optimization.  Structural and Multidisciplinary Optimization, 48(3):607–626, 2013\.
  * [46] Victor Picheny, David Ginsbourger, Yann Richet, and Gregory Caplin.  Quantile-based optimization of noisy computer experiments with tunable precision.  Technometrics, 55(1):2–13, 2013.
  * [47] Sait Cakmak, Raul Astudillo, Peter Frazier, and Enlu Zhou.  Bayesian optimization of risk measures.  arXiv preprint arXiv:2007.05554, 2020.
  * [48] Christian Sabater, Olivier Le Maitre, Pietro Marco Congedo, and Stefan Görtz.  A bayesian approach for quantile optimization problems with high-dimensional uncertainty sources.  Computer Methods in Applied Mechanics and Engineering, 376:113632, 2021.
  * [49] Warren Scott, Peter Frazier, and Warren Powell.  The correlated knowledge gradient for simulation optimization of continuous parameters using gaussian process regression.  SIAM Journal on Optimization, 21(3):996–1026, 2011.
  * [50] Léonard Torossian, Victor Picheny, and Nicolas Durrande.  Bayesian quantile and expectile optimisation.  arXiv preprint arXiv:2001.04833, 2020.
  * [51] Mélina Ribaud, Christophette Blanchet-Scalliet, Céline Helbert, and Frédéric Gillot.  Robust optimization: A kriging-based multi-objective optimization approach.  Reliability Engineering and System Safety, 200:106913, 2020.
  * [52] Janis Janusevskis and Rodolphe Le Riche.  Simultaneous kriging-based estimation and optimization of mean response.  Journal of Global Optimization, 2012.  DOI 10.1007/s10898-011-9836-5.
  * [53] J Zhang, AA Taflanidis, and JC Medina.  Sequential approximate optimization for design under uncertainty problems utilizing kriging metamodeling in augmented input space.  Computer Methods in Applied Mechanics and Engineering, 315:369–395, 2017.
  * [54] Morgane Menz, Christian Gogu, Sylvain Dubreuil, Nathalie Bartoli, and Jérôme Morio.  Adaptive coupling of reduced basis modeling and kriging based active learning methods for reliability analyses.  Reliability Engineering & System Safety, 196:106771, 2020.
  * [55] V Dubourg and Bruno Sudret.  Meta-model-based importance sampling for reliability sensitivity analysis.  Structural Safety, 49:27–36, 2014.
  * [56] Julien Bect, David Ginsbourger, Ling Li, Victor Picheny, and Emmanuel Vazquez.  Sequential design of computer experiments for the estimation of a probability of failure.  Statistics and Computing, 22(3):773–793, 2012.
  * [57] Cécile Mattrand, Pierre Beaurepaire, and Nicolas Gayton.  Adaptive kriging-based methods for failure probability evaluation: Focus on ak methods.  Mechanical Engineering in Uncertainties From Classical Approaches to Some Recent Developments, page 205, 2021.
  * [58] Maliki Moustapha and Bruno Sudret.  Surrogate-assisted reliability-based design optimization: a survey and a unified modular framework.  Structural and Multidisciplinary Optimization, pages 1–20, 2019\.
  * [59] Sylvain Lacaze and Samy Missoum.  Reliability-based design optimization using kriging and support vector machines.  In Proceedings of the 11th International conference on structural safety & reliability, New York, 2013.
  * [60] Xiaoke Li, Haobo Qiu, Zhenzhong Chen, Liang Gao, and Xinyu Shao.  A local kriging approximation method using mpp for reliability-based design optimization.  Computers & Structures, 162:102–115, 2016.
  * [61] Zhenzhong Chen, Siping Peng, Xiaoke Li, Haobo Qiu, Huadi Xiong, Liang Gao, and Peigen Li.  An important boundary sampling method for reliability-based design optimization using kriging model.  Structural and Multidisciplinary Optimization, 52(1):55–70, 2015\.
  * [62] Mathieu Balesdent, Loïc Brevault, Jérôme Morio, and Rudy Chocat.  Overview of Problem Formulations and Optimization Algorithms in the Presence of Uncertainty, pages 147–183.  Springer International Publishing, Cham, 2020.
  * [63] Gilles Pujol, Rodolphe Le Riche, Olivier Roustant, and Xavier Bay.  Optimisation multidisciplinaire en mécanique: Réduction de modèles, robustesse, fiabilité, réalisations logicielles, chapter L’incertitude en conception: formalisation, estimation.  Hermes Science Publications, Paris, 2009.
  * [64] Laura P Swiler, Patricia D Hough, Peter Qian, Xu Xu, Curtis Storlie, and Herbert Lee.  Surrogate models for mixed discrete-continuous variables.  In Constraint Programming and Decision Making, pages 181–202. Springer, 2014.
  * [65] Yulei Zhang and William I Notz.  Computer experiments with qualitative and quantitative variables: A review and reexamination.  Quality Engineering, 27(1):2–13, 2015.
  * [66] Julien Pelamatti, Loïc Brevault, Mathieu Balesdent, El-Ghazali Talbi, and Yannick Guerin.  Overview and comparison of gaussian process-based surrogate models for mixed continuous and discrete variables: Application on aerospace design problems.  In High-Performance Simulation-Based Optimization, pages 189–224. Springer, 2020.
  * [67] Qiang Zhou, Peter ZG Qian, and Shiyu Zhou.  A simple approach to emulation for computer models with qualitative and quantitative factors.  Technometrics, 53(3):266–273, 2011.
  * [68] Olivier Roustant, David Ginsbourger, and Yves Deville.  Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization.  2012\.
  * [69] Yves Deville, David Ginsbourger, Olivier Roustant Contributors, Nicolas Durrande, Maintainer Yves Deville, Depends Rcpp, Suggests DiceKriging, MASS Imports, and LinkingTo Rcpp.  Package ‘kergp’.  2015\.
  * [70] Michael JD Powell.  A direct search optimization method that models the objective and constraint functions by linear interpolation.  In Advances in optimization and numerical analysis, pages 51–67. Springer, 1994.
  * [71] Lonnie Reid and Royce D Moore.  Design and overall performance of four highly loaded, high speed inlet stages for an advanced high-pressure-ratio core compressor.  Technical report, 1978.  report NASA-TP-1337.
  * [72] John Dunham.  CFD validation for propulsion system components.  Technical report, Advisory Group for Aerospace Research and Development, Neuilly-sur-Seine (France), 1998.
  * [73] Alan Genz and Frank Bretz.  Computation of multivariate normal and t probabilities, volume 195\.  Springer Science & Business Media, 2009.

##  Appendix A Main notations and acronyms

(by alphabetical order)

  * BO : Bayesian Optimization.

  * $c(.),C(.)$ : $1-\alpha-$ the probability of satisfying all the constraints, resp. its approximation as a process. Negative values mean feasibility with confidence $1-\alpha$.

  * $\mathcal{D}^{(t)},\mathcal{D}_{f}^{(t)},\mathcal{D}_{g_{p}}^{(t)},\mathcal{D}_{\mathbf{g}}^{(t)}$ : Design of Experiment at iteration $t$, containing a collection of points $(\mathbf{x}^{i},\mathbf{u}^{i})$. The subscript $f$, $g_{p}$, $\mathbf{g}$ refer to the objective function, the $p$-th constraint or the joint constraints vector.

  * DoE : Design of Experiment.

  * I, FI, EI, EFI : Improvement, Feasible Improvement, Expected Improvement, Expected Feasible Improvement.

  * $f(.,.),F(.,.)$ : function $\mathcal{S_{X}}\times\mathcal{S_{U}}\to\mathbb{R}$ and associated GP.

  * $\mathbf{f}^{(t)}$ : vector of evaluations of the function $f$ evaluated at the DoE $\mathcal{D}^{(t)}$ or $\mathcal{D}_{f}^{(t)}$

  * $g_{p}(,)$ : $p$-th constraint function $\mathcal{S_{X}}\times\mathcal{S_{U}}\to\mathbb{R}$, $p=1,\ldots,l$.

  * $\mathbf{g}^{(t)},\mathbf{g}_{p}^{(t)}$ : vector of evaluations of the constraints functions $g_{1}(,),\ldots,g_{l},)$ (resp. only $g_{p}(,)$) evaluated at the DoE $\mathcal{D}^{(t)}$ or $\mathcal{D}_{\mathbf{g}}^{(t)}$ or $\mathcal{D}_{g_{p}}^{(t)}$.

  * $\mathbf{g}(.,.),\mathbf{G}(.,.)$ : vector of $l$ constraints and associated (not necessarily independent) GPs.

  * $GP$ : Gaussian Process.

  * $k(),\mathbf{K}()$ : covariance function and matrix of a GP. Additional subscripts indicate which GP is concerned.

  * $l$ : number of constraints.

  * $m(),\mathbf{m}()$ : mean of a GP or of a vector of GPs. Additional subscript indicates which GP is concerned.

  * MMCS : Multiple Model of the constraints and Constraint Selection algorithm.

  * MMCU : Multiple Model of the constraints and Common $\mathbf{u}$ algorithm.

  * REF : reference algorithm (independent constraints, same iterate for all constraints and objective function).

  * $S(.,.)$ : proxy to the one-step-ahead EFI variance.

  * $S_{f}(.,.)$ : one-step-ahead improvement variance.

  * $S_{g}(.,.)$ : one-step-ahead variance of the feasibility probability.

  * SMCS : Single Models of the constraints and Constraint Selection algorithm.

  * $x,x^{\prime},x^{\star},x_{targ}$ : vectors of optimization variables $\in\mathcal{S_{X}}\subset\mathbb{R}^{d}$.

  * $u,u^{\prime},\tilde{u}$ : vector of random parameters $\in\mathcal{S_{U}}\subset\mathbb{R}^{m}$.

  * $z(.),Z(.)$ : mean of $f(.,U)$ over $U$ and associated Gaussian process.

##  Appendix B Covariance kernel for a nominal input: the hypersphere
decomposition.

The hypersphere decomposition [67] is a possible choice to parameterize a
discrete covariance kernel. The underlying idea is to map each of the $l$
levels of the considered discrete variable onto a distinct point on the
surface of a $l$-dimensional hypersphere:

$\displaystyle\phi(z):\mathcal{S}_{z}\rightarrow\mathbb{R}^{l}$
$\displaystyle\phi(z=z_{m})=\sigma_{z}[b_{m,0},b_{m,1},\dots,b_{m,l}]^{\top}\quad\mbox{
for }m=1,\ldots,l$

where $b_{m,d}$ represents the $d$-th coordinate of the $m$-th discrete level
mapping, and is calculated as follows:

$\displaystyle b_{1,1}=1$ $\displaystyle
b_{m,d}=\cos\theta_{m,d}\prod_{k=1}^{d-1}\sin\theta_{m,k}$
$\displaystyle\mbox{ for }d=1,\dots,m-1$ $\displaystyle
b_{m,m}=\prod_{k=1}^{d-1}\sin\theta_{m,k}$ $\displaystyle\mbox{ for }m\neq 1$
$\displaystyle b_{m,d}=0$ $\displaystyle\mbox{ for }d\geq m\neq 1$

with $-\pi\leq\theta_{m,d}\leq\pi$. It can be noticed that in the equations
above, some of the mapping coordinates are arbitrarily set to 0. This allows
to avoid rotation indeterminacies (_i.e.,_ an infinite number of
hyperparameter sets characterizing the same covariance matrix), while also
reducing the number of parameters required to define the mapping. The
resulting kernel is then computed as the Euclidean scalar product between the
hypersphere mappings,

$k(z,z^{\prime})=\phi(z)^{\top}\phi(z^{\prime})~{}.$

The discrete kernel can then be characterized as an $l\times l$ symmetric
positive definite matrix $\mathbf{T}$ containing the covariance values between
the discrete variable levels computed as:

$\mathbf{T}=\mathbf{L}^{\top}\mathbf{L}$

where each element of $\mathbf{L}_{i,j}$ is computed as $b_{i,j}$,

$\mathbf{L}=\sigma_{z}\left[\begin{array}[]{c c c c c}1&0&\dots&\dots&0\\\
\cos\theta_{2,1}&\sin\theta_{2,1}&0&\dots&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
\cos\theta_{l,1}&\sin\theta_{l,1}\cos\theta_{l,2}&\dots&\cos\theta_{l,l-1}\prod_{d=1}^{l-2}\sin\theta_{l,d}&\prod_{d=1}^{l-1}\sin\theta_{l,d}\end{array}\right]~{}.$

In the output-as-input model of the constraints, Equation (13), the matrix
$\mathbf{T}$ contains the covariance terms related to the constraint index,
$k_{p}(i,j)=\mathbf{T}_{ij}~{},~{}i,j=\\{1,\ldots,l\\}$.

##  Appendix C Probability of feasibility with dependent constraints

The probability of satisfying the coupled constraints intervenes in Equation
(17) for the EFI acquisition criterion. It can be estimated with GPs as:

$\mathbb{P}(C^{(t)}(\mathbf{x})\leq
0)\approx\frac{1}{N}\sum\limits_{k=1}^{N}\mathbbm{1}_{\big{(}1-\alpha-\frac{1}{M}\sum\limits_{j=1}^{M}\mathbbm{1}_{\big{(}\mathbf{G}^{(t)}(\mathbf{x},\mathbf{u_{j}},\omega_{k})\leq\mathbf{0}\big{)}}\leq
0\big{)}}.$ (34)

In practice, a set of $M$ instances of the random variables $\mathbf{U}$ is
sampled from $\rho_{U}$. Subsequently, $N$ independent multi-output
trajectories of $\mathbf{G}(\mathbf{x},\cdot)$ are simulated at the
aforementioned sampled random variables. The probability of feasibility can
finally be computed by simply counting the number of trajectories for which
the ratio of samples associated to feasible constraints is larger than
$1-\alpha$. Note that a multi-output GP prediction of the constraint vector is
defined as feasible when all of its components are below their specific
threshold value (which is $0$ in this work).

In (18), the improvement is computed with respect to the incumbent best
feasible solution $z_{min}^{feas}$. However, the objective function mean is
not observed therefore $z_{min}^{feas}$ cannot be read in the data set used to
condition the GPs $F$ and $\mathbf{G}$. For this reason, in [23] the incumbent
best feasible solution is defined by taking into account the mean of the GP
$Z^{(t)}(\mathbf{x})$ and the expected value of the process
$C^{(t)}(\mathbf{x})$:

$z_{min}^{feas}=\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0~{}.$ (35)

Given that the Fubini condition holds since the value of $C(\cdot)$ is bounded
by definition, the expected value of the process $C$ can be written

$\displaystyle\mathbb{E}[C^{(t)}(\mathbf{x})]$ $\displaystyle=$ $\displaystyle
1-\alpha-\mathbb{E}_{\mathbf{U}}[\mathbb{E}[1_{\mathbf{G}(\mathbf{x},\mathbf{U})\leq\mathbf{0}}]]$
(36) $\displaystyle=$ $\displaystyle
1-\alpha-\int_{\mathbb{R}^{m}}\mathbb{P}(\mathbf{G}(\mathbf{x},\mathbf{u})\leq\mathbf{0}))\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}$
$\displaystyle=$ $\displaystyle
1-\alpha-\int_{\mathbb{R}^{m}}\Phi\left(\mathbf{0}-\mathbf{m}_{\mathbf{G}}^{(t)}(\mathbf{x},\mathbf{u}),\mathbf{K}_{\mathbf{G}}^{(t)}(\mathbf{x},\mathbf{u})\right)\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}$

where $\Phi(\cdot)$ is the cumulative distribution function of a multivariate
Gaussian distribution which, like the univariate version, is estimated
numerically [73]. The above Equation and Equation (35) are used to compute the
current feasible minimum. If Equation (35) yields no solution,
$m^{(t)}_{Z}(\mathbf{x})$ at the $\mathbf{x}$ providing the largest
probability of feasibility is taken as the incumbent optimal solution.

##  Appendix D Estimating the proxy of the one-step-ahead feasible improvement
variance

During the optimization process, the value of the coordinates
$\mathbf{u}^{t+1}$ of the point to be added to the training data set
$\mathcal{D}$ is computed by minimizing the proxy of the one-step-ahead
variance of the EFI at $\mathbf{x}_{targ}$ (see Equation (19)), which is
defined as:

$S(\mathbf{x}_{targ},\mathbf{\tilde{u}})=Var\big{(}I^{(t+1)}(x_{targ})\big{)}\int_{\mathbb{R}^{m}}Var\big{(}1_{\\{\mathbf{G}^{(t+1)}(x_{targ},\mathbf{u})\leq
0\\}}\big{)}\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u},\\\
=~{}Var(\big{(}z_{\min}^{\mbox{feas}}-Z^{(t+1)}(x_{targ})\big{)}^{+})\int_{\mathbb{R}^{m}}Var\big{(}1_{\\{\mathbf{G}^{(t+1)}(x_{targ},\mathbf{u})\leq
0\\}}\big{)}\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}.$ (37)

In this appendix, we recall some details about the computation of the first
term in the previous equation. An expression of the improvement variance
$Var\left(I^{(t+1)}(x_{targ})\right)$ which bears some similarities to the one
of expected improvement is given in [23] and can be expressed in terms of
probability and density functions of Gaussian distributions:

$\displaystyle EI^{(s)}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\big{(}z_{\min}^{\text{feas}}-m_{Z}^{(s)}(\mathbf{x})\big{)}\phi\bigg{(}\frac{z_{\min}^{\text{feas}}-m_{Z}^{(s)}(\mathbf{x})}{\sigma_{Z}^{(s)}(\mathbf{x})}\bigg{)}+\sigma_{Z}^{(s)}(\mathbf{x})\Phi\bigg{(}\frac{z_{\min}^{\text{feas}}-m_{Z}^{(s)}(\mathbf{x})}{\sigma_{Z}^{(s)}(\mathbf{x})}\bigg{)},$
$\displaystyle=$ $\displaystyle\psi_{EI}(m_{Z}^{(s)}(x),\sigma_{Z}^{(s)}(x)).$

$\displaystyle Var\left(I^{(s)}(\mathbf{x})\right)$ $\displaystyle=$
$\displaystyle
EI^{(s)}(\mathbf{x})\big{(}z_{\min}^{\text{feas}}-m_{Z}^{(s)}(\mathbf{x})-EI^{(s)}(\mathbf{x})\big{)}$
$\displaystyle+(\sigma_{Z}^{(s)})^{2}(\mathbf{x})\Phi\bigg{(}\frac{z_{\min}^{\text{feas}}-m_{Z}^{(s)}(\mathbf{x})}{\sigma_{Z}^{(s)}(\mathbf{x})}\bigg{)},$
$\displaystyle=$ $\displaystyle\psi_{VI}(m_{Z}^{(s)}(x),\sigma_{Z}^{(s)}(x)).$

At step $t+1$, the training data set $\mathcal{D}^{(t)}$ is enriched by
$(\tilde{\mathbf{x}},\tilde{u})$ on which the output
$f(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$ is unknown and represented by
$F(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$.

As a consequence, $Var\left(I^{(t+1)}(x_{targ})\right)$ cannot directly be
computed because of the randomness of $m_{Z}^{(t+1)}(x_{targ})$. Indeed,
$m_{Z}^{(t+1)}(x_{targ})$ is equal to

$E\left(Z(x)|F(\mathcal{D}^{(t)})=f^{(t)},F(\tilde{\mathbf{x}},\tilde{u})\right)$
and follows :

$m_{Z}^{(t+1)}(x_{targ})\sim\mathcal{N}\Bigg{(}m_{Z}^{(t)}(x_{targ}),\bigg{(}\frac{\int_{\mathbb{R}^{m}}k_{*F}^{(t)}(x_{targ},\mathbf{u};\tilde{\mathbf{x}},\mathbf{\tilde{u}})\rho_{\mathbf{U}}(\mathbf{u})d\mathbf{u}}{\sqrt{k_{*F}^{(t)}(\tilde{\mathbf{x}},\mathbf{\tilde{u}};\tilde{\mathbf{x}},\mathbf{\tilde{u}})}}\bigg{)}^{2}\Bigg{)}.\\\
$

We can note that $\sigma_{Z}^{(t+1)}(x_{targ})$ is not random as it depends
only on the location $(\tilde{\mathbf{x}},\mathbf{\tilde{u}})$ and not on the
function evaluation at this point. By applying the law of total variance, it
can be shown that:

$\displaystyle Var\left(I^{t+1}(x_{targ})\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[Var\left(\big{(}z_{\min}^{\text{feas}}-Z(x_{targ})\big{)}^{+}|m_{Z}^{(t+1)}(x_{targ})\right)\right]$
$\displaystyle+Var\left[\mathbb{E}\left(\big{(}z_{\min}^{\text{feas}}-Z(x_{targ})\big{)}^{+}|m_{Z}^{(t+1)}(x_{targ})\right)\right].$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\psi_{VI}(m_{Z}^{(t+1)}(x),\sigma_{Z}^{(t+1)}(x))\right]$
$\displaystyle+Var\left[\psi_{EI}(m_{Z}^{(t+1)}(x),\sigma_{Z}^{(t+1)}(x))\right].$

This calculation is performed numerically using samples of
$m_{Z}^{(t+1)}(x_{targ})$. For the sake of clarity, the reader is referred to
the previous work ([23]) for details regarding the implementation.

##  Appendix E Flow charts of the algorithms

This appendix contains the flowcharts of the four algorithms which are
implemented and compared in the body of the article. The REF, SMCS, MMCU and
MMCS methods are described in Algorithms 2, 3, 4 and 5, respectively. REF
stands for reference and was initially proposed in [23]. It is a method where
the GPs of the constraints and objective function are independent, and the
same pair $(\mathbf{x}^{t},\mathbf{u}^{t})$ is added to every GP at each
iteration. SMCS, which stands for Single Models of the constraints and
Constraint Selection, has independent GPs, like the REF algorithm, but only
one constraint is updated at each iteration. The random parameters of the
objective function and the selected constraint, $\mathbf{u}_{f}$ and
$\mathbf{u}_{g}$, are different. MMCU means Multiple Model of the constraints
and Common $\mathbf{u}$. The MMCU algorithm has a joint model of all the
constraints and the same iterate $(\mathbf{x}^{t},\mathbf{u}^{t})$ enriches
all GPs. Finally, MMCS is the acronym of Multiple Model of the constraints and
Constraint Selection. The MMCS algorithm relies on a joint model of the
constraints and identifies at each iteration a single constraint and the
associated random sample, $\mathbf{u}_{g}$, to carry out the next evaluation.

Algorithm 2 REF algorithm

Initialize $t_{init}$, $t_{max}$

Define the initial data set
$\mathcal{D}^{(t_{init})}=\\{(\mathbf{x}^{i},\mathbf{u}^{i}),\ i=1,\ldots
t_{init}\\}$

Evaluate the objective and the constraint functions :

$\mathbf{f}^{(t_{init})}=\left(f(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,f(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$\mathbf{g}_{1}^{(t_{init})}=\left(g_{1}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{1}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

⋮

$\mathbf{g}_{l}^{(t_{init})}=\left(g_{l}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{l}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$t\leftarrow t_{init}$

while $t\leq t_{max}$ do

Train the GPs $F$ and $G_{p}$ for $p=1,\ldots,l$ w.r.t. $\mathcal{D}^{(t)}$,
$\mathbf{f}^{(t)}$ and $\mathbf{g}_{p}^{(t)}$

Compute the targeted variable by maximizing the EFI (Equation (17)),

$\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}\mathbb{E}[FI^{(t)}(\mathbf{x})]$

Set $\mathbf{x}^{t+1}=\mathbf{x}_{targ}$

Compute the random parameter $\mathbf{u}^{(t+1)}$ as the minimizer

of the proxy to the one-step-ahead EFI variance (Equation (20)),

$\mathbf{u}^{t+1}=\arg\min_{\mathbf{u}}S(\mathbf{x}_{targ},\mathbf{u})$

$\mathcal{D}^{(t)}\leftarrow\\{\mathcal{D}^{(t)},(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{f},\mathbf{u}^{t+1}_{g})\\}$

$\mathbf{f}^{(t)}\leftarrow\\{\mathbf{f}^{(t)},f(\mathbf{x}^{t+1},\mathbf{u}^{t+1})\\}$

$\mathbf{g}_{p}^{(t)}\leftarrow\\{\mathbf{g}_{p}^{(t)},g_{p}(\mathbf{x}^{t+1},\mathbf{u}^{t+1})\\}$
for $p=1,\ldots,l$

$t\leftarrow t+1$

end while

return the best feasible point according to the means of the processes,

$\arg\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0$

Algorithm 3 SMCS algorithm

Initialize $t_{init}$, $t_{max}$

Define the initial data set
$\mathcal{D}^{(t_{init})}=\\{(\mathbf{x}^{i},\mathbf{u}^{i}),\ i=1,\ldots
t_{init}\\}$

Evaluate the objective and the constraint functions :

$\mathbf{f}^{(t_{init})}=\left(f(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,f(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$\mathbf{g}_{1}^{(t_{init})}=\left(g_{1}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{1}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

⋮

$\mathbf{g}_{l}^{(t_{init})}=\left(g_{l}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{l}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$t\leftarrow t_{init}$

while $t\leq t_{max}$ do

Train the GPs $F$ and $G_{p}$ for $p=1,\ldots,l$ w.r.t.
$(\mathcal{D}_{f}^{(t)},\mathbf{f}^{(t)})$ and
$(\mathcal{D}_{g_{p}}^{(t)},\mathbf{g}_{p}^{(t)})$

Compute the targeted variable by maximizing the EFI (Equation (17)),

$\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}\mathbb{E}[FI^{(t)}(\mathbf{x})]$

Set $\mathbf{x^{t+1}}=\mathbf{x}_{targ}$

Compute the random parameters of the objective function by minimizing

the one-step-ahead improvement variance (Equation (5.1))

$\mathbf{u}_{f}^{t+1}=\mbox{arg}\min_{\mathbf{u}}S_{f}(\mathbf{x}_{targ},\mathbf{u})$

Determine the next constraint to evaluate and its random parameter

by minimizing the one-step-ahead variance

of the feasibility probability (Equation (30))

$\\{\mathbf{u}_{g}^{t+1},p\\}=\mbox{arg}\min_{\mathbf{u},p^{\prime}}S_{g}(\mathbf{x}_{targ},\mathbf{u},p^{\prime})$

$\mathcal{D}_{f}^{(t)}\leftarrow\\{\mathcal{D}_{f}^{(t)},(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{f})\\}$
,
$\mathbf{f}^{(t)}\leftarrow\\{\mathbf{f}^{(t)},f(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{f})\\}$

$\mathcal{D}_{g_{p}}^{(t)}\leftarrow\\{\mathcal{D}_{g_{p}}^{(t)},(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{g})\\}$
,
$\mathbf{g}_{p}^{(t)}\leftarrow\\{\mathbf{g}_{p}^{(t)},g_{p}(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{g})\\}$

$t\leftarrow t+1$ , $t_{p}\leftarrow t_{p}+1$

end while

return the best feasible point according to the means of the processes,

$\arg\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0$

Algorithm 4 MMCU algorithm

Initialize $t_{init}$, $t_{max}$

Define the initial data set
$\mathcal{D}^{(t_{init})}=\\{(\mathbf{x}^{i},\mathbf{u}^{i}),\ i=1,\ldots
t_{init}\\}$

Evaluate the objective and the constraint functions :

$\mathbf{f}^{(t_{init})}=\left(f(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,f(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$\mathbf{g}_{1}^{(t_{init})}=\left(g_{1}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{1}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

⋮

$\mathbf{g}_{l}^{(t_{init})}=\left(g_{l}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{l}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$t\leftarrow t_{init}$

while $t\leq t_{max}$ do

Train the GPs $F$ w.r.t. $\left(\mathcal{D}^{(t)},\mathbf{f}^{(t)}\right)$ and

and $\mathbf{G}$ w.r.t.
$\left(\underbrace{\\{\mathcal{D}^{(t)},...,\mathcal{D}^{(t)}\\}}_{l\text{
times}},\\{\mathbf{g}_{1}^{(t)},\ldots,\mathbf{g}_{l}^{(t)}\\}\right)$

Compute the targeted variable by maximizing the EFI (Equation (17)),

$\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}\mathbb{E}[FI^{(t)}(\mathbf{x})]$

Set $\mathbf{x^{t+1}}=\mathbf{x}_{targ}$

Compute the random parameter $\mathbf{u}^{t+1}$ as the minimizer

of the proxy to the one-step-ahead EFI variance (Equation (20)),

$\mathbf{u}^{t+1}=\arg\min_{\mathbf{u}}S(\mathbf{x}_{targ},\mathbf{u})$

$\mathcal{D}^{(t)}\leftarrow\\{\mathcal{D}^{(t)},(\mathbf{x}^{t+1},\mathbf{u}^{t+1})\\}$

$\mathbf{f}^{(t)}\leftarrow\\{\mathbf{f}^{(t)},f(\mathbf{x}^{t+1},\mathbf{u}^{t+1})\\}$

$\mathbf{g}_{p}^{(t)}\leftarrow\\{\mathbf{g}_{p}^{(t)},g_{p}(\mathbf{x}^{t+1},\mathbf{u}^{t+1})\\}$
for $p=1,\ldots,l$

$t\leftarrow t+1$

end while

return the best feasible point according to the means of the processes,

$\arg\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0$

Algorithm 5 MMCS algorithm

Initialize $t_{init}$, $t_{max}$

Define the initial data set
$\mathcal{D}_{f}^{(t_{init})}=\\{(\mathbf{x}^{i},\mathbf{u}^{i}),\ i=1,\ldots
t_{init}\\}$

and
$\mathcal{D}_{\mathbf{g}}^{(t_{init})}=\\{(\mathbf{x}^{i},\mathbf{u}^{i},p)~{},~{}i=1,\ldots
t_{init}~{},~{}p=1,\ldots,l\\}$

Evaluate the objective and the constraint functions :

$\mathbf{f}^{(t_{init})}=\left(f(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,f(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$\mathbf{g}^{(t_{init})}=\left(g_{1}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{1}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}}),\ldots,g_{l}(\mathbf{x}^{1},\mathbf{u}^{1}),\ldots,g_{l}(\mathbf{x}^{t_{init}},\mathbf{u}^{t_{init}})\right)$

$t\leftarrow t_{init}$

while $t\leq t_{max}$ do

Train the GPs $F$ w.r.t. $(\mathcal{D}_{f}^{(t)},\mathbf{f}^{(t)})$ and
$\mathbf{G}$ w.r.t. $(\mathcal{D}_{\mathbf{g}}^{(t)},\mathbf{g}^{(t)})$

Compute the targeted variable by maximizing the EFI (Equation (17)),

$\mathbf{x}_{targ}=\mbox{arg}\min_{\mathbf{x}}\mathbb{E}[FI^{(t)}(\mathbf{x})]$

Set $\mathbf{x^{t+1}}=\mathbf{x}_{targ}$

Compute the random parameters of the objective function by minimizing

the one-step-ahead improvement variance (Equation (5.1))

$\mathbf{u}_{f}^{t+1}=\mbox{arg}\min_{\mathbf{u}}S_{f}(\mathbf{x}_{targ},\mathbf{u})$

Determine the next constraint to evaluate and its random parameter

by minimizing the one-step-ahead variance

of the feasibility probability (Equation (30))

$\\{\mathbf{u}_{g}^{t+1},p\\}=\mbox{arg}\min_{\mathbf{u},p^{\prime}}S_{g}(\mathbf{x}_{targ},\mathbf{u},p^{\prime})$

$\mathcal{D}_{f}^{(t)}\leftarrow\\{\mathcal{D}_{f}^{(t)},(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{f})\\}$
,
$\mathbf{f}^{(t)}\leftarrow\\{\mathbf{f}^{(t)},f(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{f})\\}$

$\mathcal{D}_{\mathbf{g}}^{(t)}\leftarrow\\{\mathcal{D}_{\mathbf{g}}^{(t)},(\mathbf{x}^{t+1},\mathbf{u}_{g}^{t+1},p)\\}$
,
$\mathbf{g}^{(t)}\leftarrow\\{\mathbf{g}^{(t)},g_{p}(\mathbf{x}^{t+1},\mathbf{u}^{t+1}_{g})\\}$

$t\leftarrow t+1$

end while

return the best feasible point according to the means of the processes,

$\arg\min_{\mathbf{x}}m^{(t)}_{Z}(\mathbf{x})~{}\mbox{s.t.}~{}\mathbb{E}[C^{(t)}(\mathbf{x})]\leq
0$