*2.1. Genetic Algorithms (GAs)*

Genetic algorithms (GAs) proposed by Professor J. Holland [19] are among the evolutionary algorithms (EAs). GAs operate on the whole population with individuals, and their main operators include selection, crossover and mutation. For a particular problem, GAs define the search space as the solution space, and each feasible solution is encoded as a chromosome. Before the search starts, a set of chromosomes is usually randomly selected from the solution space to form the initial population. Next, the fitness value of each individual is evaluated according to the objective function, then the selection, crossover and mutation operators are applied sequentially to generate a new generation of populations. The process is repeated until the stopping criterion is reached.

### *2.2. Cultural Algorithm (CA)*

The two-layer evolutionary mechanism used by the cultural algorithm consists of two main evolutionary spaces at the micro and macro levels, namely the population space and the belief space [28], and the basic structure of the cultural algorithm is shown in Figure 1. The evolution on the micro level refers to the internal evolution of the population space that realizes the evolution of individuals, and the evolution on the macro level refers to the evolution of the belief space that realizes the extraction and updating of knowledge sources. The evolutions between these two spaces are independent of each other, but they are connected through communication protocols (influence and acceptance functions). Figure 2 describes the basic pseudo-code of the CA. The figure shows how the process is executed in each generation. Firstly, the objective function Obj() evaluates individuals in the population space, and the Acceptance() function selects the best individuals for updating the belief space knowledge source. After that, the Influence() function influences the evolution of the next generation of populations. More details on the knowledge sources used and how they affect the population of this proposed work are given in Section 3.

**Figure 1.** Framework of the cultural algorithm (CA).

**Figure 2.** Pseudo-code of cultural algorithm.

### **3. The Hybrid Evolutionary Optimization Method Coupling CA with GAs**

Cultural systems possess the ability to incorporate heterogeneous and diverse knowledge sources into their structures. As such, they are ideal frameworks within which to support hybrid amalgams of knowledge sources and population components [29]. In order to make full use of the advantages of CA and GAs, an efficient hybrid evolutionary optimization method coupling CA with GAs (HCGA) is proposed in this paper. The cultural framework of HCGA is shown in Figure 3, which includes population space, belief space and communication protocol, whose population space is modeled using GAs, and belief space includes situational knowledge, normative knowledge and historical knowledge. In addition, HCGA introduces population entropy and population variance to judge population diversity, and a knowledge-guided *t*-mutation operator is developed based on population diversity to balance the exploration and exploitation ability of the algorithm. In the remainder of this section, we describe each part of the HCGA in detail.

### *3.1. Population Space*

In fact, the population space can support any population-based evolutionary algorithm or swarm intelligence algorithm, which can also interact and run simultaneously with the belief space. The standard cultural algorithm has only a single mutation operator in the population space, making its global convergence and exploration capability insufficient. The GAs has a strong global search capability and high robustness, which can effectively explore the search space with the increasing population convenience and global exploration capability of the algorithm, and thus the population space is evolved using the GAs in this paper. A detailed description of the genetic algorithms is given in Section 2.2 and will not be repeated here.

### *3.2. Belief Space*

In this paper, according to the characteristics of genetic algorithms, combined with the manner of extracting and updating knowledge sources in the belief space, the knowledge sources are divided into situational knowledge, normative knowledge and historical knowledge. The manner of updating the knowledge sources in the belief space every *K* generations is adopted, so that the memory consumption brought by redundant information can be reduced. Different knowledge sources have different update strategies. Taking the maximization problem as an example, the update of the knowledge sources is described as follows:

(1) Situational knowledge. Situational knowledge was proposed by Chung in 1997 [30] to record the excellent individuals with a guiding role in the evolutionary process and is structured as follows:

$$ \tag{1}$$

where *E<sup>i</sup>* is the *i*th best individual, and in this paper the best individual is selected to update the situational knowledge, that is, *s* = 1. The process of updating situational knowledge is described as follows:

$$ = \begin{cases}  & \text{if } f(x\_b) > f(E(T)) \\  & \text{else} \end{cases} \tag{2}$$

where *x<sup>b</sup>* is the best individual in the Tth generation of the population space.

(2) Normative knowledge. Normative knowledge was also proposed by Chung [30] for limiting the search space and judging the feasibility of an individual. When an individual is outside the search space described by the normative knowledge, the normative knowledge will guide the individual into the dominant search space through the influence function, thus ensuring that evolution proceeds is in the dominant region, and for the *n*-dimensional optimization problem, the structure of the normative knowledge is described as follows:

$$\tag{3}$$

where *V<sup>i</sup>* = [(*l<sup>i</sup>* , *ui*),(*L<sup>i</sup>* , *Ui*)], *i* ≤ *n*. *u<sup>i</sup>* and *l<sup>i</sup>* are the upper and lower bounds of the *i*th dimensional variables, and *U<sup>i</sup>* and *L<sup>i</sup>* are the upper and lower bounds of the fitness value, respectively. The normative knowledge is updated with the change of the dominant search region, and gradually approaches the region where the best individual is located. Therefore, when there is a better individual in the *T*th generation beyond the current search range described by the normative knowledge, the normative knowledge is updated as follows:

$$l\_i(T+1) = \begin{cases} \mathbf{x}\_j^i(T) & \text{if } \mathbf{x}\_j^i(T) \le l\_i(T) \text{) or } f(\mathbf{x}\_j(T)) > L\_i(T) \\ l\_i(T) & \text{otherwise} \end{cases} \tag{4}$$

$$L\_i(T+1) = \begin{cases} f(\mathbf{x}\_j(T)) & \text{if } \mathbf{x}\_j^i(T) \le l\_i(T) \text{) or } f(\mathbf{x}\_j(T)) > L\_i(T) \\ L\_i(T) & \text{otherwise} \end{cases} \tag{5}$$

$$u\_i(T+1) = \begin{cases} x\_j^i(T) & \text{if } x\_j^i(T) \ge u\_i(T)) \text{ or } f(x\_j(T)) > \mathcal{U}\_i(T) \\ u\_i(T) & \text{otherwise} \end{cases} \tag{6}$$

$$\mathcal{U}\_{\boldsymbol{i}}(T+1) = \begin{cases} f(\boldsymbol{x}\_{\boldsymbol{j}}(T)) & \text{if } \boldsymbol{x}\_{\boldsymbol{j}}^{i}(T) \ge \boldsymbol{u}\_{\boldsymbol{i}}(T) \text{) or } f(\boldsymbol{x}\_{\boldsymbol{j}}(T)) > \mathcal{U}\_{\boldsymbol{i}}(T) \\ \mathcal{U}\_{\boldsymbol{i}}(T) & \text{otherwise} \end{cases} \tag{7}$$

(3) Historical knowledge. Historical knowledge was introduced into the belief space by Saleem [31] to record important events that occurred during the evolutionary process, and its main role is to adjust the offset distance and direction when the optimization falls into a local optima. The historical knowledge structure is divided as

shown in Figure 4, where *e<sup>i</sup>* is the *i*th outstanding individual of historical knowledge preservation, *W* is its maximum capacity, and *dsj* and *drj* are the average offset distance and the average offset direction of the *j*th design variable. The expressions of *dsj* and *drj* are as follows:

$$d\_{sj} = \frac{\sum\_{k=1}^{w-1} |e\_{k+1}\mathbf{x}\_j - e\_k\mathbf{x}\_j|}{w-1} \tag{8}$$

$$d\_{r\bar{j}} = \text{sgn}\left(\sum\_{k=1}^{w-1} |e\_{k+1}x\_{\bar{j}} - e\_k x\_{\bar{j}}|\right) \tag{9}$$

**Figure 4.** Structure of historical knowledge.

### *3.3. Proposed t-Mutation Operator*

Evolutionary algorithms require good exploration capabilities in the early stages and good exploitation capabilities in the later stages of evolution. The t distribution contains the degree of freedom parameter n, which approaches the standard Gaussian distribution infinitely when *n* → +∞ and the t distribution is the standard Cauchy distribution when *n* = 1. That is, the standard Gaussian distribution and the standard Cauchy distribution are two boundary special cases of the t distribution. The probability density functions of the standard Gaussian distribution and the standard Cauchy distribution are shown in Figure 5. Obviously, the application of the Cauchy operator can produce a larger mutation step, which is conducive to the algorithm to guide individuals to jump out of the local optimal solution and ensure the exploration ability of the algorithm, and Gaussian distribution shows a better exploitation ability.

**Figure 5.** Probability density functions of the Cauchy and Gaussian distributions.

