A Study on the Surrogate-Based Optimization of Flexible Wings Considering a Flutter Constraint †

: Accounting for aeroelastic phenomena, such as flutter, in the conceptual design phase is becoming more important as the trend toward increasing the wing aspect ratio forges ahead. However, this task is computationally expensive, especially when utilizing high-fidelity simulations and numerical optimization. Thus, the development of efficient computational strategies is necessary. With this goal in mind, this work proposes a surrogate-based optimization (SBO) methodology for wing design using a predefined machine learning model. For this purpose, a custom-made Python framework was built based on different open-source codes. The test subject was the classical Goland wing, parameterized to allow for SBO. The process consists of employing a Latin Hypercube Sampling plan and subsequently simulating the resulting wing on SHARPy to generate a dataset. A regression-based machine learning model is then used to build surrogate models for lift and drag coefficients, structural mass, and flutter speed. Finally, after validating the surrogate model, a multi-objective optimization problem aiming to maximize the lift-to-drag ratio and minimize the structural mass is solved through NSGA-II, considering a flutter constraint. This SBO methodology was successfully tested, reaching reductions of three orders of magnitude in the optimization computational time.


Introduction
The need to account for aeroelastic phenomena, namely, flutter, during the conceptual design phase has gained relevancy as wing designs become slenderer to reduce the induced drag [1].However, this need comes with a computational challenge, especially when integrated into a Multidisciplinary Design Optimization (MDO) [2] environment to perform Fluid-Structure Interaction (FSI) simulations to estimate flutter speed.Therefore, searching for more computationally efficient solutions that are able to maximize wing performance without reaching flutter within the flight envelope is required.With this goal in mind, different approaches are considered in the literature, namely, multi-fidelity models [3,4], conventional [5] and machine learning-based surrogate models [6], and reduced-order models [7].
Even though reduced-order and surrogate models have been applied to aeroelasticity problems for some time [13], especially to study unsteady nonlinear aeroelasticity problems that arise from unconventional features such as limit-cycle oscillations [14][15][16], their usage for MDO problems considering flutter is scarce [17][18][19] in the open literature, particularly for MDO problems considering aerodynamic performance and structural weight.For instance, Sohst et al. [17] developed an MDO strategy that uses multi-fidelity solvers and surrogate models to design strut-braced and high-aspect-ratio wings considering flutter and stress constraints.They found that, although the flutter constraint was respected in the optimization process, unaccounted buckling was observed in a nonlinear assessment of the optimized strut-braced wing aircraft design.Cea and Palacios [18] implemented an optimization framework based on ROMs that couples different open-source codes, including SHARPy, for minimizing the mass of a flexible strut-braced wing while accounting for flutter speed as a constraint.Toffol and Ricci [19] developed a methodology combining in-house codes with surrogate models to optimize the structural layout of a conventional aircraft such that its mass is minimized while considering stress and flutter constraints.The application of surrogate models based on machine learning is scarcer.Nevertheless, this research field is starting to emerge, particularly for aerodynamics [6,20].
In this work, a surrogate-based optimization (SBO) strategy is presented to conceptually design flexible wings based on machine learning-based surrogate models with the aim of maximizing aerodynamic performance, minimizing structural mass, and satisfying the flutter constraint.

Methodology
The methodology employed in this work and depicted in Figure 1 adheres to three fundamental steps commonly found in SBO: (i) dataset definition and generation; (ii) surrogate model training, validation, and testing; and (iii) the optimization process.All these tasks were performed using available open-source Python codes and specially developed Python scripts.For illustrative purposes, the methodology is applied for a cruise speed of 30 m/s, but it is applicable to other speeds.Firstly, the Goland wing [21] was selected to be the baseline model for this work.The reason for this choice was two-fold: (i) experimental data exist to validate the implementation; (ii) it has been already used to validate SHARPy [22], the open-source code chosen for the aeroelastic simulations and developed by Imperial College London.SHARPy provides a coupled aeroelastic solution considering the Unsteady Vortex Lattice Method (UVLM) [23] and Geometrically Exact Beam Model (GEBM) [24] for aerodynamic and structural analyses, respectively.Krylov and modal projections are used in SHARPy to linearize aerodynamic and structural subsystems, respectively.These are then used to estimate the flutter speed by means of the iterative p-k method at the nonlinear aeroelastic equilibrium.The flutter speed was predicted in the SHARPy implementation used for this work to be under 5% of the relative error for the Goland wing, in agreement with the value obtained by SHARPy developers [25].To improve the accuracy of the drag estimation, the flat-plate theory [26] was implemented, offering an estimation of skin-friction drag multiplied by form (more adequate for small angles of attack) and interference factors.

