Adaptive Model Refinement with Batch Bayesian Sampling for Optimization of Bio-inspired Flow Tailoring

TL;DR

This paper enhances adaptive surrogate-based optimization by integrating Bayesian sampling to intelligently select sample locations, improving efficiency and effectiveness in complex design problems like bio-inspired flow control.

Contribution

It introduces AMR-PBS, a novel method combining adaptive model refinement with penalized batch Bayesian sampling for better sample placement in optimization.

Findings

AMR-PBS outperforms Bayesian EGO on benchmark functions.

The method achieves a 10% drag reduction in flow control design.

Effective sample placement improves optimization efficiency.

Abstract

This paper presents an advancement to an approach for model-independent surrogate-based optimization with adaptive batch sampling, known as Adaptive Model Refinement (AMR). While the original AMR method provides unique decisions with regards to "when" to sample and "how many" samples to add (to preserve the credibility of the optimization search process), it did not provide specific direction towards "where" to sample in the design variable space. This paper thus introduces the capability to identify optimum location to add new samples. The location of the infill points is decided by integrating a Gaussian Process-based criteria ("q-EI"), adopted from Bayesian optimization. The consideration of a penalization term to mitigate interaction among samples (in a batch) is crucial to effective integration of the q-EI criteria into AMR. The new AMR method, called AMR with Penalized Batch Bayesian Sampling (AMR-PBS) is tested on benchmark functions, demonstrating better performance compared to Bayesian EGO. In addition, it is successfully applied to design surface riblets for bio-inspired passive flow control (where high-fidelity samples are given by costly RANS CFD simulations), leading to a 10% drag reduction over the corresponding baseline (i.e., riblet-free aerodynamic surface).

Tables

Table 1. TABLE I: Problem settings for the benchmark problems; d 𝑑 d : problem dimension; N 0 subscript 𝑁 0 N_{0} : initial sample size; N f subscript 𝑁 𝑓 N_{f} : final sample size.

Benchmark Function	$d$	$N_0$	$N_f$	$N_{\text{run}}$	RAE; Median (Max)
					AMR-PBS	BEGO
Three-Hump Camel	2	20	30	20	4.39 (5.25)	23.47 (43.07)
Branin-Hoo	2	20	30	40	9.57 (30.63)	46.68 (646.07)
Hartmann6	6	60	90	20	0.16 (0.19)	0.30 (0.78)

Table 2. TABLE II: Results from Adaptive Model Refinement for 2 ∘ superscript 2 2^{\circ} Angle of Attack; N f subscript 𝑁 𝑓 N_{f} : final sample size; N refin subscript 𝑁 refin N_{\text{refin}} : number of refinements; h + superscript ℎ h^{+} : Height; s + superscript 𝑠 s^{+} : Spacing; σ 𝜎 \sigma : Distribution; RDC: Reduction in C D subscript 𝐶 𝐷 C_{D} compared to baseline.

SM	$N_0$	$N_{\text{added}}$	$N_{\text{refin}}$	$h^{+}$	$s^{+}$	$\sigma$	Cd-AMR	Cd-CFD	RAE [$\%$]	RDC [%]
Kriging	30	29	1	10.43	16.15	0.12	0.008017	0.008046	0.36	10.2
(Gaussian)

Equations

$$Y(x)=\hat{f}(x)+\epsilon$$

$$\begin{array}{rl}& H_0:{\mathbb{Q}}_{{\mathbb{P}}_{{}_{SM}}}(p_{cr})\ge {\mathbb{Q}}_{\Theta }(1-p_{cr})\end{array}$$

$$[{\mathbb{P}}_{({\mu }_{\epsilon },{\sigma }_{\epsilon })},{\Gamma }^{Infill}]=f_{{}_{PEMF}}(X,{\epsilon }^{\ast })$$

