Aerodynamic Shape Optimization    AIAA 2001-2473

Using A Real-Number-Encoded Genetic Algorithm

 

Terry L. Holst* and Thomas H. Pulliam

NASA Ames Research Center

Moffett Field, CA 94035

 

 


Abstract

 

A method for aerodynamic shape optimization using a genetic algorithm with real number encoding is presented. The algorithm is used to optimize three different problems, a simple hill climbing problem, a two-dimensional airfoil problem using an Euler/Navier-Stokes equation solver and a three-dimensional transonic wing problem using a nonlinear potential solver. Results indicate that the genetic algorithm is easy to implement, flexible in application and extremely reliable, being relatively insensitive to design space noise.

 

Background

 

Numerical methods for optimizing aerodynamic performance have been studied for many years. Perhaps the most widely used general approach involves the computation of sensitivity gradients. These methods—called gradient methods—have been utilized to produce optimal aerodynamic performance in a wide variety of different forms. The reliability and success of gradient methods generally requires a smooth design space and the existence of only a single global extremum. A good review of gradient methods used in aerodynamic design is presented by Reuther.1

 

 

–––––––––––––––––––––––

*Research Scientist, Numerical Aerospace Simulations Division, Applications Branch. Fellow AIAA.

Research Scientist, Numerical Aerospace Simulations Division, Applications Branch. Associate Fellow AIAA.

 

Copyright Ó 2001 by the American Institute of Aeronautics and Astronautics, Inc. No Copyright is asserted in the United States under Title 17, U. S. Code. The U. S. Government has a royalty-free license to exercise all rights under the copyright claimed herein for Government Purposes. All other rights are reserved by the copyright owner.

A key characteristic among all gradient methods lies is how the gradients or sensitivity derivatives are computed. The simplest approach, often called the “brute-force” or finite-difference method, consists of using CFD flow solver solutions to determine the effect of each design perturbation on the objective function. Sensitivities are then constructed using finite-difference formulas. If a design problem using the above approach has K decision variables, the sensitivity derivative computation for each design iteration will require K+1 function evaluations, i.e., one CFD solution for the unperturbed or baseline geometry and K solutions corresponding to the K decision variable perturbations. Hicks and Henne2 presents one study where the finite-difference approach is used to compute sensitivity derivatives.

 

Other methods for evaluating sensitivity derivatives that seek to reduce the large computational cost associated with the finite-difference method have been developed, e.g., the ADIFOR (Automatic Differentiation of Fortran) approach,3 the quasi-analytic method of El-banna and Carlson4,5 and Arslan and Carlson6 and the adjoint method of Jameson.7,8 

 

In the latter approach, an adjoint equation, derived from control theory, is solved to determine design space sensitivity derivatives. This approach is superior to the finite-difference approach for generating sensitivities because all the sensitivities are obtained (no matter how many) by solving one adjoint equation, with a cost on the order of one flow field solution.  Thus, for design problems containing a large number of decision variables, the cost savings for this approach over the finite-difference approach is significant. 

 

In contrast to gradient based methods, design space search methods such as genetic algorithms (GA) offer an alternative approach with several attractive features. The basic idea associated with the GA approach is to search for optimal solutions using an analogy to the theory of evolution. During solution iteration (or “evolution” using GA terminology) the decision variables or “genes” are manipulated using various operators (selection, combination, crossover or mutation) to create new design populations, i.e., new sets of decision variables.  General GA details can be found in Goldberg,9 Davis10 and Beasley, et al.11 Each design is evaluated using an objective-like “biological fitness function” to determine survivability. 

 

Constraints can easily be included in this approach either by direct inclusion of the fitness function or by preprocessing the candidate design.  For example, if a design violates a constraint, its fitness is set to zero, i.e., it doesn’t survive to the next evolution level. 

 

Because GA optimization is not a gradient-based optimization technique, it does not need sensitivity derivatives.  It theoretically works well in non-smooth design spaces containing several or perhaps many local extrema. It is also an attractive method for multi-objective design applications offering the ability to compute so-called “pareto optimial sets” instead of the limited single design point traditionally provided by other methods.

 

A disadvantage of the GA approach is expense.  In general, the number of function evaluations required for a GA algorithm exceeds the number required by a finite-difference-based gradient optimization (see the results presented in Obayashi and Tsukahara12 and Bock13).

 

Example applications utilizing potential-based flow solvers in the context of GA optimization can be found in Quagliarella and Della Cioppa14 for airfoil applications, Vicini and Quagliarella15 for multi-point and multi-objective airfoil applications and Obayashi et al.16 for multi-disciplinary optimization of transonic wings.  Examples of multi-design-point wing optimization using Euler and Navier/Stokes flow solvers can be found in Sasaki, et al.17 and Oyama.18,19

 

