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
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.
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
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.
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
.
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.
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.
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.
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.
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.
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.
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.
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.