Preference-driven Kriging-based multiobjective optimization method with a novel multipoint infill criterion and application to airfoil shape design

JID:AESCTE AID:105555 /FLA

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

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 14 15 16 17 18 19 20

Preference-driven Kriging-based multiobjective optimization method with a novel multipoint infill criterion and application to airfoil shape design Youwei He a , Jinju Sun a,b,∗ , Peng Song a , Xuesong Wang a , Asif S. Usmani c School of Energy and Power Engineering, Xi’an Jiaotong Univ., 28 West Xianning Road, Xi’an, 710049, China b Collaborative Innovation Center for Advanced Aero-Engine (CICAAE) of China, 37 Xueyuan Road, Beijing, 100191, China c Department of Building Services Engineering Faculty of Construction and Environment, The Hong Kong Polytechnic University, Kowloon, Hong Kong

27 28

i n f o

31 32 33 34 35

81 82 83 85 86

89

a b s t r a c t

90

Article history: Received 4 July 2019 Received in revised form 14 October 2019 Accepted 10 November 2019 Available online xxxx

29 30

80

88

a r t i c l e

24 26

79

87

22

25

78

84

a

21 23

77

Keywords: Kriging-based multiobjective optimization Expected angle-penalized length improvement Multipoint infill sampling criterion Preference-driven optimization Airfoil shape optimization

36 37 38 39 40 41 42 43

Multipoint infill sampling has emerged recently as a promising tool to enhance the efficiency of the adaptive Kriging-based optimization method for computational expensive multiobjective problems. However, the number of infill samples in each iteration in the open literature is up to 10. A novel multipoint infill criterion is proposed on the basis of Expected Angle-Penalized Length Improvement (EAPLI). One of its distinct features is that it makes the number of infill samples in each iteration increase tenfold along with a fair convergence speed by using multiple reference vectors to search for multiple infill points. Its implementation involves several major elements: firstly, angle-penalized length based improvement is established for each reference vector; then maximizing EAPLI is conducted for all vectors to obtain sufficient candidate infill points; finally, multiple infill points are selected from the candidate pool according to the candidate niche counts. Experimental results on benchmark test problems show that the proposed EAPLI criterion is highly competitive in comparison with other criteria. Then the performance of EAPLI using different numbers of multiple infill points, from several to tens and even a hundred, is studied. The other distinct feature of EAPLI is, the distribution of reference vectors can be manipulated to include the designer’s preferences to perform preference-driven optimization, which can reduce searches towards solutions violating the preferences and thus further enhance the optimization efficiency. A preference-driven airfoil shape optimization method is established by integrating the EAPLI, free-form deformation method, and a flow solver. It is used to optimize the NACA0012 airfoil, by which the non-dominated solutions are fairly driven into the preferred subspaces with significant gains in both objectives of lift/drag coefficient ratio and drag coefficient simultaneously. © 2019 Published by Elsevier Masson SAS.

91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109

44

110

45

111 112

46 47

1. Introduction

48 49 50 51 52 53 54 55 56 57 58

Multi-objective optimization technique has been widely used in the computational expensive shape design of aerospace and mechanical engineering applications, such as the design of high-lift airfoil [1], heat exchanger [2], transonic axial compressor rotor [3], re-entry spacecraft [4], aerodynamic and aeroacoustic airfoil [5], and centrifugal compressor [6] etc. Various multi-objective optimization methods are used in these applications ranging from nature-inspired evolutionary algorithms (EA), gradient-based and surrogate-based algorithms. But with EA algorithms, thousands of simulation calls are often required in such applications to eval-

59 60 61 62 63 64 65 66

*

Corresponding author at: School of Energy and Power Engineering, Xi’an Jiaotong Univ., 28 West Xianning Road, Xi’an, 710049, China. E-mail address: [email protected] (J. Sun). https://doi.org/10.1016/j.ast.2019.105555 1270-9638/© 2019 Published by Elsevier Masson SAS.

uate the optimization objectives, which largely restrains the usefulness of the optimization methods in solving the practical expensive aerospace/engineering design problems. Difficulties also arise in utilizing the gradient-based methods that a number of local optimizations with multiple initial points and weights are required to obtain an approximation of the complete Pareto front (PF) but the computation of the gradients is often impossible in some cases. Quite recently, the surrogate-based optimization method with adaptive sampling method has emerged as an efficient alternative to solve the expensive multiobjective problems in aerospace engineering applications with much reduced objective evaluations [7–10]. With adaptive sampling, only a portion of samples are selected initially to establish a preliminary model, and then additional samples are successively identified by adaptive sampling, it is evaluated by the real model and used to refine the surrogate model iteratively. (Rules to guide the selection of

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

JID:AESCTE

2

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:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.2 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

additional samples are termed as infill sampling criteria and the refining sample points termed as infill points [11]). There are different kinds of surrogate models to be used such as radial basis function, support vector machine, and Kriging model, among which the Kriging model is quite popular as Kriging can approximate the real complex optimization function and also provide the approximation errors [12]. Incorporation expected improvement (EI) infill criterion into the Kriging model to make it adaptive is pioneered by Mockus et al. [13,14] and then promoted by Schonlau [15] and Jones et al. [16]. It is acknowledged that EI criterion balances exploration and exploitation which helps the optimizer to search globally and efficiently, thus it is regarded as Efficient Global Optimization (EGO) [16]. Further modifications to EI have been conducted by some other investigators. For instance, it is improved by Ponweiser et al. [17] and Kleijnen et al. [18], and extended to constrained problems [19,20], variable-fidelity problems [21], and noisy problems [22–24]. In [25–30], the EI infill criterion is extended to multipoint criteria, by which the multiple infill points are selected and used in each update iteration, this improves the searching efficiency. And EI has also been used in expensive multiobjective optimizations. In general, the extension of the single-objective EI to multiobjective EI is realized by integrating a defined multiobjective improvement function over the objective space, so that the multiobjective EI-based infill criteria is obtained. For instance, Emmerich [31] defined the expected hypervolume improvement (EHVI) and used it as a prescreening criterion in multiobjective EAs. Łaniewski-Wołłk et al. used EHVI as infill criterion [32] and compared it to other criteria in [33]. Martinez-Frutos and Herrero-Perez [34] extended EHVI to deal with constraints in multi-objective optimization. Keane [35] and Bautista [36] proposed the Expected Euclidean Distance Improvement (EEDI) and Expected Maximin Distance Improvement (EMDI) infill criteria by using the Euclidean distance and maximin distance as the improvements function respectively. But in common, EHVI, EEDI, and EMDI are of heavy computational burden, which largely compromises their usefulness. To reduce computations, some modified criteria are further proposed. Zhan et al. [37] proposed three multiobjective EI criteria on the basis of Expected Improvement Matrix (EIM), namely Euclidean Distance-based EIM (EIMe ), Maximin Distance-based EIM (EIMm ) and Hypervolume-based EIM (EIMh ) criteria. Competitive approximation features of the PF are demonstrated from the comparisons of the EIM criteria with EHVI, EEDI, and EMDI and their effectiveness in time-saving has also been verified. It is noted that the aforementioned EI-based criteria of multiobjective optimization are all single-point infill criteria, and only one infill point is selected and used to refine the surrogate model for every update iteration. By contrast, multipoint infill criteria for EI-based adaptive sampling can help to further enhance the optimization efficiency, in which multiple sample points are selected and used to refine the surrogate model, where the calculations of each infill point essentially are of parallel feature, thus they can be done in parallel so as to accelerate the optimization. Interesting work has been conducted on multiple infill-point EIbased criteria by some previous investigators. For instance, Zhang et al. [38] extended the single-objective EGO algorithm into multiobjective EGO algorithms by aggregating the multiple objectives into scalar functions being termed as MOEA/D-EGO. The EI values for each weight vector are maximized to obtain sufficient candidate points, from which some are selected as infill points. In [38], the number of multiple infill points used is up to 5. Shimoyama et al. [39] combined the EHVI with the estimated objective function values (termed as EHVI+EST) and explored the nondominated solutions (NDSs) of EHVI and M objective functions in an (M+1)-dimensional objective space. Then multiple infill points

