A Comparison of Modern Metaheuristics for Multi-Objective Optimization of Transonic Aeroelasticity in a Tow-Steered Composite Wing

Abstract

This study proposes a design procedure for the multi-objective aeroelastic optimization of a tow-steered composite wing structure that operates at transonic speed. The aerodynamic influence coefficient matrix is generated using the doublet lattice method, with the steady-state component further refined through high-fidelity computational fluid dynamics (CFD) analysis to enhance accuracy in transonic conditions. Finite element analysis (FEA) is used to perform structural analysis. A multi-objective transonic aeroelastic optimization problem is formulated for the tow-steered composite wing structure, where the objective functions are designed for mass and critical speed, and the design constraints include structural and aeroelastic limits. A comparative analysis of eight state-of-the-art algorithms is conducted to evaluate their performance in solving this problem. Among them, the Multi-Objective Multi-Verse Optimization (MOMVO) algorithm stands out, demonstrating superior performance and achieving the best results in the aeroelastic optimization task.

1. Introduction

Improving aircraft design technology for performance, limitations, and sustainability is increasingly required. These are important to reach new frontiers in commercial aircraft [1,2,3,4,5] and unmanned aerial vehicle (UAV) development [6,7]. To manufacture these aircraft, ensuring aircraft design is in accordance with aviation standards and regulations, such as those of the Federal Aviation Administration (FAA) and European Union Aviation Safety Agency (EASA) [8,9], is essential. Preliminary aircraft design is crucial because it is the center of the design process. Conceptual and detailed design can designate and actualize important aircraft components. At this design stage, aeroelasticity analysis is vital in achieving reliable design and has rudimentary design requirements. Aeroelastic optimization is essential for minimizing mass and construction costs while maximizing safety. Since most commercial transport aircraft operate at transonic speeds, aeroelastic optimization in this regime is essential.

Aeroelasticity has multiple mechanics applications that require understanding the interaction between elastic solid bodies and fluid. Aeroelastic analysis considers interactions of inertial, elastic, and aerodynamic forces in aircraft design. Both aerodynamic and structural forces are considered to analyze each interaction. Aerodynamic analysis is commonly conducted using panel methods that include medium-fidelity and classical techniques. The Vortex Lattice Method (VLM) [10,11,12,13,14,15] and the Doublet Lattice Method (DLM) [16,17,18,19] are especially important. Other methods are based on the lifting-line theory [20,21]. These aerodynamic methods can perform calculations rapidly; however, they cannot capture complex phenomena such as flow separation, oscillation, and shock, which commonly occur at transonic speeds. Therefore, a higher-fidelity method is needed. CFD [22] is usually employed to model complex problems and accurately capture aerodynamic phenomena, but it is highly time-consuming. Currently, both classical and CFD methods have been employed, and they have improved aerodynamic analysis for aircraft design [14,15,23,24]. Structural aircraft analysis is regularly performed with FEA, a classical structural engineering approach. Numerous aeroelastic analyses are required in the optimization process. For this reason, using CFD as an aerodynamic solver may be too time-consuming. Previous studies [25,26,27,28] demonstrated a correction technique that uses CFD results to improve DLM accuracy at transonic speeds.

Composite materials are commonly used to reduce the mass of aircraft. Normally, they are based on carbon or fiberglass. Not only does this material reduce mass, but tailoring the sequence stacking and orientation of the material layers can produce a designed strength that can withstand expected loads. The layout of the composite needs to adapt to the direction of stress. Conventionally, the direction of fibers in a ply is likely straight, parallel, or in a specified orientation in layered structures. The tow-steered technique is continuously applied to lay fiber with varying fiber orientations to address the stress field. This can be successful by employing automated fiber-placement technology. Presently, the design and optimization of composite tow steering involving aeroelastic analysis [29,30,31,32,33,34] have primarily focused on predefined layering configurations during manufacturing. These studies have contributed to innovative aircraft design utilizing composite materials.

Optimization has been performed in engineering problems to more rapidly determine optimal solutions than can be experimentally achieved, especially with aircraft [35,36,37]. These optimization methods often employ gradient-based [38,39] and stochastic-based methods [40,41]. Gradient-based and rapidly converging methods are normally applied, while stochastic-based methods use randomization to generate solutions that converge slowly. However, real-world problems may not be convex, so the gradients of objective and constraint functions may be unacceptable because of discrete variables. Likewise, the design variables are discrete in aeroelastic optimization for aircraft design. Stochastic-based methods are appropriate for handling discrete and continuous variables. Furthermore, non-gradient problems can be solved. This development requires metaheuristics (MHs) to be implemented for performance improvement in seeking optimal solutions to higher-complexity problems. Recently developed MHs can be categorized by their algorithms and classified into several groups. Evolutionary algorithms [42,43,44,45,46,47,48] are the most common MH type. Other forms, like swarm intelligence [49,50,51,52,53,54], mathematical and physical-inspired algorithms [55,56,57,58], and human society-inspired algorithms [59,60,61,62], are further MH types. The use of MHs has been studied in a variety of applications where the problem definitions can be optimized, typically using single and multiple objectives. MH algorithms are inspired by different techniques, but they are primarily focused on population-based algorithms, and their overall procedures are identical. Generally, initial solutions, called populations, are randomly generated in the optimization procedure. Afterward, the population reproduction becomes iterative by regulating the reproducing operators, which define the algorithm and motivation. In every iteration, the iterative population proceeds with domination sorting techniques to evaluate the non-dominated population. The calculations are terminated when the maximum number of function evaluations is reached. Non-dominated solutions can be structured to support the optimal Pareto front.

