A Relative Adequacy Framework for Multi-model Management in Design Optimization

We consider the problem of selecting among different physics-based computational models of varying, and often-times not assessed, ﬁdelity for evaluating the objective and constraint functions in numerical design optimization. Typically, higher-ﬁdelity models are associated with higher computational cost. Therefore, it is desirable to employ them only when necessary. We introduce a relative adequacy framework that aims at determining whether lower-ﬁdelity models (that are typically associated with lower computational cost) can be used in certain areas of the design space as the latter is being explored during the optimization process. We implement our approach by means of a trust-region management framework that utilizes the mesh adaptive direct search derivative-free optimization algorithm. We demonstrate the link between feasibility and ﬁdelity and the key features of the proposed approach using two design optimization examples: a cantilever ﬂexible beam subject to high accelerations and an airfoil in transonic ﬂow conditions.


Introduction
Numerical engineering design optimization requires computational models to predict system behavior in large and multi-dimensional design spaces.But what constitutes an adequate model?How can we characterize it?How do we choose among different models during the optimization process depending on their fidelity level in different areas of the design space?According to Gross, fidelity is the degree to which a model or simulation reproduces the state and behavior of a real-world object, feature or condition; therefore, fidelity is a measure of the realism of a model or simulation [1].In the literature, the term "multi-fidelity" is used as an adjective for several terms: physics [2,3], surrogates [4], approximations [5], analysis [6], optimization [7], mapping [8] and multidisciplinary design optimization (MDO) [9].Therefore, the use of the term can differ substantially, which can lead to misconceptions and inappropriate methods for multi-model management.The objective of this work is to develop a framework for managing the use of physics-based models of varying fidelity regardless of their disciplinary context.The proposed methodology will support decisions related to which model(s), and at which computational cost, should be used during the design exploration of the optimization process.
Typically, when model fidelity increases, so does computational cost as argued in [10].However, this is not a law, and there can be exceptions.In this work, we may sometimes expect that a high-fidelity model has higher computational cost than a low-fidelity model, but our methodology does not depend on such assumptions.Existing multi-fidelity optimization (MFO) approaches attempt to calibrate lowfidelity models or replace low-fidelity analysis results using data from higher fidelity analyses [11].Moreover, the performance of multidisciplinary systems is obtained not only by conducting the analysis of each discipline but also by capturing and accounting for their interactions.Therefore, quantifying system interactions using high-fidelity (HF) models is a computationally intensive process, because all disciplinary models have to be repeatedly exercised at each iteration until they achieve multidisciplinary convergence.A possible remedy is to use low-fidelity (LF) models in lieu of HF ones as long as adequate accuracy is maintained.Furthermore, intermediate-fidelity models may be possible to obtain.But how should multiple models be managed to achieve conver-gent solutions for single and multidisciplinary analysis?Before we discuss how to manage multiple models, we have to define how to assign fidelity levels to the available models.
The term multi-fidelity models (MFM) is used in the literature to denote using more than one model to predict the same phenomenon.To the best of our knowledge, the most recent literature reviews on the MFO topic are the ones of Fernández-Godino et al. [12] and Peherstorfer et al. [10].Fernández-Godino et al. classified MFM into two categories: multi-fidelity surrogate models (MFSM) and multifidelity hierarchical models (MFHM).In both of these, the HF model is used as the reference to validate, calibrate, or replace lower-fidelity models or surrogate models.However, the adequacy level of the HF model may change at different points of the design space.Consequently, considering it as an absolute criterion for model selection may mislead the optimization process as the latter explores the design space.Peherstorfer et al. concluded that relative relationships among models may be more effective than the hierarchy assumed by existing MF methods, and recommended such an approach for future developments of MFO frameworks [10].
Peherstorfer et al. characterize MFM approaches as multi-fidelity management strategies that define how different models are employed during the optimization process and how outputs from different models may be combined [10].They also classify MFM methods according to adaptation, fusion, and filtering.Adaptation enhances the output of the LF model with information from the HF model while the computation proceeds and the surrogate model is updated at each iteration.Methods based on fusion evaluate LF models and HF models and then combine information.An example of fusion is the co-Kriging method [6].Filtering methods invoke the HF model following the evaluation of an LFM filter.That is, the HF model is used only if the LF model is inaccurate [6].The authors include MFHM in the adaptation and filtering management strategies.
Tables 1 and 2 list works that use MFHM, distinguishing which of them was used as LF model and which as HF model, in two different disciplines (structural mechanics and fluid mechanics).
Fernández-Godino et al. stated in their review that there is no clear relationship between the computational cost ratio of the used models to HF models and the computational cost ratio of the MFO process to the HF optimization process [12].Figure 1 depicts the cost ratio of the MFO process to the optimization process using high-fidelity reference models (RFO) vs. the cost ratio of a single analysis of the lowfidelity model (LFA) to a single analysis of the high-fidelity reference model (RFA) for 18 papers that perform optimization using MFM.It shows that employing multiple models Fig. 1: Cost ratio of the optimization process using MFM (MFO) to the optimization process using high (reference)fidelity models (RFO) vs. cost ratio of a single analysis of the low-fidelity model (LFA) to a single analysis of the high(reference)-fidelity model (RFA); reproduced from [12], which uses data reported in [2, 5, 25-35, 17, 36-39] with low (LFA/RFA) cost ratio in the optimization process does not necessarily yield a big reduction in (MFO/RFO) cost ratio.Moreover, these ratios may change for different applications or optimization problem formulations even if the models stay the same [12].The various types of models (theoretical, mathematical, and computational) and errors (observational, modeling, and discretization) involved in the verification and validation process are presented in [40].In our work, we refer to discrepancies among models as inadequacies and use the hierarchy shown in Figure 2. Different approximation models A i may be available for a physical phenomenon.These models can be solved using different numerical methods N i j .Finally, parameter values p i j k pertaining to the numerical method used define a particular model instantiation M i j k .
Model inadequacy represents a model's inability to predict the true behavior of a physical phenomenon [41].An adequacy level is assigned to each model instantiation based on absolute and relative measures.These measures require a reference model to be selected based on the highest computational cost among available models.We use the term absolute model inadequacy (AMI) to describe the difference between a model's prediction and the reference model's prediction.If more than two models are available, we use the term relative model inadequacy (RMI) to describe the difference in the output between two different models other than the reference model.
Model selection is then based on the capability to attain feasible solutions, maintaining adequate accuracy, and reducing computational cost.Oberkampf et al. discussed the relationship between prediction and validation [42].Our methodology adopts their point of view and applies it to our proposed relative adequacy framework (RAF).

