Optimum aerodynamic shape design under uncertainty by utility theory and metamodeling

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.1 (1-17)

Aerospace Science and Technology ••• (••••) ••••••

1

Contents lists available at ScienceDirect

67 68

2 3

Aerospace Science and Technology

4

69 70 71

5

72

6

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13

Optimum aerodynamic shape design under uncertainty by utility theory and metamodeling

16 17

81

Department of Aerospace Engineering, Iowa State University, Ames, IA 50011, United States of America

83

82 84 85

19

a r t i c l e

i n f o

a b s t r a c t

86 87

21 22 23 24 25

Article history: Received 28 February 2019 Received in revised form 5 October 2019 Accepted 8 October 2019 Available online xxxx

26 27 28 29 30 31

79

Xiaosong Du, Leifur Leifsson ∗

18 20

78 80

14 15

77

Keywords: Aerodynamic shape optimization Design under uncertainty Robust design Utility theory Aerodynamic metamodeling

32 33 34 35 36 37 38 39 40

In this work, utility theory is introduced to optimum aerodynamic shape design under uncertainty (also called robust aerodynamic design). Specifically, utility theory is used to formulate the objective function of the optimization problem. The advantage of the proposed approach over the commonly used weighted sum method is that it does not require the use of weighting factors or the addition of constraints on the statistical moments. The polynomial chaos expansion metamodel with the least-angle regression and the hyperbolic truncation scheme are used to accelerate the uncertainty propagation and the optimization process. The proposed approach is demonstrated on the optimum airfoil shape design under uncertainty at steady transonic conditions. Two design cases are considered. In both cases the drag coefficient is to be minimized subject to constraints on the lift coefficient and the airfoil thickness at two chordwise locations. The first case treats the Mach number as uncertainty parameter and the second case treats the Mach number and the lift coefficient at target as uncertain. The proposed method is compared with deterministic single-point optimization and multi-point optimization, and the standard robust design formulations. Results show that the proposed approach is capable of obtaining airfoil design of the lowest mean drag coefficient with the smallest standard deviation as compared to the other approaches. In the first case, the proposed approach yields a design with a mean drag coefficient comparable to the standard robust design (around 97.6 drag counts, cts), but with a significantly lower variance (0.6 cts versus 2.5 cts). In the second case, the proposed approach yielded a design with a higher mean drag coefficient than the standard approach (98.4 cts versus 97.6 cts), but with a lower standard deviation (4.1 cts versus 4.8 cts). In both cases, the deterministic approaches yield designs with comparable or lower drag coefficients but with significantly larger standard deviations. © 2019 Elsevier Masson SAS. All rights reserved.

88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106

41

107

42

108

43

109 110

44 45

1. Introduction

46 47 48 49 50 51 52 53 54 55 56 57 58 59

Simulation-based design optimization (SDO) [1–3] is important in aerodynamic design. In SDO, predictive simulations based on the governing flow physics drive the aerodynamic design. In particular, in a SDO framework the optimization technique proposes a new design and the predictive simulations are used to evaluate the performance and feasibility of the design. The process is iterated until convergence on an optimum design is achieved. The aerodynamic SDO process, typically, requires many evaluations of the predictive simulations, which is a major challenge if each evaluation is timeconsuming. A large body of work on SDO methods for aerodynamic design has been reported in the literature. The majority of the works assume that all system parameter values are explicitly known, i.e.

60 61 62 63 64 65 66

*

Corresponding author. E-mail address: [email protected] (L. Leifsson).

https://doi.org/10.1016/j.ast.2019.105464 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

any uncertainty in the parameters is neglected. Examples of such works include the ones by Jameson et al. [4–7] where the authors developed and improved the continuous and discrete adjoint methods for single-point optimum aerodynamic design in steady and unsteady flows at cruise and high-lift conditions. Here, the term ‘single-point’ refers to optimum design at a single operating condition. Lyu and Martins [8,9] integrated a gradient-based optimization algorithm with the continuous adjoint method for efficient single-point three-dimensional optimization of an adaptive morphing trailing edge wing on the blended-wing-body aircraft. Grasso et al. [10–12] implemented single-point airfoil shape design optimization of wind turbine rotor blades. Poole et al. [13, 14] investigated the influence of shape parameterization techniques based on a single-point benchmark aerodynamic optimization problem. Palacios et al. [15,16] proposed optimal single-point rolling maneuvers with flexible wings in order to identify actuation strategies for improving the vehicle aerodynamic maneuverability. Rogers et al. [17,18] optimized the single-point aerodynamic performance of guided projectiles. Reuther et al. [19] described

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105464 /FLA

1 2 3 4 5 6

A A0 A min

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

const Cd ??? Cl ??? C l∗ cts ??? d l g h H i

24 25 26 27 28

67

Nomenclature

7 8

[m5G; v1.261; Prn:14/10/2019; 15:48] P.2 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

2

l M MPC

non-dimensional cross-section area of current design, scalar [-] non-dimensional cross-section area of baseline, scalar [-] non-dimensional minimum cross-section area, scalar [-] constant value set on inequality constraints, scalar [-] drag coefficient, scalar [-] 2 d/(1/2ρ∞ V ∞ S) lift coefficient, scalar [-] 2 l/(1/2ρ∞ V ∞ S) lift coefficient at target, scalar [-] drag count C d = 1E − 4 drag force, scalar [N] lift force, scalar [N] vector containing inequality constraints vector containing equality constraints objective function, scalar index of optimal design found within each optimization step vector containing lower bounds response from actual model, scalar [-] prediction from polynomial chaos expansions, scalar [-]

68

M∞ N testing p p∞ P S u u V∞ x Ytesting

ˆ Y

α α ε εL O O ρ ρ∞ λ

ω 

free-stream Mach number, scalar [-] total number of testing points [-] static pressure at evaluated point, scalar [Pa] free-stream static pressure, scalar [Pa] total number of truncated polynomial basis terms reference area, scalar [m2 ] vector containing upper bounds on design variables speed relative to fluid [m/s] free-stream velocity [m/s] vector of design variables, vector [-] observations at testing points, vector [-] PCE predictions at testing points, vector [-] angle of attack, scalar [deg] coefficient vector of stochastic expansions, vector [-] residual in least-squares stochastic expansions, scalar [-] leave-one-out error, scalar [-] density of fluid, scalar [kg/m3 ] free-stream fluid density, scalar [kg/m3 ] penalty factor, scalar [-] weighting factor, scalar [-] general basis function of stochastic expansions, scalar [-] orthogonal basis function of polynomial chaos expansions, scalar [-]

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 96

30 32

70

95

29 31

69

the derivation, development, and practical application of adjointbased efficient aerodynamic shape optimization using control theory, involving multi-point aircraft configurations subject to geometric constraints. Here, ‘multi-point’ refers to optimum design at a given operating conditions as well as at various off-design conditions [20]. Leung et al. [21,22] applied the Newton-Krylov algorithm on aerodynamic shape optimization in three dimensions for multi-point optimization over a range of operating conditions of interest. Buckley et al. [23,24] applied multi-point optimum design to balance the aerodynamic performance of an airfoil over a range of operating conditions. Madsen et al. [25] proposed a three-dimensional computational fluid dynamics (CFD) adjointbased multi-point shape optimization of a modern 10 MW offshore wind turbine taking a range of wind speeds into consideration for a more robust design. The aerodynamic system parameters may, however, vary due to manufacturing, environmental, and operational uncertainties. The uncertainties can be classified into aleatory and epistemic types [26,27]. Aleatory uncertainty, also known as statistical uncertainty or inherent uncertainty, comes from a random process and is due to the lack of knowledge of the physical process. Aleatory uncertainty always exists in reality and cannot be eliminated. It can be described by empirical or pre-assumed probabilistic density functions [28]. In contrast, epistemic uncertainty, also known as systematic uncertainty, is due to the incomplete or imprecise understanding of the underlying physical process. The epistemic uncertainty is theoretically reducible by adding more observations or improving the physical models of the system under consideration. Epistemic uncertainty can be represented using interval analysis [29], possibility theory [30], Bayesian theory [31], or evidence theory [32]. Accounting for uncertainties, as well as incorporating an appropriate treatment, as part of the analysis of complex systems is important when reaching decisions in the design process, see, e.g., [33–35]. Considering uncertainty in optimal design is of practical importance and value and has attracted attention in the engineer-

ing design community. The concept of design under uncertainty was proposed by Taguchi [36,37] in the 1950s. It aims at reducing the effects of uncertain parameters on the system performance by considering the key figures of interest within a given uncertain parameter space [38,39]. Jafarsalehi et al. [40] demonstrated the practical value of nondeterministic optimization over deterministic optimization for a small satellite in low Earth orbit with remote sensing mission. Considering the mean of the original objective function over the uncertain space subject to constraints on the mean of the design variables provided reliable designs of probabilistic analysis. Calafiore et al. [41] implemented both expected value and worse case of original objective functions for optimum design. In addition, the authors tackled semidefinite optimization program with many constraints by means of sampling-based probabilistic relaxation, which may be extended to practical problems for large-scale applications. Eshghi et al. [42] noticed the sensitive dynamics of vibration harvesters with respect to the uncertainty of physical parameters and proposed a reliability-based design optimization containing probabilistic constrains as a solution to the problem. Zhang et al. [43] addressed a network system design problem by keeping the objective function deterministic with probabilistic constraints, which can be extended to more realistic scenarios. Aoues et al. [44] introduced a reliability-based design optimization for economical and reliable structures by decomposing the design under uncertainty into several cycles of deterministic design optimizations where each is based on the optimal system safety factors. Yi et al. [45] proposed a multidisciplinary robust design optimization method for the design of an air-cooling battery thermal management system. The proposed approach is based on the set strategy handling the coupled variables between the subsystems. The optimum design under uncertainty has also been widely implemented into aerodynamic and aero-structural systems in recent years. Papadimitriou et al. [46] propagated the probabilistically distributed flow-related and geometrical uncertainties through aerodynamic simulation models for airfoil shape optimiza-

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.3 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