nearest to the center of (M+1)-dimensional NDSs are selected and used to update the surrogate model. But the application of such an EHVI+EST infill criterion is restrained to 2-objective optimization problems, and it is unsuited for higher dimension problems due to the excessively heavy computational burden of EHVI. Zhan et al. [40] extended the three single-point multiobjective EIM criteria to three multipoint one, i.e. Pseudo EIM criteria (PEIM), in which the first infill point is selected at the maximum of an EIM criterion, and then multiply the EIs in the EI matrix with the influence function of the newly obtained (first) infill point to obtain the updated EIM criterion; the second point is selected by maximizing the approximated EIM criterion, and similarly other successive infill points are sequentially determined. In comparison with EHVI+EST, the PEIM criterion is less computationally intensive, however, the more the infill points are selected, the less reliable the landscape of the approximated EIs becomes and the less promising the refining points are. In [40], the number of multiple infill points used is up to 10. In a word, to solve expensive multiobjective problems efficiently, Kriging-based optimization with EI-based adaptive sampling is essential. The single infill point criterion of EI-based adaptive sampling is well-established but there is yet much to do to establish a generally-suited multipoint infill criteria for dealing with more infill points (i.e. larger than 10) at each Kriging-model update to further accelerate the optimization of expensive multiobjective problems, since multipoint infill criteria is of sound parallel feature, thus the evaluation process of infill points can be accelerated with parallel algorithms. Another important issue for practical multi-objective optimization is of preferences, that is, the designer prefers the solutions in one or several specific regions, beyond which the solutions are unwanted and have to be abandoned. However, all the aforementioned methods search the solutions approximating whole PF, where some of the solutions satisfying the preferences are accepted by the designer and the rest violating the preferences are wasted, which can largely compromise the optimization efficiency. Preference-driven approaches can make optimal searching cost-effective. Some efforts have been made by the previous investigators towards developing the preference-driven multiobjective evolutionary algorithms [41–43] and preference-based surrogate model building technique in single-objective optimization [44]. But with Kriging-based multiobjective optimization methods, the preference-driven feature of optimal searching hasn’t yet been realized and it calls for some in-depth work. The present study is aimed at solving expensive multiobjective problems efficiently through adaptive Kriging-based optimization. To that end, a novel EI-based multipoint infill criterion is proposed based on the Angle-Penalized Length (APL), namely Expected Angle-Penalized Length Improvement (EAPLI). Implementation of EAPLI involves several major elements: firstly, multiple reference vectors are used and APL-based improvement is established for each reference vector; then maximizing EAPLI is done for all reference vectors and sufficient candidate infill points are obtained; finally, the multiple infill points are selected from the candidate pool based on the relation of the niche counts for all the candidate points. The proposed EAPLI criterion is of two distinct features: with sufficient reference vectors, the number of infill points for each update can be up to tens and hundreds; by a proper manipulation of reference vectors, the EAPLI is well-suited for the non-preference optimization as well as for the preferencedriven optimization, where a basic uniform distribution of reference vectors are appropriate for the former while a tailored distribution of reference vectors according to the preferences must be used for the latter. The newly developed EAPLI is used to solve benchmark test problems, the obtained results are compared with that of previous investigators, from which the capacity of the pro-

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:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.3 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6 7 8 9 10 11

3

posed EAPLI in terms of the non-preference optimization is welldemonstrated. Besides, the performance of EAPLI using different numbers of infill points, from 3 to 100, in one update iteration are studied with the observations of competitive performance in terms of solution quality but remarkably reduced iterations. The capacity of EAPLI in preference-driven optimization is demonstrated on analytic problems as well as a practical airfoil shape optimization problem. In particular, the NACA0012 airfoil is optimized, NDSs are fairly driven into the preferred subspaces, and significant gains in both objectives of lift/drag coefficient ratio and drag coefficient are produced from the optimization.

67 68 69 70 71 72 73 74 75 76 77 78

12 13 14

2. Kriging-based optimization method for expensive multiobjective problems

79 80

15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

The novel infill criterion is developed for Kriging-based Multiobjective Optimization (KMOP). The KMOP method adopted in the present study is briefly introduced in this section. Fig. 1 presents the flowchart of the KMOP method. It starts with generating a small set of initial samples using the design of experiment (DOE) method, such as Latin hypercube sampling (LHS) [45]. Then the initial samples are evaluated using the computational expensive models, based on which, the Kriging models for each objective function are constructed. Next, the ideal and nadir points are defined and used to scale the objective space. The infill criterion is maximized to obtain the candidate infill points, from which the infill points are then selected and evaluated, and they are then added to the sample set. Such a procedure will iterate until the stopping rule is met.

30 31

2.1. Kriging model

32 33 34 35

Given p n-dimensional points x1 , . . . , x p ∈ R n and their associated responses f 1 , . . . , f p , the Kriging model approximates the unknown objective function, written as

36 37 38 39 40 41 42 43 44 45 46 47 48

y (x) = μ(x) + ε (x)

(1)

where μ(x) is a global model, and ε (x) represents a local deviation from the global model being defined as a Gaussian distribution of zero mean and the standard deviation σ . In general, a spatial correlation between the deviation terms of two points xi and x j is assumed and used. The commonly used correlation functions are the squared exponential function, the Matérn function and the exponential function [46]. The squared exponential function with additional hyperparameters in the field of design and analysis of computer experiments is adopted, written as

 

49 50

Corr

ε xi , x

 j



n   j p = exp − θk xki − xk k

53 54 55 56 57 58 59 60 61





s (x) = σˆ 2

2

T −1

1−C C

C+

(1 − 1T C−1 C)2 1T C−1 1



(3) (4)

62 63 64 65 66

83

2.2. Multiobjective expected improvement

84 85

Consider the continuous multiobjective optimization problem with only box constraints, given as

minimize

F (x) = ( f 1 (x), . . . , f m (x))T

subject to

x = (x1 , . . . , xn ) ∈ T

n

i =1 [l i , u i ]

(5)





E I (x) = y∈ Rm

  m

1 y i − yˆ i I (x) φ dy i i =1

si

where C denotes a matrix with entry C i j = Corr[ε (xi ), ε (x j )], c is an n-dimensional vector with entry c i = Corr[ε (x), ε (xi )], y = ( f 1 , f 2 , ... f p ) is the vector of the observed function values, and 1 is n-dimensional unit vectors.

89 90 92 93 94 95 96 97 98 99 100 101

i = 1, . . . , m

(6)

si

102 103

where yˆ i (x) and s2i (x) are the prediction and variance value of the concerned point x provided by the ith Kriging model. The multiobjective EI (corresponding to the expected value of the defined improvement over the current NDSs among the sample points) is defined as the integration of the improvement function I (x) over the whole objective space [37], written as



87

91

n

where −∞ < li < u i < +∞ for all i = 1, . . . , n. F : i =1 [li , u i ] → R m consists of m individual objective functions f i , i = 1, . . . , m. All the objectives are continuous, noise-free and expensive to evalun ate. i =1 [li , u i ] is called the design space and R m is the objective space. In the present KMOP, each objective is approximated by a Kriging model and the objective vector of a concerned point x can then be described as a set of m-dimensional Gaussian random variables, written as

Y i (x) ∼ N yˆ i (x), s2i (x)

86 88

104 105 106 107 108 109 110 111

(7)

(2)

where θk > 0 and pk ∈ {1, 2}. As a result, the Kriging model is determined by 2n + 2 parameters: μ, σ 2 , θ1 , . . . , θn , p 1 , ... pn . They are obtained by maximizing the likelihood function. The prediction and uncertainty of the unvisited point x by Kriging model are expressed respectively as (3) and (4)

ˆ + 1T C−1 (y − 1μ ˆ) yˆ (x) = μ

82

where φ is the Gaussian probability density function.

k =1

51 52

81

Fig. 1. Flowchart of the KMOP.

112 113 114 115 116

3. Expected angle-penalized length improvement criterion and its validation

117 118 119

In the present study, a novel EAPLI multipoint infill criterion is developed and incorporated with the KMOP as introduced in the last section. Fig. 2 presents the flowchart of the KMOP incorporated with EAPLI. It mainly includes 8 successive steps and starts with the generation and evaluation of initial samples, which is followed by generation and clustering of reference vectors. Then Kriging models for each objective are built and ideal and nadir points are obtained. Non-dominated samples are associated with reference vectors and the reference APL for each vector is calculated. Maximization of EAPLI for each vector is done in parallel and a large number of candidate points are obtained. The desired number of infill points are selected from the candidates in the last step of the iterative process.

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

JID:AESCTE

AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.4 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

4

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

Fig. 3. Reference vectors and sampling points in 2-objective space. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

15 16 18

Fig. 2. Flowchart of KMOP incorporated with EAPLI.

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

It is noted that maximizing EAPLI usually acts as a single-point infill criterion, but maximizing EAPLI given in the flowchart essentially acts as a multipoint infill criterion. To achieve this, efforts are made to extend EAPLI from a single-point criterion to multipoint infill one. Territories and niche counts are introduced to assess the fitness of each candidate which helps to select the desired number of infill points from the candidates. In the following subsections, for the purpose of easy understanding, the EAPLI single-point infill criterion is briefed firstly then the development of the multipoint infill EAPLI criterion is detailed and followed by its validations. 3.1. EAPLI as a single-point infill criterion In adaptive Kriging-based multiobjective optimization, the quality of infill sampling points in objective space is characterized by both high levels of convergence and diversity. How to measure the quality of infill sampling points is crucial for infill criterion, since it largely influences the convergence and diversity of the optimization solutions. In the present study, a reference vector-based APL approach is proposed, and it is well-suited for justifying the quality of points in the objective space. On such a basis, EAPLI is developed. In the following subsections, APL is depicted first and then followed by EAPLI.

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

81 82 83

17 19

80