Proposed Methodology
We will illustrate our approach by means of a running example.In this example, computational models are used to compute the dynamic free response of a cantilevered threedimensional beam under the effect of gravitational acceleration.The beam is fixed at one end as shown in Figure 3. Let us first consider modeling this problem using two approximation theories: a Timoshenko beam element with 6 degrees of freedom (DOF) and an absolute nodal coordinate formulation (ANCF) for a beam element with 12 DOF.The derived computational models are then exercised using the parameter values listed in Table 3, i.e., for varying gravitational acceleration.Initially, the beam is horizontal and at rest.Table 3: Flexible beam dimensions and material properties As can be seen in Figure 4, the free response of the beam's tip changes smoothly until g v = 40m/s 2 ; then the tip starts to suffer from high-frequency fluctuations.ANCF captures the high-frequency modes and large displacements successfully.However, the Timoshenko beam model fails to predict the response accurately throughout the parameter space.The inadequacy between the two models increases with increasing acceleration.Furthermore, Figure 5 shows how the Timoshenko beam model cannot converge to attainable solutions for gravity accelerations higher than 50m/s 2 .Consequently, if only the Timoshenko beam model is used, inaccurate predictions (if at all attainable) will misguide the optimization process.On the other hand, ANCF has high computational cost; however, it has similar fidelity with the Timoshenko beam model for predicting the dominant frequency within the acceleration range g v ∈ [−10, −20]m/s 2 as depicted in Figure 6.

Relative Adequacy Framework
Let us consider the general constrained optimization problem min For example, suppose that we have one objective function and one constraint.When a numerical method N i j is used with parameter values p i j k to evaluate a function f , the corresponding function value is denoted by f i j k .Given available models n m , a reference model is chosen and assigned the new index naught (M 0 ).New indices are also assigned to the remaining n m − 1 models.In the remainder of the text, we will use l ∈ {1, 2, ..(n m − 1)} and h ∈ {1, 2, ..(n m − 1)} with l = h to consider two models at a time.The number of possible absolute inadequacies is n a = n m − 1, while the number of possible relative inadequacies between the n m − 1 models is