There are many recent studies in aeroelastic optimization employing MH for formulating and solving wing optimization problems [63,64,65,66,67,68,69,70,71,72,73,74,75,76,77]. Nevertheless, outdated MH algorithms like Genetic Algorithms (GA) [68,71,78,79], Particle Swarm Optimization (PSO) [68,80,81,82], and the Non-dominated Sorting Genetic Algorithm (NSGA-II) [83,84,85] have been used to solve these problems. Several studies utilized more recent optimization algorithms [67] and developed new algorithms [66]. Most aircraft optimization studies focus on designing optimal wings using isotropic or composite materials at low speeds. Studies have been performed examining aeroelastic optimization of a composite wing at transonic speeds [67], but none account for tow-steering. Recent works have explored tow-steered wing design [68,69,86,87,88,89,90,91,92,93]. In these studies, the primary design variables include laminate thickness and orientation, which are used to define the tow-steering pattern [92]. The primary objective of the optimization problems addressed in these studies is to minimize mass, while some have also incorporated multi-objective formulations. Typical constraints include allowable stresses, critical aeroelastic effects (such as divergence, flutter, and gust responses [68,69,87,91]), and mechanical failure modes (such as buckling [68,88,89] and fatigue [68]). However, the optimization of composite wings with tow steering at designed transonic speeds remains unexplored. Additionally, many existing studies rely on outdated optimization algorithms.

To address these gaps, our study aims to investigate the aeroelastic optimization of a tow-steered composite wing at transonic speed, utilizing aerodynamic transonic corrections from [28]. MHs offer superior global search capabilities compared to gradient-based methods, but their slower convergence is a drawback. To mitigate this, aerodynamic transonic correction is applied to enhance the aerodynamic influence coefficient (AIC) from the DLM using CFD results, as detailed in [67]. The contributions of this study are twofold:

• The development of a design procedure for tow-steered aeroelastic optimization at transonic speed.
• A comparative performance study of recent MHs in tow-steered aeroelastic optimization at transonic speed.

The research presented here introduces a novel composite wing design procedure and evaluates the optimization of performance for solving these challenges. This design procedure can be further developed, making it more practical and applicable for aircraft design, offering a foundation for future researchers to enhance its utilization. The comparative analysis of algorithms in this study demonstrates their effectiveness in solving the formulated problems, leading to the development of new algorithms and advancements in metaheuristic optimization.

A problem formulation and design variable encoding procedure specifically for tow-steering is detailed in Section 2, while the numerical test settings and experimental setup are presented in Section 3. The results and discussion of the study are given in Section 4, and study conclusions are offered in Section 5.

2. Aeroelastic Design Problem

Multiobjective optimization aims to simultaneously identify all possible solutions for a given set of scenarios. In this study, two objectives, mass and critical speed, are considered. The former is a cost indicator, while the latter reflects reliability. They tend to conflict. Thus, multiobjective design is a tool to explore all possible optimal solutions. Notably, all the designs are feasible, but they have varying levels of construction costs and reliability.

The practical benefit of obtaining these trade-off solutions is that designers can visualize all possible options related to the objective functions during the design process. This allows them to freely select a design for construction using decision-making techniques, such as a weight rating evaluation method. Moreover, selecting multiple designs for different scenarios is also viable, as all possible solutions are readily available.

In this work, recent and well-known multi-objective metaheuristic optimization algorithms are employed to solve the aeroelastic optimization problem of a tow-steered Common Research Model (CRM) wing [94,95,96,97]. The wing has a half-span of 26.23 m. Structural and aerodynamic analyses are performed using FEA and the DLM with AIC correction [28,67], respectively. The CFD results of this CRM wing model are from [96] and used to improve the aerodynamic AIC. Mindlin shell elements with a drilling degree of freedom [98] and wrapped correction [32] are applied. The transverse shear elasticity matrix for preventing a shear-locking phenomenon is computed based on [99,100,101]. The structural components, including the front spar, rear spar, and all ribs, are made of an isotropic material (Aluminum 2024-T3 [102]). Both the upper and lower skins consist of five layers of composite material (carbon fiber) [103]. The mechanical properties of the structural and skin components are presented in