The Genetic Algorithm used in the present study is described next. Details associated with each of the operators, including selection, passthrough, random average crossover, perturbation mutation and mutation are presented. The GA is then evaluated using a number of optimization problems including a simple multi-modal hill climbing search optimization, two different types of airfoil optimizations (one transonic inviscid and the other viscous) and a transonic wing optimization.

 

Genetic Algorithm

 

Several GA variations have been used in the present study. The GA now presented is a representative example. As mentioned above the general idea behind any GA optimization is to discretely describe the design space using a number of genes, . In this notation the i-subscript is the gene number and the superscript is the generation. Each set of genes that leads to the complete specification of an individual design is called a chromosome

 

                    (1)

 

In this notation is the jth chromosome for the nth generation that consists of NG genes. A second subscript j has been added to each gene to show which chromosome it is identified with.

 

For aerodynamic shape optimization the set of genes associated with a single chromosome must fully describe the geometry being modified. In many GA applications genes are computationally represented using bit strings and the operators used are designed to manipulate bit string data. In the present approach, following the arguments of Oyama,19  Houck, et al.20 and Michalewicz,21 real-number encoding is used to represent all genes. The key reason for using real number encoding is that it has been shown to be more efficient in terms of CPU time relative to binary encoded GAs.21 In addition, real numbers are used for all parameters defining an aerodynamic surface or shape, e.g., sweep, twist, thickness, camber. Thus, using a real number encoding eliminates the need for binary-real number conversions.

 

Once the design space has been defined in terms of a set of real-number genes, the next step is to form an initial generation, ,  represented by

 

                      (2)

 

where NC is the total number of chromosomes. Each gene within each chromosome is assigned an initial numerical value using a process that randomly chooses numbers between fixed user-specified limits. For example, the tenth gene in an arbitrary chromosome is computed using

 

                  (3)

 

where  and  are the upper and lower limits for the tenth gene, respectively, and R(0,1) is a random number generator that delivers an arbitrary numerical value between 0.0 and 1.0.

 

After the initial generation is established, fitness values, , are computed for each chromosome using a suitable function evaluation. This is analogous to the objective function evaluation in gradient methods and is represented using

 

                                             (4)

 

For aerodynamic shape optimization F represents a suitable CFD flow analysis. Fitness determination is followed by a ranking process where the most fit individual is given a number one ranking (IR = 1), the second most fit individual is ranked number two (IR = 2), and so on. This is determined using

 

 

This completes the GA initialization process.  The next several subsections describe how the GA progresses from generation to generation using a variety of operators that manipulate the chromosomes allowing evolution to occur.

 

Selection

 

The first operation required to determine the (n+1)st generation is selection. The chromosomes that will be used by the other GA operators (to be discussed shortly) must be selected from the nth generation chromosomes. The selection operation used in the present study is given by

 

where the selected chromosome has been placed in a temporary holding array given by

 

 

Note how the fittest individuals in the nth generation are selected multiple times, the average individuals are selected a small number of times, and the least fit individuals are not selected at all. This biasing toward the fittest individuals is a key element in any GA. The chromosomes represented by  are used by succeeding operators to produce .

 

Passthrough

 

The simplest operator used in the present study is “passthrough.” As the name implies, a certain number of the fittest chromosomes are simply “passed through” to the next generation. Because passthrough is always performed first, it operates on the chromosomes that have the highest fitness. This guarantees that the maximum fitness never drops from generation to generation. The number of chromosomes that are passed through to the next generation is controlled by the parameter pB. For example, if pB = 0.1 then ten percent of the chromosomes will be passed through to the next generation.

 

Random Average Crossover

 

The next GA operator to be discussed is called the random average crossover operator and is implemented gene-by-gene using the following formula:

 

         (5)

 

where the j1 and j2 subscripts are randomly chosen between 1 and NC. Values for  are taken from  and the newly computed values  are stored in . The number of chromosomes modified using random average crossover is determined by the parameter pA.  For example, if pA = 0.2, then twenty percent of the chromosomes will be determined for the (n+1)st generation using random average crossover.

 

Perturbation Mutation

 

The next GA operator is called perturbation mutation and is implemented by first selecting a random chromosome  from . Then a single gene  is randomly selected from  and modified using the following formula:

 

             (6)

 

where  is a user-specified tolerance. Because this operator can cause the value of a particular gene to exceed one of its constraints (xmaxi or xmini), checks are required to make sure this doesn’t happen. The number of chromosomes modified using perturbation mutation is determined by the parameter pP. For example, if pP = 0.3, then thirty percent of the chromosomes will be determined for the (n+1)st generation using the perturbation mutation operator.

 

Mutation

 

The last GA operator used in the present study is called the mutation operator and is implemented similarly to the perturbation mutation operator. First, a random chromosome  is chosen from . Then a single gene  is randomly selected from  and supplied with a completely different value using the following formula:

 

      (7)

 