2
. The AMI of model h is given by where ε f h and ε G h are model inadequacies of the objective and constraint functions, respectively.They are defined as with respect to the reference model M 0 .At this point, model inadequacies can be represented graphically in an inadequacy space.The number of dimensions in the inadequacy space depends on the number of functions, i.e., ε ε ε h ∈ R 1+U .To keep dimensionality at two, we define one aggregate constraint function that sums constraint violations (if any) [43].Consequently, the inadequacy vector is reduced to consisting of two components ε f h and ε G h .
In the example of Figure 2, if we choose M 1 1 1 as a reference model, we have 13 remaining models; thus, n a = 13 and n r = 78.After specifying a reference model, a relative adequacy matrix (RAM) ε ε ε of size ((n m − 1) × (n m − 1)) is built to summarize the relative inadequacies among remaining models.For illustration, we will show the relative inadequacy space and the RAM for model instantiations   5 shows the RAM, whose elements are given by Each element ε(h, l) in the RAM reflects a relative rank between pairs of models.For example, the positive value of ε(1, 3) indicates higher adequacy of model On the contrary, the negative value of ε (1,6) indicates lower adequacy of model M 1 relative to tation that depicts relative inadequacies for each element in the RAM.

Penalizing Relative Inadequacy
Let us assume that M l has lower computational cost than M h .Our objective is to minimize εhl so we can improve the fidelity of M l to approach the one of M h as much as possible while keeping the low computational cost of M l .(dB) Fig. 6: Frequency response of the beam's tip at selected gravitational acceleration values Table 5: An example of the relative adequate matrix A traditional penalty function φ h penalizes constraint violations in the objective function: where r c is the penalty term parameter and F h = ∑ U u=1 max(g h u , 0) 2 is the aggregate constraint function.We adopt this approach and modify the penalty function to take into account model's fidelity: where By utilizing Eq. ( 4), this function can be simplified to We then minimize this penalty function using the derivativefree optimization algorithm described in the next section.

Mesh Adaptive Direct Search
Gradients of functions evaluated by means of computational models may not exist or may not be prone to reliable approximations at a reasonable cost [44].Moreover, these models may be prone to severe numerical noise or fail to return a value.In these cases, we need to resort to derivativefree optimization algorithms.In this work we use the mesh adaptive direct search (MADS) algorithm of [45], which is supported by convergence properties.
MADS ensures global convergence to a solution satisfying local optimality conditions based on the Clarke calculus for nonsmooth functions [46].Each MADS iteration consists of two steps: the search and the poll.The search step allows functions evaluation at any finite number of trial points generated by evaluating the objective function on a scaled discretized space of variables, called the mesh M .The search favours exploration of the design space and allows users to implement any appropriate method, given their knowledge of the problem [47].The mesh is centred around the incumbent solution d k ∈ R n where the mesh size is parametrized by ∆ m .The poll step involves local exploration around d k .During this step, the distance between the trial point and the current incumbent solution is dictated by the poll size parameter ∆ p [48].If the search fails to find an improvement, the poll proposes trial points around the incumbent solution following polling conditions.The optimization process either converges when the mesh size parameter is smaller than a user-specified threshold or terminates when a pre-described function evaluation budget is exhausted [47].

Our Search Step Strategy
We use the search step to solve a surrogate optimization problem.In this problem, the SM are not built for the output of the problem functions, but for the models' inadequacies.The purpose is to mitigate the discrepancy among available models to approach the adequacy level of the reference model while favoring the utilization of low-cost models.A Latin Hypercube (LH) is used to populate a set of design points S k to obtain a set of observations and build a surrogate model εhl which can be used to predict the value of εhl (d k ) for a given d k ∈ S k .The difference between the predicted inadequacy in the current iteration and the previous one is ASME c ; CC-BY distribution license In this work, we use polynomial response surfaces (PRSs) to build the SM [49][50][51].The surrogate problem is defined as min where ∆ r hl (d k ) is the radius of a trust region in the inadequacy space where surrogate model predictions are assumed to be reliable.For example, consider that we have four model instantiations to solve the problem in Eq. ( 10 The trust region radius is updated as where The trust region method used here is the one reported in [52] as implemented by March and Willcox in [53].The positive parameters R 1 < R 2 < 1 and c 1 < 1, c 2 > 1 are chosen to define the ranges of the ratio value for which reduction or enlargement are necessary [54].In our implementation we have used the values R 1 = 0.01, R 2 = 0.9, c 1 = 0.5 and c 2 = 1.5.
In case multiple estimated inadequacies exist within the trust region, the least expensive model is selected.If all estimated inadequacies lie outside the trust region, then all models are evaluated to update the SM.The process of evaluating candidates using estimated model inadequacy within a trust region allows for an improvement of the solution with fewer evaluations.The modified MADS pseudo-algorithm using the RAF and the trust-region approach is provided in Algorithm 1. Figure 10 depicts a flow chart of the proposed algorithm.