Table 1. Material properties.

Material	Property	Value	Unit
Aluminum	Young’s modulus ($E$)	7.31 × 1010	Pa.
	Poisson’s ratio ($\vartheta$)	0.33	-
	Density ($\rho$)	2780	kg/m3
	Yield stress (${\sigma }_y$)	4.10 × 107	Pa.
	Shear modulus ($G$)	2.70 × 1010	Pa.
Carbon fiber	Young’s modulus ($E_{11}$)	1.35 × 1011	Pa.
	Young’s modulus ($E_{22}$)	1.00 × 1010	Pa.
	Shear modulus ($G_{12}$)	5.00 × 109	Pa.
	Shear modulus ($G_{13}$)	5.00 × 109	Pa.
	Shear modulus ($G_{23}$)	5.00 × 109	Pa.
	Poisson’s ratio ($\vartheta$)	0.3	-
	Density ($\rho$)	1600	kg/m3
	Tensile strength ($S_{ty,11}$)	1.50 × 109	Pa.
	Tensile strength ($S_{ty,22}$)	5.00 × 107	Pa.
	Tensile shear strength ($S_{ty,12}$)	7.00 × 107	Pa.
	Compressive strength ($S_{cy,11}$)	1.20 × 109	Pa.
	Compressive strength ($S_{cy,22}$)	2.50 × 108	Pa.
	Compressive shear strength ($S_{cy,12}$)	7.00 × 107	Pa.
	Stress of bonding ($S_b$)	1.00 × 107	Pa.

.

The wing FEA model used for structural analysis is shown in Figure 1. Since this is a preliminary design phase, a simplified wing mesh is used to reduce computational effort. The design variables for tow-steered layup orientation are presented in Figure 2 to address variations in carbon fiber orientation. Two design variables, root (${\theta }_{in}$) and tip (${\theta }_{tip}$ or ${\theta }_{out}$) angles, are assigned. The angles of carbon fiber are linearly distributed along the spanwise direction. The ply orientation of each element can be computed using Algorithm 1.

Algorithm 1. Calculation of orientation distribution of a tow-steered wing.

Input: orientations at the root chord (${Ori}_{root}$), orientations at the tip chord (${Ori}_{tip}$), FEA element data ($element$)       - reshape the FEA element data ($element$) into groups of span elements (${ele}_{span}$)     - get the total span elements (${nele}_{span}$), minimum span location ($y_{min}$), maximum span location ($y_{max}$)       ${fori1=1:1:nele}_{span}$         Compute span location ($y$) $${Ori}_{@span}\left(i1\right)={Ori}_{root}+\left({Ori}_{tip}-{Ori}_{root}\right)\ast \left({\displaystyle \frac{y-y_{min}}{y_{max}-y_{min}}}\right)$$ (1)       $end$   Output: ${Ori}_{@span}$

The CRM wing is simulated at transonic speed at a 0.85 Mach number (corresponding to a flight speed of 251.26 m/s). The wing is analyzed at a 2° angle of attack (AOA) and a 0.0155272 kg/m³ air density. The reference length of the wing is 7.0053 m. In this proposed problem, mass and critical speed are objective functions. A lower mass results in a lighter structure and reduced construction costs, while a higher critical speed contributes to a more reliable and safer design. The critical speed is defined as the lower value of the divergence and flutter speeds, representing static and dynamic aeroelastic phenomena, respectively. The goal of this multi-objective optimization is to obtain the optimal solutions or an optimal Pareto front that demonstrates a trade-off between construction costs and safety. The aeroelastic multi-objective optimization problem for the CRM wing is as follows:

$$Minimize:f\left(\mathit{x}\right)=[f_1\left(\mathit{x}\right),f_2\left(\mathit{x}\right)]$$

(2)

$$f_1\left(\mathit{x}\right)=totalwingmass\left(W\right),f_2\left(\mathit{x}\right)=-Criticalspeed\left(V_{cr}\right)$$

subject to

Lift effectiveness:	$0.8\le {\eta }_L\le 1.2;$
Critical speed:	$V_{cr}\ge 1.2\times V_{\infty };$
Displacement:	$u_{max}\le 10\%$ of half span;
Stress:	${\sigma }_{von,max}\le {\sigma }_y$ (isotropic material);
	${TH}_{max}\le 1$ (composite material);
Buckling:	${\alpha }_{bf}\le 1;$
Boundary:	${\mathit{x}}_L\le \mathit{x}\le {\mathit{x}}_U.$

where