Population diversity is considered as the primary reason for premature convergence, which determines the search capability of the algorithm. In evolutionary algorithms, population diversity decreases over time as evolution proceeds. Therefore, population diversity can be used to determine the stage of evolution; we can thus use the population diversity to construct the degree of freedom n. By changing the degree of freedom parameter n, the mutation scale changes adaptively with evolution to balance the exploitation and exploration capabilities of the algorithm. In this paper, we introduce population variance and population entropy to determine population diversity. The expression of population variance *D<sup>T</sup>* in the *T*th generation is as follows:

$$D\_T = \frac{1}{N \times l} \left( \sum\_{i=1}^{N} \sum\_{j=1}^{l} x\_i^j - \overline{x^j} \right) \tag{10}$$

where *x j i* is the *j*th gene value of the ith individual, *N* is the number of populations and *l* is the individual coding length. The expression of *x j* is as follows:

$$\overline{\boldsymbol{x}^{j}} = \frac{1}{N} \left( \sum\_{i=1}^{N} \boldsymbol{x}\_{i}^{j} \right) \tag{11}$$

The solution space *A* of the optimization problem is divided equally into *L* small spaces, and the number of individuals belonging to the ith space *A<sup>i</sup>* in generation *T* is |*A<sup>i</sup>* |. The expression of population entropy *S<sup>T</sup>* in the *T*th generation is as follows:

$$S\_T = -\sum\_{i=1}^{L} p\_i \log(p\_i) \tag{12}$$

where

$$p\_i = \frac{|A\_i|}{N} \tag{13}$$

From the definitions of population variance and population entropy, it is clear that population variance reflects the degree of dispersion of individuals in the population and that population entropy reflects the number of individual types in the population. Therefore, the *t*-mutation operator *t*(*n*) can be constructed based on the population variance and population entropy. The degree of freedom parameter *n* is expressed as follows:

$$m = \left[1 - \ln\left(\frac{D\_T + S\_T}{D\_{\max} + S\_{\max}}\right)\right] \tag{14}$$

where [ ] is the least integer function, and *Dmax* and *Smax* are the maximum values of population variance and population entropy, respectively. Obviously, the degree of freedom parameter n of the *t*-mutation operator is 1 in the first generation and increases gradually as evolution proceeds, then the degree of freedom parameter *n* converges to positive infinity in the late evolutionary stage, and the *t* distribution becomes a standard Gaussian distribution. The *t*-mutation operator can ensure the exploration capability of the algorithm in the early evolutionary stage and the exploitation capability of the algorithm in the late evolutionary stage.

### *3.4. Communication Protocol*

The information interaction between the belief space and the population space is realized through the acceptance function and the influence function. The acceptance function passes the better individuals in the population space as samples to the belief space for knowledge sources extraction and update, and the influence function is the way to influence the population space by the belief space, which can use the knowledge sources in the belief space to guide the population space to complete and accelerate the evolution.

### 3.4.1. Acceptance Function

In this paper, a dynamic version of the acceptance function [31] is used. The number of accepted individuals is given as follows:

$$n\_{Accept} = \left(p\% + \frac{p\%}{T}\right) \times N \tag{15}$$

where [ ] is the least integer function, *T* is the current generation, *N* is the number of populations, *p*% is the preset fixed proportion and *p*% = 20%. In this paper, the acceptance function accepts *nAccept* better individuals into the belief space. The dynamic acceptance function makes the number of individuals entering the belief space decrease with the depth of evolution, which increases the global search ability of the algorithm at the early stage of evolution, and reduces the number of individuals entering the belief space at the late stage of evolution because the population tends to converge and carries mostly similar information, which can maintain the diversity of knowledge sources and avoid the consumption of memory by redundant information.

### 3.4.2. Influence Function

The core of the influence function is the manner and proportion in which each type of knowledge affects the population. Knowledge acts on each type of influence function to control the number of individuals affected by each type of influence function. Therefore, the proportion by which each type of knowledge affects the population is the relative role that each type of influence function has in the population. The proportion of the effect of the influence function is determined by the success rate of the knowledge effect, and is expressed as follows:

$$P\_k(T) = \begin{cases} 1/N\_k & \text{if } T = 1\\ \alpha + \beta \times \frac{\upsilon\_k(T-1)}{\upsilon(T-1)} & \text{if } (T \bmod K) = 0\\ P\_k(T-1) & \text{otherwise} \end{cases} \tag{16}$$