Scaling Properties of our Algorithm
We investigate the scalability of our proposed algorithm to larger number of variables by solving the following test problem from [45] min To consider multiple models of varying fidelity and computational cost, we introduce what we call inadequacy terms to the objective and constraint functions as listed in Table 6.Using the reference model (i.e., for h = 0), there is a single optimal solution to the problem: every component of the vector d is equal to − √ 3 and the optimal objective function value is − √ 3n.We solved the problem using our RAF-MADS algorithm for 6 values of n, namely 5, 10, 20, 30, 40, and 50.For each value of n, we solved the problem 5 times using different initial guesses.We used a convergence criterion of ∆ p ≤ 10 −12 and a termination criterion of a maximum of 600n function evaluations.In the search step, we used a 5n point Latin hypercube sample at the initial iteration and a n/5 random search at other iterations (if necessary).
Figure 11 depicts how the average (over the 5 runs for each n) number of function evaluations increases with n.The plot shows an increase in the reduced cost ratio: the fact that RAF-MADS utilizes relative relationships among models leverages the reduction in the number of evaluations of the reference model M 0 .We fitted an exponential to the obtained data, obtaining the expression where n f eval is the number of function evaluations.We tested this equation at n = 70 and obtained a result associated with a 4% relative residual error (see Figure 11).

Numerical Examples
We demonstrate the application of our proposed methodology on constrained design optimization problems considering both computational structural dynamics (CSD) and computational fluid dynamics (CFD) disciplines.
The CSD design problem of a flexible beam is used to investigate the impact of framework parameters on the ASME c ; CC-BY distribution license Algorithm 1 The MADS optimization pseudo-algorithm using relative model adequacies in trust regions ∈ D k Disregard models that fail to return values Calculate model cost ratios Select the reference model that has the highest computational cost Generate the RAM using Eq. ( 4) Build a surrogate model for each inadequacy εhl (d k ) in the RAM Evaluate trust ratios and trust radiuses using Eq. ( 12) and Eq.(11) cost ratio and accuracy of the overall optimization process.The beam analysis and initial configuration are adopted from [55].The problem has a low-dimensional design space but involves non-smooth functions and solutions that exhibit numerical noise.A total of 9 optimization runs are conducted by choosing six different sizes of initial trust regions, and three different sampling sizes at the initial step (k = 0).For each run, we allow a maximum of 100n black box evaluations; the optimization process can obviously stop earlier if the MADS termination criterion is satisfied.
The CFD design problem is adopted from the works of Koziel and Leifsson [56][57][58][59][60].In these works, the authors proposed a hybrid approach of MFM that utilizes both MFHM and MFSM to solve a physics-based design optimization problem [56].Moreover, they tested their proposed approach for different design space ranges of the airfoil shape optimization problem.To investigate the scalability of our proposed framework, we reproduce their approach and illustrate the different aspects between our proposed framework and their MFM method.

