# Non-Parametric Measure Approximations for Constrained Multi-Objective Optimisation under Uncertainty

Mickael Rivier, Nassim Razaaly, Pietro Marco Congedo

► To cite this version:

Mickael Rivier, Nassim Razaaly, Pietro Marco Congedo. Non-Parametric Measure Approximations for Constrained Multi-Objective Optimisation under Uncertainty. International Journal for Numerical Methods in Engineering, 2023, 10.1002/nme.7403. hal-03781832

HAL Id: hal-03781832

https://inria.hal.science/hal-03781832

Submitted on 20 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-Parametric Measure Approximations for Constrained Multi-Objective Optimisation under Uncertainty

M. Rivier<sup>a,b</sup>, N. Razaaly<sup>a,c</sup>, P.M. Congedo<sup>a</sup>

<sup>a</sup>Inria, Centre des Mathématiques Appliquées, Ecole Polytechnique, IPP, 91120 Palaiseau, France

<sup>b</sup>ArianeGroup, Le Haillan, 33185 Cedex, France

<sup>c</sup>Institut Pprime, ISAE-ENSMA, 2 Bd des Frères Lumière, 86360 Chasseneuil-du-Poitou, France

## Abstract

In this paper, we propose non-parametric estimations of robustness and reliability measures approximation error, employed in the context of constrained multi-objective optimisation under uncertainty. These approximations with tunable accuracy permit to capture the Pareto front in a parsimonious way, and can be exploited within an adaptive refinement strategy.

First, we illustrate an efficient approach for obtaining joint representations of the robustness and reliability measures, allowing sharper discrimination of Pareto-optimal designs. A specific surrogate model of these objectives and constraints is then proposed to accelerate the optimisation process.

Secondly, we propose an adaptive refinement strategy, using these tunable accuracy approximations to drive the computational effort towards the computation of the optimal area. To this extent, an adapted Pareto dominance rule and Pareto optimal probability computation are formulated.

The performance of the proposed strategy is assessed on several analytical test-cases against classical approaches. We also illustrate the method on an engineering application, performing shape optimisation under uncertainty of an Organic Rankine Cycle turbine.

## Keywords:

Optimisation under uncertainty, Gaussian processes, Robustness and reliability, Organic Rankine cycle

# 1. Introduction

There has been a growing interest in optimisation with uncertainty in recent years [1, 2, 3, 4]. We focus here specifically on cost-efficient and derivative-free strategies for tackling multiple objectives and constraints in this context. Optimisation Under Uncertainty (OUU) can be roughly classified under two categories, Reliability-Based Design Optimisation (RBDO), which deals with probabilistic and worst-case feasibility constraints, and Robust Design Optimisation (RDO), where the deterministic objectives are replaced with averaged or worst-case ones, possibly in a multi-objective context (i.e. Taguchi optimisation).
In the context of multi-objective optimisation [5, 6], a popular derivative-free optimisa-
tion strategy for constrained multi-objective problems are Evolutionary Algorithms (EA).
A comparison of multi-objective EA strategies is proposed in [7], in the bi- or tri-objective
setting. When dealing with more objectives (≥ 4), one should refer to the many-objective
optimisation literature, notably reviewed in [8]. For more details, [9] gives a thorough de-
scription of both the theoretical foundations and practical metaheuristics for application to
engineering test cases.

In this work, we focus specifically on performing optimisation in the presence of uncer-
tainties. Such a setting has received significant attention from the optimisation community,
as can be seen in several reviews [10, 11, 12]. Notably, [10] proposed four types of uncertain-
ties: (i) noise, which is an intrinsic bias on the evaluation of the fitness function, coming,
for example, from imperfect sensors; (ii) robustness, which refers to specific parameters that
may vary from their nominal value between the optimisation procedure and the real-world
application; (iii) fitness approximation, when the true fitness value is not computable or very
expensive, and less accurate but cheaper representation is used instead; (iv) time-varying
fitness function, meaning that the optimum changes through time and that the optimisation
algorithm should tackle this variability.

Much work has been devoted to tackling noise within the optimisation process. This
field is usually referred to as Noisy Optimisation (NO). In this context, some parsimonious
techniques rely on the assumption that there exist either different levels of accuracy for
computing the Quantities of Interest. Namely, multi-level techniques [13, 14] usually involve
several mesh refinements for numerical simulation, and multi-fidelity approaches [15, 16]
use different modeling strategies. Recent works have also proposed adaptive approaches for
tackling problems where accuracy can be tuned online. This approach is particularly adapted
when the fitness function is the limit of intermediate approximations of increasing fidelities,
such as iterative solvers or Monte Carlo estimations. This concept has been proposed in [17]
with online computation time allocation in a Bayesian Optimisation (BO) context. In this
work, a similar online refinement strategy is employed to achieve high parsimony.

In a multi-objective setting, the presence of noise and cohabitation of different approx-
imation accuracies raises the need for extending the classical Pareto dominance to an un-
certain environment. The EA community has proposed techniques for problems featuring
uniform [18] or Gaussian noises [19, 20, 21]. Interval and worst-case approaches are tackled
in [22, 23, 24], within the classical NSGA-II and MOEA frameworks. Recent developments
[25, 26, 27, 28] extends these computations to non-parametric noise distributions through
the construction of histogram-based approximate densities, which permits fast computations
of the probability of dominance. A straightforward extension of the Pareto dominance rule
to interval uncertainty in the presence of constraints is proposed in [29] under the name of
Bounding-Boxes. This extension has been exploited in several papers [30, 31] within an adap-
tive strategy for performing multi-objective optimisation on objective functions computed
with tunable accuracy.

In this paper, the core idea of the approach is to tune the accuracy of statistical objec-
tives and constraints approximations to focus computational power on the most promising
designs. The first contribution is the computation of non-parametric representation of the robustness and reliability measures error, permitting to estimate the Pareto front in a par- simonious way. Specifically, we aim to avoid the restriction to independent uniform [31] or Gaussian [32, 33, 34, 35] distributions, which often rely on the so-called projected process technique. We propose to use empirical sampling-based distribution approximations, that should allow for more representative error estimations given the current information, and hopefully quicker convergence. An original approach is introduced for constructing these measure samples, which relies on a coupled space Gaussian Process (GP) surrogate model for drawing joint realisations of objectives and constraints over the whole design space. These approximated measures are treated as a random vector and are compared and ranked using the concept of probabilistic Pareto dominance [36, 26, 27]. The second main contribution of this work is an adaptive refinement strategy, where the accuracy is improved only on the most promising individuals. We assess the proposed approach on several (mono and bi- objective) optimisation problems to estimate its computational cost and accuracy in terms of average distance to the exact Pareto front. Finally, we tackle the problem of an Organic Rankine turbine shape optimisation with random operating conditions.

Section 2 illustrates the formulation of multi-objective optimisation problems under uncertainties treated here within a probabilistic setting. In Section 3, we introduce the numerical ingredients permitting to compute non-parametric measure approximations. Then, Section 4 illustrates the proposed general framework to solve the OUU problem. In Section 5, we analyse the performance of the proposed approach on some algebraic test cases by comparing it with known techniques in the literature. Finally, Section 6 is devoted to showing the applicability of the framework to complex engineering test cases with the shape optimisation of Organic Rankine cycles turbines.

## 2. Uncertainty-based optimisation problem

In this section, the context of this work is made explicit. The general formulation of the tackled uncertainty-based optimisation problem is given, and classical choices in terms of robustness and reliability measures are introduced. The proposed probabilistic setting is then given, which allows for comparing aleatory values in a constrained multi-objective context.

### 2.1. Optimisation

A deterministic constrained multi-objective optimisation problem can be written as follows:

$$ \begin{array}{ll} \text{minimise:} & f(\boldsymbol{x}), \\ \text{satisfying:} & g(\boldsymbol{x}) \le \mathbf{0}, \\ \text{by changing:} & \boldsymbol{x} \in \mathcal{X}, \end{array} \qquad (1) $$

with objectives $f(\boldsymbol{x}) \in \mathbb{R}^{m_1}$ and constraints $g(\boldsymbol{x}) \in \mathbb{R}^{m_2}$, and where the design vector $\boldsymbol{x}$ lives in $\mathcal{X} \subset \mathbb{R}^n$. Note that $g(\boldsymbol{x}) \le \mathbf{0}$ stands for $\forall i, g_i(\boldsymbol{x}) \le 0$.
The computation of $f$ and $g$ can also depend on a multitude of uncontrollable parameters (e.g. environmental, material, geometrical, ...), here denoted by $\xi \in \Xi$.

The most straightforward setting is to consider these uncertain parameters as aleatoric. This work is limited to this scope and does not deal with epistemic uncertainties. In this context, $\xi$, $f(x, \xi)$ and $g(x, \xi)$ are random variables and Equation (1) is not suitable anymore.

It is classical to reformulate this problem into a measure-based Optimisation Under Uncertainty (OUU) problem, where specific statistics of the random outputs $f(x, \xi)$ and $g(x, \xi)$ are plugged as objective and constraint functions.

## 2.2. Robustness and reliability measures

Writing $\rho_f$ and $\rho_g$ some statistical measures on the objective function $f$ and constraint $g$, the OUU problem can be formulated as follows:

$$ 
\begin{array}{ll}
\text{minimise:} & \rho_f(x), \\
\text{satisfying:} & \rho_g(x) \le 0, \\
\text{by changing:} & x \in \mathcal{X}.
\end{array} 
\qquad (2) $$

As raised in the introduction, using statistical measures as objectives is usually referred to as Robust Design Optimisation (RDO), and these statistics ($\rho_f$) are called robustness measures. Similarly, Reliability-Based Design Optimisation (RBDO) deals with so-called reliability measures $\rho_g$ as constraints.

In the following, the original functions $f$ and $g$ are gathered in a vector $q \in \mathbb{R}^{m_1+m_2}$, that refers to Quantities of Interest (QoI). Namely,

$$ \forall (x, \xi) \in \mathcal{X} \times \Xi, \quad q(x, \xi) = \begin{pmatrix} f(x, \xi) \\ g(x, \xi) \end{pmatrix}. \qquad (3) $$

Similarly, $\rho_f$ and $\rho_g$ are gathered in a single statistical measures vector $\rho \in \mathbb{R}^{m_1+m_2}$,

$$ \forall x \in \mathcal{X}, \quad \rho(x) = \begin{pmatrix} \rho_f(x) \\ \rho_g(x) \end{pmatrix}. \qquad (4) $$

The most classical robustness and reliability measures for problem (2) are the following:

* Expectation: $\rho(x) = \mathbb{E}_{\xi}[q(x, \xi)]$,

* Variance: $\rho(x) = \mathbb{V}_{\xi}[q(x, \xi)]$,

* Worst case: $\rho(x) = \max_{\xi}[q(x, \xi)]$

* Quantile associated with probability $p$: $\rho(x) = q_{\xi}^{p}[q(x, \xi)]$ (also called Value at Risk, VaR),

* Superquantile: $\rho(x) = \mathbb{E}_{\xi}[q(x, \xi) | q(x, \xi) \ge q_{\xi}^{p}[q(x, \xi)]]$ (also called Expected Value at Risk, EVaR).
Note that any combination of these may be used, and one has to choose the most relevant for a given application.

In practice, robust objectives usually are either expectation or worst-case measures. The latest can also be relaxed to quantile or superquantile measures.

As for reliability measures, they can impose hard constraint, by using worst-case measures, but are usually relaxed into probabilistic constraints, that ensures $g(x, \xi) \le 0$ with a given probability $p$. This can be readily obtained by using the $p$-quantile as a reliability measure.

Given the distribution of the random input variable $\xi$, and specific measures $\rho$, problem (2) can be solved with classical techniques. However, if the unitary evaluation of $q$ is expensive, $\rho$ becomes very tedious to compute. A compromise between the accuracy and computational cost of $\rho$ must be chosen.

Specifically, we propose to consider these estimated measures as random variables, with a given Probability Density Function (PDF) representing the associated inaccuracy. Next section presents this probabilistic setting in details.

## 2.3. Probabilistic setting

In this work, we aim to tackle multi-objective optimisation problems ($\le 3$ objectives) where the statistical measures $\rho$ used as objectives and constraint functions are not exactly known and regarded as random fields $\mathcal{P}$ over $\mathcal{X}$. In this context, the classical Pareto dominance does not hold and a probabilistic one must be defined. Several propositions exist in the literature [36, 26, 27], notably differing by the assumed distribution shape.

The classical Pareto dominance is described in the case of minimisation. For two vectors $\mathbf{a}, \mathbf{b} \in \mathbb{R}^m$,