tion by taking the weighted sum of the mean and standard deviation of the drag coefficient as the objective function subject to a constraint on the reliability of lift using the first order reliability method. Vuruskan and Hosder [47,48] investigated the impact of two commonly used turbulence models with three different shape parameterization techniques on deterministic and robust aerodynamic airfoil shape optimum designs. The robust design problem formulation was to minimize the weighted sum of the mean and standard deviation of the drag coefficient subject to a constraint on the minimum mean lift coefficient. Results showed that robust design tended to reduce the impact of modeling and parameterization selections. Liem et al. [49] emphasized the importance of multi-condition aircraft optimization and noticed the heuristic choice on conditions and weights. The authors proposed a new approach minimizing expected drag value given by a probability density function in the conditional parameter space, which can be estimated based on actual flight data. Cook et al. [50,51] applied the horsetail matching method to formulate the optimization under uncertainty for aero-structural design of transonic wings. The advantage of the approach is to avoid the specification of the probability distributions on the uncertain parameters when the information is not available. Furthermore, the approach escapes stochastically dominated designs that are obtained by the existing problem formulations, which treat the statistical moments as separate objective functions. Caboni et al. [52] presented different approaches for wind turbine airfoil shape optimization under uncertainty subject to a set of aerodynamic and structural constraints. The authors applied both probability- and imprecise probabilitybased methods with the results compared against state-of-the-art design methods to demonstrate the advantages of optimum design under uncertainty. Abhiram et al. [53,54] introduced the signal/noise ratios and the Taguchi method for robust aerodynamic blade-shape design of a two-bladed small-scale unmanned helicopter for optimal hover performance. The proposed approach resulted in a systematic reduction in the number of numerical simulations while providing variance index for a robust solution. It is evident that significant progress has been made on optimum design under uncertainty in various engineering fields. The objective function, however, is not thoroughly investigated, and is commonly formulated using the weighted-sum (WS) [55,56] of the mean and the standard deviation of the figure of interest. The WS method, however, suffers from the following key weaknesses: (1) it is challenging to choose the weighting factors [57], (2) the estimate of the variance lacks an estimate of its confidence [57], and (3) the approach may not be able to reach certain Pareto optimum designs [58–61]. A few other methods, such as actual flight databased WS method [49], compromise programming (CP) [62] and the probabilistic robustness index method [63], have been developed to formulate the problem in a different manner. However, actual flight data-based WS method might require a large number of representative information, the CP method captures the Pareto optimum designs better but still has to determine weighting factors for single-object optimization, and the probabilistic robustness index method needs the estimation of variance still. Utility theory (UT) [64,65] as a typical decision-making method in finances has not been investigated for aerodynamic optimization area. This work introduced the UT method for the first time to formulate the objective function of the optimum aerodynamic shape design under uncertainty. In addition, the least-angle regression (LARS)-based polynomial chaos expansion (PCE) [66–68] metamodel is implemented as a computationally efficient model to quantify statistics of deterministic objective function. In this work, UT [69,70] is used to account for the designers’ risk preference in the objective function during the decision-making process of optimum design under uncertainty, without having to determine the weighting factors or constraining the statistical moments. The de-

3

sign constraints, both the inequality and equality constraints, are cast as nondeterministic by using their mean values. The proposed approach is demonstrated on the simulationbased optimum airfoil shape design under uncertainty at steady transonic conditions. The physics-based fluid flow simulator solves the Reynolds-Averaged Navier-Stokes (RANS) equations with the Spalart-Allmaras turbulence model. The numerical examples are transformed from Benchmark Case II developed by the Aerodynamic Design Optimization Discussion Group (ADODG) to fulfill the need for aerodynamic design optimization benchmarks where typical and novel methods can be compared. With the satisfactory comparison results obtained on the benchmark, effective and efficient application to more practical and realistic optimum engineering design scenarios under uncertainty can be expected. The paper is organized as follows. Section 2 describes the current state of the art methods for formulating the objective function optimum aerodynamic design under uncertainty. The section also derives the UT-based objective function. Section 3 gives the details of the aerodynamics modeling and the LARS metamodeling. Section 4 gives the results of two numerical airfoil robust design cases, where the proposed UT-based objective formulation is demonstrated and compared with the current state of the art deterministic optimization and standard robust design formulations. Section 5 concludes the paper.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91

2. Optimum aerodynamic shape design problem formulations

92 93

This section describes the generalized formulation of optimization problem, based on which deterministic design formulations and robust design formulations are given. A new robust aerodynamic design problem formulation using utility theory is proposed. The key terminologies and risk preferences of utility theory are also introduced in this section.

94 95 96 97 98 99 100

2.1. Deterministic optimum design formulations

101 102

Deterministic optimization problems can be generalized as follows

x∗

106

l≤x≤u

j = 1, . . . , M , k = 1, . . . , N ,

104 105

= arg min H (f(x)),

s.t. g j (x) ≤ 0, hk (x) = 0,

103

(1)

107 108 109 110

where x is the design variable vector, f is the vector with the figures of merit, H is the scalar valued objective function, g is a M-dimension vector of inequality constraints, and M is the total amount of inequality constraints, h is a N-dimension equality constraints, and N is the total amount of equality constraints.

111 112 113 114 115 116

2.1.1. Single-point deterministic design Single-point optimization (SPO) method [71–73] has the same format as (1), but only at one set of operating parameters, e.g. at one specific free-stream Mach number. The SPO formulation is as follows

119 120 121 123

l≤x≤u

j = 1, . . . , M , k = 1, . . . , N ,

118

122

x∗ = arg min f (x), s.t. g j (x) ≤ 0, hk (x) = 0,

117

(2)

where f (x) is the drag coefficient in this work.

124 125 126 127 128 129

2.1.2. Multi-point deterministic design Multi-point optimization (MPO) method [74,75] considers multiple sets of operating parameters, e.g. free-stream Mach numbers

130 131 132

JID:AESCTE

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.4 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

4

1 2 3 4

and lift coefficients at the target, which increases the robustness of optimal design. Its objective function is formulated in a weightedsum form. In this work only equal weighting factors are considered, as follows

67 68 69 70 71

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

x∗ = arg min l≤x≤u

s.t. g j (x) ≤ 0, hk (x) = 0,

N MPO

72

ωi f i (x),

73

i =1

76 77

where N MPO is the total operational conditions, f (x) is the drag coefficient, and ω is weighting factor.

78

2.2. Standard optimum design under uncertainty formulations

81

The first type of standard design under uncertainty formulation is the weighted-sum (WS) method [76–78], which takes the weighted sum of mean (μ) and the standard deviation (σ ) of model response as objective function. “μ + σ ” does not emphasize weight on σ , making the optimal design not robust enough (see the numerical examples), so only the weighting factor (equal or greater than one) on the standard deviation is considered in this work, shown as follows

x∗



 = arg min(μ( f (x)) + ω · σ f (x) ), l≤x≤u

s.t. 



μ g j (x) ≤ 0, j = 1, . . . , M , μ hk (x) = 0, k = 1, . . . , N ,

(4)

where ω is the weighting factor, μ is the mean, and σ is the standard deviation. The second type of standard design under uncertainty formulation is to minimize the mean of model response, setting constraints on the standard deviation. This formulation can guarantee the mean to be as small as possible in the feasible region as follows

x∗ = arg min μ( f (x)), l≤x≤u

s.t. 

 σ  f (x) ≤ const, μ g j (x) ≤ 0, j = 1, . . . , M , μ hk (x) = 0, k = 1, . . . , N .

50 51 52 53 54

(5)

2.3. Optimum design under uncertainty by utility theory The utility theory (UT) utilized in this work is the expected utility theory (EUT) [79–81]. It was first proposed for the solution to the Saint Petersburg paradox [79], a game in which a fair coin will be tossed until getting a head. The player can win 2n dollars, where n is how many tosses to reach the first head. The expected outcome of this game is

55 56 57 58 59 60 61 62 63 64 65 66

ESP =

∞  n =1

∞  pn · 2 ≡ (1/2)n · 2n , n

79 80 82

48 49

75

j = 1, . . . , M , k = 1, . . . , N ,

46 47

74

(3)

(6)

n =1

where pn is the probability of getting the first head at nth toss, making (6) goes infinite. The paradox is that no one would pay any substantial sum to play this game. Daniel Bernoulli suggested that the expectation was relevant, but the approximate expectation involved a measure of the personal worth, or utility [82], of the money at stake, rather than the expectation of the monetary amounts themselves. This represents the earliest suggestion that numbers can be attached to monetary amounts to represent the

83 84 85

Fig. 1. Risk aversion, risk neutral, risk loving plots, based on utility function.

86 87

personal worth of these amounts, and the decisions under risk are made based on the expectation of this personal worth. This concept of the maximization of expected utility worked as the basis for decision making under risk, then got fully developed by Von Neumann and Morgenstern [83–85] and others. Since the theory of Von Neumann – Morgenstern expected utility [83–85], it has been the dominant and normative theory [79] of decision making under uncertainty, due to the requirement on consistency of the behavior of any rational person when considering the preference. Two strict axioms, transitivity and independence, are the basic rules to start with when talking about expected utility theory. The transitivity means that if action a > action b (a is preferred to b), and action b > action c, then action a > action c. The independence means that if a ≥ b (a is preferred to or equivalent to b), then the action e offering a with a probability p and c with probability 1 − p ≥ the action f offering b with probability p and c with probability 1 − p, for any c. Once the consistency axioms are fulfilled, a numerical quantity, generally referred to as utility, can be assigned to each outcome, with the “preferred” action being which has the highest expected utility. This expected utility is indeed the weighted sum result of utilities of all possible outcomes, with the probability of each outcome as weighting factors. The utility of each outcome and preference of individual are set with the utility function [86,87], within which one can change the function parameter for various risk preferences [88–92], namely, high/low risk aversion, risk neutral and low/high risk loving. These five types of risk preferences refer to more/less curved concave, linear, and less/more curved concave curve in mathematics, respectively, as shown in Fig. 1. Utility theorists define the risk-aversion utility function as if and only if for all potential probabilities p, E p {U } ≤ U ( E p {x}), while risk-loving utility function as if and only if E p {U } ≥ U ( E p {x}), where U is the selected utility function. Taking risk aversion as example, we can represent the definition mathematically as