Structural Dynamics of Flexible Cantilever
The dynamic response of the cantilever beam described in Table 7 can be computed using three different models: absolute nodal coordinate formulation (ANCF) (A 1 ), Euler-Bernoulli (A 2 ), and Timoshenko (A 3 ).The main difference Table 7: Flexible beam dimensions and material properties among the three models is the way the cross sectional rotations are calculated, which has a considerable effect on the simulation of deformation shapes under large loads.Euler-Bernoulli beam's cross sectional rotations are assumed to occur due to bending only, since the plane section remains normal to the longitudinal axis.However, in Timoshenko's theory they are assumed to be a sum of two contributions of bending and shear deformations.ANCF makes no assumptions on the sectional rotations which are represented by interpolation functions.This allows for the cross section to deform and change shape [61].It is expected that this type of description, together with a three dimensional continuum mechanics approach, leads to more adequate results [62].
For isotropic material the vertical displacement r z satisfies the following partial differential equation under the Euler-Bernoulli beam assumption where q(x) is the distributed load which satisfies the following condition under the Timoshenko beam assumption [63,64]: where γ is the shear correction factor, and λ g is the modulus of rigidity.In ANCF the shape function is obtained using polynomials that are cubic in x and linear in y and z, where x is a material coordinate along the beam axis and y and z are the two other perpendicular coordinates.The configuration of the beam element is determined by the position and slope vectors of two end nodes N p and N q as shown in Figure 12.Each node is defined by one vector for the position r and three vectors the slopes r x , r y and r z .Thus the element has where S is the element shape function and e is the vector of nodal coordinates I is the 3 × 3 identity matrix, and with the nondimensional coordinates and l the initial length of the beam element.For Euler-Bernoulli and Timoshenko beam elements the shape functions S 1 , S 2 , S 5 and S 6 are calculated using Eq. ( 20), and S 3 = S 4 = S 7 = S 8 = 0. Using the finite-element method leads to solving the equations of motion in discretized form where [M] and [K] are the global mass and stiffness matrices respectively, and Q is the global force vector.The mass and stiffness matrices can be found in [61,65,66].The kinetic energy is calculated by

Modeling Parameters
The beam response can be obtained by means of three models named here "Euler-Bernoulli beam," "Timoshenko beam," and "ANCF."The optimization process is conducted considering a vector of parameters p for each model, p = [ne ∆t], where ne is the number of meshing elements and ∆t is the time step.The beam is under the effect of gravitational acceleration g v = −30 m/s 2 .
"ANCF" is assumed to be the reference model M 0 ; "Euler-Bernoulli beam" and "Timoshenko beam" are used as models M 1 and M 2 , respectively, where the numerical integration method is Newmark.We use the validated models reported in [55,67,68].Parameter values for each model are shown in Figure 13.

Results
Figure 14 depicts three different feasible design regions obtained using M 0 , M 1 and M 2 .The high-cost model M 0 is used as a reference to improve both lower-fidelity models M 1 and M 2 while keeping the computational cost as low as possible.Computational cost ratios are summarized in Table 8.Table 9 summarizes the MFO/RFO and reduced cost ratios for nine optimization runs; the first six runs use different initial trust region radii and the last three runs use different initial sample sizes.We can conclude that when the initial Fig. 13: Available models for predicting beam response trust radius increases, the computational cost decreases, but with an increase in the error of the obtained solution (defined as the relative deviation from the best-known solution of the problem).Figure 15 depicts the convergence history of the first six runs.Figure 16 depicts the choice of models based on the multi-model management method at each iteration of the optimization algorithm for the first six optimization runs.
We observe that the higher-fidelity reference model is preferred in the early stages of the optimization process (where data are collected to improve the prediction capability of the lower-fidelity models), but with less frequency as the trust region is decreasing.We can also observe that the higher-fidelity reference model is preferred when the predictive capability of the lower-fidelity models deteriorates as new areas of the design space are explored as depicted in Figure 17. Figure 18 illustrates that there is no clear relationship between the size of the initial trust regions and the number of evaluations of M 0 in the first six runs.For the last three runs, we observe a decrease in the number of evaluations of M 0 and an increase in the number of evaluations of M 1 .This observation explains the linear increase in the reduced cost ratio, the reason being that the increase of the number of samples at the initial design step contributes to larger trust regions in the approximated model inadequacy as depicted in Figure 19.

Airfoil Shape Parameterization
We consider this problem to investigate the scalability of our proposed framework and illustrate its differences from existing MFO approaches.Koziel and Leifsson solved this Fig.14: The feasible design space obtained using three different models problem to study the relationship between the design problem dimensionality and the computational cost of the optimization process solved using their proposed algorithm [56].NACA-four digits has a limited control to parameterize the airfoil shape with only three design variables [56].However, shape parameterization techniques (e.g., Bèzier curves) can provide better control by invoking more design parameters (control points).

