Metaheuristic Approaches to Solve a Complex Aircraft Performance Optimization Problem

The increasing demands for travelling comfort and reduction of carbon dioxide emissions have been considered substantially in the stage of conceptual aircraft design. However, the design of a modern aircraft is a multidisciplinary process, which requires the coordination of information from several specific disciplines, such as structures, aerodynamics, control, etc. To address this problem with adequate accuracy, the multidisciplinary analysis and optimization (MAO) method is usually applied as a systematic and robust approach to solve such complex design issues arising from industries. Since MAO method is tedious and computationally expensive, genetic programming (GP)-based metamodeling techniques incorporating MAO are proposed as an effective approach to minimize the wing stiffness of a large aircraft subject to aerodynamic, aeroelastic and stability constraints in the conceptual design phase. Based on the linear small-disturbance theory, the state-space equation is employed for stability analysis. In the process of multidisciplinary analysis, aeroelastic response simulations are performed using Nastran. To construct metamodels representing the responses of the interests with high accuracy as well as less computational burden, optimal Latin hypercube design of experiments (DoE) is applied to determine the optimized distribution of sampling points. Following that, parametric optimization is carried out on metamodels to obtain the optimal wing geometry shape, elastic axis positions and stiffness distribution, and then the solution is verified by finite element simulations. Finally, the superiority of the GP-based metamodel technique over genetic algorithm is demonstrated by multidisciplinary design optimization of a representative beam-frame wing structure in terms of accuracy and efficiency. The results also show that GP metamodel-based strategy for solving MAO problems can provide valuable insights to tailoring parameters for the effective design of a large aircraft in the conceptual phase.


Introduction
The increasing environmental issues, such as air pollution and global climate change, force the airline industry and aircraft designers to seek more economical airplane designs with the maximum lift-to-drag ratio and minimum structural weight as well as favorable stability. The success of such an aircraft design can provide opportunities to consider the welfare and comfort of passengers in its design phase and lead to a decrease in fuel consumption. To achieve this aim, faster and more efficient approaches to design optimization of lightweight but strong aircraft have been well studied by many researchers [1][2][3].
Finally, sensitivity analysis was conducted to identify the critical input parameters that affect aircraft performance in the conceptual design phase.

Brief Review of Aircraft Multidisciplinary Analysis
In this research, the wing jig shape was optimized using the proposed MAO method. Its aerodynamic shape, elastic axis positions, and stiffnesses were considered in the process of parameterized modeling. Then, the structural deformation was obtained by the static aeroelastic analysis and the deflections of structural nodes were transformed into the deformations of aerodynamic grids using the infinite plate spline (IPS) interpolation. Once the aerodynamic model of the wing cruise shape was constructed, the flow field was estimated and the stability was analyzed using small-disturbance theory. A brief review of aerodynamics, flight stability and aeroelasticity in the multidisciplinary analysis is introduced in this section.

Aerodynamics
Aerodynamics analysis of a typical aircraft is usually performed using high-fidelity CFD tools based on the Reynolds-averaged Naiver-Stokes (RANS) equations. Solving RANS equations is a tedious process due to the massive computational burden involved, especially when the aerodynamic and structural coupling analysis is required to determine a preliminary aircraft design. In our research, the governing Euler equations for inviscid compressible flow were discretized by equally spatial Cartesian meshes, which greatly facilitated the grid generation. A viscous-inviscid iteration technique provided by the commercial CFD solver MGAERO was used to calculate the lift-to-drag ratio. The integral form of Euler equations for the general control volume can be written as: where vol means the general control volume, → F is the vector of dependent variables, whose transpose is ρ, ρu, ρv, ρw, ρe , ρ is the air density, u, v and w are the velocity components in the Cartesian coordinate system, e is the internal energy of the unit volume,n is the outward normal vector, and → E is the flux vector. By employing the far-field drag extraction theory [30], the shock wave and the induced drag were determined by the decomposition of the total inviscid drag.