$$\begin{array}{l}{\sigma }_{von}=\sqrt{{\sigma }_1^2-{\sigma }_1{\sigma }_2+{\sigma }_2^2}\\ TH=-{\left({\displaystyle \frac{{\sigma }_{11}}{{\mathbf{X}}_{11}}}\right)}^2+\left({\displaystyle \frac{{\sigma }_{11}{\sigma }_{22}}{{\mathbf{X}}_{11}^2}}\right)-{\left({\displaystyle \frac{{\sigma }_{22}}{{\mathbf{X}}_{22}}}\right)}^2-{\left({\displaystyle \frac{{\tau }_{12}}{{\mathbf{S}}_{12}}}\right)}^2\\ {\eta }_L={\displaystyle \frac{L_F}{L_R}}={\displaystyle \frac{q\left\{S\right\}{\left[AIC\right]}_F{\left\{{\alpha }_i\right\}}_F}{L_R}}\\ \mathit{x} \mathrm{is} \mathrm{a} \mathrm{vector} \mathrm{of} \mathrm{design} \mathrm{variables}.\end{array}$$

The aeroelastic problem has a total of 46 design variables, including the thicknesses and orientations of wing components, as listed below:

${\mathit{x}}_1$	=	The root thickness of the front spar (spar No. 1);
${\mathit{x}}_2$	=	The root thickness of the rear spar (spar No. 2);
${\mathit{x}}_3$	=	The root thickness of the ribs;
${\mathit{x}}_4$	=	The tip thickness ratio of the front spar (spar no. 1);
${\mathit{x}}_5$	=	The tip thickness ratio of the rear spar (spar No. 2);
${\mathit{x}}_6$	=	The tip thickness ratio of the ribs;
${\mathit{x}}_{7-11}$	=	The thicknesses at the root chord of the upper skin (the inner layer to the outer layer);
${\mathit{x}}_{12-16}$	=	The thicknesses ratio at the tip chord of the upper skin (the inner layer to the outer layer);
${\mathit{x}}_{17-21}$	=	Orientations at the root chord of the upper skin (the inner layer to the outer layer);
${\mathit{x}}_{22-26}$	=	Orientations at the tip chord of the upper skin (the inner layer to the outer layer);
${\mathit{x}}_{27-31}$	=	The thicknesses at the root chord of the lower skin (the inner layer to the outer layer);
${\mathit{x}}_{32-36}$	=	The thicknesses ratio at the tip chord of the lower skin (the inner layer to the outer layer);
${\mathit{x}}_{37-41}$	=	Orientations at the root chord of the lower skin (the inner layer to the outer layer);
${\mathit{x}}_{42-46}$	=	Orientations at the tip chord of the lower skin (the inner layer to the outer layer).

All thicknesses are linearly distributed, with the maximum thickness at the wing root, decreasing in the spanwise direction, and reaching a minimum at the wing tip. The thickness ratio represents the tip thickness to the root thickness ratio. Most design variables are discrete. The thicknesses of the composite material are selected from the set {0.25, 0.5, 0.75, 1.0, …, 9.75, 10.0} mm, while the list of ply orientations is {0, 5, 10, …, 170, 175} degrees. The ribs and spars are made of an isotropic material, with thicknesses chosen from the set {0.25, 0.5, 0.75, …, 4.75, 5.0, 6.0, 7.0, …, 19, 20, 25, …, 45.0, 50.0} mm.

The function evaluation process for the aeroelastic wing optimization problem is outlined in Algorithms 1 and 2. Since the wing shape remains unchanged during optimization, the AIC matrix is pre-calculated using the DLM. Design variables containing all structural data are then input to perform FEA. Displacements, stresses, and buckling factors are obtained at this stage. Then, aeroelastic analysis is conducted to determine the lift effectiveness and critical speed.

Algorithm 2. The tow-steered CRM wing aeroelastic function evaluation

Input: design variable ($\mathit{x}$)       - Load the aerodynamic results data including     - Transform the design variable into the real geometric values including skin, rib and spar thicknesses and orientation of tow steer composite wing skin (detailed in Algorithm 1).     - Create the finite element model, assemble the stiffness matrix and mass matrix     - Solve the free vibration analysis and static analysis:         -     Free vibration analysis results including mode shape, lambda, omega.         -     Static analysis results including the structural deformation ($u_{max}$), Tsai-Hill criteria stress (${TH}_{max}$) and Von-Mises stress (${\sigma }_{von,max}$).     - Solve the static aeroelastic: obtain the lift effectiveness (${\eta }_L$) and divergence speed ($V_d$)     - Solve for the flutter speed ($V_f$) using v-g-method     - Compute the critical speed: $$V_{cr}=\mathrm{m}\mathrm{i}\mathrm{n}(V_d,V_f)$$ (3)     - Perform the objective function and constraints         -     $f\left(\mathit{x}\right)=\left[f_1\left(\mathit{x}\right),f_2\left(\mathit{x}\right)\right],f_1\left(\mathit{x}\right)=total wing mass,$         -     $f_2\left(\mathit{x}\right)=-Critical speed$       Subject to         -     $g_1\left(\mathit{x}\right):u_{max}\le 2.623$ m (10% of half span)         -     $g_2\left(\mathit{x}\right):{\sigma }_{von,max}\le 410$ MPa (isotropic material)         -     $g_3\left(\mathit{x}\right):{TH}_{max}\le 1$ (composite material)         -     $g_4\left(\mathit{x}\right):{\alpha }_{bf}\ge 1$ (Buckling factor constraint)         -     $g_5\left(\mathit{x}\right):V_{cr}\ge 301.512$ m/s (1.2 times flight speed)         -     $g_6\left(\mathit{x}\right):0.8\le {\eta }_L\le 1.2$         -     ${\mathit{x}}_L\le \mathit{x}\le {\mathit{x}}_U$   Output:  $f\left(\mathit{x}\right)$, $g_i\left(\mathit{x}\right)$

