Optimization of Airfoils for the Design of Long Endurance Aircrafts Using Deep Learning Models and Metaheuristics Algorithms

1. Introduction

As a key component of aircraft, an excellent aerodynamic configuration of the airfoil plays a crucial role in fuel economy, reducing polluting emissions and improving performance [1,2,3]. In the design of unmanned aerial vehicles (UAV) and gliders, the development of wing configurations with a high aspect ratio is also being opted for, because these wing configurations provide a better range and flight duration. In UAVs there is a greater interest in improving the flight duration, especially in subsonic flight missions, in missions such as: reconnaissance, target tracking, pseudo-satellites, etc [4].

Taking into account the aircraft powered by propeller engines, aerodynamically, an aircraft to achieve a large endurance its C_L^1.5/C_D (C_L, lift coefficient; C_D, drag coefficient) ratio values should be the highest possible [5]. This C_L^1.5/C_D ratio is usually called the endurance parameter. As mentioned above, this is mainly achieved by having a wing with a high aspect ratio, but the shape of the airfoils along the wingspan is also of utmost importance, which should also provide a high c_l^1.5/c_d value [4,6]. The design of this type of airfoils is obtained by optimizing the shape of the airfoil, which must produce a high value of the endurance parameter with respect to a design lift coefficient. It is important to mention that airfoils with high values of c_l^1.5/c_d and c_l usually have convex shapes that generate very negative values of pitch moment coefficient (c_m). This drawback is compensated by the empennage of the aircraft. Although in aircraft configurations that lack a horizontal empennage (such as flying wings) it is necessary that the c_m be as close to zero, for this reason in these cases an extra objective function is usually added to the optimization task, minimizing |c_m| [4,7].

Computational fluid dynamics (CFD), wind tunnel testing and theoretical analysis have gradually become the main approaches to aerodynamic design, advancing the progress of aerodynamic optimization of airfoils. However, in real engineering design, wind tunnel tests are expensive, and theoretical analysis is less applicable when dealing with complex engineering problems, so CFD has gradually become the main method for aerodynamic analysis [8].

The basic element of aerodynamic optimization of the airfoil mainly includes the design object, the design objective, the design constraints and the design method (optimization method). Design object refers to the aerodynamic geometric configuration. The design objective refers to the expected aerodynamic performance, the characteristics of the flow field, etc. The design constraints illustrate the interdependence and constraints between the design variables of the design object. And optimization methods provide the strategies and means to achieve the design goal [8].

Currently there are a large number of optimization methods that are widely used. These methods can be generally classified into gradient-based and gradient-free optimization methods. Current gradient-based approaches use discrete adjunct sensitivity analysis, scale approximately linearly with the number of design variables, and are very capable of handling thousands of design variables and constraints. The advantages of gradient-based optimization to address airfoil shape optimization are low computational costs in large design spaces and a track record of successful implementations in aerodynamics. However, there are some drawbacks to this approach, such as the tendency to converge on local optima and a high sensitivity to the starting point, poor efficiency for nonlinear cost functions and the need for continuous shapes (the shape gradient must exist at all points). The gradient-based approach may be more appropriate for a detailed aerodynamic design, as they are only able to offer a limited range of solutions. Gradient-free methods can be more complex to implement than gradient-based methods, but they do not require continuity or predictability in the design space and usually increase the probability of finding a global optimum. Optimization methods known as metaheuristics can offer robust methods for finding a solution and increase the probability of converging on a solution at the global optimum. It is known that these gradient-free methods are able to address numerically noisy optimization problems that are difficult for gradient-based methods. And unlike gradient methods, derivatives of cost functions are not necessary. In addition, no predefined baseline design or knowledge of the design space is required and gradient-free methods usually optimize several solutions in parallel. However, their convergence speeds are usually slow enough that they are unable to cooperate directly with time-consuming CFD analysis. To counteract this disadvantage that metaheuristic methods have, in the first instance, more robust computing equipment is used (which implies raising the costs of the optimization process). Another alternative is to make use of methods based on surrogate models [8,9,10,11].

Traditional methods for training surrogate models include the polynomial response surface method, the support vector regression method, and the Kriging model. The traditional substitute model has adjustable parameters; however, the main problem is that it is not suitable for handling large-scale training data, so it is usually trained with a small amount of data in a relatively limited space [8,9,10,11,12]. Recent advances in machine learning (especially in deep learning methods: artificial neural networks) have provided a new option for generating surrogate models. Compared to traditional surrogate models, the advantages of applying neural networks in aerodynamic data modeling lie in the fact that: a) neural networks, as a data-driven model, are not based on aerodynamic theories or physical models, which means that it does not require in-depth knowledge of aerodynamics; b) neural networks can be used to address high-dimensional, multi-scale and nonlinear problems that are difficult for traditional surrogate models; and c) some neural networks have the ability to process time sequence data [13].

In the present work, a methodology is developed to accelerate the process of airfoil optimization under conditions of subsonic, incompressible and turbulent flow, making use of current deep learning models and metaheuristic algorithms. This research work is specifically aimed at the design of aircraft with long endurance operating in subsonic flight regime.

2. Materials and Methods

2.1. Mathematical Formulation of the Problems

Two cases are considered:

• determination of the airfoil shape producing a long endurance parameter (c_l^1.5/c_d) with respect to a specified lift coefficient;
• determination of the airfoil shape producing a long endurance parameter (cl1.5/cd) and a near-zero value of the pitch moment coefficient (cm) with respect to a specified lift coefficient.

The purpose of the above cases is to contemplate aircraft configurations with empennage and without empennage.

For the first case, it is only necessary to maximize the endurance parameter and consider design constraints such as the design lift coefficient.

max f₁(ξ) = c_l^1.5/c_d,

such that h_a(ξ) = 0, a = 1, 2, ..., A,

(1)

g_b(ξ) ≤ 0, b = 1, 2, ..., B,

ξ ∊ Ω.

For the second case, the objective of minimizing the magnitude of c_m is added.

max f₁(ξ) = c_l^1.5/c_d,

min f₂(ξ) = |c_m|

such that h_a(ξ) = 0, a = 1, 2, ..., A,

(2)

g_b(ξ) ≤ 0, b = 1, 2, ..., B,

ξ ∊ Ω.

where ξ is the vector of design parameters; h(ξ) are constraint functions (equalities); g(ξ) are constraint functions(inequalities); Ω is the decision space of the design parameters.

To provide a solution to the above optimization tasks, the following methodology is proposed:

• Selecting a method of parameterization of the aerodynamic profile, the method should allow modeling aerodynamic profiles used for the design of aircraft with long-span wings;
• Select a robust method that allows the aerodynamic coefficients of the airfoils to be calculated. The method should take into account the viscosity and turbulence effects of the flow;
• Propose a model based on deep learning models to predict the values of the objective functions, this in order to speed up the process of calculating the objective functions;
• Finally, select the appropriate optimization algorithms for each task.

2.2. Selection of the Method of Airfoils Parameterization

To determine the design parameters, two of the best airfoil parameterization methods were evaluated: CST (see Appendix A) and Bezier-PARSEC (see Appendix B) methods. The way to evaluate the methods is to evaluate the capabilities they have when reconstructing the geometry existing airfoils (such as NACA, Eppler, TSAGI, etc.). Obtaining the original coordinates of the airfoils is carried out using the profile database of the University of Illinois at Urbana-Champaign (UIUC) [14]. Only the selection of asymmetric airfoils that are used in general aviation (aircraft operating in subsonic and/or transonic flight regime) is contemplated [15].