Good
Currently, there are some efforts being made by the SHARPy team [27,28] to include viscous considerations in the UVLM, which are typically not of major importance when the focus is on aeroelastic phenomena.

Design Space
In regard to the design variables, the choice rested on four key parameters in aeroelasticity: wing aspect ratio (AR), sweep angle (Λ), torsional stiffness (GJ), and angle of attack (α).The AR is vital to improve aerodynamic efficiency and fuel consumption, although it poses aeroelastic and structural challenges [1].The sweep angle (Λ) influences transonic flow onset and aircraft aeroelastic stability [29].GJ relates how wings respond to aerodynamic loads and is critical to averting aeroelastic instabilities [29].Lastly, α impacts lift and drag, affecting flight stability.Besides their relevancy to and impact on the aeroelastic behavior, SHARPy's ability to modify these variables without coding was fundamental in their selection.The upper and lower boundaries (UB and LB, respectively) of these parameters are listed in Table 1 alongside the initial values that correspond to the Goland wing specifications.These boundaries were later adjusted based on the preliminary optimization results.For sampling, Latin Hypercube Sampling (LHS) was adopted using the pyDOE2 package [30] due to its uniform coverage across the design space and its efficient sampling scheme.Simulations conducted in SHARPy focused on the wing metrics to be used in the optimization: structural mass, flutter speed, and aerodynamic coefficients (C L and C D ).For this study, 2000 simulation samples were generated.Despite this number appearing to be modest given the complexity of the problem, the high computational demands of SHARPy made this a significant effort.Each simulation averaged around 93.5 s.The simulations were run on a host computer with Windows 11, an Intel i7-10510U processor, 16 GB DDR4 RAM, and a 477 GB SSD.This machine also used NVIDIA GeForce MX250 graphics.Within this system, a virtual machine operated on Ubuntu 22.10, utilizing 4 cores, 11 GB RAM, and 59.2 GB of storage space.Virtualization was facilitated by Oracle VM VirtualBox 7.0.6.

Surrogate Model Training, Validation, and Testing
This subsection provides an overview of surrogate model development, from training to testing.The goal is to find a machine learning-based surrogate model that meets computational needs and offers a high predictive accuracy.Our criteria for selecting the model include accuracy, efficiency, and robustness.
We considered various regression-based machine learning models from the scikitlearn Python library [31], including Bayesian Ridge, Decision Tree Regressor, Extra Trees Regressor, Lasso, ARD Regression, and Linear Regression.For the model generation, we considered a 60%-30%-10% split for training, validation, and testing, respectively, based on experimenting with different combinations for generating surrogates of C L , C D , mass, and flutter speed.A random procedure was followed to assign the data points to each set.These machine learning models were set with random hyperparameters to gauge their baseline performance.Each was then assessed on the validation subset.
The methods were evaluated for flutter speed, mass, drag, and lift coefficients, i.e., the output parameters from SHARPy that were used for the optimization problem definition.All the considered models were able to well represent the structural mass, drag, and lift coefficients, with values of R 2 higher than 0.99 and of Root Mean Squared Error (RMSE) lower than 1%.However, the models evidenced difficulty in provide accurate predictions of flutter speed, as shown in Table 2. While, for this model, only the Decision Tree method approached the performance of Extra Trees (still remaining inferior), for the other models, the difference between the methods was minimal.With these considerations in mind and for consistency and simplicity in the workflow, we chose to use the Extra Trees Regressor as the surrogate model method for all output variables.After determining the Extra Trees Regressor's aptitude, the training and validation processes began.
Alternatively, another SBO strategy could be explored in the future, where new infill points are added alongside the optimization using an appropriate acquisition function and Bayesian Optimization [32].This approach was recently applied to an aero-structural problem by Cardoso et al. [33].

Training and Validation
The model was trained using the dataset, focusing on hyperparameter tuning.Hyperparameters were carefully adjusted to ensure model accuracy and prevent overfitting.The optimal settings were identified via a Grid Search, using RMSE and R 2 for selection.The model's hyperparameters analyzed were the number of trees in the ensemble, maximum tree depth, minimum samples at a leaf node, and minimum samples for a node split [34].The objective during training is to minimize a loss function for improved prediction based on input data, ensuring that hyperparameters are set for both known (training set) and unseen (validation set) data.