3.1.1. Reference vector-based APL For easy understanding, 2-objective optimization is used as an example and Fig. 3 depicts a 2-objective (of f 1 and f 2 ) space. Say 2 different sampling points of Point 1 and Point 2 in the objective space, and their positions are described by the position vectors denoted by the blue solid arrow lines in Fig. 3. The position vector starts from the same point, termed as the starting point, and for a detailed description of such a starting point, one can refer to section 3.2. The convergence of the sampling point can be directly described by the length of position vector (i.e. the distance of the sampling point away from the starting point): the shorter it is, the better the convergence would be; vise verse, being that the optimization goal is to search for the minimum of objectives. Obviously, the convergence of Point 1 is better than Point 2 in Fig. 3. Diversity of the sampling points can’t be readily described by a single parameter and it requests a suited scheme. As it is known, a good diversity of the sampling points means they are welldispersed over the objective space, according to which, a reference vector-based scheme is proposed to measure the sampling point diversity.

In the proposed scheme, reference vectors are used and they play significant roles. They are uniformly distributed vectors and start also from the same starting point as the position vectors do, being denoted by the dashed arrow lines in Fig. 3. For a detailed description of reference vector generation, one can refer to section 3.2. As shown in Fig. 3 by the solid squares, if the sampling points are uniformly distributed over or nearby each of the reference vectors, a good diversity of the sampling points will then be achieved. Clearly, the diversity of sampling point can be described by the angle, θ , between the position vector and its closest reference vector: the smaller the angle is, the closer the sampling point would be to the reference vector, indicating better diversity of the sampling point; vise verse. As for sampling Point 1 and Point 2 in Fig. 3, their diversity is respectively described by θ1 and θ2 . Obviously, Point 2 is of better diversity, since θ2 is smaller than θ1 . As discussed above, convergence and diversity of sampling points are separately represented by the length of the position vector and by angle θ . But how to allow for both aspects simultaneously remains unanswered. To resolve such a problem, APL is proposed and it helps to reach a good compromise between both convergence and diversity. It is the position vector length penalized with an angle-based factor and written as







APL f(x), v, z∗ , γv = 1 + α with

θ = cos−1



f(x) · v

θ

γv



  · f(x)

(8)

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

(9)

f(x)

where z∗ is the starting point; α is a penalty parameter, and it takes a value of 0.05 in the present study; θ represents the angle between position vector and reference vector v; γv is the constant angular parameter and used to normalize the angle θ ; f(x) is the length of the position vector. A 2-objective illustration of elements in the APL is presented in Fig. 4.

112 113 114 115 116 117 118 119 120

3.1.2. Definition of EAPLI The APL Improvement (APLI) for a reference vector v is defined as the difference between the reference APL and APL of a concerned sampling point for v, written as



APLI F, f, v, z∗ , γv







121 122 123 124 125

∗

= max APLref F, v, z , γv





  − APL f, v, z , γv , 0 ∗

126

(10)

where APLref is the minimum APL among the APLs of n N D S NDSs NDS for the reference vector v, and F = (f(˜x1 ), . . . , f(˜xn )) represents i the objective values of the NDSs (x˜ , i = 1, . . . , n N D S ) associated with the reference vector v.

127 128 129 130 131 132

JID:AESCTE AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.5 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

1 2 3 4 5 6

9 10 11 12 14

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

Fig. 5(a) depicts the APLref in 2-objective space, where a reference vector v is present and denoted by the black dashed arrow line. Say 3 NDSs are present and denoted by the solid green dots, and they all are associated with reference vector v. For the purpose of easy illustration, all of the NDSs are supposed to be projected onto the reference vector v being denoted by the hollow white dots, where their positions are described by the position vectors being denoted by the blue solid arrow lines and the lengths of the position vectors are angle-penalized length for the NDSs, i.e. the so-called APLs. Whereas the smallest among the APLs is used as APLref , and its corresponding contour line is given and denoted by the green dotted-dash line, on which the APL is identical. On the basis of APLref , a brief explanation of (10) is given by means of Fig. 5(b). Suppose a concerned sampling point f(x) in 2-objective space: as it falls within the region enclosed by the contour line like f(x1 ), its APLI will be the difference between the APLref and the APL of f(x1 ); as it falls outside the aforementioned enclosed region like f(x2 ), the APL takes a zero value. Substituting the APLI into equation (7), the EAPLI is obtained and expressed as



∗

EAPLI x, F, v, z , γv



= y∈ Rm





APLI F, f, v, z∗ , γv

m  1 i =1

si

 φ

y i − yˆ i si

 dy i

(11)

41 42 43 44 45 46 47 48

The value of EAPLI is numerically calculated by Monte Carlo integrations with 1000 points. A candidate infill point is obtained by maximizing the EAPLI of one reference vector. This is essentially a single-objective optimization problem, and it is then solved by Particle Swarm Optimization (PSO) algorithm. So far, the EAPLI used in the single-point infill criterion is involved, how to extend EAPLI into multipoint infill criterion is detailed below.

49 50

3.2. Multipoint infill criterion based on EAPLI

51 52 53 54 55

As shown in Fig. 2, incorporation of KMOP with EAPLI is realized through 8 steps, among which, several elements are crucial in establishing the EAPLI based multipoint infill criterion and they are detailed below.

56 57 58 59 60 61 62 63 64 65 66



70

u j = (u 1j , . . . , um ) j

vj =

Fig. 4. Elements in APL.

13

16

67

∈

{ H0 , H1 , . . . , HH },

i i =1 u j

=1

(12)

3.2.1. Generation and clustering of reference vectors As mentioned above, a candidate infill point can be obtained by maximizing EAPLI for one reference vector. To establish the multipoint infill criterion, multiple reference vectors are required. The distribution of reference vectors directly influences the distributive feature of NDSs in the objective space. To obtain well-distributed NDSs, uniformly distributed reference vectors are essential. The approach similar to that detailed in [41] is used, where all the unit reference vectors are placed uniformly with the starting point. The implementation is given below.

69

uj

(13)

u j 

72 73 74 75 76 77 78 79 80

which maps the points from a hyperplane to a hypersphere as depicted by an example in Fig. 6. In the above transformation operation, the number of points and vectors are both identical and takes a value of N. Given H and m, then N is expressed as

N=

68

71

m

where j = 1, . . . , N and N is the number of uniformly distributed points and H is a positive integer for the simplex-lattice design. Subsequently, the corresponding unit reference vectors v j are obtained by a transformation operation

8

15

First, a set of uniformly distributed points on a unit hyperplane are generated using the canonical simplex-lattice design method [47], written as

u ij

7

5

−1 Cm H +m−1

(14)

−1 Cm H +m−1

81 82 83 84 85 86 87

where is the number of alternatives for choosing (m − 1) out of ( H + m − 1) objects. In the 3D objective space as shown in −1 Fig. 6, in total N = C 43+ 3−1 = 15 reference vectors, are generated for H = 4. The reference vectors are predefined and hold unchanged during the iterative phase of KMOP (see Fig. 2). For each vector, the smallest angle between this vector and the rest vectors is used as the γV for the APL calculation by formula (8). If EAPLI is used as a multipoint infill criterion, an arbitrary number (nadd ) of sample points will be added simultaneously during one iteration, where the reference vectors can be classified into nadd clusters by the k-means method according to reference vector endpoint values.

88

3.2.2. About ideal and nadir points The ideal point and nadir point are used to normalize the objective space. The ideal and nadir points are obtained respectively by searching for the minimum and maximum objective values of the NDSs with the Kriging models by means of optimization algorithm where NSGA-II is used in this work. In such a manner, the ideal point and nadir point are updated at every iteration. Fig. 7 illustrates an ideal and nadir point in 2-objective space. The ideal point is used as the starting point for all reference vectors.

101

3.2.3. Association of NDSs and definition of reference APL An NDS x˜ i will be regarded as being associated with the reference vector v j if the angle θi , j (θi , j = cos−1 ((f(˜xi ) · v j )/f(˜xi ))) between the objective vector f(˜xi ) and the reference vector v j is minimal among all the reference vectors. An illustration of one NDS association is presented in Fig. 8, where the NDS x˜ 1 is regarded as being associated with v3 since θ1,3 is minimal. The APLs are evaluated for all NDSs of the associated reference vectors, among which the minimum APL of each reference vector is used as the reference APL value. Say x˜ i is being associated with a reference vector v j and written as x˜ i ∈  j . Then, the reference APL for v j is expressed as

111



APDre f f, v j , z∗ , γv, j



   = min APD f x˜ i , v j , z∗ , γv, j x˜ i ∈

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

(15)

j

For the reference vectors without NDS association, a larger reference APL than the maximal one among the vectors will be assigned, where 1.1 folds is adopted in the present study. As shown in Fig. 9, the reference APL for the vector v4 without NDS being associated is 1.1 APLref5 , as the reference APL of the vector v5 is the largest among the 4 reference APLs.

125 126 127 128 129 130 131 132

JID:AESCTE

AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.6 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

6

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

Fig. 5. Definition of Reference APL and APLI.

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 99

33

Fig. 8. Association of NDS.

34

Fig. 6. Generation of reference vectors.