$$ 
\begin{align} 
\mathbf{a} \succ \mathbf{b} (\mathbf{a} \text{ dominates } \mathbf{b}) &\iff \forall j \in [1, m], a_j \le b_j \quad \text{and} \nonumber \\ 
&\phantom{\iff} \exists j \in [1, m], a_j < b_j, \nonumber \\ 
\mathbf{a} \succcurlyeq \mathbf{b} (\mathbf{a} \text{ strictly dominates } \mathbf{b}) &\iff \forall j \in [1, m], a_j < b_j, \tag{5} \\ 
\mathbf{a} \sim \mathbf{b} (\mathbf{a} \text{ is indifferent to } \mathbf{b}) &\iff \mathbf{a} \not\succcurlyeq \mathbf{b} \text{ and } \mathbf{b} \not\succcurlyeq \mathbf{a}. \nonumber 
\end{align} 
$$

In the presence of constraints, we follow the constrained Pareto dominance proposed in [31]. For two vectors $\mathbf{a}, \mathbf{b} \in \mathbb{R}^m$, with $\mathbf{a}_f$ and $\mathbf{b}_f$ the objective values and $\mathbf{a}_g$ and $\mathbf{b}_g$ the constraint values,

$$ 
\begin{align} 
\mathbf{a} \succcurlyeq_c \mathbf{b} &\iff \mathbf{a} \in \mathcal{A} \text{ and } \mathbf{a}_f \succcurlyeq \mathbf{b}_f \text{ or } \mathbf{b} \in \mathcal{F}, \nonumber \\ 
\mathbf{a} \succcurlyeq_c \mathbf{b} &\iff \mathbf{a} \in \mathcal{A} \text{ and } \mathbf{a}_f \succcurlyeq \mathbf{b}_f \text{ or } \mathbf{b} \in \mathcal{F}, \tag{6} \\ 
\mathbf{a} \sim_c \mathbf{b} &\iff \mathbf{a} \not\succcurlyeq_c \mathbf{b} \text{ and } \mathbf{b} \not\succcurlyeq_c \mathbf{a}, \nonumber 
\end{align} 
$$

where $\mathcal{A}$ is the admissible set, where all constraints are satisfied, and $\mathcal{F}$ is the failure set, where at least one constraint is not verified. Using this definition, two individuals that both violate the constraints will mutually dominate each other. This allows to discard inefficient individuals even when no satisfying individual has yet been found.
Here, a probabilistic Pareto dominance between two designs $x_1$ and $x_2$ hence consists in computing the probability $\mathbb{P}_P[\mathbf{P}(x_1) \succ_c \mathbf{P}(x_2)]$ of a design to dominate the other. Recall that $\mathbf{P}$ is a random field, meaning that $\mathbf{P}(x_1)$ and $\mathbf{P}(x_2)$ are random variables. For better generality and readability, these are denoted $\mathbf{A}$ and $\mathbf{B}$ in the following.

Formally, the unconstrained probability is defined as:

$$ \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succ \mathbf{B}] = \int_{\mathbb{R}^m} \phi_{\mathbf{A}-\mathbf{B}}(\mathbf{y}) d\mathbf{y} = \Phi_{\mathbf{A}-\mathbf{B}}(\mathbf{0}), \quad (7) $$

where $\phi_{\mathbf{A}-\mathbf{B}}$ is the Probability Density Function (PDF) of $\mathbf{A}-\mathbf{B}$ and $\Phi_{\mathbf{A}-\mathbf{B}}$ is the associated Cumulative Density Function (CDF). An equivalent formula can be written using the joint PDF $\phi_{\mathbf{A},\mathbf{B}}$ of $\mathbf{A}$ and $\mathbf{B}$:

$$ \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succ \mathbf{B}] = \iint_{\mathbf{a}>\mathbf{b}} \phi_{\mathbf{A},\mathbf{B}}(\mathbf{a}, \mathbf{b}) d\mathbf{a} d\mathbf{b}. \quad (8) $$

As for the constrained case, definition (6) shows that for $\mathbf{a}$ to dominate $\mathbf{b}$, either $\mathbf{b}$ must be in the failure domain, or $\mathbf{a}$ must be in the admissible domain and dominate $\mathbf{b}$ according to the classical Pareto dominance rule. The integral in (8) can be split up depending on $\mathbf{b}$ belonging to the failure or admissible domain, and an additional dominance constraint can be added in the latest case. Hence, the domination probability reads:

$$ \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succ_c \mathbf{B}] = \int_{b \in \mathcal{F}} \phi_{\mathbf{B}}(\mathbf{b}) d\mathbf{b} + \iint_{\substack{\mathbf{a} \in \mathcal{A} \\ \mathbf{a} > \mathbf{b}}} \phi_{\mathbf{A},\mathbf{B}}(\mathbf{a}, \mathbf{b}) d\mathbf{a} d\mathbf{b}. \quad (9) $$

Note that this formulation naturally tackles the existence of dependency between the random variables $\mathbf{A}$ and $\mathbf{B}$, which greatly differs from the assumptions made in [26, 27].

We finally propose a probabilistic constrained Pareto dominance rule between two random vectors:

$$ \begin{align} \mathbf{A} \succ_c \mathbf{B} &\iff \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succ_c \mathbf{B}] = 1, \nonumber \\ \mathbf{A} \succcurlyeq_c \mathbf{B} &\iff \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succcurlyeq_c \mathbf{B}] = 1, \nonumber \\ \mathbf{A} \sim_c \mathbf{B} &\iff \mathbf{A} \not\succcurlyeq_c \mathbf{B} \text{ and } \mathbf{B} \not\succcurlyeq_c \mathbf{A}. \tag{10} \end{align} $$

This dominance rule generalises the constrained boxed Pareto dominance introduced in [31].

However, in the context of random variables with infinite support (e.g. Gaussian variable), the dominance probability cannot reach exactly 1. It has been proposed in [26] to relax this threshold. The $\epsilon$-relaxed probabilistic constrained Pareto dominance finally reads:

$$ \begin{align} \text{Domination: } \mathbf{A} \succ_c \mathbf{B} &\iff \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succ_c \mathbf{B}] \ge 1 - \epsilon, \nonumber \\ \text{Strict domination: } \mathbf{A} \succcurlyeq_c \mathbf{B} &\iff \mathbb{P}_{\mathbf{A},\mathbf{B}}[\mathbf{A} \succcurlyeq_c \mathbf{B}] \ge 1 - \epsilon, \tag{11} \\ \text{Indifference: } \mathbf{A} \sim_c \mathbf{B} &\iff \mathbf{A} \not\succcurlyeq_c \mathbf{B} \text{ and } \mathbf{B} \not\succcurlyeq_c \mathbf{A}. \nonumber \end{align} $$
The above formulation gives a general setting for comparing aleatory vectors, regardless of their distribution shape. In practice, $\epsilon$ is usually taken at $10^{-2}$. Next section proposes a sampling-based technique to get an accurate representation of the random field $P$ over the design space $\mathcal{X}$.

### 3. Non-parametric measure approximation

We aim to make as few assumptions as possible on the distribution of $P(\boldsymbol{x})$, for $\boldsymbol{x}$ in $\mathcal{X}$. To this extent, we only assume that realisations $P^{(i)}(\boldsymbol{x})$ can be drawn and exploit this empirical set of samples to characterise $P(\boldsymbol{x})$.

More precisely, because $P$ is a random field over the design space $\mathcal{X}$, we propose to draw realisations $P^{(i)}$ of the whole random field. This consists in jointly drawing $P^{(i)}(\boldsymbol{x})$ on a set of locations $\boldsymbol{x}$.

In this setting, dominance must only be verified between these joint samples. With a number $N_{samples}$ of realisations, the probability of dominance between two designs can be empirically estimated as follows:

$$ \mathbb{P}_{P}[P(\boldsymbol{x}_j) \not\subset_c P(\boldsymbol{x}_k)] \approx \frac{1}{N_{samples}} \sum_{i=1}^{N_{samples}} \mathbb{1}_{\{P^{(i)}(\boldsymbol{x}_j) \not\subset_c P^{(i)}(\boldsymbol{x}_k)\}}. \quad (12) $$

In practice, we tend to take $N_{samples}$ between $10^3$ and $10^4$.

The major task is to obtain realisations $P^{(i)}(\boldsymbol{x})$ of the robustness and reliability measures. To this extent, we assume that a probabilistic surrogate can be constructed on the original Quantities of Interest $q$ and draw realisations $P^{(i)}$ through measure estimations on functional realisations $q^{(i)}(\boldsymbol{x}, \cdot)$, simultaneously for every studied location $\boldsymbol{x} \in \mathcal{X}$. In this work, we use Gaussian process surrogate models, that are built using the GPy python package [37].

We first recall the main characteristics and formulas of Gaussian Process surrogate models. Then, an approximated approach for drawing these realisations of $P$ from a GP surrogate of $q$ in the coupled space $\mathcal{X} \times \Xi$ is proposed.

#### 3.1. Gaussian processes

This section gives a quick reminder of the GP setting, its interpretability and the main formulas that will be used in the following. Gaussian process surrogate modelling, also called Kriging, is a Bayesian kernel-based regression method [38].

The kernel setting, with the famous so-called kernel trick, allows to transpose most of the classical parametric (e.g. linear) approaches to the non-parametric world. In practice, basing on data points location $\boldsymbol{x}_i$, parametric features can be incorporated by replacing scalar product between the different parameters $\boldsymbol{x}_1^\top \boldsymbol{x}_2$ by a scalar product between chosen features $\phi(\boldsymbol{x}_1)^\top \phi(\boldsymbol{x}_2)$. Then, the kernel trick consists in choosing a kernel function $k$, and to use $k(\boldsymbol{x}_1, \boldsymbol{x}_2)$ in place of the previous scalar product. In this way, the features are not explicitly defined and can be of infinite number. Some kernels are called universal and engender a space of functions that is dense in the space of continuous function.
The Bayesian setting, in addition, tracks the distribution of all parameters (of infinite number with the kernel trick), and returns not only a predictive model of the Quantity of Interest, but its whole distribution at each location $\boldsymbol{x}$. In the context of Gaussian process, this distribution is fully determined by a predictive mean and a predictive variance for the QoI.

In practice, Gaussian processes can also simply be viewed as an infinite dimensional Gaussian distribution, which reduction to a finite set of points is also Gaussian. In this setting, the kernel plays the role of covariance function between points and the predictive mean and variance are the conditional mean and variances at each point $\boldsymbol{x}$. Given the training data $\{\boldsymbol{x}_{train}, \boldsymbol{y}_{train}\}$ where $\boldsymbol{y}_{train} = f(\boldsymbol{x}_{train}) + \epsilon$ and $\epsilon$ is a random Gaussian noise, they read:

$$ \mu_{GP}(\boldsymbol{x}) = \boldsymbol{k}_{*}^{\top}(K + \Delta)^{-1}\boldsymbol{y}_{train}, \\ \Sigma_{GP}(\boldsymbol{x}) = \boldsymbol{k}_{**} - \boldsymbol{k}_{*}^{\top}(K + \Delta)^{-1}\boldsymbol{k}_{*}, \qquad (13) $$

where $K$ is the covariance matrix between training points $\boldsymbol{x}_{train}$, $K = k(\boldsymbol{x}_{train}, \boldsymbol{x}_{train})$, $\boldsymbol{k}_{*}$ is the covariance $k(\boldsymbol{x}_{train}, \boldsymbol{x})$ between the new points $\boldsymbol{x}$ and the training points and $\boldsymbol{k}_{**}$ is the covariance between the new points, $k(\boldsymbol{x}, \boldsymbol{x})$. Matrix $\Delta$ refers to the so-called input noise or nugget effect. Usually, it regularises the problem with the same noise on each training point, $\Delta = \sigma^2 I_N$, with $N$ the number of training data. However, in the heteroscedastic case, one can readily assign a different noise to each training point, meaning $\Delta = \text{diag}(\{\sigma_i^2\}_{i=1}^N)$.

In this paper, hyperparameters of the kernel function are optimised by Maximum Likelihood. Because Automatic Relevance Determination (ARD) is exploited, a lengthscale parameter is assigned to each dimension. This heavily limits the use of these models when the problem features more than tens of dimensions.

Finally, one can draw functional realisations of a GP surrogate model, discretised on a finite set of points. A GP is a Gaussian random vector, so a realisation $\boldsymbol{f}^{(i)} \sim \mathcal{N}(\boldsymbol{\mu}_{GP}, \boldsymbol{\Sigma}_{GP})$ is drawn as follows:

$$ \boldsymbol{f}^{(i)} = \boldsymbol{\mu}_{GP} + L_{GP}Y^{(i)}, \quad Y \sim \mathcal{N}(\mathbf{0}, I_N), \qquad (14) $$

with $L_{GP}$ referring to the Cholesky decomposition of the covariance matrix, so that $\boldsymbol{\Sigma}_{GP} = L_{GP}L_{GP}^{\top}$.