The reconstruction of the airfoils is carried out through the use of the evolutionary algorithm Differential Evolution (DE), whose objective is to find the values of the design parameters of each method which provide the least geometric deviation when comparing the reconstructed profiles and the original ones [16]. The value of the geometric deviation is obtained by the error norm L₂:

$${‖{\mathrm{y}}_{\mathrm{n}}-{\mathrm{y}}_{\mathrm{t}\mathrm{a}\mathrm{r},\mathrm{n}}‖}_{{\mathrm{L}}_2}=\sqrt{{\sum }_{\mathrm{n}=1}^{\mathrm{P}}{\left({\mathrm{y}}_{\mathrm{n}}-{\mathrm{y}}_{\mathrm{t}\mathrm{a}\mathrm{r},\mathrm{n}}\right)}^2},$$

(3)

where y_tar is the set of heights of the original points of the target profile; y is the set of heights of the reconstructed points; P is the number of points that make up the geometry of the airfoil.

2.3. Creation of a Database of Airfoils and Their Aerodynamic Coefficients

2.3.1. Geometries of the Airfoils

In the first instance, it is proposed to make use of the selected profiles from the UIUC database, in case these are not enough, the use of a Generative Adversarial Network (GAN) is suggested for the creation of new profiles, which will be based on the selected profiles from the UIUC database. GANs use a learning scheme that allows the defining attributes of the probability distribution to be encoded in a neural network, making it possible to generate instances that resemble the original probability distribution [17]. The GAN architecture is formed by two constituent neural networks: one called discriminator (D) and the other generator (G). The G network is responsible for generating new instances of the same domain as that of the source dataset. The network D is in charge of discriminating whether the input data are real, that is belonging to the input data set or whether they are fictitious, that is artificially generated. Both networks are trained together in such a way that G maximizes its chances of not being detected by D and D in such a way that it makes its methods of detecting the artificially generated data by G more and more sophisticated. These two adversarial networks compete in a zero-sum game in which it is hypothesized that they eventually reach a Nash equilibrium (see Figure 1) [18].

2.3.2. Obtaining the Aerodynamic Coefficients

The aerodynamic coefficients of the profiles are obtained by means of Computational Fluid Dynamics, using the CFD package OpenFOAM 11 [19] and the GMSH mesh generator [20]. It is considered a viscous, turbulent, steady and incompressible flow.

$${\displaystyle \frac{\partial {\mathrm{u}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{i}}}}=0,$$

(4)

$${\mathrm{u}}_{\mathrm{j}}{\displaystyle \frac{\partial {\mathrm{u}}_{\mathrm{i}}}{\partial \mathrm{t}}}=-{\displaystyle \frac{1}{\mathsf{\rho }}}{\displaystyle \frac{\partial \mathrm{p}}{\partial {\mathrm{x}}_{\mathrm{i}}}}+\mathsf{\nu }{\displaystyle \frac{{\partial }^2{\mathrm{u}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}\partial {\mathrm{x}}_{\mathrm{j}}}}-{\displaystyle \frac{\partial {\mathrm{u}\mathrm{\text{'}}}_{\mathrm{i}}{\mathrm{u}\mathrm{\text{'}}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{j}}}},$$

(5)

where u_i is the average velocity; ν is the kinematic viscosity of the air; ρ is the density of the air; p is the pressure; u'_iu'_j is the Reynolds stress tensor.

To model the turbulence the k-ω SST model is considered. This is a two-equation turbulence model that is widely used to calculate the aerodynamic coefficients in subsonic and transonic flows. It combines the low Reynolds number version of the k-ω model and the high Reynolds number version of the k-ε model to obtain accurate results under a wide range of flow conditions. The k-ω SST model is computationally more expensive than other two-equation models, but can provide more accurate results for more complex flows [21,22,23]. OpenFOAM uses the version of the k-ω SST model exposed by Menter in [24]. The turbulence specific dissipation rate (ω) equation is given by:

$${\displaystyle \frac{\partial \left(\mathsf{\rho }\mathsf{\omega }\right)}{\partial \mathrm{t}}}+{\displaystyle \frac{\partial \left(\mathsf{\rho }{\mathrm{u}}_{\mathrm{i}}\mathsf{\omega }\right)}{\partial {\mathrm{x}}_{\mathrm{i}}}}=\mathsf{\alpha }\mathsf{\rho }{\mathrm{S}}^2-\mathsf{\beta }\mathsf{\rho }{\mathsf{\omega }}^2+{\displaystyle \frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}}\left[(\mathsf{\mu }+{\mathsf{\sigma }}_{\mathsf{\omega }1}{\mathsf{\mu }}_{\mathrm{t}}){\displaystyle \frac{\partial \mathsf{\omega }}{\partial {\mathrm{x}}_{\mathrm{i}}}}\right]+2(1-{\mathrm{F}}_1)\mathsf{\rho }{\mathsf{\sigma }}_{\mathsf{\omega }2}{\displaystyle \frac{1}{\mathsf{\omega }}}{\displaystyle \frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{i}}}}{\displaystyle \frac{\partial \mathsf{\omega }}{\partial {\mathrm{x}}_{\mathrm{i}}}},$$

(6)

and the kinetic energy of turbulence (k) is defined by:

$${\displaystyle \frac{\partial \left(\mathsf{\rho }\mathrm{k}\right)}{\partial \mathrm{t}}}+{\displaystyle \frac{\partial \left(\mathsf{\rho }{\mathrm{u}}_{\mathrm{i}}\mathrm{k}\right)}{\partial {\mathrm{x}}_{\mathrm{i}}}}={\tilde{\mathrm{P}}}_{\mathrm{k}}-{\mathsf{\beta }}^{\mathrm{*}}\mathsf{\rho }\mathrm{k}\mathsf{\omega }+{\displaystyle \frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}}\left[(\mathsf{\mu }+{\mathsf{\sigma }}_{\mathrm{k}1}{\mathsf{\mu }}_{\mathrm{t}}){\displaystyle \frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{i}}}}\right],$$

(7)

where F1 is the blending function; S is the invariant measure of the strain rate; ${\tilde{\mathrm{P}}}_{\mathrm{k}}$ is a production limiter. All constants are computed by a blend from the corresponding constants of the k - ε and the k - ω models via α = α₁F + α₂(1 - F). The constants for this model are: β* = 0.09, α₁ = 5/9, β = 3/40, σ_k1 = 0.85, σ_ω1=0.5, α₂ = 0.44, σ_ω2 = 0.856 [24].

The k-ω SST model allows the use of wall functions, which allows the use of meshes that are not so thin next to the wall (y+ values between 10 and 100) without losing as much precision [24]. By having this tolerance in the y+ values, greater flexibility is possible when developing the mesh.