u (x) f (x)dx ≤ u



xf (x)dx ,

88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122

(7)

123 124

where f (x) is the probability density function of x, and this mathematical expression is actually the Jensen’s inequality [80], which is the defining property of concave functions. Similarly, it is straightforward to understand that convex curve represents risk-loving preference. Fig. 2 explains the risk aversion with V 1 and V 3 having the same probability and V 2 as the mean of V 1 and V 3 . In this work, we apply a commonly used utility function, see e.g. Kannan et al. [93], which is the following

125 126 127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.5 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

5

the novel UT-based formulation proposed in this work. A summary of the formulations is given in Table 1.

1 2

3. Simulation-based aerodynamic design optimization framework

4 5

70 71 72

6

Quantifying the statistical moments, which will be required for standard robust design methods, using real physics model is computationally expensive. Here we proposed a type of stochastic metamodels to alleviate the computational costs. In this section, details on computational fluid dynamics model and the polynomial chaos expansions (PCE) are provided.

7 8 9 10 11 12

16

This section talks about the Reynolds-Averaged Navier-Stokes (RANS) equations, corresponded turbulence model, shape parameterization, solver and grid convergence study.

17 18

Fig. 2. Concave curve-based risk aversion.

19 20

25 26 27 28

U (V ) = −

EU =



31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

51 52

1 −α V e ,

(8)

α

U ( V i ) · P ( V i ),

(9)

i

where V i is the measurement associated with option i, and P ( V i ) is probability of that measurement happens. Certainty equivalent (CE) is defined as

CE = U

−1

(EU ).

(10)

The certainty equivalent is the minimum value that one would accept, so it should be maximized in typical economic problems. Since the utility function is monotone, it is the same to optimize expected utility and certainty equivalent. In this work, we will work on the CE, because it represents the quantity of interest directly. To convert it into a minimization problem on the drag coefficient, the utility function is fed with negative drag coefficient values. The negative CE needs to be minimized in optimizer, as follows

x∗ = arg min −CE(−C d ), l≤x≤u

49 50

s.t. 



(11)

μ g j (x) ≤ 0, j = 1, . . . , M , μ hk (x) = 0, k = 1, . . . , N .

53 54

2.4. Summary of problem formulations

55 56 57

Previous parts of this section describe detailed formulations of objective function in deterministic and robust design, as well as

3.1.1. Airfoil shape parameterization B-spline [94,95] is a piecewise polynomial function of the order k in a variable t defined over the range of [t 0 , tm ], m = k + 1, where the points t = t j are known as knots. A spline function of order k on a given set of knots can be expressed as follows

p (t ) =

n 

60

77 78

82 83 84

(12)

i =0

N i ,1 =

87 88 89 90 91

1, if t i ≤ t ≤ t i +1 0,

93 94

where N i ,k are B-spline basis functions, and P i , i = 0, 1, . . . , n, are control points. The basis functions are calculated from the CoxDeBoor recursion relation as



86

92

N i ,k (t ) P i ,

otherwise

95 96 97 98 99

(13)

,

100 101 102

and

N i ,k (t ) =

t − ti t i +k−1 − t i

N i ,k−1 (t ) +

t i +k − t t i +k − t i +1

103

N i +1,k−1 (t ),

(14)

104 105

for k = 2, 3, . . . , K , as well as for all needed values of i. The basis functions N i ,k can be polynomials of the order one, two, or higher. The design variables in this work are the vertical coordinates of six control points in a B-spline airfoil shape parameterization since this number can keep the relative L 2 norm error close to 1% for typical existing airfoils. Corresponding upper and lower bounds of (1 + 25%)xbaseline and (1 − 25%)xbaseline , respectively (Fig. 3).

106

3.1.2. Governing fluid flow equations The flow past the airfoil is taken to be steady, compressible, and fully turbulent. The RANS equations are taken as the governing equations. For steady compressible flow the RANS equations can be written (in Cartesian tensor form) as follows [96]

114

∂ (ρ¯ u¯ j + ρ  u j ) = 0, ∂xj

120

(15)

107 108 109 110 111 112 113 115 116 117 118 119 121 122

and

123 124

58 59

76

85

where V denotes the scalar objective function, U is the utility function of V , α is the constant controlling risk reference. Positive α represents risk aversion, while negative α represents risk loving in (8). Some basic terminologies are shown as following, for the use of decision making. Expected utility (EU) is defined as

29 30

75

81

15

24

74

80

3.1. Computational fluid dynamics model

14

23

73

79

13

22

68 69

3

21

67

125

Table 1 Summary on formulations of objective function.

126

61

Method

Objective function

Inequality constraints

Equality constraints

127

62

SPO MPO WS-based σ -constrained UT-based (this work)

f (x)

g(x) g(x) μ(g(x)) σ ( f (x)), μ(g(x))

h(x) h(x) μ(h(x)) μ(h(x)) μ(h(x))

128

63 64 65 66

n

i =1 ωi f i (x) μ( f (x)) + ω · σ ( f (x)) μ( f (x))

−CE

μ(g(x))

129 130 131 132

JID:AESCTE

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.6 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

6

1

67

Table 2 Grid convergence study based on RAE 2822.

2 3

68

Mesh

Grid size

C l (cts)

C d (cts)

Simulation timea (min)

7

1 2 3 4

9,836 38,876 154,556 616,316

0.5 0.5 0.5 0.5

101.6 99.7 103.7 103.4

1.9 6.4 36.7 193.9

8

a

4 5 6

70 71 72 73

Computed on a high-performance cluster with 32 processors. Flow solution only.

9 11

cb1 = 0.13355,

12

c w1 =

13 14 15

cb1

κ2

+



17

fw = g

18 19

cb2 = 0.622,

(1 + cb2 )

κ = 0.41,

16

σ f v1 =

1 + c 6w3 g 6 + c 6w3

 ,

χ3 χ 3 + c 3v1

22

Fig. 3. Airfoil shape parameterization.

S˜ = S +

23

26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

∂ p¯ ∂ ∂ (ρ¯ u¯ i u¯ j ) = − + (τ¯i j − ρ u i u j ), ∂xj ∂ xi ∂xj

where τ¯i j is the time averaged viscous stress tensor, defined as



  ∂ u¯ j 2 ∂ u¯ k ∂ u¯ i − δi j τ¯i j = μ + ∂xj ∂ xi 3 ∂ xk

    ∂ u j ∂ ui 2 ∂ u  − δi j k , + +μ ∂xj ∂ xi 3 ∂ xk

50

53

(17)

−ρ u i u j is commonly called the Reynolds stress which represents the effects of turbulence and needs to be modeled to close this equation. Here we selected the one equation model, Spalart-Allmaras (SA) model [97], solving a single partial differential equation for the Reynolds stress. The SA model is considered to be accurate for attached wall-bounded flows and flows with mild separation and recirculation. It employs the Boussinesq approach, providing the equation for the Reynolds stress as

−ρ u  u  i

j

 = 2μ T

51 52

 ∂u j ∂ ui 1 ∂ u¯ k 2 + − δi j − ρ¯ k¯ δi j , ∂xj ∂ xi 3 ∂ xk 3

(18)

56 57 58 59 60 61 62 63 64 65 66

μT = ρ v˜ f v1 ,

(19)

where the parameter v˜ is obtained from the solution of the transport equation

 ∂ v˜ ∂ v˜ 1 ∂ ˜ (v + v ) + cb1 (1 − f v2 ) S˜ v˜ uj = ∂xj σ ∂ xk ∂ xk  2 v˜ cb2 ∂ v˜ ∂ v˜ − c w1 f w + , d σ ∂ xk ∂ xk where the closure coefficients and functions are given by

2 2 κ d

f v2 ,

81 82 83 84

(21)

v

S=

 ∂u j ∂ ui , i j = − 2 ∂xj ∂ xi 1

80

v˜

r=

v˜ S˜ κ 2 d2

78 79

χ f v2 = 1 − , 1 + χ f v1

,

77

85 86

,



2 i j  i j ,

and d is the distance from the closest surface.

87 88 89 90 91 92 93 94

3.1.3. Solver and grid convergence study In this work, the Stanford University Unstructured (SU2) [98–100] is utilized as the numerical solver, which can accelerate the calculation when using multi-grid method built inside. Since the test cases will be drag-reduction optimization, convergence criteria is set as CAUCHY on drag, and epsilon as 1E−5, meaning that the calculation is viewed as convergence if the change on drag coefficient is less than or equal to 1E−5 within 100 iteration steps. For mesh generation, efficient hybrid mesh is transferred to the format of pointwise mesh, for the use of SU2. It follows the curvilinear body-fitted grid of C-topology, i.e. at the inlet of the grid surrounds the leading edge of the airfoil but is open at the other end (Fig. 4). To determine the size of mesh, while keeping best efficiency with the converged-level accuracy, grid converge study on RAE 2822 is made with Mach number set as 0.734, lift coefficient at target as 0.5, in fixed-C l mode [101,102], as shown in Table 2. The difference of drag coefficient between Mesh 3 and Mesh 4 is 0.3 cts. Mesh 3 is selected because it has more than 5 times higher efficiency than Mesh 4.

96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117

where k¯ is the kinetic energy of turbulence, and μ T is the turbulence viscosity which is defined by the SA turbulence model as

54 55

75

95

and u¯ i is the time averaged velocity component (i = 1, 2, 3), ρ¯ is the time averaged fluid density, p¯ is the time averaged pressure, and μ is the fluid dynamic viscosity, δi j is Kronecker delta,

48 49

(16)

v˜

σ = 2 /3 ,

c w3 = 2,

χ= ,

g = r + c w2 (r 6 − r ),

21

c v1 = 7.1,

c w2 = 0.3,

,

20

25

74 76

10

24

69

(20)

3.2. Polynomial chaos expansion metamodeling