Like the other operators used in the present study, the mutation operator is controlled by a parameter, pM. For example, if pM = 0.4, then forty per cent of the chromosomes will be determined for the (n+1)st generation using the mutation operator.

 

General Algorithm Comments

 

For consistency the parameters (pB, pA, pP, pM) must sum to one. The passthrough operator is always performed first and always passes through the top pB chromosomes from  to . Otherwise, the order in which each operation is performed is immaterial. Once all values of  have been established, the algorithm proceeds to fitness evaluation, ranking and on to succeeding generations until the optimization is sufficiently converged.

 

Results

 

Case 1: Hill-Climbing Problem

 

The first problem used to evaluate the GA is a simple hill-climbing problem. It utilizes a continuously differentiable analytic function of the form  that has several peaks and valleys, (the so called peak function from MATLAB). An isometric graphical view of this function is displayed in Fig. 1 over the range –3 ≤ x ≤ 3 and –3 ≤ y ≤ 3. The object of this exercise is to find the maximum value of z using the GA and in so doing gain insight into the workings of the GA process. Thus, the x and y values are the genes, each (x,y) pair is a chromosome and the value of z for each chromosome is the fitness.

 

Fig. 1 Isometric view of the function used in the hill-climbing problem.

 

Two typical GA convergence histories for the hill-climbing problem are presented in Fig. 2. For each curve the number of chromosomes, NC, was 20 and the  parameter in the perturbation mutation operator was 0.01. For each of these computations the GA procedure was continued until the error in the solution was reduced below 0.00001, i.e., CONV = 10-5. The P vector notation appearing in the Fig. 2 caption is defined by P = (pB, pA, pP, pM). The maximum fitness value at the end of each GA generation is plotted. As can be seen there is not that much difference between the two convergence histories, although the case with more mutation (higher values of pP and pM) and less passthrough (smaller value of pB) does have slightly superior convergence.

 

Using the hill-climbing problem to study the effect of P, NC and  on GA convergence is handy because of the speed with which each complete GA optimization can be performed. Such a study is presented in Figs. 3 and 4. For these results a value of CONV = 10-5 was used for each computation. Each curve was generated by performing a separate GA over a large number of NC values ranging from 12 to 1024 in increments of 4. The resulting data were curve fit to produce each of the curves displayed.

 

Fig. 2 Sample GA convergence histories for hill-climbing problem,  = 0.01, CONV = 10-5, NC = 20.

 

Figure 3 presents the effect of wide variations in NC and  on GA convergence. Generally, the most optimal convergence (in terms of minimizing function evaluations) occurs when the number of chromosomes is small, on the order of 12-20, and when  is between 1.0 and 0.001. Except for extremely small values (10-4 and smaller)  has little effect on GA performance. This is understandable since the perturbation mutation operator serves the role of providing small arbitrary corrections to the convergence process. As  is reduced in size the corrections become too small and GA convergence is delayed.

 

Fig. 3    Genetic algorithm convergence as a function of the number of chromosomes and the size of perturbation mutations (), CONV = 10-5, P = (0.1, 0.2, 0.3, 0.4).

 

Fig. 4    Genetic algorithm convergence as a function of the number of chromosomes and the type of operators used, CONV = 10-5,  = 0.01.

 

For optimal values of  the number of chromosomes used in each calculation has a relatively significant effect on convergence. This is due to the fact that a reasonable number of generations must be selected and operated upon in order to achieve convergence. As the number of chromosomes is increased, the number of selection operations that are performed for a fixed number of function evaluations decreases. In the limit, when the number of chromosomes is equal to the number function evaluations, no selection, crossover or mutation operations are performed. The present algorithm is relegated to a simple random search and is quite inefficient. This is the mechanism that produces slower convergence for larger values of NC.  Of course, break down will occur for any GA if the number of chromosomes becomes too small as well. Because of the present algorithm’s P vector definition, NC values below 12 are not included.

 

Figure 4 shows the effect of several different P vectors on GA performance as a function of chromosome size. As can be seen the different values of P generally have little influence on convergence. The only minor conclusion from Fig. 4 is that the curve generated when pB = 0.4 is slightly poorer for smaller values of NC.

 

Case 2: Airfoil Optimization

 

The GA is next applied to two transonic-flow airfoil problems—an inverse design problem utilizing the inviscid Euler equations and a constrained-lift drag-minimization problem utilizing the Navier-Stokes equations. The two-dimensional implicit flow solver (ARC2D) employing a central-difference/artificial-dissipation scheme is used to perform all function evaluations.22 Each airfoil geometry used in the present section is described using the PARSEC airfoil parameterization.23 The ten parameters used in this description, which are also the genes used by the GA optimization, are defined in Fig. 5. The parameter rle is the leading edge radius of curvature, xup is the x-coordinate location for which the upper surface thickness and curvature are yup and yxxup, respectively, (similar definitions apply for xlo, ylo and yxxlo on the airfoil lower surface), bte is the trailing edge included angle, yte is the trailing edge displacement away from the x axis, and ate is the angle between the trailing edge bisector and the x axis.

 