NACA and Bèzier Airfoils
In the NACA-four digits, the airfoil is defined by three design variables relative to the chord length, namely maximum ordinate of the mean camber m c , chordwise position of the maximum ordinate p c , and maximum thickness to chord ratio t c [69].
Bèzier curves of order n B are defined as where P B (i) are control points and b B is a row vector with m B components verying from zero to one.Control points location is selected to represent the airfoil shape.Koziel and Leifsson used two setups in their study: one with 6 and one with 8 control points (see Figure 20) [56].
The norm of the difference between the actual NACA0012 and the curve shapes generated by Eq. 25 is minimized to find the approximate Bèzier airfoil shapes as shown in Figure 20.Invoking more control points yields a better approximation to the actual airfoil shape.This improves the control of the airfoil shape during the optimization process.However, additional control points increase computational cost.ASME c ; CC-BY distribution license

Problem Formulation
The transonic design problem considers a free-stream Mach number Mach ∞ = 0.75 and an angle of attack α = 1 o .The objective is to maximize the lift coefficient C l subject to constraints on the wave drag coefficient C dw,max = 0.006 and the non-dimensional cross-sectional area A cs,min = 0.075.The initial designs in all cases are approximations of the NACA 2412 airfoil shape.These approximation shapes are obtained for each design case by minimizing the norm of the difference between NACA2412 and the airfoil shape obtained by parameterization.The initial designs are different in each design problem since the number of design variables differs.The design domain for each problem is defined by the lower and upper bounds d and d listed in Table 10.The Nondimensional coefficients of lift C l and wave drag C dw are calculated from The force coefficients C x and C z are calculated by integrating the pressure coefficient C p counter clockwise around the airfoil surface as where ds is the panel length on the surface of the airfoil and θ is the angle between the panel and the horizontal axis.The pressure coefficient is defined as where p is the pressure in the flow-field, p ∞ is the free-stream pressure, ρ ∞ is the free-stream density, and V ∞ is the freestream velocity.

Derivation of Available Models
Transonic flow phenomenon occurs when there are mixed subsonic and supersonic local flows over a body moving with a speed near the speed of sound.Typically, the supersonic region of the flow is terminated by a shock wave where the flow slows down to subsonic speed.Across the shock there is an abrupt increase in pressure, temperature, density, and entropy.In this work, we consider a twodimensional transonic flow past an airfoil section assuming a steady, inviscid, and adiabatic flow with no body forces.the aforementioned assumptions: the Euler equations and the transonic small disturbance equation (TSDE).The Euler equations are a set of coupled non-linear partial differential equations that represent the conservation of mass, momentum, and energy [56].The TSDE is an approximation of the Euler equations (thin airfoils at small angle of attack).

Solution Methods
The Euler equations require a timemarching solution method.In the steady case, it is assumed that time-marching proceeds until a steady-state solution is reached.In this work, the temporal discretization of the governing equations is accomplished using an implicit time-marching algorithm.For the flow spatial discretization, a second-order upwind scheme is utilized.Fluent [70] is used to solve this approximation model.The TSDE can be solved with a finite difference scheme on a rectangular grid.We also used the legacy code (TSFOIL) to solve the TSDE; however, the finite difference scheme which is embedded in the code is inaccessible [71], which makes the tool a blackbox and thus amenable to our framework since our derivative-free optimization algorithm is ideal for blackbox problems.
Modeling parameters To set up multimodels we use the grid independence studies that exist in [56].These studies were conducted for the initial airfoil configuration only; therefore, they may need to be repeated if the airfoil shape parameters are modified.Figure 21 depicts the available models that utilize modeling attributes used in the adopted independent studies; RMSE is the root-mean-square error used as the termination criterion of each model instantiation.
For the Euler model, the free-stream Mach number, static pressure, and angle of attack are prescribed at the farfield boundary.The solution domain boundaries are placed at chord lengths in front of the airfoil, 50 chord lengths behind it, and 25 chord lengths above and below it.The computational grids are of structured curvilinear body-fitted C-topology with elements clustering around the airfoil and growing in size with distance from the airfoil surface as depicted in Figure 22.The number of nodes in the discretized flow domain is denoted n nodes .For the TSDE model, the number of panels n panels represents the discretized elements of the airfoil surface.

