\ '~vnuh
- nu.

 

 

   

v‘v

HIGAN STAT

lllllll lllllllllllllllllllll

31293 01716 4223

llllllllll

 

 

 

 

 

 

 

 

                           

This is to certify that the
thesis entitled
Use of Injection Island Genetic Algorithms in the

Optimization of Composite Flywheels

presented by

David James Eby

has been accepted towards fulfillment
of the requirements for

Master's Mechanics

 

 

degree in

Mam

Major professor

Date I 0/3/97

0.7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

 

 

LIBRARY

Michigan State

University

 

 

PLACE IN RETURN BOX
to remove this checkout from your record.
TO AVOID FINES return on or before date due.

 

ME DUE

DATE DUE

DATE DUE

 

new area

 

 

 

“JV" "

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1M mus-p14

USE OF INJECTION ISLAND GENETIC
ALGORITHMS IN THE
OPTIMIZATION OF COMPOSITE FLYWHEELS

By
David J. Eby

A THESIS

Submitted to
Michigan State University
in the partial fullfillment of the requirements

for the degree of

MASTERS OF SCIENCE
Department of Materials Science and Mechanics

1997

ABSTRACT

USE OF ISLAND INJECTION GENETIC
ALGORITHMS IN THE
OPTIMIZATION OF COMPOSITE FLYWHEELS
By
David J. Eby

This thesis presents an approach to optimal design of elastic flywheels using an Island
Injection Genetic Algorithm (iiGA). An iiGA in combination with a structural finite ele-
ment code is used to search for shape variations to optimize the Specific Energy Density
(SED) of elastic flywheels while controlling angular velocity and penalizing designs that
are deemed unsafe. SED is defined as the amount of rotational energy stored per unit mass.
iiGAs seek solutions simultaneously at different levels of refinement of the problem repre-
sentation (and correspondingly different definitions of the fitness function) in separate sub-
populations (islands). Solutions are sought first at low levels of refinement with an axi-
symmetric plane stress finite element code for high speed exploration of the coarse design
space. Next, individuals are injected into populations with a higher level of resolution that
use an axi-symmetric three dimensional finite element code to “fine-tune” the structures.
Discretizing the fitness function into “sub-fitness” functions while using a computationally
inexpensive axi-symmetric plane stress finite element code allows efficient and robust ex-
ploration. The problem at hand is a combinatorial and parametric search. Since Genetic Al-
gorithms (GAS) are particularly efficient at combinatorial optimization but lack in the area
of hill climbing, hybrid iiGAs that combine Threshold Accepting (TA) and various local

optimization procedures are developed to increase the robustness of the GA.

ACKNOWLEDGMENTS

I wish to express my gratitude to Dr. Erik Goodman and Dr. William F. Punch III
for all the support, insight guidance and expertise extended throughout my research. I
would like also thank the Genetic Algorithms Research Group for their patience and advice.
Particularly Owen Mathews, Boris Gelfand and Victor Miagkikh for their hard work and
assistance during this research.

Also, I would like to thank the Computational Mechanics Research Group for their
support and assistance. Finally, I would like to thank my advisor Dr. Ronald Averill for his
incredible amount of guidance, understanding and foresight throughout the duration of my

thesis.

Table of Contents

List of Tables ............................................................................................................ iv

List of Figures .......................................................................................................... v

Chapter 1

Introduction ................................................................................................................ 1
1.1 Introduction ................................................................................................. 1
1.2 Problem Definition ...................................................................................... 2
1.3 Literature Review ........................................................................................ 5
1.4 Present Study ............................................................................................... 29

Chapter 2

Algorithm Design and Comparison ...................................................................... 30
2.1 Introduction ................................................................................................. 30
2.2 Genetic Algorithms ..................................................................................... 31
2.3 Simulated Annealing and Threshold Accepting ......................................... 41

Chapter 3

Finite Element Models of Flywheels ................................................................... 44
3.1 Introduction .................................................................................................. 44
3.2 Finite Element Model Formulations ........................................................... 44
3.3 Verification of Results .................................................................................. 51

Chapter 4

Optimal Design of Flywheels Using Genetic Algorithms .............................. 57
4.1 Introduction .................................................................................................. 57
4.2 Description of Optimization Problem .......................................................... 57
4.3 Model Description ...................................................................................... 60
4.4 Island Injection Genetic Algorithm Approach ............................................. 63
4.5 Results and Discussion ............................................................................... 70

Chapter 5

Summary and Conclusions .................................................................................... 117

References .................................................................................................................. 120

iii

List of Tables

Table 4.2. 1. Weighting Coefficient Values ..................................................................... 59
Table 4.3.1. Isotropic Material Properties ...................................................................... 62
Table 4.3.2. Composite Material Properties ................................................................... 62

Table 4.5.1 Hybrid Island Injection Genetic Algorithm Para-
meters for the Enumerable Coarse Design Space of a Solid
Isotropic Flywheel ......................................................................................................... 76

Table 4.5.2. Statistics on the Total Number of Evaluation Calls
made for Various Optimization Approaches in Average 2-Hour

Run ............................................................................................................................ 80
Table 4.5.3 Island Injection Genetic Algorithm Parameters .......................................... 88
Table 4.5.4 Genetic Algorithm Topological “Ring” Parameters .................................... 88

Table 4.5.5. Statistics on the Total Number of Evaluation Calls
for iiGA, Hybrid iiGA and Topological “Ring” GA Optimization
Approaches .................................................................................................................... 9O

List of Figures

Figure 1.1. The “knees” of a typical Pareto front .......................................................... 8

Figure 2.1. Simple GA structure .................................................................................... 34
Figure 2.2. Possible iiGA topology that searches at various levels

of geometric resolution .................................................................................................. 40
Figure 3.1. Flywheel geometry ...................................................................................... 53
Figure 3.2 Element topology .......................................................................................... 53

Figure 3.3. Typical axi-symmetric planar finite element mesh with
symmetry boundary conditions (1,134 degrees of freedom) ......................................... 53

Figure 3.4. Tangential and radial stresses as a function of radius

predicted by plane stress finite element code and analytical

solution for a flat solid flywheel. Arrow represents cross section

analyzed ......................................................................................................................... 54

Figure 3.5. Tangential and Radial Stresses as a Function of Radius

Predicted by Plane Stress Finite Element Code and Analytical

Solution for a Flat Annular Flywheel. Arrow represents cross

section analyzed ............................................................................................................. 54

Figure 3.6. Variation of radial stress through the thickness for the
MARC and present finite element code for a solid isotropic flywheel
that varies in thickness. Arrow represents cross section plotted. ................................. 54

Figure 3.7. Variation of tangential stress through the thickness

for the MARC and present finite element code for a solid

isotropic flywheel that varies in thickness. Arrow represents cross

section plotted ................................................................................................................ 54

Figure 3.8. Variation of transverse stress through the thickness

for the MARC and present finite element code for a solid isotropic

flywheel that varies in thickness. Arrow represents cross section

plotted ............................................................................................................................ 56

Figure 3.9. Variation of shear stress through the thickness for the
MARC and present finite element code for a solid isotropic flywheel

that varies in thickness. Arrow represents cross section plotted ................................... 56
Figure 4.1. Visual display of flywheel ........................................................................... 59
Figure 4.2. Plane stress finite element representation ................................................... 61

Figure 4.3. Typical planar finite element mesh with symmetry
boundary conditions used in the three-dimensional fem evaluations
(1,134 degrees of freedom) ............................................................................................ 61

Figure 4.4. iiGA topology that performs multiple refinements in
geometric resolution ....................................................................................................... 65

Figure 4.5. Island injection GA topology with various levels of
fitness evaluation efficiency ........................................................................................... 67

Figure 4.6. Island injection genetic algorithm topology with
multiple fitness definitions ............................................................................................. 68

Figure 4.7 Typical hybrid island injection genetic algorithm toplogy ........................... 69

Figure 4.8. Simple genetic algorithm exploitation of plane stress
evaluation ....................................................................................................................... 94

figure 4.9. Theoretical constant stress profile compared to simple
genetic algorithm constant stress profile ........................................................................ 94

Figure 4.10. The simple genetic algorithm evolution of a constant

stress flywheel ................................................................................................................ 95
Figure 4.11. Simple genetic algorithm non-constant stress profile

compared to the theoretical constant stress profile ........................................................ 96
Figure 4.12. Simple genetic algorithm compared to ta algorithm

and hybrid genetic algorithms for a solid isotropic flywheel in

an enumerable search space ........................................................................................... 96
Figure 4.13. Island injection GA topology ..................................................................... 97
Figure 4.14. Single computational node island injection genetic

algorithm results for a solid isotropic flywheel in an enumerable
search space ................................................................................................................... 97

vi

Figure 4.15. Hybrid island injection GA topology (TA algorithm) ............................... 98

Figure 4.16. Single computational node hybrid island injection
GA (TA algorithm) results for a solid isotropic flywheel in an
enumerable search space ................................................................................................ 98

Figure 4.17. Hybrid island injection GA topology (local search
method) .......................................................................................................................... 99

Figure 4.18. Single computational node hybrid island injection
GA (local search method) results for a solid isotropic flywheel in
an enumerable search space ........................................................................................... 99

Figure 4.19 Hybrid island injection GA topology (local search
and TA) .......................................................................................................................... 100

Figure 4.20. Single computational node hybrid island injection.
GA (local search and TA) results for a solid isotropic flywheel in
an enumerable search space ........................................................................................... 100

Figure 4.21. Twelve “ring” GA topology ...................................................................... 101

Figure 4.22. Single computational node topological “ring” GA
results for a solid isotropic flywheel in an enumerable search
space ............................................................................................................................ 101

Figure 4.23. Demonstration of a flywheel that is an “artifact”
for a three-dimensional finite element model (containing 2,482
degrees of freedom) ....................................................................................................... 102

Figure 4.24. Hybrid iiGA “raw” fitness for each island as a

function of time for a solid multi-material composite flywheel.

Islands 0-2: plane stress FEM evaluation. Islands 3-11: 3-d FEM

evaluation ....................................................................................................................... l 02

Figure 4.25. Hybrid iiGA fitness re-evaluated for each island as
a function of time. All re-evaluations performed with 3-d FEM
containing 12,272 degrees of freedom ........................................................................... 103

Figure 4.26. Hybrid iiGA quickly finding building blocks at low
levels of geometric resolution to inject into islands of higher
resolution ........................................................................................................................ 104

Figure 4.27. Best flywheel found at each level of resolution and
fitness definition ............................ . ................................................................................. 105

vii

Figure 4.28. iiGA fitness re-evaluated for each island as a

function of time. All re-evaluations performed with 3-d FEM

containing 12,272 degrees of freedom for an annular multi-material

composite flywheel ........................................................................................................ 106

Figure 4.29. Hybrid iiGA fitness re—evaluated for each island as

a function of time. All re-evaluations performed with 3-d FEM

containing 12,272 degrees of freedom for an annular multi-

material composite flywheel .......................................................................................... 107

Figure 4.30. iiGA quickly finding building blocks at a low level
of resolution to inject into islands of higher resolution for an
annular multi—material composite flywheel ................................................................... 108

Figure 4.31. Hybrid iiGA quickly finding building blocks at a low
level of resolution to inject into islands of higher resolution for an
annular multi-material composite flywheel ................................................................... 109

Figure 4.32. Twenty island “ring” GA topology ............................................................ 110

Figure 4.33. “Ring” topology GA results for an annular multi-
material composite flywheel .......................................................................................... 110

Figure 4.34. Comparison of iiGA to hybrid iiGA to “ring” topology

GA. For this problem, both the iiGA and hybrid iiGA search

extremely efficiently when compared to a “ring” topology GA.

The hybrid iiGA discovers solutions with higher fitness values

when compared to the “ring” topology GA and the iiGA alone .................................... 111

Figure 4.35. Best flywheel found by “ring” topology GA,
iiGA and hybrid iiGA for an annular multi-material composite
flywheel .......................................................................................................................... 1 12

Figure 4.36. Hybrid island injection GA topology used in the
constrained optimization of a multi-material annular flywheel ..................................... 1 13

Figure 4.37. Hybrid iiGA fitness re-evaluated for each island as

a function of time. All re-evaluations performed with 3-d FEM

containing 12,272 degrees of freedom for the constrained

optimization of an annular composite multi-material flywheel ..................................... 114

Figure 4.38. Hybrid iiGA quickly finding building blocks at a

low level of resolution to inject into islands of higher resolution

for the constrained optimization of an annular multi-material

composite flywheel ........................................................................................................ 115

viii

Figure 4.39. Best flywheel found at each level of resolution and
fitness definition for the constrained optimization of a multi-
material composite flywheel .......................................................................................... 116

Chapter 1

Introduction

1.1 Introduction

Everyday a new optimization problem is encountered. For instance, what is the
quickest path to work? Where and how congested is the road construction? Am I better off
just riding my bike? If so, what is the shortest distance? This simple problem is easily
solved, while many engineering optimization problems cannot be as cleanly approached.
Engineering involves a wide class of problems and optimization techniques. Many engi-
neering design approaches such as the “make-it-and-break-it” mentality are simple, out-of-
date concepts that have been revolutionized with computer simulations that exploit various
mathematical concepts such as the finite element method to avoid costly design iterations.
Even with the aid of high speed supercomputers, this design process can be hindersome,
producing designs that evolve extremely slowly over a long period of time. The next step
in the engineering of systems is the automation of optimization through computer simula-
tion. After all, optimization is simply engineering on a grander scale.

Optimization approaches include analytical hill climbing, random search, directed
random search and hybrid methods. Hill climbing or gradient-based methods are single
point search methods that are severely restricted in their application due their numerous
limitations. Random search methods are limited to a class of practical problems that have

small search spaces. A directed random search method such as a Genetic Algorithm (GA)

is a multiple point search method that can be an effective optimization approach to a broad
class of problems. Injection Island Genetic Algorithms (iiGAs) can help reduce the com-
putational intensity associated with typical GAs by searching at various levels of resolution
within the search space. iiGAs can also incorporate various optimization schemes to create
hybrid methods.

An optimization problem can have single or multiple goals. Multi-objective opti-
mization first arose in a natural fashion in economics. Typically, most engineering optimi-
zation problems contain multiple objectives such as: minimizing cost while maximizing
safety, manufacturability and or performance. Of all the applications of optimization in en-

gineering, structural optimization using GAs is the focus of this study.

1.2 Problem Definition

1.2.1 Flywheels

Shape optimization of flywheels for the maximization of SED is an appealing
thought that has received its fair attention by researchers. The concept of a flywheel is as
old as the axe grinder’s wheel, but could very well hold the key to tomorrow’s problems of
efficient energy storage. The flywheel has a bright outlook because of the recent achieve-
ment of high specific energy densities with composite materials. Composite materials have
high specific strength, which can be seen to govern the threshold of specific energy density
in the equation:

SED ax = K.°“”

m

 

(1)

Cult is the ultimate strength of the material and p is the density of the material. Ks is a

non-dimensional shape factor with minimum value of zero and a maximum value of unity
that occurs as the radius tends to infinity.

The overall objective in the current study was to maximize specific energy density
(SED) of flywheels, which is defined as:

11002
SED = -— (1.2)
2 m

where (DIS angular velocity, m is the total mass and I is the mass moment of inertia of the
flywheel defined by:

I = iprzdv (1.3)
V

where p is the density of the material.

A simple example of a flywheel is a solid, flat rotating disk. The SED of a flat solid
disk can be increased by varying the shape of the disk to redistribute the inertial forces in-
duced from rotation. This thesis presents an approach to optimal design of elastic flywheels
using Genetic Algorithms (GAS). A GA in combination with a structural finite element
analysis code is used to search for shape variations in the optimization of composite fly-
wheels. The variables to be optimized include: Specific Energy Density (SED) and “air
gap” growth. Failure at this angular velocity would be catastrophic. Excessive “air gap”
growth occurs when the inner most ring of an annular flywheel expands beyond acceptable
tolerances of flywheel bearings due to forces induced from rotations. It would be ideal to
maximize SED, and to reduce “air gap” growth while maintaining a reasonable angular ve-
locity of a flywheel, by means of shape optimization. A GA represents a robust, efficient

optimization technique used to solve this problem.

A flywheel stores kinetic energy by rotating a mass about a constant axis, which
makes it easy to integrate flywheels into energy conservation systems. Systems such as ve-
hicles currently use flywheels during braking for regenerating energy lost during decelera-
tion. Another practical application is energy storage in low satellite orbits where
photoelectric cells are exposed to 60 minutes of light to charge, followed by 30 minutes of
darkness where stored energy must be used. Electrochemical energy storage (e. g., in bat-
teries) is limited by low cyclic lifetimes, low longtime reliability and low specific energies,
all major concerns in satellite applications. The flywheel is well-suited for this application
due to high cyclic lifetimes, longtime reliability and high specific energies. Also, large
scale flywheels could be used in energy plants to store huge amounts of energy. Finding
practical applications for flywheels is not the problem, but optimizing the SED of fly—

wheels, given a problem specific set of parameters and constraints, provides a challenge.

1.2.2 Optimization Approach

This thesis presents and compares a set of various optimization approaches applied
to flywheels. This set includes simple Genetic Algorithms (sGA), GAs with topological
“ring” structure, Island Injection Genetic Algorithms (iiGAs), Threshold Accepting (TA)
algorithms and hybrid GAs. Hybrid GAs contain local improvement methods and/or TA al-
gorithms to help improve the performance of a typical GA.

Unlike traditional optimization techniques, GAs do not require the gradient of the
(objective function to exist. GAs only need the definition of the objective, or fitness. This
makes GAs applicable to a large set of problems. A GA is a multi-point search procedure

modeled specifically on the mechanics of natural selection. A 36A is multi-point search

procedure that occurs within isolated surroundings, interacting with nothing else. Parallel
GAs that use for an example, a topological “ring’ structure constitute a set of multi-point
searches. Structured migration allows the exchange of search points among the set. An
iiGA also uses structured migration of search points among a set, but allows exploration of
various levels of resolution of the search space. iiGAs search various levels of resolution
to quickly find building blocks to inject into higher levels of resolution. iiGAs can also use
various fitness evaluation tools and fitness definitions. TA algorithms are a single point hill
climbing search procedure with a mechanism for climbing out of local optima (based on
statistical mechanics of annealing). A local improvement method is also sometimes incor-
porated into the GA, allowing more rapid search of local areas to be explored. Hybrid GAs

combine a mixture of any of the GA approaches with local and/or TA algorithms.

1.3 Literature Review
The following is a summary of literature reviewed-for the current research. The ar-
ticles are reviewed in categories. The first is structural optimization in general, the second

is GA optimization applications and the last is the optimization of flywheels.

1.3.1 Structural Optimization

According to Jenkins [25] the process of structural optimization is an attractive goal
to a designer because there is a satisfaction in producing the “best” design for a stated ob-
jective. This is, of course, only a fraction of the motivation behind structural optimization.
Olhoff and Taylor [34] state that considerations such as limited energy resources, shortage

of economic and some material resources, strong technological competition, and environ-

mental problems are the motivations behind structural optimization. Structural optimiza-
tion often contains multiple objectives such as minimizing cost while maximizing quality.
Stadler and Dauler [49] present a broad overview and history of multicriterion optimiza-
tion. Horn [24] argues that most, if not all optimization problems have multiple objec-
tives. Both Horn [24] and Olhoff and Taylor [49] claim that the entire “trade off” set, or
Pareto front, should be sought, but they vary drastically in their approach to this Pareto
front. Olhoff and Taylor [49] make multicriterion decisions first with multiple runs using
the classical variations of calculus approach while Horn [24] claims that the Pareto front
should be found simultaneously with a Pareto “Niched” GA. Belegundu et al. [1], Sundare-
san, et a1. [50], and Khot, et a]. [26] present an alternate point of view for optimization, that
is robust design through minimum sensitivity. The goal of the optimization is to design en-
gineering systems so that they are desensitized to variations in uncertainty variables, such
as manufacturing errors. Osyczka, Kuchta and Czula [35] present an approach to optimi-
zation of multicriterion systems with computationally expensive evaluations. Jenkins [25]
presents a general paper on structural optimization using GAs while Haslinger and Jedel-
sky [22] concentrate on shape optimization via GAS.

Stadler and Dauler [49] state that multicriterion optimization is rooted in econom-
ics, and was first treated in the bobk The Wealth of Nations (1776) by Smith. The first ref-
erence to multicriterion optimization was given by Edgeworth’s book Mathematical
Physics [24]. Horn [24] states Vilfred Pareto (1896) noted that a partial ordering of solu-
tions exists before any multicriterion decisions are drawn defined by the Pareto criterion.
The Pareto criterion states that a solution is superior if it is at least as good as other solutions

in some criteria and better than other solutions in at least one criterion. More rigidly, Horn

[24] states that, given k attributes, all of which are to maximized, a solution A, with at-
tributes (a0,a1,a2,...,ak,1), and a solution B with attributes (b0,b1,b2,...,bk_1), we say that A
dominates B if and only if for all i: 3i 2 bi, and there exists a j such that aj>bj. Some of the
pairs of the solution set are not comparable in the sense that neither solution dominates the
other and are termed incompatible. But conversely there are solutions that are dominated
and can be discarded before any multicriterion decisions are made. The set of solutions that
are not dominated by any other solution constitutes the Pareto front. This is demonstrated
in the two dimensional (k=2) graph in Figure 1.1, where the objectives are to minimize cost
and maximize quality. This simple example can be extended to include more complex mul-
ticriterion objective functions (k>2) by simply increasing the dimension of the plot. Figure
1.3 also demonstrates why the Pareto front could be sought. At points E and F in Figure
1.3, the “knees” of the front exist, which represent points at which small sacrifices in one
objective causes large sacrifices in the other. In this case a rational decision maker needs to
be presented with this information for analysis.

Olhoff and Taylor [34] discuss the fundamentals of structural optimization based
on the mathematical properties of the governing equations for optimal design. They started
by reviewing the basic concepts and well defined areas of optimization. This includes the
conventional design method in which changes in the design are followed by an analysis of
the design. Typically, this procedure is continued until a satisfactory design in cost and/or

response is developed. This procedure is largely based on educated guess work by the de-

signer.

 

 

 

O
9 0
0 O
O O
O
0 o
7 0
BO

6 0 O .

t O O
O
' 51 o O
a o , Knee
O
4 0
O
3 O
1 . F o
. 0 Pareto Front
2 ‘1 o 0
O
O
1 . O
o Knee
0 fi 4 1 v r g 1
O 1 2 3 5 O 7 0 9 fl
Goat

Figure 1.1. The “knees” of a typical Pareto front.

Olhoff and Taylor [34] labeled structural optimization as a design problem in which
preassigned parameters and design variables are subdivided with an objective function to
be maximized or minimized with respect to the design variables. Typical design variables
were categorized, such as: cross-sectional, layout, material, and loading. These design
variables are roughly continuous or discrete according to the manner in which they are re-
lated to the spatial variables. Cross-sectional design variables such as area and moments
of area can be viewed as discrete or continuous, depending on what cross-sectional proper-
ties are manufacturable or available. Layout of a structure can be divided into topological

and configurational design variables. Topological design variables account for the defining

connectivity of the members in the structure and are integer variables, while configuration-
al design variables describe surfaces and or centerlines of bars and are usually considered
continuous variables. Material design variables are usually discrete, representing choices
among alternate material properties. Loading design variables can be continuous or dis-
crete, representing conditions placed on a structure.

Olhoff and Taylor [34] state that constraints are any set of design variables that can
be represented by a point in the design space, which help reduce the total number of possi-
ble designs. Geometrical constraints such as maximum and minimum thickness are im-
posed explicitly on the design variables. Behavioral constraints are broken down into
equality and inequality constraints. Equality constraints are stated as the equations govem-
ing structural response. Inequality constraints are imposed locally, such as local allowable
stress, or globally, such as compliance. Both equality and inequality constraints are often
nonlinear and implicit with respect to design variables.

The objective function must be expressed in terms of the design variables such that
its value is definable anywhere in the design space. This objective function is to be maxi-
mized or minimized with respect to the set of design variables within the design space. A
problem with a single objective function defines a single criterion optimization problem.
Olhoff and Taylor [34] present the classical calculus of variations for single criterion opti-
mization problems. A single criteria optimization problem can be rigorously stated as:

Determine Di, i=1,...,n which

minimize F(D],...,Dn) = 0,