118 119

Stochastic expansions [103] were initially developed to solve for stochastic partial differential equations. It addresses equations with stochastic inputs and is mainly applied to computational physics problems. In this field, stochastic Galerkin methods employ polynomial chaos expansions (PCE) to represent the solution and inputs to stochastic differential equations. However, for design under uncertainty, this method performs poorly, due to its inefficient requirements on intrusiveness of derived set of differential equations. Thus, another type of method, non-intrusive [104] stochastic expansions, treating existing computational models as a black box, helps make construction process more efficient and effective. In this work, the non-intrusive stochastic expansions, polynomial chaos expansions, are generated [105–108].

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.7 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84 85

19

Fig. 4. Structured curvilinear body-fitted C-topology grid: (a) solution domain; (b) mesh overview.

20

86 87

21 22

Stochastic expansions are generalized in the format of

23 24 25 26 27 28 29 30 31 32 33 34 35

M (x) =

P 

where M is dependent quantity,  is basis function corresponded with the distribution of uncertain parameters, α is coefficient of basis function, here x is a vector of deterministic and uncertain parameters, i is index of each polynomial term, P is total number of polynomial terms. In this work, Legendre polynomial basis is utilized because we assume uniformly distributed uncertain parameters. When constructing polynomial chaos expansions, P does not go to infinity, but takes a truncated form:

M P C (x) ≈

38 39 40 41 42 43 44 45 46 47 48 49 50 51

P −1 

αi i (x),

where M is the prediction from PCE model, x is a vector, containing design variables, and

(n + p )! P= , n! p !

(24)

where n is the number of design variables, p is the required order of PCE, this gives the coefficient matrix as

AM , p =

56 57 58 59



α ∈ N M : α 1 ≤ p ,

60 61 62 63 64 65 66

(25)

where λ is penalty factor, and ∗ 1 is L 1 norm. Then, the relationship between actual response and PCE prediction is

M (x) = M P C (x) + ε =

P −1 

α j j (x) + ε ≡ α T (x) + ε,

where, M is response from actual model, ε is the residual. PCE can be constructed as least-squares minimization problem, by minimizing ε as follows

2  αˆ = arg min E α (x) − M (x) . T

(27)

α

Adding a regularization penalty term to favor low-rank solution [109], due to sparsity-of-effects principle, claiming that low-order interaction among inputs are important.



2 

αˆ = arg min E α T (x) − M (x) α

εL O O =