$$\mbox{$\Gamma^{Infill}$}=\Bigg{\{}\left.\begin{aligned} \tiny&\scalebox{0.99}{$\lceil\frac{\ln(\varepsilon_{mod}^{*})-\ln(a^{CURR})}{b^{CURR}}\rceil-N_{s}^{CURR}$}&\mbox{\small Type~{}1}\\ &\scalebox{0.99}{$\lceil\exp(\frac{\ln(\varepsilon_{mod}^{*})-\ln(a^{CURR})}{b^{CURR}})\rceil-N_{s}^{CURR}$}&\mbox{\small Type~{}2}\end{aligned}\right.$$

$${\mathcal{B}}_r(\mathbf{x},{\mathbf{x}}_i)=\{\mathbf{x}\in \mathcal{X}:\Vert {\mathbf{x}}_i-\mathbf{x}\Vert \le \rho \};{\mathbf{x}}_i\in X^{Infill}$$

$$\gamma (\mathbf{x},{\mathbf{x}}_i)=1-P(\mathbf{x}\in {\mathcal{B}}_r)$$

$$\gamma (\mathbf{x},{\mathbf{x}}_i)$$

$$=P(\mathcal{N}(0,1)\le \frac{L\Vert \mathbf{x}-{\mathbf{x}}_i\Vert -M+\hat{f}({\mathbf{x}}_i)}{{\sigma }_r({\mathbf{x}}_i)})$$

$$=\frac{1}{2}\text{erfc}\left(-\frac{L\Vert \mathbf{x}-{\mathbf{x}}_i\Vert -M+\hat{f}({\mathbf{x}}_i)}{\sqrt{2{\sigma }_r^2({\mathbf{x}}_i)}}\right)$$

$$\begin{array}{rl}& L_j^{ϝ}=\min \{x_j^{min},x_j^{min}\}\\ & U_j^{ϝ}=\max \{x_j^{max},x_j^{max}\}\end{array}$$

$$z(y)=h{e}^{-y^2/{\sigma }^2}$$

$$h^{+}=\frac{h{U}_{\tau }}{\nu }\sqrt{\frac{C_f}{2}}$$

$$s^{+}=\frac{s{U}_{\tau }}{\nu }\sqrt{\frac{C_f}{2}}$$

$$\underset{\mathbf{x}=h,s,\sigma }{\mathrm{minimize}}$$

$$g_1(\mathbf{x})=6\sigma -s\le 0$$

$$g_2(\mathbf{x})=s-6h\le 0$$

$$g_3(\mathbf{x})=\sigma -0.6h\le 0$$

$${\tau }_w=\mu {\frac{dU}{dy}}_{y=0}-\rho \overline{u^{\prime }{v}^{\prime }}$$

$$Y={\hat{y}}_{{}_{SM}}+\epsilon$$

$$\mbox{$\Delta f^{~{}t}_{k}$}=\Bigg{\{}\left.\begin{aligned} \small|\frac{f^{~{}t}_{k}-f^{~{}t-\tau}_{k}}{f^{~{}t}_{k}}|~{}~{}&~{}\mbox{if~{}$f^{~{}t}_{k}\neq 0$}\quad\\ |f^{~{}t}_{k}-f^{~{}t-\tau}_{k}|~{}~{}&~{}\mbox{if~{}$f^{~{}t}_{k}=0$}\quad\\ {\rm where}~{}k=1&,2,3,...,N_{pop}\end{aligned}\right.$$

$$p_{cr}=Pr(\epsilon \le \eta )={\int }_0^{\eta }\mathbb{P}({\epsilon }^{{}^{‎}\prime })d{\epsilon }^{{}^{‎}\prime }$$

$$\begin{array}{rl}& H_0:{\mathbb{Q}}_{{\mathbb{P}}_{{}_{SM}}}(p_{cr})\ge {\mathbb{Q}}_{\Theta }(1-p_{cr})\end{array}$$

$${\mathbb{Q}}_{\Psi }(p)=inf\{x\in \mathbb{R}:p\le {\Psi }_{(c.d.f.)}(x)\}$$

$${\mathbb{Q}}_{{\mathbb{P}}_{{}_{\mathrm{S}{\mathrm{M}}^{\mathrm{CURR}}}}}\ge {\mathbb{Q}}_{\Theta }^{t=t^{\ast }}$$

$${\epsilon }_{mod}^{\ast }=(1-\frac{{\mathbb{Q}}_{\Theta }^{t=t^{\ast }}-{\mathbb{Q}}_{\Theta }^{t=t^{\ast }-\tau }}{{\mathbb{Q}}_{\Theta }^{t=t^{\ast }}})\times {\epsilon }_{mod}^{CURR}$$

$$I=0.16R{e}^{-1/8}$$

$$k=1.5(UI{)}^2$$

$$\omega =0.9^{-1/4}\sqrt{k/l}$$

Full text

I Introduction

Optimizing complex systems often involves computationally expensive simulations (e.g., CFD) to evaluate system behavior and estimate quantities of interest. While surrogate models and surrogate based optimization [1, 2, 3] provides a tractable alternative, in their native form they tend to compromise on the fidelity of the optimization process. Specifically, low-fidelity surrogate models often mislead the search process during optimization, leading to sub-optimal or even infeasible solutions. Variable-fidelity optimization approaches seek to address these issues and offer attractive trade-offs between computational efficiency and fidelity of the optimal solutions obtained, that is, the ability to quickly arrive at optimal solutions that can be relied upon. In doing so, these methods must answer the following questions depending on the interplay of model uncertainty and function improvement over iterations – when (during optimization) and where (in the design space) to re-sample in order to refine the model on the fly? There exist very few frameworks that allow flexible choice of a variety of surrogate models (i.e., model independent) and optimization processes, while also answering both these questions. In this paper, we present important fundamental extension of one such method, known as Adaptive Model Refinement [4, 5]; the new advancements extend the capability of the original method by providing a mechanism to decide the optimal location of batches of new samples, allowing effective utilization of resources. This multifidelity optimization framework is first tested on benchmark functions and then applied to design ribleted 3D airfoil surfaces for bio-inspired flow tailoring. The remaining portion of this introduction section provides a brief review of variable fidelity optimization in general, batch sampling methods in specific, and converges on the objectives of this paper.

I-A Variable Fidelity Optimization

Major surrogate modeling methods that have been used in surrogate based optimization include polynomial response surfaces [6], Kriging [7, 8], moving least square [9, 10], radial basis functions (RBF) [11], support vector regression (SVR) [12], artificial neural networks [13] and hybrid surrogate models [14]. In Variable fidelity optimization approaches, model management strategies, e.g., model selection, switching and/or refinement, adaptively integrate models of different levels of fidelity and (computational) cost into the optimization process. Different model management approaches have been reported in the literature for integrating low-fidelity models within optimization algorithms. One class of model management strategies is developed based on the Trust-Region concept [15, 16, 17, 18, 19]. The Trust-Region methods seek the agreement of the function and its gradient values in the low-fidelity model with those estimated in the high-fidelity model. However, these techniques may not be directly applicable in problems where gradients are expensive to evaluate, or where zero-order algorithms are being used for optimization.

Another class of insitu model management strategies developed for non-physics-based low-fidelity models (aka. surrogate models) involves adding infill points where additional evaluations of the high-fidelity model or experiment are desired to be performed. In these approaches, infill points are generally added in (i) the region where the optimum is located (local exploitation), and/or (ii) the entire design space in order to improve the global accuracy of the surrogate (global exploration)[20, 21, 22]. Infill points can be added in a fully sequential (one-at-a-time), or a batch sequential manner. Various criteria exist for determining the locations of the infill points including (i) index-based criteria (e.g., (integrated- and maximum) Mean Square Error (MSE) and maximum entropy criteria) and (ii) distance-based criteria (e.g., Euclidean distance, Mahalonobis distance, and Weighted distance criteria)  [23, 24, 25, 26, 27, 28, 29]. Jones et al. in 1998 [23] developed a well-known model management strategy that is based on an expected improvement (EI) criterion, and is called efficient global optimization (EGO), essentially a type of Bayesian Optimization approach. Many variations of this method exist in the literature [30, 31, 32, 33], detailed review of which is outside the scope of this paper. The EGO method and its variations, however, are generally applicable only to surrogate based optimization based on Gaussian process models.

Jin et al. [34] reviewed different one-at-a-time sequential sampling criteria, and illustrated their clear benefits over single-stage methods [35, 36]. Loeppky et al. [37] explored different batch and one-at-a-time sequential criteria for the Gaussian process model, and pointed out potential cost advantages of batch sampling approaches over one-at-a-time augmentation. One major practical advantage of batch sampling is the ability to better exploit parallel computing resources. A more fundamental advantage of batch sampling is tied with the flexible refinement needs to ensure a reliable optimization progress, which may not necessarily be met with a single sample addition. However batch sampling also brings with it new challenges with regards to computational tractability of modeling the potential interactions among the samples in the batch. A brief review of batch sampling methods specific to Bayesian Optimization, converging on the concept that we build upon, is given next.

I-B Batch Sampling in Optimization

Given the growing popularity of Bayesian optimization (BO) in solving high-dimensional optimization problems and for hyper-parameter tuning of machine learning architectures, and benefits of being able to add samples in a batch as opposed to the one-at-a-time approach in classical BO, a number of batch Bayesian optimization approaches have emerged in recent years. In the work by Ginsbourger et al. [38], a multi-point expected improvement (qEI) is generalized to the batch setting, by using the qEI acquisition function for batches of $q$ points. However, it is challenging to identify the points that jointly maximize qEI since the computational cost of evaluating the function and its derivative scales poorly with increasing $q$. Later efforts have suggested heuristics for approximating the qEI computation [39, 40, 41]. One such method is the Constant Liar approach  [42], which is a sequential batch building method, that iteratively adds a point that maximizes the single point acquisition function. This approach assumes that evaluating this point will return a particular constant ‘lie’ value, which will temporarily augment the model training set with synthetic values and refitting the GP.

In a recent work, González et al. [43] proposed an alternative batching method by approximating the repulsive effect when batching. Under a Gaussian Process prior, target values of nearby points in sample space are expected to be highly correlated. Thus, when choosing a batch of samples, we may wish for the batch members to be sufficiently far apart to maximize the information gained. To do this, the authors propose the Local Penalization method that sequentially assembles batches of samples by successively penalizing the acquisition function around points previously selected, using a penalization radius based on the estimated Lipschitz constant of the acquisition function surface. We extend this method for usage in determining optimum locations (i.e., with qEI) of the new batch of samples during refinement events invoked by the AMR method [4], where the batch size is given by AMR request to meet the needs of balancing model uncertainty with fitness improvement rate. In doing so, we integrate a lower-level BO search instance (for the purpose of sampling) within the top level AMR-based optimization, where the latter typically uses a population based algorithm such as Particle Swarm Optimization [44]. This integration is done such that the main advantages of AMR, w.r.t. surrogate model-independence, exploitation of powerful (population-based) global optimizers, and ability to decide when to refine (important for computational efficiency), are retained.

I-C Objectives of this Paper

The primary objectives and associated contributions of this paper are given below:

Integrate, into the Adaptive Model Refinement method, a request-adaptive batch sampling approach that maximizes the EI while mitigating interactions among samples (within the batch) in a tractable manner. The central contribution lies in the ability (thus enabled) to maximally exploit infill sampling resources while preserving the original benefits (e.g., model independence) of the AMR framework. 2. 2.

Test the upgraded AMR method, called AMR with Penalized Batch-Bayesian Sampling or AMR-PBS, on different benchmark functions and compare its performance with that of the Bayesian Efficient Global Optimization (BEGO) method. 3. 3.

Apply the new AMR-PBS method to perform multi-fidelity optimization of bio-inspired surface riblets on 3D airfoil sections, to achieve maximum drag reduction. This objective lends itself to the practical contribution of this paper – namely provisions evidence of the effectiveness of the upgraded AMR method is solving complex multi-fidelity optimization problems.

The remainder of the paper is structured as follows. The next section provides an overview of the Adaptive Model Refinement method, followed by a detailed description of the new batch sampling strategy in AMR. In Section III, the performance of AMR-PBS is analyzed and compared with a Bayesian EGO algorithm, using a set of test benchmark problems. For studying the effectiveness of the proposed algorithm, a CFD-based problem which is computationally very expensive to simulate is considered. Section IV-C describes the automated CFD framework and its optimization formulation, and discusses the results. Finally, concluding remarks are given in Section V.

II Variable-Fidelity Optimization with Adaptive Model Refinement (AMR)

In this section, we introduce Adaptive Model Refinement (AMR), which is a decision-making tool for the timing of surrogate model refinement and identifying the optimal batch size and location of infill points in surrogate-based optimization (SBO). Performing model refinement (i.e., increasing model fidelity) too early in the design process can be computationally expensive while wasting resources to explore undesirable regions of the design domain. On the other hand, updating the surrogate model too late might mislead the search process early on to sub-optimal regions of the design domain, that is, leading to scenarios where the global optimum is outside of the region spanned by the population of candidate solutions in later iterations. The implementation of the proposed AMR in a population-based algorithm involves the five major steps, which are illustrated within the AMR flowchart shown in Fig. 1. These five steps can be summarized as:

Step 1

This step occurs once. First, a set of initial sampling points are generated in the design space using a Latin hypercube (LH) sampling [35]. The best surrogate model is then constructed using the Concurrent Surrogate Model Selection (COSMOS) framework [45, 46], which automatically selects the best-suited model (from a pool comprising RBF, Kriging, SVR, and ANN models with different underlying kernel choices). The fidelity of the chosen surrogate model is given by PEMF [47], which provides a measure of the predicted error distribution. Finally, the initial population is generated at $t=1$, using the initial surrogate model.

Step 2

At every iteration ($t$) of the heuristic optimization algorithm, the current surrogate model is used to update the function values of the population, and then set $t=t+1$. In this paper, particle swarm optimization is the chosen optimization algorithm.

Step 3

The stopping criterion is checked. The following two different methods can be used as the stopping criteria: (i) the difference between optimum values of five consecutive iterations is less than a threshold value, and (ii) the maximum allowed number of evaluations of function is reached. In this paper, the optimization algorithm stops when the relative change in the fitness function values in five consecutive iterations is less than a predefined function tolerance, $\delta_{F}$. If the termination criterion is satisfied, the current optimum (the best global solution in case of PSO) is identified as the final optimum and the optimization process is terminated. Otherwise, go to Step 4.

Step 4

The AMR metric is evaluated in this step. If the AMR metric is satisfied, go to Step 5. Otherwise go to Step 2.

Step 5

In this step, a model refinement occurs, and the surrogate model is updated through a process that involves five sub-steps (see the expansion of “Update Surrogate Model” step in fig. 1), of which the batch sampling is the important new addition. Go to Step 2

In practice, the AMR metric (Step 4) is not necessary to be applied at every iteration; the user can specify that it be applied after every $\tau$ iteration. The first version of the AMR is presented in our previous work [5, 4], hence in this work we only focus on a new extension of the AMR. The details of Step 1 to Step 3 can be found in Appendix A. In the next subsection, we briefly talk about Step 4 and then expand on the batch sampling aspects of Step 5.

II-A The Adaptive Model Refinement (AMR) Metric

In AMR, it is assumed that the uncertainty associated with surrogate model can be evaluated in the form of an error distribution, $\mathbb{P}_{i}$. Under this assumption, the fitness function values evaluated using the current surrogate model (in the SBO process) can be related to the corresponding high-fidelity model estimation or experiment-based system evaluations as

[TABLE]

Here, $\hat{f}(.)$ and $\varepsilon$, respectively, represent the response of the current surrogate model and the stochastic error associated with it, and $Y$ is the corresponding high-fidelity model (or experiment-based system) response.

The AMR metric is then defined based on ”whether the uncertainty associated with a surrogate model response is higher than the observed improvement in the relative fitness function of the population”. The proposed model refinement metric is designed to use the stochastic global measures of surrogate model error and the distribution of solution improvement. This metric is formulated as a statistical test for surrogate model (SM) as follows

[TABLE]

where $p_{cr}$ is the critical probability ($0<p_{cr}<1$) and $\mathbb{Q}$ represents a quantile function of a distribution.

For the sake of illustration, assume the surrogate model error distribution ($\mathbb{P}_{{}_{SM}}$) and the distribution of fitness function improvement ($\Theta$)) follow a log-normal distribution, and $p_{cr}=p^{*}$. In this case, the null hypothesis will be rejected, and the optimization process will use the current surrogate model if $\mathbb{Q}_{\Theta}~{}>~{}\mathbb{Q}_{\mathbb{P}_{{}_{SM}}}$, as illustrated in Fig. 2(a). Conversely, if $\mathbb{Q}_{\Theta}~{}<~{}\mathbb{Q}_{\mathbb{P}_{{}_{SM}}}$, the null hypothesis will be accepted, the surrogate model will be updated in the optimization process, as shown in Fig. 2(b).

II-B Batch Sampling

As mentioned above, the Adaptive Model Refinement (AMR) technique is a decision-making tool for the timing of surrogate model refinement, and for identifying not only the optimal batch size, but also the location of infill points needed to refine the surrogate model during a population based optimization process (e.g., genetic algorithm or PSO based process).

II-B1 Batch Size

The strategy of tracking the variation of error with an increasing density of training points in the PEMF method [47] is then applied to determine the batch size for the samples to be added to the current sample points (

$X^{CURR}$

). The inputs and outputs of PEMF (in the AMR method) can be expressed as:

[TABLE]

where the vector $\rm X$ represents the sample data (input and output) used for training the surrogate model; and $\varepsilon^{*}$ is the desired fidelity in the model refinement process. In Eq. 3, $\mathbb{P}_{(\mu_{\varepsilon},\sigma_{\varepsilon})}$, and $\Gamma^{Infill}$, respectively, represent the distribution of the error in the surrogate model, and the batch size for the infill points to achieve a desired level of fidelity in the AMR method. To this end, the desired batch size (

$\Gamma^{Infill}$

) is estimated using the inverse of regression functions used to represent the variation of error with sample density (the Type 1 and Type 2 prediction functions used in PEMF), as shown below:

[TABLE]

where

$a^{CURR}$

and

$b^{CURR}$

are regression coefficients of the VESD function in the current surrogate model that are determined using the least square method, and

$N_{s}^{CURR}$

represents the size of the current sample points.

II-B2 Batch Location

The location of the new infill points (

$\rm X^{Infill}|N(X^{Infill})=\Gamma^{Infill}$

) in the input space is determined based on a GP-based multi-points expected improvement criterion (q-EI). For a batch sampling, it is necessary to account for (and in our case mitigate) the interaction between the batch of future samples. Direct joint optimization of future candidate samples (in the batch) [48] is expensive to compute. Recently, Gonzalez et al. [43] reported a computationally tractable approximation to model the interactions, using a local penalization term. Based on this penalization factor, we adopt a criterion based on Expected Improvement with Local Penalization (EI-LP). The implementation of this approach is summarized in Algorithm 1.

In the given algorithm, we are calculating the penalty factor, $\gamma(x;x_{i},\mathcal{L})$ (line 12 of Algorithm 1), which enables local exclusion zones based on the Lipschitz properties of the actual function ($f(x)$). This term tends to smoothly reduce the acquisition function in the neighborhood of the existing samples. To compute this penalty, we define a ball $\mathcal{B}_{r}$ with radius $\rho$ around each subsequently chosen sample point:

[TABLE]

The local penalty associated with a point $\mathbf{x}$ is defined as the probability that $\mathbf{x}$ does not belong to the ball $\mathcal{B}_{r}$, as given by:

[TABLE]

By assuming that the distribution of the ball radius $\rho$ is Gaussian with mean $(M-\hat{f}(\mathbf{x}_{i}))/L$ and variance $\sigma_{r}^{2}(\mathbf{x}_{i})/L^{2}$, we can write the following expressions for the local penalty:

[TABLE]

Here, $M=\max_{\mathbf{x}}f(\mathbf{x})$ is the actual maximum value and $L$ is a valid Lipschitz constant ($\|f(\mathbf{x}_{1})-f(\mathbf{x}_{2})\|\leq L\|\mathbf{x}_{1}-\mathbf{x}_{2}\|$).

The lower bound ($L^{\digamma}_{j}$) and the upper bound ($U^{\digamma}_{j}$) of the $\rm j^{th}$ dimension of the current sampling range ($\digamma$) are given by

[TABLE]

where, $x_{j}^{min}$ and $x_{j}^{max}$ are, respectively, the lower and the upper bounds spanned by the current population of particles in the $\rm j^{th}$ dimension of the design space; $X_{j}^{min}$ and $X_{j}^{max}$ are the lower and upper bounds of the $\rm j^{th}$ dimension of the design space.

II-C Optimization Algorithm: Particle Swarm Optimization

The AMR can be integrated with most population-based optimization algorithms, such as evolutionary and swarm intelligence algorithms, where improvement can be tracked from one generation of the population to another. Particle swarm optimization (PSO) [49] is particularly well-suited to AMR. This is because identities of population members or search agents (known as particles) are preserved across generations, allowing ready-made tracking of solution quality improvement for each agent from one generation to another.

In this paper, we use an advanced implementation of the PSO algorithm called mixed-discrete PSO (MDPSO), developed by [44]. Unlike the conventional PSO algorithm, MDPSO provides: (i) an explicit diversity preservation capability that mitigates the possibility of premature stagnation of particles, and (ii) an ability to deal with both discrete and continuous design variables. MDPSO has been used to solve a wide variety of highly non-convex (often multimodal) mixed-integer nonlinear programming problems in wind farm design [50] and design of unmanned aerial vehicles [51]. Further description of the MDPSO algorithm can be found in the following paper [44].

III Benchmark Testing

III-A Testing Overview

To evaluate the proposed AMR algorithm, two types of experiments are conducted using three distinct test benchmark functions, namely:

1. the two-dimensional Three-Hump Camel function,
2. the two-dimensional Branin-Hoo function, and 3) the six-dimensional Hartmann6 function. The detailed descriptions of these functions can be found in [52]. In Experiment 1:, AMR is run to analyse its performance over the three test functions, and compare with that of Bayesian Efficient Global Optimization (BEGO) algorithm [33]. In Experiment 2: a batch size analysis is undertaken to explore the performance of AMR in the Hartmann6 test function, across multiple initial investment (batch size). The results of the experiments are evaluated and compared in terms of relative absolute error (RAE). The problem settings used (e.g., number of initial/total samples, number of runs per problem for robustness analysis, etc.) are given in Table I