35

100 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

Fig. 9. Reference APLs for vectors. Fig. 7. Ideal and nadir points in 2-objective space.

52

118 119

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

116 117

51

3.2.4. Maximization of EAPLI in territories for each reference vector To guide the search of infill point surrounding the reference vector, exploration territories are assigned to each reference vector v, written as





T f(x), v, z∗ =

π 4H

120 121 122 123 124

−θ ≥0

(16)

where θ is the angle between the concerned point objective vector f(x) and reference vector v being calculated by the formula (9). The territories are hyper circular cones (i.e. the sectors in 2-objective space shown in Fig. 10. If the estimated objective values of the concerned point ˆf(x) does not satisfy (16) (i.e. the left side takes a negative value), the

125 126 127 128 129 130 131

Fig. 10. Territories in 2-objective space.

132

JID:AESCTE AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.7 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

1 2 3

left-hand side value of (16) will then be assigned to EAPLI, so as to guide the optimal searching into the territory of positive EAPLI, i.e. around the reference vector. It is written as

4 5 6



EAPDI x, F, v, z∗ , γv

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

=



⎧ m 1 y i − yˆ i ∗ ⎪ ⎪ ⎪ y ∈ R m APDI(F, f, v, z , γv ) i =1 si φ( si )dy i ⎪ ⎪ ⎪ ⎪ ⎨ T (fˆ(x), v, z∗ ) ≥ 0 ⎪ ⎪ T (fˆ(x), v, z∗ ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ T (fˆ(x), v, z∗ ) < 0

(17)

This approach also eliminates tuning the penalty parameter as used in (8). The N candidate infill points are obtained by maximizing EAPLI for all the N reference vectors by means of PSO. For a problem of N reference vectors, in total N single-objective optimization problems. Such a problem-solving process is of good parallel feature, it can be solved fast and efficiently with parallel computing technique. 3.2.5. Selection of infill points from candidates The multipoint infill criterion is derived from EAPLI by selecting nadd infill points from N candidates. It is much straightforward to choose candidates corresponding to the maximum EAPLI in each cluster since the reference vectors are already classified into nadd clusters. However, increasing the number of selected candidates (i.e. increasing the value of nadd ) can possibly cause infill points mal-distribution being concentrated in certain regions (where the Kriging models overestimate the function values), which reduces the diversity. The use of niche count can resolve such mal-distribution. For a reference vector v j , a niche count is assigned to each candidate and is given by

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

nc j =

N  i =1

niN D S

θi , j /θ j ,min + 1

(18)

niN D S

where is the number of NDSs for vi , θi , j denotes the angle between vi and v j , and θ j ,min represents the minimum nonzero θi , j . Obviously, a reference vector being close to NDSs has a high niche count. As the above mentioned, the niche count is proposed to prevent the infill points from clustering, and the candidate points relatively being away from the existing NDSs are favorable for maintaining the diversity of the infill points. Consequently, the candidate points of lower niche counts are chosen. The fitness of the jth candidate is defined as the maximized EAPLI of the jth candidate EAPLI j divided by the niche count nc j , which guides the selection of nadd infill points from N candidates, written as

fitness j =

EAPDI j nc j

(19)

With such a strategy, the NDSs of high-level diversity can be obtained by choosing the candidates of the highest fitness from each cluster.

59

3.3. Comparative studies

60 61 62 63 64 65 66

3.3.1. Implementation of datum multipoint EI-based criteria Two different criteria of Pseudo Expected Improvement Matrix criterion (PEIMm ) [40] and Expected Penalty-Based Boundary Intersection Improvement (EPBII) [48] are used as the datum to justify the presently developed EAPLI. Their implementations in the KMOP framework are distinctive. For the maximum distance-based PEIMm , nadd infill points are simultaneously obtained by maximizing PEIMm using PSO at every update iteration; for the PenaltyBased EPBII, N weight vectors are generated and classified into nadd clusters, then nadd solutions of the highest fitness relevant to EPBII in each cluster are added into the sample set at every update iteration. Monte Carlo integration over 1000 points is used in the calculation of EPBII.

In this section, numerical experiments are conducted for 9 distinctive test problems to study the performance of EAPLI using 5 infill points, where the same number of infill points was used by some other previous investigators [38,40,48]. The proposed EAPLI is compared with the two existing multipoint EI-based infill criteria and verified.

67 68 69 70 71 72 73 74 75 76 77 78 79 80

3.3.2. Test problems Test problems are selected from the previous MOEA/D-EGO studies [38], and they include 2-objective LZ08-F1-F4 [49], 2-objective ZDT1-ZDT3 [50], and DTLZ2 [51] with respectively 3 and 4 objectives. KMOP optimization is respectively conducted for each test problem with the presently developed multipoint infill criterion and two datum multipoint EI-based criteria.

81 82 83 84 85 86 87 88

3.3.3. Optimization settings The number of decision variables n is set as 6. The number of initial sample points generated by LHS is set as 65 (= 11n − 1) as suggested in [15]. The number of infill points in each update iteration nadd is identically 5 for PEIMm , EPBII, and EAPLI. The KMOP will stop as the total number of expensive function evaluation for initial sample points and added infill points reach 200 for 2-objective and 300 for the 3-objective and 4-objective problems. The number of reference vectors is set as N = 101( H = 100) for the 2-objective problems, N = 231( H = 20) for the 3-objective problem and N = 286( H = 10) for the 4-objective problem based on the preliminary test results. The number of independent runs is 10. The DACE toolbox [52] is used for constructing the Kriging models in the experiments with the zero-order polynomial regression function and Gaussian correlation function. Two optimizers are used: NSGA-II for exploring NDSs of Kriging models to determine the ideal and nadir points and PSO for maximizing the infill criterion. Limiting populations and generations of NSGA-II are set respectively as 500 and 100 while the crossover probability and parameter are set respectively as 0.9 and 20; the mutation probability mutation parameters are set as 1/n and 20 respectively. Swarm size and maximum generations of PSO are 150 and 100 respectively.

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

3.3.4. Performance metrics To justify the exploration capacity of KMOPs with the datum and developed criteria, the Inverted Generational Distance (IGD) [53] of the final NDSs obtained by each criterion is compared. Such a single metric describes both convergence and diversity features of the obtained NDSs [54]. Let P ∗ be a set of uniformly distributed points in the objective space over the PF and P be the obtained NDSs approximating the PF, then IGD is given as



IGD =

58

7

∗ p ∈ P d( p , P )

| P ∗|

114 115 116 117 118 119 120 121 122 123

(20)

where d( p , P ∗ ) denotes the minimum Euclidean distance from a point p to any arbitrary point of P ∗ , | P ∗ | is the number of PF solutions. The PFs of the test problems are given, and they are represented by uniformly distributed grid points over the objective space. The PF is represented respectively by 1000 points for the 2-objective problems, 5041(= 712 ) for the 3-objective problem and 9261(= 213 ) for the 4-objective problem.

124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.8 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

8

1 2 3

Problem

m

4 5 6 7 8 9 10 11 12

67

Table 1 Statistical results of IGD values.

LZ08-F1 LZ08-F2 LZ08-F3 LZ08-F4 ZDT1 ZDT2 ZDT3 DTLZ2

2 2 2 2 2 2 2 3 4

68

EAPLI

PEIMm

min

mean

max

0.0066 0.0593 0.0330 0.0173 0.0081 0.0242 0.1059 0.0454 0.1175

0.0079 0.1090 0.0422 0.0309 0.0101 0.0533 0.3589 0.0522 0.1249

0.0109 0.1975 0.0576 0.0652 0.0159 0.1252 0.5698 0.0661 0.1446

↑ ↓ ≈ ≈ ↑ ↑ ↓ ↑ ↑

EPBII

min

mean

max

0.0102 0.0376 0.0241 0.0206 0.0122 0.0482 0.0248 0.0775 0.1670

0.0124 0.0589 0.0363 0.0251 0.0159 0.0803 0.0561 0.0858 0.1708

0.0153 0.0768 0.0645 0.0309 0.0218 0.1315 0.1127 0.0951 0.1761

↑ ≈ ≈ ↑ ↑ ↑ ≈ ↑ ↑

MOEA/D-EGO

69 70

min

mean

max

mean

0.0130 0.0520 0.0289 0.0288 0.0673 0.0312 0.3085 0.0934 0.1666

0.0317 0.1030 0.0513 0.0454 0.1929 0.1807 0.4395 0.1059 0.2017

0.0755 0.2393 0.1023 0.0657 0.3007 0.6095 0.6360 0.1185 0.2659

0.0100 0.1704 0.1488 0.0844 0.0148 0.0156 0.0665 0.1398 –

15

72 73 74 75 76 77 78 79

13 14

71

80

Table 2 IGD values of NDSs obtained with different number of infill points for LZ08F1 problem.

81 82

16 17

nadd

3

5

10

15

20

30

40

50

70

100

83

18

niter

45

27

14

9

7

5

4

3

2

2

84

19

IGD

0.0056 0.0068 0.0107

0.0066 0.0079 0.0109

0.0062 0.0090 0.0141

0.0079 0.0091 0.0116

0.0079 0.0098 0.0180

0.0066 0.0093 0.0121

0.0076 0.0117 0.0200

0.0082 0.0117 0.0238

0.0096 0.0126 0.0160

0.0064 0.0126 0.0218

85

20 21

min mean max

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

87 88

22 23

86

3.3.5. Comparisons and discussions In the KMOP, random operations are largely involved in different procedures like initial sample generation and EA-based maximizing EI. To eliminate the influence of such random operations, 10 independent runs of optimization are conducted for each problem. Table 1 presents the minimum, mean and maximum of IGD obtained through the KMOP with the developed infill criterion and two datum ones, where the minimum and maximum is chosen from the 10 run results and mean is the averaged IGD. For each performance criterion and test problem, the best results of the three criteria are highlighted in boldface with underline and the second-best ones are highlighted in boldface. A Wilcoxon ranksum test is used to compare the results obtained by KMOP with EAPLI and the other two criteria at a significance level of 0.05. In Tables 1–2, the symbol ↑ indicates that the KMOP with EAPLI performs statistically better than with the datum of PEIMm and EPBII;↓ indicates inversely that the datum performs better than EAPLI, while ≈ indicates that insignificant difference is produced by the KMOP with EAPLI and the datum. Besides, the mean values of IGD obtained by MOEA/D-EGO [38] are also used as an additional comparison datum, as shown in the last column of Table 1. The NDSs nearest to the mean IGD value obtained with each infill criterion are presented in Figs. 11–13 respectively for each problem. For test problems of LZ08-F1, ZDT1, ZDT2, and 3-objective and 4-objective DTLZ2, EAPLI outperforms PEIMm and EPBII in terms of IGD metrics as shown in Table 1. The plots of NDSs in Figs. 11(a), 11(e), 11(f), 12 and 13 also show that the NDSs obtained with EAPLI are better than those with PEIMm and EPBII. In particular, for the 3-objective and 4-objective DTLZ2 problems, the NDSs obtained with EPADI are more uniformly distributed than those with PEIMm and EPBII. For the test problem of LZ08-F2, PEIMm performs best in terms of IGD values as shown in Table 1, but the NDSs obtained with the three criteria are not well converged, Fig. 11(b), which indicates this problem may be difficult for the Kriging-based method. For test problems of LZ08-F3 and LZ08-F4, EAPLI produces a comparable performance to PEIMm and EPBII in terms of IGD values shown in Table 1 and NDSs of Figs. 11(c) and 11(d). For the test problem of ZDT3 with a discrete PF, with all the three criteria, the NDSs over the sectional PFs are missed out more or less, Fig. 11(g). With PEIMm and EAPLI, 1 sectional PF and 3 sectional PFs are missed out respectively. Moreover, the NDSs obtained with EAPLI are still close to the PF. It also can be seen from Table 1 that KMOP with EAPLI outperforms MOEA/D-EGO in terms of IGD met-

ric for all the test problems of LZ08-F1-F4, ZDT1, and 3-objective DTLZ2. With the above observations, it can be concluded that the performance of EAPLI using a common number of infill points, i.e. 5, is well-demonstrated.

89 90 91 92 93

3.4. High capacity of EAPLI for multiple-point infill

94 95

It is acknowledged that the multipoint infill criteria help to enhance the optimization efficiency and accelerate the optimization process with parallel computations. That is, the multiple infill points can be evaluated simultaneously by high-fidelity simulation tools with parallel computation. Thus the time for evaluation of multiple points or one point is identical. With this, the duration of the optimization would depend on the number of model updates, hence fewer model updates enhance the optimization efficiency. With enough computational resources available, the larger the multipoint infill number is, the higher the optimization efficiency would be. But to date, the multipoint infill number reported in the open literature is up to 10. The presently developed EAPLI criterion is of distinctive feature in this respect as it supports a considerably larger number of multipoint infill. For instance, it works efficiently with 100 infill points in solving the LZ08F1 problem as detailed below, which demonstrates the great potentials in accelerating the optimization of expensive multiobjective problems. To test the EAPLI capacity for multipoint infill, the LZ08F1 problem is used, where to let nadd = 3, 5, 10, 15, 20, 30, 40, 50, 70, 100 respectively and all other settings are identical to that used by EAPLI with nadd = 5 in Section 3.3. The minimum, mean and maximum IGD values obtained over 10 independent runs are presented in Table 2, where the values for nadd = 5 are taken from Table 1 for comparison convenience. The NDSs of the run with IGD value closest to the mean IGD value obtained with the different number of infill points are presented in Fig. 14. As shown in Table 2, the mean IGD value increases with the number of infill points. The mean IGD value of NDSs obtained with EAPLI for nadd = 100 is 0.0126 being the worst of all the cases listed in Table 2, but it is yet comparable to that of the PEIMm for nadd = 5 that has a mean IGD value of 0.0124 as shown in Table 1. It can be concluded that EAPLI has permitted the well-distributed NDSs for a larger number of infill points but fewer model updates. Such distinctive features can help to significantly reduce the iterations, hence accelerate the optimization process of adaptive Kriging-based optimization of expensive multiobjective problems.

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:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.9 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

9

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

121

56

122 123

57 58

Fig. 11. NDSs with IGD nearest to the mean value obtained for 2-objective problems.

60

4. Preference-driven capacity of EAPLI

61 62

4.1. Preference-driven KMOP with EAPLI

63 64 65 66

124 125

59

Most expensive multiobjective optimizations of practical applications have owned the designer/decision maker’s preference. However, the preference-driven feature has rarely been realized for

the adaptive sampling based KMOP method. With the presently developed EAPLI, the designer/decision maker’s preference can be readily incorporated to facilitate preference-driven optimization as the distribution of vectors used can be manipulated easily. The preference-driven KMOP method with EAPLI shares the same framework of general optimization in Fig. 2. But implementations of two procedures are different which are introduced below.

126 127 128 129 130 131 132

JID:AESCTE

AID:105555 /FLA

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

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Fig. 11. (continued)

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

Firstly, the reference vector generation method is different as vectors are generated in the designer/decision-maker preferred subspaces only under this situation. The adopted preference vector generation method [55] is introduced as follows. A central vector vc and radius r are determined to set the preferred space including position and size, where vc is a unit vector and r ∈ (0, 1). With this, the reference vectors inside the subspace are generated using the following transformation

v i =

vc + r · (vi − vc )

(21)

where i = 1, . . . , N is the index of each reference vector, vi represents a uniformly distributed vector generated with (12) and (13), and v i denotes the reference vector transformed from vi lying inside the subspace specified by vc and r. Fig. 15 illustrates the above transformation procedure. 10 uniformly distributed reference vectors are generated √ inside a region of 2-objective space specified √ by vc = ( 2/2, 2/2), and r = 0.3. Notably, the above scheme for preference-based vector generation is of flexibility where an arbitrary number of preferred subspaces can be readily obtained through specifying the corresponding pairs of central vectors and corresponding radii. Next, the size of exploration territory needs to be modified since the preference-based reference vectors are not welldistributed from a point of view of the entire objective space. EAPLI exploration territory of the reference vector v i is defined as





θi ,min 2

− θ ≥ 0.

(22)

Where θi ,min represents the minimum nonzero angle between v i and other vectors, θ is the angle between the concerned objective vector f(x) and reference vector v i being calculated by the formula (9).

61 62 63 64 65 66

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

99 100

5.1. Implementation

101

KMOP with EAPLI is used as an optimizer and integrated with an evaluation module, which consists of FFD parameterization [56] and XFOIL flow solver [57], to form a multiobjective airfoil optimization method. Fig. 17 illustrates such integration. When the initial samples and infill points need to be evaluated by the real functions, the evaluation module is involved with the decision variables transferred. In the evaluation module, the airfoil is generated firstly by means of FFD; then XFOIL flow solver is evoked automatically to predict the flow; finally, the objective functions are calculated on the basis of flow prediction and feedback the objectives to the optimizer. The FFD parameterization, XFOIL flow solver, and objective function definitions are briefed below.

102 103 104 105 106 107 108 109 110 111 112 113 114 115

FFD parameterization A commonly used geometric parameterization approach of FFD is used [58–60]. With the FFD approach, the geometry shape parameterization is realized in two steps. In the first step, the parametric coordinates (s, t ) of the original airfoil are computed. As shown in Fig. 18, for a point of Cartesian coordinates (x, y ) in the deformation rectangular region, the parametric coordinates s and t are calculated by

116 117 118 119 120 121 122 123

4.2. Demonstrative examples of preference-driven optimization

59 60

68

98

5. Airfoil shape optimization

57 58

67

97

vc + r · (vi − vc )

T f(x), v i , z∗ =

number of central vectors in each case. The KMOP with EAPLI stops as 60 infill points are added. Other settings for constructing the Kriging models and optimizers hold identical to that introduced in Section 3. Fig. 16(a) and (b) present the obtained NDSs with different preferences for LZ08-F1 and DTLZ2 respectively. The optimal searching is well-driven by the preferences in both problems where most NDSs appear in the preferred corner and center regions in accord with the reference vectors. But there is a small number of NDSs unexpectedly appear outside the preferred space. The reasons may be the following: in the first several update iterations, the starting point of reference vectors are not wellestimated since the initial Kriging models constructed roughly with DOE are lack of accuracy, which makes the reference vectors being biased away from the real starting point of coordinate origin of (0, 0, 0) for these two problems. Clearly, infill points in the first several update iterations are located in the corner and center of the objective space due to poor-estimated reference vector starting points. As the iteration proceeds, the estimation of starting points are getting better and better, and they eventually approach to the real ones, thus the infill points/NDSs are found in the preferred subspaces merely. It is noted that the corner vectors are also involved in dealing with the center preferences of objective space, i.e., {(1, 0), (0, 1)} in 2-objective LZ08-F1 and {(1, 0, 0), (0, 1, 0), (0, 0, 1)} in 3-objective DTLZ2 since the determination of ideal and nadir point requires the extreme values of each objective (see Section 3.2.2). The use of corner sample points increases the accuracy of Kriging models in the corresponding region, and the estimation accuracy of ideal and nadir points can then be enhanced.

The ability of preference-driven KMOP with EAPLI is demonstrated from solving 2-objective LZ08-F1 and 3-objective DTLZ2, which is depicted below. Table 3 presents the central vectors, radii, and the number of reference vectors used to specify the preferred subspaces of the example problems. The number of decision variables and initial samples is set as 6 and 65 respectively. The number of infill points in each update iteration nadd is identical to the

s=

x − X min X max − X min

,t =

y − Y min Y max − Y min

(23)

where ( X min , Y min ) and ( X max , Y max ) are the corner points of the variation space. In the second step, the control points scattered over the variation space are determined and then the airfoil deformation is accomplished through moving these control points. The deformed airfoil coordinates are expressed by

124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105555 /FLA

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

Y. He et al. / 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 81

15 16

82

Fig. 12. NDSs with IGD nearest to the mean value obtained of 3-objective DTLZ2 problem.

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 97

31 32

Fig. 13. Parallel coordinate plot of the NDSs with IGD nearest to the mean value obtained of 4-objective DTLZ2 problem.

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 120

54 55

Fig. 14. PF of LZ08F1 problem and NDSs (filled circles) with IGD nearest to mean value obtained with a different number of infill points.

121

56

122

57

123

58 59

124

Table 3 Central vectors and radii of preference incorporation.

125 126

60 61

Problem

m

Central vectors vc

Radius r

Vector number for each subspace

nadd

LZ08-F1

2

{(1, 0), (0, 1)} √ √ {(1, 0), (0, 1), ( 2/2, 2/2)}

{0.1, 0.1} {0.05, 0.05, 0.1}

{10, 10} {3, 3, 10}

2 3

129

DTLZ2

3

{(1, 0, 0), (0, 1, 0), (0, 0, 1)} √ √ √ {(1, 0, 0), (0, 1, 0), (0, 0, 1), ( 3/3, 3/3, 3/3)}

{0.1, 0.1, 0.1} {0.1, 0.1, 0.1, 0.2}

{10, 10, 10} {3, 3, 3, 10}

3 4

131

62 63 64 65 66

127 128 130 132

JID:AESCTE

AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.12 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

12

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 81

15

Fig. 17. Flowchart of the multiobjective airfoil optimization system.

16

84

18 19

82 83

17

85

Fig. 15. Generation of reference vectors inside the preferred subspace.

20

86

21

87

22

88

23

89

24

90

25

91

26

92 93

27

Fig. 18. FFD parametric coordinates.

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 111

45

Fig. 19. Parameterization with FFD.

46 48 49 50 51

Fig. 16. NDSs obtained with reference vectors incorporating preferences for different problems, (a): 2-objective LZ08-F1 problem, (b): 3-objective DTLZ2 problem.

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

112 113

47

x( s , t ) =

n  m 

n (Pi j + Pi j ) B m i (s) B j (t )

Fig. 19 presents the deformed airfoil and baseline NACA0012 airfoil together with their control points. In addition, during the deformation, the abscissa of each control point is held while its ordinate is varied within a range [−0.05, 0.05].

114 115 116 117 118

(24)

j =0 i =0

where Pi j are the parametric coordinates of the control points,

Pi j the coordinate variation of the control points; B m (t ) and i B nj (t ) are both Bernstein polynomials, m and n are the degrees of the FFD function, indicating that m + 1 vertical columns and n + 1 horizontal rows of the control points, i.e. in total (m + 1) • (n + 1) control points are used. NACA0012 airfoil is parametrized with the FFD. As shown in Fig. 19, the diagonal corner coordinates of the rectangular deformation region are set to (0, −0.25) and (1, 0.25) while m and n are set to 2 and 2 respectively. Thus, 9 uniformly distributed control points are used in the NACA0012 airfoil parameterization.

XFOIL flow solver The open-source XFOIL code [57] is used to simulate the flow around each airfoil and calculate its lift and drag coefficient throughout the optimization. The code combines the panel method and integral boundary layer formulation for simulating the airfoil flow. To eliminate the influence of spatial discretization, a sensitivity study is performed. The panel number of 160 seems sufficient for the airfoil simulation, thus it is used throughout the present study. The incoming flow condition is given, where the angle of attack is set to 0 degrees, free stream Mach 0.4 and Reynolds number 2500000.

119 120 121 122 123 124 125 126 127 128 129

Definition of objective functions For airfoil, least drag and most lifts are always expected and its design is essentially a bi-objective optimization problem. The airfoil shape optimization aims to si-

130 131 132

JID:AESCTE AID:105555 /FLA

[m5G; v1.261; Prn:14/11/2019; 15:28] P.13 (1-15)

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

1 2 3

67

Table 4 Central vectors and radii of preference.

68

Preferred region

Central vectors vc

Radius r

Vector number for each subspace

69

Upper corner Center Lower corner

(√ 0, 1) √ ( 2/2, 2/2) (0, 1)

0.15 0.15 0.15

6 6 6

71

4 5 6 7

13

70 72 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 26 27 28 29

multaneously maximize the lift coefficient and minimize the drag coefficient, given by



30 31 32 33 34 35 36 37

91

Fig. 20. Reference vectors used for conducting preference-driven optimization.

minimize

−C L /C D

92

Fig. 21. NDSs obtained by non-preference and preference-driven optimizations.

94 95



CD

96

(25)

97 98

where C L is the lift coefficient and C D is the drag coefficient. The objective values for the baseline NACA0012 airfoil are −C L /C D = 0 and C D = 0.0054 respectively.

99

5.2. Results and discussions

103

100 101 102 104

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

93

Non-preference and preference-driven optimizations of NACA0012 airfoil by KMOP with EAPLI, and non-preference optimization of NACA0012 airfoil by PEIMm are conducted for a fair comparison, where the 45 identical initial samples are used for all the cases. In the non-preference optimization of EAPLI, 100 uniform distributed vectors are used. In the preference-driven optimization of EAPLI, three regions of the upper corner, center and lower corner are simultaneously set to the preferences. They constitute the preferred subspaces and are specified by the central vectors, radii, and reference vectors, as given in Table 4 and Fig. 20. For all three optimization runs by EAPLI or PEIMm , the number of infill points of each update iteration and the maximal iteration number is set to 3 and 40 respectively. Fig. 21 compares the obtained NDSs for both non-preference and preference-driven optimizations by EAPLI, and NDSs by PEIMm . The followings can be observed: (1) obviously, the NDSs arising from all optimizations are much better than the baseline airfoil in terms of both the lift/drag coefficient ratio and drag coefficient (for baseline airfoil, −C L /C D = 0 and C D = 0.0054); (2) regarding the two non-preference optimizations respectively by EAPLI and PEIMm , the solutions by EAPLI are better than that of PEIMm in the upper corner region but they are worse than PEIMm in the central region; (3) in contrast to the non-preference optimizations, the NDSs of preference-driven optimization by EAPLI are well driven into the preferred subspaces; (4) the NDSs of the center preference-driven optimization by EAPLI are better than that of non-preference optimization by EAPLI; (5) the NDSs of the center and corner preference-driven optimizations by EAPLI are

105 106 107 108 109

Fig. 22. Airfoil shape and surface pressure distributions of chosen NDSs and baseline airfoil.

110 111 112 113

better than that of non-preference optimization by PEIMm . For preference-driven optimization, the infill points are selected in the preferred subspaces rather than in the whole objective space, so that with the identical number of samples, the Kriging model produces higher-level accuracy prediction in such preferred subspaces than that in the whole objective region with non-preference optimization. Thus the solutions by the preference-driven optimization would be better than that of non-preference runs. Fig. 22 compares the optimized airfoil shapes and their surface pressure distributions of two extreme and one compromise NDSs with the baseline acting as the datum, of which the objective values are presented in Table 5. Obviously, in comparison with the baseline shape, the airfoils of the NDSs are thinner and asymmetric, Fig. 22(a). Such geometric derivations in optimized airfoils produce significant aerodynamic performance improvement. The baseline airfoil surface pressure distribution is overlaid, Fig. 22(b), and produces zero-lift, Table 5; with all the 3 optimized airfoils, the surface pressure distribution is improved significantly, where the Max CL /CD is on the top and followed by the Compromise, and the Min CD is moderate.

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

JID:AESCTE

AID:105555 /FLA

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

Y. He et al. / Aerospace Science and Technology ••• (••••) ••••••

14

1 2 3 4 5

References

Table 5 Objective values of chosen NDSs and baseline airfoil.

68

Objective value

Baseline

max C L /C D

Compromise

min C D

C L /C D CD

0.0 0.0054

181.7 0.0047

159.1 0.0044

112.4 0.0042

6 7 8 9

6. Conclusions

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

A novel multipoint infill criterion for Kriging-based multiobjective optimization is proposed on the basis of EAPLI, and it is well-suited for efficiently solving multiobjective optimization problems with expensive objective functions. 1) The EAPLI based optimization method is used to solve various benchmark test problems, where territory is assigned to each reference vector permitting the EAPLI to explore the objective space much efficiently so that the diversely converged NDSs are obtained approximating the complete PF. The proposed EAPLI exhibits a high capacity for balancing the diversity and convergence of NDSs. 2) The EAPLI is proven of high capacity for a large number of multipoint infill. For instance, in solving the 2-objective LZ08F1 problem, the infill point number increases from the previous several up to 100 where these infill points are simultaneously done in a single update without any significant compromise on optimization performance. Such a competitive feature of EAPLI is of great significance as it significantly reduces the iterations and hence accelerates the optimization process for adaptive Kriging-based optimization of expensive multiobjective problems. 3) The EAPLI is of sound manipulative feature in multiple reference vector distribution, and it is favorable for the designer/decision-maker to include his/her preference, so as to perform preference-driven optimization for expensive multiobjective problems. 2-objective LZ08-F1 and 3-objective DTLZ2 with various preferences are solved respectively and the capacity of EAPLI for incorporating the designer/decision-maker preference is welldemonstrated. Such a feature of EAPLI can largely reduce the wasted searching for NDSs violating the preference and enhance the optimization efficiency. 4) Through the incorporation of EAPLI, FFD parameterization, and flow solver XFOIL, the multiobjective airfoil optimization method is established and applied to optimize the NACA0012 airfoil. Significant improvements in dual objectives of lift/drag coefficient ratio and drag coefficient are achieved from the nonpreference and preference-driven optimization. For the latter, the NDSs are well driven into the preferred subspaces while these violating the preferences are effectively avoided, which demonstrates the high-capacity of the proposed methods for solving the practical preference-driven aerodynamic design problems. An extension of the proposed method to the constrained practical problems is among the authors’ major future concerns.