It satisfies the condition that *β* + *αN<sup>k</sup>* = 1, where *N<sup>k</sup>* is the number of knowledge sources types, *v<sup>k</sup>* (*T* − 1) denotes the number of individuals influenced by knowledge *k* that are better than their parents in generation *T* − 1 and *v*(*T* − 1) denotes the number of individuals influenced by all knowledge sources that are better than their parents in generation *T* − 1. The success rate of knowledge sources influenced in the previous generation determines the proportion of the effect of each knowledge source in the next generation. In order to allow each kind of knowledge source to always have the possibility of being used, we took *α* = 0.1, *β* = 0.7, ensuring that the lower bound of *P<sup>k</sup>* is 0.1 and the proportion of all knowledge sources in the first generation is the same, which is 1/*N<sup>k</sup>* .

Next, we introduced the proposed *t*-mutation operator into the influence function to develop a knowledge-guided *t*-mutation strategy.

(1) Situational knowledge. Situational knowledge has a guiding role in the evolutionary process, and the effect of situational knowledge on the population space under the action of the *t*-mutation operator is noted as follows:

$$\mathbf{x}\_{ij}^{'} = \begin{cases} \mathbf{x}\_{ij} + |(\mathbf{x}\_{ij} - E\_j)t(n)| & \text{if } \mathbf{x}\_{ij} < E\_j \\ \mathbf{x}\_{ij} - |(\mathbf{x}\_{ij} - E\_j)t(n)| & \text{if } \mathbf{x}\_{ij} > E\_j \\ \mathbf{x}\_{ij} + \gamma t(n) & \text{otherwise} \end{cases} \tag{17}$$

where *xij* is the *j*th dimensional design variable of the ith individual, *x* 0 *ij* is the *j*th dimensional design variable of the newly generated *i*th individual, *γ* is a constant and *E<sup>j</sup>* is the *j*th dimensional design variable of the situational knowledge.

(2) Normative knowledge. The normative knowledge guides the population to search in the dominant region, and the effect of the normative knowledge on the population space under the action of the *t*-mutation operator is noted as follows:

$$\mathbf{x}\_{ij}^{'} = \begin{cases} \mathbf{x}\_{ij} + |(u\_j - l\_j)t(n)| & \text{if } \mathbf{x}\_{ij} < l\_j\\ \mathbf{x}\_{ij} - |(u\_j - l\_j)t(n)| & \text{if } \mathbf{x}\_{ij} > u\_j\\ \mathbf{x}\_{ij} + \mu(u\_j - l\_j)t(n) & \text{otherwise} \end{cases} \tag{18}$$

where *µ* is a constant, and *u<sup>j</sup>* and *l<sup>j</sup>* are the upper and lower bounds of the *j*th dimensional design variables preserved by the normative knowledge of the current generation belief space, respectively.

(3) Historical knowledge. Historical knowledge is used to adjust the offset distance and direction when the optimization is trapped in a local optima, and the effect of historical knowledge on the population space under the action of the *t*-mutation operator is noted as follows:

$$\mathbf{x}\_{ij}^{'} = \begin{cases} \varepsilon \mathbf{x}\_{j} + d\_{rj} t(n) & \text{for } 45\% \text{ of the time} \\ \varepsilon \mathbf{x}\_{j} + d\_{sj} t(n) & \text{for } 45\% \text{ of the time} \\ \ random (l\_{j}^{'} \boldsymbol{u}\_{j}^{'}) & \text{for } 10\% \text{ of the time} \end{cases} \tag{19}$$

where *ex<sup>j</sup>* is the jth dimensional design variable of the best individual ex stored in the historical knowledge and *u* 0 *j* and *l* 0 *j* are the upper and lower bounds, respectively. Here a roulette wheel is used to determine how new individuals are generated, with a 45% probability that individuals produce a bias in direction, a 45% probability that individuals produce a bias in distance and a 10% probability that new individuals are generated randomly within the entire search space [32].

### *3.5. The Main Numerical Implementation of HCGA*

The main numerical implementation of HCGA is described step-by-step as follows:


**Step 9:** Stop the algorithm if the stopping criterion is satisfied; otherwise *T* = *T* + 1 and go to Step 4.

### **4. Numerical Validation and Performance of Hybrid HCGA Algorithm**

*4.1. Parameter Discussion*

Tuning parameters properly is very important for an evolutionary algorithm to achieve good performance. In HCGA, there are seven main parameters: *Pc*, *Pm*, *α*, *γ*, *µ*, *K*, *p*. In this section, we used the factorial design (FD) [16] approach in order to obtain a guideline on how to tune the designed parameters in HCGA.

Ten benchmark mathematical optimization problems were used to evaluate and compare optimization algorithms. These functions can be divided into unimodal functions and multimodal functions. Functions *F*<sup>1</sup> − *F*<sup>4</sup> are unimodal functions with only a global optimal value, which are mainly used to evaluate the exploitation ability and convergence speed of the algorithm. Functions *F*<sup>5</sup> − *F*<sup>10</sup> are multimodal functions, which have multiple local optimal values in the search space, and the number of local optima will increase with the increase of the problem size, which is an important reference for assessing the exploration capability of the algorithm. Seven of these test functions (*F*<sup>1</sup> − *F*7) are dimension-wise scalable. The details of the test functions are listed in Table 1.

**Table 1.** Details of the mathematical optimization problems.


In the experiments, the population size was set to twice the dimension for the *F*<sup>1</sup> − *F*<sup>7</sup> function and five times the dimension for the *F*<sup>8</sup> − *F*<sup>10</sup> function, and the *Tmax* was set to 5000. As shown in Table 2, we used seven parameters as factors for seven levels in an orthogonal experimental design. Table 3 shows the test results of the orthogonal parameter table with the 10*D* − *f*<sup>1</sup> function. Trials of 30 times were performed for each set of parameters. The unabridged result tables, similarly to Table 3 of other experiments, were too large, and they were omitted here.

As shown in Table 3, to estimate the effects of each set of parameters, the mean fitness of the 30 runs were obtained and listed in in the last column of the table. *K<sup>i</sup>* is the mean value of mean fitness for this column parameter at level *i* (*i* = 1, 2, . . . , 7). Std is the standard deviation of each column *K*<sup>1</sup> − *K*7. The larger the Std value is, the more this column parameter influences the algorithm performance. Furthermore, for each column, if the value of *K<sup>i</sup>* is the smallest *K* value in that column, then the best value of the parameter is the parameter value on level *i*. The best parameters (*B–P*) are listed in the last row. Tables 4 and 5 show the Std and *B–P* of all benchmark functions. The symbol ∼ in Table 5 indicates that each set of parameters enables the algorithm to optimize to the same optimal value.


**Table 2.** Factors and levels for orthogonal experiment.

**Table 3.** Orthogonal parameter table of *L*49(7) <sup>7</sup> and experimental results on *f*<sup>1</sup> (Dimension = 10).



**Table 3.** *Cont.*

**Table 4.** Standard deviations of orthogonal experiments.



**Table 5.** Best parameters (*B–P*) of orthogonal experiments.

From Table 4, it can be seen that for low-dimensional functions with unimodal functions, *P<sup>c</sup>* has a greater influence on the algorithm performance, while for high-dimensional complex functions, it is *p*, *α* and *K* that have a greater influence. This indicates that when dealing with simple functions, the population space plays a major role, and when dealing with complex functions, the belief space plays a guiding role and has an influence on the evolution of the population space.

Some rules for adjusting parameters can be obtained from analyzing the results in Table 5. For simple functions, *P<sup>c</sup>* and *P<sup>m</sup>* can be set to a lower level, while for complex functions, they need to be set to a higher level. For most functions, *α* can be set to a level of about 0.2, and for multimodal functions with many local minima, *α* should be set to 0.1. For parameter *γ*, setting it to 0.3 is enough in most cases. *µ* has roughly the same rule as *P<sup>c</sup>* and *Pm*. For high-dimensional multimodal functions, *µ* should be set to 0.1, but for unimodal or low-dimensional functions, setting it at 0.3–0.4 is enough. *K* should be set to a smaller value as the complexity of the function increases, which determines the frequency of updating the knowledge in the belief space. As for *p*, setting it at 0.1–0.2 should be enough for both unimodal and multimodal functions.

### *4.2. Validation in Numerical Experiments*

In order to verify the performance of the algorithms, cultural algorithm (CA) [16], genetic algorithms (GAs) [11], differential evolution (DE, rand/1/L) [33] and HCGA were selected for comparison with numerical experiments. Ten mathematical functions optimization test problems shown in Table 1 were used to compare the performance of HCGA with GAs, CA and DE.

The parameters in HCGA were selected based on the results of the parameter discussion in Section 4.1, and the parameters of each algorithm in the experiments were set as shown in Table 6. Since evolutionary algorithms are essentially stochastic optimization algorithms, the solution found may not be the same each time. Therefore, each benchmark function was repeated 30 times.


**Table 6.** Main parameters of the HCGA, GAs, CA and DE.

The optimal values, means and standard deviations of HCGA, GAs, CA and DE for 30 independent runs are listed in Table 7, which were used to evaluate the optimization accuracy, average accuracy and stability of the algorithms. To obtain more reliable statistical conclusions, Wilcoxon nonparametric statistical tests were performed at *α* = 0.05, and the symbols +, − or = mean that the optimization results of HCGA were significantly better, worse or similar to the comparison algorithm, respectively. Figure 6 shows the convergence curves of some of the benchmark test functions. The results are summarized as +/ − / = as the last row of each Table in Table 7.

As can be seen from Figure 6, HCGA shows a better performance for most functions. For unimodal functions, the convergence speed and accuracy of HCGA are significantly better than those of other algorithms. For multimodal functions, HCGA is able to achieve higher optimization accuracy in shorter iterations for all functions except for function *F*10, which was slightly inferior to GAs and DE in terms of convergence speed in the initial search stages. This means that HCGA not only has good search ability and fast convergence, but also moderates quite well the conflict well between convergence speed and premature convergence, which means that it has a balanced exploitation and exploration ability.

The experimental results in Table 7 show that HCGA performs better for most of the tested functions compared with other algorithms, and can obtain higher optimization accuracy, average accuracy and better stability. This indicates that HCGA is less affected by randomness and that it can maintain optimization accuracy under multiple independent runs. The results of Wilcoxon nonparametric statistical tests for CA, GAs and DEs were 23/0/1, 22/1/1, and 21/3/0, respectively, indicating that the differences between the HCGA and the other three compared algorithms are statistically significant, implying that for all test functions, HCGA shows better performance or is close to the best performance of the other algorithms, which means that it is more robust.

In addition, HCGA shows an optimization capability for high-dimensional problems that cannot be matched by CA, GAs and DE. In high-dimensional optimization problems (100 dimensional *F*<sup>1</sup> − *F*7), HCGA has significant advantages in optimization accuracy, average accuracy and stability. For the functions *F*<sup>4</sup> − *F*7, the number of local optima will increase with the increase in the problem size, and the HCGA does not fall into dimensional disaster; it also scales well with the increasing dimensionality and converges in the proximity of the global optimum, which indicates its high level of performance in solving high-dimensional functions. HCGA can still maintain strong optimization accuracy and robustness in solving high-dimensional optimization problems, which lays the foundation for the application of HCGA in practical problems.

*Appl. Sci.* **2022**, *12*, 3482



**Figure 6.** Convergence curves of the test functions.

### *4.3. Mechanistic Analysis of Improved Hybrid Algorithm Performance*

Considering the benchmark function optimization results in Section 4.2, it is obvious that HCGA is superior compared to CA and GAs. The mechanistic analyses of improved hybrid algorithm performance are as follows:


The benchmark results obtained with mathematical functions and the above analysis demonstrate that HCGA is an efficient optimization algorithm with potential applications for complex optimization problems.

### **5. Applications to Aerodynamic Design Optimization of Wing Shapes**

The aerodynamic shape optimization design of a wing is one of the important components of aircraft configuration design, and it has been the goal of researchers to design the aerodynamic shape of a wing for decades in terms of efficiency and quality to meet engineering objectives. Cruise factor is one of the most important aerodynamic characteristics that determine the performance of an aircraft. The objectives of this section are introducing and using HCGA for the aerodynamic design optimization of a wing to achieve the cruise factor optimization.

### *5.1. Parameterization Strategy*

Airfoil parameterization is a crucial step in aerodynamic optimization, and its accuracy determines the accuracy and reliability of the optimized airfoil. The commonly used parameterization methods are the free-form deformation (FFD) technique [34], Bezier curves [35], the class/shape transformation (CST) method [36], etc. In this work, a fourorder CST parameterization method is used to control the airfoil shape, and the parametric expressions of the upper and lower surface curves are defined in Equations (20) and (21).

The design variables are the leading-edge radius *Rle* of the airfoil, the inclination angles *b<sup>i</sup>* and *b* 0 *i* of the upper and lower surface curves at the trailing edge and the upper and lower surface shape function control parameters *β* and *β* 0 . For a total of nine airfoil design parameters, with the reference geometry being the RAE2822 airfoil, the design parameters and corresponding constraint ranges are shown in Table 8.

$$\frac{\frac{\pi}{c}}{\frac{\pi}{c}} = \sqrt{\frac{\pi}{c}} \left( 1 - \frac{\underline{x}}{c} \right) \left\{ \sqrt{\frac{2R\_{\ell}}{c}} \left( 1 - \frac{\underline{x}}{c} \right)^{\underline{n}} + \tan \beta \left( \frac{\underline{x}}{c} \right)^{\underline{n}} + \sum\_{i=1}^{n-1} \left[ b\_{i} \frac{\underline{n}!}{\overline{\iota}!(\underline{n}-i)!} \left( \frac{\underline{x}}{c} \right)^{i} \left( 1 - \frac{\underline{x}}{c} \right)^{\underline{n}-\bar{i}} \right] \right\} \tag{20}$$

$$\frac{\pi}{c} = \sqrt{\frac{\pi}{c}} \left( 1 - \frac{\underline{x}}{c} \right) \left\{ -\sqrt{\frac{2R\_{\parallel c}}{c}} \left( 1 - \frac{\underline{x}}{c} \right)^{\underline{n}} + \tan \beta \left( \frac{\underline{x}}{c} \right)^{\underline{n}} - \sum\_{i=1}^{n-1} \left[ b\_{i \cdot \overline{1(n-i)}}^{\prime} \left( \frac{\underline{x}}{c} \right)^{i} \left( 1 - \frac{\underline{x}}{c} \right)^{\underline{n}-i} \right] \right\} \tag{21}$$


**Table 8.** Range of design parameter.

### *5.2. Wing Shape Optimization*

In this design, eight sections were used to describe the whole wing geometry for its shape optimization, its configuration and control surface distribution, as shown in Figure 7. The parameterization method described in Section 5.1 was used to control the wing shape, with a total of 72 design variables. The optimal design of the wing for the cruise factor was considered in the cruising condition at the flow condition of Mach 0.785, a 1.92◦ angle of attack and a Reynolds number of <sup>2</sup> <sup>×</sup> <sup>10</sup><sup>7</sup> based on the aerodynamic mean chord. HCGA was used to optimize the shape with a population size of 150 and evolutionary iterations of 100. The objective was to maximize the cruise factor, and the constraints were that the maximum thickness of each control surface and the lift coefficient should not to be reduced. The mathematical optimization model is described as follows:

$$\max: f(\mathbf{x}) = Ma \frac{L}{D} \tag{22}$$

$$\text{s.t.}: t\_{\text{i}} \ge t\_{\text{Initial}} (\text{i} = 1, 2, \dots, 8), \mathcal{C}\_{\text{L}} \ge \mathcal{C}\_{\text{L0}} \tag{23}$$

where *Ma* is Mach number, *L* is lift, *D* is drag, *t<sup>i</sup>* is the maximum thickness of the *i*th section, *tinitial* is the maximum thickness of the initial control surface, *C<sup>L</sup>* is the lift coefficient and *CL*<sup>0</sup> is the lift coefficient of the initial wing. The flow was modeled by the compressible full potential flow with viscous boundary layer correction, and the total number of mesh points was about 0.5 million. The pressure coefficient contours on the upper and lower surface of the initial and optimized wing are shown in Figure 8. The respective pressure distributions at each section are shown in Figure 9 and the wing section shapes before and after optimization are compared in Figure 10. It is seen that the shock waves were significantly smeared owing to the shape modification, which resulted in a considerable reduction of wave drag on the upper surface and therefore in better aerodynamic performance. The aerodynamic parameters of the wing before and after optimization are shown in Table 9. The cruise factor *Ma* × *L*/*D* was significantly increased from 21.863 to 26.938 because the drag coefficient *C<sup>D</sup>* of the optimized wing was obviously reduced from 0.01605 to 0.01302. It can be seen that the cruise factor increased by 23.21%, while the drag coefficient decreased by 18.88% and the constraints of lift coefficient and thickness were satisfied. There was also no significant change in the induced drag coefficient *CDIND* since there was no change in the lift coefficient *CL*. The wave drag coefficient *CDwave* of the wing was reduced apparently from 0.00240 to 0.00030, and the profile drag coefficient *CDPD* was also reduced from 0.00680 to 0.00602.