## 3.2. Direct measure sampling in coupled space

We assume that a GP surrogate model can be built on $\boldsymbol{q}$ in the coupled design/uncertainty space. We here simply build a scalar GP for each output dimension. Based on training data $\{\boldsymbol{q}(\boldsymbol{x}_i, \boldsymbol{\xi}_i)\}_{i=1}^N$, the surrogate $\hat{\boldsymbol{q}}$ is constructed and yields an estimation of $\boldsymbol{q}$ at any location $(\boldsymbol{x}, \boldsymbol{\xi})$ in the coupled space $\mathcal{X} \times \Xi$.

This surrogate model can be exploited to obtain an estimation of the robustness and reliability measure through Monte Carlo integration. This integration is directly performed on the predictive mean $\hat{\boldsymbol{q}}$ of the GP. For example, with the expectation measure $\boldsymbol{\rho}(\boldsymbol{x}) =$
$\mathbb{E}_{\xi}[q(\mathbf{x}, \boldsymbol{\xi})]$, one can use the estimator :

$$ \tilde{\rho}(\mathbf{x}) = \frac{1}{N_{\xi}} \sum_{j=1}^{N_{\xi}} \hat{q}(\mathbf{x}, \boldsymbol{\xi}_j). \qquad (15) $$

More interestingly, one can draw realisations of the GP built on $\mathbf{q}$, which directly translates into realisations of $\mathbf{P}$ by integrating over the uncertain dimension.

For example, in the case of the expectation measure $\rho(\mathbf{x}) = \mathbb{E}_{\xi}[q(\mathbf{x}, \boldsymbol{\xi})]$, by evaluating $\mathbf{q}$ with $N_{\xi}$ locations in the uncertain space at $N_x$ different designs, the joint realisations of $\mathbf{P}$ over all these designs write as follows:

$$ \forall k \in [1, N_x], \quad \mathbf{P}^{(i)}(\mathbf{x}_k) = \frac{1}{N_{\xi}} \sum_{j=1}^{N_{\xi}} q^{(i)}(\mathbf{x}_k, \boldsymbol{\xi}_j), \qquad (16) $$

drawn simultaneously on a realisation $\mathbf{q}^{(i)}$ at all designs.

The same can be carried out with empirical estimators of other robustness and reliability measures.

However, one can notice that $\mathbf{q}^{(i)}$ is drawn at $N_x \times N_{\xi}$ locations, which can quickly reach a value of more than $10^5$. Looking back at Equation (14), it reveals that Cholesky decomposition should hence be performed on a matrix of size $10^5 \times 10^5$. The computational cost associated with such a task is prohibitive. This problem is quite common in the literature, and several approximations techniques have been proposed, such as a circulant embedding of the covariance matrix [39], regularisation [40] or the use of a Gibbs sampler for Gaussian Markov Random Field [41]. However, some of these approaches require equispaced points, and the others remain too expensive to allow substantial sampling of the $\mathbf{q}$ random field.

As raised in [42], the complexity is particularly high when the domain is much larger than the characteristic correlation length of the random field. Intuitively, when a value is drawn, it gives information in its surrounding at a distance proportional to the correlation length.

Based on this idea of local information, an inducing point strategy is proposed in the next section for drawing realisations $\mathbf{q}^{(i)}$ at an extremely high number of locations simultaneously. This strategy is presented in Section 3.2.1, and a quantitative quality indicator is proposed in Section 3.2.2 for use in an iterative strategy in Section 3.2.3.

### 3.2.1. Inducing points strategy

The central idea is to draw realisations $\mathbf{q}^{(i)}$ of the Gaussian random field on few well-spread points, from which a smooth functional realisation can be extended to the whole input space. These well-chosen points will be referred to as inducing points.

In the following, for conciseness, coordinates $(\mathbf{x}, \boldsymbol{\xi})$ in the coupled space are denoted with $\mathbf{z}$. With a set of training designs $\mathbf{z}_{train} = (\mathbf{x}_{train}, \boldsymbol{\xi}_{train})$ and a set of inducing points $\mathbf{z}_{ind} = (\mathbf{x}_{ind}, \boldsymbol{\xi}_{ind})$, we denote $\mathbf{z}_{tot} = \mathbf{z}_{train} \cup \mathbf{z}_{ind}$.

We assume that a GP has been constructed on $\mathbf{z}_{train}$, as plotted in Figure 1, on values $\mathbf{q}_{train}$ and with an optimised kernel $k_{train}$. Its predictive mean and covariance are written
244 $\mu_{train}$ and $\Sigma_{train}$. This surrogate model is the current estimation of $q$ over the whole coupled space.
245

<figure><img src="image_2.png" /><figcaption>Figure 1: Gaussian process model built on the training data.</figcaption></figure>

246 We aim to generate random realisations of this Gaussian field on an extensive set of points
247 $z_{draw}$. Direct Gaussian sampling would perform Cholesky decomposition on the matrix
248 $K_{draw} = \Sigma_{train}(z_{draw})$ of size $N_{draw}^2$, which quickly becomes unfeasible.

We propose to only draw realisations on the inducing points:

$$ q_{ind}^{(i)} = \mu_{train}(z_{ind}) + L_{ind}Y^{(i)}, \quad Y \sim \mathcal{N}(0, I_{N_{ind}}), \qquad (17) $$

249 with $L_{ind}$ the Cholesky decomposition of matrix $K_{ind} = \Sigma_{train}(z_{ind})$ of much smaller size
250 $N_{ind}^2$. Figure 2 gives an example of such drawn values at the inducing points.

<figure><img src="image_3.png" /><figcaption>Figure 2: Realisations drawn at the inducing points.</figcaption></figure>

Similarly to the locations $z_{tot}$, we gather training values $q_{train}$ and drawn ones $q_{ind}^{(i)}$ in $q_{tot}^{(i)}$. The realisations drawn at the inducing points are then extended to all locations $z_{draw}$.
by constructing a new GP on the $(z_{tot}, q_{tot}^{(i)})$ database. Note that the new GP kernel is not optimised as $k_{train}$ is directly used. Final realisations estimates $q_{draw}^{(i)}$ are obtained as follows:

$$q_{draw}^{(i)} = \mu_{tot}^{(i)}(z_{draw}) = k_{*}^{\top} K_{tot}^{-1} q_{tot}^{(i)}, \quad (18)$$

where $k_* = k_{train}(z_{tot}, z_{draw})$ and $K_{tot}^{-1}$ only need to be computed once as they do not depend on $i$. These realisations are represented in Figure 3, and the empirical $2-\sigma$ band of the random field is computed on $10^5$ samples. One can note that they fit well the real confidence band of the underlying GP.

Intuitively, this approach corresponds to parametrising the random field as a sum of kernels located at all training and inducing points. This is noticeable in the above Equation (18). By denoting $K_{tot}^{-1}q_{tot}^{(i)}$ as a vector $\beta^{(i)}$ of coefficients, the proposed realisations $q_{draw}^{(i)}$ are linear combinations of the kernel function $k_{train}$ located at each training and inducing point. Coefficients $\beta^{(i)}$ manage the randomness and kernel functions the spatial characteristic. The functional realisations write $k_{train}(z_{tot}, \cdot)^{\top}\beta^{(i)}$.

In practice, the proposed approach has a lot of similarities with using a Karhunen-Loève Expansion (KLE) of a Gaussian random field with Nyström approximation of the basis functions. Computational details are given in Appendix A.

<figure><img src="image_4.png" /><figcaption>Figure 3: Functional realisations from draws in Figure 2. Approximated $2-\sigma$ band of the stochastic process in dashed red.</figcaption></figure>

### 3.2.2. Quality indicator

To attain a good accuracy for this random field approximation, a key element is the number of inducing points. We propose in this section an indicator to assess the approximation quality. Intuitively, if the density of inducing points is high enough with respect to the lengthscales, the predictive means $q_{draw}^{(i)}$ (Eq. (18)) will be able to represent all scales of the random process. In such cases, the associated predictive variance remains very low over the whole design space $\mathcal{X}$. The GP covariance reads:

$$\Sigma_{tot}(z_{draw}) = K_{**} - k_{*}^{\top} K_{tot}^{-1} k_{*}, \quad (19)$$
with $K_{**} = k_{train}(\mathbf{z}_{draw}, \mathbf{z}_{draw})$. The predictive variance vector $\sigma_{tot}^2(\mathbf{z}_{draw})$ is the diagonal of $\Sigma_{tot}(\mathbf{z}_{draw})$. We propose to base the quality indicator on the mean predictive variance over all points $\mathbf{z}_{draw}$. More precisely, as $K_{tot}^{-1}$ is positive definite, $\mathbf{k}_{*}^{\top} K_{tot}^{-1} \mathbf{k}_{*} > 0$ and thus $\forall i, \sigma_{toti}^2(\mathbf{z}_{draw}) \in [0, K_{**ii}]$, we propose the following normalised indicator:

$$ \alpha = 1 - \frac{1}{N_{draw}} \sum_{i=1}^{N_{draw}} \frac{\sigma_{toti}^2(\mathbf{z}_{draw})}{K_{**ii}}, \qquad (20) $$

where $\alpha \in [0, 1]$ and $\alpha \simeq 1$ denotes a low predictive variance, meaning a good approximation quality.

Note that for most of the usual kernel functions, the diagonal of $K_{**}$ is filled with $\sigma^2$, the variance hyperparameter of the chosen kernel. This is notably the case with the exponential, squared exponential and Matern kernels. In these cases, using $\sigma^2 N_{draw} = \text{Tr}(K_{**})$ and linearity of the trace, the indicator reads:

$$ \alpha = \frac{\text{Tr}(\mathbf{k}_{*}^{\top} K_{tot}^{-1} \mathbf{k}_{*})}{\sigma^2 N_{draw}}. \qquad (21) $$

One may recognise a Frobenius norm at the numerator of the proposed indicator:

$$ \begin{aligned} \text{Tr}(\mathbf{k}_{*}^{\top} K_{tot}^{-1} \mathbf{k}_{*}) &= \text{Tr}(\mathbf{k}_{*}^{\top} (L_{tot} L_{tot}^{\top})^{-1} \mathbf{k}_{*}), \\ &= \text{Tr}((L_{tot}^{-1} \mathbf{k}_{*})^{\top} L_{tot}^{-1} \mathbf{k}_{*}), \\ &= ||L_{tot}^{-1} \mathbf{k}_{*}||_{F}^{2}, \end{aligned} \qquad (22) $$

where $L_{tot}$ refers again to the Cholesky decomposition of $K_{tot}$.

The proposed indicator is pictured in Figure 4 with nine and three inducing points, respectively. The approach gives an empirical $2-\sigma$ band very close to the real one in Figure 4(a) associated with a value of the quality indicator of 0.995. On the contrary, the approximated random field in Figure 4(b) is quite different from the underlying GP, and the quality indicator only has a value of 0.825.

For deeper understanding of this indicator, we study the behaviour of $\alpha$ on an analytical example. With $D$ the dimension of the input parameters $\mathbf{z}$, we study the following quantity of interest:

$$ q(\mathbf{z}) = \sin(\mathbf{L}^{\top}\mathbf{z}) = \sin\left(\sum_{i=1}^{D} L_i z_i\right). \qquad (23) $$

To vary the complexity of $q$, we regard parameters $D$ and $\mathbf{L}$ are random variables, which yields a random function $q_{D,\mathbf{L}}$. Here, we take $D \sim \mathcal{U}(\llbracket 1, 4 \rrbracket)$ and $\mathbf{L} \sim \mathcal{U}(\llbracket 3, 20 \rrbracket)^D$.

For each sampled function $q_{D(k),\mathbf{L}(k)}$, a GP is built on 10 training data $q_{D(k),\mathbf{L}(k)}(\mathbf{z}_{train})$. A set of points $\mathbf{z}_{draw}$ is drawn uniformly within the input space. The predictive covariance $K$ on these points is given by Equation (13). A set of well spread inducing points of size $N_{ind} \sim \mathcal{U}(\llbracket 10, 60 \rrbracket)$ is also drawn (e.g. with a maximin Latin Hypercube Sampling (LHS)).
<figure><img src="image_5.png" /><figcaption>Figure 4: Approximated random field with (a) 9 inducing points, $\alpha = 0.995$ and (b) 3 inducing points, $\alpha = 0.825$.</figcaption></figure>

strategy) and $10^5$ realisations $q_{D^{(k)}, \mathcal{L}^{(k)}}^{(i)}(\mathbf{z}_{draw})$ are computed at $10^3$ locations $\mathbf{z}_{draw}$. The empirical covariance between all these realisations gives an approximation $\tilde{K}$ of $K$ and the indicator $\alpha$ from Equation (20) is computed.