The first case is an inverse design problem where the computed  distribution for an NACA0012 airfoil (with closed trailing edge) at inviscid transonic flow conditions (M = 0.7, a = 2°) is used as a target solution.  Each gene is given a random initial value using Eq. (3) and the min/max values given in Table 1. The NACA 0012 airfoil is included within this design space. Note that the last two parameters, yte and ate, are constrained to be zero. This effectively reduces the number of genes to eight and helps the optimization recapture airfoil symmetry more efficiently. The fitness function is defined as

 

 

where the sum is over the  points which define the airfoil surface,  is the surface pressure coefficient distribution, and  is the target surface pressure coefficient distribution.

Fig. 5    Airfoil parameterization used for gene encoding.

 

 

Parameter

Min

Max

rle

0.002

0.015

xup

0.2

0.40

yup

0.03

0.09

yxxup

-0. 60

-0.20

bte

0.0

20.0

xlo

0.2

0.40

ylo

-0.09

-0.03

yxxlo

0.2

0.60

yte

0.0

0.0

ate

0.0

0.0

Table 1.  xmin/xmax values used for initialization and GA processing for the inverse airfoil design.

 

The GA was run with P = (0.1, 0.2, 0.4, 0.3),  = 0.3 and NC = 20 and 100.  Figure 6 shows selected  distributions at different points in the GA convergence history compared with the target distribution.  In particular, Fig. 6a (focusing on the upper surface shock location) shows several early time solutions, including two initial gene results (labeled: Generation 0) and the best result after two generations. At two generations the “best” current solution is still far from the target.

 

Figures 6b and 6c show later time designs where the GA has produced by generation 100 a good fitness and design.  Finally, in Fig. 6d, after 200 generations, the inverse airfoil design is fully converged and the design process has reproduced the target solution.  The last 100 generations may not be necessary in a real design process, as the iteration could be terminated earlier when the fitness reaches some prescribed tolerance.

Fig. 6    Initial and later generation pressure distributions compared to the target solution for the inverse design problem.

 

Figure 7a shows the GA fitness versus generation number, and Fig. 7b shows the GA fitness versus number of function evaluations (ARC2D calls) for the inverse design problem displayed in Fig. 6. These results represent typical GA convergence histories for an inverse design problem, e.g., Tse and Chan24 show inverse airfoil designs requiring 10,000 or more function evaluations. The 100-chromosome case (labeled Non-sym 100) converges faster in terms of generations, but is a factor of 3-4 slower in terms of function evaluations, (which is directly indicative of computer time).  High fitness for the 20-chromosome case (labeled Non-sym 20) is achieved after 5000 function evaluations. 

 

This performance can be improved by modifying the fitness measure or by including a symmetry constraint. Also shown in Fig. 7, is the GA convergence for such a case (labeled Sym 20). Imposition of symmetry reduces the airfoil definition to five genes (rle, xup, yup, yxxup and bte). As can be seen, the more restrictive symmetric design case converges more rapidly than the non-symmetric case.  Examination of the solution at various stages in the convergence process reveals that the non-symmetric case rapidly produces airfoils close to the target shape, but requires a large number of additional generations to converge to the symmetric solution.

Fig. 7    Fitness convergence for the inverse design problem.

 

The next problem is an airfoil design for viscous transonic flow (M = 0.7, a = 2°, Re = 9 x 106), which uses the Spalart-Almaras turbulence model. Twenty initial random chromosomes were generated using the min/max values shown in Table 2. The P vector and  values are the same as in the previous case.  The fitness function is given by

 

 

which produces a drag minimization optimization with nearly constant lift (). 

 

Figure 8 shows the pressure distributions (Fig. 8a), fitness convergence (Fig. 8b) and airfoil shapes (Fig. 8c) at various stages of the GA design process for the best chromosome. Clearly, one result of this optimization is reduction in the airfoil thickness, something that may not be desirable when structural considerations are included. However, the important thing to note is that the optimization performed as expected maintaining nearly constant lift and satisfying all problem constraints. The final airfoil represents a reasonable derivative from the class of airfoils possible within the given design space. The fitness function improved by 200% over the course of the optimization and converged in approximately 50 generations or 1000 function evaluations.

 

 

Parameter

Min

Max

rle

0.005

0.025

xup

0.25

0.50

yup

0.03

0.08

yxxup

-0. 50

-0.25

bte

5.0

20.0

xlo

0.25

0.50

ylo

-0.08

-0.03

yxxlo

0.25

0.50

yte

-0.02

0.02

ate

0