Testing
After training and validation, the model is tested on a dataset it has not seen before, known as the test set, to gauge the model's effectiveness.Learning curves further evaluate the model's generalization.The final surrogate models used for testing and the later optimization process are based on the optimal hyperparameters obtained for each model after tuning.The suitability of surrogate models for the aeroelastic analyses required for the optimization process was visualized using parameters from the testing dataset.The results for each model are shown in Figure 2. From the results, the surrogate models accurately predict the mass and lift coefficient.Discrepancies in the drag coefficient suggest complexities in aeroelastic interactions in certain regions of the design space.The flutter speed showed the greatest difference, with an R 2 of 0.921.The flutter model might have been overfitting, as can be seen from the learning curves in Figure 3, and more data might improve its performance.Challenges in flutter modeling could arise from estimating flutter speed using the iterative p-k method.

Optimization Process
The last step in the implemented SBO methodology is to solve the optimization problem, which is a multi-objective one, as we want to maximize aerodynamic performance in terms of the lift-to-drag ratio (C L /C D ) and minimize structural mass while ensuring that the flutter speed (V flutter ) is far from the cruise speed (V cruise ).Mathematically, the optimization problem can be stated as follows: with respect to x = (AR, Λ, GJ, α) where f , x, and c denote the classical notation for objective functions, the design variable set, and constraints, respectively.The factor 1.5 provides a safety factor against flutter, a standard in aerospace engineering [35].Since this is a multi-objective problem and the run time of the surrogate model is inexpensive, a gradient-free optimizer was chosen to better explore the design space, namely, the implementation of Non-dominated Sorting Genetic Algorithm II (NSGA-II) [36] in pygmo [37].However, as NSGA-II is an unconstrained optimization algorithm, a penalty was added to the objective function when the flutter constraint is not respected.In this work, the population size, mutation rate, and crossover rate were adjusted to improve the resulting Pareto front while ensuring a good balance between the computational cost, hypervolume, and ranking.It is worth mentioning that hypervolume measures the space occupied by solutions in the target space, while ranking helps select solutions based on levels of non-dominance [36].
After obtaining the Pareto fronts with the initial design variable boundaries (Table 1), a significant discrepancy was noted in some results.A concentration of solutions was found near the lower limit for the sweep angle and torsional stiffness, as can be observed in Figure 4, suggesting that extending these limits might produce better solutions.In addition, surrogate models showed high prediction errors for high aspect ratios.

Wi n g ma s s
Wi n g ma s s Wi n g ma s s Wi n g ma s s To solve this problem, the boundaries of the design variables were extended and the dataset was augmented, focusing particularly on the critical aspect ratio interval.An iterative analysis and modification of the dataset show a systematic approach to addressing optimization challenges, emphasizing the importance of identifying weaknesses in the model and making necessary changes.

Results
The outcomes of the optimization procedures for three cruise speeds are examined in this section.Using surrogate models for design optimization, optimal solutions are derived and then compared to reference from SHARPy simulations.

Case 1: Cruise Speed of 30 m/s
The first cruise speed considered was 30 m/s.Initially, an optimization was performed with a specific set of parameters obtained by tuning, and the solutions were compared with SHARPy simulations.Even though the results were accurate, especially for C L , C D , and mass, they evidenced discrepancies in flutter velocity estimation as mentioned previously and presented in Table 2. To solve this problem, the boundaries of the design variables were expanded, as explained in the previous section, adding more data points to improve the accuracy of the prediction.A new optimization with the updated data showed a reduced error for most parameters.The final surrogate models were built based on a dataset of 2500 points with the metrics shown in Table 3.However, to further improve accuracy, more computational resources were allocated.Increasing the population size from 180 to 280 and the number of generations from 50 to 200 improved the flutter speed estimation, confining the metric deviations within a 5% margin of error, as can be observed in Table 4. Yet, this came at the cost of computational efficiency, with a five-fold increase in computational time.Pareto front analysis revealed trade-offs between the two optimization objectives.The extended configuration covered a wider Pareto range, suggesting its ability to explore a complete solution space and produce better results, particularly toward higher C L /C D values, as can be seen in Figure 5. Throughout the study, the flutter constraint remained inactive, prompting further investigation of the Goland wing's behavior at higher cruise speeds.