All this process is repeated 300 times, and we plot in Figure 5 the evolution of the normalised Frobenius distance between $K$ and $\tilde{K}$ with respect to $\alpha$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Accuracy of the empirical covariance $\tilde{K}$ approximated with inducing points with respect to the proposed indicator.</figcaption></figure>
Figure 5 shows that a strong correlation exists between the proposed indicator $\alpha$ and the accuracy of the empirical covariance. A value of the indicator very close to 1 seems to ensure an accurate representation of the random field. Notably, the zoomed window permits to see that when $\alpha$ is greater than 0.99 or 0.995, the worst recorded empirical covariances only yield 5 to 10% relative error.

Note that several appropriate norms for covariance matrices exist. Namely, the logarithmic, Cholesky-based, and Riemannian distances between covariance matrices allow to preserve positive definiteness and have many desirable properties [43]. However, because we draw very badly conditioned empirical covariances, very small negative eigenvalues may arise, that prohibit the use of such distances. For its simplicity and interpretability, the Frobenius distance is the most reliable choice.

### 3.2.3. Adaptive strategy

To take advantage of the proposed quality indicator, we propose an iterative algorithm ensuring that $\alpha$ is sufficiently close to 1 before returning the approximated realisations. The number of inducing points is hence increased as long as a prescribed value for $\alpha$ has not been reached. At each iteration, these points are spread in the input space using the classical maximin criterion on a Latin Hypercube Sampling (LHS). The algorithm is given in Algorithm 1.

**Algorithm 1** Algorithm for choosing the number of inducing points
1: Choose threshold $\alpha_* \in [0, 1]$
2: Choose initial number of inducing points $N_{init}$ and added per iteration $N_{add}$
3: Build a GP on the training data
4: $N_{ind} = N_{init}$
5: **while** True **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;Construct a DoE of $N_{ind}$ points using a maximin LHS
7: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\alpha$ from Equation (20)
8: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\alpha \ge \alpha_*$ **then**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw realisations with Equation (18) using the DoE of inducing points
10: &nbsp;&nbsp;&nbsp;&nbsp;Return realisations
11: &nbsp;&nbsp;&nbsp;&nbsp;**else**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$N_{ind} += N_{add}$
13: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
14: **end while**

In the following, the threshold $\alpha_*$ will take values of 0.99 or 0.995, basing on Figure 5. The simple approach proposed here is very intuitive and reveals very cost-efficient.

## 4. Application to Optimisation Under Uncertainty

In this section, we give a general technique, heavily inspired from the Surrogate-Assisted Bounding-Box approach introduced in [44, 31], with the aim of making use of the non-
parametric approximations presented in the previous section in an Optimisation Under Uncertainty (OUU) context.

A general formulation of the proposed approach is presented in Section 4.1, and its application with non-parametric measure approximation is made explicit in Section 4.2.

### 4.1. General formulation

The proposed approach integrates the concept of tunable accuracy within the OUU process. This concept has notably been presented in [17] and exploited in [30, 44, 31], with the aim of tuning the accuracy of each computation of $\rho$ on the fly in order to focus computational power on the most promising designs.

In practice, robustness and reliability measures are computed using some evaluations of $q$, and their approximation error can be estimated. We propose to treat these approximated measures as random vectors, that can be compared under the $\epsilon$-relaxed probabilistic constrained Pareto dominance paradigm, proposed in Section 2.3, in order to discard all dominated approximations

Following the developments from [31], we further propose to rank all remaining measures with the $POP_{min}$ metric, which is an approximation of the Pareto Optimal Probability (POP) with desirable properties. Using the notations from Section 3, the metric writes as follows:

$$ POP_{min}(\mathbf{P}(\mathbf{x}_k)) = \min_{\substack{j \in \llbracket 1, N_x \rrbracket \\ j \neq k}} \left( \mathbb{P}_{\mathbf{P}}[\mathbf{P}(\mathbf{x}_j) \not\leq_c \mathbf{P}(\mathbf{x}_k)] \right). \quad (24) $$

The accuracy is then adaptively tuned by improving these approximations only on the most promising individuals. We refer to this approach as Measure Approximation with Tunable Accuracy (MATA), as different levels of refinement coexist, depending on the estimated performances.

In addition, we construct a surrogate model on these measure approximations in the design space. This model is referred to in the following as the Surrogate-Assisting (SA) model. It directly yields predictions of the robustness and reliability measures at any design $x$, so as to bypass $\rho$ estimations when the SA model is sufficiently accurate. In essence, it permits to further lower the required number of function evaluations when the measures $\rho$ are easier to model in $\mathcal{X}$ than $q$ in the coupled $\mathcal{X} \times \Xi$ space.

These two methods (SA and MATA) couple very efficiently to lower the cost of an uncertainty-based optimisation process. Figure 6 gives a flowchart of the global strategy, called SAMATA. The specific notations are presented in the following sections 4.1.1 and 4.1.2.

#### 4.1.1. Surrogate-Assisting strategy (SA)

Surrogate-Assisting (SA) strategies are common practice for accelerating the optimisation process. In the proposed approach, a surrogate model is constructed and updated throughout the optimisation process directly on the robustness and reliability measures $\rho$ in the design space. For each visited design $x_i$, an aleatoric approximation $\mathbf{P}(x_i)$ is computed. The SA
<figure><img src="image_7.png" /><figcaption>Figure 6: SAMATA strategy flowchart.</figcaption></figure>

model then builds on the training set $\{(x_i, P(x_i))\}_i$. When enough information is gathered at a new design $x_{new}$, the SA model bypasses the computation of $P(x_{new})$ and returns the current prediction $\rho_{SA}(x_{new})$. One can expect this approach to be extensively exploited at the end of the optimisation process, when a lot of training points are available and most new designs are gathered in the optimal area.

The SA model takes as training data aleatoric approximations of $\rho$, and must yield an accuracy metric at any location. We propose to write the accuracy of the SA model at a new design $x_{new}$ as a vector $\Delta_{SA}(x_{new})$. The SA prediction $\rho_{SA}(x_{new})$ is then directly returned to the optimisation process whenever $\Delta_{SA}(x_{new}) \le s_1$, where $s_1$ is a user-defined threshold chosen beforehand.

Similarly to [31], we propose to carry out the whole optimisation process several times with decreasing values of $s_1$. All training data $\{(x_i, P(x_i))\}_i$ are transfered from one iteration to another. This procedure allows to get several intermediate results, that compromise between computational cost and accuracy. With high $s_1$ value, the number of evaluations of $q$ is kept quite low as the SA model prediction $\rho_{SA}$ is quickly exploited. On the contrary, with a small $s_1$ value, many training data are required before employing the SA model, thus ensuring a better accuracy of $\rho_{SA}$.

The structure of this approach is depicted in Algorithm 2. The computation and refinement of the probabilistic approximations $P(x)$ (lines 4 and 12 of the algorithm) are described in the following section 4.1.2.

### 4.1.2. Measure Approximation with Tunable Accuracy (MATA)

Throughout the optimisation process, we denote $\mathcal{X}_c$ the set of designs at which an approximation $P$ has been computed. At each iteration, among the new designs $\mathcal{X}_{new}$ given by the optimisation algorithm, all designs which computation is not bypassed by the SA strategy are added to $\mathcal{X}_c$ and first approximations $P$ and $\Delta$ are computed. Within this set, we denote by $\mathcal{X}_r$ the set of designs that must be further refined. These designs satisfy two criteria: (i) they belong to the current Pareto front $\mathcal{P}(P(\mathcal{X}_c))$ and (ii) the user-defined accuracy threshold $s_2$ is not reached. Formally:

$$ \mathcal{X}_r = \{x \in \mathcal{X}_c | \Delta(x) \not\le s_2 \ \& \ P(x) \in \mathcal{P}(P(\mathcal{X}_c))\}, \quad (25) $$

where $\Delta(x)$ is a chosen measure of the variability of $P(x)$, here its standard deviation, and $\mathcal{P}$ corresponds to the non-dominated measures under the $\epsilon$-relaxed probabilistic constrained
<table>
  <thead>
    <tr>
      <th>Algorithm 2</th>
      <th>Algorithm of the Surrogate-Assisting strategy</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>1:</td>
      <td>Set initial thresholds <b>s</b><sub>1</sub></td>
    </tr>
    <tr>
      <td>2:</td>
      <td>Initialise the training set <i>T</i> = {(<b>x</b><sub>train<sub>i</sub></sub>, <b>P</b>(<b>x</b><sub>train<sub>i</sub></sub>))}<sub>i</sub> empty</td>
    </tr>
    <tr>
      <td>3:</td>
      <td><b>for</b> k from 0 to <i>N</i><sub>threshold</sub> <b>do</b></td>
    </tr>
    <tr>
      <td>4:</td>
      <td>Refine approximations in <i>T</i></td>
    </tr>
    <tr>
      <td>5:</td>
      <td>Launch optimiser</td>
    </tr>
    <tr>
      <td>6:</td>
      <td><b>while</b> optimisation running <b>do</b></td>
    </tr>
    <tr>
      <td>7:</td>
      <td>Update the SA model on <i>T</i></td>
    </tr>
    <tr>
      <td>8:</td>
      <td>Get new designs <i>X</i><sub>new</sub> to visit</td>
    </tr>
    <tr>
      <td>9:</td>
      <td><i>X</i><sub>SA</sub> = {<b>x</b> ∈ <i>X</i><sub>new</sub> | Δ<sub>SA</sub>(<b>x</b>) ≤ <b>s</b><sub>1</sub>}</td>
    </tr>
    <tr>
      <td>10:</td>
      <td><i>X</i><sub>c</sub> = <i>X</i><sub>new</sub> \ <i>X</i><sub>SA</sub></td>
    </tr>
    <tr>
      <td>11:</td>
      <td>∀<b>x</b> ∈ <i>X</i><sub>SA</sub>, return <i>ρ</i><sub>SA</sub>(<b>x</b>) to the optimiser</td>
    </tr>
    <tr>
      <td>12:</td>
      <td>∀<b>x</b> ∈ <i>X</i><sub>c</sub>, compute <b>P</b>(<b>x</b>) and return <i>ρ̃</i>(<b>x</b>) to the optimiser</td>
    </tr>
    <tr>
      <td>13:</td>
      <td>∀<b>x</b> ∈ <i>X</i><sub>c</sub>, add (<b>x</b>, <b>P</b>(<b>x</b>)) to <i>T</i></td>
    </tr>
    <tr>
      <td>14:</td>
      <td><b>end while</b></td>
    </tr>
    <tr>
      <td>15:</td>
      <td>Take next <b>s</b><sub>1</sub> value</td>
    </tr>
    <tr>
      <td>16:</td>
      <td><b>end for</b></td>
    </tr>
  </tbody>
</table>

367 *Pareto dominance* paradigm.

Given this set of designs to refine, we propose to use the POP ranking (Eq. 24) to select the most promising design within $\mathcal{X}_r$. The variability of the associated approximation $\mathbf{P}(\mathbf{x})$ is then reduced, at the cost of some evaluations of $\mathbf{q}$. We denote $\mathbf{x}^*$ the design maximising the $POP_{min}$ metric.

$$ \mathbf{x}^* = \underset{\mathbf{x} \in \mathcal{X}_r}{\operatorname{argmax}} POP_{min}(\mathbf{P}(\mathbf{x})). \qquad (26) $$

368 The MATA refinement technique, for a given set of designs $\mathcal{X}_c$ and a fixed threshold $\mathbf{s}_2$, 369 is detailed in Algorithm 3.

**Algorithm 3** Algorithm of the Measure Approximation with Tunable Accuracy technique

<table>
  <tbody>
    <tr>
      <td>1:</td>
      <td>Read threshold <b>s</b><sub>2</sub> and design set <i>X</i><sub>c</sub></td>
    </tr>
    <tr>
      <td>2:</td>
      <td>Compute first approximations <b>P</b>(<b>x</b>) or retrieve the ones that already exist</td>
    </tr>
    <tr>
      <td>3:</td>
      <td>Extract <i>P</i>(<b>P</b>(<i>X</i><sub>c</sub>))</td>
    </tr>
    <tr>
      <td>4:</td>
      <td>Compute <i>X</i><sub>r</sub> with Equation (25)</td>
    </tr>
    <tr>
      <td>5:</td>
      <td><b>while</b> <i>X</i><sub>r</sub> is non-empty <b>do</b></td>
    </tr>
    <tr>
      <td>6:</td>
      <td>Find <b>x</b><sup>*</sup> with Equation (26)</td>
    </tr>
    <tr>
      <td>7:</td>
      <td>Refine approximation <b>P</b>(<b>x</b><sup>*</sup>) with some evaluations of <b>q</b></td>
    </tr>
    <tr>
      <td>8:</td>
      <td>Update <i>P</i>(<b>P</b>(<i>X</i><sub>c</sub>)) and <i>X</i><sub>r</sub></td>
    </tr>
    <tr>
      <td>9:</td>
      <td><b>end while</b></td>
    </tr>
    <tr>
      <td>10:</td>
      <td><b>return</b> <b>P</b>(<b>x</b>) for all <b>x</b> in <i>X</i><sub>c</sub></td>
    </tr>
  </tbody>