III-B Experiment 1: A Comparative Analysis of AMR

Figure 3 shows the performance of the AMR and the BEGO algorithms, in terms of RAE. As shown in these results, the AMR outperforms the BEGO algorithm in all the test functions, both in terms of the median value and variance across multiple runs (with randomized LHS samples). The median value and the variation of the results obtained by AMR is superior to BEGO, particularly, for Branin-Hoo and Three-Hump Camel functions.

Figure 4 (a) shows the convergence history and the history of the surrogate model error over the optimization process. The contours for the actual function and the initial and final surrogate models in the AMR are illustrated in Fig. 4 (a). It can be seen from this plot that the AMR is able to reach to one of the optimums.

III-C Experiment 2: A Parameter Analysis of Initial Sample Size

In the second experiment, we study the impact of the initial investment/sample size ($N_{0}$) on the performance of the AMR technique. For this purpose, we use five different initial sample sizes varying from 60 to 180 with step size 30, and the maximum investment is set at 210 samples. For each case, we run the AMR approach 10 times, where the critical probability is set at 0.3 (i.e., $p_{cr}=0.3$) and $\tau$ is set at 2. Also, the population size ($N_{\text{pop}}$) and the maximum iteration ($\text{Iter}_{\text{max}}$) are set at 60 and 100, respectively. The obtained results are shown in Fig. 5. It can be seen from this figure that small initial sample size ($N_{0}=60$ in this problem) improves the chance to lead the optimizer to the actual optimum. For initial sample size larger than 30, the performance of the optimizer decreases (and for large values is almost similar, which can be because of the low complexity of the Hartmann6 problem). The reason for this observation is that when the initial investment ($N_{0}$) is small, the optimization is allowed to explore the design space and garner significant improvements before the refinement events get invoked.