20.0

Table 2.  xmin/xmax values for initialization and GA processing of viscous transonic airfoil design.

 

Fig. 8 Viscous airfoil optimization, a) Selected pressure distributions, b) Fitness convergence history and c) Selected airfoil shapes.

 

The present GA can produce airfoil parameters that violate geometric constraints, such as crossing surfaces. In that event the fitness is set to zero for that chromosome, causing it to drop out in the next generation. The convergence tolerance of ARC2D was not varied during GA iteration, and in fact, only a small number of flow solver iterations (400) were performed for each chromosome. Therefore, it is possible that flow solver lack-of-convergence error for some chromosomes could compromise GA convergence by producing non-realistic results such as negative lift or drag. This problem was corrected by checking the values of lift and drag after each function evaluation. If either was negative, the fitness was set to zero. For the inverse design case approximately 10% of the function evaluations failed because of physical or geometric constraint violations.  For the viscous case, no failures occurred in the optimizations attempted thus far.  Further study of this aspect will be conducted in the future.

 

The airfoil design cases described above were run on a cluster of workstations (typically SGI workstations with a mix of R10000 (250MHz) and R12000 (400MHz) processors).  The GA approach is inherently embarrassingly parallel, and for these cases, the function evaluations (calls to ARC2D) were processed in parallel producing a linear speedup. Each function evaluation required approximately 30 processor-seconds of CPU time for the inviscid cases and 1-2 processor-minutes for the viscous cases (exact times depend on the processor speed and loading) and accessed a common data area for the chromosomes and solution restart files. The grid size was 193X33 for inviscid cases and 193x49 for viscous cases. A cluster of 30 machines was able to process on the order of 10,000 function evaluations in a wall clock time of approximately 6 hours for the inviscid case and 12 hours for the viscous case. 

 

Improvement in load balancing and hardware fault tolerance issues are required to improve overall parallel performance. Other strategies, such as non-synchronized GA processing (e.g., ignoring slower function evaluations) that produce a non-deterministic approach, micro-GA’s (smaller co-evolving GA’s) or hierarchal GA’s (using heterogeneous flow codes, i.e., full-potential—Euler—Navier-Stokes) are possible and would have an impact on the GA performance and parallelization strategy.

 

 

 

 

Case 3: Wing Optimization

 

The last case presented demonstrating the use of GA optimization involves a three-dimensional transonic wing. The baseline geometry used for this set of computations is the RAE Wing A25 mounted on a symmetry plane. This geometry involves a symmetric wing with a 9% maximum thickness-to-chord ratio, a taper ratio of 0.333, a leading edge sweep of 36.65° and an aspect ratio of 6.0. The freestream flow conditions used for this set of computations were held fixed at M = 0.84, a = 4°.

 

The TOPS26,27  CFD code is utilized for this case to perform all function evaluations.  This code is a three-dimensional full potential solver that utilizes a chimera zonal grid approach for handling complex geometries. It has the capability of producing complete numerical solutions about wing configurations (~115K grid points) in 1-2 min on a single SGI Origin 2000 processor.  Each run consists of surface and volume grid generation; chimera hole cutting, donor cell search, and interpolation coefficient computation; and, finally, the flow solver step. Each of these steps is automatically coupled and executed without user intervention, making GA/TOPS implementation easy and efficient.

 

Each grid used consists of two grid zones, an inner C-H topology grid fitted to the wing surface and an outer sheared-stretched Cartesian grid. The inner grid is generated using the HYPGEN grid generation program.28 During GA iteration another option for generating the inner grid is available that first reads in a baseline-HYPGEN grid and then modifies it according to the wing geometry perturbations that have just been computed via the GA process. This saves a small amount of computer time in that the grid does not have to be recomputed from scratch at the beginning of each function evaluation.

 

All flow computations performed during GA iteration used a coarse grid consisting of 115K points. This provides adequate accuracy for the optimization to proceed while allowing efficient code operation. Once the design optimization is complete, flow solutions for the initial and final geometries are recomputed using a finer grid (494K points), thus allowing a more accurate assessment of the performance improvement actually achieved.

 

The design space discretization for this wing optimization problem consists of ten genes. Eight of the genes are associated with the wing upper surface thickness—four thickness variables at two span stations each, wing root and tip. These thickness values are implemented using bump functions2 located at x/c = 0.15, 0.35, 0.45, 0.65. The value of wing twist at the root and tip are the two remaining genes. The maximum and minimum gene values used for all thickness functions are ±0.01. The wing root twist maximum and minimum values are 3° and 0°, respectively, and the wing tip twist values are, 0° and –3°, respectively. The simple definition of the design space described above is chosen because the emphasis in this study is on GA optimization—it’s implementation strategy and efficiency—not on aerodynamic efficiency.  All computations in this section have been performed on an SGI workstation with a single R10000 (250MHz) processor using Fortran 77.

 