</table>

370 Similarly to the SA strategy, we propose to iterate over decreasing values of $\mathbf{s}_2$ in order 371 to obtain intermediate approximations of the Pareto front.
## 4.2. Application to non-parametric measure approximations

The above SAMATA strategy is presented in its most general form. In practice, specific choices must be made on the computation of approximations $P(x)$, metrics $\tilde{\rho}(x)$ and $\Delta(x)$, and SA predictions $\rho_{SA}(x)$ and $\Delta_{SA}(x)$.

The inducing point strategy presented in Section 3.2.1 permits to generate realisations of $q$ on an extensive set of points $z_{draw}$, to empirically estimate realisations of the measures $P$, e.g. using Equation (16). In short, the obtained empirical distribution of $P$ comes from the Gaussian assumption of the GP surrogate model built on $q$, composed with the specific statistical estimator.

The mean of this approximated distribution is returned to the optimisation process as $\tilde{\rho}$, and three times the standard deviation is used as accuracy metric $\Delta$. Concerning the refinements, we use the criterion proposed in [31], corresponding to a linear combination of measure-specific criteria.

Finally, the SA model is constructed in a specific way to best exploit all the available information. SA-based measure approximations should take the form of set of samples and be jointly drawn with the realisations obtained in (16). To this extent, we propose to define the SA-based realisations as random draws of a GP surrogate model conditioned on the computed joint realisations $P_c^{(i)}$. In practice, on any set of design $x_{SA}$, SA-based measure realisations $P_{SA}^{(i)}(x_{SA})$ are computed as such:

$$P_{SA}^{(i)}(x_{SA}) = \mu_{SA}^{(i)}(x_{SA}) + L_{SA}Y^{(i)}, \quad Y \sim \mathcal{N}(0, I), \qquad (27)$$

where the predictive mean writes:

$$\mu_{SA}^{(i)}(x_{SA}) = k_{*}^{\top} K_c^{-1} P_c^{(i)}, \qquad (28)$$

and $L_{SA}$ results from a Cholesky decomposition performed on the predictive covariance $\Sigma_{GP}(x_{SA})$ defined in Equation (13). In the above, $P_c^{(i)}$ refers to realisations drawn from Equation (16), on the set designs $x_c$. $K_c$ corresponds to the autocovariance matrix $k(x_c, x_c)$ and $k_*$ to the covariance $k(x_c, x_{SA})$. Note that hyperparameters are only optimised once for numerical efficiency.

This SA model construction is illustrated on an example set of samples in Figure 7, where realisations of $P$ are jointly drawn at ten locations in ascending order. Two bimodal sets of samples are notably drawn at $x = 0.15$ and $x = 0.45$ in order to emphasise the need for non-parametric reconstruction. Some SA-based functional realisations $P_{SA}^{(i)}$ form Equation (27) are represented. The two bi-modal densities are naturally well captured using the proposed approach. All computed $P_c^{(i)}$ and SA-based $P_{SA}^{(i)}$ realisations are jointly redrawn at each computed measure refinement, in order to use the joint dominance probability from Equation (12).

These formulations can be readily incorporated within the SAMATA strategy, presented in Section 4.1 and pictured in Figure 6. The specific algorithm is given in Algorithm 4.

Finally, specific numerical choices such as the optimisation algorithm and the threshold values $s_1$ and $s_2$ will be made case by case.
<figure><img src="image_8.png" /><figcaption>Figure 7: Realisations of the SA model. Training data in light blue and underlying GP model in dark blue, in the background.</figcaption></figure>

**Algorithm 4** SAMATA algorithm with joint measure approximations
1: Launch $N_{init}$ evaluations $q(\boldsymbol{x}_i, \boldsymbol{\xi}_i)$
2: Generate a first coupled space surrogate model $\hat{\boldsymbol{q}}$
3: Loop over values of $\boldsymbol{s}_1$ and $\boldsymbol{s}_2$
4: **while** optimisation running **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;Get new designs $\mathcal{X}_{new}$
6: &nbsp;&nbsp;&nbsp;&nbsp;Initialise $\mathcal{X}_c$ empty
7: &nbsp;&nbsp;&nbsp;&nbsp;Compute $P_{SA}$ for all $\boldsymbol{x} \in \mathcal{X}_{new}$ from (27)
8: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\mathcal{X}_{SA}^0 = \{\boldsymbol{x} \in \mathcal{X}_{new} | \Delta_{SA}(\boldsymbol{x}) \le \boldsymbol{s}_1\}$ and $\mathcal{X}_c^0 = \mathcal{X}_{new} \setminus \mathcal{X}_{SA}^0$
9: &nbsp;&nbsp;&nbsp;&nbsp;Deduce $\mathcal{X}_r^0$ from $\mathcal{X}_c^0$ using (25)
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k = 0$
11: &nbsp;&nbsp;&nbsp;&nbsp;**while** $\mathcal{X}_r^k \ne \emptyset$ **do**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $\forall \boldsymbol{x} \in \mathcal{X}_{new}, POP_{min}(\boldsymbol{x})$ w.r.t. all designs using (24) and (12)
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Find $\boldsymbol{x}^*$ from $\mathcal{X}_r^k$ using (26)
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $\boldsymbol{x}^* \notin \mathcal{X}_c$ **then**
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Add $\boldsymbol{x}^*$ in $\mathcal{X}_c$
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Find $\boldsymbol{\xi}^*$ that maximises the refinement criterion from [31]
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Launch evaluation of $q(\boldsymbol{x}^*, \boldsymbol{\xi}^*)$
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update the coupled space surrogate model $\hat{\boldsymbol{q}}$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update all approximations $P_c^{(i)}(\boldsymbol{x})$ for $\boldsymbol{x} \in \mathcal{X}_c$ using empirical estimators like (16)
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update $P_{SA}$ for all $\boldsymbol{x} \in \mathcal{X}_{new} \setminus \mathcal{X}_c$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute $\mathcal{X}_c^{k+1}, \mathcal{X}_{SA}^{k+1}$ and $\mathcal{X}_r^{k+1}$
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$k = k + 1$
24: &nbsp;&nbsp;&nbsp;&nbsp;**end while**
25: &nbsp;&nbsp;&nbsp;&nbsp;Return $\mathbb{E}[P(\boldsymbol{x})]$ to the optimiser for all $\boldsymbol{x} \in \mathcal{X}_{new}$
26: &nbsp;&nbsp;&nbsp;&nbsp;**end while**
27: Return all measure approximations $P(\boldsymbol{x})$ and their associated $POP_{min}$ score
## 5. Analytic comparison

In this section, the proposed SAMATA approach depicted in Alg. 4 is quantitatively compared to SABBa from [31] with Coupled-Space surrogate models, which was shown very cost efficient. As a reference, we also compare these techniques to an A Priori Metamodel strategy, where a Coupled-Space surrogate model is constructed in the coupled $\mathcal{X} \times \Xi$ space on a given number of data points and is readily employed as is for measure estimation and optimisation.

The three approaches are compared in terms of the expected modified Hausdorff distance between the approximated Pareto optimal designs and the true ones. Such a metric was proposed in [31] in the context of Bounding-Box approach. Here, by construction, joint realisations $\mathbf{P}^{(i)}(\mathbf{x})$ of the aleatory measure approximations can be drawn simultaneously on a set of designs $\mathbf{x}$ using (16) and (27). This readily permits to employ the aforementioned metric, which formula is recalled here:

$$Q_B = \mathbb{E}_{\tilde{\mathcal{X}}_P} [d'_H(\mathcal{X}_P, \tilde{\mathcal{X}}_P)], \quad (29)$$

where $d'_H$ is the modified Hausdorff distance from [45], $\mathcal{X}_P$ is the Pareto front pre-image, and $\tilde{\mathcal{X}}_P$ is an aleatory variable corresponding to the pre-image of the currently approximated aleatory Pareto front. This expected value is empirically estimated, exploiting the capability of drawing joint realisations of $\mathbf{P}$. Note that in a normalised input space, a value $Q_B = 1$ corresponds to a bad score while $Q_B = 10^{-2}$ reveals a very accurate approximation of the Pareto front.

The NSGA-II optimisation algorithm is used in the following test-cases. The number of individuals per generation is set to 32. Both the optimisation and refinement processes are sequential, and the coupled-space surrogate $\hat{q}$ is initialised with a Latin Hypercube Sampling (LHS) of 30 points.

Each technique is run ten times to get the mean and standard deviation of the convergence curves.

### 5.1. Test-case 1: Unconstrained Taguchi optimisation

This problem is a bi-objective robust optimisation proposed in [30]. There are two design variables $\mathbf{x}$ and one uncertain parameter $\xi$. The problem reads:

$$
\begin{array}{ll}
\text{minimise:} & \rho_f(\mathbf{x}) = \begin{pmatrix} \mu(\mathbf{x}) \\ \sigma^2(\mathbf{x}) \end{pmatrix} \\
\text{where:} & \mu(\mathbf{x}) = \mathbb{E}_{\xi}[q(\mathbf{x}, \xi)] \\
& \sigma^2(\mathbf{x}) = \mathbb{V}_{\xi}[q(\mathbf{x}, \xi)] \\
\text{with:} & q(\mathbf{x}, \xi) = \xi - x_1\xi^5 + \cos(2\pi x_2\xi) + 5 \\
& \xi \sim \mathcal{U}([0, 1]) \\
\text{by changing:} & (x_1, x_2) \in [1, 2]^2
\end{array}
\quad (30)
$$

The Pareto front associated with this problem is discontinuous, and the optimal set in the design space consists in segment and a point (see Fig. 8).
<figure><img src="image_9.png" /><figcaption>Figure 8: Test-case 1: a) Discretisation of the design space in red and Pareto optimal sets in black. b) Image of the discretised points in the objective ($\mu, \sigma^2$) space.</figcaption></figure>

The performance of the three compared strategies is depicted in Figure 9, showing the mean convergence over ten runs and the Bollinger Band corresponding to the associated standard deviation. Both the Bounding-Box and the proposed non-parametric SAMATA approaches show a steeper convergence than the A Priori Metamodel technique. In addition, the non-parametric approach reveals significantly more parsimonious than the other two techniques, reaching a $Q_B$ metric of $10^{-2}$ around 60 function evaluations, corresponding to a gain of nearly one order of magnitude over the reference methods.

<figure><img src="image_10.png" /><figcaption>Figure 9: Test-case 1: $Q_B$ metric for the A Priori Metamodel (solid blue), SABBa (dashed orange) and the proposed non-parametric SAMATA (dotted green) approaches.</figcaption></figure>
## 5.2. Test-case 2: Quantile-constrained mean optimisation

This second test-case, proposed in [31], assesses the performance of the proposed non-parametric approach on a mono-objective optimisation problem with a reliability based constraint, in a higher dimensional context. The objective is a robustness measure derived from a simplified Rosenbrock function and the constraint a quantile formulation of the Six-Hump Camel function. It features four design variables and three uncertain parameters, and is stated as follows:

$$ \text{minimise: } \rho_f(\boldsymbol{x}) = \mu(\boldsymbol{x}) $$

$$ \text{satisfying: } \rho_g(\boldsymbol{x}) = q^{95\%}(\boldsymbol{x}) \le 1 $$

$$ \text{where: } \mu(\boldsymbol{x}) = \mathbb{E}_{\boldsymbol{\xi}}[q_1(\boldsymbol{x}, \boldsymbol{\xi})] \\ q^{95\%}(\boldsymbol{x}) = q_{\boldsymbol{\xi}}^{95\%}[q_2(\boldsymbol{x}, \boldsymbol{\xi})] $$

$$ \text{with: } q_1(\boldsymbol{x}, \boldsymbol{\xi}) = \sum_{i=1}^{3} \left[ (1 - x_i) + 3 \left( 1 + \frac{\arctan(5(\xi_i - 0.5))}{2} \right) (x_{i+1} - x_i^2)^2 \right] \\ q_2(\boldsymbol{x}, \boldsymbol{\xi}) = \left( 4 - 2.1x_1^2 + \frac{x_1^4}{3} \right) x_1^2 + x_1x_2 + (-4 + 4x_2^2)x_2^2 \\ + \frac{\cos(2\pi\xi_1) - \sin(\frac{\pi}{2}\xi_1) - \xi_1 - \left( \cos(2\pi0.05) - \sin(\frac{\pi}{2}0.05) - 0.05 \right)}{5} \\ \boldsymbol{\xi} \sim \mathcal{U}([0, 1]^3) \\ \text{by changing: } \boldsymbol{x} \in [-0.2, 1.2]^4 \qquad (31) $$