subject to hj(D,,...,Dn) = 0, j=1,...,p

and gk(D1,...,Dn) SO, k=1,...,q.

10

Where Di is a vector of design variables, F is the objective function which is to be
minimized with respect to the design vector under p equality and q inequality constraints.
Assuming that the function F is differentiable with respect to the design variables, a calcu-
lus of variations approach can be defined. There are numerous well established single cri-
terion optimization problems in the area of static loading, including design against plastic
collapse, elastic response, optimal layout of trusses and dynamic loading under frequency
constraints using calculus of variations, Olhoff and Taylor [34].

Granhi [21] presents an overview of optimization of structures using sensitivity de-
rivatives calculation for repeated eigenvalues for frequency constraint approximations.
Various optimization algorithms and applications are also reviewed.

Horn [24] argues that most, if not all, optimization problems have multiple objec-
tives. Horn [24], Olhoff and Taylor [34] present approaches to multiple objective optimi-
zation problems in different ways. Both papers agree that the Pareto front should be sought,
but they approach it differently.

Horn [24] classifies three approaches according to how they handle the problems of
searching for the Pareto front and when the multicriterion decision is made. The first ap-
proach makes a multicriterion decision first and then performs a search. The second search-
es first, and then makes a multicriterion decision using a rational decision maker. The third
integrates search and multicriterion decision making.

When a multicriterion decision is to be made before searching, the objective func-
tion is aggregated into a single objective function. A scalar-aggregate or an order—aggre-
gate can define a multicriterion objective function. The simplest scalar aggregate is a

weighted sum of the single objective functions. The decision maker tries weighting the ob-

ll

jectives according to their merit. One drawback of such an approach is that it only accounts
for linear variations among the criteria. This weighted sum can be generalized into any
non-linear relationship, but weights are still at best guided by intuitional guesswork. This
approach can be, in general, applied to any type of optimization algorithm to find the Pareto
front by running multiple optimization runs, varying each weight from 0 to a maximum val-
ue while holding the other weights constant; but a weighted sum cannot capture concave
portions of the Pareto optimal surface. This due to the fact that linear aggregative methods
are biased towards convex portions of the Pareto front [24]. A lexicographic approach is a
non-scalar aggregate approach to multicriterion objective problems, Horn [24]. This ap-
proach requires the decision maker to order the criteria and rank the solutions best to worst
without assigning any scalar values. The approach can be tackled with evolutionary com-
putation such as GAS by using rank-based and toumament-style selection methods.

An aggregated approach is, at best, an overly simplistic view of multicriterion op-
timization problems. Is it possible to combine objectives that contrast each other into a sin-
gle objective? Horn [24] recognizes this difficulty, if not impossibility, of making
multicriterion decisions before searching. Many dominated solutions can be eliminated by
searching first and isolating the Pareto solutions in the current set defined by the Pareto cri-
terion. Horn [24] obviously limits this application to a multi-point search algorithm such
as GAS. Integrated search and decision making requires a preliminary multicriterion search
that is given to a rational decision maker to narrow the search. Additional multicriterion
- search is conducted in the reduced search area. Again the decision maker narrows the
search. This process is repeated until satisfactory results are found. Of course, any reduc-

tion of search space also reduces diversity and eliminates a full exploration of the search

12

space for a multi-point search such as a GA. Pareto ranking plus niching tries to overcome
this downfall by maintaining diversity along the Pareto front to deter middling of the Pareto
front. Middling is a phenomenon that causes multi-point searches such as GAS (that are
searching for the Pareto front) to have a population converge to the middle or center of the
Pareto front. Even with controlled middling of the front, the search space is limited and not
fully explored.

Olhoff, Taylor [25] and Bendsoe et a1. [2] present three approaches to multicriterion
optimization problems using the calculus of variations: weighted sum, goal attainment
(bound method) and constraint method. All three are explained in the following para-
graphs.

The weighted sum method is stated as:

mig(;wigi) (1.4.1)
The design D is to be minimized, where gi is a Single optimization criterion that defines the
multicriterion optimization problem and w; is the weight assigned to each criterion. Each
weight is to represent the worth of each criteria based on the designer’s preference.

The constraint method is given as:

min(ck) . .
subject to Ci — bi S 0 (i=l,2,...,n i¢k) (1.42)
where the constraint bounds bi _>_ 0 are specified by the designer.

The goal or bound method is given as:

migfli) subject to Wigi — B S 0 (i=1,...,n) (1.4.3)
The object is to minimize D with prescribed gi over [3.
Osyczka, Kuchta and Czula [35] present an approach to multicriterion optimization

with expensive function evaluations. They maintain the most efficient approach to optimi-

l3

zation of multicriterion systems with computationally expensive evaluations can be ap-
proximated by replacing the expensive evaluation with an alternate evaluation function that
reflects the overall behavior of the system and is computationally cheaper while easier to
optimize. The search for an alternate cheap accurate evaluation function can be viewed as
another complex problem.

Belegundu et al. [1], Sundaresan et a1. [50], and Khot et a1. [36] present the idea of
robust optimization through desensitizing engineering systems to uncertainties, such as
manufacturing errors and operational variances. Traditionally, the sensitivity of optimal
designs are recorded after the product is manufactured. G. Taguchi [50] developed a meth-
od to measure sensitivity of a system during the inital design process. The method had
three stages: 1.) System design, where inital design varibles are chosen, 2.) Parameter de-
sign, where design parameters are optimized to desensitize the system to external influenc-
es, 3.) Tolerance design, where designers decide on allowable tolerances to represent an
optimal deign. The method was first developed at a time when little computer simulation
was available, so most data was gathered from actual experiments. The method was also
developed to be efficient with respect to the number of evaluations needed. Due to the lack
of computer Simulation, Taguchi used orthogonal arrays for approximating expected sen-
sitivity design values. The approximations made in the statistcal model are presented along
with examples of robust design of gears and beverage cans. The optimization of robust de-
sign of gears was achieved by desensitizing gear profile to manufacturing errors while
minimizing transmission error. The robust design of the beverage can was to minimize
weight by desensitizing beverage can geometry while lessening the sensitivity to manufac-

turing error.

14

1.4.2 Application of Genetic Algorithms to Engineering Systems

This next section describes applications of GAS to engineering systems. A wide ar-
ray of subjects is covered to demonstrate the diversity of GAS in engineering applications.
Engineering systems included are: general engineering, structural, composite and aerody-
namic systems.

To begin the general engineering applications, a paperby Crossely and Laananen
[8] on the conceptual design of helicopters to reduce weight for specified helicopter mis-
sions using GAS will be reviewed. The GA was combined with an industrial sizing code
specifically developed for helicopter design. Tournament selection was used with uniform
crossover and bitwise mutation. A G-bit improvement scheme was also applied, which is
analogous to a hill climbing approach, but was considered to be more like a gradient-like
bitwise improvement approach. This scheme is applied to an individual if no improvement
has been made in three generations. The G-bit improvement scheme varies one bit at a time
and places the most fit individual found back into the population, replacing the poorest in-
dividual. This approach requires more computational effort, but previous studies found it
to improve results. The GA run was terminated if no improvement was found for five gen-
erations. The design parameters the GA searched included integer, discrete and continuous
variables. Variables that were integer include the number of engines and rotor blades. Vari-
ables that were continuous were rotor tip Speed, disk loading, wing loading, wing aspect
ratio and maximum blade loading. Various turboshafts were represented with discrete vari-
ables. The observation was made that the binary coding can cause bias in the initial gener-
ation if a design variable is not a power of two. The authors claim this bias is easily

overcome if poor fitness is associated with the biased design variable.

15

Homaifar et a1. [23] used a simple GA in the optimization of an engineering system
consisting of turbofan engines. The criterion to be optimized were thrust per unit mass flow
rate and overall efficiency. The key parameters were Mach number, compressor pressure
ratio, fan pressure ratio and bypass ratio. Single criterion runs were first found for both ob-
jectives (thrust per unit mass flow rate and overall efficiency). The multicriterion optimi-
zation approach was based on the scalar aggregated weighted sum as noted previously by
Horn [24]. Simple approximations were introduced to predict the field flow equations.

The final two general engineering papers reviewed are applications of GAS in heat
transfer. Fabbri [10] used a GA to search for two-dimensional fin surfaces that have opti-
mal thermal characteristics. The thickness of a symmetric two-dimensional fin was mod-
eled with a polynomial and the distribution of temperature was obtained by solving
Laplace’s equation using the finite element method. The fin was considered to be im-
mersed in a fluid with a constant bulk temperature with an applied temperature at one end.
The fitness was based on compared effectiveness. Compared effectiveness is the ratio of
the heat flux removed by the fin of variable shape to a rectangular fin of same volume and
length. The polynomial was varied from first to fifth order in separate runs. The results
seem to be affected by the polynomial order. AS the order of the profile increased the os-
cillation of fin Shape also increased. The results seem to be an artifact of the modeling of
fin thickness as a polynomial. In fact, the number of local maximum and minimum thick-
nesses are identical to the order of the polynomial and none are similar with respect to fin
profile.

The last general engineering paper reviewed is an application of simple GAS in the

optimization of cooling of electronic components by Queipo et al. [41]. First, the combi-

l6

natorial problem of searching for the optimal or near optimal combination of various rect-
angular cross sections placed at equally spaced locations along a ventilated two-
dimensional channel was approached using a GA. The objective was to minimize the total
thermal failure rate of a set of in-line convectively cooled electronic components. A tran-
sient finite difference flow code was used in the thermofluids-optimized packing of elec-
tronic components. The failure rate is considered to vary according to the Arrhenius
equation, which is an exponential function. It was found that the heat transfer from elec-
tronic components is very sensitive to variations in geometry when convection is the dom-
inant heat transfer mechanism. Due to excessive computational effort (30 minutes per
evaluation) in solving the transient analysis, a population Size of seven that was run for sev-
en generations yielded at best a partial solution. A good first approach to any problem is
to consider a linear static approximation, which was not considered by Queipo et a1. [41].
A second problem presented extends the first problem with additional constraints of mini-
mizing total interconectivity wire length. The effects of buoyancy were neglected in the
balance of momentum equation. The last example reviews the effects of buoyancy-induced
flow on the previous examples. Buoyancy is essentially flow induced from the tendency of
heated fluids to rise. The effect of neglecting buoyancy caused slight errors in the Simpli-
fied two-dimensional study.

There are many papers concerning applications of GAS in the optimization of struc-
tures. The first reviewed is by Nakanishi and Nakagiri [32], which used a GA with a
boundary cycle and the finite element method to search for planar and three-dimensional
topological frames that have a minimal deformation at the applied loads. The weight was

constrained to be constant. The boundary cycle is used to represent the topology of a frame

17

in which there are no members that have a free tip. The authors assume that a frame struc-
ture that has members that are not connected to at least one other member cannot be an op-
timal design. This is because members with a free tip do not transmit any forces or
contribute any stiffness to the structure while wasting mass. Two and three—dimensional
example problems demonstrate that the boundary cycle increases efficiency by exploring
only reasonable portions of the search Space, reducing the size of search space significantly.

Jenkins [25] first gives motivations to use GAS in structural optimization. This in-
cludes efficient exploration of large search spaces and robustness with respect to global 0p-
timization. An introduction to binary coding of chromosome structure for discrete and
continuous variables is given. A simple GA is introduced with definitions of genetic oper-
ations such as crossover and mutation. A simple example of the minimization of a simple
algebraic function is presented. Penalty functions are introduced in order to satisfy con-
straints in the optimal design of a trussed-beam roof structure.

Haslinger and Jedelsky [22] concentrate on shape optimization using GAS. The
main point of this paper is that when using the finite element method, if a portion of the
structure is always constant, i.e. never changing, the stiffness matrix associated with the fi-
nite element method for this portion of the structure is also constant. Instead of continu-
ously recomputing the stiffness matrix, it should be stored and reused in order to decrease
computational expense. It should be noted that Haslinger and Jedelsky [23] do not mention
that the primary variables do not need to be computed when they are not of practical value
for certain locations in the structure. This also will reduce computational intensity associ-
ated with computing gradients.

Furuya and Haftka [l4] determined optimum placement of actuators on large scale

18

space structures using simple GAS and effective indices. Actuators are used to help control
the vibrational response of the structure. Effective indicies represent the vibration-suppres-
sion characteristics of the structure. An approximate solution to predict the effective indi-
cies was used in the analysis by ignoring the fractions of elastic energy contributed by the
actuators. This reduced the computational time required to compute the modal characteris-
tics by an order of magnitude. This approximation was justified by assuming the actuators
are small, not affecting the energy distribution of the structure.

Parmee and Vekeria [36] used island injection GAS to perform shape optimization
of concrete plates to minimize weight while maintaining an allowable stress level. The fi-
nite element method was used in the evaluation of the stresses in the variable thickness
plate. An island injection Genetic Algorithm (iiGA) approach was applied to help reduce
the computational intensity associated with a refined finite element mesh. Each island
evaluated individuals with a refined finite element mesh to accurately predict the response
of the plate. The iiGA searched at various levels of resolution by representing the plate by
various levels of geometric resolution. An iiGA has structured migration of individuals
amongst islands, which injects good solutions into islands with an increased level of reso-
lution. As an example, an island with a coarse geometric resolution represented a plate by
dividing it into four rectangular areas, allowing each rectangle to have a thickness that var-
ied linearly. A refined geometric resolution simply subdivided each rectangular, increasing
the number of thickness design variables. The process of injecting solutions from low lev-
els of resolution was continued until a satisfactory level of geometric resolution was ob-
tained. Islands with a low level of resolution converged much quicker compared to islands

with a higher level of resolution due to the reduced search space. This motivated the authors

19

to consider a dynamic island injection GA, where the islands at a low level of resolution
were terminated and reseeded after convergence. The iiGA helped decrease computational
time with little or no decrease in the fitness of the best individual found. The dynamic island
injection GA showed a further decrease in computational effort with little or no effect on
the final result.

The first application of GAS in composite structures reviewed is by Todoroki and
Watanabe [51]. They used a simple GA to optimize the placement and fiber orientation of
stiffeners on a laminated composite plate that had a bolt at the center. A uniform displace-
ment across one edge of a rectangular plate caused a stress concentration near the bolt.
Only a quarter of the rectangular plate was modeled due to symmetry. The objective of the
problem was to reduce the strain energy near the stress concentration by redistributing
stresses throughout the plate with composite stiffeners. The stiffeners were to be applied
after the curing of the laminated composite plate. Evaluation of fitness was performed with
a finite element code that used three-noded triangular plate elements. A location consisted
of a square composed of two triangular finite elements with a similar face. Three different
fiber orientations were used for each location: [0]2, [45/-45] and [90]2. Also, there was a
possibility that the location lacked a stiffener. Several simple GA runs were performed to
find good GA parameters with a similar simplified problem. For this particular problem,
the authors found that a low crossover rate (20%) with a high mutation (10%) performed
the best.

Haftka et al. [13, 14] used a Simple GA to optimize the stacking sequence of a lam-
inated composite plate for buckling load maximization. First, GA parameters such as pop-

ulation size, and genetic operator probabilities were tuned by numerical experiments.

20

Often GAS require a large number of analyses, sometimes in the millions. A serious prob-
lem with GAS occurs when a single analysis is computationally demanding. In light of this,
the objective of the work was to reduce the computational intensity associated with GAS.
One step made in this direction was the use of a binary tree. A binary tree was used to store
all previous analyses performed in the optimization and retrieved later if an exact duplicate
analysis was reiterated, avoiding costly repeat analyses. A binary tree can reduce compu-
tational costs when a high computational price is paid for each evaluation (compared to the
time used to search the binary tree). A drawback of the binary tree is the need for a large
amount of available memory to store previous design information. A second step made in
the direction of reducing the computational intensity of GAS was the approximation of cer-
tain buckling loads. Buckling loads were computed exactly for each design string created
by genetic operators. The set of all possible combinations of design Strings is approximat-
ed by a linear least square fit based on bending lamination parameters. The best of these
combinations replaces the nominal design in the population. This technique searches a
small neighborhood of designs for a local optimum. This particular local search caused
difficulties with Singular optima by interfering with crossover operations. This problem
was alleviated by seeding the initial population with certain designs.

Flynn and Sherman [12] outlined a two-step procedure to locate various GA con—
figurations for multicriterion optimization of flat aircraft panels. The panels were com-
posed of 13 materials (3 isotropic and 10 composite). A GA that searched for the Pareto
front was used in the multicriterion optimization based on four object functions: panel
buckling, bay buckling, weight and the total number of frames and stiffeners. Just as Haft-

ka et a1. [14] performed pre-trial runs to tune GA parameters, so did Flynn and Sherman

21

[12] in the first phase of their two-step procedure. ‘Rule of thumb’ values were gathered
from related GA publications. The GA parameters were then tested on a simplified aircraft
panel design. The second phase identified successful configurations and placed them into
a “preferred” set. Comparing results to other optimization techniques, the “preferred” set
of GA configurations produced designs that were lighter and at least as good in all other
criteria.