The fitness function for all wing optimizations reported in this section is given by

 

 

In this equation CL and CD are the wing inviscid lift and drag coefficients, respectively. The 0.025 constant is added to the drag coefficient to approximate viscous drag. The 0.451 constant in the denominator’s second term is the coarse-grid baseline solution lift coefficient and provides a simple constraint on the optimization process keeping the lift approximately fixed. 

 

A typical result obtained using the GA optimization and design space discretization described above is presented in Figs. 9 and 10. This GA optimization was computed with NC = 20,  = 0.3 and P = (0.1, 0.2, 0.3, 0.4). Figure 9 shows pressure coefficient distributions at several span stations for both the baseline and optimized wing solutions. Figure 10 shows Mach number contours for the upper wing surface for the same two solutions presented in Fig. 9. Note that the optimization has produced a solution with significantly reduced shock strength, especially outboard of mid semi-span. In addition, the single shock characteristic of the baseline solution has been replaced with a lambda shock pattern in the optimized solution. The inviscid drag for this computation was reduced by 67 counts while the lift increased by only 0.0016.

The next several results focus on GA convergence performance. Figure 11 shows the effect of flow solver convergence level on GA convergence. There are three curves plotted in Fig. 11, each showing how the maximum fitness converges with the number of CFD flow solver function evaluations over 30 generations. The first curve shows GA convergence when each flow-solver solution maximum residual (RMAX) is reduced below 10-4. Since the initial maximum residual is on the order of 0.01, this roughly corresponds to a two order of magnitude reduction in maximum residual for each solution—a level of convergence that does not achieve plottable accuracy. As can be seen in Fig. 11, the GA process barely improves the maximum fitness above the baseline level.

 

 

Fig. 9    Pressure coefficient distributions at selected semi-span stations showing the baseline and optimized solutions, M = 0.84, a = 4°, NC = 20,  = 0.3, P = (0.1, 0.2, 0.3, 0.4), NG = 10.


a) Baseline solution

b) Optimized solution


 


Fig. 10 Mach number contours plotted on the upper wing surface for the baseline and optimized solutions, M = 0.84, a = 4°, NC = 20,  = 0.3, P = (0.1, 0.2, 0.3, 0.4), NG = 10.

 


The second curve in Fig. 11 displays GA convergence when each flow-solver solution maximum residual is reduced below 10-5. This roughly corresponds to a three order of magnitude reduction in maximum residual and produces plottable accuracy for all but the most difficult areas to converge, e.g., the solution around the shock may not be completely converged for this level of maximum residual reduction. The third curve in Fig. 11 shows GA convergence when each flow solver solution maximum residual is reduced below 10-6. This roughly corresponds to a four order of magnitude reduction in maximum residual and produces solid plottable accuracy over the entire solution. Note that for the last two curves displayed in Fig. 11, the GA convergence histories are nearly identical. However, the third curve, corresponding to the more tightly converged flow solutions, does produce a peak maximum fitness value that is 13% higher. It is interesting to note that the GA optimization performance displayed in Fig. 11 (second and third curves) is comparable to convergence of the finite-difference gradient method optimization reported in Ref. 29, which was used to optimize a similar problem.

 

Fig. 11 Genetic algorithm convergence history comparison showing the effect of flow solver convergence level on GA performance, M = 0.84, a = 4°, NC = 20,  = 0.3, P = (0.1, 0.2, 0.3, 0.4), NG = 10.

 

Another important way to view GA efficiency is to measure convergence against computer time and not against the number of function evaluations. The results presented in Fig. 11 are replotted and displayed versus CPU time in Fig. 12. Note that the level of flow solver solution convergence—as expected—has a profound effect on computational cost. The second curve, while producing 87% of the third curve’s maximum fitness, costs only about 40% as much in terms of computer time. 

 

In the hill climbing and airfoil problems the effect of NC on GA convergence was presented in Figs. 3 and 4 and Fig. 7, respectively. For the present wing optimization problem a similar comparison is presented in Fig. 13. The four curves correspond to GA optimizations that utilize 10, 20, 50 and 100 chromosomes, respectively. All computations utilized flow solutions converged to 10-6,  = 0.3 and P = (0.1, 0.2, 0.3, 0.4). The number of generations for each optimization was set so that each convergence history would consist of just under 2000 function evaluations. For this ten-gene aerodynamic shape optimization problem it is clear that the use of a smaller number of chromosomes (10-20) produces superior GA convergence. The 50-chromosome case is only slightly inferior, but the 100-chromosome case is significantly slower.

 

Fig. 12 Genetic algorithm convergence history comparisons showing the effect of flow solver convergence level on GA CPU time, M = 0.84, a = 4°, NC = 20,  = 0.3, P = (0.1, 0.2, 0.3, 0.4), NG = 10.

 