Flight Stability
It is quite common to consider momentary disturbances from non-stationary airflow, wakes, movements of passengers, or release of stores, on the aerodynamics of aircraft and stability during the steady flight mode. Such interference can be viewed as a small disturbance, thus the aircraft equations of motion can be decoupled according to the linear small-disturbance theory [31]. In the cruise condition, the normal state-space equation of longitudinal/lateral stability is written as: where A is the system matrix, B is the control matrix, x is the state vector, and u is the control vector. A and B are constructed by the aerodynamic derivatives, the mass property and the aerodynamic coefficients, x is composed by the disturbed flight parameters, and u is composed by the deflections of control surfaces.

Aeroelasticity
The basic equation of static aeroelasticity can be described as: where the subscript α denotes the degree of freedom of the equation, K aa is the structural stiffness matrix, Q aa is the matrix of aerodynamic influence coefficients, M aa is the structural mass matrix, Q ax is the matrix of unit aerodynamic loads, u a is the structural deformation vector, .. u a is the acceleration vector of rigid body motion, u x is the vector of aerodynamic trim parameters (e.g., angle of attack, elevator deflection), which is used to define the deflection of the aerodynamic control surface and the overall rigid motion of the aircraft, P a is the vector of applied loads, and q is the dynamic pressure.
To perform flutter analysis, the equation of the p-k method is used as follows [32]: where h denotes the coordinate of the normal modes, M hh is the modal mass matrix, B hh is the modal damping matrix, K hh is the modal stiffness matrix, V is the flow speed, b is the reference chord length, ρ is the air density, k is the reduced frequency, u h is the modal amplitude vector, Q hh is the generalized aerodynamic force matrix, which is the function of k, and the superscripts R and I represent the real and imaginary parts, respectively. Equation (3) can also be rewritten as: [A − pI]u h = 0, where u h denotes the vector containing modal displacements and velocities, p means the complex eigenvalue.
In terms of displacement interpolation between the structural nodes and the aerodynamic surface grid points, it is determined by the surface spline method [33].