The next application of GAS in composite structures is by Punch et al. [27, 38].
They searched for optimal laminated composite beams with a coarse grained island injec-
tion Genetic Algorithm (iiGA). The objective was to find stacking sequences for laminated
composites beam that maximized the amount of mechanical energy absorbed by the beam
before fracture. A goal of approaching the problem with an iiGA was to reduce computa-
tional effort. This was accomplished with an iiGA by searching at various levels of spatial
resolution on separate computational nodes. Structured migration amongst islands allowed
good individuals at low spatial resolution to be injected into islands of higher spatial reso-
lution. The fitness of the beam was evaluated with a new finite element model developed
for composite laminate analysis. The new finite element model accounted for layerwise
variations of displacements and stresses by assuming a piecewise linear continuous
through-the-thickness in-plane displacement distribution. A clamped-clamped graphite-
epoxy composite layered beam had a point load applied at its midspan. The GA was al-
lowed to place either a 0 or 90 degree ply in each element in the discretized beam. The
concept of thin compliant layers to modify the load path characteristics of {the structure to
increase the amount of mechanical energy absorbed was modeled by allowing the GA to

choose either an identical material or a material with a significantly reduced stiffness above

22

each composite layer. First, a symmetric beam with a small search Space (approximately
106 designs) was explored with the island injection GA. The final results from the island
injection GA were confirmed to be the global optimal solution by enumerating (exhausting)
the search space. Next, larger search spaces that did not assume beam symmetry were ex-
plored with the island injection GA with promising results.

Mallot, et al. [28] used an island injection GA to optimize an idealized airfoil. The
object of the problem was to find composite stacking sequences of sandwich plates that
maintain an appropriate opposite twist to compensate for in-flight aerodynamic loads that
cause adverse airfoil twisting, while minimizing weight under stiffness and ply clustering
constraints. Displacements and stresses were predicted with a finite element model based
on first order shear deformation plate theory. A single “optimal” finite element mesh was
employed that balanced accuracy and computational efficiency. Three different runs ob-
tained from different GA topologies were presented. A single node, a “ring” t0pology GA
containing seven islands and an island injection Genetic Algorithm topology containing
seven islands. Each GA topology contained a total of 1400 individuals. In this study, both
GA Single node and ring topology converge within 50 generations. Using the island injec-
tion GA, the behavior sought from the islands exploring smaller search spaces (lower levels
of resolutions) should initially make more progress, providing “building blocks” to be in-
jected into islands with a larger search space (higher level of resolution). The islands with
a higher resolution should receive the migrating individuals, then outperform islands with
lower levels of resolution due to the expansion of search space. The iiGA approach was
able to find designs which reduce weight by 66%, increase minimum in-plane stiffness val-

ues by 33% and increased the twist by a magnitude of 77 times the baseline design, but in

23

opposite direction. The iiGA approach was the most efficient in terms of computational
time.

Quagliarella [40], Obayah [33] and Foster [1 1] apply GAS to design Shockless tran-
sonic airfoils. Quagliarella computed aerodynamic coefficients (drag and lift) of airfoils us-
ing a potential flowfield solver. A potential flowfield solver is a static linearized form of
the Navier flow equations. This study began by using a GA to for wave drag minimization
though Shape optimization of a baseline airfoil. The results obtained were characterized by
good behavior for the minimization of wave drag but presented poor results for lift coeffi-
cients. A new objective was formulated to simultaneously optimize drag and lift coeffi-
cients.

Obayah [33] approached the same problem for Shockless plus non-Shockless wave
constraints with a GA by an inverse design method. A direct optimization method couples
analysis methods to maximize (minimize) an object function by iterating directly on the ge-
ometry. The inverse design method in this paper deals with pressure distributions rather
than geometry, to minimize, for example, drag under given lift and pitching moment. Once
optimized target distributions are found, a corresponding geometry can be determined
through the inverse method. A direct approach requires computational fluid dynamic eval-
uations of each individual in a population. This is not feasible when a Single computational
fluid dynamic evaluation could require an incredible amount of computational time. This
makes inverse design methods attractive for any optimization approach in computational
fluid dynamics because it does not require any computational fluid dynamic evaluation.
This helps transform fluid flow optimization into a computationally tractable problem,

which is vitally important for GAS.

24

Foster [11] compares two hybrid optimization techniques, gradient and GA opti-'
mization methods in aerodynamic shape optimization in hypersonic flow. This flow region
was chosen because computational costs were not extremely high using an inexpensive
Newtonian flow analysis. The first flow analysis method was based on the modified New—
tonian impact theory, a simplified technique that is computationally cheaper, but can not
capture all the physical behavior associated with complex fluid flow phenomena. The sec-
ond flow analysis method is based on a parabolized Navier-Stokes solver, which can cap-
ture complex fluid flow phenomena at a slightly higher computational cost. The hybrid
gradient technique combined hill climbing, the method of feasible directions and Rosen’s
projection method. The method of feasible directions is particularly effective for determin-
ing a search direction when inequality constraints are present. If an equality constraint is
present then it is advisable to follow the constraint, which is performed with Rosen’s pro-
jection method. The hybrid GA applied the hybrid gradient technique if the entire set of
newly generated designs had been analyzed with no improvement. Examples included op-
timization of a half-sphere-cone to minimize wave drag and maximize lift, using both hy-
brid techniques. It was found that the hybrid GA convergence rates were comparable to
standard gradient techniques. Also, the hybrid GA had the ability to search design spaces

that gradient techniques could not.

1.3.3 Design Optimization Applications Via Other Methods
Chen et al. [5] used simulated annealing to place active and passively damped mem-
bers in the optimization of truss structures. Actuators were used to help control the vibra-

tional response of the structure. Unlike Furuya and Haftka [14], the elastic energy

25

contributed by the actuators was not ignored.

Chu et al. [7] have been researching structural optimization problems with stiffness
constraints. The paper presents a simple procedure with finite element analysis to mini-
mize weight while satisfying stiffness requirements that are placed on structures such as
bridges. At the end of each finite element analysis a sensitivity number indicating the
change in stiffness for each member of the structure is calculated. Members that have little
contribution to global Stiffness of the structure are slowly removed in a step by step process.
The sensitivity number for problems that have stiffness constraints is the inverse of the
strain energy of the structure, or the overall stiffness of the structure. Maximizing the over-
all stiffness of the Structure is equivalent to minimizing the strain energy. By removing one
of the structure members the global stiffness will be decreased and the change in the strain
energy can be defined (if the higher order terms are neglected). This change in strain energy
was the sensitivity number for Stiffness constrained problems. The process for calculating
the sensitivity number for problems with displacement constraints can be defined by first
applying a unit force at each of the nodes. Next, the displacements of the nodes under the
displacement constraints are calculated. A relationship between elements and the displace-
ment of the constrained nodes is defined as the sensitivity. The process is performed on a

short cantilever beam, Mitchell truss and a plate in bending.

1.3.4 Optimization of Flywheels
The optimization of flywheels has been attempted by many researchers [3,6,15-
l7,46,48]. All of these optimization techniques, expect for Bhavikatti [3], make a plane

stress assumption.

26

Bhavikatti [3] represented the shape of the annular flywheel as a polynomial and
used a single point sequential linear programming method with various objective functions
to optimize an annular isotropic flywheel using the finite element method with quadratic
elements. The first objective function was the minimization of the differences between the
maximum and minimum tangential stresses measured at 16 points along the radius of the
flywheel. This resulted in an awkward bulging near the inner radius of the flywheel. This
can be attributed to the fact that the highest stress will be in the tangential direction near the
inner radius of the flywheel. The second objective was to level the stress along the center
of the flywheel. This resulted in a flywheel that formed a shape that mainly concentrated
on increasing stress in areas of low stress, instead of reducing stress in areas of high stress.
The last objective function involved an aggregated weighted sum that tried to minimize
volume and level stresses. This resulted in a flywheel with less volume yet the same max—
imum stress as the results from the second objective function.

Often, optimization problems require simplifications to make the problem compu-
tationally feasible and numerically tractable. The results are often artifacts of these Simpli-
fications. Seireg [48], Sandgren [46], Christensen [6], and Georgian [17] all researched
shape optimization of flywheels in a day when the speed of computers was still relatively
slow. Assuming a plane stress approximation was the only possible approach to form the
problem that was computationally feasible.

Recently, Genta [16] used a simple GA to help design flywheels using a plane Stress
finite difference model in search of the constant stress profile. Genta acknowledged the
shortcoming of the plane stress analysis by using various constraints on the fitness function

to reduce the complex three-dimensional stresses that will occur, for example, at the inner

27

most portion of an annular flywheel that has Steep gradients in ring thickness. To avoid fly-
wheel ring oscillation that affected generations of flywheels, the initial population was

seeded with flywheels with thicknesses that only varied linearly in height.

4.4.5 Paper Review for Improvement

The iiGA has the powerful ability to search at various levels of resolution with dif-
ferent evaluation tools and fitness objectives. This approach by the iiGA can be applied to
any problem that the fitness of an individual can be first approximated with an efficient,
simplified analysis. For instance, Osyyczka et al. [35] attempt an approach on a lower level
by Simply approximating expensive objective evaluations with less expensive evaluations
that can hopefully reflect the overall response of a system accurately. A flaw in this ap-
proach lies in the fact that approximations in objective evaluations often give rise to solu-
tions that are artifacts of the approximation. An iiGA approach can use approximate fitness
evaluations to quickly search for low level building blocks to inject into islands with a re-
fined analysis. The island with the refined analysis can recombine good parts of existing
solutions (building blocks) to discover better solutions while discarding portions of the so-
lution that violate the approximate fitness evaluation.

Queipo [41] used GAS to optimize the cooling of electronic components, evaluating
fitness with a transient finite difference evaluation. Due to excessive evaluation time per
individual (30 minutes), a small population size was allowed to generate a few genera-
tions. This approach could be improved by making computationally efficient approxima-
tions in the analysis. An iiGA topology can be designed that could slowly decrease the

level of approximation until satisfactory evaluations occur. The approach could be im-

28

proved further using an iiGA to search at various levels of resolution with the simplified
analysis to discover building blocks that can be injected into a refined analysis.

Furuya and Haftka [8] approximate the modal response of large scale space struc-
tures by ignoring the elastic energy contributed from small actuators. A possible downfall
in the approximation occurs if many acutators are all placed in a small area, contributing a
large portion of elastic energy locally in the structure, changing the global response of the
structure. This optimization approach could be improved with an iiGA by injecting solu-
tions found with the approximations to an island that evaluates the objective function con-
sidering the elastic energy contributed from the actuators.

Crosseley and Laananen [8] used a G-bit improvement scheme to help improve the
search efficiency of the GAS in the conceptual design of helicopters. This improvement
scheme mimics Threshold Accepting by perturbing a solution to find a better solution on a
bitwise basis. However, while the G-bit improvement scheme performs a local enumera-
tive search on an individual, a Threshold Accepting algorithm climbs a hill in a dynamic
fashion. A search that has direction should be more efficient than an exhaustive local
search.

Many GA optimization approaches [1 1-13, 27, 33 43, 49, 51] can benefit through
the use of an iiGA to search at various levels of resolution to help improve computational
efficiency.

Typical examples of optimization approaches that result in solutions which are ar-
tifacts of an approximation of the objective evaluation are shown in [6, 7, 15-17, 46, 48].
The approximations of a plane stress objective function for flywheels in [6, 7, 15-17, 46,

48] resulted in the optimization tool exploiting the conditions of plane stress with flywheels

29

that have thick gradients in ring thickness. It is completely reasonable to first attempt to
approximate the response of a flywheel with a plane stress condition, but any optimization
technique that finds better solutions by violating the assumption will only boast solutions
that are artifacts of that assumption. Understanding approximations and their potential pit-

falls can be crucial when analyzing any structure.

1.4 Present Study

An overall objective of this thesis was to deve10p tools for optimization of large
scale three-dimensional composite structures. Combining a GA with the finite element
method is by now a familiar approach in the optimization of structures. Using an iiGA with
multiple evaluation tools with different fitness functions is a new approach aimed at de-
creasing computational time while increasing the robustness of a typical GA. New hybrid
GAS were also developed to help increase computational efficiency and robustness of a typ-
ical GA. Hybrid techniques developed include GAS and iiGAs that incorporate a combina-
tion of local search methods and TA algorithms.

For reader familiarization, the mechanics of GAS and hybrid GA techniques will be
reviewed. Mathematical development of two axi-symmetric finite element models will be
presented. Results from the developed finite element codes will be verified with an indus-
trial finite element code. An enumerable search space will be designed in order to compare
and contrast optimization approaches. Next, optimum solid and annular multi-material fly-
wheels will be sought. Finally, optimization of multi-material composite flywheels will be
presented while constraining angular velocity and penalizing flywheels with large “air gap”

growths.

Chapter 2

Algorithm Design and Comparison

2.1 Introduction

According to Carlson et al. [4], a good optimization scheme should not require full
exploration of the search space, yet should progressively find better solutions without ex-
cessive computation. Unfortunately, no Single approach is best for all classes of problems.
Fortunately, certain optimization approaches work better for a Specific class of problems.
Calculus-based hill climbing search methods work well for an explicit class of problems
that have a smooth convex objective function with a unimodal search space. Hill climbing
methods assume the solution space is somewhat well behaved, which is only true for a rel-
atively small class of problems. In fact, most search spaces are non-convex and multi-mod-
a1. Furthermore, a gradient-based approach requires that the problem have continuous
design variables, which is inadequate for discrete problems. Dynamic programming is a
hill climbing technique that can overcome some discrete problems but is limited to prob-
lems that have a relatively small search space. Enumerative search methods typically just
exhaust the search space by exploring all possibilities, which iS not efficient for problems
with a large search space and or a problem associated with high evaluation times. A GA is
a directed search procedure based on abstractions of the mechanics of natural selection and
natural genetics [12].

Traditional methods can be limited in multicriteria optimization applications while

GAS tend to excel with increased problem difficulty due to the increased search space com-

30

31

plexity. Typically, a GA can converge to a global optimum, or a nearly globally optimal
local optimum, by searching only a fraction of the search space. A possible downfall of
GAS is their computational intensity. A GA is an evolutionary, combinatorial technique
that is partial to problems with discrete solution spaces.

First an extensive background on the mechanics of Parallel GAS will be presented.
GAS that use ring topologies will also be reviewed. Next, the mechanics of Island Injection
GAS will be reviewed. Simulated Annealing and Threshold Accepting methods of optimi-

zation will also be presented.

2.2 Genetic Algorithms

Genetic Algorithms (GAS) are a powerful technique for search and optimization
problems, and are particularly useful in the optimization of composite Structures. The
search space for a composite structure is generally discontinuous and strongly multimodal,
with the possibility for many local sub-optimal solutions or even singular extrema. These
facts severely limit gradient-type approaches to optimization, bringing this broad class of
problems under scrutiny for application of GAS.

A GA is a search procedure modeled on the mechanics of natural selection rather
than simulated reasoning processes. Specific knowledge is embedded in a representation
called an organism, or in our case, a flywheel. A Single chromosome will represent an or-
ganism (flywheel). Each chromosome will consist of a set of values of a certain length.
Each chromosome of length n will be represented in the form of a vector, <x1,x2,
°°°°°°°°°°xn>. Each X, is an allele, or a gene. Often, an allele is expressed through binary

numbers, although many now use GAS with real numbers represented as floating point vari-

32

ables directly on the chromosome, together with O—mean Guassian random mutation oper-
ators. Each allele depicts a characteristic of a certain flywheel. Evolution takes place
specifically on chromosomes, not the “living creature” (flywheel) itself. Living creatures
are created and destroyed by a process of decoding chromosomes. Natural selection is the
connection between chromosomes and the performance of the chromosome structure and
is a process that causes superior performing chromosomes to have a higher probability of
reproduction, while inducing poorly performing chromosomes to have lower probability
for reproduction. During reproduction, three things are commonly modeled which can
cause children to differ from their parents: crossover, mutation and inversion. This in turn
produces flywheel “children” with decoded chromosomes that can be quite different from

9”

the “parents coded structures. The amount of crossover is controlled by the user. A high
crossover rate will produce new organisms quickly, but will have more probability to dis-
card higher performing organisms, as well. A low crossover rate can Stagnate the popula-
tion, thus producing no new organisms. A large genetic pool increases diversity in a
population, which in turn improves the ultimate results of the GA. A set of co-existing or-
ganisms (flywheels) defines a population, while successive populations are termed genera-
tions. Incorporating these processes in a computer algorithm will produce an algorithm that
solves problems in a manner reminiscent of natural evolution.

GA’s have search methods that are different than other search methods. GA’S have
the property of implicit parallelism. Implicit parallelism for a GA means that each evalua-
tion Searches (incompletely) in multiple hyperplanes of a given space. A GA does this ef-

ficiently by optimizing the trade-off between exploring new Space and exploiting the

information gained simultaneously. A hyperplane is a plane in hyperspace that corresponds

33

to a plane in ordinary space. A hyperspace is any space of dimension higher than euclidean
space (three dimensions). Large diverse populations increase the number of hyperspaces
in which the GA searches for answers in hyperplanes. This improves the GA performance
and computational intensity because the GA explore a larger search Space. GA’s explore
many solutions simultaneously, and information gathered is shared concurrently among
multiple searches. All of the searches have possible access to the information from other
searches simultaneously. This is almost like a pair of hands that work separately, while
knowing what the other hand is doing and having a common goal, only the GA can have
more than 2 hands.

Figure 2.1 displays the structure of a Simple GA. The simple GA begins by creating
a Single initial population, wherein flywheels of different shapes are randomly created. At
this point a finite element program evaluates the “goodness” of each flywheel based on Spe-
cific design parameters. Biased by the evaluations obtained, the GA uses unary and binary
operators on these designs to create another population. This population maintains the pre—
viously “good” flywheels while discarding poorly performing flywheels. The program
evaluates the new population members and continues with additional rounds of generation

and selection. This is repeated until satisfactory solution(s) are obtained [12].

34

 

Create the Initial Population

i

Evaluate Fitness |

t

Create The Next Generation by
Crossover, Mutation

onvergence

l Optimal Solutions I

 

 

 

 

 

 

 

 

 

 

 

 

Figure 2.1. Simple GA structure.

2.2.1 Genetic Algorithm Coding and Structure
2.2.1.1 Design Coding and Search Space Dimensions

This section will define the design codings used to represent flywheels in this study.
The design space and constraints must first be identified. The design Space and constraints
are identified by defining a maximum and minimum thickness, inner and outer radii, and a
specified number of material choices. The fidelity of each design space is prescribed by the
number of bits in the string defining each locus, or position on the chromosome. Each locus
will represent a design variable in the flywheel. A group of loci in vector form represents
a chromosome. For instance in a very simple case, the first two strings (composed of two
bits) in the chromosome might represent the inner and outer radii of the flywheel. The third
and fifth loci (three bits each) of the chromosome represent the inner and outer thickness of

the first ring. The third loci (one bit) iS the material property associated with the first ring.

35

This is demonstrated in the simple single-ringed flywheel chromosome shown below:
[01,10,111, 1,011]

The locus that represent continuous design variables such as ring thickness are
translated from binary to base ten quite easily. In binary mode, three bits convert into 23 or
8 equally spaced choices of heights ranging between the maximum and minimum thick-
nesses. For this chromosome there are 21 material property choices and 22 maximum and
minimum radii. Binary decoding schemes have a slight drawback due to the fact that the
number of choices (alleles) at loci that represent material properties must be a power of
two or a biased solution will be achieved (this was also recognized by Furura and Haftka
[13]). The total design space dimension is 22 - 22 - 23 ~ 21 - 23 = 2, 048 . This design space
could be quickly enumerated. The design space quickly grows: a six ringed flywheel with
constant inner and outer radii with a Single material property allowing three bits to define
the search Space for the ring thickness contains over 2 million search points. A GA can
search a fraction of this design space and find the optimal solution. Increase the fidelity of
the ring thickness representation to five bits and the search Space grows to over 343 billion
search points. This search space and much larger ones can also be explored efficiently with

a GA.

2.2.1.2 Initial Population Generation

First, the GA must perform initiation and evaluation of the initial population. Fol-
lowing the above coding design, the GA creates a population of flywheels of various Shapes
composed of various materials, stochastically. Often, populations are initialized by the us-

er, as in the case of Genta [16]. A measure of “goodness” or fitness is assigned to each in-

36

dividual by the fitness function.
2.2. 1.3 Reproduction/Selection

The next process the GA must perform is selection. The selection mechanism for
GAS is Simply the process that favors the selection of fit individuals to be placed in a mating
pool. Selection pressure is the degree to which the better individuals are favored: a higher
selection pressure corresponds to an increase in favoring of highly fit individuals. The fit-
ness function provides a means of comparing individuals and the selection process deter-
mines how the individual fitness values are compared. Selection mechanisms determine
which chromosomes will be placed in the mating pool. The mating pool is comprised of
the individuals that will reproduce to create new individuals through crossover and muta—
tion. Roulette wheel, stochastic remainder, rank-based and tournament are all selection
mechanisms. Tournament selection was used exclusively in this study due to its low selec-
tion pressure. Tournament selection allows for a group of individuals that are chosen sto-
chastically to gather and compete in a tournament, the individual with the highest fitness is
placed in a mating pool. Individuals from the mating pool are used to generate new off—
spring, with the resulting offspring forming the basis of the next generation. According to
Miller and Goldberg [23], a smaller group size lessens pressure to converge by disallowing
any “super” individual in initial populations to dominate the reproduction cycle, promoting
diversity. Reproduction is the GA search catalyst that works by mimicking natural selec-

tion.

2.2.1.4 Crossover

The next genetic process performed by the GA is crossover. Pairs chromosomes that

37

were placed in the mating pool through the selection process are combined through cross-
over. Uniform order-based, cycle, partially matched, one-point, two-point and uniform
crossover are commonly used procedures. One-point crossover was used in this study and
is explained in detail. Consider the two chromosomes shown below:

[01,10,111. 1.0111

[11,00,001. 0.0101
One-point crossover only occurs between loci. One-point crossover occurring between the
second and third loci would result in the following new individuals:

[01,10,001. 0.0101

[1 1 .0 0.1_1_1_._1_.fl_ll
In this particular case, the GA would create new individuals that have swapped their asso-

ciated radii.

2.2.1.5 Mutation
Mutation is another process that aids GAS in the search procedure. The mutation op-

erator iS applied at a low rate to chromosomes in order to help introduce new genetic data
into the population. This helps maintain diversity and reduces the probability of premature
convergence. A single-bit mutation operator was used in this study. A single-bit mutation
operator is demonstrated in the following chromosome at the third bit:

[11,10,111, 1,011]
The new chromosome would be:

[11,90,111, 1,011]

This single—point mutation would cause a change in the outer radius of the flywheel. When

38

real numbers are represents as floating point numbers on the chromosome, it is common to
use an alternative form of mutation, typically an additive Guassian Operator with 0 mean
and variance determined by the feasible range of the variable or the current diversity of the

values at that locus in the current population.

2.3.1.6 Penalty Method

GAS cannot directly satisfy all constraints in many optimization problems. A GA
can indirectly satisfy constraints through the use of the penalty method. Typically, the fit-
ness of an individual will be punished for violations of constraints. Often, the global opti-
mal solution lies near or at a constraint, so the penalty method can help the GA maintain

solutions that nearly or just slightly violate the constraint.

2.2.2 Parallel Genetic Algorithms

Two problems associated with GAS are their computational intensity and their pro-
pensin to converge prematurely. An approach that ameliorates both of these problems is
parallelization of GAS (PGAS), which also produces a more realistic model of nature than
a single large population. PGAS both decrease processing time and better explore the search
space.

Unlike some specialized sequential GAS which pay a high computational cost for
' maintaining subpopulations based on similarity comparisons (niching techniques, etc.),
PGAS maintain multiple, separate subpopulations which may be allowed to evolve nearly
independently. This allows each subpopulation to explore different parts of the search

space, each maintaining its own high-fitness individuals and each controlling how mixing

39

occurs with other subpopulations, if at all.

2.2.3 Island Injection Genetic Algorithms

Island injection Genetic Algorithms (iiGAs) represent an approach to search at var-
ious levels of resolution within a given space. This includes first searching at low levels of
resolution on different nodes (islands) and then injecting the high-performance individuals
into an island of higher resolution to “fine-tune” them. Islands with a low resolution can
evaluate fitness with a simplified analysis that is computationally cheaper, while islands
with a high resolution can use a refined, computationally expensive analysis. Also, differ-
ent GA parameters can be used for each population. The rate of crossover, mutation, and
the island interaction can all vary from island to island. For example, islands at low resolu-
tion can have a larger population size to take advantage of a computationally cheap fitness
evaluation. Also, islands with a computationally cheap fitness evaluation and a low reso-
lution can be further exploited by increasing the number of generations evaluated before
injecting the information to other islands. Figure 2.2 shows a typical iiGA that can first
search at a low level of geometric resolution and inject individuals into islands with a high-
er level of geometric resolution. iiGAs can use multiple evaluation tools, various fitness
functions and performs multiple refinements in the geometric representation by increasing
the number of rings in the flywheel. For composite analysis, material properties can vary

from ring to ring.

40

Three Rings
(low resolution)

Six Rings
(medium resolution)

Twelve Rings
(high resolution)

Twenty—Four Rings
(higher resolution)

 

Figure 2.2. Possible iiGA topology that searches
at various levels of geometric resolution.

Multi-objective optimization requires combining many performance measures to
find an optimal design. Each individual fitness measure may have its own optimal or sub-
optimal solutions. By using iiGAs, each individual performance could be used as a “sub-
fitness” function since it represents “good” designs within the search space. Often, only
the best individual in a population is allowed to migrate, to allow “good” ideas to be com-
bined with other “good” ideas to find “better” ideas amongst islands of different “sub-fit-
ness” functions. Next, the individuals are injected into a final node where the evaluation
of an overall fitness function is employed. This search method ensures a robust explora-
tion of the search space for all aspects of the overall fitness. Of course, many variations on
these island injection architectures can be custom tailored for specific problems.

iiGAs have the following advantages over other PGAS. Building blocks of lower
resolution can be directly found by search at that resolution. After receiving lower resolu-
tion solutions from its parent node(S), a node of higher resolution can “fine-tune” these so—

lutions. The search space in nodes with lower resolution is proportionally smaller. This

41

typically results in finding “fit” solutions more quickly, which are injected into higher res-
olution nodes for refinement. Nodes connected in the hierarchy (nodes with a parent-child
relationship) share portions of the same search space, since the search space of parent is
contained in the search space of child. Fast search at low resolution by the parent can po-
tentially help the child find fitter individuals. iiGAs embody a divide-and-conquer and par-
titioning strategy which has been successfully applied to many problems. Homogeneous
PGAS cannot guarantee such a division Since crossover and mutation may produce individ-
uals that belong to many subspaces, i.e., the divisions cannot be maintained. In iiGAs, the
search space is fundamentally divided into hierarchical levels with well defined overlap

(the search space of the parent is contained in the search Space of the child).

2.3 Simulated Annealing and Threshold Accepting

Simulated Annealing (SA) is a combinatorial optimization technique that is based
on the statistical mechanics of annealing of solids [45]. To understand how such an ap-
proach can be used as an optimization tool, one must consider how to coerce a solid into a
low energy state. Annealing is a process typically applied to solid materials to force the
atomic Structure of the material into a highly ordered state. Atomic structures that maintain
a highly ordered state are also at a low energy state. In an annealing process, a material is
heated to a temperature that allows many atomic arrangements, then cooled slowly, mini-
mizing energy, while statistically allowing an occasional increase in atomic energy. When
the material is extremely hot, the probability of an increase in atomic energy is very high.
As the cooling continues, the probability of an increase in atomic energy decreases. Simi-

larly, SA methods use analogous set of parameters that simulate controlled cooling effects

42

found in the annealing of materials.

SA methods begin with an initial solution that is often generated randomly, and try
to perturb the solution to improve it. If the perturbation improves the solution then it is ac-
cepted and the process of perturbing continues. In this manner, SA methods are like itera-
tive methods that climb hills. As with hill climbing methods, this process of searching just
for a better solution tends to force the process to a local optimum. However, SA methods
are different in this manner: annealing occasionally allows perturbations that are harmful
to the solution to be accepted. This allows SA methods to climb out of local optima to
search for the global optima. In actual physical systems, jumps to a higher (harmful) State
of energy actually do occur. These jumps are monitored by the current temperature. As
the annealing process (cooling) continues, the probability that only better solutions will be
accepted increases. At the beginning of the annealing process (associated with a high tem-
perature), the chance that a harmful solution is accepted is higher while later in the anneal-
ing process (at a lower temperature) that chance that a harmful solution is accepted is small.

This probability of accepting harmful solutions is based on a Boltzmann distribution:
-AE
Pr[accept] = e T 2.1)

By successively lowering the temperature T, the simulation of material coming into equi-
librium at each newly reduced temperature can effectively simulate physical annealing.
Threshold Accepting (T A) is a simplified version of Simulated Annealing. The
probability of accepting a harmful solution is governed by the Boltzmann distribution for
SA applications and TA algorithm, but the TA algorithm is not dependent upon a specified
temperature. Instead the TA algorithm “cooling” is based on a specified percentage of the

current solution’s fitness. This percentage decreases over the set of generations. This

43

causes the TA in lower generations to have a higher probability of accepting a worse indi-

vidual, while later generations in the optimization are less likely to accept a worse solution.

Chapter 3

Finite Element Models of Flywheels

3.1 Introduction

Two axi-symmetric finite element models were developed to predict planar and .
three-dimensional stresses that occur in flywheels composed of orthotropic materials un-
dergoing a constant angular velocity. Both finite element models can be developed by ap-
plying the principal of minimum potential energy. An alternate but equivalent approach to
formulating the finite element model is based on the weak statement of the governing dif-
ferential equation.

To insure realistic results in an optimization problem the evaluation too] should be
verified for accuracy. The axi-symmetric plane stress finite element model is verified by
comparing results with the closed form solution [46] for constant thickness solid flywheel.
Results from the axi—symmetric three-dimensional finite element model are compared to

the commercial finite element code MARC [47] for a flywheel of variable thickness.

3.2 Finite Element Model Formulations

In this section, two axi-symmetric finite element models that predict the response
of an orthotropic flywheel undergoing a constant angular velocity will be developed. First,
the plane stress finite element model will be developed. Next, the plane stress approxima-

tion will be ignored to formulate the three-dimensional finite element model.

44

45

3.2.1 Plane Stress Finite Element Model Formulation

The finite element model that assumes a planar state of stress assumes no variation
of stress will occur through the transverse normal (thickness) direction of the flywheel. Al-
so, the model assumes that the transverse normal and in plane shearing stresses are small
when compared to the tangential and radial stresses. The plane stress finite element model
is accurate when the radial gradient of flywheel ring thickness is small.

The equilibrium equation in polar coordinates for a flywheel in a plane stress state

dtrO'r
dr

where the thickness t is a function of the radius. Figure 3.1 displays how the geometry of

 

—0’et+trBr = O (3.1)

the flywheel was modeled with concentric rings that varied in thickness linearly as a func-
tion of the radius. The body forces induced from rotations in the radial direction, Br can be
defined as:

B, = rp m2 (3.2)
where (D is the angular velocity and p is density.

Next, strains in cylindrical coordinates for plane stress axi-symmetric problems are

defined as:
e, = u,r (3.3)
u
89 = ; (34)

Hooke’s law for orthotropic material in the condition of plane stress:

[or] = [C11 C12 5r (3.5)
99 C12 C22 59

46

The governing differential equations for plane stress axi-symmetric problems with

orthotropic materials is defined by substituting equations 3.3-3.5 into equation 3.1 yielding:

ea[IT(CH%—u+C12%)}II(C12%§+szuj-“i-trB= 0 (3.6)

The plane stress weak statement:

{/LG (t r(C11%—+C12L—r‘ D—t 1(C12g—I’3+C22L;l -)+trB )rdrdedz = 0 (3.7)

Introducing finite element approximations:

L = ‘1". (3.8)
= 2 141-1in (3.9)
j=1

and integrating by parts the plane stress stiffness matrix becomes:

aw, 8‘? ‘1' aw]
Ki}. = K5. [tr(CnarJ +C12— r”)]+~il.t([cna— +sz “-3) drdOdz (3.10)
V

The plane stress forcing vector is:

{f }= in? r}rdrd9dz (3.11)

The plane stress boundary conditions are:

EBC’S NBC’s
u tr(6,nr) (3.12)

A three-noded element with one degree of freedom per node using Lagrangian qua-
dratic shape functions in the plane stress finite element model is Shown in Figure 3.2a. A
four point Guassian quadrature rule was used to numerically integrate the Stiffness matrix.

A two point Guassian quadrature rule was used to evaluate stresses.

47

3.2.2 Three-dimensional Finite Element Model Formulation

The finite element model that yields a three-dimensional stress state is truly a two-
dimensional finite element model and is accurate for all shapes. It captures the variation of
stress through the transverse normal direction, and makes no assumptions about the trans-
verse normal and in-plane shearing stresses.

Assuming that the displacement field does not change as a function of angle, the prob-
lem can be modeled with a planar mesh. Figure 3.3 displays a planar symmetric cross sec-
tion of a flywheel to be analyzed. Only half the cross Section of a flywheel needs to be
modeled due to symmetry. Symmetry boundary conditions (no displacement in the trans-
verse normal direction) were applied along the horizontal symmetry plane of the flywheel.
Forces induced from angular rotation were included as body forces. The essential boundary
condition at the center of the flywheel was not directly applied due to the coupling of the
stiffness matrix, enforcing displacements according to the physics of the problem.

The formulation of the three-dimensional axi-symmetric finite element model be-

gins by defining the equilibrium equations in cylindrical coordinates:
o,‘,+rzr,r+%(o,—oe)+Br = O (3.13)

1
17,2, , + oz, z + fr: = O (3.14)

Next, strains in cylindrical coordinates for three-dimensional axi-symmetric problems are

defined as:
E, = u,r (3.15)
u
89 — ; (3.16)
8 = w, (3.17)

48

Yr: = u’z +W9r (3.18)

Hooke’s law for three-dimensional axi-symmetric orthotropic material:

0, C11 C12 C13 0 3,
oe = C12 C22 C23 0 so (3.19)
Oz C13 C23 C33 0 8;
-rz. LO 0 0 Caitlin:

 

 

 

 

 

 

The governing differential equations for three-dimensional axi-symmetric problems with
orthotropic materials are defined by substituting equations 3.15-3.19 into equations 3.13

and 3.14 yielding:

l a a a a
(C113u+c‘2gcl3aw C123“ C22; C23T:)+
(3.20)
a 3“ a 3.. aw
ar(C“8r+Clzl;l )+ +azC55(d_z+ 87)”; = O
a a a a 8 a
a—r(rC55(d_: + 3%))4' az(C133_:l + C23? + C33aw)+
(3.21)
a flu 8w 8 Bu 3w
ar(rC55(d—z+$))+ 72(C13ar4' C23“ 7+C33a—z )th = 0
The weak statement can be written as:
a a a
JL(1(CHau+C‘ 2+Cl3a—W_C128—u—C22r—C23a_:)+
(3.22)

a 3 a a a
’a—r(Cllau+C12:l)2+‘a'-z C55(a:+ 8—2))!- B,)rdrd0dz = O

49

d Bu 8w 3 Bu SW
{324 ail’aC55(z 37D a‘Elcwar” WC 2+C333z)

a a a a Bu aw
a—(rcss 8-:+a_?))+a—2(C133_r +C23u+C33aZ )+
1(C55@_“+ g—W))+Bz)rdrd6dz = 0
Z r (3.23)
Introducing finite element approximations:
u = 2 ujwj (3.24)
j=l
w = 2 W]. (3.25)
i=1

and integrating by parts, the three-dimensional axi-symmetric stiffness matrices become:

B‘P- a‘I’- W-
11 _ i j _J
Kij — {(37 [Clly +C12 r]+

r, av}. ‘1']. _a_r 8—?

(3.26)

3‘? as! ‘1’,- awj aw aw].
2:“? icl: 1332 ]+ 7[C23'& ]+ a—Z [(3555: Drdrdfldz (3.27)
V
3‘? B‘P- ‘1']. ‘11,. aw].
”=[bz [C13a—Z’+C23—]+ 7[C55'a—z Drdrdedz (3.28)
V
aw air! ‘11 am
22:...(32 iCl: 3382]]... 'r_i[C558r JDrdrdOdz (3.29)
V

The axi-symmetric three-dimensional forcing vectors are:

50

{f1} = j {Brr}rdrd0dz (3.30)
V

2
{f } = 0 (3.31)
Where Br is the body forces in the radial direction induced from rotations.

The axi-symmetric three-dimensional essential and natural boundary conditions

are:
EBC’S NBC’S
r(0',n + T, n )
u ' z Z (3.32)
w r(O'Znz + trznr)

Equations 3.26-3.31 can be presented to form the complete axi-symmetric three-dimen-
sional finite element model in following matrix form:

11 12 l
K K H = f (3.33)

K21 K22 w 2

The three-dimensional axi-symmetric finite element model was implemented into a
finite element code that used four noded elements with two degrees of freedom per node
(Figure 3.2b) with linear Lagrangian interpolation functions. A two point Guassian quadra-
ture rule was used to numerically integrate the stiffness matrix. A one point Guassian
quadrature rule was used to evaluate stresses. An automated mesh generator was written to
allow for element refinement through the transverse normal and the radial directions.
Therefore, the finite element code that predicts three-dimensional stresses can have various
levels of refinement. A coarse mesh with a small number of degrees of freedom will be less

accurate but more computationally efficient compared to a refined mesh containing a larger

number of degrees of freedom. The mesh was also generated in such a way as to minimize

51

the time required to solve the set of banded linear equations created by the finite element

code.

3.3 Verification of Results

Figure 3.4 displays results comparing the analytical solution to the plane stress fi-
nite element computational solution for the tangential and radial stresses in a flat solid fly-
wheel as a function of flywheel radius. Ten quadratic elements are used in the plane stress
finite element model. Stress contour plot obtained using the MARC analysis code displays
an arrow which represents the cross-section where the stresses are compared. There is little
difference between the plane stress FEM solution and analytical solution for the radial and
tangential stresses. Figure 3.5 displays results comparing the analytical solution to the
plane stress FEM computational solution for radial and tangential stresses for a flat annular
flywheel composed of 20 quadratic elements as a function of flywheel radius. Again, the
plane Stress FEM solution agrees with the analytical solution.

To compare three-dimensional results with a solution obtained using MARC, the
variation of stresses were plotted through the thickness near the inner radius, being very
close to zero. The flywheel shape has moderately Steep ring gradients, and was chosen to
demonstrate the capabilities of capturing the variation of transverse normal and in-plane
shearing stresses with the axi-symmetric three-dimensional finite element code. Figures
3.6-3.9 compare the solutions from MARC and the present FEM code as a function of fly-
wheel thickness. Again the MARC analysis stress contour plot displays an arrow which
represents the cross-section where the stresses are compared. All components of the stress

agree with the solution obtained by MARC. It should be noted that the solution obtained

52

by MARC measures the stresses at the center of each element and then interpolates to the
nodes of each element. The present FEM code measures stresses (with no interpolation) at
the center of each element. Interpolating the stresses from the center of the element to the
nodes will cause slight variations in solutions obtained from MARC and the present FEM
code. This is more evident near the surface of the flywheel, where a stress free boundary
condition exists.

The three-dimensional analysis convergence was sought by refining the mesh from
30, to 120, 480 and finally 960 linear elements. The difference in maximum Stress in the
problem converged from 15.1% to 10.1% to 2.1% respectively. The transverse normal
stress converged much more Slowly. This could come from the domination of the tangential
and radial stresses in the minimization of potential energy. Often, the minimization of po-
tential energy sacrifices the accuracy of a non-dominating stress component in order to ac-

curately capture the dominating stress components.

53

Z
R ;-_—7- ‘p v j R
e

   

 

 

 

 

   
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Top View Side View
Figure 3.1. Flywheel geometry.
11
w
u
u
o—o—o—>
3.2a. 1 Dimensional Element 32b 2 Dimensional Element
(1 DOF per node) (2 DOF per node)
Figure 3.2 Element topology.
"‘— ‘ r Transverse Normal
L J T“ - - Tan ta]
«3 D t gen 1
Av ijflw ‘1 "cc ion Direction
_ 12.1..1~__,______ -
3 H“ “ Radial
. M“ Direction
Axrs of “M
Rotation
\iAxis of Symmetry LSymmetry Boundary Conditions

Figure 3.3. Typical axi-symmetric planar finite element mesh with symmetry
boundary conditions (1,134 degrees of freedom).

54

 

 

 

 

 

 

 

 

 

0006000 VP
\\ —.— FEM 3.61.! Ctr...

7.005000 0
——.— An.lytlo.l H.dl.l 0"...
—.— FEM T.ng.ntl.l Ctr...
—.— An.lytlc.l T.ng.ntl.l 8"...

0.005000 «

lace.“ 4L

3 C 005.00 4»

3.00E000 «v

2.00:.00 4

1.005006 0

0006.00 . . 4 ‘

0005.00 1005.00 0005.00 0002.00 0002.00 0002.00 0005.00 7005.00 0005.00 0.00:.00 v.005.01

Hutu.

Figure 3.4. Tangential and Radial Stresses as a Function of Radius Predicted
by Plane Stress Finite Element Code and Analytical Solution for a Flat Solid
Flywheel. Arrow Represents Cross Section Analyzed.

‘ ..°EO°7 ‘P

 

 

 

 

 

1.005907 0
—.— FEM ll.¢llnl Btr... ‘
-—.—-— An.lyflo.l R.dl.l
. .1 an...
‘ ”5.07 —.-- FEM 1’.li
an...
—.— AIS-Iva“! T.n¢.ntl.l
an...
1.20E007 1
" LOOEOO‘I 4»
3 0006.00 4
0.005000 0-
4,005.00 .1
2005-00 4 A
0.005000 4

 

000E000 LOOEOOO 2.00E000 8.006.900 ‘.006000 0.005000 0.0050” 7.0059!” 0.00I000 0.00E000 1.00600‘
mm.

Figure 3.5. Tangential and Radial Stresses as a Function of Radius Predicted
by Plane Stress Finite Element Code and Analytical Solution for a Flat Annular
Flywheel. Arrow Represents Cross Section Analyzed.

55

24:05.07 .

v “007
—.— MARC
+ Pro-0n! FEM Coda

Radial 51ml (PI)

 

 

 

0.!!!“ 1.0060! {WE—CR sane-o: 4.8““ 5.0064. tax-m [“02 0006-0:
Thlckn...

Figure 3.6. Variation of Radial Stress Through the Thickness for the MARC
and Present Finite Element Code for a Solid Isotropic Flywheel that Varies
in Thickness. Arrow Represents Cross Section Plotted.

2 (new? .—

1.05.07 .. ‘,
‘4

+ MARC
4~ Present FEM
Code

1.20E007 4i

Tangentail Sims (Pa)

i
t

 

 

0.00E000
0.0015000

 

LOGS“ 1.005“ 3.0(1-02 ‘ DOE-02 5.0064” OVOOEM 7.00602 0.00502
Thlckn...

Figure 3.7. Variation of Tangential Stress Through the Thickness for the
MARC and Present Finite Element Code for a Solid Isotropic Flywheel
that Varies in Thickness. Arrow Represents Cross Section Plotted.

Radial Strut (Pa)

56

   

—.— MARC
~- ~ Present FEM Cod.

tons.“ »

 

monsoon
0.0015000 Loos-oz zone-la aooi—m 0.00:4): 5.00500 anoint 7.00:0: tonic:

Figure 3.8. Variation of Transverse Stress Through the Thickness for the
MARC and Present Finite Element Code for a Solid Isotropic Flywheel
that Varies in Thickness. Arrow Represents Cross Section Plotted.

 

 

   

—.— MARC
3-- Pro-om FEM Cod.

Shearing Sims (Pa)

 

OWE.” INS-4H acacia: Sll‘m 4mm Iu-DZ 000:4: 7““ 0““
Thick"...

 

Figure 3.9. Variation of Shear Stress Through the Thickness for the
MARC and Present Finite Element Code for a Solid Isotropic Flywheel
that Varies in Thickness. Arrow Represents Cross Section Plotted.

Chapter 4
Optimal Design of Flywheels Using Genetic Algorithms

4.1 Introduction

Various GA approaches are presented in this thesis to search for flywheels with
varying shape to maximize SED, while controlling angular velocity and “air” gap growth.
A simple Genetic Algorithm (SGA) will be compared to a set of optimization techniques.
This set includes parallel GAS with topological ring structures, Island Injection Genetic Al-
gorithms (iiGAs), local improvement methods, Threshold Accepting (TA) algorithms and
hybrid GAS. First, the optimization of solid isotropic flywheels will be presented. Next,
comparisons of the set of optimization techniques for an enumerable search space will be
presented. Optimization of annular and solid composite flywheels are presented. Finally,
“safer” annular flywheels (while constraining angluar velocity) will be sought by penaliz-

ing flywheels that have a large “air gap” growth.

4.2 Description of Optimization Problem

Structural optimization problems have been around for centuries and many have
been reviewed in this thesis. The formulation of this optimization problem begins with the
definition of overall fitness, in this case the maximization of Specific Energy Density
(SED). Conflicting interests arise when performing optimization of flywheels. The criteria
considered are: maximize SED, minimize “air gap” growth while controling high angular

velocities. The minimization of “air gap” growth is considered to find “safer” flywheel de-

57

58

signs. SED can be defined as:

 

(4.1)

The definition of SED contains angular velocity squared, so as we search for flywheel
shapes with a controlled angular velocity we are also limiting SED. The GA will then be
able to optimize the fitness function by placing more mass near the end of the flywheel, in-
creasing the mass moment of inertia. “Air gap” growth (which is the maximum displace-
ment of the inner radius of an annular flywheel due to forces induced from angular
rotations) will typically be higher for flywheels that have a higher SED value due to the
association of higher failure angular velocities.

Constraints, such as maximum angular velocity and air gap growth, are here en-
forced in the GA by penalizing individuals that violate specific constraints. In this case, a
constraint on maximum angular velocity was not practically known. The goal of decreas-
ing angular velocity is to search for safer flywheels. Therefore, an initial GA run was per-
formed without constraining angular velocity. Next, the values for the constraints on
maximum angular velocity were formed based on the maximum angular velocity for the
best individual found in the initial GA run. This approach was also taken with constraints
on the “air gap” growth in flywheels.

The fitness was defined as:

C] SED _lefailure_PzGAPfailure
SEDmax GAPmax

subject to

Fitness(maximize) = 4.2)

(Omax

or, min 5 0'r 3 Cr» max, 0'0, min 5 00 5 0'0, max

CZ, min S OZ S OZ, max , Trz, min STI’Z S TrZ’ max

59,

where C] is the value used to weight the normalized fitness and Pi are penalty coefficients
used if an individual violates a constraint. This definition of fitness is based on the maxi-
mum stress failure criteria for isotropic materials which can easily be transformed into
terms of strain for composite analysis. If the angular velocity was lower than a prescribed
amount then the coefficient P1 was set to zero. Otherwise, the values for the coefficients
in equation 4.4.1 are given in Table 4.2.1.

Table 4.2.1. Weighting Coefficient Values.

 

C1 P1 P2

 

100.0 40.0 20.0

 

 

 

 

 

As mentioned, an initial GA run was performed to detemiine the maximum SED,
angular velocity and “air gap” values without penalizing for constraints and normalizing
equation 4.2 with unity. The constraints were determined by taking a percentage of the

maximum angular velocity while penalizing based on “air gap” growth.

Top View Side View

I I ’

   

Figure 4.1. Visual Display of Flywheel.

60

4.3 Model Description

Figure 4.1 is a visual display of an annular flywheel that shows how the flywheel is
modeled as a series of concentric rings. The thickness of each ring varies linearly in the ra-
dial direction with the possibility for a diverse set of material choices for each ring. When
using the axi-symmetric plane stress analysis, a three-noded “line” element was used to
represent the flywheel (see Figure 4.2). The axi-symmetric three-dimensional finite ele-
ment model required a planar mesh. Figure 4.3 displays a typical planar mesh generated,
using symmetry to increase computational efficiency for the axi-symmetric three-dimen-
sional analysis.

The flywheel dimensions were chosen based on research provided by [44]. The
maximum and minimum thickness, inner and outer radii were taken from [44] to ensure re-
alistic dimensions.

An aluminum alloy material was chosen for isotropic flywheels due to its high
strength low, density characteristics. Table 4.3.1 contains all isotropic material properties
used in this thesis. For composite flywheels, quasi-isotropic material properties were cho-

sen and are given in Table 4.3.2.

61

 

Three noded plane stress “line” element.

Figure 4.2. Plane stress finite element representation.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

“NT Transverse Normal
, ~~~~E1TN\\ Direction .
4‘ .__ "ix“ hm Tangential
THNHN NNLN \\ . .
4.. ”NFL N~D~N Direction
NNHNTM~H NENN~H~WN\DN
”-11-...“ Ni~~~ “‘NNN ~~\\ Rad.al
~~M-Hhi H~~~ Hfi~:\\\\ l
i~~~NHu~- NFN‘H ~~~~ WQHHQQK D. .
Ansof ~~~4~.-.-:~+~L~ :::::::: 2:355:13»
.' .‘..;.'.A'. 1‘1"“? 5""),11 ":‘r‘.")“1,f :3}.;.‘I‘,'-‘.“_‘L.I‘,"."11,",‘L‘;~i{~ gm‘ “11"“[‘2"“"::‘1r'1u KENT-j=—
1
Axis of Symmetry \—Symmetry Boundary Conditions

Figure 4.3. Typical planar finite element mesh with symmetry boundary
conditions used in the three-dimensional fem evaluations (1,134 degrees
of freedom).

62

Table 4.3.1. Isotropic Material Properties.

 

 

 

 

 

- Tensile Compressive
De t
Material E (GPa) v “813), Strength Strength
(Kg/m ) (MPa) (MPa)
Steel 200.0 0.30 8.86E3 400 400
Aluminum 720- 0.285 2.80E3 480 480

 

 

 

 

 

 

Table 4.3.2. Composite Material Properties.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Property E-Glass0 E-Glass 0
Epoxy [0 ] Epoxy [90 ]
E (GPa) 18.96 18.96
Density (Kg/m3) 1.6E3 1.6E3
v 0.27 0.27
G(GPa) 7.47 7.47
e,(tensile) 8.99E-3 15.9E-3
3: (tensile) 15.9E-3 8.99E-3
ez (tensile) 8.99E-3 8.99E-3
an 31.37E-3 31.37E—3
e, (compressive) 22.47E-3 8.21E-3
e, (compressive) 8.21E-3 22.47E—3
t:Z (compressive) 22.47E-3 22.47E-3

 

 

 

 

 

63

The maximum stress failure criterion was used to predict the maximum failure an-
gular velocity in the analysis of isotropic flywheels, while the maximum strain criterion was
used for composite flywheels. The absolute value of the maximum stress (or strain) can be
used to scale the angular velocity and “air gap” growth. By assuming an initial angular ve-
locity, the stresses and strains were calculated by the finite element method. A failure index
was defined as the largest ratio of the absolute maximum stress (or strain) to the maximum
allowable stress (or strain) produced by the initialized angular velocity. The failure index
was used to scale the initialized angular velocity to the maximum failure angular velocity.
The failure index was also used to scale the “air gap” growth to the “air gap” growth at fail-

ure. Also, the origin of the maximum stress (or strain) predicted the flywheel failure ring.

4.4 Island Injection Genetic Algorithm Approach

This next section will describe the advantages of first exploring the search space ef-
ficiently at various levels of resolution using an efficient and simple calculation performed
by a finite element code that assumes that the flywheel is in a plane stress state, basing fit-
ness on a sub-fitness function using an Island Injection Genetic Algorithm (iiGA). The
iiGA can discretize an objective function by breaking it down into separate sub-fitness
functions used to evaluate fitness on various islands. Building blocks are first sought using
a low level of geometric resolution with a plane stress assumption basing the fitness on a
sub-fitness function. Good solutions found at a low level of geometric resolution based on
a plane stress assumption with the sub-fitness function are then injected into islands that
measure fitness with a refined three-dimensional finite element code at a higher level of

geometric resolution, basing fitness on an overall fitness definition.

64

4.4.1 Island Injection GA Geometric Resolution

Figure 4.4 displays an iiGA topology that first searches at a low level of geometric
resolution (three rings) for building blocks needed to inject into islands of medium geo-
metric resolution. At the end of a prescribed number of generations per cycle a file is writ-
ten containing good solutions from the existing population. The iiGA then searches at a
medium level of geometric resolution (six rings) to find building blocks to inject into a high
level of geometric resolution (twelve rings). One complete cycle ends by generating and
evaluating the populations on the islands with the highest level of geometric resolution. The
next cycle begins with the islands with the lowest geometric resolution (three rings) search-
ing to improve the existing building blocks in the island. Now, the islands with a medium
level of geometric resolution (six rings) receives migrating solutions from the file contain-
ing good building blocks from the island with the lowest level of geometric resolution. The
process of migrating building blocks from a lower level of geometric resolution to a higher
level of geometric resolution continues. Islands with the highest level of geometric resolu-
tion receive building blocks from the islands with a medium level of geometric resolution.
This process of searching at low levels of resolution to improve building blocks is contin-

ued and is exploited to accelerate the convergence of a typical GA.

65

Three Rings
(low resolution)

Six Rings
(medium resolution)

Twelve Rings
(high resolution)

Figure 4.4. iiGA t0pology that performs multiple refinements in geometric resolution.

 

 

- . . llllllll

4.4.2 Efficient Evaluation Tools

The finite element codes can evaluate fitness at various levels of efficiency. The
axi-symmetric plane stress finite element evaluation is extremely efficient, but less accu-
rate for odd shapes. The axi-symmetric three-dimensional finite element code will also
have various levels of efficiency. A three-dimensional finite element model with a large
number of degrees of freedom will be extremely accurate, but computationally less effi—
cient. Using these varying levels of efficiency in fitness evaluations can be used to increase
the computational efficiency of the iiGA. A unique dilemma occurs when comparing the
plane stress and three-dimensional finite models. The plane stress evaluation is efficient,
but not entirely accurate while the three-dimensional model is computationally hindersome
yet accurate. Therefore it would be satisfactory to first quickly approximate the response

of a flywheel with the plane stress evaluation to discard poor individuals and then fine tune

66

the design with a refined analysis. Figure 4.5 displays an iiGA topology that searches at
various levels of geometric resolution and uses various objective evaluations governed by
two factors. The first factor would be whether or not to evaluate with a plane stress assump-
tion. If the evaluation is not based on a plane stress analysis then the level of accuracy, or
the total number of degrees of freedom must be chosen for the three—dimensional evalua-
tion. Any migration of individuals requires re-evaluation of fitness due to the different lev-
els of accuracy in measuring fitness.

Further inspection of the iiGA topology in Figure 4.5 reveals that at each low level
of geometric resolution, both plane stress and three-dimensional evaluations exist. The is-
lands that evaluate fitness based on plane stress quickly evaluate individuals to be injected
into the islands that evaluate fitness with the a more refined analysis tool. Also, the islands
that evaluate fitness with the more refined three-dimensional model are allowed to evolve
independently, sharing only portions of the search space with other islands at the same level

of geometric resolution.

67

   
    
    
  
 
 
 
 
  

 

KEY: Evaluation

 
    

Tool
Three Rings 0 O Plane
(low resolution) Stress FEM
7. .
Six Rings ( 13 25 DOF)
(medium resolution) ./ 0 3-D FEM
(150 DOF)
Twelve Rings O.“ .~. 49:9} 3-D FEM
(high resolution) 333:: @533 (590 DOF)

@ 3-D FEM
(1,500 DOF)
. 3-D FEM

(10,000 DOF)

Twenty-Four Rings
(higher resolution)

Figure 4.5. Island injection GA topology with various
levels of fitness evaluation efficiency.

 

4.4.3 Sub-Fitness Functions

Each island could possibly evaluate fitness based on a sub-fitness function. For in-
stance, it was found earlier by Eby [9] that angular velocity was an important factor in the
optimization of Specific Energy Density (SED) of flywheels when limited to a plane stress
evaluation. Figure 4.6 displays an iiGA topology that searches at various levels of geomet-
ric resolution evaluating fitness with plane stress and three-dimensional analyses basing fit-
ness on angular velocity and SED respectively. This allows for efficient search at low
levels of geometric resolution with a plane stress fitness evaluation, which bases fitness on
a “sub-fitness” function, while also using a refined analysis to fine tune the flywheel at

higher levels of geometric resolution.

68

KEY: Evaluation Fitness
Tool Definition

0 Plane Angular
Stress FEM Velocity
(7, 13, 25 DOF)

0 3-D FEM SED
(190 DOF)

Three Rings 0
(low resolution)
Six Rings '

(medium resolution) Vi/ ‘ /////

Twelve Rings
(high resolution)

     

3-D FEM SED
(962 DOF)

@ 3-D FEM

(2,774 DOF)

. 3—D FEM
(12,272 DOF)

Figure 4.6. Island injection genetic algorithm topology
with multiple fitness definitions.

Twenty—Four Rings
(higher resolution)

 

4.4.4 Hybrid Injection Island Genetic Algorithms

GAs excel at searching on a global level to quickly find near-optimal solutions,
however they lack in the abilities to quickly climb hills in a local convex search space.
Threshold Accepting (TA) is an efficient hill climbing method and seems ideal to incorpo-
rate into the GA. The TA algorithm can be initiated with the best solution found by the GA,
allowing the TA to quickly climb from a near-optimal solution to the optimal solution. This
was done to climb hills that are close to the nearly globally optimal solutions found by the
GA. The new solution found with the TA approach can be reinserted into the iiGA replac-
ing the worst individual in the population.

A local search procedure was also developed to use with the iiGA. This local search

69

varied the thickness of the failure ring of the flywheel. This was done with the best individ-
ual at the end of each generation to discover improvements through slight variations in the
local failure ring thicknesses. Each ring has an inner and outer thickness, so the local search
first held the inner thickness constant and varied the outer thickness by increasing and de-
creasing by small increment. Next, the outer thickness was varied in a similar fashion. If a
better solution was found it replaced the worse solution in the population, and if only worse
solutions were found, then no individuals were replaced. Each local search contained four
evaluation calls working on a specific failure ring. Either a local or a TA technique could
be applied to an island, but not coincidentally. Figure 4.7 displays a hybrid iiGA topology

that benefits from local and TA methods.

KEY: Evaluation Fitness
_ _ .1. _ _ 1. _ _. .1 Tool Definition

0 Plane Angular
Stress FEM Velocity
(7, 13, 25 DOF)

0 3-D FEM SED
(190 DOF)

....... :O ::. '0
(low resolution)I

Six Rings |
(medium resolution]

  

Twelve Rings

3-D FEM SED
(high resolution)

(962 DOF)

3-D FEM

Twenty-Four Ringsl (2,774 DOF)

(higher resolution) I

I

(12,272 DOF)

 

Hybrid Local Search Plain GA Hybrid Threshold Accepting Algorithm

Figure 4.7 Typical hybrid island injection genetic algorithm toplogy.

70

4.5 Results and Discussion

First, the solid constant stress flywheel shape will be found and compared to the
plane stress closed fOrmed solution [1]. Next, the natural boundary condition at the end of
the plane stress constant stress flywheel will be made homogeneous to remove the assump-
tion that the optimal design is one with equal stresses throughout, an assumption that facil-
itates the analytical solution.

The global optimum solution for a relatively small design search space (2,097,152
designs) was found by enumerating (exhausting) the search space. The known global opti-
mum solution was then used to benchmark a set of algorithms that includes: simple Genetic
Algorithms (sGAs), Genetic Algorithms that compose topological “rings”, island injection
Genetic Algorithms (iiGAs), Threshold Accepting (TA) algorithms and hybrid Genetic Al—
gorithms. The hybrid Genetic Algorithms incorporate TA algorithms and or a local im-
provement method into either a 86A or an iiGA.

Next, multi-material solid and annular composite flywheels will be optimized with
no constraints on angular velocity or air gap growth. Finally, multi-material annular com-
posite flywheels will be optimized constraining angular velocity while minimizing “air

gap” growth.

4.5.1 Results of Optimized Solid Isotropic Flywheels

The first attempt to optimize the SED of solid isotropic flywheels with no con-
straints on angular velocity or “air gap” growth involved using a simple GA that evaluated
fitness with the plane stress finite element code to search for the well known constant stress

profile. First the fitness function was defined as SED with no constraints; the final result

71

is shown in Figure 4.8. The GA is simply placing mass at the outer edge of the flywheel to
increase the mass moment in the fitness function. The GA is exploiting the simplified plane
stress fitness evaluation. A stress concentration will be located at the point of inflection
near the end of the flywheel. This flywheel design violates the simplified plane stress as-
sumption that stress does not vary through the thickness.

The fitness function was next defined as angular velocity due to the fact that SED
is directly proportional to the square of the angular velocity of the flywheel. Figure 4.9
compares the analytical solution of a constant stress flywheel to the simple GA’s constant
stress flywheel. The difference in SED is less than 1%. All simple GA runs had a 75%
crossover rate, 1% mutation rate and initial population size of 300. The search for the con-
stant stress profile converged within 5 minutes on a Sun Sparc Ultra workstation. Figure
4.10 displays the evolution of the constant stress flywheel. Each fitness evaluation using
the plane stress evaluation lasted less than 0.001‘h of a second.

To solve the boundary value problem of a flywheel of constant stress a force (nat-
ural boundary condition) must be applied in the radial direction at the end of the flywheel.
To search for the optimal shape of a flywheel with stress free edges the natural boundary
condition at the end of the flywheel will be made homogeneous (zero). Again the fitness
was evaluated on angular velocity. Figure 4.1 1 compares the design of a solid isotropic fly-
wheel with stress free edges to the constant stress flywheel. The flywheel with stress free
edges has a 55% increase in SED when compared to the constant stress flywheel with the
axi-symmetric three-dimensional finite element code. The commonly used rotor shape of

the constant stress flywheel is not optimal when a flywheel has stress free edges.

72

4.5.2 Comparison of Optimization Methods Through Enumeration

To prepare a problem that is computationally enumerable, a relatively small search
space was designed. A six-ringed solid single isotropic material flywheel that used a coarse
spatial resolution (three bits per ring height) represented a design search space of
(2#bits)#Ri"gs+1= (23)7 =2,097,152. An algorithm was written to search exhaustively for the
global solution, taking 49 hours on a Sun Sparc Ultra with the axi-symmetric three-dimen-
sional finite element analysis containing 962 degrees of freedom in each analysis. Each
evaluation lasted approximately a 0.084th of a second. Now, comparisons for a set of opti-
mization techniques can be made with a known global solution. The optimization tech-
niques include simple Genetic Algorithms (sGAs), island injection Genetic Algorithms
(iiGAs), Threshold Accepting (TA) algorithms and hybrid iiGAs that incorporate TA and
local improvement methods. Statistics that pertain to the number of evaluations needed to
find the global solution are relevant measures of algorithm exploration efficiency. But after
exploring a search space with an iiGA using various levels of computationally efficient
evaluation tools, a better measure of algorithm performance can be conceived as the total
computation time (or effort) spent before discovering the global solution. Of course, this
does not measure the additional gains possible using the iiGA on multiple processors. Sta-
tistics for the total number of evaluation calls will also be given to satisfy the different ways
to evaluate algorithm exploration efficiency. For each optimization approach, five random-
ly generated initial populations were allowed to run for a total of two hours each to compute
the average time to find the global optimal solution. The fitness was based only on SED

for this simple problem.

The first optimization techniques to explore the enumerable search space were the

73

sGA and the TA methods. Figure 4.12 displays the normalized fitness as a function of time
for typical runs for the sGA, TA, and a hybrid GAs. Hybrid GAs incorporated either a TA
algorithm or local search method at the end of each generation. In a hybrid approach, the
TA algorithm was initiated with the best solution in the population and tried to improve the
solution within 10 calls, replacing the worse solution in the population. As previously de-
scribed, the local search method explored a small search space associated with the failure
ring. This was done by slightly varying the failure ring heights. In Figure 4.12, the sGA
alone quickly converges to a sub-optimal solution, never discovering the global solution in
any of the five runs. The TA algorithm starts at a lower fitness because it starts the problem
with a single, randomly generated flywheel. The TA algorithm quickly climbs to a sub-op-
timal solution. Figure 4.12 displays a typical run in which the TA algorithm seems to be
taking small steps in climbing to a local solution. Figure 4.12 also displays typical runs
where the hybrid iiGAs that incorporate either the local search method or TA algorithm dis-
cover solutions that are extremely close to the global optimal solution.

The iiGA topology shown in Figure 4.13 was run on a single computational node.
The iiGA parameters are shown in Table 4.5.1 demonstrate how the iiGA topology in Fig-
ure 4.15 can explore a design space with various levels of resolution and evaluation accu-
racy.‘ Islands 0-2 can quickly search the design space with larger populations due to the
efficient plane stress finite element code which bases fitness on angular velocity. These
islands will converge quickly due to the low level of geometric resolution and the efficient
plane stress finite element code. The iiGA completes a high number of generations (termed
generations per cycle) before solutions migrate to neighboring islands to further increase

the rate of convergence. A slightly higher crossover rate will encourage the iiGA to intro-

74

duce new designs into the population. Islands 3 and 4 search at a low level of geometric
resolution, evaluating individuals with an efficient axi-symmetric three-dimensional finite
element code (with a low number of degrees of freedom) basing fitness on SED. Islands
5-8 have lower population sizes with a lower number of generations per cycle and are in-
creasingly refined in geometric resolution. They measure fitness with various levels of ac-
curacy with an axi-symmetric three-dimensional finite element code. The iiGA completes
a lower numbers of generations before migrating solutions to neighboring islands. Islands
9-11 have the highest level of geometric resolution and three-dimensional finite element
evaluation accuracy. They also have the lowest number of generations completed before
migrating any solutions to neighboring islands. The iiGA topology and parameters were
specifically chosen to reduce the total number of expensive finite element evaluations
needed in the optimization of large scale engineering structures. Due to the fact the GAs
often require a large number of evaluations, decreasing the number of expensive finite el-
ement evaluations will help increase the computational efficiency of a typical GA. The TA
algorithm was called at most 10 times at the end of each generation. The TA algorithm was
initiated with the best solution in the population and replaced the worst solution in the pop-
ulation.

Figure 4.14 displays the normalized fitness (recomputed with the axi-symmetric
three-dimensional 962 degree of freedom analysis) of each island as a function of time.
Figure 4.14 displays the iiGA slightly “jumping” to higher values in fitness due to the in-
jection of high fitness low resolution individuals. The second and third islands (which are
at a low level of resolution) quickly converge (less than 10 seconds) to near global optimal

solutions. Figure 4.14 conveys the benefits of searching at various levels of resolution to

75

quickly find building blocks to be injected into islands of higher resolution. Unfortunately,
for this particular small search space, the difference in geometric resolution plays a small
role in the actual fitness. A near global optimal solution (within 2% of the global Optimal
solution) can be found using only three rings, causing difficulties to fully realize the impact
of searching at various geometric levels with the iiGA. Most of the work done by the iiGA
occurs quickly in the first few islands (at a low level of resolution). This makes it difficult
to realize the full benefits of an iiGA for such a simple problem. The global optimal solu-
tion is discovered first on the fourth island for this optimization run in less than 200 sec-
onds. The average time for the iiGA in Figure 4.13 to find the global optimal solution in
any island of high resolution measuring fitness with the axi-symmetric three-dimensional
finite element model (containing 962 degrees of freedom) in five runs was 768 seconds.
The sixth, seventh and eighth islands discover the global solution within 900 seconds in
Figure 4.14. These iiGA results illustrate the efficiency of searching at various levels of
resolution to quickly find building blocks to increase robustness and help decrease compu-

tational demands of a typical GA.

76

Table 4.5.1 Hybrid Island Injection Genetic Algorithm Parameters for the
Enumerable Coarse Design Space of a Solid Isotropic Flywheel.

 

 

 

 

 

 

 

 

 

 

 

Island Pop 6:8 Crossover Mutation Analysis Fitness #
# Size €3,016 Probability Probability Tool Definition Mi grated
0 300 12 75% 1% Plane (1) 3

Stress
1 300 12 75% 1% Plane 0) 3
Stress
2 200 6 70% 1% 3-D SED 3
3 200 6 70% 1% 3-D SED 3
4 200 4 65% 1% 3-D SED 3
5 200 4 65% 1% 3-D . SED 3
6 100 2 60% 1% 3-D SED 3
7 100 2 60% 1% 3-D SED 3
8 100 2 60% 1% 3-D SED 3

 

 

 

 

 

 

 

 

 

 

The iiGA topology shown in Figure 4.15 was run on a single computational node
that combined a TA algorithm to create a hybrid iiGA. iiGA parameters are shown in Table
4.5.1. Figure 4.16 displays the normalized fitness (recomputed with the axi-symmetric
three-dimensional 962 degree of freedom analysis) for a typical optimization run for the hy-
brid iiGA for each island as a function of time. The second and third islands quickly con-
verge to near global optimum solutions. The global optimum is discovered first on the fifth
island in less than 200 seconds. The average time for the iiGA in Figure 4.15 to find the

global optimum in any island of high resolution measuring fitness with the axi-symmetric

77

three-dimensional finite element model (containing 962 degrees of freedom) in five runs
was 674 seconds. The sixth, seventh and eight islands converge to the global optimum
within 500 seconds for this optimization run. When comparing Figure 4.14 to 4.16, it can
be seen that the hybrid iiGA climbs from local optimum solutions that are near by the global
optimum to the global optimum through help from the TA algorithm. Also when compar-
ing Figure 4.14 to 4.16, a the hybrid iiGA makes noticeably smaller steps to the global op-
timum when compared to the iiGA. These hybrid iiGA results demonstrate how a local hill
climber such as a TA algorithm can help the iiGA to find a global optimum in an efficient
manner.

The iiGA topology shown in Figure 4.17 was run on a single computational node
that combined a local search procedure to create a hybrid iiGA. iiGA parameters are shown
in Table 4.5.1. The local search procedure was called at most 10 times at the end of each
generation. The local search procedure took the best solution in the population and replaced
the worst solution in the population if a better solution was discovered. Figure 4.18 dis-
plays the normalized fitness (recomputed with the axi-symmetric three-dimensional 962
degree of freedom analysis) for the hybrid iiGA for each island as a function of time. The
second and third islands quickly converge to near global optimums. The global optimum
is discovered first on the fifth island in less than 200 seconds. The average time for the
hybrid iiGA in figure 4.17 to find the global optimum in all islands of high resolution mea-
suring fitness with the axi-symmetric three-dimensional finite element model (containing
962 degrees of freedom) in five runs was 715 seconds. The sixth, seventh and eight islands
converge to the global solution within 750 seconds for this optimization run. When com-

paring Figure 4.14 to 4.18, it can be seen that the hybrid iiGA that uses a local search pro-

78

cedure to help climb from local solutions that are near by the global solution to the global
solution. These hybrid iiGA results demonstrates how a local search procedure can help the
iiGA find the global optimum in an efficient manner.

The iiGA topology shown in Figure 4.19 was run on a single computational node
that combined local search procedure and TA algorithms to create a hybrid iiGA. iiGA pa-
rameters are shown in Table 4.5.1. The TA algorithm and local search procedure was
called at most 10 times at the end of each generation in each island. The local search pro-
cedure took the best solution in the population and replaced the worst solution in the pop-
ulation if a better solution was discovered. Figure 4.20 displays the normalized fitness
(recomputed with the axi-symmetric three-dimensional 962 degree of freedom analysis) for
the hybrid iiGA for each island as a function of time. The second and third islands quickly
converge to near global optimums. The global optimum is discovered first on the fifth is-
land in less than 200 seconds. The average time for the iiGA in Figure 4.19 to find the glo-
bal optimum in any island of high resolution measuring fitness with the axi-symmetric
three-dimensional finite element model (containing 962 degrees of freedom) in five runs
was 417 seconds.

The GA ring topology shown in Figure 4.21 was run on a single computational
node. GA parameters were as follows: 75% crossover and 1% mutation. Figure 4.22 dis-
plays the normalized fitness computed with the axi-symmetric three-dimensional 962 de-
gree of freedom analysis for the GA ring topology for each island as a function of time. All
the islands slowly converged to a near optimal solution. The algorithm was allowed to run
5 hours, never discovering the global optimum.

When comparing all the optimization techniques presented above, some general

79

conclusions can be drawn. The best solution discovered with the TA algorithm had a lower
fitness when compared to all other approaches for all optimization runs. The sGA approach
was improved in terms of computational efficiency when combining with either the TA al-
gorithm or the local search method. The iiGA approach outperformed the TA algorithm
and all single island and ring topology GA approaches in terms of a lower computational
effort to locate the global optimum. All hybrid iiGA approaches required less computa-
tional effort (averaged over five runs) compared to other iiGA approaches. The hybrid iiGA
that used a combination of the TA algorithm and the local search procedure required less
computational effort (when averaged over five runs) compared to a hybrid iiGA that used
either the TA algorithm or the local search procedure independently. It appears the hybrid
iiGA that incorporates both TA and the local search method is more robust in terms of ef-
ficient exploration of the search space, at least for the problem considered. Table 4.5.2 con-
tains the statistics on the average number of calls made over the five runs. The hybrid iiGA
that used the TA algorithm and local search method searched less than 5% of the entire
search space, taking less than 0.5% of the time to enumerate the entire search space, mea-
suring more than half of the evaluations with the plane stress finite element model to find

the global optimum.

Table 4.5.2. Statistics on the Total Number of Evaluation Calls
made for Various Optimization Approaches in Average 2-Hour Run.

 

 

 

 

 

 

 

 

 

 

Optimization Total FEM $31538 StrglanEM Total TA Total Local
Approach Calls Calls FEM Calls FEM Calls
TA 12,201 0 12,201 12,201 0
sGA 26,877 0 26,877 0 0
SGA plus TA 19,902 0 19,902 980 0
SGA plus 23,543 0 23,543 0 3,110
Local
iiGA 122,900 52,828 70,072 0 O
iiGA plus 117,102 49,051 68,051 0 12,150
Local
iiGA plus 1 11,521 46,051 65,470 11,340 0
TA
iiGA plus TA 98,853 43,771 55,082 3,880 7,800
and Local

 

 

 

 

 

 

 

 

4.5.3 Optimization of Multi-Material Composite Flywheels

Finally, optimum multi-material composite flywheel shapes will be sought. First,
solid multi—material flywheel shapes will be optimized with no constraints on maximum
angular velocity with an iiGA, hybrid iiGA and topological “ring” GA. Next, annular
multi-material composite flywheels will be optimized while controlling the maximum al-
lowable angular velocity while reducing “air gap” growth. It was found that the iiGA and
hybrid iiGA search was extremely efficient when compared to a “ring” topology GA. Also,

it was found that the hybrid iiGA found solutions with slightly higher fitness values when

81

compared to the iiGA or “ring” topology GA.

It was discovered in previous iiGA optimization runs that an extremely refined
three-dimensional finite element mesh was needed to avoid “artifacts” in the solutions. The
solutions that contained “artifacts” of the evaluation typically placed a large amount of ma-
terial near the outer radius of the flywheel, creating a sharp inflection point near the ring
interfaces. Figure 4.23 displays a solution that is an artifact of the analysis due a stress con-
centration that would be associated with the sharp inflection point. The three-dimensional
finite mesh contained 2,482 degrees of freedom. This solution was re-evaluated multiple
times with an increasingly refined three-dimensional finite element model until the maxi-
mum stress in the model ceased to increase by more than 5%. The three-dimensional finite
element model required 12,272 degrees of freedom to converge. As expected, the origin of
failure in the converged solution was located near the sharp point of inflection. Further re-
finements in the three-dimensional model would be fruitless due to the nature of the specific
flywheel geometry. In other words, a further refinement in the analysis would only accen-
tuate the stress concentration located near the sharp point of ring inflection. The maximum
stress in the flywheel will be dominated by the stress concentration. The iiGA needs only
realize that designs of this nature are not optimal. This level of refinement helped the GA
realize that flywheel designs of this sort are not optimal.

The optimization of solid and annular composite multi-material flywheels will be
approached with the hybrid iiGA topology shown in Figure 4.7 using 10 bits per ring height
in spatial resolution while ranging the geometric resolution from a 3 to 24 rings. This type
of hybrid iiGA topology performed better on average compared to other hybrid iiGA topol-

ogies in the previous study of the coarsely designed solid isotrOpic flywheel. Notice islands

82

9-11 are the most refined in geometric resolution and are evaluated with the three-dimen-
sional finite element code with 12,272 degrees of freedom to avoid “artifacts” in the final
solution.

iiGA parameters are contained in Table 4.5.3. The iiGA parameters were chosen to
help decrease computational costs as previously outlined in the optimization of an enumer-
able search space of a solid isotropic flywheel. Table 4.3.3 contains one quasi-isotropic ma-
terial, but the GA can choose the direction of “strength” in either the radial or tangential
direction. Also, the GA was given the choice of two isotropic materials (Steel and Alumi-
num, Table 4.3.1). It is hypothesized that the GA could place materials that have a higher
density at the outer radius of the flywheel to increase the mass moment of inertia term in
the definition of SED. Next, a solid multi-material flywheel will be optimized with respect
to SED. I

Figure 4.24 shows the “raw” fitness for solid multi-material flywheels as a function
of time for each island in an iiGA depicted in Figure 4.7. iiGA parameters are contained in
Table 4.5.3. The “raw” fitness is the actual fitness measured by each island’s specific finite
element evaluation. Islands 0-2 measure “raw” fitness with an approximate but efficient
evaluation based on angular velocity. The axi-symmetric plane stress finite element code
predicts fitness accurately for flywheel shapes that have small ring gradients, producing re-
sults in Figure 4.24 that are excessively optimistic. The axi-symmetric three-dimensional
finite element code evaluated islands 3-8 with a reduced number of degrees of freedom
when compared to islands 9-11. Therefore we can expect some discrepancy in fitness val-
ues for islands 3-8 when re-evaluating fitness with the most refined axi-symmetric three-

dimensional finite element model. From this point forward, all plots that display fitness as

83

a function of time will be recomputed with the most accurate evaluation used in the run.

Figure 4.25 displays the fitness of solid multi-material flywheels as a function of
. time (re-evaluated fitness measured at the highest level accuracy with the three-dimensional
axi-symmetric finite element model containing 12,272 degrees of freedom). Hybrid iiGA
topological structure is shown in Figure 4.7 while all parameters are contained in Table
4.5.3. In Figure 4.7, islands 08 evaluated individuals with either an axi-symmetric plane
stress or three-dimensional finite element code that is less accurate when compared to the
axi-symmetric three-dimensional finite element code that evaluated fitness with 12,272 de-
grees of freedom in islands 9-11. The inaccuracies of the evaluation in islands 0-8 are
shown in Figure 4.25. Figure 4.25 displays an expected response; islands 0-8 initially find
good solutions but begin to find worse solutions as time progresses. These solutions con-
tain building blocks that are used to help evolve islands at higher levels of resolution
through injection and therefore cannot be discarded even though they have a low fitness
when evaluated with the most accurate finite element model. We can expect, but cannot
discard what appears to be “noise” in Figure 4.25 for islands of lower resolution. This ex-
pectation would be a snafu for islands that were initially evaluated using the finite element
model with the highest degree of accuracy. An explanation of the origin of noise and noise
effects are given next.

The iiGA explores the search space defined by the fitnesses of the individuals in a
population through successive generations. Noise occurs when the iiGA cannot decipher
the differences between a solution that does or does not violate an assumption of a fitness
evaluation (for example a plane stress finite element evaluation). If a high fitness is asso-

ciated with solutions that violate a fitness evaluation then the iiGA will sooner or later ex-

84

ploit this to improve the existing solutions in the population. Initially, the intended fitness
definition dominates the iiGA search procedure for islands that evaluate fitness with an ef—
ficient, less accurate finite element evaluation. After a certain amount of time this dominat-
ing effect lessens, allowing the iiGA to explore areas of the search space that have high
actual fitness values, but violate evaluation assumptions. The intended fitness can be con-
ceived as the fitness the engineer or user ideally defines, assuming that the fitness can be
evaluated absolutely in closed form. The actual fitness can be thought of as the fitness that
can also be a function of violations in the fitness evaluation. The iiGA searches the actual
fitness definition, often finding artifacts in the final solution. This can force an engineer or
the user to use analyses that are more refined and computationally demanding. Any in-
crease in computational time will amplify the total time required to effectively explore the
search space. An analysis refinement can severely limit the practical application of a typical
GA. This is truly one area where the iiGA can shine by searching at various levels of com-
putational efficiency to quickly explore the search space to find potentially fruitful areas.
These fruitful areas can be re-evaluated at a higher level of accuracy to double check the
solution.

Figure 4.26 is Figure 4.25 “zoomed” in to help accent the effects of searching at var-
ious levels of resolution while exploiting less accurate, computationally cheap finite ele-
ment evaluations with the hybrid iiGA shown in Figure 4.7. All islands that did not evaluate
fitness with the most refined three-dimensional finite element evaluation contain noise as
previously discussed. Islands 0, 3 and 4 begin with high fitness values, initially finding fit
building blocks due to a low flywheel geometric resolution (three rings) and an extremely

efficient finite element evaluation. Island 0, 1 and 2 were evaluating fitness with a plane

85

stress finite element evaluation basing fitness on angular velocity. Figure 4.26 shows that
basing fitness on a “sub-fitness” function was fruitful. Islands 5 and 6 begin with a high
initial fitness and initially converge slower to a higher fitness when compared to islands 0
through 4. Islands 7 and 8 initially converge slower but converge to a higher fitness when
compared to islands 5 and 6. After a period of time the iiGA begins to find worse solutions
due to inaccuracies in the finite evaluation in islands 0-8. Islands 9, 10 and 11 begin with
the lowest inital fitness, converging the slowest when compared to all other islands while
evolving to the most fit individuals. This was due to the fashion that the hybrid iiGA
searched.

The technique begins by first quickly finding building blocks at a low (three rings)
then medium-low (six rings) then medium (twelve rings) and high (24 rings) levels of res-
olutions with various evaluation tools. Islands that do not search at the highest level of res-
olution use two different evaluation tools and fitness definitions. The fitness definitions for
these islands were measured with an axi-symmetric plane stress finite element code based
on angular velocity or an axi-symmetric three-dimensional finite element code basing fit-
ness on SED. Next, the iiGA migrates solutions containing good building blocks in an hi-
erarchical order (from lower to higher level of resolution). Each island receives the
solutions that contain building blocks in the next cycle of the algorithm. Figure 4.26 shows
that islands 9, 10 and 11 evolve individuals with a higher level of fitness (compared to all
other islands) due to the increase in geometric resolution.

Figure 4.27 shows the best solid multi-material composite flywheel found at each
level of geometric refinement for the axi-symmetric plane stress (basing fitness as a func-

tion of angular velocity) and the three-dimensional finite element evaluation (basing fitness

86

on SED). The plane stress evaluations allow the iiGA to efficiently explore the search space
to discover fruitful solutions by sacrificing evaluation accuracy. Often, many optimization
approaches use a similar level of accuracy to approximate expensive evaluations, giving
rise to solutions that are artifacts of the simplified analysis. Comparing the best flywheel
found by the axi-symmetric plane stress finite element code to the three-dimensional finite
element code in Figure 4.27 (at each level of resolution) displays solutions that are vividly
artifacts of a simplified analysis tool (plane stress approximation). Completely dismissing
the solutions found by the plane stress evaluation is a mistake because they are just exag-
gerated shapes that have helped evolve the shapes found with the refined axi-symmetric
three-dimensional finite element analysis. The latter statement can be better understood by
re-examining Figure 4.25, where the best flywheel found at each level of resolution and ac-
curacy is displayed as a function of time. Figure 4.27 also conveys how the iiGA “fine
tunes” the solutions in islands of higher resolution and higher accuracy. The iiGA places
only the E-Glass material throughout out the flywheel at various orientations, never using
the more dense isotropic material. Evidently, the iiGA did not find that exploiting the mass
moment of inertia (with either aluminum or steel) of an unconstrained solid composite fly-
wheel to increase the SED when compared to the associated sacrifices that would have been
made in the failure angular velocity.

Figures 4.28 and 4.29 display the fitness of an annular multi-material flywheel as a
function of time (re-evaluated fitness measured at the highest level accuracy with the three-
dimensional axi-symmetric finite element model containing 12,272 degrees of freedom) for
the iiGA and hybrid iiGA displayed in Figures 4.6 and 4.7, respectively. All iiGA param-

eters are contained in Table 4.5.3. Figure 4.30 and Figure 4.31 are Figures 4.28 and 4.29

87

that are “zoomed” in to help emphasize the effects of searching at various levels of resolu-
tion while exploiting less accurate, computationally cheap finite element evaluations.

As previously stated regarding the iiGA results for a solid multi-material composite
flywheel, islands with a coarse geometric resolution initially find solutions containing
building blocks extremely fast. This is due to the combination of a small search space and
computationally efficient finite element evaluations. Islands with a refined geometric reso-
lution converge slower due to a larger search space and a computationally demanding three-
dimensional finite element model. Islands at a lower resolution quickly find solutions that
contain good building blocks that were injected in islands in an hierarchical fashion (from
lower to higher level of resolution). Islands with higher levels of geometric resolution con-
verge slower to solutions that have an increase in fitness (when compared to islands of low-

er geometric resolution) due to the increase in search space.

88

Table 4.5.3 Island Injection Genetic Algorithm Parameters.

 

Gens

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Island Pop Prob. Prob. Analysis Fitness # #
# Size 6;:th Xover Mutation Tool Def. Migrated Rings
0 300 12 75% 1% Plane Stress (a) Best 3
l 300 12 75% 1% Plane Stress (0 Best 6
2 300 12 70% 1% Plane Stress to Best 12
3 200 8 70% 1% 3-D SED Best 3
4 200 8 65% 1% 3-D SED Best 3
5 200 6 65% 1% 3-D SED Best 6
6 100 6 60% 1% 3-D SED Best 6
7 100 4 60% 1% 3-D SED Best 12
8 100 4 60% 1% 3-D SED Best 12
9 100 2 60% l % 3-D SED Best 24
10 100 2 60% 1% 3-D SED Best 24
1 1 100 2 60% 1% 3-D SED Best 24
Table 4.5.4 Genetic Algorithm Topological “Ring” Parameters.
Island Pop Gens Prob. Prob. Analysis Fitness # #
Size C336 Xover Mutation Tool Def. Mi grated Rings
All 105 3 65% 1 % 3—D SED Best 24

 

 

 

 

 

 

 

 

 

 

 

89

Figure 4.32 depicts 20 islands that form a topological “ring” GA used to obtain the
fitness as a function of time presented in Figure 4.33 for an annular multi-material compos-
ite flywheel. All islands operated on the identical GA parameters given in Table 4.5.4 mea-
suring fitness with the axi-symmetric three-dimensional finite element model containing
12,272 degrees of freedom. The topological “ring” GA contained the same total number
of individuals found in the hybrid iiGA, equally divided amongst the islands. Figure 4.33
shows that each island slowly finds better solutions in the “ring” topology GA, taking near-
ly six days to converge.

Figure 4.34 compares the results obtained from the iiGA, hybrid iiGA and topolog-
ical “ring” GA. For the iiGA approach, Figure 4.25 displays that the best solution found in
each island (that evaluated fitness with the highest level of finite element accuracy) is nearly
identical as a function of time. This is also true for the hybrid iiGA and topological “ring”
GA. For these reasons, the results shown in Figure 4.34 are drawn from typical islands that
evaluated fitness with the most refined three-dimensional finite element model for both the
iiGA, hybrid iiGA and topological “ring” GA. Both the iiGA and the hybrid iiGA converge
over five times faster when compared to the topological “ring” GA. The hybrid iiGA also
converges to a higher fitness compared to both the iiGA and topological “ring” GA. Table
4.5.5 contains statistics for the total number of evaluation calls made during each GA opti-
mization approach. Since each topological “ring” GA run takes nearly six days, no attempt
was made to compute averages over a set of independent runs for this optimization ap-
proach. However, a second set of independent runs for both the iiGA and hybrid iiGA were
done for two reasons: they took just over one day to converge and to display the fact that

these results were not biased by the initial population. The second runs produced nearly

90

identical results in terms of final fitness values and total time until convergence for both the
iiGA and hybrid iiGA.

Figure 4.35 compares the best solution found the iiGA, hybrid iiGA and the topo-
logical “ring” GA. They are all similar in shape, with some slight variations. The best so-
lution found by the “ring” topology GA is not completely smooth in shape, placing more
material at the end of the flywheel to increase the mass moment of inertia. All solutions
contain similar material placement near inner radius of the flywheel. This can be expected
since the maximum stress found typically in an annular flywheel is located at the inner ra-
dius in the tangential direction. All the optimization approaches placed the strongest direc-
tion of E-Glass epoxy in the tangential direction. The hybrid iiGA and iiGA found
solutions that are very similar in shape, but vary slightly in material placement away from
the inner radius of the flywheel. The hybrid iiGA found a solution that has the least amount

of mass placed near the end of the flywheel.

Table 4.5.5. Statistics on the Total Number of Evaluation Calls for
iiGA, Hybrid iiGA and Topological “Ring” GA Optimization Approaches.

 

 

 

 

 

 

Optimization Total FEM 11.123152138 StrgiznlgEM Total TA Total Local
Approach Calls Calls FEM Calls FEM Calls
Topological 12,201 12,201 0 0 0
“Ring” GA .
iiGA 1,244,717 383,074 861,643 0 0
Hybrid iiGA 861,889 306,544 555,345 43,750 34,980
(local and TA)

 

 

 

 

 

 

 

 

91

To use the iiGA to design flywheels that have a reasonable angular velocity, a “tar-
get” angular velocity of a flywheel must be specified. The iiGA should also penalize fly-
wheels that have a large amount of “air gap” growth to discover “safer” designs. A lower
“target” angular velocity can be specified by the designer to help control the high angular
velocities associated with the solutions found in an unconstrained optimization approach of
a flywheel (as just previously presented). A new “target” angular velocity was determined
by reducing the angular velocity that was found in best solution in the previous uncon-
strained optimization problem. A reduction of 25% in angular velocity was chosen so the
iiGA could place materials and form flywheel shapes while still maintaining a reasonable
SED. For the constrained problem, the fitness was based on equation 4.2. In equation 4.2,
the SED term was normalized with SED results ascertained previously in the unconstrained
optimization approach. This was done similarly to normalize the “air gap” growth at fail-
ure. The angular velocity was normalized with the “target” velocity. Flywheels were pe-
nalized only if the angular velocity surpassed the “target” angular velocity. This means the
iiGA can avoid all penalties associated with high angular velocities by maintaining fly-
wheels that have angular velocities below the “target” value.

For the constrained optimization of multi-material annular composite flywheels a a
slight variation in the hybrid iiGA topology is shown in Figure 4.36. The only difference
from the previous hybrid iiGA topology is the elimination of the plane stress evaluations.
Substituted instead of a plane stress evaluation was the three-dimensional evaluation. This
change in approach was made due to the fact that the optimization of an annular multi-ma-
terial composite flywheel requires control of the maximum angular velocity of a flywheel,

not simply maximizing the “sub-fitness” function (angular velocity). The islands that pre-

92

viously evaluated with the efficient plane stress evaluation were replaced with three-dimen-
sional evaluations that were reduced in the total number of degrees of freedom in the
analysis.

Figure 4.37 displays the fitness of the constrained optimization of an annular multi-
material flywheel as a function of time (re-evaluated fitness measured at the highest level
of accuracy with the three-dimensional axi-symmetric finite element model containing
12,272 degrees of freedom). Hybrid iiGA parameters are contained in Table 4.5.3. Figure
4.38 shows the hybrid iiGA quickly finding solutions that contain building blocks at various
levels of geometric resolution. Some “noise” is present due to the multiple levels in the
three-dimensional finite element fitness evaluation.

Figure 4.39 shows the best solution found a various levels of geometric resolution
for the constrained optimization of multi-material annular flywheels. All show interesting
patterns in shape an material placement. The next few lines include reasonings for specific
material placement throughout the flywheel, starting at the inner radius and ending at the
outer radius. To overcome penalties due to large “air gap” growth the iiGA placed the stiff-
er material (aluminum) near the inner portion of the flywheel. To maintain a relatively high
SED the hybrid iiGA next placed composite material in the flywheel with an orientation of
strength in the tangential direction. Next, the hybrid iiGA placed composite material in the
flywheel with an orientation of strength in the radial direction. The angular velocity con-
straint forced the iiGA to design a flywheel that placed a higher density material (alumi-
num) at the outer radius to increase the mass moment of inertia of the flywheel.

The angular velocity in the best solution found at the highest level of geometric res-

olution was extremely close to the “target” angular velocity. The “air gap” growth found

93

in the best solution found at the highest level of geometric resolution in the constrained op-
timization problem was cut in half when compared to the unconstrained single criteria op-
timization approach. The SED for the best solution found at the highest level of geometric
resolution in the constrained optimization problem was 21.6% lower than the unconstrained
single criteria optimization approach. This was due to the constraint on angular velocity,

which is squared in the fitness definition.

Thickness

94

 

0 .25

0.20 '-

0.15

V‘T I Y T Y V

0.10 r

0 .05

 

 

A

 

A l A_.

 

0 .00
0 .00

0.10 0.15 0.20 0.25
Radius

0.05

Figure 4.8. Simple genetic algorithm exploitation of plane stress evaluation.

Thickness

 

0.25

0.20 r

0.10 '-

0 .05

 

0 .00

 

0.00

 

0—-0 Theoretical Canahnt Shea
B—El GA Constant Streee

 

 

 

 

AA. 1

0.20

A A A _A_ A L

0.05
Radius

Figure 4.9. Theoretical constant stress profile compared to
simple genetic algorithm constant stress profile.

95

 

Generation 0 Generation 34

   

_ Generation 13

. Generation 38

he)

 
 

Generation 21 Generation 47
Generation 28 Generation 650

Figure 4.10. The simple genetic algorithm evolution of a constant stress flywheel.

Normalized Fitness

96

 

 

 

 

 

 

 

 

 

 

 

 

  

 

 

 

1000.0

0.25 - f . , . - - .
0—0 GA Non-Constant Stress
0 20 13—8 Theoretical Constant 91058 1
0.15 -
ill . .
S .
is .
«C t. q
.—
010 -
0.05
0 '00 a . . a I . A g g I A - L - l . A . l - . A -
0 .00 0 .05 0.10 0 .15 0.20 025
Radius
Figure 4.11. Simple genetic algorithm non-constant stress
profile compared to the theoretical constant stress profile.
1 U ' r T 1 ' T n f 1 fl
:3 / V {FF-Pf v
av > ‘d— a" f,
0.8 ~ ' -
_.,_f H SGA 1
G—El sGA and TA
0 6 _| ————-:~TA _
' V—V sGA and Local
04 - -
0.2 - -
DU 1 a l ‘_ l A 1
0.0 200 .0 400.0 600.0 800.0
Time (seconds)

Figure 4.12. Simple genetic algorithm compared to TA
algorithm and hybrid genetic algorithms for a solid isotropic
flywheel in an enumerable search space.

97

  

   

 

KEY: Evaluation Fitness
Tool Definition

   

 

    
  

      
   
 
   
   

 

 

Normalized Fitness

 

 

 

Three Rings
(low resolution) “0 Plane Angular
Stress FEM Velocity
Six Rings (7, 13 DOF)
(high resolution)
0 3-D FEM SED
(56 DOF)
Six Rings
(high resolution) 3-D FEM SED
(350 DOF)
. 3-D FEM SED
(962 DOF)
Figure 4.13. Island injection GA topology.
1.0 '1 I i T r
‘ / __ .. —~—~
/ ‘ fl , H--.
0.8 ‘
G—O Island 0
0.6 - B—El Island 1 -
v--‘.~. Island 2
H Island 3
Island 4
' Island 6
-+——+- Island 7
Island 8
0.2 l- ..
00 a l a I - l . l .
0 .0 200 .0 400 .0 600.0 800.0 1000.0

Time (seconds)

Figure 4.14. Single computational node island injection
genetic algorithm results for a solid isotropic flywheel
in an enumerable search space.

 

KEY: Evaluation Fitness

     
 

(350 DOF)

l' "" — "l' — _ 1" —’ — 'l . .
. I I I I Tool Definltlon
Three Rlngs ® .
(low resolution) I I O Plane Angular
. . I I I I Stress FEM Velocity
Slx Rings 1 I I I (7, 13 DOF)
(high resolution) I I
I I I | 0 3-D FEM SED
l I (56 DOF)
Six Rings I I
(high resolution) I 1 3-1) FEM SED
|
I

|
'--—-L—--'----' .3-DFEM
(962DOF)

\/

Threshold Accepting Algorithm

 

Figure 4.15. Hybrid island injection GA topology (TA algorithm).

 

 

 

 

e——e Island 0
H Island 1
«‘«-—————.,; Island 2 d
0'8 I A—Alsland 3

Island 4
"5‘ Island 5

Island 8
0-4 ’ -l——-I-|sland 7 ‘

- .. Island 8

Normalized Fitness

 

 

 

0.0 A l . I a I . .
0.0 200 .0 400.0 600.0 800.0 1000.0

Time (seconds)

Figure 4.16. Single computational node hybrid island
injection GA (TA algorithm) results for a solid isotropic
flywheel in an enumerable search space.

Tool Definition

     

r. ._ _ .I. _ ._ .I. _ ._ 1 KEY: Evaluation Fitness
(9”
(low resolution) I O Plane Angular
. . ' ' ' Stress FEM Velocity
Slx Rlngs I I I (7, 13 DOF)
(high resolution) I
I 0 3-D FEM SED
I 56 DOF)
Six Rings I (
(high resolution) 1 3-1) FEM 513])
I (350 DOF)
I
L

. 3-D FEM
(962 DOF)

 

Local Search Method

Figure 4.17. Hybrid island injection GA topology (local search method).

 

 

 

Normalized Proteus

 

 

 

G—O Island U
H Island 1
____._.§. Island 2 _
“I ' Hlsland 3
Island 4
- Island 5
Island 8
0.4 . -I—I-Island 7 ‘
— Island 8
r
02 l- -
00 . l . J . l a I .
0 .0 200 .0 400 .0 600.0 800.0 1000.0

Time (seconds)
Figure 4.18. Single computational node hybrid island

injection GA (local search method) results for a solid
isotropic flywheel in an enumerable search space.

100

     
     

 

KEY: Evaluation Fitness
Tool Definition

0 Plane Angular
Stress FEM Velocity

§®IGIO'

     
 

 

Three Rings
(low resolution)I

  

  
  

  
   

Six Rings (7, 13 DOF)
(high resolution)
0 3-D FEM SED
(56 DOF)

 

(high resolution) 3-D FEM SED

(350 DOF)

 
   
 
  
 

I
|
I
I
Six Rings I
I
I
I
L

  

3-D FEM
(962 DOF)

  

Hybrid Local Search None Hybrid Thershold Accepting Algorithm

Figure 4.19 Hybrid island injection GA topology (local search and TA).

 

 

 

Normalized Fitness

 

 

 

1'0 -f. _-
0.8 .
G——€) Island 0
B—EI Island 1
0 6 _ <:~~-~<.~~ Island 2 4
' H Island 3
Island 4 I
" Island 5
Island 6
0.4 - >e——+< Island 7 ‘
< Island 8
0.2 a .
00 L t . l . l a l I
0.0 200.0 400.0 600.0 600.0 1000.0

Time (seconds)
Figure 4.20. Single computational node hybrid island
injection. GA (local search and TA) results for a solid
isotropic flywheel in an enumerable search space.

Normalized Thickness

101

 

 

KEY: Evaluation Fitness
Tool Definition

. 3-D FEM SED
(962 DOF)

Figure 4.21. Twelve “ring” GA topology.

 

 

 

 

 

 

1.0 , . A. . , ,
A, up, ~t':J-'»"+-""* ‘El —-‘--
H- _‘__,_.E-—— H {.7— ‘M A
foe—M “F.
0.8 -
o/
. * ” G—{llsland 0
H Island 1
0 6 '.=.~—-—-.~ Island 3 _
' H Island 4
Island 5
I ~—-—-“ Island 6
Island 7
0-4 ” +——+Island s -
0.2 - a
L
00 l A 1 . 1 L . l A
0.0 1000.0 2000.0 3000.0 4000.0 5000.0 “00.0

Time (seconds)

Figure 4. 22. Single computational node topological “ring” GA
results for a solid isotropic flywheel 1n an enumerable search

space.

102

enter

Location of Stress Concentration

 

 

oooooooooooooooooooooooo
Figure 4.23. Demonstration of a flywheel that is an “artifact” for a
three-dimensional finite element model (containing 2,482 degrees
of freedom).

 

“- rr 7” _

G—e Island 0
El—El Island 1
«m Island 2
H Island 3 ‘

Island 4

:~-~~ Island 5

Island 6

xz—x Island 7

Island 8

G—6 Island 8
-‘ Island 10
-—‘ Island 11

 

 

on I I l l l
0.0 30000.0 60000.0 90000.0 1200000 1500000

Time (seconds)
Figure 4.24. Hybrid iiGA “raw” fitness for each island as a function
of time for a solid multi-material composite flywheel. Islands 0-2:
plane stress fern evaluation. Islands 3-11: 3—D FEM evaluation.

 

103

 

 

 

 

 

 

 

40.0 - _. -
.. . . W; x xwml,‘ \
f [I It?" X—XK 1
30.0111; ' . G—Olslando
'- . B—EI Island 1
(If-IP- .. —-———~ . 2.... Island2
33.1 WM ' A—A Island 3
gII'.“ Isl and 4
20-0 :1 ‘ - ----~—{~\ Island 5
595 . Island 0
E. .4." "i , x—x Island 7
.2‘.‘ ‘ -——--—« Island a
I e——o Island 9
I” . -- lsland10
t: .5: -——- Island 11
El ”‘ ' I
0.0 ‘ l ‘ ‘ I ‘ ‘ I - L l - a L a
0.0 30000.0 60000.0 900000 1200000 1500000

Time (seconds)

Figure 4.25. Hybrid iiGA fitness re-evaluated for each island
as a function of time. All re-evaluations performed with 3-D
FEM containing 12,272 degrees of freedom.

104

 

40.0 . . ‘

G—O Island 0
[El—El Island 1
9—--~--< lsI and 2
A—A Island 3

Island 4
I“~M~~~~i:> Island 5

Island 8
x——>< Island 7
,. ,. ~ v Island 8
151—9 Island 9
Wm...» Island 10

 

 

Fitness

 

 

 

 

 

0.0 2000.0 4000.0 6000.0 8000 .0 10000.0
Tlme (seconds)

Figure 4.26. Hybrid iiGA quickly finding building blocks at low
levels of geometric resolution to inject into islands of higher
resolution.

105

     
 
  
 
  
         
  

KEY: Evaluation Fitness
Tool Definition

0 Plane Angular
Stress FEM Velocity
(7, 13, 25 DOF)
0 3—D FEM SED
(190 DOF)

3-D FEM SED
(962 DOF)

@ 3-D FEM SED
(2,774 DOF)

. 3-D FEM SED
I (12,272 DOF) I

Material

I E-Glass Epoxy [0°]

 

 

 

 

   

Fitness .
Definition Angular Velocrty

Evaluatio 3-1) Plane Stress
Tool

? .- E-Glass Epoxy [90°]

 

Figure 4.27. Best flywheel found at each level of resolution
and fitness definition.

 

 

 

0.0

106

 

 

 

 

 

 

 

 

 

 

 

 

 

' t l - l
: ’ JMM .
-._.:.._..._-_.._.1,~.*......‘.,-...,.‘+,.C_._ H __ ._ ,, ~ -.. 4M
" ' " ' ‘ ’ ‘ " G—elslando
{4,1,1 . B—Ellsland1
. I C , 9 _4;\n . +———+ Island2
' I . . ‘1 V V . (is—Alslanda
an 33 E . 3 35—05—- Island4
~ A . ,, .3? - g ' “Eass......_g,_....._§ law '0‘:- Island 5
. .- - a... .. .. Island 8
E‘ " " “7‘ ’W‘ '“ " ‘ x——X Island?
» Islands
I* . H Islands
...._.._. Island 10
r.—. Island 11
0.0 30000.0 I000.0 90000.0 1200000
Tlme (seconds)

Figure 4.28. iiGA Fitness re-evaluated for each island as
a function of time. All re-evaluations performed with 3-D
FEM containing 12,272 degrees of freedom for an annular

multi-material composite flywheel.

Khan.” mnr—w
. l

0.0

107

 

 

 

 

L A A l

A

 

‘ G—Olsland 0
‘ B—El lsland1
. 'tH'Isand 2

H Island 3
Island 4

‘ 1 Island 5

Island 6

. x-—>< Island 7

~ Island 8
H Island 9
Island 1!

‘ -—-—1 Island 11

 

0.0 30000.0 “000.0 90000 .0 120000 .0

Time (seconds)

Figure 4.29. Hybrid iiGA fitness re-evaluated for each island
as a function of time. All re-evaluations performed with 3-D
FEM containing 12,272 degrees of freedom for an annular
multi-material composite flywheel.

 

 

 

108

 

 

 

G—O Island 0

 

‘ B—El Island 1
. >———s Island 2
_ Hlslanda
Island 4

 

 

 

A A A 4L

0.0

 

, ‘ i:~»-~-~.'= Island 5
Island 6
. x———x Island 7
. Island 8
H Island s
- - Island 10
‘ .__. Island 11

 

 

6000 .0 9000 .0 120000
Tlrne (seconds)

0.0 3000.0

15000.0

Figure 4.30. iiGA quickly finding building blocks at a low
level of resolution to inject into islands of higher resolution

for an annular multi-material composite flywheel.

Jr— ‘

 

—._”. :1. mi.-- ‘3‘ ‘5:

 

109

 

 

 

 

 

 

 

 

 

 

 

15.0 , , 4
." I
w. . {I
I”... G—O Island 0
_‘ . ; - - --"7:"‘f~-~ , 1Hlyand1
g! g 31' ‘ g” .l; \ 0 - Island2
‘ r. _‘ '- ~ ... -=- ' Hlslanda
10.0 " .1 ‘ ' I "
g A of "’ Island 4
g 3 . .j Ida—Mr 1'", - -.._..- Tl? “Quhwg’n ”f :“"“"'t Island 5
u. I: 91/4" ‘ Islands
5;; ' . X——-X Island 7
4.- Islands
5.0 s d H Island 9
Island 10
* .—. Island 11
I .
00 . l . l . l A l 4
0.0 3000.0 8000.0 9000.0 12000.0 15100.0

Time (seconds)

Figure 4.31. Hybrid iiGA quickly finding building blocks at a
low level of resolution to inject into islands of higher resolution
for an annular multi-material composite flywheel.

Fitness

150

5.0

0.0

Figure 4.32. Twenty island “ring” GA topology.

110

 

 

r

   

 

 

KEY: Evaluation
Tool

  

3-D FEM
(12,272 DOF)

 

 

 

 

4 L A

9—0 Island 0
H Island 1
w Island 2
H Island 3
Island 4
—~i Island 5
Island 6
x———x Island 7
. Island 8
H Island 9
WW. Island 1U
:— Island 11

4

L

 

l

 

0.0

1000000

20 0000.0

Time (seconds)

300000 .0

4000000

5000000

Figure 4.33. “Ring” topology GA results for an annular multi-material

composite flywheel.

Fitness
Definition

 

 

15.0

5.0

0.0

111

 

 

—— Hybrid iiGA *
——- iiGA «
.......... Ring Topology GA

 

I' «

 

 

1 4L

0.0 1 00000.0

A

2000000 300000 .0 4000000 5000000
Time (seconds)

 

Figure 4.34. Comparison of iiGA to hybrid iiGA to “ring” topology GA.
For this problem, both the iiGA and hybrid iiGA search extremely
efficiently when compared to a “ring” topology GA. The hybrid iiGA
discovers solutions with higher fitness values when compared to the
“ring” topology GA and the iiGA alone.

112

Best Flywheel
Found by “Ring” —>
Topology GA

Best. Flywheel

Found by iiGA —>

 

Best Flywheel
Found by Hybrid —>
iiGA

 

 

Material

I E-Glass Epoxy [00]

I E-Glass Epoxy [90°]

 

 

\flflflfllcfllnflflalnnlnllnl

 

IiIIIOOIOODOOODIIIDOIIIO

D-l

IIIIOODOIIIIOOOIOOOOOOiI

Figure 4.35. Best flywheel found by “ring” topology GA, iiGA
and hybrid iiGA for an annular multi-material composite flywheel.

 

 

113

KEY: Evaluation Fitness
r__.,.__1.__1 Tool Definition

Three Rings I O I 3-D FEM SED
(low resolution) : l I
l

(54, 190, 260 DOF)

0 3-D FEM SED
(340 DOF)

 
     

 

Six Rings |
(medium resolution)

      

Twelve Rings
(high resolution)

3-D FEM SED
(962 DOF)

, 3-D FEM SED
Twenty-Four ngsl (2,774 DOF)

(higher resolution) I

L _..L _..L__J .3-DFEM SED
] “T \ (12,272 DOF)

Hybrid Local Search Plain GA Hybrid Thershold Accepting Algorithm

 

Figure 4.36. Hybrid island injection GA toplogy used in the
constrained optimization of a multi-material annular flywheel.

114

 

eoo . . . . . .'.

 

60.0

 

 

 

 

G—e Island 0
G—a Island 1

4 .. ..’Island2

 

A——-A Island 3
Island 4

' » ..-¢--‘——'§;,- Idand 5

Island 6

‘ >6—x Island 7

- .. Island 8
H Island 9
20.0 - I W... Island 10
-—- lsland 11
0.0 ‘ ‘ * ‘ ‘ ’ m
0.0 30000.0 600000 90000 .0 1200000
Time (seconds)

Figure 4.37. Hybrid iiGA fitness re-evaluated for each island
as a function of time. All re-evaluations performed with 3-D
FEM containing 12,272 degrees of freedom for the constrained
optimization of an annular composite multi-material flywheel.

115

80.0 - u . u

 

 

60.0 - A ....,_.”" :_ .. p . .
r :: ‘l - ' ' - G—Olslandu
I /"' B—Ellsland1
" I] f ‘ ‘-‘---‘—L> Island 2
. A—fl Island 3
3...... ~--‘ "51.---.. ”w Idand 4
~ ' 7»*'-*-—‘:5" Island 5
Island 6
. H Island 7
-~~-- Island 8
B—EI Island 9
20.0 b - ----~' Island10
-—' Island 11

 

   

 

 

0.0 ‘ ' ‘ ‘ ‘
0.0 10000.0 20000.0 30000.0

Time (seconds)

 

Figure 4.38. Hybrid iiGA quickly finding buildin blocks at a
low level of resolution to inject into islands of hig er resolution
for the constrained optimization of an annular multi-material
composite flywheel.

116

 

 

 

Best GEY: Evaluation Fitness‘

ver Tool Definition

3-D FEM SED
(54, 190, 260 DOF)
0 3-D FEM SED

(342 DOF)

est 3-D FEM SED

(962 DOF)
.ver

 

@ 3—D FEM SED
(2,774 DOF)

. 3-D FEM SED
| (12,272 DOF) ;
f Material ‘

I E-Glass Epoxy [00]

 

Best
Ever

E

 

E-Glass Epoxy [90°]

Aluminum

 

 

 

J

 

 

 

»1\ucaoocooo

 

 

 

Figure 4.39. Best flywheel found at each level of resolution
and fitness definition for the constained optimization of a
multi-material composite flywheel.

Chapter 5

Summary and Conclusions

This thesis presents an approach to the optimization of composite flywheels and
was broken down into four major segments. To familiarize the reader, a broad range of op-
timization topics were first reviewed. These topics included past research involving a
broad class in the optimization of engineering systems. The second segment presented the
mechanics of GAs. The third segment contained the development and accuracy confirma-
tion of the finite element models used to model flywheels. The fourth segment contained
the optimization of flywheels using GAS.

The simple GA found the optimal shape of a flywheel with stress free edges that has
a 55% increase in SED when compared to a constant stress flywheel. The commonly used
rotor shape of the constant stress flywheel is not optimal when a flywheel has stress free
edges.

The iiGA can approach optimization problems in a unique way. For many prob-
lems, the iiGA can be used to break down a complex fitness function into “sub-fitness”
functions, which represent “good” aspects of the overall fitness. Also, the iiGA can use dif-
fering evaluation tools. A simplified analysis tool can be used to quickly search for good
building blocks. This, in combination with searching at various levels of resolution, makes
the iiGA efficient and robust. Mimicking a smart engineer, the iiGA can first quickly eval-

uate the overall response of a structure with a coarse representation of the design and finish

117

118

the job off by slowly increasing the levels of refinement until a “finely tuned” structure has
been evolved. This approach allows the iiGA to decrease computational time and increase
the robustness of a typical GA. The iiGA seems to have a bright future in the optimization
of large scale three-dimensional composite structures.

The hybrid iiGA and iiGA were found to search more efficiently than the “ring” to-
pology GA for the case of the enumerable search space of a solid isotropic flywheel. In
fact, the iiGA found the global solution in all five independent runs while the “ring” topol-
ogy GA never found the global solution for any run. The hybrid iiGA that used the com-
bination of a TA algorithm and a local search method found the global optimum in less time
than all other optimization techniques, searching less than 5% of the total search space. The
hybrid iiGA found the global solution on average in less than 0.5% of the time used in the
enumeration of the search space, by measuring more than half of the fitness evaluations
with the plane stress finite element model.

The hybrid iiGA and iiGA were also found to search more efficiently when com-
pared to the “ring” topology GA for the case of the annular multi-material composite fly-
wheel. Both the hybrid iiGA and iiGA converged five times faster when compared to the
“ring” topology GA. The hybrid iiGA converged to a slightly higher fitness when com-
pared to both the iiGA and “ring” topology GA.

A hybrid iiGA designed a “safer” composite flywheel by penalizing flywheels with
large “air gap” growths. The angular velocity was constrained to a lower “target” angular
velocity while maximizing SED and penalizing flywheels that have a large “air gap”
growth. The angular velocity in the best solution found at the highest level of geometric

resolution was extremely close to the “target” angular velocity. The “air gap” growth in

119

the best solution found at the highest level of geometric resolution in the constrained opti-
mization problem was cut in half when compared to the unconstrained optimization ap-
proach. The SED for the best solution found at the highest level of geometric resolution in
the constrained optimization problem was 21.6% lower than the unconstrained single cri-

terion optimization approach.

120

REFERENCES

1.

10.

11.

12.

13.

Belegundu A. D., and Zhang S., 1992, “Robustness of Design Through Minimum
Sensitivity” Journal of Mechanical Design, Vol 114, pp213—217.

Bendsoe M. P., Olholf N., Taylor J. E., 1984, “A Variational Formulation of Mul-
ticriterion Structural Optimization”, Journal of Structural Mechanics, Vol. 11, No.
4, pp. 523-544.

Bhavikatti K., Ramakrishnan C., 1980, “Optimum Shape Design of Rotating
Disks”, Computers and Structures, Vol. 11, pp. 397-401.

Carlson S.E., Shonkwiler R., and Ingrim M.E., 1994, “Comparison of Three Non-
derivative Optimization Methods with a Genetic Algorithm For Component Selec-
tion”, Journal of Engineering Design, Vol. 5, No. 4, pp 367-378.

Chen G., Bruno R., Salama M., 1991, “Optimal Placement of Active/Passive Mem-
bers in Truss Structures Using Simulated Annealing”, AIAA Journal, Vol. 29, No.
8, pp. 1327-1334.

Christensen R., Wu E., 1977, “Optimal Design of Anisotropic (Fiber-Reinforced)
Flywheels”, Journal of Composite Materials, Vol. 11, pp395-404.

Chu D. N., Xie Y. M., Hira A., Steven G. P., “Evolutionary Structural Optimization
for Problems with Stiffness Constraints,” Finite Elements in Analysis and Design,
Vol. 22, 1996, pp. 239-251.

Crosseley W. A., and Lannanen D. H., 1996, “Conceptual Design of Helicopters via
Genetic Algorithms”, Journal of Aircraft, Vol. 33, No. 6, pp. 1062-1069.

Eby D., 1996, “Optimization of Anisotropic Flywheels Via Genetic Algorithms”,
Senior Research Project, Michigan State University, E. Lansing, MI, 48823.

Fabbri G., 1997, “A Genetic Algorithm for Fin Profile Optimization”, International
Journal of Heat and Mass Transfer, Vol. 40, No. 9, pp2l65-2172.

Foster N., and Dulikravich G., 1997, “Three-Dimensional Aerodynamic Shape op-
timization Using Genetic And Gradient Search Algorithms”, Journal of Spacecraft
and Rockets, Vol. 34, No. 1, pp. 36-42.

Flynn R., and Sherman P., 1995, “Multicriterion Optimization of Aircraft Panels:
Determining Viable Genetic Algorithm Configurations”, International Journal of
Journal of Intelligent Systems, Vol. 10, pp. 987-999.

Furura H., and Haftka R. T., 1993, “Locating Actuators for Vibration Suppression
on Space Trusses by Genetic Algorithms”, Structures and Controls Optimization,

14.

15.

16.

17.

l8.

19.

20.

21.

22.

23.

24.

25.

26.

121

ASME, AD Col. 38, pp.1-11.

Furuya H., and Haftka R., 1995, “Placing Actuators in Space Structures by Genetic
Algorithms and Effectiveness Indicies”, Structural Optimization, Vol. 9, pp. 69-75.

Genta G., 1989, “Some Considerations on the Constant Stress Disc Profile”, Mec-
canica, Vol. 24, pp. 235-248.

Genta G., Domenico D., 1995, “Use of Genetic Algorithms for the design of R0-
tors”, Meccanica, Vol. 30, pp. 707-717.

Georgian J ., 1989, “Optimum Design of Variable Composite Flywheel”, Journal of
Composite Materials, Vol. 23, pp.2-10.

Goldberg D., Genetic Algorithms in Search and Machine Learning, 1989, Addison-
Wesley Pub. Co. Inc.

Goodman E. D., 1994, GALOPPS, The Genetic Algorithm Optimized for Portabil-
ity and Parallelism System, Tech. Rept. #94-5, MSU GARAGe, Michigan State

University, 100pp.

Goodman E. D., Averill R. C., Punch III W. F., and Eby D. J., 1997, “Parallel Ge-
netic Algorithms in the Optimization of Composite Structures”, Second World Con-
ference on Soft Computing (WSCZ), www.egr.msu.edu/~ebydavid/ver5.html.

Grandi R., 1993, “Structural Optimization with Frequency Constraints-A Review”,
AIAA Journal, Vol. 31, No. 12, pp. 2296-2301.

Haslinger J .,and Jedelsky D., 1996, “Genetic Algorithms and Fictitious Based Ap-
proaches in Shape Optimization”, Structural Optimization, Vol. 12, pp 257-264.

Homaifar A., Lai H. Y., and McCormick E., 1994, “System Optimization of Turbo-
fan Engines Using Genetic Algorithms”, Applications in Mathematical Modeling,
Vol. 18, pp. 72-83.

Horn J., 1997, “Overview: Evolutionary computation for Multipleobjective Opti-
mization”, Proceedings of the Seventh International Conference on Genetic Algo-
rithms, Michigan State University.

Jenkins W. M., 1991, “Towards Structural Optimization Via the Genetic Algo-
rithm”, Computers & Structures, Vol. 40, No. 5, pp 1321-1327.

Khot N., and Veley D., 1991, “Use of Robustness Constraints in the Optimum De-

sign of Space Structures”, Journal of Intelligent Material System and Structures,
Vol. 2. PP. 161-176.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

37.

38.

122

Kogiso N., Watson L. T., Gurdal Z., and Haftka R. T., 1994, “Genetic Algorithms
with Local Improvement for Composite Laminate Design”, Structural Optimiza-
tion, Vol. 7, pp. 207-218.

Mallot B., Averill R., Goodman E., Ding Y., and Punch 111 W., 1996, “Use of Ge-
netic Algorithms for Optimal Design of Laminated Composite Cantilever Sand-
wich Panels with Bending-Twisting Coupling”, presented at AIAA SDM
(Structures, Dynamics and Materials).

Mares C., and Surace C., 1996, “An application of Genetic Algorithms to Identify
Damage in Elastic Structures”, Journal of Sound and Vibration, Vol. 195, No. 2,
pp. 195—215.

MARC Users Manual, Volume B, Part 1, MARC Analysis Research Corporation.

Miller B., and Goldberg D., 1995, “Genetic Algorithms, Tournament Selection, and
the Effect of Noise”, IlliGAL Report No. 95006, University of Illinois at Urbana-

Champaign, Department of General Engineering.

Nakanishi Y., and Nakagiri S., 1996, “Optimization of Frame Topology using
Boundary Cycle and Genetic Algorithm”, Japanese Society of Mechanical Engi-
neering International Journal, Series A, Vol. 39, No. 2, pp. 279-285.

Obayahi S., and Takanashi S., 1996, “Genetic Optimization of Target Pressure Dis-
tributions for Inverse Design Methods”, AIAA Journal, Vol. 34, No. 5, pp. 881-886.

Olhoff N., and Taylor J. E., 1983,“On Structural Optimization”, Journal of Applied
Mechanics, Vol 50, pp. 1139-1151.

Osyczka A., Kuchta W., and Czula R., 1994, “Computer Aided Multicriterion Sys-
tem for Computationally Expensive Functions”, Structural Optimization, Vol. 8,
pp. 37-41.

Parmee I., and Vekeria H., 1997, “Co-Operative Evolutionary Strategies for Single
Component Design”, Proceedings of the Seventh International Conference on Ge-
netic Algorithms, Michigan State University.

Punch 11] W. F., Averill R. C., Goodman E. D., Lin S. C., and Ding Y., February
1995,“Design Using Genetic Algorithms - Some Results for Laminated Composite
Structures”, IEEE Expert, Vol. 10(1), pp. 42-49.

Punch III W. F., Averill R.C., Goodman E.D., Lin S. C., Ding Y., and Yip Y.C.,
1994, “Optimal Design of Laminated Composite Structures using Coarse-Grain

Parallel Genetic Algorithms,” Computing Systems in Engineering, Vol. 5, No. 4-6,
pp. 414-423.

 

39.

40.

41.

42.

43.

45.

46.

47.

48.

49.

50.

51.

123

Punch III W., Averill R., Goodman B, Lin 5., and Yip Y., 1994, “Optimal Design
of Laminated Composite Structures using Coarse-Grain Parallel Genetic Algo-
rithms”, Computing Systems in Engineering, Vol. 5, pp. 415-423.

Quagliarella D., and Cioppa A., 1994, “Genetic Algorithms applied to the Aerody—
namic Design of Transonic Airfoils”, Journal of Aircraft, Vol. 32, No. 4, pp. 889-
891.

Queipo N., Devarakonda R., and Humphrey J., 1994, “Genetic Algorithms for the
Thermosciences Research: Application to the Optimized Cooling of Electronic

Components”, International Journal of Heat and Mass Transfer, Vol 37, No. 6, pp.
893-908.

Rajan S. D., 1995, “Sizing, Shape and Topology Design Optimization of Trusses
using Genetic Algorithm”, Journal of Structural Engineering, Vol. 121, No. 10, pp.
1480-1487.

Richef R., and Haftka R., 1993, “Optimization of Laminate Stacking Sequence for
Buckling Load Maximization by Genetic Algorithm”, AIAA Journal, Vol. 31, No.
5. pp. 951-955.

Ries D., Kirk J., 1992, “Design and Manufacturing for a Composite Multi-Ring Fly-

wheel”, Proceedings of the 27th Intersociety Energy Conversion Engineering Con-
ference, Vol. 2, pp 4.43-4.48.

Ruthenbar R., 1989, “Simulated Annealing Algorithms: An Overview”, IEEE C ir-
cuits and Devices Magazine, Vol. 5, No. 1, pp. 19-26.

Sangren E., Ragsdell K. M., 1983, “Optimal Flywheel Design with a General thick-
ness Form Representation”, Journal of Mechanisms, Transmissions, and Automa-
tion in Design, Vol. 105, pp. 425-433.

Schraudolph N., and Grefenstette J ., 1992, A User's Guide to GAucsd 1.4, July.

Seireg A., Surana K., 1970, “Optimum Design of Rotating Disks”, Journal of En-
gineering for Industry, Vol. 92, pp. 1-10.

Stadler W., and Dauer J ., 1992, “Multicriterion Optimization in Engineering: A Tu-
torial and Survey”, Structural Optimization: Status and promise, Progress in Astro-
nautics and Aeronautics, Vol. 150, AIAA Inc., pp 209-249.

Sundaresan S., Ishii K., and Houser D., 1992, “Design Optimization For Robust-
ness Using Performance Simulation Programs”, Engineering Optimization, Vol.

20, pp.163-178.

Todoroki A., and Watanabe K., 1995, “Applications of Genetic Algorithms to Stiff-

 

52.

124

ness Optimization of Laminated Composite Plates with Stress-Concentrated Open
Holes, Japanese Society of Mechanical Engineering International Journal, Series
A, Vol. 38, No. 4, pp. 458-404.

Ugural A., Fenster S., 1995, Advanced Strength and Applied Elasticity, pp. 339-
344.

"lllllllilllllllli