The optimum design is $\boldsymbol{x}_* \approx (0.7033, 0.7035, 0.6212, 0.3859)$ with $q^{95\%}(\boldsymbol{x}_*) = 1$ and $\mu(\boldsymbol{x}_*) \approx 0.4981$.

The convergence curves are depicted in Figure 10. Similarly to Figure 9, the A Priori Metamodel is the less efficient strategy, and the proposed non-parametric SAMATA approach is the most parsimonious. The standard deviation over the ten runs are greater in this test-case, but the proposed method permits to reach a metric of $10^{-1}$ for less than 70 function evaluations, against 110 for the Bounding-Box approach and more than 200 with the A Priori Metamodel.

## 5.3. Benefit of sampling-based approximations

The so-called projected processes techniques [32, 33, 34, 35] permit to solve Optimisation Under Uncertainty (OUU) problems in a Bayesian Optimisation setting (also called Efficient Global Optimisation, EGO). They however rely on independent Gaussian approximations of the robustness and reliability measures. Such Gaussian approximations are cheap to estimate and particularly adapted for expectation-like measures, but might fail to accurately capture more complex or correlated distributions. For deeper insights on this method, see [33].

We run ten optimisations of the first test-case using the proposed approach, but replacing the jointly drawn realisations of $\boldsymbol{P}$ with samples from a Gaussian distribution which mean and variance correspond to the empirical moments from the initial samples. In this manner,
<figure><img src="image_11.png" /><figcaption>Figure 10: Test-case 2: $Q_B$ metric for the A Priori Metamodel (solid blue), SABBa (dashed orange) and the proposed non-parametric SAMATA (dotted green) approaches.</figcaption></figure>

we aim to investigate solely the impact of using a different density estimation, with either
the proposed sample-based strategy or a simple moment-based Gaussian fit.

<figure><img src="image_12.png" /><figcaption>Figure 11: Convergence curves of the non-parametric SAMATA strategy (solid blue) against its Gaussianised equivalent (dashed orange). Variability represented with the translucent $\pm\sigma$ band.</figcaption></figure>

Figure 11 shows the two convergence curves, with their Bollinger bands. Although the convergence seems to plateau around the same metric value, the Gaussianised approach
shows slightly inferior accuracy compared to the original non-parametric strategy along the whole study. Intuitively, we argue that accurate distribution shapes and dependencies are lost when Gaussian density is forced, which may deteriorate the intermediate approximations of the Pareto optimal set.

# 6. ORC application

We perform here a robust shape optimisation of a typical converging-diverging turbine nozzle for ORC applications. It was extensively studied in the context of deterministic [46, 47, 48] and multi-point optimisation [49]. Some recent works tackled such design under epistemic uncertainties [50] due to turbulence modelling, and investigated the robustness of the blade design under aleatoric uncertainties in the operating conditions, thermodynamic model parameters and geometric tolerances [51]. We propose here to minimise the mean of the performance function while the mean mass-flow rate is constrained to lie within a prescribed range centered on the nominal value.

The Biere is a reference two-dimensional benchmark geometry to test the design of devices operating with the siloxane fluid MDM (Octamethyltrisiloxane, $C_8H_{24}O_2Si_3$). This blade profile is meant to obtain a convergent-divergent cascade passage which serves to accelerate the fluid up to supersonic speed. Compressibility effects play a crucial role and a typical fish-tail shock pattern is generated downstream the trailing edge. Strong shocks may cause significant losses, design of the trailing edge region is thus critical for the turbine efficiency. The Mach field is depicted in Figure 12 with the baseline blade profile.

<figure><img src="image_13.png" /><figcaption>Figure 12: Mach contours at nominal conditions for the baseline profile and computational grid of 36k cells.</figcaption></figure>

## 6.1. Problem setting

The objective function $\Delta P$ is defined as the standard deviation of the azimuthal distribution of static pressure half an axial chord downstream of the blade Trailing Edge (TE).
Minimising $\Delta P$ should yield a significant reduction of the shock strength. As mentioned above, a constraint is imposed on the mass-flow rate per unit span $\dot{m}$, normalised with respect to the nominal value.

Uncertainties on the operating conditions are considered both at the inlet and outlet of the turbine, through changes of boundary conditions. Following [52], we model the operational variability as independent and uniform uncertain variables, gathered in $\xi = [P_{in}^t, T_{in}^t, P_{out}^s]$. The associated ranges are reported in Table 1, alongside Nominal and Mean conditions. Uncertainties on the parameters of the thermodynamical model are neglected since previous studies [53] provided evidence of their limited impact.

<table><caption>Table 1: Operating Conditions: Nominal, Mean and Random.</caption><thead><tr><th>ξ</th><th>P<sub>in</sub><sup>t</sup> [Pa]</th><th>T<sub>in</sub><sup>t</sup> [K]</th><th>P<sub>out</sub><sup>s</sup> [Pa]</th></tr></thead><tbody><tr><td>Nominal (ξ<sub>0</sub>)</td><td>8 × 10<sup>5</sup></td><td>545.15</td><td>1.072 × 10<sup>5</sup></td></tr><tr><td>Mean (ξ<sub>μ</sub>)</td><td>8 × 10<sup>5</sup></td><td>545.15</td><td>1.5 × 10<sup>5</sup></td></tr><tr><td>Random</td><td>U[7.95 × 10<sup>5</sup>, 8.05 × 10<sup>5</sup>]</td><td>U[541.15, 549.15]</td><td>U[1 × 10<sup>5</sup>, 2 × 10<sup>5</sup>]</td></tr></tbody></table>

The optimisation problem follows a multi-objective Taguchi formulation:

$$
\begin{array}{ll}
\text{minimise:} & \rho_f(\boldsymbol{x}) = \begin{pmatrix} \mu_1(\boldsymbol{x}) \\ \sigma^2(\boldsymbol{x}) \end{pmatrix} \\
\text{satisfying:} & \rho_g(\boldsymbol{x}) = \mu_2(\boldsymbol{x}) \in [0.98\dot{m}_0, 1.02\dot{m}_0] \\
\text{where:} & \mu_1(\boldsymbol{x}) = \mathbb{E}_{\boldsymbol{\xi}}[\Delta P(\boldsymbol{x}, \boldsymbol{\xi})] \\
& \sigma^2(\boldsymbol{x}) = \mathbb{V}_{\boldsymbol{\xi}}[\Delta P(\boldsymbol{x}, \boldsymbol{\xi})] \\
& \mu_2(\boldsymbol{x}) = \mathbb{E}_{\boldsymbol{\xi}}[\dot{m}(\boldsymbol{x}, \boldsymbol{\xi})] \\
& \dot{m}_0 = 344.843 \text{ kg.s}^{-1} \\
\text{with:} & \boldsymbol{\xi} \sim \mathcal{U}(\Xi) \\
\text{by changing:} & \boldsymbol{x} \in \mathcal{X}
\end{array}
\qquad (32)
$$

where $\Xi = [7.95, 8.05] \times [541.15, 549.15] \times [1, 2]$. The design space $\mathcal{X}$ allows to modify the shape of the blade as much as possible while ensuring convergence of the numerical simulation. In practice, it consists of 9 of the 30 Control Points (CP) of a B-spline curve of degree 3, that can be displaced in the direction normal to the baseline geometry. Such parametrisation was proposed in [54] and is illustrated in Figure 13.

A mono-objective version of (32) is also treated for comparison purposes. Results will also be compared with two deterministic optima obtained at the Nominal and Mean operating conditions from Table 1.

Simulation is performed using the Non-Ideal Compressible-Fluid Dynamics solver included in the SU2 [55, 56, 57, 58] suite on unstructured grids, generated using an in-house meshing tool. A dedicated mesh deformation tool based on Radial Basis Functions (RBF) allows high flexibility and robustness while maintaining the grid connectivity throughout the optimisation steps.
<figure><img src="image_14.png" /><figcaption>Figure 13: B-splines parametrisation. Fixed CP in black, moving CP in red.</figcaption></figure>

## 6.2. Mono-objective results

Because nominal conditions are not centered within the uncertain range, we conduct the reference deterministic optimisations both at Nominal and Mean conditions. Similarly to [59], a state-of-the-art Bayesian Optimisation strategy is employed, with an Expected Improvement (EI) criterion multiplied by the Probability of Feasibility (PF) [60]. Convergence is reached after roughly **100** function evaluations.

The robust optimisation problem is solved using the proposed non-parametric SAMATA strategy. The NSGA-II optimisation procedure is again exploited, and constraints are handled through penalisation. The initial DoE consists of 30 evaluations drawn with a maximin Latin Hypercube Sampling (LHS) in the coupled space. Each population contains 32 individuals, and measure refinements are performed sequentially. We iteratively take threshold values of 10%, 5%, 3%, 2% and 1%.

Convergence values are gathered in Table 2, with a final optimum centered on $\rho_f = 8087$, with a standard deviation of the approximation error of 44.26. This optimum satisfies the 1% accuracy imposed by the normalised thresholds and only required **145** simulation calls, which is not significantly higher than the deterministic optimisation. The proposed approach permits to reduce by around one order of magnitude the number of required function evaluations with respect to [59] that employed a more classical nested approach.

<table><caption>Table 2: Optimal outputs associated to the different threshold values.</caption><thead><tr><th rowspan="2">N<sub>eval</sub></th><th colspan="2">ρ<sub>f</sub> [Pa]</th><th colspan="2">ρ<sub>g</sub> [kg.s<sup>-1</sup>]</th></tr><tr><th>Mean</th><th>Std. dev.</th><th>Mean</th><th>Std. dev.</th></tr></thead><tbody><tr><td>76</td><td>8.328 × 10<sup>3</sup></td><td>6.355 × 10<sup>2</sup></td><td>3.464 × 10<sup>2</sup></td><td>1.065 × 10<sup>-1</sup></td></tr><tr><td>91</td><td>8.854 × 10<sup>3</sup></td><td>3.986 × 10<sup>2</sup></td><td>3.514 × 10<sup>2</sup></td><td>3.149 × 10<sup>-2</sup></td></tr><tr><td>106</td><td>8.869 × 10<sup>3</sup></td><td>1.124 × 10<sup>2</sup></td><td>3.387 × 10<sup>2</sup></td><td>1.534 × 10<sup>-2</sup></td></tr><tr><td>118</td><td>8.288 × 10<sup>3</sup></td><td>1.068 × 10<sup>2</sup></td><td>3.500 × 10<sup>2</sup></td><td>1.144 × 10<sup>-2</sup></td></tr><tr><td>145</td><td>8.087 × 10<sup>3</sup></td><td>4.426 × 10<sup>1</sup></td><td>3.505 × 10<sup>2</sup></td><td>1.812 × 10<sup>-2</sup></td></tr></tbody></table>

The point-wise ($\xi_0$ and $\xi_\mu$) and mean values of the objective function are computed at the Nominal $x_0^*$, Mean $x_\mu^*$ and Robust $x_{OUU}^*$ operating conditions. Results are gathered in
Table 3 and the full distributions of $\Delta P$ over $\xi$ at these given designs are plotted in Figure 14. These distributions are obtained with an extensive Monte Carlo on a GP built over $\Xi$ at a given $\boldsymbol{x}$, with 50 training points.

<table><caption>Table 3: Nominal values and statistics for different blade profiles.</caption><thead><tr><th>Blade profile</th><th>ΔP(<b>x</b>, <i>ξ</i><sub>0</sub>) [Pa]</th><th>ΔP(<b>x</b>, <i>ξ</i><sub>μ</sub>) [Pa]</th><th>E<sub>ξ</sub>[ΔP(<b>x</b>, <i>ξ</i>)] [Pa]</th></tr></thead><tbody><tr><td>Nominal optimum <b>x</b> = <b>x</b><sub>0</sub><sup>*</sup></td><td>1.958 × 10<sup>3</sup></td><td>1.526 × 10<sup>4</sup></td><td>1.624 × 10<sup>4</sup></td></tr><tr><td>Mean optimum <b>x</b> = <b>x</b><sub>μ</sub><sup>*</sup></td><td>1.079 × 10<sup>4</sup></td><td>1.773 × 10<sup>3</sup></td><td>8.472 × 10<sup>3</sup></td></tr><tr><td>Robust optimum <b>x</b> = <b>x</b><sub>OUU</sub><sup>*</sup></td><td>1.210 × 10<sup>4</sup></td><td>3.082 × 10<sup>3</sup></td><td>8.084 × 10<sup>3</sup></td></tr></tbody></table>