Case 2: Cruise Speed of 60 m/s
For Case 2, the optimization was conducted at a cruise speed of 60 m/s using the same methodology as the one applied for 30 m/s.The Extra Trees Regressor was identified as the best surrogate model, with the performance metrics presented in Table 5.For a cruise speed of 60 m/s, from 9000 combinations, the flutter constraint was active only 0.322% of the time.The design variables for which the flutter constraint was active fell within the following specific ranges: AR of 15.2-15.4;Λ of 10.4-11.3degrees; GJ between 1.62 × 10 6 and 1.70 × 10 6 N.m 2 M; and α of 6.86-6.95degrees.
Nevertheless, the optimization algorithm found 180 optimal solutions where the flutter constraint was never active.These solutions differed in design variables from those where the flutter phenomenon was a concern.A subsequent, more detailed analysis aimed to reduce the relative error to under 5% for all outcomes.
The results in Table 6 show that all percentage errors were below the threshold after a deeper exploration of the design space and an increase in the optimization parameters, as was implemented for the 30 m/s analysis.
However, this intensive optimization resulted in a 491% increase in computational time.Nevertheless, the potential for greater accuracy seems to justify the additional computational effort.Comparing the Pareto fronts in Figure 6, an improvement was evident, especially at higher C L /C D values.In these cases, the advanced optimization provided better aerodynamic efficiency with a lower mass increase.However, for some design choices, there was a significant overlap in results between the initial and advanced optimizations, indicating that increasing iterations and population sizes do not always lead to different results.

Case 3: Cruise Speed of 130 m/s
For the final case, considering a cruise speed of 130 m/s, the flutter phenomenon emerged for most of points in the dataset, which made optimization difficult.The performance of the surrogate models on which the optimization was based is shown in Table 7.The initial optimization process, despite using an expanded dataset, produced results with high percentage errors compared with SHARPy, especially in terms of C L /C D and flutter velocity.This was primarily due to the active flutter constraint in a vast majority (82.2%) of the examined design combinations that led to the introduction of a penalty in the objective function.Consequently, in the Pareto front of 180 optimal solutions, 24.4% had active flutter issues.To address this, the optimization was reevaluated by increasing the population size to 280 and the number of generations to 200, as was implemented for the previous cruise speeds.
The modified approach yielded a substantial reduction in the number of solutions with an active flutter constraint, down to 7.5%.However, this improvement came at a significant computational cost, with a 690% increase in time compared to the initial optimization.Despite the improvement, the percentage errors were still remarkably high, and the Pareto fronts showed several infeasible solutions due to the active flutter problem and the penalties associated with the results.
However, a closer examination of the results, focusing particularly on the solutions not affected by flutter, showed commendable accuracy.For the most exhaustive optimization (with 200 generations and a population of 280), errors were consistently less than 5%, as listed in Table 8.These results, in line with those obtained for the cruise speeds previously examined, clearly reveal that the high penalties were the main cause of the high percentage errors previously observed.
Two distinct groups of solutions can be immediately identified from the graph in Figure 7.In both cases, the optimized design space is significantly reduced, but the second offers a more diverse set of solutions.Moreover, considering only the feasible solutions, we obtained results that are aligned with the previous cases.This further strengthens the viability of the implemented SBO strategy and shows that, despite the obstacles represented by penalties, the algorithm is able to offer very good solutions.

Comparison and Discussion of the Results
A final comparison is presented in Table 9, which details the results obtained using consistent parameters and dataset sizes across the three analyzed cruise speeds.The following optimization parameters were set: a population size of 280; 200 generations; a mutation rate of 0.5; and a crossover rate of 0.8.For the last analysis, using a cruise speed of 130 m/s, the results considered are just the feasible ones, i.e., those without the flutter constraint active.Firstly, one can notice that, as expected, when improving one parameter, the other worsens accordingly, typical of the nature of multi-objective optimization.High values of C L /C D consequently correspond to higher mass values.
It is possible to observe that the aspect ratio is consistently higher when maximizing C L /C D compared to minimizing the mass across all cruise speeds.This observation aligns with established aerodynamic principles: wings with a higher AR have longer spans relative to their chord, which reduces the induced drag and consequently improves the lift-to-drag ratio.However, the flip side of this advantage is that longer, slenderer wings tend to be more flexible.Indeed, it is possible to observe how the flutter phenomenon in this case occurs at lower speeds.To counterbalance this issue, these wings are heavier and require a higher torsional stiffness to comply with the flutter speed constraint.Deviating from the usual trend is the relatively high aspect ratio observed at the highest cruise speed for the mass minimization objective.This unexpected outcome may be attributed to the narrower feasible design space identified by the surrogate model, and the optimizer that might have been trapped in a local minima region.
At the lower cruise speeds (30 m/s and 60 m/s), when the primary objective is to minimize mass, a positive and more pronounced sweep angle is observed.This might be associated with a structural benefit, redistributing the lift toward the root, which can result in a lighter wing structure.However, when maximizing C L /C D , a lower sweep may favor aerodynamic efficiency at these speeds.
Variations in torsional stiffness and the angle of attack across the different objectives and velocities are subtler compared to other parameters.However, it is worth noting that higher C L /C D values are associated with lower angles of attack.This trend aligns with the expectation that a reduced angle of attack would lead to decreased aerodynamic drag.Lowering the angle of attack while still achieving adequate lift is a strategy to enhance the aerodynamic efficiency of the wing.
Another important observation about the aspect ratio is that for both 30 m/s and 60 m/s cruise speeds, the optimal solutions approach the boundaries of the defined design space, which was set between 6 and 16.This is significant, as it indicates that the design space might not fully encompass the optimal regions for these objectives at these speeds.
The plot presented in Figure 8 delineates the Pareto fronts for the three analyzed cruise speeds, 30 m/s, 60 m/s, and 130 m/s, derived using the same optimization parameters.It is important to note that the results at 130 m/s represent only the feasible solutions.As the cruise speed increases, the Pareto front tends to be narrower and shift toward lower C L /C D values, as expected.The results at 130 m/s present an evident distinction.Unlike the relatively distributed Pareto fronts at 30 m/s and 60 m/s, the solutions at 130 m/s are remarkably more clustered.This localized concentration indicates that, at this higher cruise speed, design solutions that avoid flutter are limited, especially within the limits of our design space.Moreover, these solutions at 130 m/s are characterized by relatively high masses and moderate C L /C D values.At higher speeds, the overly stringent constraint of flutter speed within our design space leads to results in which, to avoid flutter, the wing structure needs to be heavier, consequently compromising its aerodynamic efficiency.