(i ) P C \i (i ) 2 (x )) i =1 ( M (x ) − M . N ( i ) P C (x(i ) ))2 i =1 ( M (x ) − M

+ λ α 1 .

(29)

In practical use, there is no need to construct the PCE metamodel N times, for the calculation of ε L O O . (29) can be simplified as

εL O O =

N  

(28)

88 89 90 91 92 93 94 95

2 N 2 M (x(i ) ) − M P C (x(i ) )    (i )  ˆ Y , (30) −μ M x 1 − hi i =1

ˆ Y is the estimated mean of model response. where μ The general process of this basis-adaptive method can be summarized as

96 97 98 99 100 101 102 103 104 105 106 107 108 109

(1) (2) (3) (4)

Set error threshold ε T , and maximum iteration NImax , Start with low-order PCE, e.g. 2nd order, Calculate PCE coefficients using least-squares method, Calculate ε L O O and check the convergence on ε T and NImax .

110 111 112 113 114

Least-angle regression (LARS) [110] is an iteration-based linear regression method, based on correlation between residual and orthogonal PCE basis function. There are two main steps, initialization and iteration, in LARS method, as discussed below.

115 116 117 118 119

(26)

j =0



N

i =1

PC

54 55

(23)

i =0

52 53

(22)

i =0

36 37

αi i (x),

3.2.1. Basis-adaptive Least-Angle Regression Sparse (LARS) method Basis-adaptive technique [66–68] allows the maximum degree of PCE effectively from available data. It constructs PCE, starting from low order, which only asks for a very small number of sample points, then keeps increasing the maximum degree of PCE until the leave-one-out (LOO) residual is small enough or all points have been used. The leave-one-out residual is

(1) Initialization: (a) Set all PCE coefficients α as 0, (b) Put all preset orthogonal PCE basis into a set called candidate set, (c) Define an empty set, called active set, (d) Current residual at the end of initialization is actually the observations. (2) Iteration: (a) Find the PCE basis, i , which is the most correlated with current residual, (b) Move the corresponding coefficient, αi , of i , towards the resulting coefficient of least squares, until another basis,  j has the same correlation with current residual,

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.8 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

8

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84 85

19

Fig. 5. Hyperbolic truncation scheme on PCE: (a) q = 1 (standard truncation); (b) q = 0.5.

20

86 87

21 22 23 24 25

(c) Repeat Step 2(b), until another basis is found to have the same correlation, (d) Keep iterating, until the active set gets the required number of samples or the total number of observations.

88 89 90 91 92

26 27 28 29 30 31 32 33

The two families of adaptivity are not mutually exclusive, and the combination, basis-adaptive LARS method [68–70], has already been developed for application. The process starts from low-order PCE, then calculates coefficients using LARS method in turn to get leave-one-out (LOO) residual, finally check this residual with preset error threshold, keep increasing the maximum degree of PCE until LOO residual is small enough.

93 94 95 96 97 98 99 100

34 35 36 37 38 39 40 41

A M , p ,q =

44

47

α q =

52 53 54 55 56 57 58 59

106

Fig. 6. Flowchart of the PCE-based design under uncertainty algorithm.

1/q

αiq

after further transformation, the formula becomes

 A M , p ,q =

α ∈ A M,p :

N   i =1

αi p

q

 ≤1 .

(33)

Fig. 5 shows the detailed explanation of hyperbolic truncation method, setting the order of a 2-dimensional PCE as p = 3 and the norm order q = 1 and 0.5. Circles represent selected orders of PCE terms, and stars are those terms truncated by hyperbolic schemes. It is straightforward to see that when q = 1, the hyperbolic truncation corresponds the same as the standard truncation (25), while q = 0.5 gives a further truncation. In this work, q = 0.75 for both robust airfoil design cases.

[112,113] and the adjoint information solved by SU2. The optimum design under uncertainty problems (cf. Sections 2.2 and 2.3) are solved by an optimization framework based on off-line generated metamodel. A flowchart of our framework is given in Fig. 6. The six B-spline control points are the geometry design variables and operational parameters (Mach in Case 1, Mach and C l at target in Case 2) are selected as uncertain parameters. In this work, these variables are assigned with Uniform distributions. An initial sampling plan is generated using Latin hypercube sampling (LHS) [114]. Then, the observations are gathered using the CFD simulation model (cf. Section 3.1). Given the sampled data, a PCE-based metamodel is constructed (cf. Section 3.2) and validated using the normalized root mean squared error (NRMSE) given by

3.3. Workflow

62 63 64 65 66

The deterministic single- and multi-point optimization problems considered in the work (cf. Section 2.1) are solved using the sequential least squares programming (SLSQP) algorithm in pyOpt [111], which uses the Han-Powell quasi-Newton method

=

110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

NRMSE

60 61

108 109

(32)

,

i =1

50 51

105

(31)

107

48 49

α ∈ A M , p : α q ≤ p ,

 N 

45 46

103 104



102

where

42 43

101

3.2.2. Hyperbolic truncation method The hyperbolic method [66–68] utilizes the q-norm scheme to define the truncation, mathematically can be expressed as

  N testing i  i 2  i =1 (Yˆ − Y testing ) 

125

N testing

128

 max(Y testing ) − min(Y testing ) ,

126 127

(34)

129

where N testing is the total number of testing points, Y testing is the

131

observation of testing points, and Yˆ is the PCE prediction. If the

130 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.9 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7

9

NRMSE value is above 2%, a new set of samples are generated using LHS with 50 more training points and the loop of data collection and metamodel construction is repeated. Once the metamodel is accurate enough, i.e. the NRMSE is less than or equal to 2%, the optimum design under uncertainty is performed using the metamodel and the sequential quadratic programming (SQP) algorithm in MATLAB [115].

67 68 69 70 71 72 73 74

8 9

75

4. Numerical examples

76

10 11 12 13 14 15

77

In this section, the UT-based objective formulation based on LARS-based PCE metamodel is applied to robust design cases. The results are compared with the deterministic optimization formulations and standard robust design formulations as described in Section 2.

78 79 80 81 82

16 17 18

4.1. Case 1: RAE 2822 airfoil in transonic viscous flow with uncertain Mach number

Fig. 7. Metamodel validation of Case 1.

21 22 23 24 25 26 27 28 29

4.1.1. Problem formulation This case aims at obtaining an optimal airfoil shape, which is parameterized by six B-spline control points and is insensitive to the uncertain Mach number, statistically distributed as U (0.70, 0.75). Computational model is run at fixed-C l mode [116], with two thickness constraints at 20% and 80% chord length station. The original problem setup is solved at a free-stream Mach number of 0.734 as

min C d (x),

30

l≤x≤u

31

s.t.

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

h ≡ C l∗ = 0.5, g 1 ≡ t /c x/c =0.2 ≥ t /c R A E2822,x/c =0.2 , g 2 ≡ t /c x/c =0.8 ≥ t /c R A E2822,x/c =0.8 .

(35)

For deterministic optimization, all operating parameters are pre-fixed values. Single-point optimization is run three times, at a M ∞ of 0.700, 0.725, and 0.750, respectively. Multi-point optimization runs at a range of M ∞ , namely, 0.700, 0.715, 0.725, 0.735, 0.745, and 0.750. For robust design, M ∞ is set as U (0.70, 0.75). Standard robust design formulations are set up with increasing weighting factors or decreasing minimum required constraints on standard deviation. The UT-based robust design focus on the negative CE value of the utility function (8) as introduced in Section II.A.4. 4.1.2. Metamodeling Metamodeling process follows the workflow as described in Section 3.3. Testing data set containing 100 LHS sampling points are used to validate the LARS-based PCE metamodel. The total number of training data starts with 50, then keeps increasing by 50 until the NRMSE value is around 2%. A detailed parametric study of NRMSE with respect to the number of training data is shown in Fig. 7, based on which the PCE metamodel using 250 training points fulfills the accuracy level and is used for optimum design

under uncertainty. This PCE metamodel has a degree of 5 after the adaptive-basis process. Since only uniform distributions are considered, Legendre orthogonal basis is used.

58

86 87 88 89

4.1.3. Results Deterministic SPO and MPO designs are completed using pyOpt directly on physics-based simulation model, since the statistical moments are not needed. While the optimum design under uncertainty requires Monte Carlo evaluations to estimate the metric of interest, i.e. mean and standard deviation of drag coefficient for standard methods and CE for the proposed UT-based method. So the design under uncertainty is completed on the LARS-based PCE metamodel driven by SQP in MATLAB. Details of key information are provided in Table 3. Results in Table 3 show that deterministic optimization design is converging faster with fewer number of iteration steps because adjoint is provided as extra physics by SU2 solver. The pyOpt SLSQP algorithm reaches the vicinity of the optimal design after six to seven direct CFD evaluations and needs another six to seven evaluations to reach the final design because of the tight convergence criteria (1e−6). MPO design costs 96 CFD evaluations in total because six evaluations are needed for one objective function calculation. Following this tendency, optimum design under uncertainty would be computationally impractical, therefore, the LARS-based PCE metamodel is introduced in this work to keep the computational costs affordable. Comparison among deterministic optimizations and robust designs is shown in Fig. 8 and Table 4. It can be seen that the SPO has the best performance the corresponded specific condition parameters. The MPO has lowest mean drag coefficient of 96.5 cts, but it also has a relatively high standard deviation of 1.4 cts. The standard robust formulation (μ(C d (x)) + σ (C d (x))) reduced the mean and standard deviation to 97.7 cts and 2.5 cts, respectively. The high-risk-aversion utility theory reduced the mean and standard deviation to 97.8 cts and 0.6 cts, respectively. All the optimized

90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122

56 57

84 85

19 20

83

Table 3 Comparison of optimization process, Case 1.

123 124

59

Method

Optimization iterations

Direct CFD evaluation

Adjoint evaluation

125

60

1 1 2 2 16 16

12 14 14 16 × 6 = 96 250a 250a

1 2 3 3 N/A N/A

126

64

SPO M ∞ = 0.700 SPO M ∞ = 0.725 SPO M ∞ = 0.750 MPO Method Standard Method (μ + σ ) Utility Theory (high risk aversion)

65

a

61 62 63

66

The number of training points to construct PCE metamodel, optimization does not require any more model evaluations.

127 128 129 130 131 132

JID:AESCTE

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.10 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

10

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20 21

Fig. 8. Baseline and optimized designs of Case 1: (a) shapes; (b) drag divergence plot. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

24 25 26 27 28 29 30 31 32 33 34

89

Table 4 Results from various methods, for Case 1.

90

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

Initial shape SPO M ∞ = 0.700 SPO M ∞ = 0.725 SPO M ∞ = 0.750 MPO Method Standard Method (μ + σ ) Utility Theory (high risk aversion)

105.4 105.4 100.3 99.6 96.5 97.7 97.8

12.7 12.7 6.6 0.7 1.4 2.5 0.6

125.0 125.0 111.3 100.3 98.6 101.6 98.4

94.2 94.2 94.6 98.8 94.4 95.4 96.8

N/A 0 −1.9E−3 −2.2E−3 −4.5E−3 −5.5E−6 −7.9E−6

N/A 0 −2.17E−4 0 0 −2.4E−4 −1.7E−5

∗ All statistics of C have the units of drag counts, and statistics of SPO and MPD methods are approximated by d

six Mach numbers shown in Fig. 7(b).

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

91 92 93 94 95 96 97 98 99 100 101

35 36

87 88

22 23

86

designs have active inequality constraints, as well as the active equality constraint. In terms of risk preferences, we ranged the utility theory from high risk loving to high risk aversion as shown in Fig. 9 and Table 5. With the utility function mentioned in (8), α is set as −1,000, −100, 5,000 and 50,000 for high risk loving, low risk loving, low risk aversion and high risk aversion, respectively. Risk neutral is 2 ∗ C d in this work. It can be seen when varying risk preference from high risk loving to high risk aversion, the standard deviation of drag coefficient keeps going down from 6.2 cts to 0.6 cts. The mean value also goes down from 99.8 cts to 97.7 cts, then remains at a low level since risk neutral. This shows that utility theory can increasing the robustness of the optimal design while keeping mean drag coefficient in a low level. In terms of increasing weighting factors on the standard deviation, results of WS formulation are shown in Fig. 10 and Table 6. It can be seen when increasing weighting factors on standard deviation, the WS method reduced the standard deviation from 2.4 cts to 0.3 cts. However, mean drag coefficient keeps going up from

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119

Fig. 9. Results comparison among risk preferences, varying from high risk loving, low risk loving, risk neutral, to low risk aversion and high risk aversion.

120 121 122 123

Table 5 Case 1 results from risk preferences.

124

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

Initial shape High risk loving Low risk loving Risk neutral Low risk aversion High risk aversion

105.4 99.8 99.0 97.9 97.7 97.8

12.7 6.2 4.8 4.0 2.6 0.6

125.0 109.8 106.9 104.5 102.0 98.4

94.2 94.6 94.9 94.6 95.2 96.8

N/A −5.0E−6 −1.2E−4 −5.9E−6 −5.5E−6 −7.9E−6

N/A −1.2E−3 −1.7E−3 −3.3E−5 −3.1E−5 −1.7E−5

∗ All statistics of C have the units of drag counts. d

125 126 127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.11 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

11

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83 84

18 19 20

Fig. 10. Results comparison among weighted-sum methods, varying the weighting factor on standard deviation from 1, 10, 50, 100 to 1000.

Fig. 11. Results comparison among σ -constrained methods, varying the constraint value on standard deviation from 10 cts, 1 cts, 0.5 cts, to 0.36 cts.

23

Table 6 Case 1 results from weighted-sum methods.

Table 7 Case 1 results from

24

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

25

Initial shape

105.4 97.7 98.3 98.5 98.6 99.3

12.7 2.4 0.4 0.3 0.3 0.3

125.0 101.6 98.6 98.9 98.9 99.5

94.2 95.4 97.7 98.1 98.2 98.8

N/A −5.5E−6 −2.8E−5 −2.8E−5 −2.8E−5 −5.01E−4

N/A −3.4E−5 −1.9E−4 −1.9E−4 −1.9E−4 −1.1E−3

26 27 28 29 30

μ+σ μ + 10σ μ + 50σ μ + 100σ μ + 1000σ

∗ All statistics of C have the units of drag counts. d

88 89

σ -constrained methods.

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

Initial shape σ < 10 cts σ < 1 cts σ < 0.5 cts σ < 0.36 cts

105.4 97.5 97.6 97.7 98.9

12.7 4.0 1.0 0.5 0.4

125.0 102.5 101.6 98.5 98.9

94.2 94.6 95.4 96.6 98.1

N/A −5.2E−6 −5.9E−6 −4.7E−6 −2.4E−6

−3.6E−5 −3.4E−5 −3.3E−5 −3.0E−5

N/A

∗ All statistics of C have the units of drag counts. d

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

97.7 cts to 99.3 cts. And the maximum drag coefficient goes up to 99.53 cts. Besides, it is hard to decide the value of weighting factor. In terms of adding constraints on the standard deviation, results are shown in Fig. 11 and Table 7. When decreasing the constraint value on standard deviation, the standard deviation was reduced from 4.0 cts to 0.3 cts. However, mean drag coefficient increased from 97.54 cts to 98.85 cts, due to insufficient guarantee on μ. Moreover, the contour plots of Mach number at various freestream Mach numbers within the uncertainty range are given in Fig. 12 to get more insights on the optimized results from high risk aversion-based utility theory. It is obvious that at very low free-stream Mach number, the optimized shape has a very weak shock, which makes the drag coefficient slightly higher than the initial shape RAE 2822. However, with the increase of free-stream Mach number, the strength of the weak shock and drag coefficient almost remain unchanged, compared with the emergence of strong shock at high free-stream Mach numbers of RAE 2822. Similarly, we can see the same tendency through the plots of pressure distributions, which are shown in Figs. 13–15. Based on the plots of pressure distribution, the utility theory keeps the weak shock in an even lower strength than those from standard robust design method.

54 55 56

4.2. Case 2: RAE 2822 airfoil in transonic viscous flow with uncertain Mach and uncertain lift coefficient

57 58 59 60 61 62 63 64 65 66

90 91 92 93 94 95 96 97

31 32

86 87

21 22

85

4.2.1. Problem formulation The original problem setup of this case is the same as (34), except that both C l and M ∞ are set as uncertain parameters. In this case, the two types of deterministic optimizations: single-point optimization (SPO) is run at C l∗ = 0.5 and M ∞ = 0.725, multi-point optimization (MPO) is run at nine pairs of C l∗ and Mach, namely, (0.4, 0.700), (0.5, 0.700), (0.6, 0.700), (0.4, 0.725), (0.5, 0.725), (0.6, 0.725), (0.4, 0.750), (0.5, 0.750) and (0.6, 0.750). In addition, standard robust design formulations and the UT-based robust design

are also applied here, following: M ∞ ∼ U (0.70, 0.75) and C l∗ ∼ U (0.4, 0.6).

98 99 100 101

4.2.2. Metamodeling Following the same process as IV.A.2, another LARS-based PCE metamodel is constructed and validated as shown in Fig. 16. The metamodel with 250 training points has a NRMSE of 2.2%, so this metamodel is selected for optimization under uncertainty. The adaptive-basis scheme selects a 5th degree PCE as the best metamodel that can be obtained from training data set. And the PCE metamodel still only needs Legendre orthogonal bases since all variables still follow Uniform distributions.

102 103 104 105 106 107 108 109 110 111

4.2.3. Results Similar optimization approaches as Case 1 are implemented. Deterministic SPO and MPO designs are completed using pyOpt directly on physics-based simulation model, while optimum design under uncertainty is completed based on the constructed PCE metamodel. Details of key information are provided in Table 8. The deterministic optimization designs convergence faster because of the adjoint information provided by SU2, but still takes 12 to 14 CFD evaluations because of the tight convergence criteria (1e−6). Comparison between deterministic optimization methods and robust design methods are shown in Fig. 17–19 and Table 9. The standard method is used as μ(C d (x)) + σ (C d (x)). It can be seen that it outperforms initial shape, SPO and MPO on both mean (97.8 cts) and standard deviation (4.8 cts) of drag coefficient within the uncertain design space. However, the difference between maximum and minimum drag coefficient is 19.9 cts. The UT-based robust design provided a mean and standard deviation of 98.4 cts and 4.1 cts, respectively. In addition, the difference between maximum and minimum drag coefficient is 14.9 cts. All the optimized designs have active inequality constraints, as well as the active equality constraint.

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

12

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.12 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55 56 57

Fig. 12. Mach contours for RAE2822 and Optimized (Utility Theory) airfoil: (a) RAE2822, M ∞ = 0.70; (b) Optimized, M ∞ = 0.70; (c) RAE2822, M ∞ = 0.725; (d) Optimized, M ∞ = 0.725; (e) RAE2822, M ∞ = 0.75; (f) Optimized, M ∞ = 0.75.

60 61 62 63 64 65 66

122 123 124

58 59

121

In terms of risk preference, we varied α in utility function as −1,000, −100, 500, 10,000 for high risk loving, low risk loving, low risk aversion and high risk aversion, respectively. Risk neutral still comes from U = 2∗ C d . The results are shown in Table 10. It can be seen when increasing the aversion on risk, the standard deviation was reduced from 8.1 cts to 4.1 cts, and mean drag coefficient can still be kept in a low level. In the meantime, maximum drag coefficient keeps decreasing down to 109.3 cts, and minimum drag

coefficient keeps going up until 94.4 cts, which is a good strategy to increase robustness. In terms of increasing weighting factors on standard deviation, Table 11 shows that optimal result reduced the standard deviation from 4.83 cts to 2.88 cts. However, the mean drag coefficient increased from 97.8 cts to 100.0 cts. Besides, the difference between maximum and minimum drag coefficient is still 17. 4 cts in the best case.

125 126 127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.13 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

13

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82 83

17 18 19

Fig. 16. Metamodel validation of Case 2.

Fig. 13. Case 1 pressure distributions at M ∞ = 0.70 for RAE2822, Standard Method, and Utility Theory.

84 85

Table 8 Comparison of optimization process, Case 2.

20 21 22

25 26 27

87

Method

Optimization iterations

Direct CFD evaluation

Adjoint evaluation

88

SPO Method MPO Method Standard Method (μ + σ ) Utility Theory (high risk aversion)

1 1 32 30

14 12 × 5 = 60 250a 250a

2 1 N/A N/A

90

23 24

86

a

The number of training points to construct PCE metamodel, optimization does not require any more model evaluations.

28 29

89 91 92 93 94 95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106 107

41 42 43

Fig. 14. Case 1 pressure distributions at M ∞ = 0.725 for RAE2822, Standard Method, and Utility Theory.

108 109

44

110

45

111

46

112

47

113 114

48

Fig. 17. Shape comparison of baseline, robust design with standard method and utility theory.

49 50

In terms of decreasing constraint value on standard deviation, Table 12 shows that the standard deviation dropped from 5.8 cts to 2.79 cts. However, the mean drag coefficient increased from 96.6 cts to 100.2 cts. And the difference between maximum and minimum drag coefficient is still 18.0 cts in the best case.

52 53 54 55 56

118 119 120 121 122 123

57

5. Conclusion

58

124 125

59 60 61 62 63 64 66

116 117

51

65

115

Fig. 15. Case 1 pressure distributions at M ∞ = 0.75 for RAE2822, Standard Method, and Utility Theory.

The paper introduces utility theory for the objective function formulation of optimum aerodynamic shape design under uncertainty. The advantage of the proposed approach over the standard aerodynamic design under uncertainty formulation is that the objective function can be written without prescribing weighting coefficients of the different statistical quantities, introducing additional constraints, or estimating the statistical moments of objective func-

126 127 128 129 130 131 132

JID:AESCTE

14

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.14 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37 38

103

Fig. 18. C d contours versus C l∗ and Mach number: (a) RAE 2822; (b) MPO; (c) Standard Method; (d) Robust Design with Utility Theory.

104 105

39

Table 9 Case 2 results from various methods.

40 41

106 107

42

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

108

43

Initial shape SPO Method MPO Method Standard Method Utility Theory

102.6 102.3 101.8 97.8 98.4

10.7 7.5 5.7 4.8 4.1

163.3 135.4 118.3 123.4 109.3

91.4 91.6 92.3 93.5 94.4

N/A −2.53E4 −1.9E−3 −6.1e−4 −4.8E−6

N/A −5.2e−5 −8.8E−3 −5.4E−3 −6.5E−3

109

44 45 46 47

∗ All statistics of C have the units of drag counts. d

48

110 111 112 113 114 115

49

116

Table 10 Case 2 results from different risk preferences.

50 51

117

52

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

118

53

Initial shape High risk loving Low risk loving Risk neutral Low risk aversion High risk aversion

102.6 101.1 97.8 96.6 97.8 98.4

10.7 8.1 5.9 5.9 4.5 4.1

163.3 147.6 129.7 121.6 113.6 109.3

91.4 93.5 92.4 91.6 93.6 94.4

N/A −4.0E−3 −8.3E−4 −6.5E−6 −5.4E−6 −4.8E−6

N/A −4.4E−3 −2.9E−3 −8.6E−6 −6.5E−3 −6.5E−3

119

54 55 56 57 58

∗ All statistics of C have the units of drag counts. d

59 60

Fig. 19. Comparison of Cd contours between RAE288 and Utility Theory.

61 62 63 64 65 66

tion. The proposed approach, therefore, overcomes the difficulty of choosing the weighting coefficients. The proposed approach is demonstrated using transonic airfoil shape optimization cases and compared with deterministic single- and multi-point design, and standard robust design. Results show that the proposed method

120 121 122 123 124 125 126

with high-risk aversion formulation outperforms the deterministic methods in the sense of yielding design that are insensitive to variations in the uncertain problem parameters. Compared to the standard aerodynamic design under uncertainty method, the proposed approach yields design of comparable efficiency but with significantly lower standard deviation, i.e. the designs obtained

127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.15 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8 9 10

Table 11 Case 2 results from weighted-sum methods. Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

Initial shape

102.6 97.8 99.6 99.9 100.0 100.0

10.7 4.8 3.0 2.9 2.9 2.9

163.3 123.4 109.7 110.5 110.7 110.8

91.4 93.5 93.3 93.4 93.4 93.4

N/A −6.1E−4 −1.5E−6 −1.9E−6 0 0

N/A −5.4E−3 −5.2E−3 −5.2E−3 −5.2E−3 −5.2E−3

μ+σ μ + 10σ μ + 50σ μ + 100σ μ + 1000σ

∗ All statistics of C have the units of drag counts. d

11 12 13

Table 12 Case 2 results from

σ -constrained method.

14

Method

μ(C d )

σ (C d )

max(C d )

min(C d )

g1

g2

15

Initial shape σ < 10 cts σ < 5 cts σ < 3 cts σ < 2.8 cts

102.6 96.6 96.2 98.7 100.2

10.7 5.8 5.0 3.0 2.8

163.3 120.1 119.8 110.4 111.0

91.4 91.6 92.6 93.0 93.0

N/A −8.3E−4 −8.2E−4 −6.6E−6 −6.6E−6

N/A −3.4E−3 −3.7E−3 −5.5E−3 −5.5E−3

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

∗ All statistics of C have the units of drag counts. d

through the proposed approach are less sensitive to the parameter variations than the designs obtained with standard approach. And the optimum designs with constraints on standard deviation of drag coefficient cannot guarantee the mean of drag coefficient well. Future work will investigate problems of higher complexity, i.e. aerodynamic cases requiring models of higher simulation degrees of freedom and larger parameter spaces. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] M. Giles, M. Drela, W.T. Thompkins Jr., Two-dimensional transonic aerodynamic design method, AIAA J. 25 (9) (September 1987) 1199–1206. [2] K.D. Lee, S. Evi, Aerodynamic design via optimization, J. Aircr. 29 (6) (1992) 1012–1019. [3] A. Viswanath, A.I.J. Forrester, A.J. Keane, Dimension reduction for aerodynamic design optimization, AIAA J. 49 (6) (June 2011). [4] S. Kim, K. Hosseini, K. Leoviriyakit, A. Jameson, Enhancement of a class of adjoint design methods via optimization of parameters, AIAA J. 48 (6) (June 2010) 1072–1076. [5] S. Nadarajah, A. Jameson, Optimum shape design for unsteady flows with time-accurate continuous and discrete adjoint methods, AIAA J. 45 (7) (July 2007) 1478–1491. [6] S. Kim, J. Alonso, A. Jameson, Multi-element high-lift configuration design optimization using viscous continuous adjoint method, J. Aircr. 41 (5) (2004) 1082–1097. [7] A. Jameson, Aerodynamic design via control theory, J. Sci. Comput. 3 (1988) 223–260. [8] Z. Lyu, J. Martins, Aerodynamic shape optimization of an adaptive morphing trailing-edge wing, J. Aircr. 52 (6) (2015) 1951–1970. [9] Z. Lyu, J. Martins, Aerodynamic design optimization studies of a blendedwing-body aircraft, J. Aircr. 51 (5) (2014) 1604–1617. [10] F. Grasso, Usage of numerical optimization in wind turbine airfoil design, J. Aircr. 48 (1) (2011) 248–255. [11] F. Grasso, Design and optimization of tidal turbine airfoil, J. Aircr. 49 (2) (2012) 636–643. [12] F. Grasso, Development of thick airfoils for wind turbines, J. Aircr. 50 (3) (2013) 975–981. [13] D.A. Masters, D.J. Poole, N.J. Taylor, T.C.S. Rendall, C.B. Allen, Influence of shape parameterization on a benchmark aerodynamic optimization problem, J. Aircr. 54 (6) (2017) 2242–2256. [14] D.J. Poole, C.B. Allen, T.C.S. Rendall, Global optimization of wing aerodynamic optimization case exhibiting multimodality, J. Aircr. 55 (4) (2018) 1576–1591. [15] S. Maraniello, R. Palacios, Optimal rolling maneuvers with very flexible wings, AIAA J. 55 (9) (2017) 2964–2979. [16] Y. Wang, A. Wynn, R. Palacios, Nonlinear aeroelastic control of very flexibile aircraft using model updating, J. Aircr. 55 (4) (2018) 1551–1563.