Results
The three design optimization problems with different number of variables (denoted DP#1, DP#2, and DP#3) are solved using our proposed RAF algorithm where the problem formulation is the same, except that the number of design variables is different as listed in Table 10.While all initial designs approximate the same NACA2412 airfoil shape parameters, there is a difference among them due to different number of design variables.Consequently, there is a considerable variation in the lift and drag coefficients of the initial designs.Figure 23 depicts some similarities among the optimized airfoil shapes; however, there are differences at the tailing-edge, where additional control points provide more degree of flexibility to modify the airfoil shape.The flexible shape parameterization impacts the lift coefficient as the optimized airfoil in DP#3 has the highest lift compared to DP#1 and DP#2, as summarized in Table 11.
Model instantiations, M 1 , M 2 and M 3 , have low execution CPU time (about 1 to 3 seconds depending on the selected number of panels) and for M 0 , M 4 , M 5 and M 6 evaluation takes higher CPU time (about 10 to 30 minutes depending on the selected number of nodes and the threshold value of the RMSE).The total CPU time of evaluating the low cost models in each optimization run is approximately from half to one evaluation of the reference model.and DP#2 needs 247 iterations, giving an increase of 98 iterations with the additional 4 design variables.The number of iterations increases to 374 for DP#3, an increase of 127 for the additional 4 design variables.The overall computational cost rises accordingly.However, the reduced cost ratio of the optimization process increases by 2.4% from DP#1 to DP#2 and by 0.3% from DP#2 to DP#3. Figure 24 depicts a slope change which indicates that any further increase in the number of design variables may lead the reduced cost ratio to approach a plateau.A possible reason is that the number of the initial samples, which increases accordingly with the increase of the number of design variables, is part of the search step that contributes to reducing the total number of design iterations of the MFO compared to the RFO process.Further increase in the number of design variables may be accompanied by a significant increase in the cost ratio if the number of evaluations of the reference model approaches the number of design iterations of the RFO. Figure 25 depicts the increase in the number of model evaluations for each model instantiation with increasing Fig. 24: Change in number of model evaluations and process cost ratio for increasing number of design variables of evaluations they obtained along with ours.While curves have similar trends, our algorithm requires considerable less (approximately one fifth) total number of evaluations.It should be noted that we present this plot for information only, and not as a comparison.Koziel and Leifsson consider only one high-fidelity model and one surrogate model while our approach considers several numerical instantiations of physics-based models (even though surrogate models can be considered as well).
Figure 26 depicts the number of selections of each model instantiation during the optimization process of the three problems.We observe that the reference model is selected only at the beginning of the optimization process; thereafter more inexpensive models are selected.
Independent studies are conducted for the same boundary and initial conditions to tune modeling parameters that ensure numerical stability.Typically, these studies are conducted for the initial design guess; however, design variables are altered during the optimization process, which may change the boundary conditions of the analysis problem.This may require repeating the independent studies to ensure the adequacy of the selected modeling parameters.Our proposed framework assess these studies for the selected samples at the initial step and during the optimization process (if the approximated model inadequacies are mistrusted).The exhausted cost for that assessment is compensated as the possibility of misleading the search process is minimized which is reflected on reducing the total number of design iterations.On the other hand, that assessment contributes to reducing the number of evaluations of the reference model during the optimization process.The main difference between our methodology and other existing MFHM or MFSM approaches is that we utilize relative relationships among available models regardless of their expected fidelity to enhance

Concluding Summary
This work introduces a framework for managing the use of models (whether physics-or data-based) of varying fidelity regardless of their disciplinary context.The framework considers a set of available models and assigns an adequacy level to each model based on its absolute and relative error relative to the other models and a reference model, respectively.This requires a reference model which is selected based on the highest computational cost among available models.Model inadequacies are represented graphically in a two-dimensional inadequacy space.Relative model inadequacies are concatenated in a relative adequacy matrix.The elements of this matrix are then used in an optimization problem formulation where a penalty function is minimized using the mesh adaptive direct search (MADS) algorithm.The MADS search step is utilized to solve a surrogate of the original optimization problem within trust regions where polynomial response surfaces are built to predict model inadequacies.This enables multiple models to make considerably faster predictions and attain solutions when non-linear phenomena emerge.
We demonstrate the proposed framework by solving two design optimization problems that have different disciplinary structures.The first problem is the optimization of a cantilevered flexible beam using three models.The presented results of this numerical example showed a significant reduction in the computational cost of the optimization process.The second problem is the airfoil shape optimization that utilizes Bèzier airfoils to enhance its aerodynamic characteristics.To investigate the scalability of our proposed framework, the size of the design space is scaled up by invoking more control points to the Bèzier airfoil geometry.The results show that both reduced cost ratio and optimization computational cost increases linearly with increased number of design variables.
The initial cost ratio (at the initial design step) increases with increasing number of available models and/or design variables.However, the reduced cost ratio does not decrease sharply.A possible reason is that the increase in the number of initial model evaluations increase the trust in approximated model inadequacy.Moreover, the evaluation of models is conducted at the search step.This is reflected in the reduction in number of reference model evaluations during the optimization process and the reduction in total number of design iterations.Independent studies are conducted for the initial design guesses and may be repeated (if necessary) during the optimization process to ensure numerical stability.We found that utilizing independent studies can assist the RAF to alleviate the initial cost ratio when it is used to solve design optimization problems with larger number of variables.In this work, we demonstrated the proposed framework using physics based models; however, it can be applied to any types of models, including (data-based) surrogate models.
Ongoing and future research focuses on extending the relative adequacy framework for solving multidisciplinary design optimization problems.

Fig. 9 :
Fig. 9: Trust regions of the relative inadequacy space

[ 1 ]
Initialization (k = 0) Set maximum number of evaluations k max .Set convergence tolerance (typically machine precision) Set initial poll and mesh sizes Set initial trust radiuses ∆ r hl (d k ) = 0 and inadequacy distance ||Λ hl (d k )|| 1 = +inf Populate initial guesses using Latin Hypercube (LH) and store them in the cache matrix D k Evaluate φ k using Eq.(5) ∀d k ∈ D k Set φ min = φ k for d k that returns minimum values for f k and F k [2] Iteration [2.1] Search If k > 0, then Populate possible points using LH and store them in S k Update the RAM using built SM Compute the difference among model inadequacies of the current iteration and the previous one End

Fig. 10 :Fig. 11 :
Fig. 10: Flow chart of the algorithm integrating mesh adaptive direct search with the relative adequacy framework using trust region principles

) 3 . 1 . 1
Problem Formulation The beam kinetic energy (KE) T (d) is minimized considering the beam cross section area C and the beam length L as design variables where d = [L C] T .The maximum absolute displacement in the vertical direction of the beam tip for t ∈ [0, 1.2s] must be greater than or equal to 1.5m.The design domain is defined by the lower and upper bounds d and d of d, where d = [1 0.1] T and d = [4 1] T .min d∈R n T (d)