Concluding Remarks
An SBO methodology for designing wings considering both performance and aeroelasticity constraints was successfully implemented using four different open-source Python codes.This methodology uses Extra Tree Regressors to mitigate the computational cost associated with this wing design problem, and it was tested to solve a multi-objective optimization problem considering three cruise speeds.
Although the surrogate models exhibited suitable performance, certain discrepancies appear, especially in drag coefficient predictions and flutter speed, attributed to observed aeroelastic nonlinearities.Nevertheless, the adopted surrogate models, despite the challenges faced, yielded acceptable results, with errors generally lower than 10%.
To improve the accuracy of the surrogate models in the three cases analyzed, the original design space was revised for each of these cases individually.After this treatment, the optimized designs that maximize C L /C D and minimize mass while keeping flutter outside the prescribed boundary presented differences under 5% compared to the corresponding SHARPy simulations.Although this refinement process improved results, it also increased the computational time substantially.
Since the goal was to efficiently integrate aeroelastic analysis, the computational gains were also significant.Direct simulations with SHARPy take about 93.49 s each, while surrogate models offer predictions in under 0.01 s.In terms of optimization, the following observations are made:

•
Optimization with 50 generations and 180 population averages 321.76 s, totaling 9000 evaluations.Direct SHARPy simulations would take nearly 10 days.• A scenario with 200 generations and 280 population takes around 1825.03 s, equaling 56,000 evaluations.Using SHARPy would mean approximately 60 days.
These results highlight the substantial computational time savings when using surrogate models, particularly in optimization scenarios requiring multiple evaluations.

Figure 2 .
Figure 2. Visual comparison of the results obtained using the SHARPy simulations (in blue) and the surrogate model (in orange) in the testing phase of the SMs.

Figure 3 .
Figure 3. Learning curves of flutter velocity model.

Figure 4 .
Figure 4. Distribution of optimization solutions along design variable boundaries for both objectives: C L /C D maximization (above, in green) and mass minimization (below, in red).Units of mass are given in kg.

Table 1 .
The design variables and respective boundaries used for the dataset generation and optimization problem.

Table 3 .
Final RMSE and R 2 metrics for Case 1 (cruise speed of 30 m/s).

Table 4 .
Optimized solutions obtained with the surrogate models compared to the corresponding SHARPy results for Case 1 (cruise speed of 30 m/s).

Table 5 .
Final RMSE and R 2 metrics for Case 2 (cruise speed of 60 m/s).

Table 6 .
Optimized solutions obtained with the surrogate models compared to the corresponding SHARPy results for Case 2 (cruise speed of 60 m/s).
L D

Table 7 .
Final RMSE and R 2 metrics for Case 3 (cruise speed of 130 m/s).

Table 8 .
Optimized solutions obtained with the surrogate models compared to the corresponding SHARPy results for Case 3 (cruise speed of 130 m/s).

Table 9 .
Optimization results across the analyzed cruise speeds.