Fig. 13 Genetic algorithm convergence history comparison showing the effect of number of chromosomes on GA performance, M = 0.84, a = 4°,  = 0.3, P = (0.1, 0.2, 0.3, 0.4), RMAX ≤ 10-6, NG = 10.

For all computations presented in this section a special solution initialization feature has been used in the CFD flow solver. Instead of using a freestream initial solution, the solution is initialized from a file that contains the most recent maximum fitness solution. Whenever the GA procedure encounters a new chromosome having a fitness that exceeds the previous maximum, the flow solution is saved and used to initialize each succeeding solution until it is replaced with a superior solution. Because changes in many chromosomes are small—being perturbations away from the several fittest individuals—use of this philosophy greatly reduces the amount of computer time required for a GA computation. A quantification of this improvement is presented in Fig. 14. The first curve shows a GA convergence history that uses the solution restart option, and the second curve shows an identical GA convergence history that does not use solution restart. Each flow solution for both curves is converged to RMAX ≤ 10-6. For this optimization the use of restart saves a factor of 2.7 in computer time.

 

The successful use of solution restart and less well-converged solutions during GA iteration is a testament to the robustness of the GA approach for optimization. Both these shortcuts introduce error into the GA process, but, if suitably controlled, do not hamper overall GA convergence.

 

Fig. 14 Genetic algorithm convergence history comparisons showing the effect of solution restarts on GA CPU time, M = 0.84, a = 4°, NC = 20,  = 0.3, P = (0.1, 0.2, 0.3, 0.4), RMAX ≤ 10-6, NG = 10.

Conclusions

 

A genetic algorithm (GA) procedure suitable for performing optimizations, including aerodynamic shape optimizations, has been presented. It uses real-number encoding to represent all geometric parameters as genes in the GA. Four operators are utilized to advance from one generation to the next. They include passthrough, random average crossover, perturbation mutation and mutation.

 

The GA performed surprisingly well on all cases presented, demonstrating wide applicability. Essentially the same GA code was used for all problems presented herein—hill climbing, inviscid inverse design for airfoils, viscous L/D maximization for airfoils and inviscid L/D maximization for wings.

 

GA/flow solver coupling was quite easy to setup and required only a few hours for the simplest cases. A wide variety of optimization problems could theoretically be solved using the same GA with a minimal implementation effort providing an appropriate function evaluation mechanism was available.

 

For the problems presented herein, which were limited to a maximum of ten genes, the optimal number of chromosomes was small—on the order of 20. Reasonable levels of convergence were achieved for most cases in several hundred function evaluations with the slowest cases requiring several thousand. Utilization of solution restarts and reduced levels of convergence per function evaluation provided a total reduction in computer time of six to eight.

 

The results indicate that the GA approach is robust. It continues to converge even with significant amounts of lack-of-convergence error in each function evaluation or when as many as ten per cent of all function evaluations failed due to constraint violation or solution divergence. (Higher failure rates, up to 50%, probably would not significantly affect GA convergence, but have not been encountered).

 

Despite good performance from the GA, it is recognized that gradient methods would probably provide solutions to the single-objective, smooth problems described herein in less computer time. The strength of genetic algorithms for optimization lies in their ability to provide solutions for non-smooth design spaces and for multi-objective problems. This will be the focus of future work.

 

References

 

1.     Reuther, J., “Aerodynamic Shape Optimization Using Control Theory,” RIACS Rep. 96.09, 1996. 

 

2.     Hicks, R. and Henne, P., “Wing Design by Numerical Optimization,” J. Aircraft, Vol. 15, 1978, pp. 407-412.

 

3.          Bischof, C., Carle, A., Corliss, G., Griewank, A. and Hovland, P., “ADIFOR—Generating Derivative Codes from Fortran Programs,” Scientific Programming, No. 1, 1992, pp. 1-29.

 

4.          El-banna, H. M. and Carlson, L. A., “Determination of Aerodynamic Sensitivity Coefficients Based on the Transonic Small Perturbation Formulation,” AIAA J., Vol. 27, 1990, pp. 507-515.

 

5.          El-banna, H. M. and Carlson, L. A., “Aerodynamic Sensitivity Coefficients Using the Three-Dimensional Full Potential Equation,” J. Aircraft, Vol. 31, 1994, pp. 1071-1077.

 

6.          Arslan, A. E. and Carlson, L. A., “Integrated Determination of Sensitivity Derivatives for an Aeroelastic Transonic Wing,” J. of Aircraft, Vol. 33, 1996, pp. 224-231.

 

7.          Jameson, A., “Aerodynamic Design via Control Theory,” J. of Sci. Comp., Vol. 3, 1988, pp. 233-260.

 