As expected, the Nominal optimum yields very high variability in Fig. 14(a) but very low nominal output. The Mean and Robust optima are quite similar, Figs 14(b) and 14(c). However, the Robust optimal blade manages to reach a lower value for the expectation of $\Delta P$, with a slightly modified distribution shape.

<figure><img src="image_15.png" /><figcaption>Figure 14: Distribution of $\Delta P$ at (a) $x_0^*$, (b) $x_\mu^*$ and (c) $x_{OUU}^*$. $\Delta P(\boldsymbol{x}, \boldsymbol{\xi}_0)$ represented as a dashed black line, $\Delta P(\boldsymbol{x}, \boldsymbol{\xi}_\mu)$ as a dotted black line and $E_\xi[\Delta P(\boldsymbol{x}, \boldsymbol{\xi})]$ as a plain red line.</figcaption></figure>
### 6.3. Bi-objective results

We then tackle the complete Taguchi formulation (32) of the Optimisation Under Uncertainty of the ORC turbine blade. This formulation aims to find the Pareto front between good expected performance and low variability. The proposed non-parametric SAMATA algorithm is applied, with again the NSGA-II optimisation method and normalised thresholds of 10%, 5%, 2% and 1%.

The Pareto fronts associated with these threshold values are depicted in Figure 15, and reveal great improvements in the Pareto front accuracy throughout the iterations. The total computational cost for attaining the chosen thresholds is of **72**, **122**, **276** and **419** simulation calls respectively. Each cloud of samples associated with a given design has its transparency driven by its POP value (1: fully visible, 0: transparent). The mono-objective optimum found in the preceding section is represented as a red dot and lies on the Pareto front, as expected.

<figure><img src="image_16.png" /><figcaption>Figure 15: Pareto front of the bi-objective OUU problem with normalised threshold values of (a) 10% (72 evaluations), (b) 5% (122 evaluations), (c) 2% (276 evaluations) and (d) 1% (419 evaluations). Transparency driven by POP value and mono-objective optimum from previous section as a red dot.</figcaption></figure>
# 7. Conclusions

This paper describes an efficient strategy for solving multi-objective and constrained Optimisation Under Uncertainty problems. It consists in coupling two primary techniques: i) a Measure Approximation with Tunable Accuracy (MATA) approach that only refines promising individuals up to a given threshold, ii) a Surrogate-Assisting (SA) strategy that constructs a surrogate model of the objective and constraint functions in the design space to bypass further evaluations when it is sufficiently accurate.

To apply this framework in practice, we propose non-parametric approximations of the error distributions using a sampling-based technique. The robustness and reliability measures for a given design are represented as a set of realisations that may exhibit complex shapes and dependencies. These realisations usually reveal more representative and can be jointly drawn between several designs. Issues associated with drawing such joint realisations are highlighted, and an adaptive algorithm is proposed to keep the associated cost manageable. This strategy shows similarities with Nyström approximation for Karhunen-Loève Expansion (KLE) but allows for intuitive control of the overall accuracy and adaptive sampling. A Surrogate-Assisting (SA) model that takes advantage of these joint realisations is introduced to apply the SAMATA strategy. Quantitative comparisons with state-of-the-art approaches on several algebraic test cases reveal a drastic increase in the overall parsimony.

Finally, the proposed approach is applied to an engineering problem dealing with the geometric optimisation of the blade profile in an Organic Rankine Cycle turbine. Convergence is attained for a cost of 100 to 150 evaluations in a mono-objective setting, which is one order of magnitude less than the state-of-the-art approach. The Pareto front of a bi-objective formulation is accurately recovered at the cost of 100 to 400 evaluations.

Several perspectives can be drawn. Regarding overall parsimony, cost improvement could be achieved by improving the Uncertainty Propagation (UP) techniques associated with specific statistical measures. For example, extreme quantile estimation and refinement should be performed with appropriate techniques such as [61] or [62]. Such approaches can be readily incorporated into the proposed framework. Similarly, we only performed sequential refinement in this work, while techniques such as Kriging believer can be exploited to return a set of new locations. Multi-point refinement is critical to take advantage of nowadays' high parallelisation capacities of computational clusters. Finally, the proposed non-parametric error estimation relies on coupled-space surrogate models that may become unfeasible in high dimensional uncertain space or when the simulator is intrinsically stochastic. Specific developments for using SAMATA in this context would effectively widen its range of application.

# Appendix A. Link between inducing points and Nyström method for KLE

In practice, the proposed approach is equivalent to constructing a Karhunen-Loève Expansion (KLE) of a Gaussian random field, when computing the basis functions with Nyström approximation. Computational details are given hereafter to show the equivalence.
To this extent, note that we can write the eigendecomposition (also called spectral decomposition) of the autocovariance matrix $K_{tot}$ from Equation (18), which is Symmetric Positive Definite (SPD) by definition:

$$K_{tot} = U_{tot} \Lambda_{tot} U_{tot}^{\top}, \qquad (A.1)$$

where $\Lambda_{tot}$ is the diagonal matrix containing the eigenvalues associated with the eigenvectors stored in the columns of $U_{tot}$. To make this decomposition unique, the eigenvalues are sorted in decreasing order, $\Lambda_{tot11} \ge \Lambda_{tot22} \ge \dots \ge 0$, and the eigenvectors are orthonormal, so that $U_{tot}^{\top} = U_{tot}^{-1}$, which gives:

$$U_{tot,i}^{\top} U_{tot,j} = \sum_{k} U_{tot_{ki}} U_{tot_{kj}} = \delta_{ij}, \qquad (A.2)$$

with $\delta_{ij}$ the Kronecker delta.

The inverse of $K_{tot}$ can thus be written:

$$K_{tot}^{-1} = U_{tot} \Lambda_{tot}^{-1} U_{tot}^{\top}, \qquad (A.3)$$

and the square root $K_{tot}^{-1/2}$ of this matrix is:

$$K_{tot}^{-1/2} = U_{tot} \Lambda_{tot}^{-1/2}, \qquad (A.4)$$

579 which verifies $K_{tot}^{-1} = K_{tot}^{-1/2} K_{tot}^{-1/2\top}$.

KLE and Nyström approximation. A realisation of a random field described with a truncated Karhunen-Loève Expansion reads as follows:

$$\mathbf{q}_{KLE}^{(i)}(\mathbf{z}) = \sum_{j=1}^{N_{KLE}} \sqrt{\lambda_j} \phi_j(\mathbf{z}) \theta_j^{(i)}, \qquad (A.5)$$

where the eigenvalues $\lambda_j$ and eigenvectors $\phi_j$ are solution of the following homogeneous Fredholm integral equation of the second kind:

$$\int k(\mathbf{z}, \mathbf{z}') \phi_j(\mathbf{z}) d\mathbf{z} = \lambda_j \phi_j(\mathbf{z}') \qquad (A.6)$$

$$\text{with} \quad \int \phi_i(\mathbf{z}) \phi_j(\mathbf{z}) d\mathbf{z} = \delta_{ij}.$$

As presented in [63], the Nyström approximation consists in computing an empirical estimation of the above integral. This Monte Carlo approximation based on points $\mathbf{x}_k$ gives:

$$\frac{1}{N_{MC}} \sum_{k=1}^{N_{MC}} k(\mathbf{z}_k, \mathbf{z}') \phi_j(\mathbf{z}_k) d\mathbf{z} \approx \lambda_j \phi_j(\mathbf{z}') \qquad (A.7)$$

$$\text{with} \quad \frac{1}{N_{MC}} \sum_{k=1}^{N_{MC}} \phi_i(\mathbf{z}_k) \phi_j(\mathbf{z}_k) d\mathbf{z} \approx \delta_{ij}.$$
Replacing $z'$ with all $z_{tot_k}$, and setting $N_{MC} = N_{tot}$ so that the Monte Carlo estimation is performed on locations $z_{tot}$, Equation (A.7) yields:

$$ \frac{1}{N_{tot}} K_{tot} \Phi \approx \Lambda \Phi \qquad (A.8) $$
$$ \text{with } \frac{1}{N_{tot}} \Phi_{\cdot i}^{\top} \Phi_{\cdot j} \approx \delta_{ij}. $$

Matching (A.8) against (A.1) and (A.2) permits to derive:

$$ \Phi \approx \sqrt{N_{tot}} U_{tot}, $$
$$ \Lambda \approx \frac{\Lambda_{tot}}{N_{tot}}, $$

which gives the Nyström approximation of the j-th eigenfunction from Eq. (A.7):

$$ \phi_j(\mathbf{z}) \approx \frac{\sqrt{N_{tot}}}{\Lambda_{totjj}} \sum_{k=1}^{N_{tot}} k(\mathbf{z}, \mathbf{z}_{tot_k}) U_{tot_{kj}}. \qquad (A.9) $$

Using these approximations, KLE-based realisations write:

$$ \begin{aligned} \mathbf{q}_{KLE}^{(i)}(\mathbf{z}) &= \sum_{j=1}^{N_{tot}} \left( \sqrt{\frac{\Lambda_{totjj}}{N_{tot}}} \frac{\sqrt{N_{tot}}}{\Lambda_{totjj}} \sum_{k=1}^{N_{tot}} (k(\mathbf{z}, \mathbf{z}_{tot_k}) U_{tot_{kj}}) \theta_{tot_j}^{(i)} \right) \\ &= \sum_{j=1}^{N_{tot}} \left( \frac{\theta_{tot_j}^{(i)}}{\sqrt{\Lambda_{totjj}}} \sum_{k=1}^{N_{tot}} (k(\mathbf{z}, \mathbf{z}_{tot_k}) U_{tot_{kj}}) \right). \end{aligned} \qquad (A.10) $$

By definition, in this context, $\mathbf{q}_{tot}^{(i)} = K_{tot}^{1/2} \theta_{tot}^{(i)}$, with $\theta_{tot}^{(i)} \sim \mathcal{N}(0, I)$. A realisation at location $\mathbf{z}$ then writes:

$$ \begin{aligned} \mathbf{q}^{(i)}(\mathbf{z}) &= k(\mathbf{z}, \mathbf{z}_{tot}) K_{tot}^{-1} \mathbf{q}_{tot}^{(i)} \\ &= k(\mathbf{z}, \mathbf{z}_{tot}) K_{tot}^{-1/2} \mathbf{\theta}_{tot}^{(i)} \\ &= k(\mathbf{z}, \mathbf{z}_{tot}) U_{tot} \Lambda_{tot}^{-1/2} \mathbf{\theta}_{tot}^{(i)} \\ &= \sum_{j=1}^{N_{tot}} \left( \frac{\mathbf{\theta}_{tot_j}^{(i)}}{\sqrt{\Lambda_{totjj}}} \sum_{k=1}^{N_{tot}} (k(\mathbf{z}, \mathbf{z}_{tot_k}) U_{tot_{kj}}) \right), \end{aligned} \qquad (A.11) $$

<footnote n="580">580 *Inducing points.* We proposed to draw realisations using Equation (18). To simplify the following, we assume here that $N_{train} = 0$ and $\mathbf{q}_{tot} = \mathbf{q}_{ind}$. This implies that matrix $L_{ind}$ corresponds to the square root matrix $K_{tot}^{1/2}$.</footnote>

<footnote>581</footnote>

<footnote>582</footnote>

<footnote n="583">583 which corresponds exactly to Equation (A.10).</footnote>

<footnote>584</footnote>

<footnote>585</footnote>

<footnote>586</footnote>
[1] V. Keshavarzzadeh, R. G. Ghanem, D. A. Tortorelli, Shape optimization under uncertainty for rotor blades of horizontal axis wind turbines, Computer Methods in Applied Mechanics and Engineering 354 (2019) 271–306.

[2] P. Seshadri, P. Constantine, G. Iaccarino, G. Parks, A density-matching approach for optimization under uncertainty, Computer Methods in Applied Mechanics and Engineering 305 (2016) 562–578.

[3] J. Zhang, M. Xiao, P. Li, L. Gao, Quantile-based topology optimization under uncertainty using kriging metamodel, Computer Methods in Applied Mechanics and Engineering 393 (2022) 114690.

[4] V. Keshavarzzadeh, R. M. Kirby, A. Narayan, Stress-based topology optimization under uncertainty via simulation-based gaussian process, Computer Methods in Applied Mechanics and Engineering 365 (2020) 112992.

[5] X. Liu, T. Li, Z. Zhou, L. Hu, An efficient multi-objective reliability-based design optimization method for structure based on probability and interval hybrid model, Computer Methods in Applied Mechanics and Engineering 392 (2022) 114682.