**Table 9.** Comparison of airfoil aerodynamic parameters.


**Figure 7.** Distribution of the wing sections.

**Figure 8.** (**a**) Pressure coefficient contours on the wing upper surface of the initial wing; (**b**) Pressure coefficient contours on the wing upper surface of the optimized wing; (**c**) Pressure coefficient contours on the wing lower surface of the initial wing; (**d**) Pressure coefficient contours on the wing lower surface of the optimized wing.

**Figure 9.** Comparison of surface pressure coefficient distributions for initial and optimized wings.

**Figure 10.** Comparison of initial and optimized wing-section shapes.

For a better comparison of values with the proposed HCGA algorithm in this engineering application, the commonly used GAs [11] and PSO [37] of the engineering optimization field were selected for comparison with HCGA. The parameter settings of HCGA were the same as for the numerical experiment, and all parameters of GAs and PSO were default parameters. The population size and maximum number of iterations for all three algorithms were 150 and 100, respectively. The cruise factor convergence curve is shown in Figure 11.

**Figure 11.** Convergence curves of cruising factor.

It can be seen that the optimization results of HCGA were significantly better than those of GAs and PSO for the same number of iterations. It can be observed that HCGA is obviously a more efficient algorithm for aerodynamic optimization problems, which can achieve better-quality optimized results with fewer flow field calculations and can significantly improve the efficiency of aerodynamic optimization.