Fig. 17 :Fig. 18 :
Fig. 17: Cumulative number of model evaluations in the first six runs

Fig. 25 :
Fig.25: Total number of model evaluations for each of the three cases required by the MFSM algorithm of Koziel and Leifsson and by our RAF-MADS algorithm.The MFSM curve is reproduced to the best of our ability based on the information available in[56]

Fig. 26 :
Fig. 26: Selected models during the optimization problem

Table 2 :
Papers in fluid mechanics organized by level of fidelity of models used (AN : Analytical, EM : Empirical, LI : Linear, PF : Potential flow, EU : Euler, RANS : Reynoldsaveraged Navier Stokes)

Table 6 :
, respectively If φ k < φ min , then Update φ min = φ k and d min = d k Go to 3 End if Otherwise, Select the model inadequacy that lies inside the trust region adjacent to the model with the lowest cost ratio Enhance the model output using the selected model inadequacy If φk < φ min , then Update φ min = φk and d min = d k Poll Determine the polling directions and frame Update the RAM using the SM built in 2.1 Compute the difference between model inadequacies of the current iteration and the previous one Select the model inadequacy that lies inside the trust region adjacent to the model with the lowest cost ratio If all inadequacies lie outside their adjacent trust region, then go to 2.1 Evaluate φk at each design point in the polling set using the model adjacent to the selected inadequacy Enhance the model output using the selected model inadequacy If φk < φ min , then Update φ min = φk and d min = d k If no stopping criteria is met, then go to 2.2 Otherwise terminate the algorithm; return d min and φ min Models used in scalability studyModelInadequacies ε f h and ε g h m , ∆ P , ∆ r k , d k and D k

Table 8 :
Cost ratios among different models

Table 9 :
(MFO /RFO) cost ratios of the conducted optimization runsRunsCost ratio Reduced cost n int ∆ r hl (d 0 ) Error%

Table 10 :
Design variables of the three optimization design problems

Table 11 :
The total number Results of the three design problems for 3, 7, and 11 optimization variables using either only the high-fidelity reference model M 0 or multiple models managed by our relative adequacy framework