For the creation of the mesh, a control volume with a type C configuration [25] (see Figure 2(a)) and the use of hexahedral elements is considered (OpenFOAM does not work with two-dimensional meshes). Each mesh has 21360 hexahedral elements, which allows y+ values close to 10 when using turbulent flows with Re values between 1 million to 6 million (see Figure 2(b)). Most of the boundary conditions used are shown in Figure. 2(c) . The "wall" condition is applied to the faces that limit the geometry of the profile; the "outlet" condition is assigned to the vertical face located to the right of the control volume; the "inlet" condition is assigned to the arc faces and the upper and lower horizontal faces of the control volume; finally, the "empty" condition is applied to the front and rear faces of the control volume. This last condition is what indicates to OpenFOAM that the simulation is two-dimensional [19]. The entire control volume is defined as "internal field".

In the simulations, five physical fields are mainly evaluated: pressure (p) (in incompressible flow simulations, OpenFOAM, uses a specific pressure, p/ρ [19]), velocity (U), turbulent kinematic viscosity (ν_t), turbulent kinetic energy (k) and specific turbulence dispersion velocity (ω). Table 1 shows the initial conditions for each physical field at each boundary condition and for the internal field.

The values of ω and ν_t at the input, output and internal field are specified by default to work together with the wall function implemented on the profile wall [19,26]. The initial values of k are determined by the following equation [19]:

$${\mathrm{k}}_0=1.5{\left(\mathrm{I}‖{\mathrm{U}}_{\mathrm{\infty }}‖\right)}^2,$$

(8)

where I is the intensity of the turbulence and U_∞ is the velocity of the free flow.

The indicated solver for this type of simulations is the SIMPLE, in the version 11 of OpenFOAM this solver is called incompressibleFluid [19]. The recommended schemes for working with the SIMPLE solver and the k-ω SST turbulence model are as shown in Table 2.

2.3.3. Encoding of the Data

When training neural networks, it is often useful to make sure that the data is normalized at all stages of the network. Normalization helps to stabilize and accelerate the training of the network by means of a downward gradient. If the data has been scaled badly, the loss function value may become an undefined value (NaN) and the network parameters may diverge during training. Common ways to normalize data include rescaling data so that its interval is [0, 1] or that it has a mean of zero and a standard deviation of one [27,28].

For normalization of the input vectors, which are the design parameters of the airfoils, the MinMax scaler is used [29], where the transformation is given by:

$${\mathsf{\xi }}_{\mathrm{j},\mathrm{s}\mathrm{t}\mathrm{d}}={\displaystyle \frac{{\mathsf{\xi }}_{\mathrm{j}}-{\mathsf{\xi }}_{\mathrm{j},\mathrm{m}\mathrm{i}\mathrm{n}}}{{\mathsf{\xi }}_{\mathrm{j},\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathsf{\xi }}_{\mathrm{j},\mathrm{m}\mathrm{i}\mathrm{n}}}},$$

(9)

