Abstract
Optimisation is aimed at enhancing aircraft design by identifying the most promising wing planforms at the early stage while discarding the least performing ones. Multiple disciplines must be taken into account when assessing new wing planforms, and a model-based framework is proposed as a way to include mass estimation and longitudinal stability alongside aerodynamics. Optimisation is performed with a particle swarm optimiser, statistical methods are exploited for mass estimation, and the vortex lattice method (VLM) with empirical corrections for transonic flow provides aerodynamic performance. Three surrogates of the aerodynamic model are investigated. The first one is based on radial basis function (RBF) interpolation, and it relies on a precomputed database to evaluate the performance of new wing planforms. The second one is based on an artificial neural network, and it needs precomputed data for a training step. The third one is a hybrid model which switches automatically between VLM and RBF, and it does not need any precomputation. Its switching criterion is defined in an objective way to avoid any arbitrariness. The investigation is reported for a test case based on the common research model (CRM). Reference results are produced with the aerodynamic model based on VLM for two- and three-objective optimisations. Results from all surrogate models for the same benchmark optimisation are compared so that their benefits and limitations are both highlighted. A discussion on specific parameters, such as number of samples for example, is given for each surrogate. Overall, a model-based implementation with a hybrid model is proposed as a compromise between versatility and an arbitrary level of accuracy for wing early-stage design.
1. Introduction
Multidisciplinary optimisation represents the new frontier of aircraft design [1]. It promises to reduce cost and time for developing new products by providing better configurations at an earlier stage [1, 2]. This can be achieved with gradient-free techniques which do not rely on differentiable objective functions [3] and can be applied to aircraft design as shown in [4] by taking into account a large number of design variables. Their main limitation is the repeated evaluation of the objective function which can lead to high computational cost. Regarding aircraft, an example of gradient-free optimisation for early-stage design is available in [5] where genetic algorithms and multifidelity aerodynamics are exploited to identify the optimal location for propellers.
Remedies have been proposed in literature to lower the computational cost of optimisation by exploiting surrogate models, also known as metamodels. In general, surrogate modelling is widely employed in engineering design loops to speed up the process since they reduce the time for performance assessment [6]. It can be applied to any branch of engineering [7, 8], and, for example, application to aerofoil shape optimisation is described in [9]. A summary of surrogate modelling methods for evolutionary optimisation is available in [10, 11] while an application to aircraft conceptual design is given in [12]. Regarding aerodynamics, multiple techniques have been proposed in literature to build surrogates for both steady and unsteady flows. They usually rely on methods such as Kriging interpolation [13, 14], radial basis function (RBF) interpolation [15], proper orthogonal decomposition [16], eigenvalue decomposition [17], and artificial neural network [18, 19]. An additional class of models, called herein hybrid models, are able to exploit information coming from both the original model and its surrogate. Applications of this approach to include multiple fidelity levels can be found in literature for wing design optimisation [20] as well as airfoil design [9, 15] using Kriging interpolation. In such cases, information from the original model is employed to locate the next sampling point which improves the prediction of the underlying Kriging and a reduction up to one order of magnitude in computational time is achieved [20]. Hybrid models can exploit information in a dynamic way too by selecting, for example, the basis function of the surrogate during the optimisation process as proposed in [21]. Another approach is a first estimation of the error with the surrogate model for the new point in the parametric space, and, if the error is larger than a given threshold, the assessment of the same new point is performed with the original model. This approach is employed in [22] to speed up turbomachinery design, and the threshold value is chosen arbitrarily. It affects the number of evaluations performed with the surrogate and, in the end, the level of error in the final results. However, an objective way to switch between models, which does not rely on engineering experience, is still needed. In addition, a comparison of multiple surrogates applied to the optimisation of transonic wing planforms is missing in literature. This paper addresses both needs by suggesting an objective criterion for the switching between the surrogate and its original model as well as benchmarking multiple surrogates.
Besides aerodynamics, additional disciplines must be included to perform a successful wing planform optimisation and a summary of possible architectures is presented in [23]. An example of a multidisciplinary framework is given in [24] where a genetic algorithm is adopted to optimise general aviation aircraft at the conceptual design stage. In addition, the industrial designer requires a flexible framework which incorporates existing tools [25] and possibly interchange them at any time. Frameworks based on model-based principles [26] satisfy both requirements, i.e., the possibility of multifidelity aerodynamics levels as well as software units for distinct disciplines. Although model-based principles have been available in literature since the beginning of the 90s, their application to aircraft design is a recent development. For example, a model-based framework is employed to study the performance of general aviation aircraft in [25]. Besides aerodynamics, the system includes models for physical parts such as engines and fuel tanks.
In this paper, surrogate modelling for multidisciplinary optimisation of wing planform for transonic aircraft at early-stage design is investigated. In particular, three surrogate models are discussed, specifically RBF interpolation, neural network, and a hybrid approach which switches automatically between surrogate and its original model. Regarding the latter, a novel, objective criterion for the switching is proposed. The investigation is made possible by the versatile model-based framework which is suggested as ideal architecture for multidisciplinary and multifidelity problems. The paper proceeds with a description of framework and surrogate modelling methods in Section 2. The infrastructure is presented in detail with focus on the various models employed for this paper. Surrogate modelling techniques, normalisation of data, and error quantification are described in that section too. The wing planform optimisation problem is formulated in Section 3. A generic transonic wing planform is employed to formulate the problem so that the optimisation is not tied to any specific test case. Besides the mathematical form, a graphical interpretation of both parameters and constraints is provided. Results from the surrogate modelling investigation are given in Section 4 for the particular case of optimising the common research model (CRM) wing. First, results are produced with the model-based framework using the aerodynamic model based on VLM for two- and three-objective optimisations. Such results are the benchmark for the performance of surrogate models. Secondly, surrogates based on RBF or neural network are assessed in Section 4.1. A comparative analysis is provided between their results and the reference solution. Thirdly, the application of the hybrid model is described in Section 4.2 where numerical details are given about the criterion for switching between the surrogate and its original model. A final comparison is provided between all surrogates and the reference. Conclusions and a description of future developments are given in Section 5.
2. Model-Based Framework and Surrogate Modelling
The optimisation loop is implemented using a model-based approach [26], and the implementation is coded in an object-oriented way [27]. Interfaces declare input and output capabilities for each distinct class of models. In turn, a model belonging to a given class must meet the requirements set by the interface. Communication between models is only allowed through channels defined by interfaces. As a result, the software is composed of models and definitions of their interactions. Each model can use one or more tools to undertake its task so that multifidelity approaches and reusing of existing data and software are possible.
The optimisation workflow is depicted in Figure 1 where four blocks are shown, specifically a baseline configuration, an optimiser, a performance model for the assessment of new wing planforms, and the final set of optimal configurations. The optimiser is based on the particle swarm optimisation (PSO) algorithm implemented in parallel [28]. It performs a search for global optimality with a gradient-free method. Its objective function is a weighted sum of single objectives. Specifically, any combination of drag coefficient , operating empty weight (OEW), and aerodynamic efficiency can be chosen as objective function. Their values are normalised with respect to the performance of the baseline, as usually done in industrial practice at early-stage design, before computing the weighted sum. Constraints are geometric as well as related to aircraft stability. The optimiser interrogates the performance model to compute the objective function for new configurations. Optimal solutions are returned when the maximum number of iterations is reached or a convergence criterion is met.
In Figure 1, a close-up on the performance model is provided which shows inputs and outputs placed on the left and right sides of the block, respectively. Four models concerning geometry, aerodynamics, mass estimation, and centre of gravity are employed by the performance model. This represents a typical pattern for the model-based framework since each model relies on one or more models to compute its output. The computation starts from parameters given as input by the optimiser, and a new geometry is produced by the parametric computer-aided design (CAD) model. It is an in-house code developed at Airbus Operations Ltd. for the purpose of early-stage design. Based on the CAD model, the mass estimation model provides a value for the aircraft OEW and its spacial distribution by exploiting semi-empirical methods discussed in [29]. Specifically, contributions to the total aircraft mass from wing, tail, fuselage, and engines are taken into account. Equations linking directly the components’ masses to their geometric properties, such as volume, height, and width, are derived from statistical analysis of existing aircraft [29], Ch. 8. They are employed to compute the masses for the new configurations produced by the optimiser. As an example, requirements in terms of structural stresses for the wing are taken into account by assuming a conventional wing structure, i.e., made of ribs, stringers, skin, etc., and computing the amount of material required to withstand the forces [29], Ch. 11. The centre of gravity (CG) position is obtained from the mass distribution by a model written on purpose since such information is needed to compute stability derivatives accurately. Consecutively, the aerodynamic model uses the CAD model, the estimated OEW value, and the CG position to trim the aircraft for the specific flight condition. The design lift coefficient is computed by balancing the lift force and maximum take-off weight. At the end, the aerodynamic model provides the lift, drag, and moment coefficient, angle of attack, horizontal tail rotation angle, lift distribution, and longitudinal stability derivative. A set of these quantities describes the aerodynamic performance for a given configuration.
Regarding the aerodynamic model, its reference implementation is based on the vortex lattice method (VLM) [30]. Specifically, the open-source VLM code AVL [31] is exploited as a tool to compute the lift distribution and moment coefficient, evaluate the stability derivatives, and predict the induced drag coefficient. The angle of attack and rotation of horizontal tail are both computed with a trimming procedure which is in charge of obtaining the desired lift coefficient and a null pitching moment. When the trimming procedure returns values of angle of attack and horizontal tail rotation outside the ranges and , respectively, results from AVL are considered unreliable due to intrinsic limitations of VLM at a large angle of attack [30] and the corresponding wing planform is discarded. Empirical methods from ESDU [32, 33] are employed to enhance the estimation of the drag coefficient by taking into account transonic flow and viscous effects. This aerodynamic model constitutes the reference for the surrogate models which are described next.
When it comes to the implementation, models communicate using a Python infrastructure; specifically, all data is stored in the memory belonging to the Performance model. This is the preferred way to store data since it minimises the time spent in transferring information. In fact, when a new model is executed, data available in the Performance model can be reused for further calculations. In the specific case of external tools such as AVL, inputs and outputs are performed using files. A Python wrapper is in charge of writing the input file, executing the external tool, and loading the results from files. Data is then stored in the Performance model. Parallel execution is performed using the multiprocessing module, and results are stored in the Performance model by exploiting interprocess communication.
2.1. Surrogate Modelling for Aerodynamics
Three surrogate modelling techniques are investigated in this paper, specifically RBF interpolation, an artificial neural network, and a hybrid model. A normalisation procedure is applied to input and output data of all three surrogate models, and it is described next. Denote as a vector containing the same scalar quantity for multiple configurations. It is then expressed as where is its average value and is the vector of fluctuations. Minimum and maximum values of the fluctuation are identified and denoted and , respectively. The normalised vector is defined as
Its components are in the range . Reconstruction of the original dataset is achieved by isolating in equation (1) and adding the average value to the result.
2.1.1. Radial Basis Function Interpolation
The RBF model implements interpolation of aerodynamic data using radial basis function (RBF) [34, 35]. The process is based on three steps. First, a database of output aerodynamic data corresponding to input wing planform configurations is populated with either experimental or precomputed data. Second, the RBF coefficients are obtained by exploiting the input/output relations in the database. Third, the actual interpolation is performed by computing sums of RBFs. Although the process is quite common [34], two main parameters affect directly the quality of results and they depend on the problem investigated, specifically number of samples and formulation of the radial basis function (also known as RBF kernel). Regarding the latter, three types are considered in this paper and reported here:
They all assign different importance to configurations in the parameter space according to their distance from the interpolation point. This is done on purpose to assess the quality of RBF interpolation when only nearby configurations are employed. In fact, if considering nearby configurations is sufficient for an accurate interpolation, a reduction in computational cost could be achieved by focusing on a smaller area of the parameter space. The thin spline plate (TPS) in equation (2) takes into account all configurations in the database, and it assigns greatest importance to the ones far away. Conversely, Gaussian function in equation (3) focuses on configurations closer to the interpolation point. Its value is 1 when , and it halves at so that almost no importance is given to far configuration. Modified Gaussian takes into account just few configurations around the interpolation point since its value halves exactly at .
2.1.2. Artificial Neural Network
An artificial neural network [36] based on multilayer perception algorithm [37] and a logistic activation function is employed as a surrogate for the aerodynamic model. Weights for the neurons are computed with backpropagation [38]. One hidden layer is adopted, and its number of neurons is the object of investigation. Guidelines for its estimation are given in [39] where a value between the number of inputs and outputs is suggested. The implementation relies on an open-source machine learning framework implemented in Python [40].
2.1.3. Hybrid Model
A third approach to modelling for aerodynamics is proposed as a trade-off between the surrogate and VLM. The method, called hybrid model, switches automatically between the underlying VLM model and its RBF surrogate when aerodynamic data is requested. This avoids the precomputation of data which is needed by both RBF and neural network. A first approximation of aerodynamic performance is obtained with the surrogate model based on RBF. Quality of results is assessed using a criterion, which is discussed in the next paragraph, and two scenarios are then possible. If the criterion is satisfied, aerodynamic data is returned and no further action is taken. If not, a VLM simulation is performed and its results are returned as well as stored in a database to improve the surrogate model. Note that only VLM solutions populate the database for RBF interpolation. When a number of new configurations are assessed with VLM, RBF coefficients are recomputed by inverting the RBF matrix and stored for successive interpolations. In practice, this happens at every iteration of the optimisation loop. Until the first iteration is complete, only VLM aerodynamics is adopted since a database of aerodynamic data is only available starting from the second iteration.
Regarding the criterion for switching between VLM and RBF, it exploits information about configurations previously assessed using VLM. The process is illustrated in Figure 2(a) with a simplified representation of the parameter space. Assume is the new configuration to be assessed. At first, its aerodynamic performance is computed with RBF interpolation and stored in a vector . A number of neighbours is then identified based on their Euclidean distances. Their aerodynamic performance are stored in vectors . The standard deviation is then computed with the differences with .
(a) Parameter space
(b) Cumulative function
When the value of where is a threshold, the estimation from RBF interpolation is accepted and returned. Otherwise, a VLM computation is performed and is updated. Although other criteria can be adopted, the standard deviation is proposed herein as a way to quantify the variation of aerodynamic performance in a set of neighbouring configurations. This approach mimics a rule of thumb which usually relies on the engineer’s experience: if the new data is reasonably similar to the existing one, it can be assumed reliable.
Regarding the threshold value , the arbitrariness is avoided by choosing its value based on statistical analysis performed a priori. Specifically, when a number of configurations assessed with VLM are available, each one is considered in turn and the standard deviation is computed considering its closest neighbours. Thus, a cumulative plot similar to the one depicted in Figure 2(b) can be produced. It shows the percentage of total configuration whose standard deviation is smaller than the value reported on the horizontal axis. This provides an objective way to define the threshold value. For example, the value can be chosen to include 95% of configurations as depicted in Figure 2(b). Please note that the hybrid model presents itself to the Performance model as one, monolithic model, i.e., using the same programming interface of the other two surrogate models, according to the model-based principles given in Section 2.
2.1.4. Error Quantification
In the context of optimisation, the error resulting from using the surrogate model is quantified by focusing on the optimal configurations only. The objective function is computed with VLM for all optimal configurations. The difference between the value from the surrogate model and the one from VLM is scaled with respect to the latter. Thus, average error and standard deviation can be computed and they can be then expressed as a percentage.
3. Formulation of Wing Planform Optimisation Problem
A wing planform is unequivocally defined by a sequence of seven parameters, which is herein called configuration, as shown in Figure 3(a). They represent displacements from the baseline wing planform to be optimised. Thus, when all the parameters are zero, the baseline wing is obtained. Precisely, parameters P1, P2, and P3 define the streamwise displacement of the trailing edge points with P1 defining the root chord. Note that the first, straight, part of the wing in Figure 3(a) is assumed to be inside the fuselage and its spanwise length is assumed to be constant. Parameters P4 and P5 define the leading-edge location so that chord lengths at the kink and at the tip are given by P5 - P2 and P4 - P3, respectively. Apart from the root chord P1 for which a maximum variation of from the reference is allowed, the parameters are not bounded. The set of geometrical constraints are reported in Figure 3(b). They mainly avoid planforms with negative sweep angles or M-shaped wings by requiring a positive clockwise angle between the three segments which define each edge. In addition, a minimum length constraint is imposed on the tip chord to avoid pointy wings. Besides geometrical constraints, a requirement for a stable configuration is formulated. The corresponding inequality is based on the value of the stability derivative which must be negative for longitudinally stable aircraft. Overall, the parametrisation allows the planform to assume a variety of shapes ranging from rectangular straight to tapered swept wings. Regarding the objective function, it can be chosen among OEW, aerodynamic efficiency, and drag coefficient. Its value is normalised with respect to the performance of the baseline geometry according to the practice usually adopted in early-stage aircraft design. The mathematical formulation of the problem is given in Table 1.
(a) Parameters
(b) Constraints
4. Results
The model-based software architecture was exploited to optimise the wing planform of the common research model (CRM) [41] flying at Mach number , altitude of 10000 m, take-off weight of , and maximum thrust of . The baseline geometry is illustrated in Figure 4 using the parametric CAD. It is composed of wing, fuselage, and horizontal tail, and the latter can be rotated around its mid-chord axis. The geometry of the torsion box, size, and distribution of ribs are defined parametrically depending on wing geometry.
Optimisations were first performed using VLM with empirical corrections for transonic flows as a tool for the aerodynamic model. They represent the reference solutions for the subsequent application of surrogate modelling. Optimisations were performed using PSO with a swarm size of 128 for 48 iterations. Successful convergence of the optimisation process is assumed when the maximum distance between positions of swarm particles and the swarm best particle is smaller than at two successive iterations. Velocity is updated for each particle every iteration by taking into account 50% of its velocity at the previous iteration, 50% of particle’s change of position, and 50% of swarm best particle’s velocity.
Results for three two-objective optimisations are discussed. The first one minimises OEW and maximises aerodynamic efficiency, and the Pareto front is depicted in Figure 5(a). Configuration A is chosen as a representative optimal solution. The wing planform for A is shown in Figure 5(b) where it is compared to the baseline. OEW is minimised by 3.7% with a reduction of wingspan. The kink moves outward so that the trailing edge sweep angle of the inner wing is reduced. The angle of attack is 6.22 deg, and it can be compared with the reference value of 4.15 deg which is obtained by trimming the baseline geometry. Overall, the wing is more slender and this factor improves aerodynamic efficiency by 5.2%. The quarter-chord sweep angle is larger when compared to the baseline. However, the drag coefficient, which is not included in the optimisation, increases by 30%.
(a) Pareto front for minimum OEW and maximum aerodynamic efficiency
(b) Wing planform of config. A
(c) Pareto front for minimum OEW and minimum drag coefficient
(d) Wing planform of config. B
(e) Pareto front for maximum aerodynamic efficiency and minimum drag coefficient
(f) Wing planform of config. C
Regarding the second optimisation, it minimises both OEW and drag coefficient. Configuration B is located on the resulting Pareto front as depicted in Figure 5(c). Its wing planform is very similar to the CRM baseline, Figure 5(d). The kink is slightly moved outward while few changes were found for the sweep angle, and the angle of attack is 4.10 deg. In particular the trailing edge sweep angle does not change significantly. Results for configuration B suggest that little modifications to the baseline improve both drag coefficient and OEW by 2.5% and 1%, respectively, while aerodynamic efficiency decreases by only 0.25%. Note that the baseline configuration is very close to the Pareto front and little margin for improvement was available to the optimiser.
A third optimisation was performed which maximises aerodynamic efficiency and minimises drag coefficient, and its results are reported in Figure 5(e). Configuration C is an optimal solution on the Pareto front. The wing planform depicted in Figure 5(f) is very slender. The wingspan is increased by more than 5 meters while the tip chord is kept at the minimum. The root chord is slightly reduced, and the kink is moved outward. The angle of attack is 4.03 deg and sweep angle for the inner wing decreases whereas it is unaltered for the outer part. Overall, aerodynamic efficiency increases by 4.5% and drag coefficient is reduced by 7.3%. However, OEW is not taken into account and it increases by 2%.
Initial values of location and velocity for swam particles in the PSO algorithm were generated randomly. Thus, each of the three optimisations described in Figure 5 was repeated 20 times. The Pareto front approximations were compared using two methods, specifically the -indicator [42] and the hypervolume () indicator [43]. As an example, results are shown for the optimisations which minimise both OEW and drag coefficient in Figure 6. The indicator is computed referring to the point at and for all 20 approximations. An average value of 0.34105 with a standard deviation 0.0011884 was found as depicted in Figure 6(a). The indicator was computed as well and its values for the 20 approximations are shown in Figure 6(b). The average value is 0.0049434, and the standard deviation is 0.0010499. Both indicators suggest that Pareto front approximations are all very close and the random initialisation of the PSO algorithm does not affect the quality of the final results. Statistical analyses were performed for all optimisations, but for the sake of brevity, graphs are not included here. When OEW is minimised and aerodynamic efficiency is maximised, Figure 5(a), the indicator has an average value of 0.025126 with a standard deviation of 0.00069854 with respect to the point at and . The indicator has an average value of 0.0034692 with a standard deviation of 0.00071947 for the same case. Regarding the optimisation which maximises aerodynamic efficiency and minimises drag coefficient (Figure 5(e)), the indicator has an average value of 0.090507 with a standard deviation of 0.0021697 with respect to the point at and . The indicator has an average value of 0.0058330 with a standard deviation of 0.0013592.
(a) Hypervolume indicator
(b) indicator
Apart from the two-objective optimisations, the framework can perform multiobjective optimisations as well. In Figure 7, results are reported for a three-objective optimisation aiming at minimising both OEW and drag coefficient as well as maximising aerodynamic efficiency simultaneously. The PSO optimiser runs with a swarm size of 128 for 48 iterations. The resulting three-dimensional Pareto front is shown in Figure 7(a). Identifying the reference configuration for CRM in the plot is quite difficult. It is dominated and, thus, located in the cloud of points. A configuration named D is chosen on the Pareto front and depicted in Figure 7(b). Regarding the inner wing, its sweep angle is smaller than the reference’s one. Conversely, the outer wing shows a slightly increased sweep angle. Overall, the wing span increases by 10%, the wing kink moves outward, and the tip chord is slightly smaller as well. Aerodynamic efficiency is increased by 5.6%, the drag coefficient is reduced by 4.1%, and OEW decreased by 0.8%.
(a) Single-objective optimisation, max E
(b) Three-objective
The three-objective optimisation performed with the model-based framework is a useful tool to provide the designer with a set of optimal configuration which improves multidisciplinary objectives simultaneously. However, visually comparing Pareto fronts resulting from different models is difficult since three-dimensional plots are involved. A more immediate comparison is available for the Pareto front involving two objectives only. Without loss of generality, the two-objective optimisation aiming at minimising both OEW and drag coefficient is chosen as benchmark simulation for the surrogate models which are described next. Specifically, the Pareto front in Figure 5(c) is assumed to be the reference solution to be reproduced. Although benchmark results are available for other two-objective optimisations, for example aiming at maximising aerodynamic efficiency and minimising either OEW or drag coefficient, they are not reported in the paper for the sake of brevity. Please note that the test case represents an academic case which is suitable to assessing the proposed methodology. For example, the structural model, which is described in Section 2, is limited to traditional wing structures (i.e., composed of ribs, panels, stringers, etc.) for which it can estimate the OEW based on considerations about structural integrity. Reliability of the structural model for any industrial configuration is not guaranteed. Inclusion of requirements which are usually adopted for aircraft design in the industrial context is not attempted in this paper, and as a consequence, optimal configurations might not represent the best configurations from an industrial point of view.
4.1. Surrogate Modelling for Aerodynamics
The capability of the model-based framework of replacing models at any time was exploited to replace the aerodynamic model based on VLM with surrogate models. Results from two types of surrogate models are presented in this section, specifically RBF interpolation and neural network. They both rely on a database of precomputed aerodynamic data to evaluate performance of new configurations. Inputs consist in the seven geometrical parameters which define the wing planform. Outputs are 8, specifically lift, drag, and moment coefficient, angle of attack, horizontal tail rotation, lift distribution, longitudinal stability derivative, and aerodynamic efficiency.
Regarding the aerodynamic database, it can be populated with preexisting experimental data when it is available. For this paper, it will be generated using VLM instead. The 7 parameters which define the wing planform (Table 1) are now limited to the range . An exception is made for the root chord which spans a smaller range . Aerodynamic data for a number of configurations uniformly distributed in the parameter space was computed with VLM. It is normalised in a range of so that all inputs and outputs have the same order of magnitude as described in Section 2.1. In addition, only interpolation is allowed in order to avoid inaccuracy due to extrapolation.
The surrogate model based on RBF is described first. Its application is composed of three steps. First, an RBF function is chosen. Secondly, the RBF matrix is assembled, inverted, and stored in memory. Thirdly, the interpolation is computed by exploiting the inverted matrix only when aerodynamic data is requested by the performance model. Two choices have to be made when using RBF interpolation, specifically the type of function and the number of samples. Concerning the radial basis function, it should be chosen according to the data to be interpolated. Three radial functions are compared in Figure 8(a) for values of distance . Note that the maximum distance between any two configurations is because of the input data normalisation. The three functions assign different weights to configurations in the parameter space as already explained in Section 2.1. In the specific case, Gaussian function covers the half parameter space and its modified version spans a quarter of parameter space. Performances of the three functions were compared by running a two-objective optimisation aiming at minimising both OEW and drag coefficient. Results are depicted in Figure 8(b) where a comparison with VLM is included as well. Three approximations of the Pareto front were produced. Using the Gaussian function leads to results which reproduce the general trend with an average error of 10.0%. Results produced with the modified Gaussian function match accurately the reference for configurations located in the central region of the Pareto front, and an average error of 11.5% is found. Conversely, TPS, which focuses on far configurations instead, provides the best results in any region of the Pareto front with an average error of 7.62%.
(a) RBF function
(b) Pareto front
Regarding the number of samples, the benchmark optimisation was performed using RBF interpolation with databases containing 5, 10, 50, 100, 200, and 300 samples uniformly distributed in the parameter space. The average error was computed as described in Section 2.1 for each Pareto fronts corresponding to a different sample number. The resulting convergence curve is presented in Figure 9(a). Using databases with less than 400 configurations returns an error larger than 9% with respect to the VLM reference. The average error is smaller than 8% (precisely 7.62%) when 438 samples are employed. Although a large number of configurations is impractical for industrial application, it could be useful for academic investigation. Hence, an additional simulation with a database containing 15220 configurations was performed and the average error becomes 3.70%. Note that the TPS function consistently provided the least error regardless of the number of samples. For example, when 15220 configurations are employed, an average error of 3.70% is found with TPS whereas Gaussian and modified Gaussian functions provide an average error of 4.0% and 9.34%, respectively. In Figure 9(b), Pareto fronts identified using TPS for three databases are reported. Increasing the number of samples by a factor of 43 lowers the average error from 21.43% (obtained with 10 samples) to 7.62% (computed for 438 samples). The optimisation based on TPS and 438 samples was repeated 20 times in order to perform a statistical analysis on the resulting Pareto front approximations. They were compared with the indicator, whose average value was evaluated in 0.0091440 with a standard deviation of 0.0044769, and the indicator, which showed an average value of 0.36016 with a standard deviation of 0.0022778 with respect to the point at and . The reference point was chosen for consistency with the reference results reported in Section 4 and summarised in Figure 6(a). Similar results were obtained for the other optimisations involving a different number of samples as well as RBF functions. It was found that the variability related to the stochastic initialisation of the PSO algorithm is smaller than the error between results produced with VLM and RBF for any given combination of function and number of samples.
(a) Average error with respect to the number of samples using TPS
(b) Pareto front produced with TPS using 438 and 15220 samples
The surrogate model based on the neural network is described next. Note that the same database with 438 samples and the same scaling of input and output variables employed for the RBF model were used to train the neural network too. A brief investigation was performed to choose the size of the network. Regarding inner and outer layers, the amount of neurons is defined by the number of inputs and outputs. According to the number of inputs (7) and outputs (8), guidelines proposed in [39] suggest a number of hidden neurons ranging from 7 to 20. Thus, three two-objective simulations aiming at minimising OEW and drag coefficient were performed with a size of the hidden layer ranging between 4 and 16. In addition, three large values (specifically 64, 256, and 512) were included in the investigation as well. Results are reported in Figure 10(a) and they are compared to VLM reference. Using 4 neurons produces configurations which dominate the VLM ones. However, a quantification of the error showed that the neural network overestimated the performance and an average error of 10.4% was found. Increasing the number of neurons to 8 improves the results. Using neural networks, the smallest average error of 9.6% is found for 8 neurons. Further increasing their number to 16, 64, 256, and 512 leads to inaccuracies, and the average error becomes 11.2%, 12.7%, 11.2%, and 14.6%, respectively. In Figure 10(b), a comparison of Pareto fronts produced with neural network and RBF using the same database of 438 samples is presented. Both surrogate models are able to identify the Pareto front accurately. However, a better agreement is found near the central region of the Pareto front when it comes to neural network.
(a) Neural network performance
(b) Pareto front comparison
Results from a statistical analysis are provided here for the surrogate based on neural network and trained with 438 samples. It was performed on 20 Pareto front approximations which were produced by running the optimisation with 8 neurons. An average value of 0.0047656 and a standard deviation of 0.0010489 were found for the indicator. The indicator has an average value of 0.36429 and a standard deviation of 0.0033477 with respect to the point at and . This is the same reference point which was chosen for analysing the reference results in Section 4 as well. Similar results were obtained for the optimisation involving a different number of neurons and the same training set. They confirm that the error introduced by approximating VLM with a neural network is larger than the variability associated to the stochastic initialisation of the PSO algorithm.
A summary concerning the computational cost of surrogate modelling is given. Statistics is reported here for the RBF model which uses TSP and the neural network with 8 neurons. Note that only the aerodynamic model is replaced, and when evaluating the objective function, part of computational time is still employed to run the mass estimation model and stability model as well as to transfer data. Computations were performed using a single core of an Intel i7 4810MQ CPU. The total computational cost is split into two contributions. The first one is the computation of the aerodynamic database which was used by both surrogate models. It took 438 evaluations of the objective function using VLM for a total of 17 m 50 s which includes 1 s for each evaluation to execute external models and transfer the data. Regarding RBF, the matrix inversion needed to compute RBF coefficients took roughly 5 s and it was performed only once. Secondly, interpolation of aerodynamic data for one configuration took 0.23 s and this task had to be repeated for 48 iterations and a swarm size of 128. Thus, the total time needed to run the MDO simulation with RBF model is 41 m 20 s split between building the surrogate model (17 m 50 s) and using it (23 m 30 s). Regarding the surrogate model based on neural network, it needs 25 m 51 s to perform the optimisation for a total of 43 m 41 s, which is comparable to the amount of time needed by the RBF model. These numbers are compared to the VLM reference. The cost of evaluating the objective function using aerodynamics based on VLM is 1.5 s on an Intel i7 4810MQ CPU, and this task was repeated for 48 iterations and a swarm size of 128. The total time employed to perform the optimisation is about 2 h 33 m using a single core.
4.2. Hybrid Modelling for Aerodynamics
The computation of the threshold for the hybrid model was performed by applying the procedure described in Section 2.1 to the reference results. Specifically, two sets of results were considered. The first one is composed of all configurations assessed with VLM during the two-objective optimisation aiming to minimise OEW and drag coefficient. They were depicted in Figure 5(c). The second one contains data produced with RBF model based on the TPS function when performing the same optimisation. Results were already shown in Figure 8. Each set is considered in turn, and the cumulative distribution is given in Figure 11(a). Such distribution is shown for both databases and for three numbers of neighbours, specifically 2, 4, and 8. Values of ranging from 0 to 0.35 are found. Note that aerodynamic data for 95% of configurations has a standard deviation of less than when 2 or 4 neighbours are considered. A close-up in the region around is depicted in Figure 11(b). It is shown that results produced with VLM and RBF converge to when the number of neighbouring configurations is 8. Information in Figure 11 was exploited to set a threshold and a number of neighbours as a criterion to switch between VLM and RBF model. It means that data from the RBF model is considered reliable when its standard deviation from the closest 5 neighbours is .
(a) Global distribution
(b) Close-up
Performing a two-objective optimisation to minimise OEW and drag coefficient with the hybrid model using a swarm size of 128 for 48 iterations, a threshold and a number of neighbours took 44 m 17 s on an Intel i7 4810MQ CPU using a single core. Results are shown in Figure 12. The ratio between the number of evaluations of the objective function using VLM and RBF is depicted in Figure 12(a). The first iteration is performed with VLM since no database is available to perform RBF interpolation. The following 3 iterations show no contribution by the surrogate model as well. This is due to values of for interpolated aerodynamic data which are above the threshold . Starting from the fifth iteration, results from the surrogate model are considered accurate with . The number of evaluations performed with RBF increases, and by the end of the simulation, a total of 43% of configurations was assessed exclusively with the surrogate model. The Pareto front obtained with the hybrid model is compared to the ones computed with RBF and VLM in Figure 12(b). Overall, a very good agreement is found between the hybrid model and the VLM reference. Results for the lower region of the Pareto front match accurately and the upper region, for which fewer configurations are available, are identified too. In addition, a comparison is provided between the hybrid model and the RBF one. The latter is able to reproduce the upper region of the Pareto front properly. Concerning the lower part, the hybrid model provides more accurate results. This result was confirmed by performing the hybrid optimisation 20 times and comparing the resulting Pareto fronts. The average value of the indicator is 0.0061991 with a standard deviation of 0.0026108 whereas the indicator has an average value of 0.36016 and a standard deviation of 0.002277 with respect to the reference point at and . This is chosen for consistency with the reference results reported in Section 4.
(a) Type of simulation
(b) Pareto front
A quantitative summary of the investigation concerning surrogate modelling for aerodynamics using RBF, neural network, and hybrid model is provided in Table 2. Results for the same two-objective optimisation aiming at minimising the OEW and drag coefficient are compared. The surrogate model based on RBF interpolation identifies the Pareto front with an average error of 7.62%. The computational cost is reduced to almost a quarter of the VLM, and this figure is similar for all surrogate models. The neural network is the least performing of the surrogate models since its error is 9.60% on average and the corresponding computational time is decreased to 28.5%. Conversely, the hybrid model provides the best results with an error smaller than 1% at the same computational cost. The total time needed by the hybrid optimisation is 44 m 17 s since evaluations took 1.5 s and 0.23 s using VLM and RBF, respectively. A higher number of VLM evaluations (around 1500) is performed for the hybrid model when compared to the RBF surrogate, which used 438 VLM simulations for the sampling instead. However, the computational time for data-transferring and execution of external models is negligible for the hybrid optimisation since both VLM and RBF are now part of the same model and no additional overheads are needed. This is the key element of the cost saving which is reported here. Note also that the computational time and accuracy of the hybrid model change as a function of the threshold. For example, further reduction of time was obtained by performing the same two-objective optimisation using a threshold . An average error of 1.18% was found, and the total computational time was 27 m 41 s which is 18% of the reference.
5. Conclusions
The paper describes an investigation into surrogate modelling for wing planform multidisciplinary optimisation which exploits a novel model-based framework. The workflow is composed of a particle swarm optimiser and a performance model which is in charge of computing the objective function for each configuration to be assessed. It employs a mass estimation model to estimate the mass, a stability model for the computation of longitudinal stability derivative, and an aerodynamic model.
Four models were exploited for the calculation of aerodynamic performance. The original model relies on the vortex lattice method, and it provided the reference results. Three models provide surrogates for the calculation of aerodynamic data. The first one is based on radial basis function interpolation whereas the second one employs an artificial neural network. They both require the precomputation of an aerodynamic database prior of their usage. A third, hybrid approach based on both VLM and RBF is proposed to avoid the precomputation. Switching between the two underlying models is based on a novel criterion which is defined in an objective way. When assessing a new configuration, results provided by the surrogate model are accepted if their standard deviation with respect to neighbouring configurations in the parameter space is below a given threshold. In turn, the threshold value is obtained by analysing results of previous optimisations only once. Results and performance of surrogate models were compared by performing the same benchmark optimisation. It minimises two objectives, specifically operating empty weight and drag coefficient.
Key findings are summarised as follows. Regarding the RBF model, employing the thin plate spline function provides the most accurate results. This is valid for all numbers of samples which were investigated. The convergence curve is almost flat since an average error of 10% is obtained with 50 samples, but 438 are needed to go below 8%. The computational cost is reduced to a third of the original one. Regarding the artificial neural network, the investigation showed that a number of 8 neurons provide results which are comparable to the ones from the RBF model when using the same precomputed database. The computational cost is similar too, and increasing the number of neurons does not improve the final results in terms of average error. Regarding the hybrid model, the results show an excellent agreement with the original model. The threshold affects both average error and computational cost. Modifying its value, it is possible to act directly on the percentage of evaluations performed with the original model or its surrogate. Thus, the level of accuracy, and the corresponding computational cost, can be arbitrarily decided.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors want to acknowledge the support received from Airbus Operations Ltd., in particular from the Future Projects Office. The research is part of the Agile Wing Integration project sponsored by Innovate UK.