### **6. Conclusions**

In this paper, an efficient hybrid evolutionary optimization method coupling CA with GAs (HCGA) was proposed to improve the efficiency of the optimization procedure for the aerodynamic shape of an aircraft. HCGA aims to improve the ability to solve complex problems and increase the efficiency of optimization. To improve the robustness of the algorithm, HCGA uses GAs as an evolutionary model of the population space. HCGA constructs the belief space using three kinds of knowledge: situational knowledge, normative knowledge and historical knowledge. Meanwhile, the knowledge-guided *t*mutation operator was developed to dynamically adjust the mutation step and balance the exploitation and exploration ability of the algorithm. The optimization performance of HCGA was demonstrated on many benchmark functions for which the global optima are known a priori. The optimization results obtained with the benchmark functions show that HCGA provides a better global convergence, a better convergence speed and a better optimization accuracy compared to CA and GAs. In particular, HCGA shows the potential for solving large-scale design variable optimization problems.

By combining HCGA with a CFD solver, an efficient decision-maker design tool for aerodynamic shape design optimization was developed to find the best aerodynamic shape to satisfy the design requirements. For the three-dimensional wing design problem, the proposed HCGA optimizer successfully reduced the wing drag computerized design, thus significantly improving the wing cruise factor. Compared with the baseline wing, the drag coefficient was reduced by 18.88%, which resulted in a 23.21% improvement in the cruise factor. This proves the capability and potential of HCGA for solving complex engineering design problems in aerodynamics. As a practical engineering application of the super-heuristic algorithm, the potential and value of such algorithms for engineering applications are further validated.

However, this study is only preliminary and further testing is needed to evaluate the performance of HCGA in complex engineering optimization. In addition, the practical application of HCGA only considered single-objective optimization. Multi-objective optimization problems should thus be the next step for investigation in future research.

**Author Contributions:** The contributions of the five authors in this paper are conceptualization, X.Z., Z.T. and J.P.; methodology, X.Z. and Z.T.; formal analysis, X.Z., Z.T. and F.C.; investigation, X.Z. and Z.T.; resources, F.C. and C.Z.; data curation, X.Z. and Z.T.; writing—original draft preparation, X.Z. and Z.T.; writing—review and editing, X.Z. and Z.T.; supervision, Z.T.; project administration, C.Z.; funding acquisition, Z.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China (NSFC-12032011, 11772154) and the Fundamental Research Funds for the Central Universities (NP2020102).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The figures, tables and data that support the findings of this study are mentioned in the corresponding notes, with reference numbers and sources, and are publicly available in the repository.

**Conflicts of Interest:** The authors declare no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

## **References**