IV Application of AMR-PBS: Optimization of Riblet-based Flow Tailoring

IV-A Riblet-based Flow Tailoring

Lulekar et al. [53] presented a new approach to represent and optimize surface riblets inspired by passive surface features observed in marine animals, for the purpose of reducing drag coefficients. While understanding of a broad range of riblet dimensions that could provide aerodynamic benefits existed, actual benefits might be realizable only in a very narrow range, which in the case of nature, is an outcome of millions of years of evolution. Hence, investigating (hypothesized) aerodynamic efficiency gains with minute surface riblets can become a challenging endeavor. The only ways to capture the complex flow mechanisms of interest induced by $O(10^{-1})$ millimeter scale riblets on a surface of a length scale, $O(10^{2})$ mm, are to conduct high-resolution computational fluid dynamic (CFD) simulations or physical (wind tunnel) experiments. The prohibitive time/cost of both approaches limits the potential number of trials needed to identify the right desired riblet geometry. In this paper, we are using a set of open-source tools for running high-fidelity CFD simulation, where each run needs approximately 16 hours of computational time, executed on the UB CCR academic clusters using 2 compute nodes with Intel Xeon Processor E5-2660 (25M Cache, 2.20GHz), 16 cores per node, and 48GB RAM, indicating the expense of one single CFD simulation. One solution for overcoming this computational cost is minimizing the number of high-fidelity evaluations. For this purpose, we use AMR, where the surrogate model is trained over CFD or physical experiments, can address this challenge, and effectively lead us to the very small range of aerodynamically-promising riblet geometry.

