ERGO-II: An Improved Bayesian Optimization Technique for Robust Design With Multiple Objectives, Failed Evaluations, and Stochastic Parameters

Wauters, Jolan

doi:10.1115/1.4064674

Abstract

In this work, the efficient robust global optimization (ERGO) method is revisited with the aim of enhancing and expanding its existing capabilities. The original objective of ERGO was to address the computational challenges associated with optimization-under-uncertainty through the use of Bayesian optimization (BO). ERGO tackles robust optimization problems which are characterized by sensitivity in the objective function due to stochasticity in the design space. It does this by concurrently minimizing the mean and variance of the objective in a multi-objective setting. To handle the computational complexity arising from the uncertainty propagation, ERGO exploits the analytical expression of the surrogate model underlying BO. In this study, ERGO is extended to accommodate multiple objectives, incorporate an improved predictive error estimation approach, investigate the treatment of failed function evaluations, and explore the handling of stochastic parameters next to stochastic design variables. To evaluate the effectiveness of these improvements, the enhanced ERGO scheme is compared with the original method using an analytical test problem with varying dimensionality. Additionally, the novel optimization technique is applied to an aerodynamic design problem to validate its performance.

Graphical Abstract Figure

View large Download slide

1 Introduction

Driven by the goal of reducing the computational burden associated with conventional approaches to optimization-under-uncertainty (OUU), recent years have witnessed the emergence of various surrogate-assisted (SA) techniques aimed at addressing this challenge [1–3]. This paper focuses on addressing this computational challenge within the context of robust design optimization (RDO), which seeks to find designs that exhibit robustness against small variations in uncertain parameters [4]. More specifically, the performance of the efficient robust global optimization (ERGO) scheme is enhanced, which is positioned in the OUU research field below.

The existing SA-RDO approaches can be categorized in two ways: first based on the manner by which the surrogate is introduced and second on the formulation of the RDO problem. Persson and Ölvander divided the former in two based on whether the surrogate is updated during the optimization or not [5]. In reference to the latter, Kanno distinguished three RDO formulations: worst-case optimization (minimization of the maximum objective evaluation), discrepancy optimization (minimization of the difference between objective value evaluations), and variance minimization (minimization of the moments of the objective probability density functions) [6]. Persson and Ölvander’s study illustrated the outperforming ability of the actively updating approach [5]. Ryan et al. put a multi-objective (MO) reformulation of the variance minimization forward as the most competitive formulation of the RDO problem [7].

Among the various surrogate-modeling techniques available, Gaussian process interpolation (GPI), also known as Kriging, has received significant attention [8]. This can be attributed to GPI’s flexibility as a parameter-free method, as well as its capability to quantify prediction uncertainty. Bayesian optimization (BO) builds on GPI’s through the introduction of an acquisition functions that are typically formulated in such a way that the design space is explored and the objective is exploited. The efficient global optimization (EGO) framework is the most well known in which the acquisition function corresponds to the expected improvement (EI) [9,10].

The MO–BO approach to RDO was pursued by Keane in which the uncertainty propagation was performed by means of sampling the surrogate to obtain the quantities of interest and optimization these by means of NSGA-II [11,12]. This approach was also followed by Cheng et al. [13]. The approach of Ribaud et al. [14] relies on a Taylor expansion in order to propagate mean and variance through the surrogate. However, none of these methods make active use of the underlying model uncertainty to ensure objective space exploration and therefore, the end result might strongly depend on the initial set of function evaluation on which the surrogate is built. Most recently, the ERGO scheme was devised to account for this problem by formulating an acquisition function that actively leverages on the surrogate-induced uncertainty [15]. The ERGO scheme builds on the EGO algorithm and further increases the efficiency of solving RDO problems by performing the uncertainty propagation (UP) to assess the robustness measures without additional experiments or sampling. In this paper, the ERGO scheme is improved upon by permitting the treatment of multiple objectives, by improving upon the assessment of the predictive uncertainty of the stochastic measures, by treating function evaluation failures, and by treating the occurrence of stochastic parameters along with stochastic design variables. The contributions are summarized as follows:

The original ERGO scheme relied on a reformulation of the Keane’s Euclidean EI [16] as the acquisition function that drives the selection of designs to be evaluated. However, the Euclidean EI has been proven to be sensitive to magnitude variations of the objectives and is furthermore limited to two objectives [17], corresponding to the mean and variance of a single function output. As such, the hypervolume EI is introduced in the ERGO framework to account for both shortcomings [18].
The original ERGO scheme made use of a first-order second moment (FOSM) approach [19] to assess the predictive variance (epistemic uncertainty [20]) of the predictions of mean and variance of the objective due to input variability (aleatoric uncertainty [20]). This low-order approach drove the acquisition function to under emphasize certain regions in the objective space. To account for this shortcoming, a second-order second moment (SOSM) approach is introduced to assess the predictive variance.
The original ERGO scheme did not treat the possible occurrence of objective function failures. The SAMURAI framework [21], which incorporated the ERGO scheme in the context of asynchronous reliability-based robust design optimization, accounted for this particularity by means of adding the prediction of failed evaluations to the training set, mimicking Forrester et al.’s missing-at-random (MAR) approach [22]. The introduction of the MAR approach in the original scheme can introduce excessive non-linearities in the objective space when the predictive uncertainty of the surrogate is substantial. This can lead to the prediction of non-existent multi-model objective spaces and regions with high predictive uncertainty. Consequently, the acquisition function may remain elevated throughout the optimization process. To account for this, a classifier is introduced to predict the feasible domain and guide the acquisition function, following the approach of Backov et al. [23].
The original ERGO scheme only permitted stochastic design variables. As such, the impact of parameter variability on the objective could not be directly assessed. To account for this defect, the treatment of stochastic parameters is examined in the form of a two-step acquisition function optimization.

The remainder of this paper is organized as follows. In Sec. 2 of this paper, a brief overview of the notations used throughout this work is introduced. In Sec. 3, the ERGO framework is briefly summarized after which the three improvements upon the strategy are presented in Sec. 4, namely the novel acquisition function which relies on the hypervolume expected improvement, the higher-order predictive variance assessment of the stochastic objective measures, the treatment of crash constraints by means of Gaussian process classification, and the treatment of stochastic parameters. In Sec. 5, the aforementioned improvements are brought together making up ERGO-II and compared with the original scheme on an analytical test problem. In Sec. 6, the improved approach is examined on an aerodynamic test case inspired by the aerodynamic design optimization discussion group [24]. Concluding remarks and outlooks are formulated in Sec. 7.

2 Notation

2.1 Sub- and Superscripts.

$x_{i}$ refers to the $i th$ sample of $X$ (⁠ $i = 1, \dots, n (X)$ ⁠) with $n (X)$ the number of elements in set $X$ ⁠, with $X$ a proper subset of vector space $X$ and $X$ a subset of $R^{d (X)}$ ⁠: $x \in X \subset X \subseteq R^{d (X)}$ ⁠. $d (X)$ is used to denote the dimensionality of the vector space, such that $x^{(i)}$ refers to the $i th$ component (⁠ $i = 1, \dots, d (X)$ ⁠) of design vector $x$ ⁠. Furthermore, with $x_{\sim i}$ the vector $x$ is denoted from which the $i th$ component has been omitted.

The superscripts ${\cdot}^{i}$ and ${\cdot}^{c}$ are used to differentiate the interpolator from the classifier. Similarly, the subscripts ${\cdot}_{x}$ and ${\cdot}_{u}$ are used to differentiate the design variable, entities that influence the model that the engineer can change within a specified range during the optimization design study, and the operation parameters, entities that influence the model and might vary during operation, but can not be changed by the engineer.

2.2 Operators.

The concept of vectorizing a matrix, denoted as $v e c (A) ≜ A^{V}$ ⁠, involves the process of stacking its columns. Similarly, the matricization of a tensor, expressed as $m a t (A) ≜ A^{M}$ ⁠, entails the stacking of matrices. In the case of a rank 3 tensor, $A \in R^{d} \times R^{d} \times R^{d}$ ⁠, this operation yields $m a t (A) \in R^{d} \times R^{d^{2}}$ ⁠, while for a rank 4 tensor, $A \in R^{d} \times R^{d} \times R^{d} \times R^{d}$ ⁠, it results in $m a t (A) \in R^{d^{2}} \times R^{d^{2}}$ ⁠.

Furthermore, the use of the double dot product, denoted as “:”, serves to represent the Frobenius inner product: $A : B =$ $t r (A^{⊤} \cdot B)$ ⁠, with $t r$ signifying the trace operator of a matrix. Additionally, a noteworthy observation is made: $A \cdot B : B \cdot A = t r ((A \cdot B)^{2})$ holds true when both $A$ and $B$ are symmetric. The Hadamard product is aptly denoted by “ $\circ$ ”, and the notation “ ${\cdot}^{\circ \cdot}$ ” is introduced to represent the Hadamard power. The Kronecker product is represented by “ $\otimes$ ”. Moreover, an “outer power” of a vector is introduced as $x^{\otimes} ≜ x \otimes x = x x^{⊤}$ ⁠, providing a means to simplify notation. Similarly, an “outer power” of a matrix is defined as $A^{\otimes} ≜ A \otimes A = v e c (A) v e c (A)^{⊤}$ ⁠. Based on these definitions, equivalent expressions are formulated, such as $x^{\otimes} : A = x^{⊤} A x =$ $| x ‖_{A}^{2} = t r (x^{\otimes} \cdot A^{⊤})$ ⁠, which are employed interchangeably for the sake of brevity and clarity.

Finally, use is made of $ϕ$ and $Φ$ to respectively denote the standard normal and cumulative distribution function. The notation of both is eased out by introducing $ϕ_{j} (x) ≜ ϕ (y_{j} | μ_{j} (x), Σ_{j} (x))$ in the context of Gaussian processes.

2.3 Fonts.

Epistemic uncertainty denotes non-deterministic behavior arising from a degree of ignorance or a deficiency in knowledge concerning the system or the environment and is depicted in this work using the calligraphic font, symbolized as $y$ ⁠. In contrast, aleatory uncertainty pertains to the inherent variability linked to the physical system or environment under investigation, which is in this work represented by the Fraktur typeface, denoted as $x$ [20].

3 Efficient Robust Global Optimization

The ERGO framework, upon which this paper builds, is a SA-RDO technique in which a MO–BO routine is combined with a second-order Taylor-based UP approach [15]. In the ERGO framework, the mean and variance of the objective function are minimized simultaneously through the formulation of an innovative acquisition function, referred to as the robust expected improvement (REI). The UP approach leverages the analytical formulation of the surrogate model to obtain a closed expression of the mean and variance and their respective predictive uncertainties. As such, no additional sampling is required, which further decreases the computational cost of solving the problem. The RDO problem solved by ERGO can formally be written as

X_{*} = \arg \underset{X \in X}{Par} {E_{x}, V_{x}} [y (X)]

(1)

where

E x

and

V_{x}

respectively denote the expected value and variance operators over the stochastic design variables

x

⁠, which are treated as aleatoric uncertainties. The RDO problem translates itself to finding the Pareto front

P ≜ Par {E_{x}, V_{x}} [y] ≜ {E_{x} [y_{*}], V_{x} [y_{*}]} : {E_{x} [y_{*}], V_{x} [y_{*}]} ≼ {E_{x} [Y], V_{x} [Y]}

with

Y

the objective space:

y_{*} \subset Y \subseteq R

⁠. Furthermore, it is assumed that the input variability is normally distributed

x \sim N (x, Φ_{x})

⁠, such that

x

can be fully described by its mean

x

and variance