3. Numerical Setup

In the aeroelastic analysis procedure, the FEA mesh is completely fixed in both static and dynamic analyses. A grid-independence test is conducted. In this study, we deliberately chose a coarser mesh over a fully converged one since this research centers on crafting a design approach specialized for enhancing the tow-steered composite wing structures. Thus, achieving absolute precision in analytical outcomes is not the primary aim.

The integration of the function evaluation and an optimization algorithm is the key to this study. The common procedure for multi-objective metaheuristics (MHs) is illustrated in Figure 3. At the beginning of the optimization process, a population or a set of design solutions is randomly generated, while the Pareto archive, a matrix for collecting non-dominated solutions, is generated. Then, a new population is iteratively derived using reproduction operators, which vary between algorithms. Subsequently, function evaluations are performed. The members in the Pareto archive and those in a new population are sorted using a non-dominated sorting technique to identify non-dominated solutions and saved to the Pareto archive.

When the archive size is too large, a clustering technique is employed to remove some solutions to conserve computer memory. This clustering technique preserves the diversity of the members in the archive while addressing memory constraints. The updated Pareto archive at the last iteration is regarded as the Pareto set of optimal solutions. An equal number of function evaluations is used to terminate all algorithms to ensure that this comparison is fair.

During the optimization process, the wing shape remains fixed. The FEA mesh, aerodynamic loadings, and AIC are pre-calculated and imported for structural and aeroelastic analysis. Aerodynamic loadings are derived using the VLM, while the AIC is constructed using the DLM with steady-state corrections based on transonic CFD results.

Each reproduced population contains solution vectors used during function evaluations. These vectors contain structural data such as material thicknesses and orientations. The structural parameters, combined with the imported FEA mesh, are used to generate the FEA model. Pre-calculated aerodynamic loadings are then applied to perform structural static analysis. Subsequently, aeroelastic analysis is conducted using the FEA model and the imported AIC to determine divergence and flutter speeds. The critical speed is the lower of these two values.

The performance of the state-of-the-art algorithms on multi-objective aeroelastic optimization of the CRM composite wing is investigated. They are evaluated based on the convergence rate, consistency, and versatility. The hypervolume (HV), which can measure both the spread and advancement of a Pareto front, is used as a performance indicator in this study [104,105]. The statistical parameters of the HV, including the average, standard deviation, maximum, minimum, and rank of average, are evaluated for performance comparison. The average HV, maximum, and rank are used as convergence rate indicators, while the standard deviation and minimum HV are consistency indicators. Higher values of average, maximum, and minimum HV are better, while lower values of standard deviation and rank are more favorable. Ten independent runs of each algorithm are evaluated in the aeroelastic optimization problem. The population size and maximum function evaluations are 100 and 20,000, respectively. The parameter settings of all algorithms are detailed in

Table 2. Parameter settings of MHs.

Algorithm	Parameter Name	Abbreviation	Value
Multobjective Mayfly Optimization Algorithm (MOMA)	Inertia weight	$g$	0.8
	Inertia weight damping ratio	$g_{damp}$	1
	Personal learning coefficient	$a_1$	1
	Global learning coefficient	$a_2$ & $a_3$	1.5
	Distance sight coefficient	$\beta$	2
	Mutation coefficient	$d$	0.77
	Mutation coefficient damping ratio	$d_{damp}$	0.99
	Random flight	$f$	0.77
	Random flight damping ratio	$f_{damp}$	0.99
	Mutation rate	$\mu$	0.02
Multiobjective Salp Swarm Algorithm (MSSA)	-	-	-
Multiobjective Grasshopper Optimization Algorithm (MOGOA)	Exploration and exploitation proportional	$C$	[0.00004, 1]
Multiobjective Multi-Verse Optimization (MOMVO)	Wormhole existence probability	$WEP$	[0.2, 1]
Multiobjective Lichtenberg Algorithm (MOLA)	Reference Lichtenberg point	$ref$	0.4
	Number of particles	$N_P$	100,000
	Creation radius	$R_c$	150
	Stick probability	$S_c$	1
	Number of grids in each dimension	$NGrid$	30