[6] C. Comis Da Ronco, R. Ponza, E. Benini, Aerodynamic shape optimization of aircraft components using an advanced multi-objective evolutionary approach, Computer Methods in Applied Mechanics and Engineering 285 (2015) 255–290.

[7] E. Zitzler, D. Kalyanmoy, L. Thiele, Comparison of multiobjective evolutionary algorithms: Empirical results, Evolutionary computation 8 (2) (2000) 173–195.

[8] H. Ishibuchi, N. Tsukamoto, Y. Nojima, Evolutionary many-objective optimization: A short review, in: 2008 IEEE Congress on Evolutionary Computation (IEEE World Congress on Computational Intelligence), IEEE, 2008, pp. 2419–2426.

[9] X.-S. Yang, Engineering optimization: an introduction with metaheuristic applications, John Wiley & Sons, 2010.

[10] Y. Jin, J. Branke, et al., Evolutionary optimization in uncertain environments-a survey, IEEE Transactions on evolutionary computation 9 (3) (2005) 303–317.

[11] C.-K. Goh, K. C. Tan, Evolutionary multi-objective optimization in uncertain environments, Issues and Algorithms, Studies in Computational Intelligence 186 (2009) 5–18.

[12] K. C. Tan, K. C. Goh, Handling uncertainties in evolutionary multi-objective optimization, in: IEEE World Congress on Computational Intelligence, Springer, 2008, pp. 262–292.
[13] C. P. Ho, P. Parpas, Multilevel optimization methods: Convergence and problem structure (2016).

[14] M. El-Beltagy, A. Keane, Optimisation for multilevel problems: A comparison of various algorithms, in: Adaptive computing in design and manufacture, Springer, 1998, pp. 111–120.

[15] A. March, K. Willcox, Provably convergent multifidelity optimization algorithm not requiring high-fidelity derivatives, AIAA journal 50 (5) (2012) 1079–1089.

[16] M. K. Zahir, Z. Gao, Variable-fidelity optimization with design space reduction, Chinese Journal of Aeronautics 26 (4) (2013) 841–849.

[17] V. Picheny, D. Ginsbourger, Y. Richet, Noisy expected improvement and on-line computation time allocation for the optimization of simulators with tunable fidelity (2010).

[18] J. Teich, Pareto-front exploration with uncertain objectives, in: International Conference on Evolutionary Multi-Criterion Optimization, Springer, 2001, pp. 314–328.

[19] E. J. Hughes, Evolutionary multi-objective ranking with uncertainty and noise, in: International Conference on Evolutionary Multi-Criterion Optimization, Springer, 2001, pp. 329–343.

[20] P. Limbourg, D. E. S. Aponte, An optimization algorithm for imprecise multi-objective problem functions, in: 2005 IEEE Congress on Evolutionary Computation, Vol. 1, IEEE, 2005, pp. 459–466.

[21] H. Eskandari, C. D. Geiger, R. Bird, Handling uncertainty in evolutionary multiobjective optimization: Spga, in: 2007 IEEE Congress on Evolutionary Computation, IEEE, 2007, pp. 4130–4137.

[22] D.-w. Gong, N.-n. Qin, X.-y. Sun, Evolutionary algorithms for multi-objective optimization problems with interval parameters, in: 2010 IEEE Fifth International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA), IEEE, 2010, pp. 411–420.

[23] G. L. Soares, F. G. Guimarães, C. A. Maia, J. A. Vasconcelos, L. Jaulin, Interval robust multi-objective evolutionary algorithm, in: 2009 IEEE Congress on Evolutionary Computation, IEEE, 2009, pp. 1637–1643.

[24] M. Li, R. Silva, F. Guimarães, D. Lowther, A new robust dominance criterion for multiobjective optimization, IEEE Transactions on Magnetics 51 (3) (2015) 1–4.

[25] F. Khosravi, M. Müller, M. Glaß, J. Teich, Uncertainty-aware reliability analysis and optimization, in: 2015 Design, Automation & Test in Europe Conference & Exhibition (DATE), IEEE, 2015, pp. 97–102.
[26] F. Khosravi, M. Borst, J. Teich, Probabilistic dominance in robust multi-objective optimization, in: 2018 IEEE Congress on Evolutionary Computation (CEC), IEEE, 2018, pp. 1-6.

[27] F. Khosravi, A. Raß, J. Teich, Efficient computation of probabilistic dominance in robust multi-objective optimization, arXiv preprint arXiv:1910.08413 (2019).

[28] F. Khosravi, A. Rass, J. Teich, Efficient computation of probabilistic dominance in multi-objective optimization, ACM Trans. Evol. Learn. Optim. 1 (4) (oct 2021). doi:10.1145/3469801.
URL https://doi.org/10.1145/3469801

[29] M. Mlakar, T. Tušar, B. Filipič, Comparing solutions under uncertainty in multiobjective optimization, Mathematical Problems in Engineering 2014 (2014).

[30] F. Fusi, P. M. Congedo, An adaptive strategy on the error of the objective functions for uncertainty-based derivative-free optimization, Journal of Computational Physics 309 (2016) 241-266.

[31] M. Rivier, P. Congedo, Surrogate-assisted bounding-box approach applied to constrained multi-objective optimisation under uncertainty, Reliability Engineering System Safety 217 (2022) 108039. doi:https://doi.org/10.1016/j.ress.2021.108039.
URL https://www.sciencedirect.com/science/article/pii/S0951832021005445

[32] S. Da Veiga, F. Delbos, Robust optimization for expensive simulators with surrogate models: Application to well placement for oil recovery, Safety, Reliability, Risk and Life-Cycle Performance of Structures and Infrastructures 11 (2013) 3321-3328.

[33] V. Baudoui, Optimisation robuste multiobjectifs par modèles de substitution, Ph.D. thesis, Toulouse, ISAE (2012).

[34] J. Janusevskis, R. Le Riche, Simultaneous kriging-based estimation and optimization of mean response, Journal of Global Optimization 55 (2) (2013) 313-336.

[35] B. J. Williams, T. J. Santner, W. I. Notz, Sequential design of computer experiments to minimize integrated response functions, Statistica Sinica (2000) 1133-1152.

[36] R. F. Coelho, Probabilistic dominance in multiobjective reliability-based optimization: Theory and implementation, IEEE Transactions on evolutionary computation 19 (2) (2014) 214-224.

[37] GPy, GPy: A gaussian process framework in python, http://github.com/SheffieldML/GPy (since 2012).

[38] C. E. Rasmussen, C. K. I. Williams, Gaussian processes for machine learning., Adaptive computation and machine learning, MIT Press, 2006.
[39] C. E. Powell, et al., Generating realisations of stationary gaussian random fields by circulant embedding, matrix 2 (2) (2014) 1.

[40] M. Pourahmadi, Covariance estimation: The glm and regularization perspectives, Statistical Science (2011) 369–387.

[41] C. Lantuéjoul, N. Desassis, Simulation of a gaussian random vector: a propagative version of the gibbs sampler, in: The 9th international geostatistics congress, 2012, pp. 174–181.

[42] A. M. Panunzio, R. Cottereau, G. Puel, Large scale random fields generation using localized karhunen–loève expansion, Advanced Modeling and Simulation in Engineering Sciences 5 (1) (2018) 20.

[43] M. Herdin, N. Czink, H. Ozcelik, E. Bonek, Correlation matrix distance, a meaningful measure for evaluation of non-stationary mimo channels, in: 2005 IEEE 61st Vehicular Technology Conference, Vol. 1, IEEE, 2005, pp. 136–140.

[44] M. Rivier, P. M. Congedo, Surrogate-assisted bounding-box approach for optimization problems with tunable objectives fidelity, Journal of Global Optimization 75 (4) (2019) 1079–1109.

[45] M. P. Dubuisson, A. K. Jain, A modified hausdorff distance for object matching, in: Proceedings of 12th International Conference on Pattern Recognition, Vol. 1, 1994, pp. 566–568 vol.1. doi:10.1109/ICPR.1994.576361.

[46] M. Pini, G. Persico, D. Pasquale, S. Rebay, Adjoint method for shape optimization in real-gas flow applications, ASME Journal of Engineering for Gas Turbines and Power 137 (3) (2015).

[47] S. Vitale, T. A. Albring, M. Pini, N. R. Gauger, P. Colonna, Fully turbulent discrete adjoint solver for non-ideal compressible flow applications, Journal of the Global Power and Propulsion Society 1 (2017) Z1FVOI.

[48] G. Persico, P. Rodriguez-Fernandez, A. Romei, High-fidelity shape-optimization of non-conventional turbomachinery by surrogate evolutionary strategies, Journal of Turbomachinery 141 (8) (2019) 081010.

[49] M. Pini, G. Persico, V. Dossena, Robust adjoint-based shape optimization of supersonic turbomachinery cascades, in: ASME Turbo Expo 2014: Turbine Technical Conference and Exposition, American Society of Mechanical Engineers, 2014, pp. V02BT39A043–V02BT39A043.

[50] N. Razaaly, G. Gori, G. Iaccarino, P. Congedo, Optimization of an orc supersonic nozzle under epistemic uncertainties due to turbulence models, in: GPPS 2019, 2019.
722

723

724

[51] N. Razaaly, G. Persico, P. M. Congedo, Impact of geometric, operational, and model uncertainties on the non-ideal flow through a supersonicorc turbine cascade, Energy 169 (2019) 213-227.

725

726

727

728

[52] M. Pini, G. Persico, V. Dossena, Robust adjoint-based shape optimization of supersonic turbomachinery cascades, in: ASME Turbo Expo 2014: Turbine Technical Conference and Exposition, American Society of Mechanical Engineers, 2014, pp. V02BT39A043-V02BT39A043.

729

730

731

[53] P. Congedo, G. Geraci, R. Abgrall, V. Pediroda, L. Parussini, Tsi metamodels-based multi-objective robust optimization, Engineering Computations (Swansea, Wales) 30 (8) (2013) 1032-1053.

732

733

[54] G. Farin, Curves and Surfaces for CAGD: A Practical Guide, 5th Edition, Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2002.

734

735

736

737

738

[55] F. Palacios, M. F. Colonna, A. C. Aranake, A. Campos, S. R. Copeland, T. D. Economon, A. K. Lonkar, T. W. Lukaczyk, T. W. R. Taylor, J. J. Alonso, Stanford University Unstructured (SU2): An open-source integrated computational environment for multi-physics simulation and design, in: 51st AIAA Aerospace Sciences Meeting and Exhibit, 2013.

739

740

741

742

743

[56] T. D. Economon, D. Mudigere, G. Bansal, A. Heinecke, F. Palacios, J. Park, M. Smelyanskiy, J. J. Alonso, P. Dubey, Performance optimizations for scalable implicit {RANS} calculations with {SU2}, Computers & Fluids 129 (2016) 146-158.
doi:http://dx.doi.org/10.1016/j.compfluid.2016.02.003.
URL http://www.sciencedirect.com/science/article/pii/S0045793016300214

744

745

746

747

[57] M. Pini, S. Vitale, P. Colonna, G. Gori, A. Guardone, T. Economon, J. Alonso, F. Palacios, Su2: the open-source software for non-ideal compressible flows, in: NICFD 2016: 1st International Seminar on Non-Ideal Compressible-Fluid Dynamics for Propulsion & Power, Varenna, Italy, 2016.

748

749

750

751

[58] S. Vitale, G. Gori, M. Pini, A. Guardone, T. D. Economon, F. Palacios, J. J. Alonso, P. Colonna, Extension of the SU2 open source CFD code to the simulation of turbulent flows of fluids modelled with complex thermophysical laws, no. AIAA Paper 2015-2760, 2015.

752

753

[59] N. Razaaly, Rare event estimation and robust optimization methods with applications to orc turbine cascade, Ph.D. thesis, Université Paris-Saclay (2019).

754

755

[60] J. R. Gardner, M. J. Kusner, Z. E. Xu, K. Q. Weinberger, J. P. Cunningham, Bayesian optimization with inequality constraints., in: ICML, 2014, pp. 937-945.

756

757

[61] B. Echard, N. Gayton, M. Lemaire, Ak-mcs: an active learning reliability method combining kriging and monte carlo simulation, Structural Safety 33 (2) (2011) 145-154.
758

759

760

[62] N. Razaaly, P. M. Congedo, Novel algorithm using active metamodel learning and importance sampling: application to multiple failure regions of low probability, Journal of Computational Physics 368 (2018) 92–114.

761

762

[63] C. K. Williams, M. Seeger, Using the nyström method to speed up kernel machines, in: Advances in neural information processing systems, 2001, pp. 682–688.

763

764

765

[64] R. Tipireddy, D. A. Barajas-Solano, A. M. Tartakovsky, Conditional karhunen-lo\eve expansion for uncertainty quantification and active learning in partial differential equation models, arXiv preprint arXiv:1904.08069 (2019).