54 55

Declaration of competing interest

56 57 58 59

The authors declare that they have no conflict of interest in this work.

60 61

Acknowledgements

62 63 64 65 66

67

The work is supported by the National Natural Science Foundation of China under Grant No. 51606141, and China Postdoctoral Science Foundation under Grant No. 2016M602817, both are highly acknowledged by the authors.

[1] M. Kanazaki, K. Tanaka, S. Jeong, K. Yamamoto, Multi-objective aerodynamic exploration of elements’ setting for high-lift airfoil using Kriging model, J. Aircr. 44 (2007) 858–864, https://doi.org/10.2514/1.25422. [2] R. Hilbert, G. Janiga, R. Baron, D. Thévenin, Multi-objective shape optimization of a heat exchanger using parallel genetic algorithms, Int. J. Heat Mass Transf. 49 (2006) 2567–2577, https://doi.org/10.1016/j.ijheatmasstransfer.2005.12.015. [3] E. Benini, Three-dimensional multi-objective design optimization of a transonic compressor rotor, J. Propuls. Power 20 (2004) 559–565, https://doi.org/10.2514/ 1.2703. [4] A. Viviani, L. Iuspa, A. Aprovitola, Multi-objective optimization for re-entry spacecraft conceptual design using a free-form shape generator, Aerosp. Sci. Technol. 71 (2017) 312–324, https://doi.org/10.1016/j.ast.2017.09.030. [5] L. Leifsson, S. Koziel, Y.A. Tesfahunegn, Multiobjective aerodynamic optimization by variable-fidelity models and response surface surrogates, AIAA J. 54 (2016) 531–541, https://doi.org/10.2514/1.J054128. [6] S. Guo, F. Duan, H. Tang, S.C. Lim, M.S. Yip, Multi-objective optimization for centrifugal compressor of mini turbojet engine, Aerosp. Sci. Technol. 39 (2014) 414–425, https://doi.org/10.1061/(ASCE)AS.1943-5525.0000338. [7] E. Iuliano, Global optimization of benchmark aerodynamic cases using physicsbased surrogate models, Aerosp. Sci. Technol. 67 (2017) 273–286, https:// doi.org/10.1016/j.ast.2017.04.013. [8] N. Bartoli, T. Lefebvre, S. Dubreuil, R. Olivanti, R. Priem, N. Bons, J.R.R.A. Martins, J. Morlier, Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design, Aerosp. Sci. Technol. 90 (2019) 85–102, https://doi.org/10.1016/j.ast.2019.03.041. [9] S. Cheng, H. Zhan, Z. Shu, H. Fan, B. Wang, Effective optimization on Bump inlet using meta-model multi-objective particle swarm assisted by expected hyper-volume improvement, Aerosp. Sci. Technol. 87 (2019) 431–447, https:// doi.org/10.1016/j.ast.2019.02.039. [10] J. Yu, Z. Wang, F. Chen, J. Yu, C. Wang, Kriging surrogate model applied in the mechanism study of tip leakage flow control in turbine cascade by multiple DBD plasma actuators, Aerosp. Sci. Technol. 85 (2019) 216–228, https://doi.org/ 10.1016/j.ast.2018.11.055. [11] A.I.J. Forrester, A.J. Keane, Recent advances in surrogate-based optimization, Prog. Aerosp. Sci. 45 (2009) 50–79, https://doi.org/10.1016/j.paerosci.2008.11. 001. [12] G. Matheron, Principles of geostatistics, Econ. Geol. 58 (1963) 1246–1266, https://doi.org/10.2113/gsecongeo.58.8.1246. [13] J. Mockus, V. Tiesis, A. Zilinskas, The application of Bayesian methods for seeking the extremum, in: Towar. Glob. Optim., vol. 2, 1978, pp. 117–129. [14] J. Mockus, Bayesian Approach to Global Optimization: Theory and Applications, Springer, Dordrecht, 1989. [15] M. Schonlau, Computer Experiments and Global Optimization, Ph.D. diss, Dept. Stat. Actuarial Sci., Univ. at Waterloo, Waterloo, 1997. [16] D.R. Jones, M. Schonlau, W.J. Welch, Efficient global optimization of expensive black-box functions, J. Glob. Optim. 13 (1998) 455–492, https://doi.org/10.1023/ a:1008306431147. [17] W. Ponweiser, T. Wagner, M. Vincze, Clustered multiple generalized expected improvement: a novel infill sampling criterion for surrogate models, in: Evol. Comput. 2008, CEC 2008 (IEEE World Congr. Comput. Intell. IEEE Congr.), IEEE, 2008, pp. 3515–3522. [18] J.P.C. Kleijnen, W. Van Beers, I. Van Nieuwenhuyse, Expected improvement in efficient global optimization through bootstrapped Kriging, J. Glob. Optim. 54 (2012) 59–73, https://doi.org/10.1007/s10898-011-9741-y. [19] J.M. Parr, A.J. Keane, A.I.J. Forrester, C.M.E. Holden, Infill sampling criteria for surrogate-based optimization with constraint handling, Eng. Optim. 44 (2012) 1147–1166, https://doi.org/10.1080/0305215X.2011.637556. [20] S. Koullias, D.N. Mavris, Methodology for global optimization of computationally expensive design problems, J. Mech. Des. 136 (2014) 81007, https:// doi.org/10.1115/1.4027493. [21] Y. Zhang, Z.-H. Han, K.-S. Zhang, Variable-fidelity expected improvement method for efficient global optimization of expensive functions, Struct. Multidiscip. Optim. 58 (2018) 1431–1451, https://doi.org/10.1007/s00158-0181971-x. [22] J. Sacks, W.J. Welch, T.J. Mitchell, H.P. Wynn, Design and analysis of computer experiments, Stat. Sci. 4 (4) (1989) 409–423, https://doi.org/10.1214/ss/ 1177012413. [23] D. Huang, T.T. Allen, W.I. Notz, N. Zeng, Global optimization of stochastic blackbox systems via sequential Kriging meta-models, J. Glob. Optim. 34 (3) (2006) 441–466, https://doi.org/10.1007/s10898-005-2454-3. [24] V. Picheny, T. Wagner, D. Ginsbourger, A benchmark of Kriging-based infill criteria for noisy optimization, Struct. Multidiscip. Optim. 48 (3) (2013) 607–626, https://doi.org/10.1007/s00158-013-0919-4. [25] A. Sóbester, S.J. Leary, A.J. Keane, A parallel updating scheme for approximating and optimizing high fidelity computer simulations, Struct. Multidiscip. Optim. 27 (2004), https://doi.org/10.1007/s00158-004-0397-9. [26] D. Zhan, J. Qian, Y. Cheng, Balancing global and local search in parallel efficient global optimization algorithms, J. Glob. Optim. 67 (2017) 873–892, https://doi. org/10.1007/s10898-016-0449-x.

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:105555 /FLA

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