Φ_{x}

⁠.

Making a prediction for a noisy input corresponds to the integration of the epistemic predictive distribution of the Gaussian process

p (y (x) | S^{i})

⁠, created from

S^{i} = {X, Y}

with

Y

the objective matrix defined by

Y^{(i, j)} = y_{i} (x_{j})

and aleatory input distribution

p (x | x, Φ_{x})

leading to mixed epistemic and aleatory distribution

p (y (x) | S^{i}, x, Φ_{x})

⁠. Evaluation of the mean and variance over the aleatoric variability leads to again epistemic stochastic processes.

\begin{aligned} y_{E}^{} (x | Φ_{x}) & ≜ E_{x} [y (x) | S^{i}, x, Φ_{x}] \\ = \int p (y (x) | S^{i}, x, Φ_{x}) \cdot p (x | x, Φ_{x}) \cdot d x \end{aligned}

(2a)

\begin{aligned} y_{V}^{} (x | Φ_{x}) & ≜ V_{x} [y (x) | S^{i}, x, Φ_{x}] \\ = \int (p (y (x) | S^{i}, x, Φ) - E_{x} [y (x) | S^{i}, x, Φ_{x}])^{2} \\ \cdot p (x | x, Φ_{x}) \cdot d x \end{aligned}

(2b)

as presented in Ref. [25].

A closed-form expression of $y_{E}^{}$ and $y_{V}^{}$ is not readily available, therefore approximations are resorted to. In the original framework, Girard et al.’s approach is followed for the universal Kriging interpolator [15], where the trends are obtained by means of Isserlis’ theorem and the processes are estimated through the assumption that they are again Gaussian, appropriately referred to as the Gaussian approximation [25]. Furthermore, the covariance is locally approximated by a second-order polynomial.

The more widely used ordinary GPI is presented here, in which the trend function

f

is assumed constant. In this case, the approximation becomes mathematically equivalent to locally approximating the predictive mean of the objective with a second-order Taylor function. In case of normally distributed input noise, the uncertainty propagation becomes analytically tractable and corresponds to the second-order first moment (SOFM) and SOSM approach for respectively the predictive means of the mean and variance of the objective due to input variability. These values are now fully defined by the predictive mean

μ_{y}

⁠, Jacobian

μ_{\nabla y}

and Hessian

μ_{Δ y}

of the stochastic process with deterministic input and the mean

x

and covariance

Φ_{x}

of the stochastic input.

\begin{aligned} μ_{y}^{T} (x) & = μ_{y}^{} (x) + (x - x)^{⊤} \cdot μ_{\nabla y}^{} (x) \\ + \frac{1}{2} t r ((x - x)^{\otimes} \cdot μ_{Δ y} (x)) \end{aligned}

(3a)

\begin{aligned} μ_{E}^{} (x | Φ_{x}) & ≜ E_{y} {y_{E}^{} (x | Φ_{x})} \\ \approx E_{x}^{} [μ_{y}^{T} (x) | x, Φ_{x}] \\ \approx μ_{y}^{} (x) + \frac{1}{2} Φ_{x} : μ_{Δ y}^{} (x) \end{aligned}

(3b)

\begin{aligned} μ_{V}^{} (x | Φ_{x}) & ≜ E_{y} {y_{V}^{} (x | Φ_{x})} \\ \approx V_{x}^{} [μ_{y}^{T} (x) | x, Φ_{x}] \\ \approx t r (μ_{\nabla y}^{\otimes} (x) \cdot Φ_{x}) + \frac{1}{2} t r ((μ_{Δ y}^{} (x) \cdot Φ_{x})^{2}) \end{aligned}

(3c)

as derived in Ref. [15].

To ensure exploration during the optimization process, the uncertainty in the prediction of the afore-derived stochastic metrics needs to be known. In a similar fashion as above, a Taylor-series based error propagation technique is employed. The original approach made use of a first-order curve and neglected the influence of covariance terms.

\begin{aligned} Σ_{E} (x | Φ_{x}) & ≜ V_{y} {y_{E}^{} (x | Φ_{x}) | S^{i}} \\ \approx V_{y} {E_{x} [μ_{y}^{T} (x) | x, Φ_{x}] | S^{i}} \\ \approx Σ_{y}^{} (x) + \frac{1}{2} Σ_{Δ y}^{M} (x) : Φ_{x}^{\otimes} \end{aligned}

(4a)

\begin{aligned} Σ_{V} (x | Φ) & ≜ V_{y} {y_{V}^{} (x | Φ_{x}) | S^{i}} \\ \approx V_{y} {V_{x} [μ_{y}^{T} (x) | x, Φ_{x}] | S^{i}} \\ \approx (μ_{\nabla y} (x)^{⊤} \cdot Φ_{x})^{\otimes} : Σ_{\nabla y} (x) \\ + (μ_{Δ y}^{V} (x)^{⊤} \cdot Φ_{x}^{\otimes})^{\otimes} : Σ_{Δ y}^{M} (x) \end{aligned}

(4b)

as derived in Ref. [15].

By means of these metric, a novel acquisition function can be formulated, namely the REI. The REI builds on Keane’s Euclidean EI [15,16]. In this context, the Pareto front

P

is redefined as the Pareto front of the aleatoric mean and variance as predicted by the surrogate in the evaluated infills:

P ≜ Par [μ_{{E, V}} (X | Φ_{x})]

⁠.

\begin{aligned} {REI}_{E u} (x | Φ_{x}, P) & ≜ P_{y} [y_{{E, V}} (x | Φ_{x}) < P] \\ \cdot ‖ {\bar{y}}_{{E, V}} (x | Φ_{x}) - P ‖_{2} \\ = ‖ {1^{⊤} q_{E}, 1^{⊤} q_{V}} - 1^{⊤} p \cdot P ‖_{2} \end{aligned}

(5)

where the first term denotes the probability of finding the infill point outside of the current predicted Pareto front and the second term denoting the distance of the centroid of the infill in the objective space to the nearest point on the predicted Pareto front.

p^{s}

⁠,

q_{E}^{s}

⁠, and

q_{V}^{s}

correspond to the sum over three vector functions

Algorithm 1

Require: Evaluated sampling plan $S^{i} = {X, y}$

1: While computational budget is not exhausted do

2: $θ_{*} \leftarrow \arg max_{θ} L^{i} (θ | S^{i})$ ▷Eq.(A2)

3: $P \leftarrow Par [μ_{{E, V}} (X | θ_{*}, Φ_{x})]$

4: $z_{*} \leftarrow max_{x} {REI}_{E u} (x | θ_{*}, P)$ ▷Eq.(5)

5: If $z_{*} > k_{R E I}$ then

6: $x_{*} \leftarrow \arg max_{x} {REI}_{E u} (x | θ_{*}, P)$

7: $y_{*} \leftarrow y (x_{*})$

8: $S^{i} \leftarrow S^{i} \cup {x_{*}, y_{*}}$

9: else

10: $break$

11: $X_{*} \leftarrow \arg {Par}_{x} [μ_{{E, V}} (X | θ_{*}, Φ_{x})]$

\begin{aligned} p & ≜ {[Φ_{E}^{(1)}, {(Φ_{E}^{(i + 1)} - Φ_{E}^{(i)}) \cdot Φ_{V}^{(i + 1)}}, (1 - Φ_{E}^{(p)}) \cdot Φ_{V}^{(p)}]}^{⊤} \\ q_{E} & ≜ {[Ω_{E}^{(1)}, {(Ω_{E}^{(i + 1)} - Ω_{E}^{(i)}) \cdot Φ_{V}^{(i + 1)}}, Ω_{E}^{(p)} \cdot Φ_{V}^{(p)}]}^{⊤} \\ q_{V} & ≜ {[Ω_{V}^{(1)}, {(Ω_{V}^{(i + 1)} - Ω_{V}^{(i)}) \cdot Φ_{E}^{(i + 1)}}, Ω_{V}^{(p)} \cdot Φ_{E}^{(p)}]}^{⊤} \end{aligned}

with

i = 1, \dots, p - 1

and

p = n (P)

⁠. Furthermore,

Φ_{E}^{(i)}

and

Ω_{E}^{(i)}

are introduced as

\begin{aligned} Φ_{E}^{(i)} & ≜ Φ (P_{E}^{(i)} | μ_{E}^{} (x | Φ_{x}), Σ_{E} (x | Φ)) \\ Ω_{E}^{(i)} & ≜ μ_{E}^{} (x | Φ_{x}) \cdot Φ_{E}^{(i)} - Σ_{E}^{} (x | Φ_{x}) \cdot ϕ_{E}^{(i)} \end{aligned}

as proposed in Ref. [15].

Plugged in EGO (cf. Appendix) the ERGO framework is obtained (Algorithm 1). Notice in Algorithm 1: line 3 that $θ_{*}$ is added as an argument. This substitution is made to indicate that a maximum likelihood estimation (MLE) approach is adopted instead of a fully Bayesian approach to evaluate the interpolator. Furthermore, $k_{R E I}$ is introduced as an additional constant that serves as a stopping criterion: if no further significant improvement can be made, stop (⁠ $b r e a k$ ⁠) the optimization.

4 Improved Efficient Robust Global Optimization

4.1 Robust Hypervolume Expected Improvement.

The initial ERGO scheme depends on adapting Keane’s Euclidean EI [16] as the acquisition function guiding the selection of designs for evaluation. Nevertheless, the Euclidean EI has been shown to be responsive to changes in the magnitudes of objectives and is also constrained to handling only two objectives [17], representing the mean and variance of a singular function output. Consequently, the ERGO framework introduces the hypervolume EI to address both of these limitations [18].

The current standard to assess the Pareto set quality is the hypervolume indicator

H (P)

⁠, which corresponds to the Lebesque measure of the hyperspace dominated by the Pareto front bounded by a reference point in the objective space

r ≜ {y_{j}^{max} + ϵ | j = 1 \dots d (Y)}

⁠. The indicator can be used to define a scalar improvement function

I_{h v} (x | P, r)

which measures the contribution (or improvement) of the point

y (x)

to the Pareto set

P

⁠. Identical to the single-objective case, the stochastic nature of the interpolator can be used to assess the uncertainty in the prediction and evaluated to define the expected hypervolume improvement

E_{y} [I_{h v} (x | Φ_{x}, P, r)]

⁠. To realize this, the exclusive hypervolume indicator

H^{e} (x | P, r) ≜ H (P \cup y (x) | r) - H (P | r)

is introduced, such that

I_{h v} (x | P, r) = max [H^{e} (x | P, r), 0]

(6a)

E_{y} [I_{h v} (x | P, r)] = \int I_{h v} (x | P, r) \prod_{j = 1}^{d (y)} ϕ_{j} (x) \cdot d y_{j}

(6b)

Extensive research has been conducted toward the analytical evaluation of this integral. In this work, the approach by Couckuyt et al. is built upon [18], in which

E_{y} [I_{h v} (x | P, r)]

is written in closed form as the sum of improvement contributions of disjoint cells

K = {κ_{i} | i = 1, \dots, n (K)}

with

κ_{i} = {[{\underline{κ}}_{i}^{(j)}, {\bar{κ}}_{i}^{(j)}] | j = 1, \dots, d (Y)}

making up the non-dominated region

I_{C} (K | x)

⁠. The overline and underline are introduced to respectively denote the lower and upper bounds of the cell, such that

κ_{i}

contains the dimensions of the cell.

E_{y} [I_{h v} (x | P, r)] = \sum_{i = 1}^{n (K)} I_{C} (κ_{i} | x)

(7)

The improvement contribution of each cell

κ_{i}

is by itself the sum over the disjoint cells