In order to exploit AMR, an automated modeling and design framework is needed that connects the various processes from design generation/CAD modeling to solution post processing. A first of its kind framework in the riblet design (for passive flow control) domain was developed by Lulekar et al. [53]. Abrief descriptin of the riblet geometry and the automated framework are given next.

IV-B The Geometry of Riblets

Different type of riblet shapes have been considered over decades to obtain drag reduction. We adopt our smooth riblet, which can be termed as patterns of Gaussian-shaped ridge lines [53]. Figure 6 shows the ridge-lines as a continuous extraction of a 2D Gaussian curve.

The geometry of the ridge is parameterized in terms of peak height ($h$), spacing ($s$) and the distribution ($\sigma$) of the curve. Thus, any point $z$, on the ridgeline can be expressed as

[TABLE]

where $-3\sigma<y<3\sigma$. As illustrated in Fig. 6, the pattern of ridgelines covers the entire top surface of the 3D airfoil section, and the cross-sectional shape/size of the ridge does not change in the streamwise direction in our current implementation. In most of the studies are conducted by changing the protrusion height ($h$) and the spacing ($s$) between the riblets to impede the streamwise vortices and their performance is measured based on a specific set of non-dimensional parameters, to account for the change in size of the flow structures like vortex diameter, where $h^{+}$ and $s^{+}$ is given as follows,