8.          Jameson, A., “Automatic Design of Transonic Airfoils to Reduce the Shock Induced Pressure Drag,” 31st Israel Annual Conf. on Aviation and Aeronautics, Tel Aviv, 1990, pp. 5-17.

 

9.          Goldberg, D. E., “Genetic Algorithms in Search, Optimization and Machine Learning,” Addison-Wesley, Reading, MA, 59-88, 1989.

 

10.      Davis, L., “Handbook of Genetic Algorithms,” Van Nostrand Reinhold, New York, 1991.

 

11.      Beasley, D., Bull, D. R. and Martin, R. R., “An Overview of Genetic Algorithms: Part 1, Fundamentals,” University Computing, Vol. 15, No. 2., 1993, pp. 58-69.

 

12.      Obayashi, S. and Tsukahara, T., “Comparison of Optimization Algorithms for Aerodynamic Shape Design,” AIAA J., Vol. 35, 1997, pp. 1413-1415.

 

13.      Bock, K.-W., “Aerodynamic Design by Optimization,” Paper 20, AGARD CP-463, 1990.

 

14.   Quagliarella, D. and Della Cioppa, A., “Genetic Algorithms Applied to the Aerodynamic Design of Transonic Airfoils,” AIAA Paper 94-1896-CP, 1994.

 

15.   Vicini, A. and Quagliarella, D., “Inverse and Direct Airfoil Design Using a Multiobjective Genetic Algorithm,” AIAA J., Vol. 35, 1997, pp. 1499-1505.

 

16.   Obayashi, S., Yamaguchi, Y. and Nakamura, T., “Multiobjective Genetic Algorithm for Multidisciplinary Design of Transonic Wing Planform,” J. of Aircraft, Vol. 34, 1997, pp. 690-693.

 

17.   Sasaki, D., Obayashi, S., Sawada, K. and Himeno, R., “Multiobjective Aerodynamic Optimization of Supersonic Wings Using Navier-Stokes Equations,” European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS 2000, Barcelona, Spain, Sept. 11-14, 2000.

 

18.   Oyama, A., “Multidisciplinary Optimization of Transonic Wing Design Based on Evolutionary Algorithms Coupled with CFD Solver,” European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS 2000, Barcelona, Spain, Sept. 11-14, 2000.

 

 

 

19.   Oyama, A., “Wing Design Using Evolutionary Algorithms,” PhD Thesis, Dept. of Aeronautics and Space Engineering, Tohoku University, Senadi, Japan, March 2000.

 

20.   Houck, G. R., Joines, J. A. and Kay, M. G., “A Genetic Algorithm for Function Optimization: A Matlab Implementation,” North Carolina State University-IE, TR 95-09, 1995.

 

21.   Michalewicz, Z., Genetic Algorithms + Data Structures = Evolution Programs, AI Series, Springer-Verlag, New York, 1994.

 

22.   Pulliam, T. H., Implicit Finite-Difference Methods for the Euler Equations, Advances in Computational Transonics, Fourth Volume in the Series on Recent Advances in Numerical Methods in Fluids, Editor W. G. Habashi, Pineridge Press Ltd., 1983.

 

23.      Sobieczky, H., “Parametric Airfoils and Wings,” Recent Development of Aerodynamic Design Methodologies-Inverse Design and Optimization, Friedr. Vieweg & Sohn Verlagsgesellschaft mbH, Braunschweig/Wiesbaden, Germany, 1999, pp. 72-74.

 

24.      Tse, D. and Chan, L., “Transonic Airfoil Design Optimization using Soft Computing Methods,” Canadian Aeronautics and Space J., Vol. 46, No. 2, June 2000, pp. 65-73.

 

25.      Treadgold, A. F., Jones, A. F. and Wilson, K. H., “Pressure Distribution Measured in the RAE 8ft X 6ft Transonic Wind Tunnel on RAE Wing ‘A’ in Combination with an Axi-Symmetric Body at Mach Numbers of 0.4, 0.8 and 0.9,” Experimental Data Base for Computer Program Assessment, AGARD-AR-138, May 1979.

 

26.      Holst, T. L., “Full Potential Equation Solutions Using a Chimera Grid Approach,” AIAA Paper No. 96-2423, June, 1996.

 

27.      Holst, T. L., “Multizone Chimera Algorithm for Solving the Full-Potential Equation,” J. of Aircraft, Vol. 35, No. 3, May-June 1998, pp. 412-421.

 

28.      Chan, W. M., Chiu, I., and Buning, P. G., “User’s Manual for the HYPGEN Hyperbolic Grid Generator and the HGUI Graphical User Interface,” NASA TM 108791, Oct. 1993.

 

29.   Cheung, S. and Holst, T., “Aerodynamic Shape Optimization Using a Combined Distributed/Shared Memory Paradigm,” Proceeding of the Fourth CAS Workshop, Moffett Field, Calif., 1998.

 

 



http://www.mathworks.com