Multiobjective Jellyfish Search (MOJS)	Hypercube limitation grid point	$NGrid$	20
Multiobjective Crystal Structure Algorithm (MOCryStAl)	Grid inflation parameter	$\alpha$	0.1
	Number of grids per each dimension	$NGrid$	30
	Leader selection pressure parameter	$\beta$	4
	Extra repository member selection pressure	$\gamma$	2
Multiobjective Manta Ray Foraging Optimizer (MOMRFO)	Grid inflation parameter	$\alpha$	0.1
	Number of grids per each dimension	$NGrid$	30
	Leader selection pressure parameter	$\beta$	4
	Extra archive member selection pressure	$\gamma$	2

.

In the aeroelastic analysis procedure, the FEA mesh is fully fixed for both static and dynamic analyses. A grid independence test is conducted to ensure accuracy. In this study, we intentionally selected a coarser mesh rather than a fully converged one because our primary focus is on developing a design methodology to enhance tow-steered composite wing structures. So, achieving absolute precision in analytical results is not the primary objective. The performance of state-of-the-art algorithms for multi-objective aeroelastic optimization of the CRM composite wing is investigated. Algorithms are evaluated based on the convergence rate and consistency. Hypervolume (HV), which quantifies both the spread and advancement of the Pareto front, is used as a performance metric. Statistical results, including the average, standard deviation, maximum, minimum, and rank of HV, are analyzed for comparative assessment.

The average HV, maximum HV, and rank indicate the convergence rate, while the standard deviation and minimum HV reflect consistency. Higher average and maximum HV values, along with a lower rank, indicate better performance. Lower standard deviation and higher minimum HV values signify greater consistency and robustness. Each algorithm is evaluated over 10 independent runs for the aeroelastic optimization problem. The population size is uniformly set to 100 across all algorithms. Each algorithm is subjected to 20,000 function evaluations to ensure a fair comparison. Parameter settings, adopted from pre-tuned configurations specified in earlier research, are detailed in Table 2. These settings represent optimal values tested across various problems in the respective studies.

4. Results and Discussion

In this section, the optimum solutions from several algorithms are presented and discussed. HV is used as an indicator of the diversity and spread of the Pareto solution. Figure 4 presents HV convergence plots for all algorithms, comparing the average hypervolume results from all optimization runs. MOMVO outperforms all other algorithms, while MOGOA, MOMA, and MSSA are the next best. Early in the optimization, MOGOA demonstrates the fastest convergence compared to other approaches. However, its HV shows minimal improvement during the second half of the optimization runs. This suggests that MOGOA excels in exploitation but suffers from poor exploration, leading to premature convergence at local optima. In contrast, MOMVO exhibits outstanding performance in the latter half of the optimization, with its rapidly increasing hypervolume. By the end of the process, it achieves nearly double the HV of the worst-performing algorithm, MOCryStAl. These results highlight MOMVO as the best-performing algorithm, striking a well-balanced trade-off between exploitation and exploration. Another notable observation from the convergence plots is the inconsistent convergence pattern of MOMA. This inconsistency suggests that the solutions provided by MOMA cannot be considered optimal, as the algorithm did not fully converge.

A comparison of the best Pareto fronts for this test problem is shown in Figure 5, which corresponds to the best HV results in

Table 3. Statistical results of a tow-steered aeroelastic problem considering a hypervolume indicator.

Algorithm	MOCryStAl	MOGOA	MOJS	MOLA	MOMA	MOMRFO	MOMVO	MSSA
average	0.42729	0.47462	0.52441	0.43754	0.56144	0.40638	0.65495	0.44987
standard deviation	0.06488	0.23092	0.03094	0.06618	0.05834	0.15415	0.07778	0.18921
maximum	0.50795	0.67949	0.56871	0.52411	0.65702	0.52966	0.76748	0.62961
minimum	0.29419	0	0.46314	0.29234	0.470095	0	0.46145	0
rank	7	4	3	6	2	8	1	5

Bold text presents the average, best, worst, standard deviation, and rank of hypervolume indicator.

. The best Pareto front in this study is obtained by MOMVO. Upon examining the Pareto fronts in Figure 5, the results can be categorized into two groups based on the critical speed. The first group includes Pareto fronts with critical speeds ranging from 500 m/s to approximately 600 m/s, while the second group features critical speeds greater than 700 m/s. However, regardless of the group, the Pareto front obtained by MOMVO exhibits significantly greater spread and advancement compared to all other examined approaches. Considering all aspects of the results, MOMVO is the best-performing algorithm for this problem, while MOGOA and MOMA provide comparable results and show potential for further improvement.