[TABLE]

[TABLE]

Further details of the physics and geometry of the riblets in the near wall flow, which affect the overall drag coefficient of the riblet surface can be found in our previous paper [53]. Here, instead, we focus on the modeling and optimization framework, which its central component is the presented AMR-based algorithm. A summary of this framework is described in the next.

IV-C Automated CFD Framework

In order to be able to study and solve the optimization problem to find the best drag reduction due to the riblets with the minimum human interaction, we developed an automated framework [53], which integrates disparate computational tools in batch mode. The overall framework is illustrated in Fig. 7 and it contains seven main blocks that can be divided into three categories: 1) Design-of-Experience (DoE) and metamodel training and selection (COSMOS) blocks; 2) CFD simulation blocks; and 3) multi-fidelity Optimization block. Each block is described in details in the next sections. Further descriptions of the CAD modeling, mesh generation, CFD flow solver, and post processing can be found in Appendix B.

IV-D Optimization Formulation and Design of Experiments

The characteristics of the boundary layer are greatly affected by the type of the riblet used, spacing between two riblets and the protrusion into the boundary layer. Hence, identifying the suitable riblet shape boils down to three parameters, namely height, spacing and the distribution of the Gaussian curve. Given the objective of reducing drag ($D$) on the 3D airfoil section, the optimization problem can be stated as:

[TABLE]

where $h\in[0.2,0.6]$, $s\in[0.72,3.6]$ and $\sigma\in[0.12,0.46]$, with the dimension of $h$ and $s$ are in millimeters. The first constraint, $g_{1}(\mathbf{x})$, is used to prevent adjacent Gaussian ridges from overlapping with each other. The second constraint, $g_{2}(\mathbf{x})$ in Eq. (12), that restricts the inter-ridge spacing to 6 times the height of the ridge. The third constraint, $g_{3}(\mathbf{x})$ in Eq 12, restricts the Gaussian curve to flatten out where the advantage of the riblets impeding the cross flow momentum is lost. Equations (10), (11), (12) can be used to set the bounds for $h$, $s$ and $\sigma$.

IV-E Results

The simulation and optimization is performed for a case study with angle of attack $2^{o}$. This particular angle of attack is chosen, to observe the effect of riblet on drag reduction at higher pressure gradients. For this problem, the following settings are used: initial number of sample points (initial investment) used for training the surrogates is $N_{0}=30$, a population size of PSO is set at $N_{\text{pop}}=30$, and the maximum iteration is set at $MaxIter=100$. In addition, the critical probablity in AMR is set at $p_{cr}=0.2$, and refinement events are allowed every 6th iteration during the optimization process. The convergence history of the AMR for this problem is shown in Fig. 8. For this problem, the AMR requested only one refinement event at iteration 12th with 29 infill points. Table II summarizes the results obtained for this problem. Kriging with Gaussian kernel is identified as the optimal surrogate model. The final results show that the optimal design variables (the chosen dimensions of the riblets) are $h^{+}=10.43$, $s^{+}=16.15$ and $\sigma=0.12$. The coefficient of drag (Cd) at the optimum is 0.008046, which the surrogate model estimated it with 0.36% accuracy. These dimensions provide a notable drag reduction (10.2%) in compare with barefoil (airfoil without riblets).