Genetic Algorithm (GA) Based Optimization Strategy
GA is considered a popular tool for locating the global optimum solution [16,[34][35][36]. GA is a population-based evolutionary algorithm and inspired by genetic evolutions. It employs genetic operators of selection, mutation, and crossover to generate the populations for improving the fitness function to achieve its evolutionary strategy. Thus, the mechanism of GA can be described as a combination of an artificial survival of the fittest and genetic operators from nature to iteratively improve the fitness of the population.
The flowchart of GA is shown in Figure 1. The implementation of the GA procedure for solving the multidisciplinary analysis and optimization problem in this paper can be explained as consisting of three steps: Step 1: Generate the initial population using genetic parameters and each individual represents a specific aircraft model.
Step 2: Perform the aforementioned multidisciplinary analysis to evaluate the responses of the interests for the current population.
Step 3: Apply genetic operators to create the next generation of population and repeat Step 2 until the number of generations is met.
To apply GA-based strategy for the aircraft design, a parameterized model of the representative aircraft configuration was developed and then was used to perform the simulations required by GA. The developed model included 19 design parameters depicted in Table 1  To apply GA-based strategy for the aircraft design, a parameterized model of the representative aircraft configuration was developed and then was used to perform the simulations required by GA. The developed model included 19 design parameters depicted in Table 1

Genetic Programming Aided Optimization of Conceptual Aircraft Design
Based on the multidisciplinary analysis described in Section 2, a preliminary aircraft design can be conducted using conventional methods. To overcome the low computational efficiency of traditional sequence approaches, multidisciplinary analysis and optimization (MAO) was used to apply the above theoretical principles in a multidisciplinary context across the disciplines of aerodynamics, structure, aircraft flight dynamics, etc. for efficient aircraft designs in this paper. However, the aircraft multidisciplinary design is actually a multi-parameter optimization problem, it is very difficult to predict and evaluate the influence of design variables on aircraft performance because there are many parameters. Based on this discussion, genetic programming, as a new artificial intelligence technique, is proposed to build mathematical expressions representing the relationship between the input parameters and output responses of the interests. Once this metamodeling process is completed, the optimal solution describing the aircraft configurations will be obtained with less computational cost by replacing the full model with a GP metamodel. The process of this metamodel-based MAO for the optimal design consisted of four steps and they were: Design of experiments 2.
Optimization formulation, design variables and constraints These four steps were integrated to perform the multi-parameter, multidisciplinary optimization of the conceptual aircraft design and explained in the following sub-sections.

Design of Experiments
Design of experiments (DoE) is a statistical technique to simultaneously study the effect of multiple variables with high efficiency and it can successfully deal with design optimization projects of various engineering applications [37]. A survey of different kinds of DoE is given by Simpson et al. [38]. However, the quality of the metamodel strongly depends on an appropriate choice of DoE type and sampling size. To improve the quality of the conventional random Latin hypercube DoE, a uniform Latin hypercube DoE based on the use of the Audze-Eglais optimality criterion [39], was employed in this paper. The main principles in this approach are as follows: The number of levels of factors (same for each factor) is equal to the number of experiments and for each level, there is only one experiment.
The points corresponding to the experiments are distributed as uniformly as possible in the domain of factors, which is controlled by Equation (5): where P is the number of points, L pq is the Euclidean distance between the points p and q (p q) in the system. According to this extended optimal Latin hypercube design of numerical experiments (DoE), a 1997-point DoE has been developed for the FE simulation to be performed at each point. As an example, the bar chart of the minimum distances between the sampling points is shown in Figure 2 indicating a good uniformity of 250-point DoE.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 16 uniform Latin hypercube DoE based on the use of the Audze-Eglais optimality criterion [39], was employed in this paper. The main principles in this approach are as follows: The number of levels of factors (same for each factor) is equal to the number of experiments and for each level, there is only one experiment.
The points corresponding to the experiments are distributed as uniformly as possible in the domain of factors, which is controlled by Equation (5): where P is the number of points, is the Euclidean distance between the points p and q (p ≠ q) in the system.
According to this extended optimal Latin hypercube design of numerical experiments (DoE), a 1997-point DoE has been developed for the FE simulation to be performed at each point. As an example, the bar chart of the minimum distances between the sampling points is shown in Figure 2 indicating a good uniformity of 250-point DoE.

Metamodel Building by Genetic Programming (GP)
GP methodology is a systematic way of selecting a structure of high-quality global approximations, which can reflect the relationship between input parameters and output responses of the interests [24]. As other nature-inspired methods, a population is generated in GP by three most important genetic operators: reproduction, mutation and crossover. The initial population is constructed by a blind random search of the space defined by the problem. The formula obtained by

Metamodel Building by Genetic Programming (GP)
GP methodology is a systematic way of selecting a structure of high-quality global approximations, which can reflect the relationship between input parameters and output responses of the interests [24].
As other nature-inspired methods, a population is generated in GP by three most important genetic operators: reproduction, mutation and crossover. The initial population is constructed by a blind random search of the space defined by the problem. The formula obtained by each individual in a population is composed of a random combination of functions (mathematical operators, such as square and addition) and terminals (variables and constants [40,41]) and these mathematical approximations are actually computer programs with a tree structure. Selection of the structure of an analytical expression is a problem of empirical model building. A tree structure-based typical program, representing the expression (x 1 × x 2 − x 3 ) 1/2 , is shown in Figure 3.

Metamodel Building by Genetic Programming (GP)
GP methodology is a systematic way of selecting a structure of high-quality global approximations, which can reflect the relationship between input parameters and output responses of the interests [24]. As other nature-inspired methods, a population is generated in GP by three most important genetic operators: reproduction, mutation and crossover. The initial population is constructed by a blind random search of the space defined by the problem. The formula obtained by each individual in a population is composed of a random combination of functions (mathematical operators, such as square and addition) and terminals (variables and constants [40,41]) and these mathematical approximations are actually computer programs with a tree structure. Selection of the structure of an analytical expression is a problem of empirical model building. A tree structure-based typical program, representing the expression (x1 × x2 − x3) 1/2 , is shown in Figure 3.
These randomly generated approximations are general and hierarchical, varying in size and shape. The main goal of GP is to solve a problem by searching a highly fit individual computer program in the space of all possible programs. GP methodology tries to find near-global solutions by keeping many solutions that may potentially be close to minima (local or global). A schematic description of the GP technique is sketched in Figure 4 to demonstrate the process of metamodel building. First, initializing the population was performed to create various tree structures. Followed by parameter insertion to generate a family of mathematical expressions with coefficients, optimization by sequential quadratic programming (SQP) algorithm was conducted to determine the optimal locations and numbers of coefficients (or numerical terminals). After tuning, the individual that had the lowest error was chosen to represent the original individual. Finally, the fitness function was evaluated to assess the quality of each individual until the termination criterion was satisfied throughout the iterations of the whole process. More details about the procedure of GP methodology have been introduced in References [24,25,42].
To encourage the evolution of smooth mathematical expressions, the fitness values G(i, t) of individual i at generation t have been defined as a weighted sum of different terms or objectives, following an approach used for multi-objective optimization in evolution-based algorithms: where G 1 is the root mean square error (RMSE) of the i-th individual in the t-th generation evaluated on the given data set, divided by the average RMSE of the archive individuals at the previous generation, N is the total number of training data, y i is the actual value, y i is the estimated value; G 2 is the square of the number of numerical coefficients (parameters) present in the individual; G 3 is the number of operations not defined (i.e., division by zero) in the individual at any of the sample points; G 4 is the number of nodes that the individual is made of and a 1 , a 2 , a 3 and a 4 are weighting factors (that add up to 1) determined by the exhaustive testing and tuning of the GP algorithm. To robustly solve multi-criterion optimization problems, a Pareto set method can be applied to select optimal weighting factors [43]. In this paper, the weighting factors were a 1 = 0.8989, a 2 = 0.001, a 3 = 0.1, and a 4 = 0.0001.
building. First, initializing the population was performed to create various tree structures. Followed by parameter insertion to generate a family of mathematical expressions with coefficients, optimization by sequential quadratic programming (SQP) algorithm was conducted to determine the optimal locations and numbers of coefficients (or numerical terminals). After tuning, the individual that had the lowest error was chosen to represent the original individual. Finally, the fitness function was evaluated to assess the quality of each individual until the termination criterion was satisfied throughout the iterations of the whole process. More details about the procedure of GP methodology have been introduced in References [24,25,42]. To encourage the evolution of smooth mathematical expressions, the fitness values G( , ) of individual i at generation t have been defined as a weighted sum of different terms or objectives, following an approach used for multi-objective optimization in evolution-based algorithms: Two fit indices were computed to evaluate the accuracy of predictive models constructed by GP: one was the aforementioned RSME, and the coefficient of determination (R 2 ) was the other fit criterion used to control the predictive performance of the metamodels: where N is the total number of training data, y i is the actual value, y i is the estimated value, and y is the mean of the values y i . Obviously, the predicted model will be quite perfect if R 2 is equal to one and RMSE is zero. It is noted that there are 17 mathematical operators developed for GP-based metamodel building. They are {+, −, *, /, square, cube, constant, sine, cosine, sinh, cosh, tanh, exponential, negative exponential, natural logarithm, absolute, reciprocal}.

Parameterized Finite Element Model
A representative beam-frame wing structure and its finite element model for static aeroelasticity and flutter analyses are shown in Figure 5. The geometric configuration of this high-aspect-ratio wing model was decided by seven parameters in Figure 5a, which are the inboard/outboard taper ratios (η in , η out ), the inboard/outboard aspect ratios (λ in , λ out ), the sweep angle of the leading edge (χ), the inboard/outboard dihedral angles (ψ in , ψ out )). There are three more parameters (p r , p k , p t ) for determining the elastic axis positions at the root, the kink and the tip, respectively. The wing root length (c r ) was kept constant. The wing structural stiffness properties were modelled by a variable cross-section beam from the root to the tip and the structural mass was represented by lumped mass elements. Material properties of aluminum were assigned to the primary beam of the wing, and the components such as the body, tail and fin, were modelled using rigid elements. The wing lifting surface was modelled with 282 aerodynamic elements. A fixed maximum take-off weight of 84.5 tons and cruise Mach of 0.785 at a height of 11,200 m were applied to perform the analysis. The Reynolds number was 5.71 million based on the mean aerodynamic chord. The root chord remains a constant length with 6.17 m, and the reference area was varied with the geometric design variables.
A representative beam-frame wing structure and its finite element model for static aeroelasticity and flutter analyses are shown in Figure 5. The geometric configuration of this high-aspect-ratio wing model was decided by seven parameters in Figure 5a, which are the inboard/outboard taper ratios (η , η ), the inboard/outboard aspect ratios (λ , λ ), the sweep angle of the leading edge (χ), the inboard/outboard dihedral angles ( ψ , ψ )). There are three more parameters ( p , p , p ) for determining the elastic axis positions at the root, the kink and the tip, respectively. The wing root length (c ) was kept constant. The wing structural stiffness properties were modelled by a variable cross-section beam from the root to the tip and the structural mass was represented by lumped mass elements. Material properties of aluminum were assigned to the primary beam of the wing, and the components such as the body, tail and fin, were modelled using rigid elements. The wing lifting surface was modelled with 282 aerodynamic elements. A fixed maximum take-off weight of 84.