15

[17] L. Fowler, J. Rogers, Airframe performance optimization of guided porjectiles using design of experiments, J. Spacecr. Rockets 52 (6) (2015) 1603–1613. [18] A. Leonard, J. Rogers, J. Sahu, Aerodynamic optimization of micorspoiler actuator for guided projectiles, J. Spacecr. Rockets 54 (6) (2017) 1216–1227. [19] J.J. Reuther, A. Jameson, J.J. Alonso, M.J. Rimlinger, D. Saunders, Constrained multipoint aerodynamic shape optimization using an adjoint formulation and parallel computers, part 1, J. Aircr. 36 (1) (1999) 51–60. [20] F. Gallard, M. Meaux, M. Montagnac, B. Mohammadi, Aerodynamic aircraft design for mission performance by multipoint optimization, in: AIAA Computational Fluid Dynamics Conference, 2013. [21] T.M. Leung, D.W. Zingg, Single- and multi-point aerodynamic shape optimization using a parallel Newton-Krylov approach, in: 19th AIAA Computational Fluid Dynamics, San Antonio, TX, 2009. [22] T.M. Leung, D.W. Zingg, Aerodynamic shape optimization of wings using a parallel Newton-Krylov approach, AIAA J. 50 (3) (2012) 540–550. [23] H. Buckley, B. Zhou, D. Zingg, Airfoil optimization using practical aerodynamic design requirements, J. Aircr. 47 (5) (2010) 1707–1719. [24] H. Buckley, D. Zingg, Approach to aerodynamic design through numerical optimization, AIAA J. 51 (8) (2013) 1972–1981. [25] M.H.A. Madsen, F. Zahle, N.N. Sorensen, J.R.R.A. Martins, Multipoint highfidelity CFD-based aerodynamic shape optimization of a 10 MW wind turbine, Wind Energy Sci. (2018), submited for publication. [26] J.C. Helton, J.D. Johnson, W.L. Oberkampf, C.J. Sallaberry, Representation of Analysis Results Involving Aleatory and Epistemic Uncertainty, SANDIA Report, 2008. [27] J. Mullins, Y. Ling, S. Mahadevan, L. Sun, A. Strachan, Separation of aleatory and epistemic uncertainty in probabilistic model validation, Reliab. Eng. Syst. Saf. 47 (2016) 49–59. [28] G.W. Parry, P.W. Winter, Characterization and evaluation of uncertainty in probabilistic risk analysis, Nucl. Saf. 22 (1) (1981) 28–42. [29] L. Jaulin, M. Kieffer, O. Didrit, E. Walter, Applied Interval Analysis, SpringerVerlag, New York, NY, 2001. [30] D. Dubois, H. Prade, Possibility Theory: An Approach to Computerized Processing of Uncertainty, Plenum Press, New York, NY, 1986. [31] A. O’Hagan, J.E. Oakley, Probability is perfect, but we can’t elicit it perfectly, Reliab. Eng. Syst. Saf. 85 (1–3) (2004) 239–248. [32] G. Shafer, A Mathematical Theory of Evidence, Princeton University Press, Princeton, NJ, 1976. [33] J.C. Helton, R.J. Breeding, Calculation of reactor accident safety goals, Reliab. Eng. Syst. Saf. 39 (2) (1993) 129–158. [34] S.R. Baker, N. Bloom, S.J. Davis, Measuring economic policy uncertainty, Q. J. Econ. 131 (4) (2016) 1593–1636. [35] J.W. Dyer, Aerospace Safety Advisory Panel, NASA Annual Report, 2012. [36] G. Taguchi, M.S. Phadke, Quality engineering through design optimization, in: IEEE Global Telecommunications Conference, GLOBECOM ’84, Atlanta, GA, November 26–29, 1984. [37] G. Taguchi, S. Chowdhury, Y. Wu, Taguchi’s Quality Engineering Handbook, John Wiley & Sons, Hoboken, N.J., 2005. [38] A. Chaudhuri, G. Waycaster, N. Price, T. Matsumura, R.T. Haftka, NASA uncertainty quantification challenge: an optimization-based methodology and validation, J. Aerosp. Inform. Syst. 12 (1) (2015) 10–34. [39] S. Han, M. Tao, U. Topcu, H. Owhadi, R.M. Murray, Convex optimal uncertainty quantification, SIAM J. Optim. 5 (3) (2015) 1368–1387. [40] A. Jafarsalehi, H.R. Fazeley, M. Mirshams, Spacecraft mission design optimization under uncertainty, J. Mech. Eng. Sci. 230 (16) (2016) 2872–2887. [41] G.C. Calafiore, F. Dabbene, Optimization under uncertainty with application to design of truss structures, Struct. Multidiscip. Optim. 35 (2008) 189–200. [42] A.T. Eshghi, S. Lee, M.K. Sadoughi, C. Hu, Y.C. Kim, J.H. Seo, Design optimization under uncertainty and speed variability for a piezoelectric energy harvester powering a tire pressure monitoring sensor, Smart Mater. Struct. 26 (2017) 1–18. [43] X. Zhang, S. Mahadevan, S. Sankararaman, K. Goebel, Resilience-based network design under uncertainty, Reliab. Eng. Syst. Saf. 169 (2018) 364–379. [44] Y. Aoues, A. Chateaueuf, D. Lemosse, A. El-Hami, Optimal design under uncertainty of reinforced concrete structures using system reliability approach, Int. J. Uncertain. Quantificat. 3 (6) (2013) 487–498. [45] Y. Yi, W. Li, M. Xiao, L. Gao, A set strategy approach for multidisciplinary robust design optimization under interval uncertainty, Adv. Mech. Eng. 11 (1) (2019) 1–17. [46] D.I. Papadimitriou, C. Papadimitriou, Aerodynamic shape optimization for minimum robust drag and lift reliability constraint, Aerosp. Sci. Technol. 55 (2016) 24–33. [47] A. Vuruskan, S. Hosder, Investigation of the impact of turbulence models on robust aerodynamic shape optimization, in: AIAA Aerospace Sciences Meeting, AIAA SciTech Forum, 2018. [48] A. Vuruskan, S. Hosder, Impact of turbulence models and shape parameterization on robust aerodynamic shape optimization, J. Aircr. 56 (3) (2019) 1–17. [49] R.P. Liem, J.R.R.A. Martins, G.K.W. Kenway, Expected drag minimization for aerodynamic design optimization based on aircraft operational data, Aerosp. Sci. Technol. 63 (2017) 344–362.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.16 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