IV-F Flow Characteristics of the Optimum Design

Riblets are known to impede the cross-stream translation of the streamwise vortices. The flow very close to the wall explains the impending physics of the drag reduction. To understand the impending physics behind the reduction, we need to visualize the flow very close to the wall. There is not a particular cross-stream velocity present but due to tip effect of the airfoil, a small amount of cross-flow stream is induced as visualized in Fig. 9(b). Although a cross-flow is induced we can see a decrease in momentum transfer. It can be seen from Fig. 9(a) that the protrusion of the riblets from the surface into the boundary layer create zones of high-speed and low-speed fluid regions. The low-speed fluid region comes in contact with the riblet valleys which constitute a higher portion of the surface area whereas the high-speed fluid comes in contact only with the riblet peaks.

Figure 10(a) shows the boundary layer in the ridge valley and at the ridge peak at $x/c=0.3$, compared with the boundary layer profile for the barefoil. This figure clearly shows the presence of different velocity gradients, attributed to low speed and high speed fluid regions. Velocity gradients are quantified by measuring the shear stress, where the wall shear stress is given by,

[TABLE]

Riblets demarcate the fluid into low speed and high speed regions, where the shear stress is simultaneously redistributed. In this study, the riblets are aligned with the free-stream velocity and hence, a shift in the wall shear stress distribution can be seen in the spanwise direction. Figure 10(b) shows the distribution of the wall shear stress when compared to the barefoil.

V Concluding Remarks

In this paper, we proposed a fundamental extension of a variable-fidelity optimization approach, so-called Adaptive Model Refinement (AMR), by incorporating a multi-point expectation improvement mechanism to decide the optimal location of batches of new samples. For evaluation, we used three benchmark functions and compared with an EGO-based approach (Bayesian EGO). The results show that the AMR-PBS performs 2-8 times better than the Bayesian EGO with significantly higher robustness. A parameter analysis of initial sample size is performed, which shows that starting optimization with small sample size and utilizing resources once you need further accuracy can improve the optimization success.

In addition, we applied the proposed variable-fidelity optimization (AMR-PBS) to design ribleted 3D airfoil surfaces for bio-inspired tailoring. While the results were not necessarily representative of converged results, since optimization were stopped due to time restrictions, the optimal design do show a promising drag reduction of 10.2% in comparison with an airfoil without riblets. A large number of high and low speed regions are created due to the presence of the riblets balancing the adverse effects of increase in wetted surface area. Ideally, the riblets should be thin with sharp peak having minimal surface area to minimize the effect of high shear stress due to high-speed fluid contact. Although, a single high-resolution CFD simulation takes 7-16 hours, the given automated CFD framework, by utilizing the AMR technique and the parallel computing, was able to provision the optimal design in less than 40 hours. Future work entails investigation of how the interaction among samples in the batch impact the optimization process, and role played by the batch size in balancing potential interaction (and thus knowledge gain gap) and the need for reduction in model uncertainty for reliable continuation of the optimization process. Further application to an extended optimization of the passive flow surfaces to convergence will provide more comprehensive evidence regarding the effectiveness of the new AMR method.

Appendix A The Adaptive Model Refinement (AMR) Algorithm

The Adaptive Model Refinement (AMR) Metric

In this paper, it is assumed that the uncertainty associated with surrogate model can be evaluated in the form of an error distribution, $\mathbb{P}_{i}$. Under this assumption, the fitness function values evaluated using the current surrogate model (in the SBO process) can be related to the corresponding high-fidelity model estimation or experiment-based system evaluations as

[TABLE]

In Equation 14, $\widehat{y}_{SM}$ and $\varepsilon$, respectively, represent the response of the current surrogate model and the stochastic error associated with it, and $Y$ is the corresponding high-fidelity model (or experiment-based system) response. The relative improvement in the fitness function value ($\Delta f$) can be considered to follow an unknown distribution, $\Theta$, over the population of solutions. Here, $\Delta f$ in the $t^{th}$ iteration ($t\geq 2$) can be expressed as

[TABLE]

where, $\tau$ ($\in~{}\mathbb{Z}_{<t^{*}}$) is a user-defined parameter which regulates the occurrences of the “surrogate model refinement” in the proposed SBO approach. Numerical experiments exploring different values for the $\tau$ indicate that the $3\leq\tau\leq 5$ can be the suitable choice.

The AMR is then defined based on ”whether the uncertainty associated with a surrogate model response is higher than the observed improvement in the relative fitness function of the population”. The proposed model refinement metric is designed to use the stochastic global measures of surrogate model error and the distribution of solution improvement. Based on prior experience or practical design requirements, the designer is likely to be cognizant of what levels of global surrogate model error, $\eta$, is acceptable in an optimization process. Hence, $\eta$ can be perceived as a user preference. The critical probability, $p_{cr}$ for the current surrogate model (in the SBO process) with an error distribution $\mathbb{P}$ is then defined as the probability of the model error to be less than $\eta$. This definition can be expressed as

[TABLE]

The critical probability ($p_{cr}$) essentially indicates a critical bound in the error distribution $\mathbb{P}$ ($0\leq\varepsilon\leq\eta$). If the predefined cut-off value ($\beta$) of the $\Theta$ distribution lies inside this region, the current surrogate model is considered to be no more reliable for use in the optimization process.