Y. He et al. / 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

[27] D. Zhan, J. Qian, Y. Cheng, Pseudo expected improvement criterion for parallel EGO algorithm, J. Glob. Optim. 68 (2017) 641–662, https://doi.org/10.1007/ s10898-016-0484-7. [28] Z. Feng, Q. Zhang, Q. Zhang, Q. Tang, T. Yang, Y. Ma, A multiobjective optimization based framework to balance the global exploration and local exploitation in expensive optimization, J. Glob. Optim. 61 (2015) 677–694, https:// doi.org/10.1007/s10898-014-0210-2. [29] J. Liu, W.P. Song, Z.H. Han, Y. Zhang, Efficient aerodynamic shape optimization of transonic wings using a parallel infilling strategy and surrogate models, Struct. Multidiscip. Optim. 55 (2017) 925–943, https://doi.org/10.1007/s00158016-1546-7. [30] D. Ginsbourger, R. Le Riche, L. Carraro, Kriging is well-suited to parallelize optimization, in: Comput. Intell. Expens. Optim. Probl., Springer, Berlin, 2010, pp. 131–162. [31] M. Emmerich, Single- and Multi-Objective Evolutionary Design Optimization Assisted by Gaussian Random Field Metamodels, Ph.D. diss, FB Inform., Univ. Dortmund, Dortmund, 2005. [32] Ł. Łaniewski-Wołłk, S. Obayashi, S. Jeong, Development of expected improvement for multi-objective problem, in: Proc. 42nd Fluid Dyn. Conf. Aerosp. Numer. Simulat. Symp., 2010. [33] K. Shimoyama, K. Sato, S. Jeong, S. Obayashi, Comparison of the criteria for updating Kriging response surface models in multi-objective optimization, in: 2012 IEEE Congr. Evol. Comput., 2012. [34] J. Martínez-Frutos, D. Herrero-Pérez, Kriging-based infill sampling criterion for constraint handling in multi-objective optimization, J. Glob. Optim. 64 (2016) 97–115, https://doi.org/10.1007/s10898-015-0370-8. [35] A.J. Keane, Statistical improvement criteria for use in multiobjective design optimization, AIAA J. 44 (4) (2006), https://doi.org/10.2514/1.16875. [36] D.C. Bautista, A Sequential Design for Approximating the Pareto Front Using the Expected Pareto Improvement Function, Ph.D. diss, Dept. Stat., Ohio State Univ, Columbus, 2009. [37] D. Zhan, Y. Cheng, J. Liu, Expected improvement matrix-based infill criteria for expensive multiobjective optimization, IEEE Trans. Evol. Comput. 21 (2017) 956–975, https://doi.org/10.1109/TEVC.2017.2697503. [38] Q. Zhang, W. Liu, E. Tsang, B. Virginas, Expensive multiobjective optimization by MOEA/D with Gaussian process model, IEEE Trans. Evol. Comput. 14 (2010) 456–474, https://doi.org/10.1109/TEVC.2009.2033671. [39] K. Shimoyama, K. Sato, S. Jeong, S. Obayashi, Updating Kriging surrogate models based on the hypervolume indicator in multi-objective optimization, J. Mech. Des. 135 (9) (2013) 094503, https://doi.org/10.1115/1.4024849. [40] D. Zhan, J. Qian, J. Liu, Y. Cheng, Pseudo expected improvement matrix criteria for parallel expensive multi-objective optimization, Adv. Struct. Multidiscip. Optim. (2018), https://doi.org/10.1007/s10898-016-0484-7. [41] J. Molina, L.V. Santana, A.G. Hernández-Díaz, C.A. Coello Coello, R. Caballero, g-dominance: reference point based dominance for multiobjective metaheuristics, Eur. J. Oper. Res. 197 (2) (2009) 685–692, https://doi.org/10.1016/j.ejor. 2008.07.015. [42] L. Ben Said, S. Bechikh, K. Ghedira, The r-dominance: a new dominance relation for interactive evolutionary multicriteria decision making, IEEE Trans. Evol. Comput. 14 (5) (2010) 801–818, https://doi.org/10.1109/TEVC.2010.2041060. [43] R. Wang, R.C. Purshouse, P.J. Fleming, Preference-inspired coevolutionary algorithms for many-objective optimization, IEEE Trans. Evol. Comput. 17 (4) (2013) 474–494, https://doi.org/10.1109/TEVC.2012.2204264.