Statistical HV results are presented in Table 3 and summarized as a box plot in Figure 6 for a deeper analysis of performance in terms of both speed and consistency. In terms of speed, MOMVO achieved the best average HV and ranking, as well as the highest hypervolume in this study, confirming its position as the fastest algorithm. MOJS achieved the lowest standard deviation, while MOMA delivered the best results in the worst-case scenario. The poorest MOMA run still had the greatest hypervolume, slightly better than the worst run of MOMVO. Although MOJS exhibited a lower standard deviation, both its average hypervolume and that of MOMA were significantly lower than those of MOMVO. When comparing the average hypervolumes, MOMVO, MOJS, and MOMA have standard deviations of 11.9, 5.9, and 10.4%, respectively. So, MOJS demonstrates the best consistency, while MOMVO and MOMA are comparable. However, in the worst-case scenario, algorithms such as MOGOA, MOMRFO, and MSSA performed poorly, with a minimum hypervolume value of zero, indicating their failure to capture any feasible solutions for the test problem.

Overall, MOMVO remains the best algorithm based on these statistical results, as it offers the fastest speed and is one of the most consistent algorithms in this study. However, it is noteworthy that MOGOA and MOMA produce comparable results. Based on the observations in this study, the performance of these algorithms can be enhanced in future research.

MOGOA exhibits rapid convergence during the early stages of the optimization, indicating strong local search capabilities. However, its poor global search capability causes convergence at local optima. This algorithm could be improved by enhancing its global search capabilities. Potential improvements include introducing additional reproduction operators to increase population diversity or integrating adaptive strategies to dynamically readjust reproduction parameters during optimization.

MOMA demonstrates comparable performance to MOGOA. It maintains a good balance between exploration and exploitation, as it continuously improves the hypervolume throughout the optimization. However, MOMA encounters an issue in its clustering process. When the number of non-dominated solutions exceeds the maximum Pareto archive limit, a clustering technique is used to eliminate excessive solutions. If the algorithm employs an effective clustering technique, the hypervolume should consistently increase rather than fluctuate, as observed in the hypervolume history of MOMA. Its performance could be enhanced by addressing this issue.

The details of the selected solutions from the best Pareto front, highlighted in Figure 7, are provided in

Table 4. Aeroelastic results of three optimal solutions sampled from the Pareto front.

Optimum Solution Number	No. 1	No. 2	No. 3
1st objective total mass [kg.]	32,639.14439	44,057.84840	69,370.10817
2nd objective critical speed [m/s]	511.90569	832.25059	876.06638
total mass [kg.]	32,639.14439	44,057.84840	69,370.10817
divergent speed [m/s]	65,535	65,535	65,535
flutter speed [m/s]	511.90569	832.25059	876.06638
maximum deflection [m.]	0.02992	0.02639	0.01727
maximum Von-Misses stress [Pa.]	200,433,275.012	142,468,641.420	58,506,362.433
maximum Tsai Hill [-]	0.00084	0.00063	0.00028
lift effectiveness [-]	1.03190	1.02289	1.01704
minimum buckling factor [-]	1.15896	1.02324	1.73896

and

Table 5. Design variables of three optimal solutions sampled from Pareto front.

Optimum Solution Number		Solution No. 1	Solution No. 2	Solution No. 3
1st objective total mass [kg.]		32,639.14439	44,057.84840	69,370.10817
2nd objective critical speed [m/s]		511.9056913	832.2505855	876.0663793
thickness at root chord	leading spar [mm.]	8	7.5	7
	tailing spar [mm.]	8.5	17	29
	ribs [mm.]	1	1.5	4.5
thickness at tip chord	leading spar [mm.]	2.3	0.7	1.4
	tailing spar [mm.]	6.6	13.6	7.6
	ribs [mm.]	1	1.4	4.3
thickness at upper skin root chord	layer 1, inner layer [mm.]	1	1	2
	layer 2 [mm.]	9	9	8
	layer 3, middle layer [mm.]	0.2	0.2	35
	layer 4 [mm.]	0.4	1	11
	layer 5, outer layer [mm.]	20	11	23
thickness at upper skin tip chord	layer 1, inner layer [mm.]	0.2	0.2	0.3
	layer 2 [mm.]	6.1	7.2	5.6
	layer 3, middle layer [mm.]	0.1	0.1	6.4
	layer 4 [mm.]	0.2	0.5	5.8
	layer 5, outer layer [mm.]	12.8	10.1	15.4
orientation at upper skin root chord	layer 1, inner layer [mm.]	100	105	105
	layer 2 [mm.]	135	130	140
	layer 3, middle layer [mm.]	125	155	150
	layer 4 [mm.]	165	165	120
	layer 5, outer layer [mm.]	30	45	35
orientation at upper skin tip chord	layer 1, inner layer [mm.]	145	145	120
	layer 2 [mm.]	120	130	155
	layer 3, middle layer [mm.]	145	160	165
	layer 4 [mm.]	45	35	45
	layer 5, outer layer [mm.]	20	0	0
thickness at lower skin root chord	layer 1, inner layer [mm.]	0.2	1	8
	layer 2 [mm.]	18	30	29
	layer 3, middle layer [mm.]	12	16	16
	layer 4 [mm.]	19	31	40
	layer 5, outer layer [mm.]	7	7	12