[50] L. Cook, J.P. Jarrett, K. Willcox, Extending horsetail matching for optimization under probabilistic, interval, and mixed uncertainties, AIAA J. 56 (2) (2018) 849–861. [51] L. Cook, J.P. Jarrett, Optimization using multiple dominance criteria for aerospace design under uncertainty, AIAA J. 56 (12) (2018) 4956–4976. [52] M. Caboni, E. Minisci, A. Riccardi, Aerodynamic design optimization of wind turbine airfoils under aleatory and epistemic uncertainty, J. Phys. Conf. Ser. 1037 (2018) 1–10. [53] D.R. Abhiram, R. Ganguli, D. Harursampath, Robust design of small-scaled unmanned helicopter for hover performance using Taguchi method, in: 55th AIAA Aerospace Sciences Meeting, 2017, pp. 1–25. [54] D.R. Abhiram, R. Ganguli, D. Harursampath, Robust design of small unmanned helicopter for hover performance using Taguchi method, J. Aircr. 55 (4) (2018) 1745–1752. [55] A. Stoebner, S. Mahadevan, Robustness in reliability-based design, in: Proceedings of the 41st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2000. [56] S. Wang, Q. Li, G. Savage, Reliability-based robust design optimization of structures considering uncertainty in design variables, Math. Probl. Eng. 2015 (2015). [57] T. Zang, M. Hemsch, M. Hilburger, S. Kenny, J. Luckring, P. Maghami, S. Padula, W. Stroud, Needs and Opportunities for Uncertainty-Based Multidisciplinary Design Methods for Aerospace Vehicles, Technical Report, National Aeronautics and Space Administration, 2002. [58] Z. Mourelatos, An Efficient Unified Approach for Reliability and Robustness in Engineering Design, 2004. [59] I. Das, J. Dennis, A closer look at drawbacks of minimizing sums of objectives for Pareto set generation in multicriteria optimization problems, Struct. Optim. 14 (1) (1997) 63–69. [60] N. Vidanovic, B. Rasuo, G. Kastratovic, S. Maksimovic, D. Curcic, M. Samardzic, Aerodynamic-structural missile fin optimization, Aerosp. Sci. Technol. 65 (2017) 26–45, https://doi.org/10.1016/j.ast.2017.02.010. [61] G. Ocokojic, B. Rasuo, A. Bengin, Aerodynamic shape optimization of guided missile based on wind tunnel testing and CFD simulation, Therm. Sci. 21 (3) (2017) 1543–1554, https://doi.org/10.2298/TSCI150515184O. [62] W. Chen, M. Wiecek, J. Zhang, Quality Utility – A Compromise Programming Approach to Robust Design, 1998. [63] S. Jeong, G. Park, Reliability-based robust design optimization using the probabilistic robustness index and the enhanced single loop single vector approach, in: 10th World Congress on Structural and Multidisciplinary Optimization, 2013. [64] P.C. Fishburn, Utility Theory for Decision Making, John Wiley & Sons, Inc., Research Analysis Corporation, 1970. [65] D. Dalalah, W. Al-Rawabdeh, Benchmarking the utility theory: a data envelopment approach, Benchmarking 24 (2) (2017) 318–340. [66] G. Blatman, B. Sudret, An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probab. Eng. Mech. 25 (2) (2010) 183–197. [67] G. Blatman, B. Sudret, Adaptive sparse polynomial chaos expansion based on least angle regression, J. Comput. Phys. 230 (6) (2011) 2345–2367. [68] G. Blatman, Adaptive Sparse Polynomial Chaos Expansion for Uncertainty Propagation and Sensitivity Analysis, Doctor of Philosophy Thesis, Blaise Pascal University, 2009. [69] T. Glazier, B. Schmerl, J. Camara, D. Garlan, Utility Theory for Self-Adaptive Systems, Research Report, 2017. [70] A. Abbas, A. Cadenbach, On the use of utility theory in engineering design, IEEE Syst. J. 12 (2) (2018) 1129–1138. [71] J. Ren, A. Thelen, A. Amrit, X. Du, L. Leifson, Y. Tesfahunegn, S. Koziel, Application of multifidelity optimization techniques to benchmark aerodynamic design problems, AIAA-2016-1542, in: AIAA Aerospace Sciences Meeting, San Diego, CA, Jan. 4–8, 2016. [72] S.E. Cliff, J.J. Reuther, D.A. Saunders, R.M. Hicks, Single-point and multipoint aerodynamic shape optimization of high-speed civil transport, J. Aircr. 38 (6) (2001) 997–1005. [73] Z. Lyu, G. Kenway, J. Martins, Aerodynamic shape optimization investigations of the common research model wing benchmark, AIAA J. 53 (4) (2015) 968–985. [74] A. Jameson, K. Leoviriyakit, S. Shankaran, Multi-point aero-structural optimization of wings including planform variations, in: 45th Aerospace Sciences Meeting and Exhibit, Reno, Nevada, Jan. 8–11, 2007. [75] G. Kenway, J. Martins, Multi-point high-fidelity aerostructural optimization of a transport aircraft configuration, J. Aircr. 55 (1) (2014) 144–160. [76] Y. Zhang, S. Hosder, L. Leifsson, S. Koziel, Robust airfoil optimization under inherent and model-form uncertainties using stochastic expansions, AIAA2012-0056, in: AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Nashville, TN, Jan. 9–12, 2012. [77] M. Shahbaz, Z. Han, W. Song, M. Aizud, Surrogate-based robust design optimization of airfoil using inexpensive Monte Carlo method, in: IEEE, Applied Sciences and Technology (IBCAST), 13th International Bhurban Conference, Jan. 12–16, 2016.