15

[44] T. Shao, S. Krishnamurty, G.C. Wilmes, Preference-based surrogate modeling in engineering design, AIAA J. 45 (2007) 2688–2701, https://doi.org/10.2514/ 1.27777. [45] M.D. McKay, R.J. Beckman, W.J. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 42 (1) (2000) 55–61, https://doi.org/10.1080/00401706. 2000.10485979. [46] C.K. Williams, C.E. Rasmussen, Gaussian Processes for Machine Learning, MIT Press, Massachusetts, 2006. [47] J.A. Cornell, Experiments with Mixtures: Designs, Models, and the Analysis of Mixture Data, Wiley, New York, 2002. [48] N. Namura, K. Shimoyama, S. Obayashi, Expected improvement of penaltybased boundary intersection for expensive multiobjective optimization, IEEE Trans. Evol. Comput. 21 (6) (2017) 898–913, https://doi.org/10.1109/TEVC.2017. 2693320. [49] H. Li, Q. Zhang, Multiobjective optimization problems with complicated Pareto sets, MOEA/D and NSGA-II, IEEE Trans. Evol. Comput. 13 (2009) 284–302, https://doi.org/10.1109/TEVC.2008.925798. [50] E. Zitzler, K. Deb, L. Thiele, Comparison of multiobjective evolutionary algorithms: empirical results, Evol. Comput. 8 (2) (2000) 173–195, https://doi.org/ 10.1162/106365600568202. [51] K. Deb, L. Thiele, M. Laumanns, E. Zitzler, A. Abraham, L. Jain, R. Goldberg, Scalable test problems for evolutionary multiobjective optimization, in: Evol. Multiobjective, 2005. [52] S.N. Lophaven, H.B. Nielsen, J. Søndergaard, DACE—A MATLAB Kriging Toolbox, Version 2.0 [Online]. Available: http://www2.imm.dtu.dk/projects/dace/, Aug. 2002. [53] C.A. Coello Coello, M. Reyes-Sierra, A study of the parallelization of a coevolutionary multi-objective evolutionary algorithm, in: Proc. Third Mex. Int. Conf. Artif. Intell, 2004. [54] D.A. Van Veldhuizen, G.B. Lamont, Multiobjective Evolutionary Algorithm Research: a History and Analysis, Tech. Rep. TR-98-03, Dep. Electr. Comput. Eng. Grad. Sch. Eng. Air Force Inst. Technol. Wright-Patterson AFB, Ohio, 1998, http://dx.doi.org/10.1.1.35.8924. [55] R. Cheng, Y. Jin, M. Olhofer, B. Sendhoff, A reference vector guided evolutionary algorithm for many-objective optimization, IEEE Trans. Evol. Comput. 20 (2016) 773–791, https://doi.org/10.1109/TEVC.2016.2519378. [56] T.W. Sederberg, S.R. Parry, Free-form deformation of solid geometric models, in: Proceedings of the 13th Annual Conference on Computer Graphics and Interactive Techniques, Dallas: SIGGRAPH, 1986. [57] M. Drela, XFOIL: an analysis and design system for low Reynolds number airfoils, in: T.J. Mueller (Ed.), Low Reynolds Number Aerodynamics, SpringerVerlag, Berlin, Heidelberg, 1989. [58] M. Li, J. Bai, L. Li, X. Meng, Q. Liu, B. Chen, A gradient-based aero-stealth optimization design method for flying wing aircraft, Aerosp. Sci. Technol. 92 (2019) 156–169, https://doi.org/10.1016/j.ast.2019.05.067. [59] J. Li, S. He, J.R.R.A. Martins, Data-driven constraint approach to ensure lowspeed performance in transonic aerodynamic shape optimization, Aerosp. Sci. Technol. 92 (2019) 536–550, https://doi.org/10.1016/j.ast.2019.06.008. [60] J. Tao, G. Sun, J. Si, Z. Wang, A robust design for a winglet based on NURBS-FFD method and PSO algorithm, Aerosp. Sci. Technol. 70 (2017) 568–577, https:// doi.org/10.1016/j.ast.2017.08.040.

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

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