K

of the product over the dimensionality of the objective space of a function

J

I_{C} (κ_{i} | x) = \sum_{j = 1}^{n (K)} \prod_{l = 1}^{d (Y)} J^{(l)} (κ_{i}, κ_{j} | x)

(8)

depending on the relative position of

κ_{i}

versus

κ_{j}

⁠: if the cell that contains the prediction dominates the cell in the dimension considered,

J_{1}

is invokes. When this domination is only partial,

J_{2}

is called. In all other cases, the contribution is zero. This can be mathematically presented as

\begin{aligned} J^{(l)} (κ_{i}, κ_{j} | x) \\ = {\begin{matrix} J_{1}^{(l)} & if {\underline{κ}}_{j}^{(l)} \geq {\bar{κ}}_{i}^{(l)} \\ J_{2}^{(l)} & if {\bar{κ}}_{j}^{(l)} > {\underline{κ}}_{i}^{(l)} {\bar{κ}}_{j}^{(l)} < {\underline{κ}}_{i}^{(l)} \\ 0 & otherwise \end{matrix} \end{aligned}

(9)

The mathematical derivation of

J_{1}

and

J_{2}

can be found in Couckuyt et al. [18] and is given below. The superscript

\cdot^{(l)}

has been omitted from

J

and

κ

for clarity.

\begin{aligned} J_{1} & = ({\bar{κ}}_{j} - {\underline{κ}}_{j}) \cdot {Φ_{l} ({\bar{κ}}_{i}) - Φ_{l} ({\underline{κ}}_{i})} \end{aligned}

(10a)

\begin{aligned} J_{2} & = ({\bar{κ}}_{j} - {\underline{κ}}_{j}) \cdot {Φ_{l} (max [{\underline{κ}}_{i}, {\underline{κ}}_{j}]) - Φ_{l} ({\underline{κ}}_{i})} \\ + ({\bar{κ}}_{j} - μ_{l} (x)) \cdot {Φ_{l} (min [{\bar{κ}}_{i}, {\bar{κ}}_{j}]) - Φ_{l} (max [{\underline{κ}}_{i}, {\underline{κ}}_{j}])} \\ + Σ_{l} (x) \cdot {ϕ_{l} (min [{\bar{κ}}_{i}, {\bar{κ}}_{j}]) - ϕ_{l} (max [{\underline{κ}}_{i}, {\underline{κ}}_{j}])} \end{aligned}

(10b)

To calculate the statistical functions, the integral bounds of the region that improves the Pareto set need to be identified. However, this region is often irregularly shaped, making analytical calculation complex for a higher number of objective functions. Instead, Couckuyt et al. proposed a branch-and-bound approach. The algorithm constructs a pseudo Pareto set by adding two extra points representing the ideal and anti-ideal points. The objective space is divided into hyperrectangles, and a binary partitioning strategy is used to determine the smallest set of hyperrectangles that dominate or augment the Pareto set. The process involves testing a single point within each hyperrectangle and recursively partitioning the space until all hyperrectangles are accepted or rejected based on a condition function. Once the integral bounds are identified, the statistical criteria can be evaluated [18]. This is schematically presented in Fig. 1.

Fig. 1

View large Download slide

Schematic representation of the partitioning of the area outside of the Pareto front (bottom-left) of two contesting objectives in disjoint cells. The light and dark areas (bottom-left) respectively correspond to the area where the Pareto front can be improved and the exclusive hypervolume for infill $y (x)$ [18].

When the expected value and variance of the objective functions are used to drive the expected hypervolume improvement the REI is obtained

{REI}_{h v} (x | Φ_{x}, P, r) ≜ E_{y} [I_{h v} (x | Φ_{x}, P, r)]

(11)

As indicated by Wagner et al., the activation function that drives the Bayesian optimization routine should decrease to zero on evaluated points [17]. This typically corresponds to the predictive variance decreasing to zero in evaluated points. However, the predictive (epistemic) variances of the (aleatoric) mean and variance do not disappear in the infills due to contribution of the predictive variance of the Jacobian and Hessian. To avoid reevaluation of infills and stagnation of the optimizer an activation function

A (x | X)

is introduced and multiplied with the acquisition function. The activation function takes on the form of a Gaussian radial basis function, which is reformulated such that it is equal to zero in evaluated points and continuously increases to one in unevaluated regions.

{REI}_{h v} (x | Φ_{x}, P, r) ≜ E_{y} [I_{h v} (x | Φ_{x}, P, r)] \cdot A (x | X)

(12a)

A (x | X) ≜ 1 - \sum_{i = 1}^{n (X)} \exp (- w | x_{i} - x ‖_{2}^{2})

(12b)

with the weighting coefficient

w

empirically defined to be equal to

(N_{S^{i}} + 50)^{2}

⁠.

4.2 Improved Epistemic Uncertainty Propagation.

The original ERGO scheme made use of an FOSM approach to assess

Σ_{E}

and

Σ_{V}

⁠. This low order approach drove the acquisition function to under emphasize certain regions in the objective space. To account for this shortcoming, an SOSM approach is introduced to assess the predictive variance [26].

\begin{aligned} Σ_{E}^{} (x | Φ_{x}) & \approx Σ_{y} (x) + Σ_{y, δ y} (x) : Φ_{x} + \frac{1}{2} Σ_{δ y}^{M} (x) : Φ^{\otimes} \end{aligned}

(13a)

\begin{aligned} Σ_{V}^{} (x | Φ_{x}) & \approx (μ_{\nabla y} (x)^{⊤} \cdot Φ_{x})^{\otimes} : Σ_{\nabla y} (x) + \frac{1}{2} t r ((Σ_{\nabla y} (x) \cdot Φ_{x})^{2}) \\ + (μ_{Δ y}^{V} (x)^{⊤} \cdot Φ_{x}^{\otimes})^{\otimes} : Σ_{Δ y}^{M} (x) + \frac{1}{2} t r ((Φ_{x}^{\otimes} \cdot Σ_{Δ y}^{M} (x))^{2}) \\ + 2 (μ_{\nabla y} (x)^{⊤} \cdot Φ_{x})^{⊤} \cdot (μ_{Δ y}^{V} (x)^{⊤} \cdot Φ_{x}^{\otimes})^{⊤} : Σ_{\nabla y, Δ y} (x) \\ + Σ_{\nabla y, Δ y} (x) \cdot Φ_{x}^{\otimes} : Φ_{x} \cdot Σ_{\nabla y, Δ y} (x) \end{aligned}

(13b)

When comparing the first-order uncorrelated (Eq. ) and second-order correlated (Eq. ) predictive variance of the aleatoric mean

Σ_{E}

the additional covariance term can be observed. Similarly, when comparing the first-order uncorrelated (Eq. ) and second-order correlated (Eq. ) predictive variance of the aleatoric variance

Σ_{V}

additional second-order terms can be observed on top of the additional covariances terms. A comparison of the two approximations is presented in Fig. 2. While

Σ_{V}

remains small between 0.2 and 0.7, the second-order correlated approach predicts a non-zero value as opposed to the first-order uncorrelated approach. These new approximations can now be introduced in Eqs. and to improve upon the prediction of REI.

Fig. 2

View large Download slide

Comparison of the first-order uncorrelated and second-order correlated predictive variance of the aleatoric mean and variance of the Gaussian process predictor trained from eight equidistant samples for test function $y (x) = (6 x - 2)^{2} \sin (12 x - 4) / 45 + 0.03$ with $x \sim N (x, 10^{- 2})$

4.3 Crash Constraint Treatment.

The original ERGO scheme did not treat the possible occurrence of objective function failures. The SAMURAI framework, which incorporated the ERGO scheme in the context of reliability-based robust design optimization, accounted for this particularity by means of adding the prediction of failed evaluations to the training set. However, this strongly affects the convergence of the scheme [21].

By storing information on whether (⁠ $c_{i} = + 1$ ⁠) or not (⁠ $c_{i} = - 1$ ⁠) the design vector (⁠ $x_{i}$ ⁠) could be evaluated, a binary classification problem is obtained. In this work, use is made of a discriminative Gaussian process classifier (GPC) to model the feasible space directly $p (c (x) | S^{c})$ with $S^{c} = {X, c} = {(x_{i}, c_{i}) | i = 1, \dots, n (S^{c})}$ ⁠. Rasmussen and Williams’ textbook is strongly drawn from to describe the theory behind GPC and its implementation [27]. The approach corresponds to turning the output of a regression model into a class probability using a response function $λ (z) : R \mapsto [0, 1]$ ⁠, guaranteeing a valid probabilistic interpretation. As such, the regression model becomes a latent or hidden process.

Inference can as such be understood as a two-step process, in which first the distribution over the latent output is determined

p (y (x) | S^{c})

⁠, after which this distribution is used to make a prediction of the probability of a class

p (c (x) | S^{c})

⁠. To ease out the notation, the arguments of the joint distribution due to the evaluation of the training set in the latent stochastic process are dropped:

y (X^{c}) ≜ y

⁠,

p (c (x) | S^{c}) = \int λ (y (x)) \cdot p (y (x) | S^{c}) \cdot d y (x)

(14a)

p (y (x) | S^{c}) = \int p (y (x) | y) \cdot p (y | S^{c}) \cdot d y

(14b)

with

p (y | S^{c}) = p (c | y) \cdot p (y) / p (c)

⁠. The evaluation of these integrals is intractable and must be approximated. Conventional methods rely on sampling, using for example, Markov Chain Monte Carlo sampling, or analytical approximations of the integrals, such as the Laplace approximation, variational inference (VI) or expectation propagation (EP). In this work, use is made of the latter, which is briefly described here. Through Bayes theorem the following can be observed

p (y | S^{c}) \propto p (y) \prod_{i = 1}^{n (S^{c})} p (c_{i} | y_{i})

(15)

with

p (y)

the Gaussian prior. Note furthermore that the proportionality is exactly equal when a normalization term, corresponding to the marginal likelihood, is taken into account. The likelihood

p (c_{i} | y_{i})

is taken equal to the probit likelihood for binary classification

Φ (c_{i} y_{i})

⁠. The latter is approximated by means of an unnormalized Gaussian function

t_{i} (y_{i})

⁠, such that posterior

p (y | S^{c})

can be approximated by

q (y | S^{c})

q (y | S^{c}) ≜ \frac{1}{k_{E P}} p (y) \prod_{i = 1}^{n (S^{c})} t_{i} (y_{i} | {\tilde{k}}_{i}, {\tilde{μ}}_{i}, {\tilde{σ}}_{i}^{2})

(16a)

= N (μ, Σ)

(16b)

with

\tilde{Σ} = d i a g ({\tilde{σ}}_{i}^{2})

⁠,

μ = Σ^{- 1} \tilde{Σ} \tilde{μ}

⁠,

Σ = (Ψ^{- 1} + {\tilde{Σ}}^{- 1}),

and

k_{E P} = q (S^{c})

⁠. The latter corresponds to the normalization term of the posterior and is equal to the EP approximation of the marginal likelihood, all of which are predicted by the EP algorithm. The marginal log likelihood

\log L^{c} (θ | S^{c}) ≜ \log p (S^{c} | θ) \approx \log q (S^{c} | θ)

is given by

\begin{aligned} \log L^{c} (θ | S^{c}) & = - \frac{1}{2} \log | Ψ + \tilde{Ψ} | - \frac{1}{2} {\tilde{μ}}^{⊤} (Ψ + \tilde{Ψ})^{- 1} \tilde{μ} \\ + \sum_{i = 1}^{n (S^{c})} \log Φ (\frac{c_{i} μ_{- i}}{\sqrt{1 + σ_{\sim i}^{2}}}) + \sum_{i = 1}^{n} \log (σ_{\sim i}^{2} + {\tilde{σ}}_{i}^{2}) \\ + \sum_{i = 1}^{n (S^{c})} \log \frac{(μ_{\sim i} - {\tilde{μ}}_{i})^{2}}{2 (σ_{\sim i}^{2} + {\tilde{σ}}_{i}^{2})} \end{aligned}

(17)

The EP algorithm relies on the minimization of the Kullback–Leibler divergence of the exact posterior over the approximate

KL (p (y | S^{c}) ‖ q (y | S^{c}))

⁠, as opposed to the VI algorithm which minimizes

KL (q (y | S^{c}) ‖ p (y | S^{c}))

⁠. This translates itself in a four-step iterative process during which

t_{i}

is sequentially updated. First, an approximate posterior is constructed by leaving out the current

t_{i}

⁠. Second, this distribution is combined with the exact likelihood to approximate the marginal distribution. This marginal is approximated by a Gaussian distribution in the third step, which boils down to a moment matching step. Finally,

t_{i}

is reevaluated in the fourth step. Details on the algorithmic implementation are given in Rasmussen and Williams [27].

The availability of the approximate Gaussian posterior distribution obtained by means of the EP algorithm now permits the evaluation of the approximate predictive distribution.

q (c (x) = 1 | S^{c}) = Φ (\frac{ψ (x)^{⊤} (Ψ + \tilde{Σ})^{- 1} \tilde{μ}}{\sqrt{1 + ψ (x, x) - ‖ ψ (x) ‖_{(Ψ + \tilde{Σ})^{- 1}}^{2}}})

(18)

An illustration of the classifier is presented in Fig. 3. The surrogate-assisted SOFM (Eq. ) and SOSM (Eq. ) predictions are cut off by

q (c (x) = 1 | S^{c}) < 0.5

⁠. This ensures that optimizer is not motivated to explore the non-evaluable region. Notice that the prediction of the unfeasible domain only slightly moves into the exact one, which indicates the effectiveness of the GPC.

Fig. 3

View large Download slide

Comparison of the surrogate-assisted SOFM (Eq. (3a)) and SOSM (Eq. (3b)) predictions with the exact mean and variance of the problem presented in Fig. 2 for which the evaluations fail at $x < 0.2$ ⁠. The predictions in by means of a GPC modeled non-evaluable region are omitted.

The treatment of objective function evaluation failure in the context of ERGO can subsequently be understood as follows. First, a GPC is trained in parallel with the GPIs of the objective functions evaluations. Subsequently, the acquisition function is conditioned to only probe the design space which is predicted to be evaluable. This is realized by taking the product of the acquisition function and the predictive classifier. The result is a “safe” robust expected improvement (SREI):

\begin{aligned} SREI (x | Φ_{x}, P, r) & ≜ E_{y} [I_{h v} (x | Φ_{x}, P, r)] \cdot A (x | X) \\ \cdot q (c (x) = 1 | S^{c}) \end{aligned}

(19)

The effect of the activation function

A (x | X)

(Eq. ) and the classifier

q (c (x) = 1 | S^{c})

(Eq. ) on the REI is illustrated in Fig. 4. It can be observed that dotted line is high in the unfeasible region, even for designs that have been attempted to be evaluated. However, by enforcing the activation function (dashed line), regions of (attempted) evaluations are forced to zero. Secondly, the full line shows that by applying the classifier, the most valuable next point to be evaluated lies outside the non-evaluable region. Proving the effectiveness of the proposed new acquisition function.

Fig. 4

View large Download slide

Illustration of the effect of the activation function $A (x | X)$ (Eq. (12a)) and the classifier $q (c (x) = 1 | S^{c})$ (Eq. (18)) on the REI for the problem presented in Figs. 2 and 3

The construction of the Gaussian process classifier is performed using the publicly available matlab® implementation of Rasmussen and Williams’ work using a conjugate gradient approach to maximize the concentrated log marginal likelihood function [27].

4.4 Stochastic Parameter Treatment.

The original ERGO scheme only permitted stochastic design variables. As such, the impact of parameter variability on the objective could not be directly assessed. To account for this defect, the treatment of stochastic parameters was included in the improved ERGO formulation. The problem solved now takes on the following form:

X_{*} = \arg \underset{X \in X}{Par} {E_{u, x}, V_{u, x}} [y (X | u)]

(20)

The Gaussian processes built from the different objective functions now include the stochastic parameters

u \in R^{d (U)}

with the assumption that

u \sim N (u, Φ_{u})

⁠. More concretely, a space filling design of experiments including

u

is constructed where

u \in [u - 3 σ_{u}, u + 3 σ_{u}]

with

d i a g [Φ_{u}] = σ_{u}

⁠. This ensures that

99.6 %

of the effect of stochastic parameters is captured. The robust Pareto front is now updated to include the effect of the stochastic parameters, evaluated at their mean value

P ≜ Par [μ_{{E, V}} (X | u, Φ)]

⁠.

To select a new infill, the robust expected improvement is optimized around the mean of the stochastic parameters. The corresponding new infill design variable vector

x_{*}

is subsequently held fixed while

u_{*}

is determined. The latter corresponds to the maximizer of the “safe” predictive variance product (SPVP) as to ensure that the predictions of aleatoric metrics is as accurate as possible. This results in the following two-step optimization process:

x_{*} = \arg max_{x \in X} SREI (x | u, P, r)

(21a)

\begin{aligned} u_{*} = \arg max_{u \in U} SPVP (u | x_{*}) \end{aligned}

(21b)

SPVP (u | x_{*}) ≜ \prod_{i = 1}^{d (y)} Σ_{i} (u | x) \cdot q (c (u | x) = 1 | S^{c})

(21c)

4.5 Algorithmic Development.

Algorithm 2

Require $S = {S^{i}, S^{c}}$ with $S^{i} = {U^{i}, X^{i}, Y^{i}}$ and $S^{c} = {U^{c}, X^{c}, c}$

1: While computational budget is not exhausted do

2: For $k = 1, \dots, n (Y)$ do

3: ${Θ^{i}}_{*} (k) \leftarrow \arg max_{θ} L^{i} (θ | {S^{i}}_{k})$ ▷Eq.(A2)

4: $θ_{*}^{c} \leftarrow \arg max_{θ} L^{c} (θ | S^{c})$ ▷Eq.(17)

5: $P \leftarrow Par [μ (X | u, Θ_{*}, Φ)]$

6: $z_{*} \leftarrow max_{x} SREI (x | u, Θ_{*}, P, r)$ ▷Eq.(21a)

7: If $z_{*} > k_{R E I}$ then

8: $x_{*} \leftarrow \arg max_{x} SREI (x | u, Θ_{*}, P, r)$

9: $u_{*} \leftarrow \arg max_{u} SPVP (u | x_{*}, Θ_{*})$ ▷Eq.(21b)

10: $c_{*} \leftarrow c (u_{*}, x_{*})$

11: If $c_{*} = 1$ then

12: $y_{*} \leftarrow y (u_{*}, x_{*})$

13: $S^{i} \leftarrow S^{i} \cup {u_{*}, x_{*}, y_{*}}$

14: $S^{c} \leftarrow S^{c} \cup {u_{*}, x_{*}, c_{*}}$

15: else

16: $break$

17: $X_{*} \leftarrow \arg {Par}_{x} [μ (X | u, Θ_{*}, Φ)]$

If these different aspects are brought together, the improved ERGO-II algorithm is obtained (Algorithm 2 and Fig. 5) with $Θ = {Θ^{i}, θ^{c}}$ and $Φ = {Φ_{x}, Φ_{u}}$ ⁠. Upon inspection of the new algorithm, it can be observed that in comparison with the original ERGO scheme (or EGO upon which it was built as presented in Algorithm 3 in Appendix) multiple GP interpolators are built for each of the different function output sets and a GP classifier to model the evaluable region (lines 2-4). Upon assessment of whether the contribution of the design that is expected to contribute most strongly to the robust Pareto is significant enough (lines 6 and 7), the design is evaluated (line 10). When the design evaluation is success, the output is added to the sets (lines 11-14). This process is subsequently repeated until the computational budget is depleted (line 1) or until the improvement decreases beyond a threshold value $k_{R E I}$ (lines 7 and 16).

Fig. 5

View large Download slide

Flowchart of the EGO algorithm in gray [10], with the additions that lead to the ERGO algorithm in blue [15] and the additions of this paper in orange resulting in the current ERGO-II scheme (Color version online.)

5 Performance Assessment

5.1 Problem Description.

A comparative study between ERGO (Algorithm 1) and ERGO-II (Algorithm 2) is setup to examine how the improved methodology performs in comparison to the original scheme. This is realized by applying the methodologies on a stochastic reformulation of Zitzler et al.’s first test problem (ZDT1) for varying dimensionality [28]. The test problem is restricted to a single stochastic function, to which two robustness objectives correspond, since the original ERGO scheme is restricted to this condition. Second, the original ERGO scheme is enhanced with the same routine to treat stochastic parameters as the ERGO-II method. Third, the MAR approach introduced in the SAMURAI framework [21] is also applied here to the original routine to permit the occurrence of evaluation failures in the comparative study, which corresponds to the constraint introduced here.

\begin{aligned} X_{*} = \arg \underset{x \in X}{Par} & {E_{u, x}, V_{u, x}} [y (X | u)] \\ s . t . & x^{(1)} \geq 0.2 \\ y (X | u) = 1 + 1.5 \sum_{i = 2}^{d (x)} x^{(i)} \\ - {u x^{(1)} (1 + 1.5 \sum_{i = 2}^{d (x)} x^{(i)})}^{0.5} \end{aligned}

(22)

with both the design variables and parameters being treated as being normally distributed

x \sim N (x, d i a g [10^{- 2}])

with

x \in [0, 1]^{d (X)}

and

u \sim N (1, d i a g [10^{- 2}])

⁠. The feasible domain is chosen such that it contains

80 %

of the input space.

5.2 Methodology.

The evaluation of the performance of the multi-objective optimization algorithms is conducted through the utilization of the hypervolume metric $H (P)$ [29]. The inherent stochastic nature of the optimization process, influenced by the random seed of the optimizer, can result in variable performance. To account for this aspect, the optimization is repeated ten times after which the average, min and max $H (P)$ in each iteration after the design of experiments (DoE) are computed.

The determination of the appropriate size for the DoE adheres to Jones et al.’s recommendation, which suggests a scaling factor of $11 \cdot (d (X) + d (U)) - 5$ [10]. For the DoE, a Latin hypercube sampling approach [30] is employed, and Morris and Mitchell’s maximin-criterion is utilized to quantify the space-filling property. This criterion aims to maximize the distances between pairs of points in ascending order while simultaneously minimizing the number of corresponding pairs [31]. The stopping criteria employed include a maximum of 200 function evaluations after the DoE or the termination of the process if the normalized expected improvement falls below 0.1% (corresponding to $k_{R E I}$ ⁠). The optimization of the REI is accomplished using a genetic algorithm.

5.3 Results and Discussion.

Figure 6 illustrates the ERGO versus ERGO-II performance assessment as a function of number of iterations for increasing dimensionality over the different plots from top to bottom with on the left the hypervolume indicator and on the right the normalized maximum expected improvement. The full line denotes the mean values, while the shades area are minimum, and maximum values of the respective performance measures. Optimal performance of an optimizer would be characterized by an early climb of $H (P)$ and a decrease to zero of the REI once $H (P)$ stagnates. The latter would ensure that the optimizer terminates if no further improvement can be obtained.

Fig. 6

View large Download slide

Performance assessment of ERGO versus ERGO-II for increasing dimensionality with on the left, the hypervolume indicator and on the right, the normalized maximum expected improvement

A notable observation is the enhanced performance of ERGO-II compared to the original scheme as the dimensionality increases (⁠ $d_{x} \in [5, 8]$ ⁠). This can be observed based on the fact that there is a strong increase in $H (P)$ early on and a decrease of REI to zero once $H (P)$ stagnated. ERGO maintains a high REI during the optimization, without further improvements of $H (P)$ ⁠. The significant performance fluctuations of the original scheme are attributed to the treatment of crash constraints, since this approach might lead to the constant generation of regions of high predictive variance, as discussed in the introduction.

At lower dimensionality (⁠ $d_{x} \in [2, 4]$ ⁠) ERGO-II also performs better than ERGO, but fails to stop once $H (P)$ stagnates. The reason for which might be attributed to the size DoE, which is now taken to be proportional to the dimensionality of the input space.

6 Validation

To illustrate the method in a more practical context, the original ERGO scheme is compared with the improved ERGO-II methodology by applying them to a computational fluid dynamics (CFD) design problem. More precisely, the stochastic reformulation of the second benchmark case introduced by the AIAA Aerodynamic Design Optimization Discussion Group (ADODG) is solved [24].

6.1 Problem Description.

The second benchmark case has the objective of minimizing the drag coefficient

C_{D}

of a modified RAE2822 profile in transonic viscid flow by altering its geometry, while a minimum cross-sectional area

S

⁠, lift coefficient

C_{L}

⁠, and pitching moment coefficient

C_{M}

constraint is satisfied [32].

\begin{aligned} x = \arg min_{\begin{matrix} x \end{matrix} \in X} & C_{D} (x) \\ s . t . & C_{L} (x) = 0.824 \\ C_{M} (x) \geq - 0.092 \\ S (x) \geq S_{R A E 2822} \end{aligned}

(23)

The problem is formulated at a freestream Mach number Ma

= 0.734

and Reynolds number Re

= 6.5 \times 10^{6}

⁠. This problem is extended upon by including stochasticity in the design and parameter space in order to account for off-point importance as was emphasized by LeDoux and colleagues [33], echoing Drela’s remarks on airfoil design [34].

\begin{aligned} X_{*} = \arg {Par}_{X \in X} & {E_{u, x}, V_{u, x}} [C_{D} (X | u)] \\ s . t . & S (X) \geq S_{R A E 2822} \\ C_{L} (X | u) \geq 0.824 \\ C_{M} (X | u) \geq - 0.092 \end{aligned}

(24)

The stochasticity in the parameter space is introduced in the form of variability in the Mach and Reynolds numbers

u = [Ma, Re]

such that these numbers are normality distributed around the original reference operation conditions with a coefficient of variation

σ_{u} / u

equal to 0.1.

6.2 Parameterization.

Kulfan and Bussoletti’s class-shape function transformation (CST) is used to parametrize the airfoil [35]. This corresponds to the thickness distribution

t (x_{c})

as the product of a class function

C (x_{c})

⁠, a generalized formulation of a circle, and a shape function

S (x_{c})

⁠, corresponding to the Bézier curve, plus a trailing edge (te) thickness correction

δ t_{t e}

to give the thickness distribution, with

x_{c} = x / c

where

x

and

c

respectively denote the chordwise position and chord length

t^{'} (x_{c}) = C (x_{c}) \cdot S (x_{c}) + x_{c} \cdot δ t_{t e}

(25a)

C (x_{c}) = {(x_{c})}^{0.5} \cdot (1 - x_{c})

(25b)

S (x_{c}) = \sum_{i = 0}^{d (b)} [b^{(i)} \cdot C_{i}^{d (b)} \cdot {(x_{c})}^{i} \cdot {(1 - x_{c})}^{d (b) - i}]

(25c)

with

C_{i}^{d (b)}

the binomial coefficient and

b

the Bézier coefficients. Since the airfoil is not symmetric, a separate curve is considered for both the pressure and suction side. Three second-order curves for both sides are considered, leading to six Bézier coefficients: [

b_{p}, b_{s}

]. Additionally, the angle of attack (AoA)

α

is added as a design variable:

x = [b_{p}, b_{s}, α]

⁠. The design space

X

is determined by estimating the CST coefficients corresponding to the RAE2822 (given in Table 1, obtained by means of an MSE minimization) around which a variation of

\pm 0.1

is permitted. The AoA is permitted to vary in

[0 \deg, 10 \deg]

⁠.

Table 1

Optimal designs CST description as predicted by ERGO (middle) and ERGO-II (bottom) in comparison with the default (RAE2822) profile at the top

	$b_{p}^{(1)}$	$b_{p}^{(2)}$	$b_{p}^{(3)}$	$b_{s}^{(1)}$	$b_{s}^{(2)}$	$b_{s}^{(3)}$
Default	$0.122$	$0.184$	$0.207$	$- .108$	$- .248$	$.035$
Design 1	$0.075$	$0.158$	$0.256$	$- .120$	$- .217$	$.001$
Design 2	$0.090$	$0.183$	$0.222$	$- .066$	$- .207$	$.010$
Design 1	$0.090$	$0.186$	$0.256$	$- .138$	$- .218$	$.007$
Design 2	$0.097$	$0.195$	$0.254$	$- .147$	$- .213$	$.002$

	$b_{p}^{(1)}$	$b_{p}^{(2)}$	$b_{p}^{(3)}$	$b_{s}^{(1)}$	$b_{s}^{(2)}$	$b_{s}^{(3)}$
Default	$0.122$	$0.184$	$0.207$	$- .108$	$- .248$	$.035$
Design 1	$0.075$	$0.158$	$0.256$	$- .120$	$- .217$	$.001$
Design 2	$0.090$	$0.183$	$0.222$	$- .066$	$- .207$	$.010$
Design 1	$0.090$	$0.186$	$0.256$	$- .138$	$- .218$	$.007$
Design 2	$0.097$	$0.195$	$0.254$	$- .147$	$- .213$	$.002$

Production tolerances are translated in the form stochastic design variables such that the design vector is normally distributed around the parametrization of the airfoil with a coefficient of variation equal to 0.1.

Self-intersecting airfoils that would appear in the DoE are accounted for by the crash constraint treatment presented in Sec. 4.3.

6.3 Methodology.

The treatment of constraints is performed through a conventional $p r o b a b i l i t y o f f e a s i b i l i t y$ approach [36]. To model the transonic flow around the airfoil, CFD simulations are applied. The governing equations for flows correspond to the unsteady Reynolds-averaged Navier–Stokes equations.

6.3.1 Turbulence Modeling.

The introduction of Reynolds’ decomposition of the flow variables in a mean and fluctuating part in the Navier–Stokes equations and the ensemble-averaging of this leads to the RANS equations, which are much cheaper to solve. However, by decomposing and averaging an additional term is introduced in the equations, which has the form of a stress-term and is typically addressed as the Reynolds stress. This leads to a set of equations that is no longer closed. Different manners exist by which this stress term is solved, such as by relating the Reynolds stress to the mean strain rate, building forth on the Boussinesq eddy viscosity assumption, by means of the eddy viscosity. Again different approaches exist to determine this term; the $k - ω$ SST model, which is used in this work, does this by means of two additional transport equations for turbulent kinetic energy $k$ and specific turbulence dissipation rate $ω$ and assuming a linear relation in the boundary layer between the turbulent shear stress and the turbulent kinetic energy [37].

6.3.2 Solver-Specific Settings.

To solve the set of equations presented in the previous sections, a pressure-based solver with a fully-coupled approach was applied. A second-order scheme was used for the pressure and second-order upwind discretization was applied for density and momentum while turbulent kinetic energy, specific dissipation rate, and energy were discretized using a first-order upwind scheme. A least-squares cell-based approach is employed to calculate the gradients. This implies that a least-squares method is used to solve the system of equations relating the cell gradient to the difference in cell values. All residuals of the equations were enforced to drop below $10^{- 5}$ ⁠, including for the continuity equation. The residuals of the former are expressed as the ratio of the sum of the absolute errors in each cell to the sum of the absolute values of the flow variable calculated in each cell. The residual of the latter is expressed as the ratio of the sum of the rate of mass creation in each cell to the largest absolute value of this sum in the first five iterations.

6.3.3 Mesh Settings.

A c-shaped grid is constructed around the airfoil. A grid independence study was conducted to determine the grid and domain sizing in accordance to the approach in [38]. The results correspond to a computational of 20 chord lengths in front, above, and below the airfoil and 40 chord lengths behind it as presented in Fig. 7. Furthermore, a maximum wall normal expansion ratio of 1.1 is imposed on 120 layers surrounding the airfoil to obtain a $y^{+}$ of maximum 1 near the stagnation point, on average 0.35, and 200 nodes are placed on the chord.

Fig. 7

View large Download slide

Computational domain for the high-fidelity solver built in accordance to the sensitivity study performed in Ref. [38]

6.3.4 Boundary Conditions.

Surrounding the computational domain a pressure far field boundary condition is imposed with a Mach number determined by the ERGO-II scheme. The values of the turbulent intensity and the turbulent viscosity ratio are respectively set to $5 %$ and $10 %$ ⁠. The wall is modeled as a no-slip boundary.

6.3.5 Computational Solver.

All calculations are performed using the CFD-code ansys®fluent®23.1.0 The computational effort of a single simulation using five cores of an Intel Xeon Gold 6240 2.6GHz CPU amounted to roughly 1 min.

6.4 Results and Discussion.

The resulting Pareto fronts of the design optimization-under-uncertainty as presented by ERGO and ERGO-II are visualized in Fig. 8. In light and dark gray, the feasible designs of respectively the DOE and the infills are presented. In black the Pareto fronts, as predicted by the surrogates of the output (⁠ $C_{D}$ ⁠) after the final iteration are visualized.

Fig. 8

View large Download slide

Pareto front of the ADODG2 aerodynamic design problem under uncertainty as predicted by ERGO (a) and ERGO-II (b))

Inspection shows that the Pareto front obtained by ERGO-II extends deeper in the bottom left than the one obtained by ERGO, conform the conclusions drawn in the previous section. Additionally, the number of feasible infills is larger when using ERGO-II, as can be seen by the larger number of dark gray dots in the bottom plot. This can again be attributed to the effect of the crash constraint treatment as pinpointed in the introduction and elaborated in Sec. 4.3: using the MAR approach might lead to distorted surrogates if the prediction is excessive or might mislead the optimizer if the MAR approach predicts the unfeasible region to be the optimal one. The behavior is avoided when using the GP classifier approach as presented in this work.

Another observation that can be made is the fact that the default prediction is different in the top and bottom plots. This can be attributed to the fact that the uncertainty propagation relies on the GPI. Theoretically, if the output function displays the smoothness characteristics enforced by the GPI and UP approach, the GP-assisted UP will converge to the exact solution. However, this is prohibited by the infills added using the MAR approach.

The two designs on the two Pareto fronts are selected for closer inspection. Their corresponding CST-parametrization values are given in Table 1 and visualized along with their characteristics for Ma $= 0.734$ and Re $= 6.5 \times 10^{6}$ as obtained using the CFD modeling approach in Fig. 9. A first observation that can be drawn is the fact that ERGO’s design 2 (in yellow in the top plots) has a smaller area than the default design (in blue), while this constraint was enforced (Eq. (24)). This is once more a consequence of the MAR approach which has added this unfeasible design to the design set with a prediction that outperformed the other evaluated designs. A possible approach to avoid this is driven up the contribution of the model uncertainty in the MAR prediction. However, experience has shown that this destabilizes the surrogate training set, as denoted above.

Fig. 9

Optimal designs visualization and characteristics prediction at Ma =0.734 and Re =6.5×106 as predicted by ERGO (a and b) and ERGO-II (c and d). (Color version online.)

View large Download slide

Optimal designs visualization and characteristics prediction at Ma $= 0.734$ and Re $= 6.5 \times 10^{6}$ as predicted by ERGO (a and b) and ERGO-II (c and d). (Color version online.)

On the other hand, the designs presented by ERGO-II meet the area constraint. It can be noted that the optimal designs have a more rounded pressure side near the leading edge and rounded suction side near the trailing edge. Inspection of the characteristics indicates a higher stall angle for the optimal designs, but a lower value of $C_{L}$ -gradient in the linear region in comparison with the baseline design. Second, a lower drag and moment coefficient is observed.

In regards to the computational cost, ERGO and ERGO-II are built on the EGO algorithm, a form of BO as discussed in the introduction, which is known to be a data-efficient method and therefore used specifically for methods where the computational cost of a function evaluation is far more substantial than the cost of training the GPs and finding the next infill. Both methods required roughly the same computational budget, corresponding to 250 function calls, of which 150 succeeded. This added up to approximately 4–5 h for each optimization. If an alternative surrogate-free approach had been attempted to solve the problem, the computational cost of the uncertainty quantification would have prohibited solving the problem due to the computational cost of the CFD calls; for example, performing a single UP of a design by means of a sampling approach or a spectral method such as the Taylor-based method or Polynomial Chaos would surmount to the same computational cost as the whole ERGO-II optimization.

The successful outcome of the validation case using ERGO-II proves the effectiveness of the improved scheme.

7 Conclusion and Outlook

In this paper improvements upon ERGO, a surrogate-assisted optimization-under-uncertainty methodology for robust design optimization, was presented. This was realized in four steps:

The hypervolume expected improvement was introduced in the framework as an alternative to the Euclidean EI. The Euclidean EI exhibited sensitivity to variations in objective magnitudes and was restricted to only two objectives. By incorporating the hypervolume EI, these limitations were mitigated, enabling more robust and versatile design selection based on the acquisition function.
A more advanced approach for assessing predictive variance was introduced. The initial implementation relied on a FOSM technique to estimate the predictive uncertainty of the stochastic objective measures. However, this lower-order approach led to the underemphasis of certain regions in the objective space. To address this drawback, the adoption of a SOSM approach was proposed, enabling a more accurate assessment of the predictive variance within the acquisition function.
A classifier to model the feasible space was introduced. In the original ERGO scheme, the consideration of objective function failures was neglected. To address this issue, the utilization of a Gaussian process classifier was proposed, which, when combined with the acquisition function, minimizes the disruption to the convergence of the optimization process.
The treatment of stochastic parameters was examined. The original ERGO scheme had a limitation in that it only allowed for stochastic design variables. Consequently, the direct assessment of the impact of parameter variability on the objective was not possible. In order to address this limitation, a two-step acquisition function optimization approach was undertaken to specifically handle stochastic parameters.

The scheme was compared with the original ERGO approach using one of the ZDG test function for varying dimensionality and validated on the second ADODG problem, the drag minimization of a RAE2822 airfoil in a viscous transonic flow while being subjected to stochasticity in both operating conditions and design variables. The outcome illustrated the success of the scheme in obtaining a robust solution within an affordable computational budget.

Future research tracks are currently focused on two aspects of the proposed ERGO-II scheme. First, alternative approaches to the uncertainty propagation routine are explored. For highly non-linear objective spaces, the second-order approximation might leave wanting. Therefore, exact propagation methods are explored for specific covariance function and input distribution combinations. A second research track is focused on dealing with the limitations inherent to Bayesian optimization: namely the assumption of a smooth output space and a limited number of variables. Some recent advances in the machine learning community, namely manifold GPs, are currently being explored.

Acknowledgment

The author gratefully acknowledges the funding by the Research Foundation-Flanders (FWO), Belgium, through the Junior Post-Doctoral fellowship of Jolan Wauters. The computational resources (Stevin Supercomputer Infrastructure) and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by Ghent University, Belgium, FWO and the Flemish Government—department EWI, Belgium.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Nomenclature

Greek Symbols

$α$ =: angle of attack
$α$ =: GP predictive mean coefficients
$β$ =: polynomial coefficients
$θ$ =: covariance function hyperparameters
$κ$ =: hypervolume cell
$ν$ =: Matérn smoothness parameter
$σ^{2}$ =: process variance
$Φ$ =: input space covariance matrix
$Ψ$ =: Gaussian process covariance matrix
$c$ =: chord length
$d$ =: dimensionality
$b$ =: Bernstein coefficients
$y$ =: vector of objective function evaluations
$F$ =: model matrix
$P$ =: pareto front
$S$ =: training set
$x$ =: deterministic/mean design vector
$x$ =: stochastic design vector
$y$ =: Gaussian process of the objective
$X$ =: design set
$X$ =: design space
$C_{D}$ =: drag coefficient
$C_{L}$ =: lift coefficient
$C_{M}$ =: moment coefficient
Ma =: Mach number
Re =: Reynolds number

Subscripts

$p$ =: pressure side
$s$ =: suction side
$h v$ =: hypervolume
$E u$ =: Euclidean
$*$ =: optimal value

Superscripts

$M$ =: matricization
$V$ =: vectorization
$c$ =: classifier
$i$ =: interpolator
$T$ =: Taylor approximation
$⊤$ =: transpose

Operators

$A (\cdot)$ =: activation function
$ϕ (\cdot)$ =: standard normal probability density function
$Φ (\cdot)$ =: standard normal cumulative distribution function
$ψ (\cdot, \cdot)$ =: covariance function
$E_{{\cdot}} [\cdot]$ =: expected value operator
$μ_{{\cdot}} (\cdot)$ =: predictive mean
$f (\cdot)$ =: trend vector function
$H [\cdot]$ =: hypervolume measure
$K_{ν} (\cdot)$ =: Bessel function
$L (\cdot)$ =: likelihood function
$N (\cdot, \cdot)$ =: normal distribution
$Σ_{{\cdot}} (\cdot)$ =: predictive variance
$V_{{\cdot}} [\cdot]$ =: variance operator
$y (\cdot)$ =: objective vector function
$y (\cdot)$ =: Gaussian process

Abbreviations

ADODG =: aero design optimization discussion group
BO =: Bayesian optimization
DoE =: design of experiments
E(R)GO =: efficient (robust) global optimization
((S)R)EI =: ((safe) robust) expected improvement
EP =: expectation propagation
(F/S)OFM =: (first/second)-order first moment
(F/S)OSM =: (first/second)-order second moment
GA =: genetic algorithm
GPI/C =: Gaussian process interpolator/classifier
KL =: Kullback–Leibler
LHS =: Latin hypercube sampling
MAR =: missing-at-random
MC(I/S) =: Monte Carlo (integration/sampling)
MLE =: maximum likelihood estimation
NSGA-II =: non-dominated sorting genetic algorithm II
OUU =: optimization-under-uncertainty
PDF =: probability density function
RANS =: Reynolds-averaged Navier–Stokes
RDO =: robust design optimization
SA =: surrogate-assisted
(S/M)OO =: (single/multi) objective optimization
SQP =: sequential quadratic programming
U(Q/P) =: uncertainty (quantification/propagation)

Appendix: Bayesian Optimization

This appendix presents a concise introduction to GPI and BO, extended with a number of metrics essential to the ERGO framework, making use of the notation introduced in Sec. 2. For more information on GPI and BO, the reader is referred to Refs. [8,27].

BO is motivated by the problem of finding a global minimizer (or the mathematically equivalent maximizer) of an unknown objective function

y

in a data-efficient manner.

x_{*} = \arg min_{\begin{matrix} x \in X \end{matrix}} y (x)

(A1)

In this context,

X

denotes the designated space of interest, typically a confined subset of

R^{d (X)}

⁠, with

d (X)

typically taken smaller than 20 [16]. It is further assumed that the function

y : X \mapsto R

is a continuous black-box function. Evaluating

y

is an expensive process devoid of noise or gradients, making it impossible to directly obtain the gradient without resorting to approximations. The core concept of BO is to construct a surrogate model of the objective function. This surrogate model can be updated and queried to guide the optimization decisions [8].

The GPI is the surrogate model commonly used in BO. From a function-space perspective, a Gaussian process $y$ can be defined as a distribution over functions, where the values of $y (x)$ at any set of points ${x_{i} | i = 1, \dots, n}$ follow a joint Gaussian distribution. A Gaussian process is fully characterized by its second-order statistics $y (x) = G P (m (x), ψ (x, x^{'} | θ))$ ⁠. Here, $m (x) : X \mapsto R$ represents the mean function, and $ψ (x, x^{'} | θ) : X \times X \mapsto R$ represents the covariance function fully described by its parameters $θ$ ⁠. Different interpolators can be obtained based on the formulation of the mean function. When a constant value is used, the method is known as “simple Kriging.” Moreover, if a multivariate polynomial $f (x)$ is employed, the interpolator is referred to as “universal Kriging.”

In a fully Bayesian approach, a prior distribution is placed over the hyperparameters. However, customary a point estimate of the hyperparameters is used in the shape of the MLE:

\arg max_{θ} L^{i} (θ | S^{i}) ≜ p (S^{i} | θ)

with

S^{i} ≜ {X, y}

the evaluated data. The likelihood function can be simplified by taking the logarithm and marginalizing out the hyerparameters

\begin{aligned} \log L^{i} (θ | S^{i}) & = \frac{n (S^{i})}{2} \log 2 π + \frac{n (S^{i})}{2} σ^{2} + \frac{1}{2} \log | Ψ | \\ + \frac{1}{2 σ^{2}} (y - F β)^{⊤} Ψ^{- 1} (y - F β) \end{aligned}

(A2)

with

F

the model matrix defined by

F^{(i, j)} = f_{i} (x_{j})

⁠,

Ψ

the correlation matrix defined by

Ψ^{(i, j)} = ψ (x_{i}, x_{j})

⁠,

β

the coefficients of the multivariate polynomial trend, and

σ^{2}

the process variance. This can be further simplified by enforcing the stationarity with respect to

β

and

σ^{2}

⁠, which corresponds to a generalized least squares approach such that

β = (F^{⊤} Ψ^{- 1} F)^{- 1} F^{⊤} Ψ^{- 1} y

(A3a)

σ^{2} = N^{- 1} (y - F β)^{⊤} Ψ^{- 1} (y - F β)

(A3b)

When leaving out the constant terms the concentrated log marginal likelihood function is obtained which can be directly optimized. Note that the resulting function on the one hand quantifies the fit of the surrogate through

σ^{2}

⁠, while on the other hand, it penalized complexity through

| Ψ |

⁠.

Building forth on the definition of a Gaussian process and the theorem of Bayes the predictive distribution

p [y (x) | S^{i}]

can be directly evaluated and gives again a normal distribution of which the mean

E [y (x | S^{i})] ≜ μ_{y} (x)

and variance

V [y (x | S^{i})] ≜ Σ_{y} (x)

can be directly evaluated

μ (x) = β^{⊤} \cdot f (x) + α^{⊤} \cdot ψ (x)

(A4a)

Σ (x) = σ^{2} {ψ (0) - | ψ (x) ‖_{Ψ^{- 1}}^{2} + | g (x) ‖_{Γ}^{2}}

(A4b)

with

ψ (x) = [ψ (x_{1}, x), \dots, ψ (x_{n}, x)]^{⊤}

α = Ψ^{- 1} (y - F β)

⁠,

Γ = (F^{⊤} Ψ^{- 1} F)^{- 1}

⁠,

g (x) = F^{⊤} Ψ^{- 1} ψ (x) - f (x)

⁠.

The availability of a closed formulation of the predictive mean and variance permits the derivation of the Jacobian of the predictive mean

\nabla μ (x)

and the corresponding variance

Σ_{\nabla μ} (x)

[39–42]

μ_{\nabla y}^{} (x) = β^{⊤} \cdot \nabla f (x) + α^{⊤} \cdot \nabla ψ (x)

(A5a)

Σ_{\nabla y}^{} (x) = σ^{2} {Δ ψ (0) - ‖ \nabla ψ (x) ‖_{Ψ^{- 1}}^{2} + ‖ \nabla g (x) ‖_{Γ}^{2}}

(A5b)

The formulation is extended even further upon for the Hessian of the predictive mean

μ_{Δ y}^{} (x)

and the corresponding predictive variance of the Hessian

Σ_{Δ y}^{} (x)

μ_{Δ y}^{} (x) = β^{⊤} \cdot Δ f (x) + α^{⊤} \cdot Δ ψ (x)

(A6a)

Σ_{Δ y}^{} (x) = σ^{2} {D^{(4)} ψ (0) - ‖ Δ ψ (x) ‖_{Ψ^{- 1}}^{2} + ‖ Δ g (x) ‖_{Γ}^{2}}

(A6b)

Additionally, the covariance between the mean prediction, gradient, and hessian can be derived

\begin{aligned} Σ_{y, \nabla y}^{} (x) & = σ^{2} {ψ (0) - ψ (x)^{⊤} \cdot Ψ^{- 1} \cdot \nabla ψ (x) \\ + g (x)^{⊤} \cdot Γ \cdot \nabla g (x)} \end{aligned}

(A7a)

\begin{aligned} Σ_{y, Δ y}^{} (x) & = σ^{2} {Δ ψ (0) - ψ (x)^{⊤} \cdot Ψ^{- 1} \cdot \nabla ψ (x) \\ + g (x)^{⊤} \cdot Γ \cdot \nabla g (x)} \end{aligned}

(A7b)

\begin{aligned} Σ_{\nabla y, δ y}^{} (x) & = σ^{2} {D^{(4)} ψ (0) - ψ (x)^{⊤} \cdot Ψ^{- 1} \cdot \nabla ψ (x) \\ + \nabla g (x)^{⊤} \cdot Γ \cdot Δ g (x)} \end{aligned}

(A7c)

In this paper, exclusive use is made of the Matérn covariance function, common in the aerospace community [27,43]. It is often preferred over the conventional exponential covariance function, which infinitely differentiable nature leads to non-physical levels of smoothness [44].

ψ (x, x^{'} | θ, ν) = \frac{2^{1 - ν}}{Γ (ν)} {(\sqrt{2 ν} ‖ x - x^{'} ‖_{θ})}^{ν} K_{ν} (\sqrt{2 ν} ‖ x - x^{'} ‖_{θ})

(A8)

with

Γ (ν)

the Gamma function and

K_{ν}

the modified Bessel function of the second kind

Γ (ν) = \int_{0}^{\infty} t^{ν - 1} e^{- t} d t

(A9a)

K_{ν} (x) = \frac{π}{2} \frac{I_{- ν} (x) - I_{ν} (x)}{\sin (ν π)}

(A9b)

I_{ν} (x) = \sum_{m = 0}^{\infty} \frac{1}{m! Γ (m + ν + 1)} {(\frac{x}{2})}^{2 m + ν}

(A9c)

In the Matérn function

θ

and

ν

respectively correspond to the width parameters and the smoothness parameter. A strongly simplified expression can be obtained when

ν = 1 / 2 + p

with

p \in N

is used. The result is a

p

-times mean square differentiable function. In this work, exclusive use is made of the case where

p = 2

\begin{aligned} ψ (x, x^{'} | θ, 2) & = (1 + \sqrt{5} ‖ x - x^{'} ‖_{θ} + \frac{5}{3} ‖ x - x^{'} ‖_{θ}^{2}) \\ \cdot \exp (- \sqrt{5} ‖ x - x^{'} ‖_{θ}) \end{aligned}

(A10)

Using the chain rule the Jacobian

\nabla_{x} ψ (x, x^{'})

and Hessian

δ_{x x} ψ (x, x^{'})

can be obtained with the former being a

d (X) \times 1

vector and the latter a

d (X) \times d (X)

matrix given by

\begin{aligned} \nabla_{x} ψ (x, x^{'} | θ, 2) & = - \frac{5}{3} \cdot θ \circ (x - x^{'}) \cdot (1 + \sqrt{5} ‖ x - x^{'} ‖_{θ}) \\ \cdot \exp (- \sqrt{5} ‖ x - x^{'} ‖_{θ}) \end{aligned}

(A11a)

\begin{aligned} δ_{x x} ψ (x, x^{'} | θ, 2) & = - \frac{5}{3} \cdot [d i a g (θ) (1 + \sqrt{5} ‖ x - x^{'} ‖_{θ}) \\ + 5 \cdot (θ \circ (x - x^{'}))^{\otimes}] \cdot \exp (- \sqrt{5} ‖ x - x^{'} ‖_{θ}) \end{aligned}

(A11b)

The Gaussian process interpolator is built using the object-oriented design and analysis of computer experiments (ooDACE) toolbox, which is an open-source matlab® software [45]. The construction process involves employing a multi-start sequential quadratic programming approach to optimize the concentrated log marginal likelihood function [46,47].

Now the enhancement in performance of a predicted point concerning the current best-evaluated point, denoted as

y_{min} =

min [y (X)]

⁠, can be formally expressed as follows:

I (x | y_{min}) =

max [y_{min} - y (x), 0]

⁠. By considering the stochastic nature of the interpolator, it becomes possible to assess the uncertainty in the prediction. This assessment is used to define the expected improvement

E_{y} [I (x | S^{i})]

[9]. It is important to note that the first term aims to minimize the objective (exploitation), while the second term aims to minimize the uncertainty in the prediction (exploration). This infill criterion serves as the foundation of the well-known EGO algorithm (Algorithm 3) [10].

\begin{aligned} EI (x | y_{min}) & ≜ E_{y} [I (x | y_{min})] \\ = (y_{min} - μ (x)) \cdot Φ (y_{min} | μ (x), Σ (x)) \\ + Σ (x) \cdot ϕ (y_{min} | μ (x), Σ (x)) \end{aligned}

(A12)

Notice in Algorithm 3: line 6 that

θ_{*}

is added as an argument. This substitution is made to indicate that an MLE approach is adopted instead of a fully Bayesian approach to evaluate the interpolator. Furthermore,

k_{E I}

is introduced as an additional constant that serves as a stopping criterion: if no further significant improvement can be made, stop (⁠

b r e a k

⁠) the optimization. Otherwise, the state of practice is to ensure the training set does not contain more than 500 designs [16].

EGO [10]

Algorithm 3

Require Evaluated sampling plan $S^{i} = {X, y}$

1: While computational budget is not exhausted do

2: $θ_{*} \leftarrow \arg max_{θ} L^{i} (θ | S^{i})$ ▷Eq.(A2)

3: $y_{m i n} \leftarrow min y$

4: $z_{*} \leftarrow max_{x} EI (x | θ_{*}, y_{m i n})$

5: If $z_{*} > k_{E I}$ then

6: $x_{*} = \arg max_{x} EI (x | θ_{*}, y_{m i n})$ ▷Eq.(A12)

7: $y_{*} \leftarrow y (x_{*})$

8: $S^{i} \leftarrow S^{i} \cup {x_{*}, y_{*}}$

9: else

10: $break$

11: $x_{*} \leftarrow \arg min_{x} y$

References

1.

Chatterjee

,

T.

,

Chakraborty

,

S.

, and

Chowdhury

,

R.

,

2019

, “

A Critical Review of Surrogate Assisted Robust Design Optimization

,”

Arch. Comput. Methods Eng.

,

26

(

1

), pp.

245

–

274

.

Google Scholar

Crossref

2.

Parnianifard

,

A.

,

Azfanizam

,

A.

,

Ariffin

,

M. K. A. M.

, and

Ismail

,

M. I. S.

,

2018

, “

An Overview on Robust Design Hybrid Metamodeling: Advanced Methodology in Process Optimization Under Uncertainty

,”

Int. J. Ind. Eng. Comput.

,

9

(

1

), pp.

1

–

32

.

3.

Balesdent

,

M.

,

Brevault

,

L.

,

Morio

,

J.

, and

Chocat

,

R.

,

2020

,

Overview of Problem Formulations and Optimization Algorithms in the Presence of Uncertainty

,

Springer International Publishing

,

Cham

, pp.

147

–

183

.

4.

Zang

,

T. A.

,

Hemsch

,

J.

,

Michael

,

H.

,

Luckring

,

J. M.

,

Peiman

,

M.

,

Padula

,

S. L.

, and

Stroud

,

W. J.

,

2002

, “Needs and Opportunities for Uncertainty-Based Multidisciplinary Design Methods for Aerospace Vehicles,” NASA, Report™-2002-211462.

5.

Persson

,

J. A.

, and

Ölvander

,

J.

,

2017

, “

How to Compare Performance of Robust Design Optimization Algorithms, Including a Novel Method

,”

Artificial Intell. Eng. Design, Anal. Manuf.

,

31

(

3

), pp.

286

–

297

.

Google Scholar

Crossref

6.

Kanno

,

Y.

,

2020

, “

On Three Concepts in Robust Design Optimization: Absolute Robustness, Relative Robustness, and Less Variance

,”

Structural Multidisciplinary Optim.

,

62

(

2

), pp.

979

–

1000

.

Google Scholar

Crossref

7.

Ryan

,

K. M.

,

Lewis

,

M. J.

, and

Yu

,

K. H.

,

2015

, “

Comparison of Robust Optimization Methods Applied to Hypersonic Vehicle Design

,”

J. Aircraft

,

52

(

5

), pp.

1510

–

1523

.

Google Scholar

Crossref

8.

Shahriari

,

B.

,

Swersky

,

K.

,

Wang

,

Z.

,

Adams

,

R. P.

, and

Freitas

,

N. D.

,

2016

, “

Taking the Human Out of the Loop: A Review of Bayesian Optimization

,”

Proc. IEEE

,

104

(

1

), pp.

148

–

175

.

Google Scholar

Crossref

9.

Mockus

,

J.

,

Tiesis

,

V.

, and

Zilinskas

,

A.

,

1978

, “The Application of Bayesian Methods for Seeking the Extremum,”

Towards Global Optimization 2: Proceedings of a Workshop at the University of Cagliari

, Vol.

2

,

L. C. W.

Dixon

, and

G. P.

Szego

, eds.,

Italy

, pp.

117

–

129

.

10.

Jones

,

D. R.

,

Schonlau

,

M.

, and

Welch

,

W. J.

,

1998

, “

Efficient Global Optimization of Expensive Black-Box Functions

,”

J. Global Optim.

,

13

(

4

), pp.

455

–

492

.

Google Scholar

Crossref

11.

Keane

,

A. J.

,

2012

, “

Cokriging for Robust Design Optimization

,”

AIAA J.

,

50

(

11

), pp.

2351

–

2364

.

Google Scholar

Crossref

12.

Keane

,

A. J.

, and

Voutchkov

,

I. I.

,

2020

, “

Robust Design Optimization Using Surrogate Models

,”

J. Comput. Design Eng.

,

7

(

1

), pp.

44

–

55

.

Google Scholar

Crossref

13.

Cheng

,

J.

,

Liu

,

Z.

,

Wu

,

Z.

,

Li

,

X.

, and

Tan

,

J.

,

2015

, “

Robust Optimization of Structural Dynamic Characteristics Based on Adaptive Kriging Model and CNSGA

,”

Struct. Multidiscipl. Optim.

,

51

(

2

), pp.

423

–

437

.

Google Scholar

Crossref

14.

Ribaud

,

M.

,

Blanchet-Scalliet

,

C.

,

Helbert

,

C.

, and

Gillot

,

F.

,

2020

, “

Robust Optimization: A Kriging-Based Multi-objective Optimization Approach

,”

Reliab. Eng. Syst. Saf.

,

200

(

1

), p.

106913

.

Google Scholar

Crossref

15.

Wauters

,

J.

,

2022

, “

ERGO: A New Robust Design Optimization Technique Combining Multi-objective Bayesian Optimization With Analytical Uncertainty Quantification

,”

ASME J. Mech. Des.

,

144

(

3

), p.

031702

.

Google Scholar

Crossref

16.

Keane

,

A. J.

,

2006

, “

Statistical Improvement Criteria for Use in Multiobjective Design Optimization

,”

AIAA J.

,

44

(

4

), pp.

879

–

891

.

Google Scholar

Crossref

17.

Wagner

,

T.

,

Emmerich

,

M.

,

Deutz

,

A.

, and

Ponweiser

,

W.

,

2010

, “On Expected-Improvement Criteria for Model-Based Multi-objective Optimization,”

Parallel Problem Solving From Nature, PPSN XI

,

Springer Berlin Heidelberg

,

Berlin, Germany

, pp.

718

–

727

.

18.

Couckuyt

,

I.

,

Deschrijver

,

D.

, and

Dhaene

,

T.

,

2014

, “

Fast Calculation of Multiobjective Probability of Improvement and Expected Improvement Criteria for Pareto Optimization

,”

J. Global Optim.

,

60

(

3

), pp.

575

–

594

.

Google Scholar

Crossref

19.

Haldar

,

A.

, and

Mahadevan

,

S.

,

2000

,

Probability, Reliability, and Statistical Methods in Engineering Design

,

Wiley

,

New York

.

20.

Oberkampf

,

W. L.

,

Trucano

,

T. G.

, and

Hirsch

,

C.

,

2004

, “

Verification, Validation, and Predictive Capability in Computational Engineering and Physics

,”

Appl. Mech. Rev.

,

57

(

5

), pp.

345

–

384

.

Google Scholar

Crossref

21.

Wauters

,

J.

,

Degroote

,

J.

,

Couckuyt

,

I.

, and

Crevecoeur

,

G.

,

2022

, “

SAMURAI: A New Technique for Efficient Optimization Under Uncertainty

,”

AIAA J.

,

60

(

11

), pp.

6133

–

6156

.

Google Scholar

Crossref

22.

Forrester

,

A.

,

Sóbester

,

A.

, and

Keane

,

A.

,

2006

, “

Optimization With Missing Data

,”

Proc. R. Soc. A: Math. Phys. Eng. Sci.

,

462

(

2067

), p.

935

.

Google Scholar

Crossref

23.

Bachoc

,

F.

,

Helbert

,

C.

, and

Picheny

,

V.

,

2020

, “

Gaussian Process Optimization With Failures: Classification and Convergence Proof

,”

J. Global Optim.

,

78

(

3

), pp.

483

–

506

.

Google Scholar

Crossref

24.

Nadarajah

,

S.

,

2013

, “Aerodynamic Design Optimization: Drag Minimization of the RAE 2822 in Transonic Viscous Flow,” McGill Computational Aerodynamics Group, Report ADODG Case 2.

25.

Girard

,

A.

,

2004

, “Approximate Methods for Propagation of Uncertainty with Gaussian Process Models,” Thesis.

26.

Xue

,

J.

,

Leung

,

Y.

, and

Ma

,

J.-H.

,

2015

, “

High-Order Taylor Series Expansion Methods for Error Propagation in Geographic Information Systems

,”

J. Geograph. Syst.

,

17

(

2

), pp.

187

–

206

.

Google Scholar

Crossref

27.

Rasmussen

,

C. E.

, and

Williams

,

C. K. I.

,

2006

,

Gaussian Processes for Machine Learning

,

MIT Press

,

Cambridge, MA

.

28.

Zitzler

,

E.

,

Deb

,

K.

, and

Thiele

,

F.

,

2000

, “

Comparison of Multiobjective Evolutionary Algorithms: Empirical Results

,”

Evolutionary Comput.

,

8

(

2

), pp.

173

–

195

.

Google Scholar

Crossref

29.

Fonseca

,

C. M.

, and

Fleming

,

P. J.

,

1996

, “On the Performance Assessment and Comparison of Stochastic Multiobjective Optimizers,”

Parallel Problem Solving From Nature – PPSN IV

,

Springer Berlin Heidelberg

,

Berlin, Germany

, pp.

584

–

593

.

Google Scholar

Crossref

30.

McKay

,

M. D.

,

Beckman

,

R. J.

, and

Conover

,

W. J.

,

1979

, “

A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output From a Computer Code

,”

Technometrics

,

21

(

2

), pp.

239

–

245

.

31.

Morris

,

M. D.

, and

Mitchell

,

T. J.

,

1995

, “

Exploratory Designs for Computational Experiments

,”

J. Statist. Planning Inference

,

43

(

3

), pp.

381

–

402

.

Google Scholar

Crossref

32.

Nadarajah

,

S.

,

2017

, “Aerodynamic Design Optimization: Drag Minimization of the NACA 0012 in Transonic Inviscid Flow,” McGill Computational Aerodynamics Group, Report ADODG Case 1.

33.

LeDoux

,

S. T.

,

Vassberg

,

J. C.

,

Young

,

D. P.

,

Fugal

,

S.

,

Kamenetskiy

,

D.

,

Huffman

,

W. P.

,

Melvin

,

R. G.

, and

Smith

,

M. F.

,

2015

, “

Study Based on the AIAA Aerodynamic Design Optimization Discussion Group Test Cases

,”

AIAA J.

,

53

(

7

), pp.

1910

–

1935

.

Google Scholar

Crossref

34.

Drela

,

M.

,

1998

, “Pros & Cons of Airfoil Optimization,”

Frontiers of Computational Fluid Dynamics 1998

,

World Scientific

,

Singapore

, pp.

363

–

381

.

Google Scholar

Crossref

35.

Kulfan

,

B. M.

, and

Bussoletti

,

J. E.

,

2006

, “

‘Fundamental’ Parametric Geometry Representations for Aircraft Component Shapes

,”

11th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference

,

Portsmouth, VA

,

Sept. 6–8

.

36.

Forrester

,

A.

,

Sóbester

,

A.

, and

Keane

,

A.

,

2008

,

Engineering Design Via Surrogate Modelling: A Practical Guide

,

Wiley

,

Hoboken, NJ

.

Google Scholar

Crossref

37.

Menter

,

F. R.

,

1994

, “

Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications

,”

AIAA J.

,

32

(

8

), pp.

1598

–

1605

.

Google Scholar

Crossref

38.

Wauters

,

J.

,

Degroote

,

J.

, and

Vierendeels

,

J.

,

2019

, “

Comparative Study of Transition Models for High-Angle-of-Attack Behavior

,”

AIAA J.

,

57

(

6

), pp.

2356

–

2371

.

Google Scholar

Crossref

39.

Han

,

Z.-H.

,

Görtz

,

S.

, and

Zimmermann

,

R.

,

2013

, “

Improving Variable-Fidelity Surrogate Modeling Via Gradient-Enhanced Kriging and a Generalized Hybrid Bridge Function

,”

Aerospace Sci. Technol.

,

25

(

1

), pp.

177

–

189

.

Google Scholar

Crossref

40.

Han

,

Z.-H.

,

Zhang

,

Y.

,

Song

,

C.-X.

, and

Zhang

,

K.-S.

,

2017

, “

Weighted Gradient-Enhanced Kriging for High-Dimensional Surrogate Modeling and Design Optimization

,”

AIAA J.

,

55

(

12

), pp.

4330

–

4346

.

Google Scholar

Crossref

41.

McHutchon, A., 2014, “Nonlinear Modelling and Control using Gaussian Processes,” Ph.D. thesis, University of Cambridge, https://mlg.eng.cam.ac.uk/pub/pdf/Mch14.pdf

42.

Lophaven

,

S. N.

,

Nielsen

,

H. B.

, and

Sondergaaard

,

J.

,

2002

, “DACE: A Matlab Kriging Toolbox,” Technical University of Denmark, Report IMM-TR-2002-12.

43.

Palar

,

P. S.

,

Zuhal

,

L. R.

,

Chugh

,

T.

, and

Rahat

,

A.

,

2020

, “On the Impact of Covariance Functions in Multi-Objective Bayesian Optimization for Engineering Design,”

AIAA SciTech Forum

,

Orlando, FL

,

Jan. 6–10

,

American Institute of Aeronautics and Astronautics

, pp.

1

–

13

.

44.

Stein

,

M. L.

,

1991

, “

A Kernel Approximation to the Kriging Predictor of a Spatial Process

,”

Ann. Inst. Statist. Math.

,

43

(

1

), pp.

61

–

75

.

Google Scholar

Crossref

45.

Couckuyt

,

I.

,

Dhaene

,

T.

, and

Demeester

,

P.

,

2014

, “

ooDACE Toolbox: A Flexible Object-Oriented Kriging Implementation

,”

J. Mach. Learn. Res.

,

15

(

1

), pp.

3183

–

3186

.

46.

Ugray

,

Z.

,

Lasdon

,

L.

,

Plummer

,

J.

,

Glover

,

F.

,

Kelly

,

J.

, and

Martí

,

R.

,

2007

, “

Scatter Search and Local NLP Solvers: A Multistart Framework for Global Optimization

,”

INFORMS J. Comput.

,

19

(

3

), pp.

328

–

340

.

Google Scholar

Crossref

47.

Nocedal

,

J.

, and

Wright

,

S. J.

,

2006

,

Numerical Optimization

,

Springer New York

,

New York

.

ERGO-II: An Improved Bayesian Optimization Technique for Robust Design With Multiple Objectives, Failed Evaluations, and Stochastic Parameters

Abstract

1 Introduction

2 Notation

2.1 Sub- and Superscripts.

2.2 Operators.

2.3 Fonts.

3 Efficient Robust Global Optimization

4 Improved Efficient Robust Global Optimization

4.1 Robust Hypervolume Expected Improvement.

4.2 Improved Epistemic Uncertainty Propagation.

4.3 Crash Constraint Treatment.

4.4 Stochastic Parameter Treatment.

4.5 Algorithmic Development.

5 Performance Assessment

5.1 Problem Description.

5.2 Methodology.

5.3 Results and Discussion.

6 Validation

6.1 Problem Description.

6.2 Parameterization.

6.3 Methodology.

6.3.1 Turbulence Modeling.

6.3.2 Solver-Specific Settings.

6.3.3 Mesh Settings.

6.3.4 Boundary Conditions.

6.3.5 Computational Solver.

6.4 Results and Discussion.

7 Conclusion and Outlook

Acknowledgment

Conflict of Interest

Data Availability Statement

Nomenclature

Greek Symbols

Subscripts

Superscripts

Operators

Abbreviations

Appendix: Bayesian Optimization

EGO [10]

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only