Optimization Formulation, Design Variables and Constraints
Designing aircraft structures is a challenging procedure because the weight as an objective should be minimized and the stiffness provided should supply the required structural rigidity. The structural stiffnesses are expressed as (flexural rigidity in the x-direction), (flexural rigidity in the y-direction), GJ (torsional rigidity), and the mass is denoted by , where is the density of the air and represents a fixed cross-section. For the constants (Young's modulus) and (polar second moment of area), the stiffness is proportional to the mass: where is the reference length of the wing. Thus, in the aircraft conceptual design phase, the optimization formation can be defined as: Objective: minimize ∑(( + + ) ) Design Constraints: Four aeroelastic constraints and one aerodynamic constraint: The maximum z-direction displacement at the wingtip was less than 7% of the magnitude at the semi-span. The wing deformation of the large aircraft will still obey the linear theory of elasticity in the cruise flight condition; The maximum torsion at the wingtip was less than 2.14°. The elastic torsional deformation that might lead to a decrease of the local angle of attack; The aileron efficiency was greater than 60%; The flutter speed at sea level was greater than 320 m/s;

Optimization Formulation, Design Variables and Constraints
Designing aircraft structures is a challenging procedure because the weight as an objective should be minimized and the stiffness provided should supply the required structural rigidity. The structural stiffnesses are expressed as EI xx (flexural rigidity in the x-direction), EI yy (flexural rigidity in the y-direction), GJ (torsional rigidity), and the mass is denoted by ρS, where ρ is the density of the air and S represents a fixed cross-section. For the constants E (Young's modulus) and J (polar second moment of area), the stiffness is proportional to the mass: where l is the reference length of the wing. Thus, in the aircraft conceptual design phase, the optimization formation can be defined as: Objective : minimize EI xx + EI yy + GJ l Design Constraints: Four aeroelastic constraints and one aerodynamic constraint: The maximum z-direction displacement at the wingtip was less than 7% of the magnitude at the semi-span. The wing deformation of the large aircraft will still obey the linear theory of elasticity in the cruise flight condition; The maximum torsion at the wingtip was less than 2.14 • . The elastic torsional deformation that might lead to a decrease of the local angle of attack; The aileron efficiency was greater than 60%; The flutter speed at sea level was greater than 320 m/s; The total drag was less than the drag obtained by the baseline design. Design Variables: The lower and upper bounds of design variables were determined by 90% and 110% of the baseline values [44], respectively, which are shown in Table 1.
The optimization flowchart shown in Figure 6 demonstrates GP metamodel-assisted multidisciplinary optimization process for an aircraft conceptual design. Firstly, the Latin optimal hypercube design of experiments was used to determine the best distribution of sampling points for metamodel building. According to these design points, various jig shape models of the representative wing were generated for aeroelastic and stability analyses. Thus, the training data for GP metamodel building was obtained. By comparison of the predicted values from metamodels with the true results, more sampling points will be used to improve the accuracy of metamodels in the next iteration if the accuracy quality is not satisfactory. Finally, parametric optimization on the metamodels using SQP was performed to seek the optimal aircraft design.

Case Study
Both GA and GP techniques were applied to solve the above MAO problem for an aircraft conceptual design. Six responses related to stiffness, wingtip displacement, wingtip torsion, aileron efficiency, flutter and drag were built using GP. A cross-validation check of GP and GA approaches was carried out by comparison of the obtained optimal designs.
To construct the metamodels using GP, 1997 sampling points were selected as training data by the optimal Latin hypercube DoE to maximally represent information in the whole design space. Two fit indices were computed to evaluate the accuracy of predictive models constructed by GP, namely root mean square error (RSME) and coefficient of determination (R2) in Equations (6) and (9), respectively. In this paper, main parameters used in GP metamodel building are shown in Table 2.
As an example, the mathematical expression for the wingtip torsion is given as below:

Case Study
Both GA and GP techniques were applied to solve the above MAO problem for an aircraft conceptual design. Six responses related to stiffness, wingtip displacement, wingtip torsion, aileron efficiency, flutter and drag were built using GP. A cross-validation check of GP and GA approaches was carried out by comparison of the obtained optimal designs.
To construct the metamodels using GP, 1997 sampling points were selected as training data by the optimal Latin hypercube DoE to maximally represent information in the whole design space. Two fit indices were computed to evaluate the accuracy of predictive models constructed by GP, namely root mean square error (RSME) and coefficient of determination (R2) in Equations (6) and (9), respectively. In this paper, main parameters used in GP metamodel building are shown in Table 2. As an example, the mathematical expression for the wingtip torsion is given as below: where X 1 to X 19 are design variables detailed in Table 1. It should be noted that results of two fit indices, for example, RSME = 1.27 × 10 −2 and R 2 = 0.995 for the predicted wingtip torsion response, indicate the metamodels by GP are quite accurate to estimate the nonlinear responses of the interests. Furthermore, design variables regarding horizontal bending stiffness and torsional stiffness have slight effects on the wingtip torsion, and this can be proved by sensitivity analysis performed in Section 5.2.
To validate the accuracy of generated metamodels, a validation data set of 100 sampling points was also used to check the overall performance capacity of the metamodels. In this paper, the relative error was used to measure the discrepancy between the actual response and the predicted model, e = r a − r p /r a (13) where r a is the actual response value and r p is the prediction value. The relative errors of 100 validation points are given in Figure 7 and their relative mean errors for all six responses were 0.083%, 0.008%, 0.024%, 0.119%, 0.106%, and 0.598%. On average, the predictions of all the responses by GP technique were quite accurate. However, the most limitation of GP methodology was its complexity, discontinuity, and consumption of the model construction. The higher the dimension of the design space was, the more complicated the constructed GP expression was, the higher the number of step discontinuities of the response functions occurred, and more RAM was required. For the design space composed of 19 variables, it took more than three days for the GP methodology to construct the six responses on a personal computer with eight cores (E5-2670, 2.6 GHz) and 24 G RAM.  On average, the predictions of all the responses by GP technique were quite accurate. However, the most limitation of GP methodology was its complexity, discontinuity, and consumption of the model construction. The higher the dimension of the design space was, the more complicated the constructed GP expression was, the higher the number of step discontinuities of the response functions occurred, and more RAM was required. For the design space composed of 19 variables, it took more than three days for the GP methodology to construct the six responses on a personal computer with eight cores (E5-2670, 2.6 GHz) and 24 G RAM.

Optimal Designs by GP and GA
Metamodels were used to compute the responses of the interests during the MAO process for the representative aircraft conceptual design. The comparison of the optimal result by GP-based and GA-based optimizations and finite element validations is given in Table 3 and the corresponding design variables are shown in Table 4. The genetic parameters used by GA include Population 100, Generation 16, Tournament Selection, Crossover 0.7, and Mutation 0.4.

Optimal Designs by GP and GA
Metamodels were used to compute the responses of the interests during the MAO process for the representative aircraft conceptual design. The comparison of the optimal result by GP-based and GA-based optimizations and finite element validations is given in Table 3 and the corresponding design variables are shown in Table 4. The genetic parameters used by GA include Population 100, Generation 16, Tournament Selection, Crossover 0.7, and Mutation 0.4.  Since the optima obtained by GP metamodel-based and GA-based optimizations satisfy the design constraints, both results are feasible solutions. It is noted that the superiority of the GP method over GA for MAO was demonstrated by the overall performance of the optimal design except for the drag response shown in Table 3. The stiffness by GP based optimization reduced from 0.806 to 0.535 by 33.6%, while the drag just deteriorated by 0.77% and no constraints were violated. In addition, it is observed that the predictions of responses agree well with finite element validation due to the high-accuracy models.
The most remarkable differences among the two cruise shapes shown in Figure 8 are that the optimum by GA-based method had the largest sweep angle of the leading edge, while the optimal design using GP-based method had a largest outboard spanwise length. The increase of the sweep angle contributed to the decrease of the local shock wave drag and the higher aspect ratio led to the induced drag in reduction. The induced drag and the local shock wave drag are two significant types of a large aircraft flight drag under the cruise Mach of 0.785. The increase of the sweep angle also results in the increase of the vertical bending-torsional coupling effect as well as the structural weight. It is the reason why the optimum by GA method had a heavier structural weight and a larger wingtip torsional deformation than the optimal design by GP method.
Since the optima obtained by GP metamodel-based and GA-based optimizations satisfy the design constraints, both results are feasible solutions. It is noted that the superiority of the GP method over GA for MAO was demonstrated by the overall performance of the optimal design except for the drag response shown in Table 3. The stiffness by GP based optimization reduced from 0.806 to 0.535 by 33.6%, while the drag just deteriorated by 0.77% and no constraints were violated. In addition, it is observed that the predictions of responses agree well with finite element validation due to the highaccuracy models.
The most remarkable differences among the two cruise shapes shown in Figure 8 are that the optimum by GA-based method had the largest sweep angle of the leading edge, while the optimal design using GP-based method had a largest outboard spanwise length. The increase of the sweep angle contributed to the decrease of the local shock wave drag and the higher aspect ratio led to the induced drag in reduction. The induced drag and the local shock wave drag are two significant types of a large aircraft flight drag under the cruise Mach of 0.785. The increase of the sweep angle also results in the increase of the vertical bending-torsional coupling effect as well as the structural weight. It is the reason why the optimum by GA method had a heavier structural weight and a larger wingtip torsional deformation than the optimal design by GP method.   Figure 9 shows the variations of vertical, horizontal and torsional stiffnesses along the wingspan. Since the vertical stiffness is the most important parameter to reflect the aeroelastic character, it is reasonable to identify the slight variation of the vertical bending performance by these two methods. The GP method helped to remarkably reduce the horizontal stiffness in the region near the wing root and decrease the torsional stiffness shown in Figure 9b,c, respectively.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 13 of 16 Figure 9 shows the variations of vertical, horizontal and torsional stiffnesses along the wingspan. Since the vertical stiffness is the most important parameter to reflect the aeroelastic character, it is reasonable to identify the slight variation of the vertical bending performance by these two methods. The GP method helped to remarkably reduce the horizontal stiffness in the region near the wing root and decrease the torsional stiffness shown in Figure 9b,c, respectively. The transient behavior of the state variables in Equation (2) for the fourth-order model is shown in Figure 10. It can be found that the disturbance quantities ultimately vanish and the responses converge towards a stable steady state. The natural response to perturbations consists of two damped modes: the short period and heavily damped mode (a and b) and the long period and lightly damped mode (phugoid mode: c and d). The results are quite typical for a large aircraft, and the stability of uncontrolled motion considered as a design constraint can be guaranteed during the optimization process. The transient behavior of the state variables in Equation (2) for the fourth-order model is shown in Figure 10. It can be found that the disturbance quantities ultimately vanish and the responses converge towards a stable steady state. The natural response to perturbations consists of two damped modes: the short period and heavily damped mode (a and b) and the long period and lightly damped mode (phugoid mode: c and d). The results are quite typical for a large aircraft, and the stability of uncontrolled motion considered as a design constraint can be guaranteed during the optimization process. The transient behavior of the state variables in Equation (2) for the fourth-order model is shown in Figure 10. It can be found that the disturbance quantities ultimately vanish and the responses converge towards a stable steady state. The natural response to perturbations consists of two damped modes: the short period and heavily damped mode (a and b) and the long period and lightly damped mode (phugoid mode: c and d). The results are quite typical for a large aircraft, and the stability of uncontrolled motion considered as a design constraint can be guaranteed during the optimization process.

Sensitivity Analysis
Finite difference method was used to determine the derivatives of the responses with respect to each design variable. The variation of each design variable was set at 5% of its design interval indicated in Table 1. The occurrence of flutter was determined by the additional damping at a set of

Sensitivity Analysis
Finite difference method was used to determine the derivatives of the responses with respect to each design variable. The variation of each design variable was set at 5% of its design interval indicated in Table 1. The occurrence of flutter was determined by the additional damping at a set of speeds, so the sensitivity of the flutter speed was not taken into consideration. In Figure 11, the results show that sensitivity analysis is used to explore how the variations of design variables have an impact on the response such as the weight, wingtip displacement, wingtip torsion, aileron efficiency, and drag. It can be also found that the more important the design variable was to the response, the larger the size of the slice displayed in the pie chart. X depicts the design variable, which is detailed in Table 1. The symbol "-" means a negative correlation. Obviously, the horizontal bending stiffness (X14) had a great influence on the weight. The vertical bending stiffness (X11), the sweep angle of the leading edge (X5), and the outboard aspect ratio (X4) were the top three key parameters to have important impacts on the other responses. Hence, aircraft designers should pay more attention to these parameters in the conceptual phase of designing a large aircraft.
the size of the slice displayed in the pie chart. X depicts the design variable, which is detailed in Table  1. The symbol "-" means a negative correlation. Obviously, the horizontal bending stiffness (X14) had a great influence on the weight. The vertical bending stiffness (X11), the sweep angle of the leading edge (X5), and the outboard aspect ratio (X4) were the top three key parameters to have important impacts on the other responses. Hence, aircraft designers should pay more attention to these parameters in the conceptual phase of designing a large aircraft. Figure 11. Sensitivity of the optimal design obtained by the GP method.

Research Significance
Research on the conceptual design of a large aircraft has revealed the immense benefits of improving overall performance, e.g., safety, comfort, and mitigation on environmental impacts. However, conventional sequence approaches can hardly adapt to the design of a modern aircraft, due to their iterative design methodology and low computational efficiency. To address this issue, the metaheuristic algorithm assisted aerostructural optimization method has been proposed in this paper and its superiority over the sequential design methods was examined by a representative beam-frame wing structure. The results give useful insights to tailoring parameters for the effective design of a large aircraft in the conceptual phase and also provide aircraft designers with a wealth of information for design practice.

Conclusions
Taking into account the interest in the integration of various disciplines for solving complex engineering design problems being rapidly growing within industries, GP metamodel-based multidisciplinary analysis and optimization (MAO) strategy has been proposed in this paper to improve the efficiency of an aircraft conceptual design. The correctness of this approach was examined by MAO of a representative beam-frame wing structure and its better strength and level of accuracy for seeking the optimum were demonstrated by comparison with the result obtained using the genetic algorithm technique. The results of the sensitivity analysis indicate that the horizontal bending stiffness had a great impact on the weight and three critical parameters significantly affecting the other responses of interest are determined: the vertical bending stiffness, Figure 11. Sensitivity of the optimal design obtained by the GP method.

Research Significance
Research on the conceptual design of a large aircraft has revealed the immense benefits of improving overall performance, e.g., safety, comfort, and mitigation on environmental impacts. However, conventional sequence approaches can hardly adapt to the design of a modern aircraft, due to their iterative design methodology and low computational efficiency. To address this issue, the metaheuristic algorithm assisted aerostructural optimization method has been proposed in this paper and its superiority over the sequential design methods was examined by a representative beam-frame wing structure. The results give useful insights to tailoring parameters for the effective design of a large aircraft in the conceptual phase and also provide aircraft designers with a wealth of information for design practice.

Conclusions
Taking into account the interest in the integration of various disciplines for solving complex engineering design problems being rapidly growing within industries, GP metamodel-based multidisciplinary analysis and optimization (MAO) strategy has been proposed in this paper to improve the efficiency of an aircraft conceptual design. The correctness of this approach was examined by MAO of a representative beam-frame wing structure and its better strength and level of accuracy for seeking the optimum were demonstrated by comparison with the result obtained using the genetic algorithm technique. The results of the sensitivity analysis indicate that the horizontal bending stiffness had a great impact on the weight and three critical parameters significantly affecting the other responses of interest are determined: the vertical bending stiffness, the sweep angle of the leading edge, and the outboard aspect ratio. It is concluded that GP metamodel-based strategy for solving MAO problems is more efficient and accurate as compared with GA and has demonstrated the capability of becoming a systematic and robust approach to the aircraft design in the conceptual phase.