[78] X. Li, P. Nair, Z. Zhang, L. Gao, C. Gao, Aircraft robust trajectory optimization using nonintrusive polynomial chaos, J. Aircr. 51 (5) (2014) 1592–1603. [79] C. Brian, Is expected utility theory normative for medical decision making, Med. Decis. Mak. 16 (1) (1996) 1–6. [80] L. Jonathan, Choice under uncertainty, Stanford Lecture Notes, October 2006. [81] G. Farus, P. Wolfgang, Expected uncertainty utility theory, Econometrica 82 (1) (2014) 1–39. [82] I. Moscati, How economists came to accept expected utility theory: the case of Samuelson and Savage, J. Econ. Perspect. 30 (2) (2016) 219–236. [83] E. Diecidue, U. Schmidt, P. Wakker, The utility of gambling reconsidered, J. Risk Uncertain. 29 (3) (2004) 241–259. [84] K. Kontek, Decision utility theory: back to von Neumann, Morgenstern, and Markowitz, in: Social Science Research Network (SSRN), December 2010. [85] H. Stefansson, R. Bradley, What is risk aversion, Br. J. Philos. Sci. (2017). [86] M. Bombardini, F. Trebbi, Risk aversion and expected utility theory: an experiment with large and small stakes, J. Eur. Econ. Assoc. 10 (6) (2012) 1348–1399. [87] N.J. Devlin, K.K. Shah, B.J. Mulhern, K. Pantiri, B. Van Hout, A new method for valuing health: directly eliciting personal utility functions, Eur. J. Health Econ. (2018) 1–14. [88] M. Rabin, Risk aversion and expected-utility theory: a calibration theorem, Econometrica 68 (5) (2000) 1281–1292. [89] R. Chetty, A. Szeidl, Consumption Commitments and Risk Preferences, The National Bureau of Economic Research (NBER) Working Paper, No. 12467, 2006. [90] G. Charness, U. Gneezy, A. Imas, Experimental methods: eliciting risk preferences, J. Econ. Behav. Organ. 87 (2013) 43–51. [91] C. Gollier, J.K. Hammitt, N. Treich, Risk and choice: a research saga, J. Risk Uncertain. 47 (2) (2013) 129–145. [92] R. Pettigrew, Risk, rationality and expected utility theory, Can. J. Philos. 45 (5–6) (2015) 798–826. [93] H. Kannan, C. Bloebaum, B. Mesmer, Incorporation of risk preferences in a value-based systems engineering framework for satellite system, in: AIAA SciTech, 18th AIAA Non-Deterministic Approaches Conference, San Diego, California, Jan. 4–8, 2016. [94] L. Leifsson, S. Koziel, P. Kurgan, Automated low-fidelity model setup for surrogate-based aerodynamic optimization, in: Solving Computationally Expensive Engineering Problems, vol. 97, 2014, pp. 87–111. [95] J.A. Samareh, Survey of shape parameterization techniques for high-fidelity multidisciplinary shape optimization, AIAA J. 39 (5) (2001) 877–884. [96] J.C. Tannehill, D.A. Anderson, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, 2nd edition, Taylor & France, 1997. [97] C. Nguyen, Turbulence modeling incompressible case 1.1 Reynolds-Averaged Navier-Stokes equations, MIT Lecture Notes, 1980. [98] F. Palacios, M. Colonno, Stanford university unstructured (SU2): an open source integrated computational environment for multi-physics simulation and design, in: AIAA 2013-0287, 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Grapevine, Texas, Jan. 7–10, Jan 2013. [99] F. Palacios, T. Economon, Stanford university unstructured (SU2): open-source analysis and design technology for turbulent flows, in: AIAA 2014-0243, 52nd AIAA Aerospace Sciences, National Harbor, Madrid, Jan. 13–17, 2014. [100] T. Economon, F. Palacios, Towards High-Performance Optimizations of the Unstructured Open-Source SU2 Suite, in: AIAA 2015-1949, AIAA Infotech@Aerospace, Kissimmee, Florida, Jan. 5–9, 2015. [101] M. Nemec, D. Zingg, T. Pulliam, Multipoint and multi-objective aerodynamic shape optimization, AIAA J. 42 (6) (2004) 1057–1065. [102] M. Drela, Pros & cons of airfoil optimization, in: Frontiers of Computational Fluid Dynamics, 1998, pp. 363–381. [103] M. Walter, A Stochastic Expansion-Based Approach for Design Under Uncertainty, PhD Thesis, School of Aerospace Engineering, Georgia Institute of Technology, May 2013. [104] M. Eldred, C. Webster, P. Constantine, Design under uncertainty employing stochastic expansion methods, in: AIAA 2008-6001, AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria, British Columbia, Canada, Sep. 10–12, 2008. [105] L. Leifsson, S. Koziel, Y. Zhang, S. Hosder, Low-cost robust airfoil optimization by variable-fidelity models and stochastic expansions, in: 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Grapevine, TX, Jan. 7–10, 2013. [106] H. Shah, S. Hosder, S. Koziel, Y.A. Tesfahunegn, L. Leifsson, Multi-fidelity robust aerodynamic design optimization under mixed uncertainty, Aerosp. Sci. Technol. 45 (2015) 17–29. [107] Sriram, A. Jameson, Robust optimal control using polynomial chaos and adjoints for systems with uncertain inputs, in: AIAA 2011-3069, 20th AIAA Computational Fluid Dynamics Conference, Honolulu, HI, Jun. 27–30, 2011. [108] S. Marelli, B. Sudret, UQLab: a framework for uncertainty quantification in Matlab, in: Vulnerability, Uncertainty, and Risk, 2014. [109] M. Udell, C. Horn, R. Zadeh, S. Boyd, Generalized low rank models, Found. Trends Mach. Learn. 9 (1) (2016) 1–118. [110] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, Least angle regression, Ann. Stat. 32 (2) (2004) 407–499.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105464 /FLA

[m5G; v1.261; Prn:14/10/2019; 15:48] P.17 (1-17)

X. Du, L. Leifsson / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8

[111] R.E. Perez, P.W. Jansen, J.R.R.A. Martins, PyOpt: a python-based objectiveoriented framework for nonlinear constrained optimization, Struct. Multidiscip. Optim. 45 (1) (2012) 101–118. [112] S.P. Han, Superlinearly convergent variable metric algorithm for general nonlinear programming problems, Math. Program. 11 (3) (1976) 263–282. [113] M.J.D. Powell, A fast algorithm for nonlinearly constrained optimization calculations, in: G.A. Watson (Ed.), Numerical Analysis Dundee 1977, in: Lecture Notes in Mathematics, vol. 630, Springer Verlag, Heidelberg, Berlin, New York, 1978, pp. 144–157.

17

[114] J.C. Helton, F.J. Davis, Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems, Reliab. Eng. Syst. Saf. 81 (1) (2003) 23–69. [115] MATLAB “MATLAB 2017b”, Version 2017b, The Mathworks Inc., Natick, Massachusetts, 2017. [116] A.M. DeGennaro, C.W. Rowley, L. Martinelli, Uncertainty quantification for airfoil icing using polynomial chaos expansions, J. Aircr. 52 (5) (2015) 1404–1411.

67 68 69 70 71 72 73 74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132