The Adaptive Model Refinement (AMR) metric is formulated as a hypothesis for testing that is defined by a comparison between (a) the distribution of the relative fitness function improvement($\Theta$) over the entire population; and (b) the distribution of the error associated with the current surrogate model ($\mathbb{P}$) over the entire design space. This statistical test for one of the intermediate surrogate model (SM) in the SBO process can be stated as

[TABLE]

where $0<p_{cr}<1$ and $\mathbb{Q}$ represents a quantile function of a distribution; The $p$-quantile, for a given distribution function, $\Psi$, is defined as

[TABLE]

In Eq. 17, the critical probability, $p_{cr}$, is an Indicator of Conservativeness (IoC). The IoC is based on user preferences, and regulates the trade-off between optimal solution reliability and computational cost in SBO. Generally, the higher the IoC (closer to 1), the higher the solution reliability and the greater the computational cost; under these conditions, model switching events will occur early on in the optimization process.

For the sake of illustration, assume $\Theta$ and $\mathbb{P}_{{}_{SM}}$ follow a log-normal distribution, and $p_{cr}=p^{*}$. In this case, the null hypothesis will be rejected, and the optimization process will use the current surrogate model if $\mathbb{Q}_{\Theta}~{}>~{}\mathbb{Q}_{\mathbb{P}_{{}_{SM}}}$, as illustrated in Fig. 11(a). Conversely, if $\mathbb{Q}_{\Theta}~{}<~{}\mathbb{Q}_{\mathbb{P}_{{}_{SM}}}$, the null hypothesis will be accepted, the surrogate model will be updated in the optimization process, as shown in Fig. 11(b).

In this paper, kernel density estimation (KDE) is adopted to model the distribution of the relative improvement in the fitness function over consecutive ${\rm\Delta t}$ iterations. The uncertainty associated with surrogate models, and the batch size for the samples to be added in model refinement process of AMR are determined using an advanced surrogate error estimation method, called Predictive Estimation of Model Fidelity (PEMF) [54]. The PEMF method is derived from the hypothesis that “the accuracy of approximation models is related to the amount of data resources leveraged to train the model” [54]. A brief description of the KDE and PEMF methods are provided in the following sections.

Model Refinement Based on PEMF

In the previous section, the formulations and concepts of the AMR metric are defined. Based on the AMR metric, the model refinement will be performed at the $t^{*}$-th iteration of SBO under the condition

[TABLE]

where $\rm SM^{{}_{CURR}}$ represents the current surrogate model in the optimization process. Model refinement is performed to efficiently improve the fidelity of the current surrogate model to achieve the desired fidelity for the upcoming iterations of SBO. In this paper, the desired fidelity, $\varepsilon_{mod}^{*}$, is determined using the history of the fitness function improvement in the optimization process, which is given by:

[TABLE]

where, $\varepsilon_{mod}^{CURR}$ is the predicted modal error value associated with the current surrogate model.

Appendix B Automated CFD Framework

Geometry and Meshing

A finite span of NACA$0012$ airfoil is considered with riblets placed on the top and bottom surface of the airfoil. The airfoil contains a constant chord with an aspect ratio of 2. The chord length is 0.127m (5 inches) and is visualized in Fig. 6 along with the riblets. The riblets are a continuous extraction of a 2D Gaussian curve, where the riblet valleys and the peaks are aligned in the flow direction 6. Three-dimensional CAD model of symmetric NACA$0012$ airfoil is build using the $SALOME-8.3.0$ and imported as an unstructured triangulated surface (.stl) with a precision of $O(10^{-2})mm$. The surface mesh is imported for volumetric grid generation carried out using $SnappyHexMesh$. To promote a greater stability with reasonable mesh count, while capturing the boundary layer profile more accurately, hex cells are preferred. Ten inflation layers are used with a growth rate of 1.2 and the smallest cell height is determined from $y^{+}\in[5,10]$.

To reduce the computational cost of the analysis, a symmetry plane is used and hence all the results are visualized on a half wing. The mesh generated in the domain approximately consists of 8 million cells. Fig. 12(a) shows the prism layer around the NACA$0012$ airfoil and a sliced section of the mesh accurately capturing the Gaussian riblet curve can be seen in Fig. 12(b).

Boundary Conditions

The CFD simulations are carried out in a bounded box.The domain has a velocity inlet ($free-stream$) with a constant velocity input. At the outlet, $pressureoutlet$ for the pressure is provided and $zeroGradient$ for the velocities. At the top, bottom and side wall we have $slip$ condition for the velocities and $zeroGradient$ for the pressure. A $symmetryplane$ boundary condition is used to reduce the computational time, and $noSlip$ boundary conditions is applied on the airfoil. The flow being transitional region, a %3.2 of turbulent intensity is given at the inlet. The turbulent kinetic energy (TKE) and the omega ($\omega$) at the inlet are calculated using the following equations:

[TABLE]

[TABLE]

[TABLE]

Solver Settings

An open source tool, $OpenFOAM$, is used to perform the CFD simulations. A pressure-based solver is used with $k-omega$ SST [55] turbulence models for wall bounded wall flow. The incoming $Re=3.38x10^{5}$, is calculated with a chord as the characteristic length, which is in the transitional ranger from laminar to turbulent. The fluid considered here is air, and treated as incompressible, as the Mach number is less than $0.4$. The pressure and velocity are coupled using the SIMPLE scheme. For spatial discretization and to calculate gradient of velocity ($\nabla U$), linear Upwind scheme is used. As stated above, the whole process is automated to evaluate multiple designs, OpenFoam is used for CFD evaluations. Simulations are carried out with a termination criteria of $O(10^{-6})$. To shorten the computation time, MPI library is used to parallelize the whole computation, and the simulations are performed using the distributed computing cluster (CCR) at the University of Buffalo. The primary quantities of interest obtained and derived from the CFD simulations are the shear stress data and co-efficient of drag.