Metaheuristic Approaches to Solve a Complex Aircraft Performance Optimization Problem

Abstract

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.

1. 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].

Aircraft design is an inherently multidisciplinary, iterative effort. However, conventional design methods weaken the collaborative effect of various disciplines so that aerodynamic, structural, and stability analyses have to be performed in a specific sequence to achieve an optimum of a system. The major limitation of this method relates to an iterative design process and will lead to less efficiency of computations, so it can hardly adapt to the design of a modern aircraft, especially the conceptual design with great potential for innovation in the early stage. Instead, multidisciplinary analysis and optimization (MAO) could take disciplines containing aerodynamics, structure, aircraft flight dynamics, etc. into consideration simultaneously and has been employed in large modern aircraft designs [4,5,6,7,8,9]. In light of its capability to overcome the limitations of the conventional method, the MAO technique with the integration of novel algorithms has been applied to solve complex problems in other industry sectors, such as automobile [10,11], energy [12], mechatronic systems [13], etc. Obviously, a crucial challenge of MAO is its massive computation burden, especially for computational fluid dynamics (CFD) simulations based on Reynolds-averaged Naiver-Stokes (RANS) equations. Though RANS method can simulate intricate turbulence and shedding vortex, the steady flight conditions make up most of the flight time of a large aircraft, hence some viscous-inviscid iteration methods are more suitable for cruise conditions with moderate flow separation [14]. As the disturbance in the cruise condition is weak compared with the steady-state, the small-disturbance theory is valid for the stability analysis.

Although a gradient-based algorithm has the advantages of rapidity and effectiveness for solving an optimization problem, it is apt to converge to a local optimal solution and sometimes, it fails in special cases where calculating the partial derivatives is hard or the derivatives are not continuous. Considering these facts, the evolutionary algorithms are suitable to seek the globally optimal solution. Genetic algorithm (GA) is based on Darwin’s theory and is used most widely as a global algorithm for this purpose [15,16]. However, function evaluations in GA are the most time-consuming part as they are usually executed by finite element analysis (FEA) or CFD with high computational demands, for example, a single function evaluation could be of the order of hours or days. Hence, the computational burden of FEA or CFD for design evaluations is the main limitation of GA for solving multidisciplinary optimization problems.

To reduce the computational burden and benefit to a great extent of better exploration and exploitation in the design space, metamodel-based optimization has recently found widespread use owing to its promising potential [17,18,19,20,21,22,23]. As a response surface method, genetic programming (GP) is a systematic way of selecting a structure of high-quality global approximations [24]. The program can represent an empirical model to be used for approximation of a response function, and the explicit expression [25], or even algorithms [26], can be entirely attained for optimization. With the high accuracy, approximation models generated by GP have already been used for solving various design optimization problems [18,21,27,28,29].

In light of the exploitation of the immense potential from a large aircraft conceptual design process and the advantages of GP methodology in support of MAO, an integrated GP metamodel-based optimization method to design a full-scale aircraft subject to aerodynamic, structural and stability constraints is proposed in this paper. Static aeroelasticity and flutter analyses are performed with the MSC/Nastran doublet-lattice method. The stability of flight dynamics is evaluated by the linear small-disturbance method. The lift-to-drag ratio is estimated by a viscous-inviscid iteration technique implemented in the CFD solver MGAERO. To address this MAO problem, a two-step procedure for designing a conceptual aircraft is proposed. First, GP methodology was applied to generate metamodels based on optimal Latin hypercube design of experiments (DoE). With these metamodels, multidisciplinary analysis can be efficiently executed due to the explicit mathematical approximations of the predictive models. Subsequently, parametric optimization was performed for the minimum stiffness of the wing structure subject to aerodynamic, aeroelastic and stability constraints. To demonstrate the strength and accuracy of GP methodology for solving complex design problems, both GA and GP methods were employed to find the optimal solutions for a representative beam-frame wing structure. Finally, sensitivity analysis was conducted to identify the critical input parameters that affect aircraft performance in the conceptual design phase.

2. 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.

2.1. 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:

$$\frac{\partial }{\partial \mathrm{t}}{\iiint }_{vol}\overrightarrow{F}dv+{\iint }_s\overrightarrow{E}·\widehat{n}dS=0$$

(1)

where $vol$ means the general control volume, $\overrightarrow{F}$ is the vector of dependent variables, whose transpose is $\left\{\rho ,\rho u,\rho v,\rho w,\rho e\right\}$, $\mathsf{\rho }$ 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, $\widehat{n}$ is the outward normal vector, and $\overrightarrow{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.

2.2. 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:

$$\dot{x}=Ax+Bu$$

(2)

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.

2.3. Aeroelasticity

The basic equation of static aeroelasticity can be described as:

$${\mathrm{K}}_{aa}-\overline{q}{\mathrm{Q}}_{aa}{u}_a+M_{aa}{\ddot{u}}_a=\overline{q}{Q}_{ax}{u}_x+P_a$$

(3)

where the subscript $\mathsf{\alpha }$ denotes the degree of freedom of the equation, ${\mathrm{K}}_{aa}$ is the structural stiffness matrix, ${\mathrm{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, ${\ddot{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 $\overline{q}$ is the dynamic pressure.

To perform flutter analysis, the equation of the p-k method is used as follows [32]:

$$\left[{\left(\frac{V}{b}\right)}^2{p}^2{M}_{hh}+\frac{V}{b}p{B}_{hh}+K_{hh}-\frac{1}{2}\rho V^2\left(Q_{hh}{}^R+\frac{P}{K}{Q}_{hh}{}^I\right)\right]u_h=0$$

(4)

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, $\rho$ 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: $\left[A-pI\right]{\overline{u}}_h=0$, where ${\overline{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].

3. 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. Bounds of 19 design variables.

Design Variables	Lower Bound	Upper Bound
X1, Inboard taper ratio	1.457	1.781
X2, Outboard taper ratio	3.658	4.471
X3, Inboard aspect ratio	0.937	1.146
X4, Outboard aspect ratio	2.509	3.067
X5, Sweep angle (degree)	24.651	30.130
X6, Inboard dihedral angle (degree)	4.346	5.312
X7, Outboard dihedral angle (degree)	0.645	0.789
X8, Axis position (root)	0.306	0.374
X9, Axis percent (tip)	0.338	0.413
X10, Axis percent (kink)	0.389	0.476
X11, Vertical bending (a)	0.600	2.000
X12, Vertical bending (b)	0.600	2.000
X13, Vertical bending (c)	0.600	2.000
X14, Horizontal bending (a)	0.600	2.000
X15, Horizontal bending (b)	0.600	2.000
X16, Horizontal bending (c)	0.600	2.000
X17, Torsional (a)	0.800	1.200
X18, Torsional (b)	0.800	1.200
X19, Torsional (c)	0.800	1.200

in Section 4.3 and Section 4.4.

4. 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
• Metamodel building by genetic programming (GP)
• Parameterized finite element model
• 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.

4.1. 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):

$$U={\displaystyle \sum }_{p=1}^P{\displaystyle \sum }_{q=p+1}^P\frac{1}{L_{pq}^2}\to min$$

(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.

4.2. 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₁ × x₂ − x₃)^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 $\mathrm{G}\left(i,t\right)$ 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:

$$RMSE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(y_i-y_i^{\prime }\right)}^2}$$

(6)

$$\mathrm{G}\left(i,t\right)=a_1{G}_1\left(i,t\right)+a_2{G}_2\left(i,t\right)+a_3{10}^6{G}_3\left(i,t\right)+a_4{G}_4\left(i,t\right),$$

(7)

$$a_1+a_2+a_3+a_4=1$$

(8)

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^{\prime }$ 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$.

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²) was the other fit criterion used to control the predictive performance of the metamodels:

$$R^2=1-\frac{{\sum }_{i=1}^N{\left(y_i-y_i^{\prime }\right)}^2}{{\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}$$

(9)

where N is the total number of training data, $y_i$ is the actual value, $y_i^{\prime }$ is the estimated value, and $\overline{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}.

4.3. 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 (${\mathsf{\eta }}_{in},{\mathsf{\eta }}_{out}$), the inboard/outboard aspect ratios (${\mathsf{\lambda }}_{in},{\mathsf{\lambda }}_{out}$), the sweep angle of the leading edge ($\mathsf{\chi }$), the inboard/outboard dihedral angles (${\mathsf{\psi }}_{in},{\mathsf{\psi }}_{out}$)). There are three more parameters (${\mathrm{p}}_r,{\mathrm{p}}_k$,${\mathrm{p}}_t$) for determining the elastic axis positions at the root, the kink and the tip, respectively. The wing root length (${\mathrm{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.

4.4. 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 $E{I}_{xx}$ (flexural rigidity in the x-direction), $E{I}_{yy}$ (flexural rigidity in the y-direction), GJ (torsional rigidity), and the mass is denoted by$\rho S$, where $\rho$ 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:

$$\left(E{I}_{xx}+E{I}_{yy}+GJ\right)l\propto \rho Sl$$

(10)

where $l$ is the reference length of the wing.

Thus, in the aircraft conceptual design phase, the optimization formation can be defined as:

$$\mathrm{Objective}:\mathrm{minimize}\sum \left(\left(E{I}_{xx}+E{I}_{yy}+GJ\right)l\right)$$

(11)

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.

5. 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. Main evolution parameters used in GP metamodel building.

Parameters	Values
Population size	500
Generations	100
Number of fitness cases	1997
Initialisation	50% full + 50% grow
Maximum depth initialised trees	4
Minimum depth initialised trees	2
Maximal depth for a tree	50
Reproduction rate	20%
Crossover rate	40%
Mutation rate	40%
Termination of the evolution if RMSE is less than	1.0 × 10−7

. As an example, the mathematical expression for the wingtip torsion is given as below:

$$\begin{gathered} F_{torsion}=-0.238478+0.04X_1-0.71435X_{10}+0.14171X_{11}-\frac{0.252627X_{11}}{X_{12}}- \\ 0.0317857X_{13}+\frac{0.12126X_{13}}{X_{12}}-762.093X_{15}+0.197786X_{18}-0.01X_2+0.682686X_3- \\ 0.297X_4+\frac{0.511492X_4}{X_{12}}-0.014X_5-\frac{0.007X_5}{X_{11}}+\frac{0.00172534X_{11}{X}_5}{X_{12}{}^3}+\frac{0.00003781X_4{X}_5{}^2}{X_{11}{X}_{12}}+ \\ \frac{0.0071\left(-\frac{16.6482X_{13}}{X_{12}}+0.106X_{15}-110.845X_3+1.73842X_5\right)}{X_{12}}-0.00245223X_7-0.261464X_{18}{X}_9+ \\ (0.194348X_{12}+762.093X_{15}+\frac{0.000126113X_4{X}_5{}^2}{X_{11}{X}_3})\mathrm{Tan}\mathrm{h}\left(64742.1X_7{}^3\right). \end{gathered}$$

(12)

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⁻² and R² = 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=\left|r_a-r_p\right|∕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.

5.1. 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. Comparison of responses.

Responses	GP	Validation	GA	Constraint
Stiffness	0.533	0.535	0.806	-
Displacement	5.61%	5.65%	6.16%	<7%
Torsion (deg)	1.39	1.44	2.05	<2.14
Efficiency	0.677	0.680	0.655	>0.6
Flutter (m/s)	>320	>320	>320 m/s	>320
Drag (N)	65,158.2	63,090.8	64,661.2	<65,438

and the corresponding design variables are shown in

Table 4. Geometric parameters and elastic axis position parameters.

Geometric Design Variables	GP	GA
Inboard taper ratio	1.74	1.60
Outboard taper ratio	4.25	3.52
Inboard aspect ratio	0.96	0.91
Outboard aspect ratio	2.90	2.83
Sweep angle (deg)	25.63	29.90
Inboard dihedral angle (deg)	5.18	4.50
Outboard dihedral angle (deg)	0.69	0.82
Axis position (root)	0.31	0.33
Axis percent (tip)	0.41	0.45
Axis percent (kink)	0.46	0.44

. 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.

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.

5.2. 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.

6. Summary

6.1. 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.

6.2. 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.