$${\mathsf{\xi }}_{\mathrm{j},\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}\mathrm{d}}={\mathsf{\xi }}_{\mathrm{j},\mathrm{s}\mathrm{t}\mathrm{d}}\left({\mathrm{d}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{d}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+{\mathrm{d}}_{\mathrm{m}\mathrm{i}\mathrm{n}},$$

(10)

where ξ_j,min and ξ_j,max are the lower and upper limits of the calculated interval of the parameter ξ_j; d_min and d_max are new design intervals, 0 and 1, respectively.

For the case of the output data, the aerodynamic coefficients are represented in graphs, that is, images, which can be represented in matrices (three matrices representing the RGB layers), and whose components (pixels) have values between 0 and 1 (the original ranges are between 0 and 255, but the Python graphers also detect the range between 0 and 1 as valid) [30]. The graphs of interest are: endurance parameter vs angle of attack, pitch moment coefficient vs angle of attack and lift coefficient vs angle of attack. Each graph was normalized in images of 256 x 256 pixels, the axes and labels are eliminated, only the curve representing the values is left, the pixels that make up the curve have values of 1 (white color), while the pixels that make up the background have values of 0 (black color). Finally, the three images become layers of an RGB-encoded image (see Figure 3).

2.4. Design of the Neural Network for the Prediction of Aerodynamic Coefficients of Airfoils

The method for the design of the neural network is depicted in Figure 4, this method can be divided into the following three parts. 1 - the Principal Component Analysis (PCA) is used to perform the reduction of representative features of the output images and to obtain an approximation of the optimal dimension of the features involved in a Variational Autoencoder (VAE). 2 - VAE is used to extract representative features from the output images by unsupervised reconstruction-based learning. 3 - an Multi-layer Perceptron (MLP) network is developed to construct a non-linear mapping of the shapes of the aerodynamic surfaces to the features extracted by the VAE. The graphs of the aerodynamic coefficients of the profiles can be predicted by the composite network connecting the MLP and the decoder of the VAE [31].

2.4.1. Approximation of the Representative Features of the Output Images

The main goal of PCA is to find the principal component (PC) space, which is used to transform information from the original higher-dimensional space to the lower-dimensional space. The PC space consists of k components that are orthonormal, uncorrelated and represent the directions of the maximum variance of the data. The first principal component of the PCA space represents the direction of the maximum variance of the data, the second component has the second largest variance, and so on [32]. For image compression, the PCA method is used together with the information loss method (LOSSY). In the LOSSY method, some of the redundant information of the original file is lost, that is, when the file is restored, only the most important information is restored [29,32].

2.4.2. VAE Configuration

A VAE is a model consisting of an encoder and a decoder (see Figure 5). The encoder, also known as the output model or the recognition model, approximates the posterior distribution of the hidden variables of the decoder. The decoder is a deep hidden variable model and is a variation of the generative model. Both components, the encoder and the decoder, are directional graphic models that are totally or partially parameterized by deep neural networks [33].

In the case of an automatic encoder (AE), with a given input of ζ, the encoder converts ζ into a hidden representation of z that has a smaller dimension compared to ζ. The decoder then uses z to reconstruct ζ'. The recovery loss function L is calculated by comparing the differences between ζ and ζ'. VAEs are probabilistic modifications of AEs that combine Bayesian variational inference and deep learning. The encoder serves as an inference model that attempts to infer the latent variables z from the input data ζ through a learned distribution q(z|ζ); the decoder serves as a generative model that generates similar input data ζ‘ under the given latent vector z by sampling from learned distribution p(ζ|z). Assuming that the posterior q(z|ζ) of the m latent variables is a diagonal Gaussian with mean μ and standard deviation σ, and the prior p(ζ|z) is a normal Gaussian with μ=0 and σ=1 [31,33]. The Kullback-Leiber divergence is used to reduce the gap between the prior and posterior distributions, i.e.,

$$\mathrm{K}\mathrm{L}\left(\mathrm{q}\right(\mathrm{z}\left|\mathsf{\zeta }\right)\left|\right|\mathrm{p}\left(\mathsf{\zeta }\right|\left|\mathrm{z}\right))={\displaystyle \frac{1}{2}}{\sum }_{\mathrm{i}=1}^{\mathrm{m}}({\mathrm{e}}^{{\mathsf{\sigma }}_{\mathrm{i}}}-({\mathsf{\sigma }}_{\mathrm{i}}+1)+{\mathsf{\mu }}_{\mathrm{i}}^2),$$

(11)

where the vector z can be assembled as follows:

$$\mathbf{z}=\mathsf{\mu }+\mathsf{\sigma }\odot \mathsf{\epsilon },$$

(12)

where ε is a random vector of R(0, 1) and ⊙ denotes component multiplication.

The binary cross entropy recovery loss function BCE(ζ', ζ) is used to minimize the differences between the generated image and the original image, and is expressed as follows [30]:

$$\mathrm{B}\mathrm{C}\mathrm{E}\left(\mathsf{\zeta }\mathrm{\text{'}},\mathsf{\zeta }\right)=-\mathsf{\zeta }\log \left(\mathsf{\zeta }\mathbf{\text{'}}\right)-(1-\mathsf{\zeta })\log (1-\mathsf{\zeta }\mathrm{\text{'}}).$$

(13)

Therefore, the loss function to be minimized has the form of:

$${\mathrm{L}}_{\mathrm{V}\mathrm{A}\mathrm{E}}=\mathrm{B}\mathrm{C}\mathrm{E}+\mathrm{K}\mathrm{L}.$$

(14)

To minimize the loss function, network parameters are adjusted through back-forward propagation in every optimization iteration.

The proposed configuration for the encoder is composed of 6 convolutional layers and 6 MaxPooling layers (see Figure 6). After several blocks of convolutional residuals, the mean and standard deviation are extracted for each input image, in order to apply (12) and determine the z vector. The first 5 convolutional layers use the rectified linear unit (ReLU) activation function and the sixth layer uses the sigmoid activation function.

The architecture of the decoder is similar to that of the encoder, as shown in Figure 7, with the only difference that the convolutional layers are replaced by deconvolutive layers (or transposed convolution layer) and the layers of the MaxPooling operation are replaced by UpSamplin layers. The first five deconvolutive layers use ReLU activation function, and the last layer uses the sigmoid activation function.

To determine the performance in image reconstruction by the VAE, the mean absolute error (MAE) metric is used. The MAE describes the difference between the output data and the original input data [31]:

$$\mathrm{M}\mathrm{A}\mathrm{E}={\displaystyle \frac{1}{{\mathrm{n}}_{\mathrm{P}}}}{\sum }_{\mathrm{p}=1}^{{\mathrm{n}}_{\mathrm{P}}}\left|{\mathsf{\zeta }\mathrm{\text{'}}}_{\mathrm{p}}-{\mathsf{\zeta }}_{\mathrm{p}}\right|.$$

(15)

where n_P is the number of pixels of the images; ζ’_p is the value of the pixel p of the reconstructed image; ζ_p is the value of the pixel p of the real image.

Having the reconstructed images of the graphs of the aerodynamic coefficients, it is required to evaluate how efficient these images are to perform the readings of each coefficient in comparison to the original graphs. For this case, the use of the determination factor (R²) as a metric is proposed. R² measures the fraction of variance of a dependent variable explained by independent variables. The values of R² vary between 0 and 1, the closer the values of R² to 1, indicates a better fit of the model [34,35,36].

2.4.3. MLP Configuration

The MLP is in charge of performing a non-linear mapping of the vector that indicates the characteristics associated with the graphs derived from the VAE and the normalized parameters of the airfoil parameterization method. The loss function to be minimized is the mean square error (MSE) between the difference between the z vectors provided by the encoder and the z' vectors predicted by the MLP [31]:

$${\mathrm{L}}_{\mathrm{M}\mathrm{L}\mathrm{P}}=\mathrm{M}\mathrm{S}\mathrm{E}\left(\mathbf{z}\mathrm{\text{'}},\mathbf{z}\right)={\displaystyle \frac{1}{\mathrm{m}}}{\sum }_{\mathrm{i}=1}^{\mathrm{m}}\left({\displaystyle \frac{1}{2}}{‖\mathrm{z}{\mathrm{\text{'}}}_{\mathrm{i}}-{\mathrm{z}}_{\mathrm{i}}‖}^2\right).$$

(16)

The proposed MLP architecture has an input layer with j neurons (corresponding to the number of parameters of the airfoil parameterization method), an output layer with m neurons (corresponding to the size of the z vector), and three hidden layers. To determine the number of neurons and the activation function of each hidden layer an optimization algorithm is used, where the optimization task is defined as:

min f₁(η) = MSE(z’, z),

such that η ∊ Ω.

(17)

In this case η represents the vector formed by the parameters n₁, n₂, n₃ representing the number of neurons in each of the hidden layers respectively, and by a₁, a₂, a₃, representing the activation function in each hidden layer. The optimization algorithm chosen to solve this optimization task is the Integer Encoding Differential Evolution (IEDE) [34]. IEDE is a version of the standard DE algorithm where only the mutation operator is modified, so that it can operate with integer type values [37].

$${\mathrm{v}}_{\mathrm{i},\mathrm{j}}=\left\{\begin{array}{c}\begin{array}{cc}{\mathsf{\eta }}_{\mathrm{r}1,\mathrm{j}}+\mathrm{i}\mathrm{n}\mathrm{t}\uparrow \left({\mathsf{\delta }}_{\mathrm{i},\mathrm{j}}\right)& \mathrm{i}\mathrm{f}\mathrm{r}<{\mathsf{\delta }}_{\mathrm{i},\mathrm{j}}-\mathrm{i}\mathrm{n}\mathrm{t}\uparrow \left({\mathsf{\delta }}_{\mathrm{i},\mathrm{j}}\right)\end{array}\\ \begin{array}{cc}{\mathsf{\eta }}_{\mathrm{r}1,\mathrm{j}}+\mathrm{i}\mathrm{n}\mathrm{t}\downarrow \left({\mathsf{\delta }}_{\mathrm{i},\mathrm{j}}\right)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\end{array}\right.,$$

(18)

$${\mathsf{\delta }}_{\mathrm{i}}=\mathrm{F}({\mathsf{\eta }}_{\mathrm{r}2}-{\mathsf{\eta }}_{\mathrm{r}3}).$$

(19)

where F is the mutation factor; r is a randomly selected number based on the normal distribution; η_r1, η_r2, η_r3 are three different vectors randomly selected in the current population; int↑(), int↓() are functions that round a number to the nearest upper and lower integer, respectively.

Having defined the architecture of the MLP, it is connected with the decoder of the VAE in order to have a neural network that is capable of predicting the graphs of the aerodynamic coefficients of a profile from a vector of design parameters (see Figure 8).

2.5. Optimization Algorithms

2.5.1. Single-Objective Optimization Algorithm

To solve single-objective optimization tasks, we opted to use the Success-History based Adaptive Differential Evolution (SHADE) algorithm [38]. This variant of the DE algorithm has the advantage that it adapts its evolutionary operators as it finds better values of the mutation and crossover factors, that is, it evolves as it searches for the optimal solution to the ongoing problem. Another advantage of the SHADE algorithm is that it can be adapted to the use of population size reduction (PSR) methods [38,39,40]. These methods help to solve part of the disadvantage of evolutionary algorithms, which is to reduce the number of objective functions to be evaluated in each generation, without having to lose performance in the search for an optimal global value. The variant of the SHADE algorithm to be used in this work incorporates the - Continuous Adaptive Population Reduction (CAPR) [41] method. The performance of the CAPR-SHADE algorithm has been tested in aerodynamic optimization tasks [42,43].

In CAPR-SHADE, the mutation factor F ∊ [0, 1] controls the magnitude of the differential mutation operator and CR ∊ [0, 1] is the crossover rate. This algorithm maintains a historical memory (M_F and M_CR) with H entries for both control parameters F and CR. In the beginning, the contents of M_F,k, M_CR,k (k = 1, ..., H) are initialized to 0.5. In each generation g, the control parameters CR_i and F_i used by each individual ri are generated by randomly selecting the index ri from [1; H], and then applying the formulas below:

$${\mathrm{C}\mathrm{R}}_{\mathrm{i}}=\left\{\begin{array}{cc}0& \mathrm{i}\mathrm{f}{\mathrm{M}}_{\mathrm{C}\mathrm{R},\mathrm{r}\mathrm{i}}=-1\\ {\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{n}}_{\mathrm{i}}({\mathrm{M}}_{\mathrm{C}\mathrm{R},\mathrm{r}\mathrm{i}},0.1)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.,$$

(20)

$${\mathrm{F}}_{\mathrm{i}}={\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{c}}_{\mathrm{i}}({\mathrm{M}}_{\mathrm{F},\mathrm{r}\mathrm{i}},0.1).$$

(21)

Here, randn_i(M, 0.1), randc_i(M, 0.1) are randomly selected values of normal and Cauchy distributions respectively, with mean of M (M_CR or M_F) and variance 0.1. In equation (20), in case a value for CR_i outside [0, 1] is generated, it is replaced by the limit value (0 or 1) closest to the generated value. When F_i > 1, F_i is rounded to 1, and when F_i ≤ 0, (21) is applied repeatedly until a valid value is generated. In (20), if MCR,ri has been assigned the -1 value, CR_i is set to 0.

The mutation strategy used by CAPR-SHADE is current-to-pbest/1 which is a variant of the current-to-best/1 strategy in which the operator is adjustable using a parameter p:

$${\mathbf{v}}_{\mathrm{i},\mathrm{g}}={\mathsf{\xi }}_{\mathrm{i},\mathrm{g}}+{\mathrm{F}}_{\mathrm{i}}\left({\mathsf{\xi }}_{\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t},\mathrm{g}}-{\mathsf{\xi }}_{\mathrm{i},\mathrm{g}}\right)+{\mathrm{F}}_{\mathrm{i}}({\mathsf{\xi }}_{\mathrm{r}1,\mathrm{g}}-{\mathsf{\xi }}_{\mathrm{r}2,\mathrm{g}}).$$

(22)

The individual ξ_pbest,g is randomly selected from the upper p·NP members in the g-th generation (p ∊ [0, 1]). F_i is the parameter F used by individual ξ_i. The elitism of the current-to pbest/1 depends on the control parameter p, which provides a balance between exploitation and exploration (small values of p indicate greater discrimination in the elite set). In case any component j of the mutated vectors v_i is out of the decision range, the value will be adjusted to the nearest extreme value.

The crossover operator uses a binary strategy, similar to the one used in DE standard:

$${\mathrm{u}}_{\mathrm{j},\mathrm{i},\mathrm{g}}=\left\{\begin{array}{cc}{\mathrm{v}}_{\mathrm{j},\mathrm{i},\mathrm{g}}& \mathrm{i}\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\left(\mathrm{0,1}\right)\le {\mathrm{C}\mathrm{R}}_{\mathrm{i}}\\ {\mathsf{\xi }}_{\mathrm{j},\mathrm{i},\mathrm{g}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right..$$

(23)

The selection operator compares each individual ξ_i,g with its corresponding trial vector u_i,g, where the best one passes to the next generation:

$${\mathsf{\xi }}_{\mathrm{i},\mathrm{g}+1}=\left\{\begin{array}{cc}{\mathbf{u}}_{\mathrm{i},\mathrm{g}}& \mathrm{i}\mathrm{f}\mathrm{f}\left({\mathbf{u}}_{\mathrm{i},\mathrm{g}}\right)\le \mathrm{f}\left({\mathsf{\xi }}_{\mathrm{i},\mathrm{g}}\right)\\ {\mathsf{\xi }}_{\mathrm{i},\mathrm{g}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right..$$

(24)

In each generation, in (24), CR_i and F_i values that succeed in generating a trial vector u_i,g which is better than the parent individual ξ_i,g are recorded as S_CR, S_F , and at the end of the generation, the M_CR and M_F memories content are updated using Algorithm 1.

Algorithm 1. Memory update algorithm in SHADE [39].
	Input: SCR, SF, MCR,k,g, MF,k,g, k , H
	Output: MCR,k,g+1, MF,k,g+1, k
1	if SCR ≠ ∅ and SF ≠ ∅ then
2	if MCR,k,g = -1 or max(SCR) = 0 then
3	MCR,k,g+1 = -1;
4	else
5	MCR,k,g+1 = meanWL(SCR);
6	MF,k,g+1 = meanWL(SF);
7	k++;
8	if k > H then
9	k = 1;
10	else
11	MCR,k,g+1 = MCR,k,g;
12	MF,k,g+1 = MF,k,g;
13	return Output

In Algorithm 1, the index k determines the position in the memory to be updated. In generation g, the element k in the memory is updated. At the beginning of the search, k is initialized to 1. k is incremented every time a new item is inserted in the history. If k > H, k is set to 1. In Algorithm 1 it is observed that when all individuals of generation g do not generate a trial vector that is better than the parent, i.e. S_CR = S_F = ∅, the memory is not updated. The weighted mean of Lehmer mean_WL(S) is calculated as follows:

$${\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}_{\mathrm{W}\mathrm{L}}\left(\mathrm{S}\right)={\displaystyle \frac{\sum _{\mathrm{m}=1}^{\left|\mathrm{S}\right|}{\mathrm{w}}_{\mathrm{m}}{\mathrm{S}}_{\mathrm{m}}^2}{\sum _{\mathrm{m}=1}^{\left|\mathrm{S}\right|}{\mathrm{w}}_{\mathrm{m}}{\mathrm{S}}_{\mathrm{m}}}},$$

(25)

$${\mathrm{w}}_{\mathrm{m}}={\displaystyle \frac{∆{\mathrm{f}}_{\mathrm{m}}}{\sum _{\mathrm{l}=1}^{\left|\mathrm{S}\right|}{∆\mathrm{f}}_{\mathrm{l}}}},$$

(26)

$${∆\mathrm{f}}_{\mathrm{m}}=\left|\mathrm{f}\left({\mathbf{u}}_{\mathrm{m},\mathrm{g}}\right)-\mathrm{f}\left({\mathsf{\xi }}_{\mathrm{m},\mathrm{g}}\right)\right|,$$

(27)

The amount of improvement ∆f_m is used to influence the adaptation parameter S.

The CAPR method gradually reduces the size of the population according to the change in the gradient of the fitness value [41]:

$${\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}=\left\{\begin{array}{cc}\sqrt[\mathsf{\gamma }]{{\displaystyle \frac{{\mathsf{\Delta }}_{\mathrm{g}}}{{\mathsf{\Delta }}_{\mathrm{g}-1}}}}& \mathrm{i}\mathrm{f}0<{\displaystyle \frac{{\mathsf{\Delta }}_{\mathrm{g}}}{{\mathsf{\Delta }}_{\mathrm{g}-1}}}<1\\ {\mathrm{N}\mathrm{P}}_{\mathrm{g}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.,$$

(28)

$${\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}=\left\{\begin{array}{cc}{\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}& \mathrm{i}\mathrm{f}{\mathrm{N}\mathrm{P}}_{\mathrm{g}+1}>{\mathrm{N}\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\\ {\mathrm{N}\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.,$$

(29)

where

$$\begin{array}{cc}{\mathsf{\Delta }}_{\mathrm{g}}={\displaystyle \frac{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathsf{\xi }}_{\mathrm{g}}\right)-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathsf{\xi }}_{\mathrm{g}-1}\right)}{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathsf{\xi }}_{\mathrm{g}}\right)}},& {\mathsf{\Delta }}_{\mathrm{g}-1}={\displaystyle \frac{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathsf{\xi }}_{\mathrm{g}-1}\right)-{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathsf{\xi }}_{\mathrm{g}-2}\right)}{{\mathrm{f}}_{\mathrm{a}\mathrm{v}\mathrm{g}}\left({\mathsf{\xi }}_{\mathrm{g}-1}\right)}}\end{array},$$

(30)

The above equations are used from the third generation, where f_avg indicates the averaging of the values of the objective functions of all individuals in each generation.

The CAPR-SHADE algorithm as well as the standard DE algorithm are created to be used in unconstrained optimization tasks. A penalty function is used to convert a constrained optimization task to a unconstrained one. The penalty function is expressed as follows [44]:

$$\mathrm{L}\left(\mathsf{\xi }\right)=\left\{\begin{array}{cc}\mathrm{f}\left(\mathsf{\xi }\right)& \mathrm{i}\mathrm{f}\mathsf{\psi }\left(\mathsf{\xi }\right)=0\\ \mathrm{R}\mathsf{\psi }\left(\mathsf{\xi }\right)+{\mathrm{U}}^{\mathrm{*}}& \mathrm{i}\mathrm{f}\mathsf{\psi }\left(\mathsf{\xi }\right)>0\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{f}\left(\mathsf{\xi }\right)\le {\mathrm{U}}^{\mathrm{*}}\\ \mathrm{R}\mathsf{\psi }\left(\mathsf{\xi }\right)+\mathrm{f}\left(\mathsf{\xi }\right)& \mathrm{i}\mathrm{f}\mathsf{\psi }\left(\mathsf{\xi }\right)>0\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{f}\left(\mathsf{\xi }\right)>{\mathrm{U}}^{\mathrm{*}}\end{array}\right.,$$

(31)

where

$$\mathsf{\psi }\left(\mathsf{\xi }\right)={\sum }_{\mathrm{k}=1}^{\mathrm{K}}\mathrm{m}\mathrm{a}\mathrm{x}\left\{0,{\mathrm{g}}_{\mathrm{k}}\left(\mathsf{\xi }\right)\right\},$$

(32)

where U* is the minimum value of the values of the objective functions of those vectors that satisfied the constraints. Initially, a large value is assigned to U*, but the value must be updated with the current best-known function value at feasible points. Therefore, the initial U* is not updated until a feasible solution has been found. R is a penalty parameter that allows to regularize the obtained values of ψ(ξ) and U*. The constraints expressed as equalities (h(ξ)) have to be reformulated as inequalities (g(ξ)).

Using the penalty function, optimization tasks as expressed in (1) are expressed as:

min L(ξ),

such that ξ ∊ Ω.

(33)

2.5.2. Multi-Objective Optimization Algorithm

The Non-Dominated Classification Genetic Algorithm (NSGA-II) is a powerful decision space exploration engine based on Genetic Algorithm (GA) to solve multi-objective optimization problems [45]. The philosophy of NSGA-II is based on four fundamental principles, which are: non-dominated classification, elite preservation operator, crowding distance and selection operator [45,46].

The non-dominated classification process begins with the assignment of the first rank to the non-dominated members (based on Pareto dominance) of the initial population. These first classified members are then placed on the first front and removed from the initial population. After that, the non-dominant classification procedure is performed on the remaining members of the population. In addition, the non-dominated members of the remaining population are assigned the second rank and placed on the second front. This process continues until all the members of the population are placed on different fronts according to their ranks, as shown in Figure 9(a). Algorithm 2 shows the process of the non-dominated classification.

Algorithm 2. Non-dominated classification [45].
	Input: Pt, f1(Pt), f2(Pt)
	Output: Fi
1	for p in Pt do;
2	Sp = ∅;
3	np = 0;
4	for q in Pt do
5	if p ≺ q then
6	q→Sp;
7	else
8	np++;
9	if np == 0 then
10	prank = 1;
11	p→F1;
12	i = 1
13	while Fi ≠ ∅ do
14	Q = ∅;
15	for p in Fi do
16	for q in Sp do
17	nq = nq -1;
18	if nq == 0 then
19	qrank = i+1;
20	q→Q;
21	i++;
22	Fi = Q;
23	return Output

The elite preservation strategy retains the elite solutions of a population by transferring them directly to the next generation. In other words, the non-dominated solutions found in each generation pass to the next generations until some new solutions dominate them.

The crowding distance (CD) is calculated to estimate the density of the solutions surrounding a particular solution. It is the average distance of two solutions on either side of the solution along each of the objectives (see Figure 9(b)). The CD of the k-th individual is defined as the average distance of two closest solutions on each side:

$${\mathrm{C}\mathrm{D}}_{\mathrm{k}}={\sum }_{\mathrm{l}=1}^2{\displaystyle \frac{{\mathrm{f}}_{\mathrm{l},\mathrm{k}+1}-{\mathrm{f}}_{\mathrm{l},\mathrm{k}-1}}{{\mathrm{f}}_{\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{f}}_{\mathrm{l},\mathrm{m}\mathrm{i}\mathrm{n}}}},$$

(34)

where f_l,max and f_l,min are the maximum and minimum values, respectively, of the l-th objective function among all vectors.

Algorithm 3 shows the process for determining the crowding distance on each front.

Algorithm 3. Crowding distance [45].
	Input: F
	Output: CD
1	nk = |F|;
2	CD = ∅;
3	CD1 = CDnk = ∞;
4	for k = 2 to nk-1 do
5	Apply (34);
6	CDk → CD;
7	return Output

The encoding of the individuals to be used in NSGA-II is floating-type. The genetic operators selected to operate with floating-type coding are: binary tournament selection, simulated binary crossover (SBX) operator and polynomial mutation.

Algorithm 4 describes the binary tournament selection.

Algorithm 4. Binary tournament selection.
	Input: Pt, pt
	Output: p
1	Select two vectors at random from Pt, ξr1≠ξr2;
2	Select a random number r; from a normal distribution;
3	if r < pt then
4	if ξr1 ≺ ξr2 then
5	p = ξr1;
6	else
7	p = ξr2;
8	else
9	if ξr1 ≺ ξr2 then
10	p = ξr2;
11	else
12	p = ξr1;
13	return Output

The SBX operator simulates the principle of operation of the single-point crossing operator on binary chains [47].

$$\begin{array}{c}{\mathrm{c}}_{1,\mathrm{j}}={\displaystyle \frac{1}{2}}\left(\left(1-{\mathsf{\beta }}_{\mathrm{j}}\right){\mathrm{p}}_{1,\mathrm{j}}+\left(1+{\mathsf{\beta }}_{\mathrm{j}}\right){\mathrm{p}}_{2,\mathrm{j}}\right)\\ {\mathrm{c}}_{2,\mathrm{j}}={\displaystyle \frac{1}{2}}\left(\left(1+{\mathsf{\beta }}_{\mathrm{j}}\right){\mathrm{p}}_{1,\mathrm{j}}+\left(1-{\mathsf{\beta }}_{\mathrm{j}}\right){\mathrm{p}}_{2,\mathrm{j}}\right)\end{array},$$

(35)

where c₁ and c₂ are the child vectors of the parent vectors p₁ and p₂; the subscript j represents the components of each vector; and β is obtained from:

$${\mathsf{\beta }}_{\mathrm{j}}\left(\mathrm{u}\right)=\left\{\begin{array}{cc}{\left(2\mathrm{u}\right)}^{{\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{c}}+1}}}& \mathrm{i}\mathrm{f}\mathrm{u}<0.5\\ {\displaystyle \frac{1}{{\left(2\right(1-\mathrm{u}\left)\right)}^{\frac{1}{{\mathsf{\eta }}_{\mathrm{c}}+1}}}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.,$$

(36)

where u is a randomly selected number from a normal distribution from 0 to 1; η_c is the distribution index for the crossover.

The mutation operator uses the polynomial mutation strategy [48]:

$${\mathrm{v}}_{\mathrm{i},\mathrm{j}}={\mathrm{c}}_{\mathrm{i},\mathrm{j}}+({\mathsf{\Omega }}_{\mathrm{j}}^{\mathrm{u}}-{\mathsf{\Omega }}_{\mathrm{j}}^{\mathrm{l}}){\mathsf{\delta }}_{\mathrm{j}},$$

(37)

where v_i is the mutated vector of the child vector c_i; Ω^l_j and Ω^u_j are the minimum and maximum values of the design parameter j, respectively; δ_j is a small deviation that is calculated from:

$${\mathsf{\delta }}_{\mathrm{j}}\left(\mathrm{u}\right)=\left\{\begin{array}{cc}{\left(2\mathrm{u}\right)}^{{\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{m}}+1}}}& \mathrm{i}\mathrm{f}\mathrm{u}<0.5\\ {\displaystyle \frac{1}{1-{\left(2\right(1-\mathrm{u}\left)\right)}^{\frac{1}{{\mathsf{\eta }}_{\mathrm{m}}+1}}}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.,$$

(38)

where η_m is the index of the mutation distribution.

The NSGA-II procedure is generally shown in Figure 10.

3. Results and Discussion

The calculations that were performed in each of the stages were performed with a computer with the following characteristics:

• hardware model - Gigabyte Technology Co., Ltd. B760 GAMING X;
• processor-13th Gen Intel® CoreTM i5-13400f x 12;
• RAM memory - 64 GB;
• GPU - NVIDIA GeForce RTX 4090;
• OS - Ubuntu 22.04.4 LTS.

3.1. Selection of the Method of Airfoils Parameterization

In the airfoil reconstruction tests using the CST and Bezier-PARSEC methods, 322 airfoils extracted from the UIUC database were used. The DE variant “current-to-best/1” was selected as the optimization algorithm [49]. In the optimization tests a value of F = 0.85 was used, the crossover operation was omitted, an initial population of NP = 10D was taken (D is the size of the vector of design parameters) and a maximum of 200 generations were evaluated for each airfoil [16]. The decision spaces for each method are defined in Table 3 and Table 4.

Figure 11 shows a histogram where the results obtained from the reconstruction tests of existing airfoils can be visualized.

It can clearly be seen that the CST method was able to reconstruct a greater number of profiles with the smallest possible geometric deviation. Figure 12 shows 4 examples of profiles reconstructed with both methods, comparing them with the original profile. The graphs show the local deviations in the upper curve and in the lower curve of the profile. In all cases it is visualized that the CST Method achieves a better reconstruction of the profile.

Based on the results obtained, the CST method was selected to demonstrate better flexibility than the CST method, the only disadvantage that this model has is the lack of intuitiveness, since the parameters used in this method lack an aerodynamic meaning of the profile.

3.2. New Geometries of the Airfoils Using a GAN

The architecture of the GAN that was used for the generation of new profiles is shown in Figure 13 and its specifications are shown in Table 5 and Table 6. Because the profiles to be generated can be encoded as ξ_real vectors (each vector consisting of the profile parameters CST), then the generator and the discriminator can be MLP. To achieve a better training of the GAN the design intervals of the CST parameters (see Table 3) were normalized between the values of -1 to 1.

1200 epochs were required for training. Due to the fact that at this point of the research there are only 322 profiles, a convergence could not be guaranteed, for this reason the results were monitored every 200 epochs. In Figure 14 it can be seen that between 800 and 1000 or at the end of the 1200 evaluated epochs, the values of the loss function of the generator and the discriminator tend to be the same, in these epochs it is where it is required to save the weights of the neural network.

Figure 15 shows geometries that are clearly similar to the airfoils used for the GAN training, of course there are airfoils that have no resemblance to the real airfoils, but still retain the geometric characteristics that can identify them as an airfoil.

When analyzing 678 new airfoils, it was found that the values of the maximum profile thickness (y_t,max) are in the range from 0.066 to 0.17, and the values of the maximum curvature of the profile (y_c,max) are in the range from 0.0117 to 0.065. The ranges of the above parameters, for the real profiles, are y_t,max from 0.06 to 0.19 and y_c,max from 0.0055 to 0.065.

3.3. Validation of CFD Simulations

To support the configuration of the simulation in OpenFOAM and the mesh generated in GMSH, the aerodynamic coefficients of the NACA 0012 profile with different Reynolds and Mach numbers were obtained, comparing them with the coefficients obtained experimentally in the laboratories of the National Aeronautics and Space Administration of the United States [50]. At the suggestion of OpenFOAM, the comparison between the experimental results and those obtained by CFD, the data obtained in the tests with fixed transition of the boundary layer should be used, the tests with grain size of 80, 120 and/or 180 are suggested [51]. The characteristics of the tests and the coefficients obtained at different Re and M are shown in detail in Table 7 and Figure 16.

The results obtained using OpenFOAM are suitable for generating the aerodynamic coefficient database of the profiles. It is important to mention that in this work the CFD simulations will be limited until the maximum value of the endurance parameter is found, this because for the moment the evaluation of the profile loss region will be avoided due to its complexity.

3.4. Approximation of the Representative Features of the Output Images

Figure 17, Figure 18 and Figure 19 show 4 tests with different numbers of PC (5, 10, 15, 20) applied to the three layers of one of the output images. From 10 PC, the shape of the curves can be recreated, only having variation in the color saturation in the pixels, which is not of great importance for reading values.

3.5. Evaluation of the VAE

Considering a database of 1000 profiles, 90% of the data was selected for the VAE training and the remaining 10% was used for testing. In the PCA tests it was determined that at least 10 main components are required to be able to reconstruct the images. A VEA being a more robust model requires fewer components, for this reason 5 different sizes of the z vector will be tested, 5, 6, 7, 8 and 9. The results of the reconstruction of the images in both stages are shown in Figure 20 and Figure 21 .

In the testing stage, the generalization capacity of a neural network is defined. The average MAE values of the test data for each z size are shown in Table 8. The best values are highlighted in green. The z vectors of size 7 and 8 come to generalize in a similar way, in this case the vector of smaller size is chosen.

An example of graph reconstruction using the VAE with a 7-component z vector is shown in Figure 22. Figure 23 shows a linear fit analysis between the read values of the aerodynamic coefficients obtained from the real graphs (those obtained from the CFD simulations) and the graphs reconstructed by the VAE. The values of R² are indicated for each case, the ideal fit line is indicated (when the predicted values (P) are equal to the target values (T)) and the line of the best linear fit (BLF) of the model.

3.6. Design of the MLP and Performance Evaluation of AZTLI-NN

To solve the optimization task shown in (17), the decision space Ω shown in Table 9 is considered.

In the optimization task the output layer has a sigmoidal activation function, the optimizer is ADAM, out of 90% of the data assigned for training 10% is assigned for validation during training, and 500 epochs were considered. The IEDE algorithm has a F = 0.85, the population size of 30, and 30 generations were evaluated. The optimal architecture is shown in Figure 24. In the input layer there are 14 neurons, corresponding to the parameters of the CST method, the first hidden layer has 161 neurons, the second with 320 and the third layer with 511, finally the output layer has 7 neurons corresponding to the size of the z vector of the VAE.

The results of R² of the coefficient reading with the graphs predicted by the AZTLI-NN assembled network are shown in Figure 25.

Statistically, it can be considered that the neural network achieves a good generalization, but at present, many of the neural networks dedicated to the prediction of aerodynamic coefficients achieve values of R²~0.99. To analyze whether it is possible to increase the accuracy of the neural network, the number of data was increased to 1500 with the help of the GAN. As the number of training data increased, the predictions improved, the R2 value for c_l^1.5/c_d improved by 2.2%, for the c_m improved by 3.2%, and for the c_l improved by 1.2%. In all three cases the R² values were closer to 0.99, and even for the prediction of c_l a value greater than 0.99 was obtained (see Figure 26).

From the profiles used for testing, information was sought about experimental tests with similar characteristics to which the neural network was trained. Information on experimental tests of the FX 66-S-161 profile was found in [52]. The experimental tests were developed in a laminar wind tunnel with Re = 1.5x10⁶ and M = 0.25. The aerodynamic coefficients corresponding to the FX 66-S-161 profile are shown in Figure 27.

3.7. Solving the First Optimization Case

The first proposed optimization case is described in detail as follows: during the conceptual design process of an aircraft it is necessary to find an aerodynamic profile that has the following characteristics: design lift coefficient, c_l,d = 0.59; minimum maximum thickness of the airfoil, y_t,min = 11%; maximum permissible angle of attack, α_max = 4°. The flow conditions are: Re = 1.5x10⁶ and M = 0.15. The optimization algorithm used to solve this task is shown in Appendix C.

In the first instance, the convergence of the algorithm was evaluated based on the number of generations, for this an initial population size NP₀ = 10|ξ|, a minimum population NP_min = 4 (recommended in [41]) was proposed. 500 generations were evaluated. The test was performed 5 times to analyze the repeatability of the results. The results of the 5 tests performed are shown in Figure 28 .

As can be visualized in Figure 29 the algorithm requires less than 100 generations to converge. In the next series of tests, a value of 100 generations was maintained to be evaluated. In the following series of tests the size of the initial population was varied, with sizes of 10|ξ|, 20|ξ| and 50|ξ|. The results of this series of tests are shown in Table 10. The 15 airfoils obtained in the different tests are shown in Figure 29 .

The optimal airfoils that obtained the maximum and minimum value of y_tmax are shown in Figure 30 . The aerodynamic coefficients of both profiles obtained with AZTLI-NN and with OpenFOAM are compared in Table 11.

3.8. Solving the Second Optimization Case

The second proposed optimization case is described in detail as follows: during the conceptual design process of an aircraft wing, it is necessary to provide profile variants that provide high values of c_l^1.5/c_d and close to zero values of c_m; and a design lift coefficient (c_l,d) of 0.59 is considered. The flow conditions are: Re = 1.5x10⁶ and M = 0.15.

The default values for NSGA-II of the population size (NP = 100) and the number of generations (T = 250) were used [45]. This test was performed 5 times to analyze the repeatability of the algorithm shown in Appendix D.

Figure 31 shows how the Pareto front of each of the 5 tests evolves. After evaluating 250 generations, it can be observed that the Pareto boundaries tend to have the same shape, which implies that there is repeatability in the results.

The computation time of each test is approximately 120 s.

The best Pareto boundary obtained is shown in Figure 32 (In test 5, a boundary with the highest number of dominant individuals was obtained). In addition, the geometries of the airfoils that are part of the Pareto front are shown in Figure 32.

Figure 33 shows in more detail the comparison of the geometries of the airfoils shown in Figure 32 . It can be seen that the position of the maximum thickness (∼25% of the chord) and the position of the maximum camber (∼0.20 of the chord) of the airfoils are similar. The main difference between the airfoils is shown in the camber at the trailing edge of the airfoil. As |c_m|decreases (or c_l^1.5/c_d decreases) the camber decreases more rapidly at the trailing edge of the airfoil.

The graphs of the aerodynamic coefficients of each of the airfoils of the Pareto front are shown in Figure 34. Among the most remarkable things that can be observed in Figure 34 is the following: in the graph of c_l vs α it can be observed that all the airfoils have a linear behavior within the range of 0 to 10 degrees; the maximum values of c_l^1.5/c_d of the airfoils are between 8 and 9 degrees; from 0 to 7 degrees the c_m values of each airfoil do not show a significant change.

In the graphs of Figure 34 the design condition is indicated with a dot (c_l,d = 0.59), in Table 12 the results obtained under this design condition are shown in more detail, in addition it is compared with the aerodynamic coefficients obtained by OpenFOAM.

4. Conclusions

As part of this work, a methodology was developed to generate optimization algorithms, aimed at the aerodynamic optimization of airfoils under conditions of subsonic, incompressible and turbulent flow. In the development process, metaheuristic algorithms (specifically evolutionary algorithms) and different artificial neural network architectures were evaluated. The purpose of the algorithms developed is that they can be used in the early stages of aircraft design with high elongation and long flight duration. To achieve the established goal, the following research and developments were carried out:

• An artificial neural network architecture, AZTLI-NN, (composed of a multilayer perceptron and a variational autoencoder) was developed for aerodynamic response prediction. A new image-based coding was proposed for the output parameters of the neural network, the neural network has the ability to generate the diagrams of the aerodynamic coefficients as a function of the angle of attack. For the training of the neural network, a database of wing airfoils was generated with their respective aerodynamic coefficients (lift coefficient, drag coefficient, pitch moment coefficient) using computational aerodynamics. The procedure of the numerical simulations was validated with experimental cases. Airfoil parameterization methods were evaluated to determine which provided the best performance when reconstructing wing airfoils.
• Evolutionary algorithms of mono-objective and multi-objective optimization were adapted so that they could be used in conjunction with the AZTLI-NN network.
• The performance of each of the algorithms was evaluated. Several tests were carried out to evaluate the repeatability of the results and the consistency in the computation times. In both cases, repeatability of the results was obtained, and the computation times are suitable for the algorithms to be considered in early-stage design processes.

Therefore, the work carried out achieves the practical importance raised, since not only was it possible to meet the established objective, but also provided a background to achieve the development of neural network architectures and adapt evolutionary algorithms to be used in optimization tasks of airfoils aimed at different tasks. For example, the proposed methodology can be used in the design of aircraft operating in different flight regimes, in the design of propellers, in the design of wind turbines, etc.