thickness at lower skin tip chord	layer 1, inner layer [mm.]	0.2	1	7.8
	layer 2 [mm.]	7.7	29.6	27.8
	layer 3, middle layer [mm.]	8.7	11.4	9.1
	layer 4 [mm.]	18.7	31	40
	layer 5, outer layer [mm.]	2.9	4.2	5
orientation at lower skin root chord	layer 1, inner layer [mm.]	130	125	85
	layer 2 [mm.]	115	110	145
	layer 3, middle layer [mm.]	10	15	10
	layer 4 [mm.]	150	150	150
	layer 5, outer layer [mm.]	165	160	115
orientation at lower skin tip chord	layer 1, inner layer [mm.]	145	145	145
	layer 2 [mm.]	75	50	45
	layer 3, middle layer [mm.]	15	0	0
	layer 4 [mm.]	115	115	120
	layer 5, outer layer [mm.]	30	30	5

. Table 4 presents the values of all objective and constraint functions for the selected solutions, while Table 5 provides the corresponding design variable values. These details are provided to demonstrate the trade-off between mass and critical speed. In Table 5 and Figure 8, the first solution has less than half the mass of the third, but it has a significantly lower critical speed. The first solution is lighter but less reliable than the third solution, while the second solution is intermediate.

All selected solutions are feasible, as all constraints—including displacement, von Mises stress (for aluminum), Tsai–Hill (for carbon fiber), buckling factor, critical speed, and lift effectiveness—are within the allowable ranges. Additionally, the trend shows that overall thickness increases with greater mass and critical speed. This trade-off reflects the goal of multi-objective design optimization. Thicker structures may provide greater safety but require higher manufacturing costs.

In Table 5 and Figure 8, the details of the thicknesses and orientations of the selected solutions are displayed. Notably, the color scale in each sub-figure is different, as each solution thickness varies considerably. The solutions, from the first to third, range from the lightest to the heaviest. The optimum orientation of each solution is not perfectly matched. However, the orientations of each solution show a similar trend in more than half of the layers. This suggests that the optimum orientation may vary slightly depending on differences in ply thickness.

A comparison of composite wing aeroelastic results [67] with tow-steered aeroelastic results is presented using the best Pareto fronts from both studies, as shown in Figure 9. This figure highlights three solutions: (1) the tow-steered composite wing from the current study with 20,000 FEs, (2) the tow-steered composite wing from the current study with 5000 FEs, and (3) the Pareto front of composite wing aeroelastic results from previous studies [67] with 5000 FEs. The Pareto front at 5000 FEs with tow-steer design is more advanced than that obtained from constant ply orientation. The Pareto front at 20,000 FEs of the two-steer design dominates that of a previous study.

Since all optimal trade-off solutions have been obtained, any design can be selected for construction using a decision-making technique, such as a weight rating evaluation method. This selection can be made by designers or customers. Moreover, a comparison between the original transonic design problem presented in previous work [67] and the current study indicates significant improvements, as illustrated in Figure 9. The tow-steered composite wing design solution demonstrates superior performance to the original problem. As expected, more optimal solutions are achievable from the tow-steered design concept, which provides greater design diversity.

5. Conclusions

A design procedure for tow-steered aeroelastic optimization is developed in the current study. A comparative analysis of recent multi-objective algorithms for the multi-objective aeroelastic optimization of a tow-steered composite wing is conducted. These recent algorithms, which are rarely used in the aerospace field, are employed to solve a complex, real-world aircraft design optimization problem. Their performance and consistency are investigated.

Overall, MOMVO is considered the best algorithm in this study, with MOGOA, MOJS, and MOMA having slightly lower performance. MOMVO is the fastest algorithm, while MOGOA ranks second in this aspect. MOJS is the most consistent, slightly outperforming MOMVO and MOMA. MOMA performs best in the worst-case scenario, slightly better than MOMVO. While MOGOA demonstrates good speed, it fails in the worst-case scenario as it cannot obtain feasible solutions in some cases.

MOMVO is the best-performing algorithm in this study, with outstanding speed and best results in both consistency and the worst-case scenario. Optimal solutions showcasing the trade-off have been identified, allowing designers or customers to select a design for construction using decision-making techniques. Additionally, the tow-steered composite wing design from this study outperforms the original transonic design problem [67], significantly increasing the potential for achieving optimal solutions.

For future research, we recommend three key steps: (1) incorporating nonlinearity effects for high aspect ratio wings and considering additional flight conditions, for example, the maneuver load and gust load cases, among others; (2) dealing with randomness in linear/nonlinear computational simulation and manufacturing employing automated fiber placement and material properties using reliability optimization; (3) developing a faster optimizer by integrating machine learning or deep learning models with metaheuristic algorithms to address more detailed problems; and (4) conducting further experiments to confirm the reliability of the optimized designs. Ground vibration tests and updates to FEA models can be experimentally performed to apply the optimized design to a physical wing in future work.