ABSTRACT
Title of dissertation: ROBUST MULTI-OBJECTIVE OPTIMIZATION
OF HYPERSONIC VEHICLES UNDER
ASYMMETRIC ROUGHNESS-INDUCED
BOUNDARY-LAYER TRANSITION
Kevin Michael Ryan,
Doctor of Philosophy, 2014
Dissertation directed by: Professor Mark J. Lewis
Department of Aerospace Engineering
The effects of aerodynamic asymmetries on hypersonic vehicle controllability
and performance were investigated for a wide range of geometries. Asymmetric
conditions were introduced by an isolated surface roughness that forces boundary-
layer transition resulting in a turbulence wedge downstream of the disturbance. The
disturbance simulates the effects of physical deformations that may exist on a vehicle
surface or leading edge, such as protruding edges of thermal protection system tiles
or non-uniform surface roughness.
Both multi-objective and robust multi-objective optimization studies were per-
formed. Traditional multi-objective optimization methods were used to identify
vehicle designs that are best suited to withstand spanwise asymmetric boundary-
layer transition while retaining its performance and payload requirements. Trade-
offs between vehicle controllability and performance were analyzed. A novel multi-
objective based robust optimization method to solve single-objective optimization
problems with environmental parameter uncertainty was proposed and tested. Un-
like commonly used robust optimization methods, the multi-objective method for-
mulates an optimization problem such that post-optimality data handling techniques
can identify multiple robust designs from a single solution set. This allows for com-
parisons to be made between different types of robust designs, thus providing more
information about the design space. Comparisons were made between the robust
multi-objective optimization formulation and conventional robust regularization-
and aggregation-based methods. The results, performance, and philosophies of each
method are discussed.
Design trends were identified for classifying the optimum and robust opti-
mum designs of hypersonic vehicle shapes under boundary-layer transition uncer-
tainties. Traditional multi-objective optimization results show that two types of
vehicle shapes bound the set of Pareto-optimal solutions: wedge-like and cone-like.
The L2-norm optimum design, representing a compromise between the competing
shapes, was a hybrid wedge-cone shape. The robust optimization results show that
a flat wedge-like vehicle design is best for a worst-case scenario, while a pyramidal
shaped vehicle design minimizes the expected detrimental effects on vehicle control-
lability. The analyses prove that the novel robust optimization method can provide
a range of robust optimum results, while also capturing trade-offs within the design
space, providing capabilities not available in state-of-the-art robust optimization
methods.
ROBUST MULTI-OBJECTIVE OPTIMIZATION OF
HYPERSONIC VEHICLES UNDER ASYMMETRIC
ROUGHNESS-INDUCED BOUNDARY-LAYER TRANSITION
by
Kevin Michael Ryan
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2014
Advisory Committee:
Professor Mark J. Lewis, Co-Chairman/Advisor
Associate Professor Kenneth H. Yu, Co-Chairman/Advisor
Associate Professor Christopher P. Cadou
Professor Roberto Celi
Professor Ashwani K. Gupta
c? Copyright by
Kevin Michael Ryan
2014
Dedication
To my girlfriend, Suzanne. Thank you for all your love and support.
ii
Acknowledgments
I would like to thank my advisors Dr. Mark Lewis and Dr. Kenneth Yu for
giving me the opportunity to conduct this research. Dr. Lewis, your knowledge
and passion for the field of aerospace has inspired me throughout the completion
of this dissertation and continues to do so. Dr. Yu, your support, patience, and
advice enabled me to push forward through meetings, conferences, and on to the
completion of my Ph.D. I would also like to thank my committee members Dr.
Cadou, Dr. Celi, and Dr. Gupta whose insight and help guided me through my
research. Their time and effort is much appreciated.
I would like to acknowledge the Air Force Office of Scientific Research funding
that supported this work through the UMD/AEDC Hypersonic Workforce Revital-
ization Program. Their excitement and encouraging words at our technical reviews
would always revitalize my motivation towards my research.
I would also like to thank my family for their constant support through my
nine year academic adventure at the University of Maryland. Finally, I would like
to thank my fellow co-workers in the Hypersonics Research Group and Advanced
Propulsion Research Lab whose help and support I could not have done without. A
special thanks goes out to Neal Smith, Vijay Ramasubramanian, and Jeremy Knittel
for their insightful discussions about my work, help with editing my papers, and their
intellectual interrogations to help me prepare for my Comprehensive Exam.
iii
Table of Contents
List of Tables vii
List of Figures viii
List of Symbols and Abbreviations xi
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Objectives and Contributions . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . 7
1.4 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.1 The Design, Aerodynamics, and Optimization of Hypersonic
Vehicles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2 Roughness-Induced Boundary-Layer Transition . . . . . . . . 15
1.4.3 Robust Optimization . . . . . . . . . . . . . . . . . . . . . . . 17
1.4.3.1 Formulations . . . . . . . . . . . . . . . . . . . . . . 18
2 Hypersonic Aerodynamics 23
2.1 Vehicle Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2 Infinitely Sharp Leading Edge . . . . . . . . . . . . . . . . . . . . . . 27
2.2.1 Pressure Forces . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2.2 Viscous Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.3 Moment Calculations . . . . . . . . . . . . . . . . . . . . . . . 29
2.3 Blunted Leading Edge . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1 Pressure Forces . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.2 Viscous Forces and Leading Edge Heat Flux . . . . . . . . . . 33
2.3.3 Moment Calculations . . . . . . . . . . . . . . . . . . . . . . . 35
3 Asymmetric Boundary-Layer Transition 36
3.1 Roughness-Induced Transition . . . . . . . . . . . . . . . . . . . . . . 36
3.1.1 Physical Effects . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.2 Parametric Trends . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.2.1 Roughness Height and Effective Roughness . . . . . 38
iv
3.1.2.2 Roughness Shape . . . . . . . . . . . . . . . . . . . . 40
3.1.2.3 Turbulence Wedge Spreading Angle . . . . . . . . . . 41
3.1.3 Turbulence Wedge Model . . . . . . . . . . . . . . . . . . . . 43
4 Optimization Methods, Algorithms, and Definitions 45
4.1 Traditional Optimization Methods . . . . . . . . . . . . . . . . . . . . 46
4.1.1 Genetic Algorithms . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.2 Multi-Objective Genetic Algorithms (MOGA) . . . . . . . . . 50
4.2 Robust Optimization Methods . . . . . . . . . . . . . . . . . . . . . . 54
4.2.1 Robust Regularization . . . . . . . . . . . . . . . . . . . . . . 56
4.2.2 Aggregation Approach . . . . . . . . . . . . . . . . . . . . . . 57
4.2.3 Novel Multi-Objective Approach . . . . . . . . . . . . . . . . . 58
4.2.3.1 Robust MOOP Formulation . . . . . . . . . . . . . . 60
4.2.3.2 Convergence Criteria and Data Handling . . . . . . . 62
4.3 Illustrative Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.3.1 Original Optimization Problem . . . . . . . . . . . . . . . . . 67
4.3.2 Analytical Robust Optimum Solutions . . . . . . . . . . . . . 70
4.3.2.1 Optimal Expected Value . . . . . . . . . . . . . . . . 70
4.3.2.2 Optimal Worst-Case Scenario . . . . . . . . . . . . . 71
4.3.3 Robust MOGA Formulation & Data Handling Solutions: Base-
line Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.4 Accuracy and Convergence Tests . . . . . . . . . . . . . . . . . . . . 81
4.4.1 Effects of Changes in Initial Population Size . . . . . . . . . . 84
4.4.2 Effects of Changes in Number of Histogram Bins . . . . . . . . 86
4.4.3 Effects of Changes in Number of Uncertainty Subsets . . . . . 87
5 Vehicle Design Optimization Methodology 93
5.1 Optimization Design Point, Constraints, and Assumptions . . . . . . 94
5.2 Baseline Multi-Objective Optimization . . . . . . . . . . . . . . . . . 97
5.3 Influence of Asymmetric Boundary-Layer . . . . . . . . . . . . . . . . 101
5.3.1 Multi-Objective Optimization . . . . . . . . . . . . . . . . . . 101
5.3.2 Robust Optimization Formulations . . . . . . . . . . . . . . . 105
5.4 MOGA Method and Parameters . . . . . . . . . . . . . . . . . . . . . 109
5.4.1 Multi-Objective Studies MOGA Parameters . . . . . . . . . . 109
5.4.2 Robust Optimization Studies MOGA Parameters . . . . . . . 112
6 Results 114
6.1 Baseline Multi-Objective Optimal Geometries . . . . . . . . . . . . . 116
6.1.1 Infinitely Sharp Case . . . . . . . . . . . . . . . . . . . . . . . 116
6.1.2 Blunted Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.2 Optimal Geometries under Asymmetric Boundary-Layer Conditions . 124
6.2.1 Multi-Objective Results . . . . . . . . . . . . . . . . . . . . . 124
6.2.1.1 Infinitely Sharp Case . . . . . . . . . . . . . . . . . . 124
6.2.1.2 Blunted Case . . . . . . . . . . . . . . . . . . . . . . 133
6.2.2 Robust Optimization Results . . . . . . . . . . . . . . . . . . 136
v
6.2.2.1 Robust Regularization Results . . . . . . . . . . . . . 136
6.2.2.2 Explicit Averaging Results . . . . . . . . . . . . . . . 137
6.2.2.3 Perturbation Results . . . . . . . . . . . . . . . . . . 138
6.2.2.4 Robust MOGA Results . . . . . . . . . . . . . . . . 141
6.2.2.5 Result Comparisons and Vehicle Analysis . . . . . . 146
6.2.2.6 Design Space Analysis . . . . . . . . . . . . . . . . . 153
7 Conclusions 157
Bibliography 164
vi
List of Tables
4.1 MOGA options used for solving Equation 4.30. . . . . . . . . . . . . . 75
4.2 Worst-case scenario comparison settings in Figure 4.7. . . . . . . . . . 81
4.3 Accuracy and convergence input settings test matrix. . . . . . . . . . 83
5.1 Design point specifications. . . . . . . . . . . . . . . . . . . . . . . . . 95
5.2 MOGA parameters for the multi-objective optimization studies. . . . 111
5.3 Genetic algorithm options for the robust optimization studies. . . . . 112
6.1 Pareto-optimal point objective function values for the infinitely sharp
leading baseline edge case. . . . . . . . . . . . . . . . . . . . . . . . . 120
6.2 Pareto-optimal point design vector values for the infinitely sharp lead-
ing edge baseline case. . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.3 Pareto-optimal point objective function values for the blunted leading
edge baseline case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.4 Pareto-optimal point design vector values for the blunted leading edge
baseline case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.5 Pareto-optimal point objective function values for the ABLT infinitely
sharp leading edge case. . . . . . . . . . . . . . . . . . . . . . . . . . 132
6.6 Pareto-optimal point design vector values for the ABLT infinitely
sharp leading edge case. . . . . . . . . . . . . . . . . . . . . . . . . . 132
6.7 Pareto-optimal point objective function values for the ABLT blunted
leading edge case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
6.8 Pareto-optimal point design vector values for the ABLT blunted lead-
ing edge case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
6.9 Robust optimum designs. . . . . . . . . . . . . . . . . . . . . . . . . . 147
vii
List of Figures
2.1 Example convex parametric waverider shape. . . . . . . . . . . . . . . 25
2.2 Example concave parametric waverider shape. . . . . . . . . . . . . . 25
2.3 Example blunted concave parametric waverider shape. . . . . . . . . 26
2.4 Profile cross-section view at some given local span location. . . . . . . 28
2.5 Profile cross-section view at some given local span location. . . . . . . 32
2.6 Shear stress components in the x and z directions. . . . . . . . . . . . 35
3.1 Example turbulence wedge generated from isolated surface roughness. 39
3.2 Effect of roughness shape on transition. . . . . . . . . . . . . . . . . . 41
3.3 Variation of turbulent spreading with boundary-layer edge Mach num-
ber. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Example of an isolated roughness element applied to upper surface of
a generic hypersonic vehicle. . . . . . . . . . . . . . . . . . . . . . . . 44
4.1 Generalized genetic algorithm process. . . . . . . . . . . . . . . . . . 49
4.2 fi.e. curves as a function of p and ~x. . . . . . . . . . . . . . . . . . . . 68
4.3 fi.e. curves as a function of p and ~x. . . . . . . . . . . . . . . . . . . . 69
4.4 Mean ?PM generation history. . . . . . . . . . . . . . . . . . . . . . 76
4.5 Baseline analytical and predicted optimum expected value design vec-
tors for p values of 0.2, 0.3, 0.4, and 0.5, respectively. . . . . . . . . . 77
4.6 Baseline analytical and predicted optimum expected value design vec-
tors for p values of 0.6, 0.7, and 0.8, respectively. . . . . . . . . . . . . 78
4.7 Analytical and predicted optimum worst-case scenario design vectors
for various [a, b] values. . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.8 Speed, quality, and accuracy changes with change in initial population
size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.9 Speed, quality, and accuracy changes with change in number of his-
togram bins used in the data handling procedure. . . . . . . . . . . . 91
4.10 Speed, quality, and accuracy changes with change in number of uncer-
tainty subsets in the multi-objective optimization problem formulation. 92
5.1 Parametric locations of leading edge isolated roughness (planform view).102
viii
5.2 Example locations of leading edge isolated roughnesses as function of
p (planform view). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.1 Two objective function example showing extreme Pareto-optimal points
and L2-norm optimum Pareto-optimal point. . . . . . . . . . . . . . . 116
6.2 Extreme Pareto-optimal point vehicle designs with best f1, f2, f3
values for the infinitely sharp leading edge baseline case (isometric
views). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
6.3 L2-norm Pareto-optimal point vehicle design for the infinitely sharp
leading edge baseline case. . . . . . . . . . . . . . . . . . . . . . . . . 119
6.4 Extreme Pareto-optimal point vehicle designs with best f1, f2, f3, f4
values for the blunted leading edge baseline case (isometric views). . . 122
6.5 L2-norm Pareto-optimal point vehicle design for the blunted leading
edge baseline case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.6 Extreme Pareto-optimal point vehicle designs with best f1, f2, f3,
values and the design with the worst f3 value for the ABLT infinitely
sharp leading edge case (isometric views). . . . . . . . . . . . . . . . . 125
6.7 L2-norm Pareto-optimal point vehicle design for the ABLT infinitely
sharp leading edge case. . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.8 Isometric views of the best f3 Pareto-optimal design with a spanwise
disturbance at p = 1, 2, 3, and 4 for the ABLT infinitely sharp case. . 127
6.9 Isometric views of the worst f3 Pareto-optimal design with a spanwise
disturbance at p = 1, 2, 3, and 4 for the ABLT infinitely sharp case. . 127
6.10 Isometric views of L2-norm Pareto-optimal design with a spanwise
disturbance at p = 1, 2, 3, and 4 for the ABLT infinitely sharp case. . 127
6.11 Spanwise distribution of F3,scaled for best, worst, and L2-norm vehicle
designs for the ABLT infinitely sharp leading edge case. . . . . . . . . 128
6.12 Spanwise distribution of the y and z components of F3 for best, worst,
and L2-norm vehicle designs for the ABLT infinitely sharp leading
edge case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
6.13 Extreme Pareto-optimal point vehicle designs with best f1, f2, f3, f4
values and the design with the worst f3 value for the ABLT blunted
leading edge case (isometric views). . . . . . . . . . . . . . . . . . . . 134
6.14 L2-norm Pareto-optimal point vehicle design for the ABLT blunted
leading edge case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
6.15 Robust regularization-based method generation history. . . . . . . . . 137
6.16 Explicit averaging generation history. . . . . . . . . . . . . . . . . . . 138
6.17 Generation history of the mean ~x values over the population for the
perturbation method. . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.18 Robust MOGA mean ?PM generation history. . . . . . . . . . . . . 141
6.19 Average FIV2 vs. average FIV1 over all generations for each of the 96
Pareto-optimal solutions. . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.20 Four-dimensional Pareto-optimal front. . . . . . . . . . . . . . . . . . 143
6.21 Three-dimensional projections of the Pareto-optimal front. . . . . . . 145
6.22 Optimal vehicle design for the robust regularization-based method. . 148
ix
6.23 Optimal vehicle design for the robust MOGA method with robust
regularization equivalent data handling. . . . . . . . . . . . . . . . . . 148
6.24 Optimal vehicle design for the explicit averaging method. . . . . . . . 149
6.25 Optimal vehicle design for the perturbation method. . . . . . . . . . . 149
6.26 Optimal vehicle design for the robust MOGA method with explicit
averaging and perturbation equivalent data handling. . . . . . . . . . 150
6.27 Original objective function values of each optimum design vector as
a function of environmental parameter p. . . . . . . . . . . . . . . . . 151
6.28 Two-dimensional projections of Pareto-optimal front. . . . . . . . . . 155
6.29 Vehicle designs A through F from Figure 6.28(b). . . . . . . . . . . . 156
x
List of Symbols and Abbreviations
A planform power law scaling parameter
a uncertainty neighborhood lower bound
B base upper surface power law scaling parameter
b uncertainty neighborhood upper bound
C base lower surface power law scaling parameter
CM moment coefficient
Ch Chapman-Rubesin parameter
Drag total vehicle drag
d upper and lower surface separation distance
dA panel area
E conditional expected value
F robust counterpart objective function
f (original) objective function
g constraint function
gw wall-to-inviscid stagnation temperature ratio
H stagnation enthalpy parameter (transformed plane)
h all generations from 1 to t
J number of constraint functions
K penalty factor
L vehicle length
Lift total vehicle lift
M Mach number, or number of objective functions
~M moment vector
Mx moment about the x-axis
My moment about the y-axis
Mz moment about the z-axis
m base & lower surface power law exponent
N number of panels in vehicle mesh
N normal distribution
n planform surface & base surface power law exponent, or number
of samples
n? outward unit normal
P pressure
P bounded parameter uncertainty set
Pr probability
p environmental parameter
p? random parameter value
p? discrete set of p values
Q number of objective functions
q heat flux
R gas constant
Re Reynolds number
~r position vector from CG to panel
xi
ry distance from CG to panel in y direction
rz distance from CG to panel in z direction
Sp vehicle planform area
Sw vehicle wetted area
T temperature
t generation index
U velocity magnitude
U uniform distribution
V vehicle internal volume
v number of design variables in design vector
w vehicle width, or post-optimality weights
X bounded design variable set
x body fixed longitudinal dimension, or design variable
x? true Pareto-optimal design set
x? approximated Pareto-optimal design set
~x design variable vector
y body fixed lateral dimension
z body fixed vertical dimension
z? ideal point
? oblique shock inclination angle
?? pressure gradient parameter
? small perturbation
?PM MOGA Pareto metric function
? leading edge inclination angle, or environmental parameter
perturbation
 uncertainty neighborhood
? volumetric efficiency
? vehicle wedge angle
? dynamic viscosity
? density
? standard deviation
? shear stress
? turbulence wedge half-angle, or probability density function
? boundary-layer stream function (transformed plane)
? ?h arithmetic mean over generations 1 through t
Subscripts
I robust regularization formulation
II explicit averaging formulation
III perturbation formulation
IV robust MOGA formulation
0 stagnation value
ab uncertainty neighborhood index
bad bad point in objective space
base vehicle base value
xii
CG about the center of gravity
cc concave
con after applying constraints
cv convex
e boundary-layer edge
ev minimized expected value
good good point in objective space
i panel index, or averaging index
i.e. illustrated example
j constraint function index
l lower surface
lam laminar
lb lower bound
local at a given spanwise location
m MOGA objective function index
max maximum allowed for attached shock
p planform
press due to pressure
q MOGA objective function index
scaled after scaling
span at a given value of p
t at generation t
tot total
turb turbulent
u upper surface
ub upper bound
visc due to viscosity
w wall
wcs minimized worst-case scenario
with with a turbulence wedge present
without without a turbulence wedge present
x based on x distance, or x component
y y component
z z component
? freestream conditions
Superscripts
? average boundary-layer conditions
? derivative in the transformed plane
? global optimum
Acronyms
ABLT asymmetric boundary-layer transition
AEDC Arnold Engineering Development Complex
CFD computational fluid dynamics
xiii
CG center of gravity
DNS direct numerical simulation
DOE design of experiments
GADS genetic algorithm and direct search
GNC guidance, navigation, and control
HTV-2 Hypersonic Technology Vehicle-2
L/D lift-to-drag ratio
MATLAB Matrix Laboratory
MAXWARP Maryland Axisymmetric Waverider Program
MOGA multi-objective genetic algorithm
MOOP multi-objective optimization problem
NASA National Aeronautics and Space Administration
NSGA-II Non-Dominated Sorting Algorithm-II
PDF probability density function
PLIF planar laser-induced fluorescence
SAEA Surrogate-Assisted Evolutionary Algorithm
STABL Stability and Transition Analysis for hypersonic Boundary-Layers
STS Space Transportation System
TPS thermal protection system
UMD University of Maryland
xiv
Chapter 1: Introduction
1.1 Motivation
Physical deformations, surface imperfections, or surface roughness on a hyper-
sonic vehicle can promote early localized boundary-layer transition. If transition
occurs asymmetrically, leading to drag asymmetry, it is essential that these effects
be small enough to be compensated by the control system of the vehicle [1]. If the
effects are too great, the asymmetries can perturb the vehicle to a state which is
uncontrollable, leading to vehicle loss.
After the space shuttle Columbia accident, as a part of the Return to Flight
effort, a focus was set on predicting hypersonic boundary-layer transition onset from
damage (e.g., tile impact, gap fillers) [2]. A lack of quality flight data and complete
understanding of the phenomena prompted a spacewalk to remove protruding gap
filler before the reentry of STS-114. Post STS-114, many space shuttle flight experi-
ments have been performed to investigate the impact of asymmetric boundary-layers
due to roughness-induced transition. During the flight experiments, it was found
that the shuttle?s flight control system had to respond to a roll/yaw moment pro-
duced by the experimentally induced asymmetric boundary-layer. Numerous firings
of control jets and deflections of control surfaces were necessary to return the shuttle
1
to the nominal attitude [2].
The importance of vehicle surface conditions and their effects on vehicle con-
trollability has also been identified for hypersonic cruise vehicles, as exemplified by
the HTV-2 (Hypersonic Technology Vehicle-2) flight tests. The second flight test of
HTV-2 ended in early vehicle termination. An independent review board concluded
that the most probable cause of the HTV-2 premature flight termination was an un-
expected aeroshell degradation, creating ?multiple upsets of increasing severity? [3].
History has shown that it is important to include this phenomena in the models and
methods used during the design process of hypersonic vehicles and projectiles.
The design process begins by listing all the mission requirements and con-
straints that are needed to be satisfied by the aerospace vehicle. The design space is
typically wide open, leaving the designers with the task of performing a preliminary
analysis in hopes to narrow the scope and provide a rough estimate of the vehicle
shape and performance characteristics. This is often accomplished by the use of
optimization methods. Traditional optimization methods tend to locate isolated,
specific design points with high precision. The resulting design can be sensitive to
uncertainty in the design variables or environmental parameters [4]. The process
of finding optimal solutions that are the least sensitive to uncertainty is referred to
as robust optimization. The problem of optimizing a hypersonic vehicle with asym-
metric roughness considerations inherently contains uncertainty. For example, the
position or severity of surface imperfections or roughnesses are not known a pri-
ori. These parameters may change due to high stress flight conditions. This design
problem is well suited for robust optimization.
2
Previous studies have used aerospace applications to evaluate robust optimiza-
tion methods. Ong et al. [5] used an airfoil design problem to demonstrate their
robust optimization scheme max-min SAEA. In a hypersonic framework, Ridolfi et
al. [6] used a capsule design problem as a test-bed for their archival-based robust
optimization scheme. Generally, robust optimization methods find an optimum by
constructing a robust counterpart to the original objective function. The robust
counterpart function can be formulated based on different philosophies and there-
fore may result in different categories of robust optima. For example, methods based
on robust regularization [7], are based on a worst-case scenario approach: given an
objective function f to be minimized, the robust optimum using robust regulariza-
tion is the design that minimizes the maximum value of f over the uncertainty [4].
The most common approach is to minimize the expected value of f over the uncer-
tainty. This is typically done through a Monte Carlo integration approximation to
the expected value of f . This method is an aggregation method sometimes referred
to as explicit averaging [6, 8, 9].
Both robust regularization and explicit averaging methods converge to a single
robust optimum design. Engineering design problems often benefit from comparing
the results from different robust optimization methods to help analyze the design
space and/or aid in a final design selection. To accomplish this, a separate optimiza-
tion problem needs to be formulated and solved for each method. If the designer
wants to analyze the problem for a different uncertainty probability density function
(PDF), then a separate optimization problem needs to be formulated and solved.
This is tedious as well as computationally expensive. This reveals a missing capa-
3
bility in state-of-the-art robust optimization methods, specifically when applied to
engineering problems. There is a need for an all encompassing robust optimization
method that can provide multiple types of robust optima for many uncertainty dis-
tributions without having to repeatedly resolve the problem. This will in turn pro-
vide engineers, scientists, or designers with information to fully analyze the design
space and understand the fundamental physics or principles behind the optimization
problem at hand.
Successfully optimizing hypersonic vehicles under asymmetric roughness-induced
boundary-layer conditions is essential to the progress of reentry vehicle technology
and hypersonic projectile technology. Advancing the state of the art requires a full
understanding of the design space. This can only be accomplished through accurate
and inexpensive modeling of roughness-induced boundary-layer phenomena coupled
with versatile optimization and robust optimization methods.
1.2 Objectives and Contributions
This work characterizes the aerodynamic sensitivities of hypersonic vehicles
to asymmetric conditions and explores concepts that balance aerodynamic capabil-
ity with control. To accomplish this task, an analytical waverider geometry model
capable of generating a wide range of concave and convex shapes was used. An
aerodynamic model was then used to approximate the aerodynamic forces and mo-
ments on the vehicles. Asymmetric flow properties were modeled to simulate the
effects of an isolated surface roughness event which may occur on a vehicle (e.g.,
4
protruding thermal protection system (TPS) tiles, ablation). Multi-objective opti-
mization studies were performed in order to identify vehicle shapes most robust to
the moment-induced effects of an asymmetric boundary-layer. Trade-offs between
controllability and vehicle performance were investigated.
In this research, a multi-objective robust counterpart formulation is proposed,
where post-optimality data handling techniques can identify, among a set of Pareto-
optimal solutions, the optimal worst-case scenario and expected value solutions for
any uncertainty PDF. Solving the robust multi-objective formulation once produces
a set of designs, thus giving the designer more information about the design space.
The ability to investigate the design space and compare multiple robust optimum
designs is a particularly important contribution to the field of engineering opti-
mization, where providing context with an optimum design is often necessary and
useful. This capability is a unique advancement in the state of the art and is not
currently available in modern robust optimization methods. The proposed robust
multi-objective optimization method was compared to three commonly used robust
optimization methods. The foundations behind the proposed method, the con-
vergence criteria, and data handling techniques are discussed in detail. A purely
mathematical example problem was used to explore the capability and nuances of
the proposed method.
The robust optimum results found using the data handling techniques were
compared to the true robust optimum results found analytically. The proposed
robust multi-objective method and the three common robust optimization methods
were all applied to the vehicle design optimization problem to identify hypersonic
5
shapes that are robust to uncertain surface roughness conditions. The results and
performance of each method were then compared and discussed. In particular, the
potential for the multi-objective based robust optimization method to reveal design
trends is illustrated through the hypersonic vehicle design problem.
In summary, the following were the main objectives of the study:
? Characterize the effects of asymmetric roughness-induced boundary-layer tran-
sition on hypersonic vehicle controllability.
? Identify hypersonic shapes best fit to withstand asymmetric flow field effects
through optimization techniques.
? Explore the hypersonic vehicle design space through state-of-the-art multi-
objective and robust optimization methods using genetic algorithms.
By accomplishing the objectives listed above, the following contributions have
been made through this work:
? An analytical roughness-induced boundary-layer transition model was devel-
oped to explore the unknown design space of hypersonic vehicles with uncer-
tain surface roughness conditions.
? Hypersonic vehicle shapes are classified for maximum aerodynamic perfor-
mance, vehicle controllability, vehicle usability, as well as hybrid shapes acting
as a compromise between these objectives.
? A novel robust multi-objective optimization method to solve single-objective
optimization problems with environmental parameter uncertainty by use of
6
post-optimality data handling techniques was developed. The method pro-
vides a range of robust optimum results, while also capturing trade-offs in the
design space, providing capabilities not available in current robust optimiza-
tion methods.
1.3 Organization of the Dissertation
The remainder of this chapter provides a review of the previous work per-
formed in the fields of hypersonic vehicle design, roughness-induced boundary-layer
transition, and state-of-the-art robust optimization methods. Background knowl-
edge and important concepts are discussed. An in-depth analysis of the hypersonic
aerodynamic prediction model and pertinent assumptions used in the simulations
are outlined in Chapter 2. Chapter 3 presents a detailed discussion on asymmetric
boundary-layer transition (ABLT) and discusses the implementation of this phe-
nomena into the aerodynamic model. The hypersonic aerodynamic model and the
ABLT model were applied to the vehicle design problem using the methods de-
scribed in Chapter 4. The design point, structure, and formulations of the opti-
mization problems solved in this work are given in Chapter 5. The optimum designs
and explanation of the results are fully outlined in Chapter 6 and summarized in
Chapter 7.
The content in this dissertation is a comprehensive piece of literature pulling
together previous publications by the current author while providing more detail.
References [10?12] form the foundation of this dissertation.
7
1.4 Previous Work
1.4.1 The Design, Aerodynamics, and Optimization of Hypersonic
Vehicles
The classic highly blunted, low lift-to-drag ratio (L/D), hypersonic vehicle
shapes have reliably served NASA for the past 50 years. Yet, looking forward, there
are particular missions of interest which may benefit from a fundamentally different
aeroshell shape. For example, the human exploration of Mars have advantages
from the use of mid-L/D vehicle designs due to the large landed mass requirements
[13]. Aerocapture missions to Neptune would benefit greatly from high-L/D shapes
by utilizing aerodynamic forces during a pass through the Neptune atmosphere to
capture into orbit, instead of a large propulsive delta V maneuver.
One of the most promising classes of high-L/D hypersonic vehicle design con-
cepts are waveriders. Waveriders are hypersonic vehicles that are designed such that
the entire bow shock is on the underside of the vehicle emanating from the leading
edge of the body. As a result, there is no flow spillage from the lower surface to the
upper surface of the vehicle. The shockwave creates a high pressure flow which is
contained beneath the vehicle, thus producing a high lift-to-drag ratio.
NASA has identified these high-L/D shapes as key research areas within the
Entry, Descent, and Landing Technology Roadmap [13] and are currently funding
research which are based of the waverider concept [14]. The goals of this proposed
research parallel the recommended areas of investment defined in the roadmap;
8
the primary category of shapes investigated in this work is waveriders (see Section
2.1). However, the analysis isn?t too restrictive as canonical hypersonic shapes
such as wedge-like, cone-like, convex, and concave vehicle shapes are all possible
candidates in the studies performed in this research. That being said, the previous
work presented in this section will mainly focus on previous literature pertaining
specifically to waveriders.
The concept of a waverider was first conceived in a paper by Nonweiler in
1959 [15]. His first waverider concept was derived by carving out a geometry based
on the streamlines and shockwave generated from a two-dimensional hypersonic
flow over a wedge. These designs were originally coined as caret waveriders. The
waverider concept was further extended to include shapes that are ?carved? out of
more complex generating flow fields, such as flow over a non-lifting cone [16]. The
generation of waverider shapes grew in complexity once again in the 1980s and 1990s,
fueled by a renewed interest in hypersonic flight. Sobieczky et al. [17] introduced
the method of osculating cones to generate waverider shapes with three-dimensional
shock structures. Most of the previous work on the generation of waverider shapes
are based of an inverse design approach: a generating flow field is first established,
then the waverider design is carved out of the resulting streamlines.
In 1998, Starkey and Lewis [18] took the opposite approach. An analytical,
power-law derived hypersonic waverider model was developed using two or three
power-law equations to define the upper, lower, and planform shapes of the wa-
verider. From this, five variables can be manipulated to generate a wide range of
waverider styles including: caret waveriders, constant wedge angle waveriders with
9
planar shocks, constant wedge angle waveriders with three-dimensional shocks, and
variable angle waveriders with three-dimensional shocks.
Once the waverider concept and shape generation methods were conceived,
the next logical task is to evaluate the aerodynamic characteristics of waveriders.
This task has been the subject of many research papers spanning from the late
1990s to the 2010s. The previous literature predicted waverider aerodynamics by
using numerical, analytical, and experimental techniques. Many papers compare
and contrast two or more of these approaches. The papers discussed in the next
few paragraphs are a select sample of the vast number of contributions to waverider
aerodynamics. These papers were selected to highlight the diversity in the aeropre-
diction approaches. An organized view of the theory behind most of these methods
is available in Anderson [19].
Many of the legacy research papers on predicting the aerodynamic character-
istics of waveriders took an inviscid analytical approach. Within the overarching
category of an inviscid analytical approach existed numerous theories and methods.
One of these theories is Piston Theory. This theory calculates pressure distribu-
tions along a waverider by assuming an analogous scenario of a piston moving into
a column of fluid. This theory can be seen applied to the calculation of stability
derivatives in references [20?23] as well as more recently in Oppenheimer et al. [24].
Corrections/additions to this theory have been provided by Liu et al. [25]. Liu
provides a hypersonic lifting surface theory, adding the capability to account for
the effects of wing thickness and/or flow incidence, upstream influence, and three-
dimensionality for an arbitrary lifting surface system in an unsteady flow.
10
Another subset of inviscid analytical methods are surface inclination meth-
ods. These methods theoretically show that surface pressures can be approximated
using only the information regarding the angle of incidence of the vehicle surfaces
with respect to the freestream. These methods include Newtonian Theory, the
tangent-wedge method, and the tangent-cone method. Rasmussen wrote several pa-
pers [26?28] outlining the approach and capabilities of using Newtonian Theory to
analytically calculate the stability derivatives of waveriders as a function of angle of
attack and sideslip angle. He also combines these efforts within the framework of
hypersonic-small-disturbance theory [29].
Efforts by Starkey et al. [30] were taken to use shock expansion theory to an-
alytically describe the aerodynamics of waveriders. His method and results were vali-
dated with numerical simulations using the waverider aeroprediction code MAXWARP
(Maryland Axisymmetric Waverider Program) [31]. Comparing analytical methods
with experimental and numerical results was a common theme throughout the lit-
erature.
O?Brien and Lewis [32] compared the analytical waverider aerodynamic predic-
tions of Newtonian Theory to a three-dimensional Reynolds-Averaged Navier Stokes
solver and found reasonable agreement. Cockrell et al. [33, 34] performed extensive
experimental work on two cone-derived waverider shapes with engine integrated
vehicle components. Comparisons were made to inviscid numerical computational
fluid dynamics (CFD) solutions. It was found that the aerodynamic performance
of the fully integrated configuration is severely degraded when compared to a pure
waverider. Gillum and Lewis [35] used hypersonic wind tunnel results to examine
11
the prediction capabilities of Newtonian theory and the tangent wedge method, and
to examine the effect of blunted leading edges.
The effects of blunting waverider leading edges is also seen as a topic of inter-
est in the literature. In 1994, Blosser et al. [36] presented waverider leading edge
design concepts. Blosser examines the techniques for blunting as well as the impact
on leading edge heat flux, via analytical and numerical aerothermodynamic calcula-
tions. Blunting concepts and the corresponding aerothermodynamic effects has also
been studied more recently by Starkey [37] and Chen et al. [38].
Research comparing analytical inviscid methods with numerical and experi-
mental results drew attention to the importance of viscous effects. Viscous effects
were found to be an important aspect to model, even in the high Reynolds num-
ber flow of the hypersonic flight regime. This implication is still emphasized in
the more recent waverider publications [39?41]. The importance of modeling shear
stress brought forth the task of designing waveriders with minimized viscous drag
characteristics. These waveriders were coined as viscous optimized waveriders [42].
The ability to accurately (and hopefully inexpensively) predict waverider aero-
dynamic characteristics helped guide the research towards optimizing waveriders for
particular missions. One of the first and foremost computer simulation tools used to
optimize waveriders was MAXWARP [31]. The MAXWARP code has been used to
generate aerodynamically optimized waverider shapes for maximum cruise mission
performance, reentry missions, chemically reacting flow fields, and interplanetary
missions [30]. Waverider optimization has been prevalent from the 1990s onward
and is still a hot topic of research. For example, Burt et al. [43] used a simulated
12
annealing optimization method to optimize waverider geometries for high altitude
flight, calculating the aerodynamics numerically using a Direct Simulation Monte
Carlo method.
An ongoing theme stressed in this study is the fact that traditional optimiza-
tion techniques tend to find optimum designs which are sensitive to off-design con-
ditions (either in design parameters or environmental parameters). The literature
shows that waverider optimization follows this tendency. However, it is important
to note that this sensitivity is not only attributed as a side-effect of optimization,
but also as a performance characteristic of waveriders in general. Much of the wa-
verider literature in the late 1990s and early 2000s studied the changes in waverider
aerodynamic performance for off-design conditions. Two main off-design categories
were investigated: off-design Mach number and off-design attitude (i.e. changes in
angle of attack and/or sideslip angle).
Tarpley et al. [22,44] presented a hybrid method to predict the aerodynamics
of waveriders for both on- and off-design conditions (including both Mach number
and attitude changes). Linear piston theory was used to calculate on-design aero-
dynamics coupled with tangent wedge and tangent cone methods for the off design
regime. The main assumption made was that the shockwave is still attached to
the leading edge, even in the off-design conditions. This assumption is somewhat
limiting, due to the fact that this would only remain true for small perturbations in
freestream Mach number and attitude angles.
For larger perturbations in freestream Mach number, Takashima et al. [42]
found that lift-to-drag ratio performance didn?t suffer significantly. The lift-to-drag
13
ratio for a waverider with a design Mach number of six changed by a maximum of
13.5% over a freestream Mach number spanning from four to eight. When extrap-
olating this effect to a hypothetical waverider mission, Takashima et al. [45] found
that when a waverider accelerates through an off-design Mach number to reach
cruise conditions, mission parameters like maximum range aren?t greatly affected.
Even though the overall parameters like lift-to-drag ratio and range aren?t greatly
affected, Mazhul and Rakchimov [46] found that the shock structures produced by
the waverider can in fact be vastly different in configuration for off-design Mach
numbers and angles of attack.
Changes in shock structure and flow physics around a waverider are even more
exaggerated when considering off-design sideslip conditions. This is seen in Shi et al.
[47,48], where CFD simulations were run to analyze waveriders at sideslip conditions.
Shi found that a waverider at sideslip conditions produces a complex asymmetric
flow field with intricate shock structures. Shockwave boundary-layer interaction was
found to occur for Mach numbers faster than the design Mach number. At an angle
of attack, the shock detaches and the presence of a small separation bubble and
a vortex on upper surface was observed. These flow structures strongly affect the
waverider pressure and skin-friction coefficient distributions.
Once the research on waverider aerodynamics simulation and experimentation
matured, the community started ramping up hypersonic flight tests [49]. Real-world
scenarios revealed more concern about the stability and control of waveriders and
hypersonic vehicles. The cause of concern is asymmetric boundary-layer transition
[1]. In particular, roughness-induced boundary-layer transition was identified as one
14
of the most important aspects to study [2, 50,51]. Previous literature pertaining to
roughness-induced boundary-layer transition is reviewed in Section 1.4.2.
1.4.2 Roughness-Induced Boundary-Layer Transition
Boundary-layer transition is a complex phenomena influenced by many con-
tributing factors including, but not limited to, Reynolds number, crossflow, freestream
turbulence levels, wall temperature, and vibrations [52]. Among all factors, it has
been found that transition can become completely dominated by roughness or abla-
tion effects [53]. As expected with any topic pertaining to turbulence, the research
community has been studying roughness-induced boundary-layer transition for a
long time, without any sign of slowing down. This may be contributed to the com-
plex nature of the problem coupled with the continuous improvements of ground
test facilities, instrumentation, and numerical simulations.
This section has no intention of providing a complete review on all of the
previous research involving this topic. The history is simply too long and too vast.
Highlights of several historically significant papers will be given and an emphasis
will be placed on discussing the status of the state-of-the-art. For a more complete
background and history of the literature, a survey paper on hypersonic roughness-
induced boundary-layer transition has been written by Schneider [53]. Also, Zhong
and Wang [54] provide a thorough discussion on the background and literature
pertaining to hypersonic boundary-layer receptivity, instability, and transition.
An analysis of roughness elements disturbing laminar flow and inducing tran-
15
sition can be seen published in the literature as far back as the early 1930s, where
Schiller suggested that ?the critical height of the excrescence that disturbs the lami-
nar boundary-layer and causes its transition to turbulent flow is the height at which
vortices form behind the excrescence? [55]. Early experimental efforts to visualize
and explain the physical phenomena responsible for forming the ?wedge? of tur-
bulence behind the roughness were completed by Gregory and Walker [56]. These
visualization-based studies were extrapolated to hypersonic conditions for conical
reentry shapes in 1968 by Canning et al. [57]. Ablative visualization techniques
were used, thus by design, highlighting the importance of turbulence wedge local
heat flux implications.
The flow structures produced by an isolated roughness in a laminar flow field
have stemmed vast amounts of research in of itself. For example, Acarlar and
Smith [58] performed experimental studies to exhaustively illustrate and explain
the standing horseshoe vortex and shedding hairpin vortex structures produced by
an isolated roughness element. Current research studies are still focused around
explaining these structures and their subsequent breakdown into a turbulence wedge
[59].
The previous literature has also focused on the turbulence wedge downstream
of the roughness element. In 1972, Fischer [60] compiled lateral spreading angle
data from a wide variety of turbulence wedge experiments. His turbulence wedge
spreading data correlation with freestream Mach number is often cited in the previ-
ous literature as a validating reference. Even 40 years later, current numerical and
experimental papers cite Fischer to validate their measurements [53, 61?63].
16
Visualizing the turbulence wedge structure is a common approach in both
experimental and numerical studies, due to the inherent three-dimensional nature
of the flow physics. Experiments performed by Danehy et al. [64, 65] employ pla-
nar laser-induced fluorescence (PLIF) techniques to visualize the wall normal and
planform planes of flat plate experiments. Numerical validation of the results pub-
lished by Danehy et al. using direct numerical simulation (DNS) was performed by
Iyer et al. [66], further illustrating the interest in roughness-induced boundary-layer
transition from both the experimental and numerical research communities.
A common approach taken in the most current literature, is to analyze the nu-
merical stability of the hypersonic boundary-layers. A state-of-the-art suite of codes
called STABL and STABL-3D was developed by Candler et al. at the University
of Minnesota [67]. This computer code solves hypersonic boundary-layer transi-
tion problems using both linear stability theory and parabolized stability equation
solutions.
1.4.3 Robust Optimization
The concepts behind robust optimization have been developed independently
in many scientific fields, yet can arguably be traced back the furthest in engineering,
specifically the field of ?quality engineering? [4]. The roots of robust optimization in
quality engineering were most strongly planted by Taguchi in the mid 1980s. Taguchi
laid the foundation of the techniques and philosophies for quantifying robustness.
Taguchi followed a Design of Experiments (DOE) approach, rather than formal
17
optimization methods, to carry out the process and find robust designs [68].
Since the 1980s, many disciplines have built off of the idea of robustness and ap-
plied formal analytical and numerical optimization techniques, making great strides
advancing the field. This section serves as a select literature review of the state-of-
the-art methods and techniques in robust optimization. Comprehensive literature
reviews and survey papers have been written by Beyer and Sendhoff [4], Jin [8] and
Kruisselbrink [9] and are recommended for further review.
1.4.3.1 Formulations
Generally, robust optimization methods incorporate uncertainty by construct-
ing robust counterpart objective functions F to the original objective function f [4].
This counterpart function is based off of a particular ?philosophy? or statistical rep-
resentation of uncertainty. For example, worst-case scenario methods use a minimax-
type formulation in order to construct a robust counterpart function. Aggregation
methods form robust counterparts based on the mathematical moments of f over
the uncertainty.
The worst-case scenario form of accounting for uncertainty is directly analo-
gous to the minimax approximation in the field of applied mathematics [69]. This
method appears in research spanning multiple disciplines and may be referred to as
minimax, robust regularization [7], or simply as a robust counterpart approach [70].
For a minimization optimization problem, methods based on robust regulariza-
tion minimize the maximum f value over the uncertainty. The robust regularization
18
technique is nominally used when uncertainty is modeled deterministically. For op-
timization problems where the uncertainty is modeled probabilistically, probability
measures are provided describing the likelihood for a particular value of the un-
certain variable to occur in the domain of interest [4]. Probabilistically modeled
uncertainty problems can also be solved with this method, as long as an upper and
lower bound on the uncertain variable can be defined.
The robust regularization approach is considered as a worst-case scenario and
may be too conservative for some applications. Sometimes the worst-case scenario
robust result is infeasible for a particular problem. In such applications, a proba-
bilistic aggregation approach to forming a robust counterpart function is often used.
Utility aggregation functions such as expected value, variance, and other higher
order moments of f can be used. For many optimization problems, the objective
function f isn?t in analytical form, and therefore numerical representations of these
utility functions need to be used. Seemingly the most common utility function used
in robust optimization is the conditional expectation of f [5, 6, 8, 9, 71?75]. The
simplest way of calculating the expected value of the original objective function for
a particular distribution of uncertainty is through explicit averaging [6, 8, 9].
In the context of applying evolutionary or genetic algorithms, Tsutsui et al.
[72,73] developed a more efficient variation of the explicit averaging method. Explicit
averaging was reduced from taking an average fitness over n samples to adding a
stochastic perturbation to each design in the population for a given generation.
Much of the previous work using this perturbation method covered uncertainty in
the design variables, while a few others [5, 74] have found applications using other
19
types of uncertainties. For an infinite population size, it has been shown that the
perturbation method is equivalent to the explicit averaging method [4,8,74,75]. Also
note that if n = 1 and the perturbation (?p, ?x, ?f , etc.) is chosen randomly
per generation, then the explicit averaging method is equivalent to the perturbation
method.
Another option investigated in the literature is a robust meta-model approach.
This approach ties in robustness by constructing a meta-model approximation and
solving for an estimated robust design. Response surface methodology, neural net-
works, and Kriging models have all been proposed. Two good reviews on meta-
modeling techniques have been published by Simpson et al. and Jin [76,77].
Beyer and Sendhoff [4] note several problems about the effectiveness of the
meta-modeling approach:
1. These techniques are not well suited for large-scale robust optimization prob-
lems, when the number of design variables is large.
2. For a fully quadratic response surface, the meta-model must be repeatedly
applied in order to get closer to the robust optimum. It is unlikely that the
quadratic meta-model represents the robust counterpart adequately.
3. The computational efficiency of the meta-model approach is at least question-
able. It is not clear what the response surface really produces as iteration
time approaches infinity. Up to now, no proof has been presented that this
methodology yields a robust optimum.
Even with these weaknesses, the method is still applied in the literature. If the
20
number of design variables are small, and a simple response surface is used, the
meta-modeling technique can be computationally efficient as seen by Elster and
Neumaier [78].
The robust counterpart approach (with or without meta-modeling) form an
optimization problem aimed at minimizing a utility function (i.e. expected value)
accounting for the uncertainty. It has been identified in the literature that optimizing
the expected value may lead to a design with a large variance; this is not always
desirable [4]. Some of the previous literature address this problem by using an a
priori aggregation scheme [79?81]. These methods assign weights of importance
to variance and expected value measures, aggregating them into a single-objective
problem. The current state-of-the-art robust optimization methods address this
problem by forming multi-objective problems using the variance and expected value
as conflicting objective functions [4, 9, 74].
The multi-objective formulations have been solved using both deterministic
and heuristic methods. For example, Das [82] solves for a Pareto-optimal front
between variance and expected value objective functions using a deterministic non-
linear programming method. Chen [83] and Rolander et al. [84] formulate their
multi-objective problem in a compromise programming (or sometimes called goal
programming) structure.
Solving the robust multi-objective formulations using genetic algorithms is a
budding research area. Jin and Sendhoff [74] solve a two objective function optimiza-
tion problem using genetic algorithms. The objective function f and the variance of
f over the uncertainty act as the two conflicting objectives. Ray [85] adds one more
21
objective forming a three objective function robust optimization problem, solved
using evolutionary algorithms. Here the objective function f , the mean value of f
and the variance of f act as the three conflicting functions. An interesting extension
to this category of heuristic-based solvers can be seen in Sorensen [86]. Sorensen
uses a meta-heuristic method, called tabu search, to search for robust optimum so-
lutions. These state-of-the-art multi-objective formulations and methods are only
a first step towards providing the designer with more insight into the nature of the
robust design space. These methods are still limited to one particular uncertainty
PDF for one nominal value of the uncertain variable. The methods are also lim-
ited to providing robust optimum answers based on moments of the uncertainty
(expected value, variance, etc.) or worst-case scenario robust optimum answers.
The literature review completed for the current study found no methods that can
concurrently provide both of these types of robust optima.
The current work developed a novel robust multi-objective optimization for-
mulation and method which allows for post-optimality data handling techniques to
find robust optimum designs that follow the worst-case scenario and best mean fit-
ness approaches of robust regularization and explicit averaging. This can be seen as
an a posteriori aggregation method to contrast the aforementioned a priori aggre-
gation methods. Ideally, this method has the ability to post-optimally pick designs
that are most robust for any uncertainty PDF and for various nominal values or
the uncertain parameter. Therefore, many types of results can be identified and
analyzed by solving only one multi-objective optimization problem.
22
Chapter 2: Hypersonic Aerodynamics
The aerodynamic model was constructed using multiple different aerodynamic
approximation theories and methods. Certain methods are only applicable to vehi-
cles with an infinitely sharp leading edge, while others must be used for a blunted
leading edge. The model is discussed in two subsections based on the leading edge
conditions. For both leading edge conditions, the flow over the waverider is assumed
to be two-dimensional with streamlines flowing straight back over the vehicle. This
assumption was also made by Starkey in his analytical calculations [30] and vali-
dated with CFD results [87]. The aerodynamic model for waveriders with infinitely
sharp leading edges is based off the work presented by Starkey [30]. The analysis
of blunted leading edges, however, required several modifications and extensions
in order to render the approach applicable to a blunted leading edge. Gillum and
Lewis [35], using experimental measurements, found that blunted waverider aerody-
namic coefficients can be accurately predicted using tangent-wedge and Newtonian
theory approximations. The results found by Gillum and Lewis serve as validation
to the approach taken in this work. The overall aerodynamic analysis in this work
is implemented as a panel method, where the waverider geometry is constructed as
a mesh and the aerodynamic calculations are applied to differential area elements.
23
2.1 Vehicle Geometry
Starkey [30] developed an analytical description of waverider geometries that
uses three, two-dimensional power-law equations to generate a parametric hyper-
sonic vehicle with a three-dimensional shock. The curvature of the planform (p)
and upper surfaces (u) are defined by
yp = Ax
n , yu = Bz
n
u . (2.1)
The curvature of the lower surface (l) can be defined either for a convex vehicle (cv)
or a concave vehicle (cc) by
yl,cv = Ccv (?zl,cv + x tan ?)
m , yl,cc = Ccc (zl,cc ? x tan ?)
m . (2.2)
These equations allow for the generation of a wide variety of vehicle geometries. An
example of a convex vehicle derived using these equations is shown in Figure 2.1,
and a concave vehicle in Figure 2.2. By varying the constants A, B, C, m, and n
in Equations 2.1 and 2.2, geometries from highly convex shapes to highly concave
shapes, representing a range of low lift-to-drag ratio shapes to high lift-to-drag ratio
shapes can be modeled. This methodology for defining a vehicle geometry with
infinitely sharp leading edges was used in the current study.
A modification to Starkey?s method was made to allow for a parametric de-
scription of leading edge bluntness. There are two ways to blunt the leading edge of
an initially sharp vehicle: one is to remove or ?sandpaper? the leading edge down
to a desired shape, the second is to separate the upper and lower surfaces and add
material to create the desired leading edge. The latter method was chosen here
24
x y
z
(a) Isometric view
y = Axn
w
L
(b) Planform view
y = Bzn
y = C (x tan ? ? z)m
??
(c) Rear View
?
?
(d) Side View
Figure 2.1: Example convex parametric waverider shape.
x y
z
(a) Isometric View
y = Axn
w
L
(b) Planform View
y = Bzn
y = C (z ? x tan ?)m
??
(c) Rear View
?
?
(d) Side View
Figure 2.2: Example concave parametric waverider shape.
25
as suggested by Starkey [37]. An upper and lower surface separation distance, d,
was added to the list of variables defining the vehicle geometry. The side cross sec-
tion of the blunted leading edge was prescribed as a semicircle. The radius of the
semicircle varies with spanwise location (y direction) such that the leading edge is
always tangent to the lower surface. An example blunted waverider shape is shown
in Figure 2.3. The profile shapes of the infinitely sharp and blunted waveriders at
each spanwise location are a wedge and a blunted wedge as seen in Figure 2.2(d)
and Figure 2.3(d) respectively.
x y
z
(a) Isometric View
y = Axn
w
L
(b) Planform View
y = Bzn
y = C (z ? x tan ?)m
d
(c) Rear View
d
(d) Side View
Figure 2.3: Example blunted concave parametric waverider shape.
26
2.2 Infinitely Sharp Leading Edge
2.2.1 Pressure Forces
The vehicle geometry modeling in Section 2.1 generates a wedge profile at each
spanwise location on the vehicle, with the top surface oriented parallel with the
freestream. When assuming two-dimensional flow, the top surface flow properties
are equal to the freestream. A constraint was placed in the optimization schemes
(see Chapter 5) limiting the wedge deflection angle of the bottom surface at any
given spanwise location ?local, to be less than the maximum deflection angle ?max
allowed for an attached shockwave. With these assumptions, the lower surface
properties can be calculated using oblique shock theory. Flow properties on the
top and bottom surfaces of the waverider are approximated using perfect gas and
isentropic relations.
The vehicle geometry is discretized into a mesh with a finite number of differ-
ential areas or panels. Examples of these meshes are seen in Figures 2.1 and 2.2. A
local deflection angle and a local shockwave angle exists at each spanwise location.
Oblique shock relations are used to find the properties on the lower surface of the
waverider. The lift, drag, and side forces are calculated by multiplying panel surface
pressure by the corresponding differential area, splitting the force into components
as determined by the unit normal direction of the panel, and summing over all pan-
els. The base pressure is assumed to be equal to the freestream pressure. A profile
cross-section view of the pressures at a given span can be seen in Figure 2.4.
27
Oblique shock
U?, P?
Pu = P?
Pl
Pbase = P?
?local < ?max
?local
x
z
Figure 2.4: Profile cross-section view at some given local span location.
2.2.2 Viscous Forces
Laminar and turbulent viscous forces for the top and bottom surfaces of the
waverider are calculated using the reference temperature method after applying a flat
plate assumption to the compressible boundary-layer equations [88]. The laminar
shear stress at the wall is given by
?w,lam ? 0.332?eUe
2
?
Ch?
Rex,e
, (2.3)
where subscript e represents boundary-layer edge conditions and superscript ? rep-
resents average boundary-layer conditions. The boundary-layer edge conditions are
given by the flow properties calculated in Section 2.2.1. The Chapman-Rubesin
parameter Ch? is found using the power-law viscosity approximation given by
Ch? =
(
T ?
Te
)? 13
, (2.4)
28
where T
?
Te
is found using Eckert?s empirical estimate for the average boundary-layer
temperature given by
T ?
Te
? 0.5 + 0.039M2e + 0.5
Tw
Te
. (2.5)
The wall temperature Tw is approximated as 1200K. This is a reasonable estimate
for hypersonic vehicles at cruise conditions assuming active wall cooling [30].
The turbulent shear stress at the wall is found by utilizing the power-law
approximation for flat plates [88] given by
?w,turb ?
0.0592??U2e
2 (Re?)0.2
, (2.6)
where
Re? =
??Uex
??
, (2.7)
?? =
Pe
RT ?
. (2.8)
The viscosity ?? is calculated using Sutherland?s formula at temperature T ?. Once
the wall shear stress is calculated for each panel, the viscous force can be found by
multiplying the shear stress with the corresponding differential area. A small angle
approximation is used, assuming the friction force only contributes to drag in the x
direction and therefore has zero contribution to lift or side force.
2.2.3 Moment Calculations
Contributions to the roll, pitch, and yaw moment coefficients are given by both
the pressure forces and the viscous forces. The pressure force moment coefficient
vector is defined as
~CMCG,press ?
~MCG,press
1
2??U
2
?SpL
, (2.9)
29
where Sp is the planform area, L is the length of the vehicle, and ~MCG,press is the
moment vector about the center of gravity due to pressure forces given by
~MCG,press =
N?
i=1
(PdA)i (~ri ? n?i) , (2.10)
where subscript i represents a given panel in the waverider mesh, N is the total
number of panels making up the entire surface (top and bottom) of the waverider,
dA is the area of the panel, n? is the outward unit normal of the panel, and ~r is a
position vector from the vehicle center of gravity to the centroid of the panel.
The viscous force contribution to the moment coefficient vector is found by
applying a small-angle approximation. The wedge angles are assumed to be small
enough that the shear stress on the lower surface can be approximated as pointing
purely in the drag direction. Therefore, the moment coefficient is defined as
~CMCG,visc ?
~MCG,visc
1
2??U
2
?SpL
, (2.11)
where ~MCG,visc is the moment vector about the center of gravity due to viscous
forces given by
~MCG,visc =
[
MxCG,visc MyCG,visc MzCG,visc
]
, (2.12)
where
MxCG,visc = 0 ,
MyCG,visc =
N?
i=1
(?wdArz)i ,
MzCG,visc =
N?
i=1
(?wdAry)i .
(2.13)
The variables ry and rz are the distances from the center of gravity to the panel
centroid in the y and z directions respectively. MxCG,visc is equal to zero because the
30
shear stress points strictly in the x direction. The total moment coefficient vector
(x, y, and z components are roll, pitch, and yaw respectively) is the sum of the
pressure and the viscous contributions
~CMCG,tot = ~CMCG,press + ~CMCG,visc . (2.14)
2.3 Blunted Leading Edge
2.3.1 Pressure Forces
The geometry model in Section 2.1 for a blunted leading edge waverider gen-
erates a blunted wedge profile at each spanwise location on the vehicle, with the top
surface oriented parallel with the freestream, as illustrated in Figure 2.5. As op-
posed to the infinitely sharp leading edge, where an attached shockwave is present,
a blunted leading edge produces a bow shock. Predicting the bow shock shape and
calculating the flow field downstream of the shock is a ?daunting task? and is of-
ten calculated using CFD [19]. The goal of the current work was to construct a
more computationally inexpensive aerodynamic prediction model so that optimiza-
tion schemes can be implemented. Therefore, a more appropriate approach was to
use surface inclination methods.
Modified Newtonian Theory is most suited for blunt surfaces where the deflec-
tion angle is greater than the max deflection angle allowed for an attached oblique
shock wave. For more slender shapes, tangent-wedge, tangent-cone, or shock expan-
sion method is more accurate [19]. For the blunted wedge shapes of interest in this
work, a hybrid approach was taken. For the blunted leading edge section, region
31
Bow shock
U?, P? Pu = P?
Pbase = P?
?local
x
z
I
II
Figure 2.5: Profile cross-section view at some given local span location.
I in Figure 2.5, Modified Newtonian Theory is applied to approximate the surface
pressure. For the lower surface, region II in the figure, the local inclination angle
is checked to see if it is greater than or less than the max deflection angle ?max. If
the local wedge angle is greater than the max deflection angle, Modified Newtonian
Theory is applied to approximate the surface pressure. If the local deflection angle
is less than the maximum deflection angle, then a tangent-wedge approximation is
used. This hybrid approach was coupled with the perfect gas and isentropic rela-
tions to approximate additional flow properties on the leading edge and the bottom
surface of the waverider. The lift, drag, and side forces are calculated by multiplying
each panel surface pressure by the corresponding differential area, splitting up the
force into components as determined by the unit normal direction of the panel, and
then summing over all of the panels. As in the infinitely sharp leading edge case,
the top and base surface flow properties are assumed to be equal to the freestream
conditions.
32
Using one or many different surface inclination methods to approximate pres-
sures is common in the preliminary design of hypersonic vehicles [19]. Yet, surface
inclination methods do not model flow spillage or leaking that is known to occur for
blunted waverider bodies [38]. Past studies have shown that even with an apprecia-
ble bow shock stand-off distance, significant spillage does not occur [36]. For this
study, flow spillage was assumed negligible.
2.3.2 Viscous Forces and Leading Edge Heat Flux
The laminar and turbulent viscous drag calculations for the upper and lower
surface of the blunted waverider are performed using the methods outlined for the
sharp leading edge waverider, but now the boundary-layer edge conditions (denoted
by subscript e in the equations) are calculated using the hybrid surface inclination
method described in the previous section.
The flat plate assumption inherent in Equations 2.3 and 2.6 no longer holds for
the curved blunt leading edge surface. Solutions to the boundary-layer equations for
the blunted leading edge will have a dependence on a pressure gradient parameter
??. Calculations of shear stress and heat flux along the leading edge are performed
using a method formulated by Bae and Emanuel [89, 90]. This method enables all
the boundary-layer properties to be determined by a set of two-dimensional lookup
tables as a function of ?? and gw. The parameter gw is the temperature ratio defined
by
gw =
Tw
T0e
. (2.15)
33
The method uses laminar compressible similarity equations under the assumptions of
steady flow of a perfect gas with unity values for the Prandtl number and Chapman-
Rubesin parameter given by
???? + ???? + ??
[
gw + (1? gw)H ? ?
?2
]
= 0 , (2.16)
H ?? + ?H ? = 0 . (2.17)
The bounding wall is assumed impermeable and the flow is assumed two-dimensional.
Numerical solutions for Equations 2.16 and 2.17 were calculated by Emanuel for a
range of ?? and gw values. Resulting data for ???w and H
?
w was compiled into two-
dimensional tables [89].
?? and gw values are calculated for each panel on the blunted leading edge of
the waverider. By interpolating from the tables described above, ???w and H
?
w values
are calculated. Shear stress at the wall, ?w, and heat flux at wall, qw, were then
found by transforming the ???w and H
?
w values back to the physical plane. The viscous
drag is found by calculating the ?w component in the x direction and multiplying
by the corresponding differential area of each panel, then summing over all panels.
The sum of the viscous forces in the z direction is assumed to be negligible due to
near symmetry above and below the stagnation point on the leading edge. Figure
2.6 below exemplifies the breakdown of ?w into the x and z components ?wx and ?wz ,
respectively.
34
??
?
Bow shock
?local
?w
?wx
?wz
?w
?wx
?wz
?w
?wx
?wz
Figure 2.6: Shear stress components in the x and z directions.
2.3.3 Moment Calculations
Contributions to the roll, pitch, and yaw moment coefficients from pressure
forces for the top, bottom, and leading edge surfaces are given by Equation 2.9. The
viscous contribution for the top and bottom surface is given by Equation 2.11. The
leading edge viscous contribution to the moment coefficients are calculated using
Equations 2.11 through 2.13, but using ?wx instead of ?w.
35
Chapter 3: Asymmetric Boundary-Layer Transition
In hypersonic vehicle design, modeling boundary-layer transition is important
for prediction and control of heat transfer and skin friction on the vehicle. By
extension, predicting boundary-layer transition is crucial in designing a vehicle?s
TPS and GNC (guidance, navigation, and control) system. Transition is a complex
phenomena influenced by many contributing factors including, but not limited to,
Reynolds number, crossflow, freestream turbulence levels, wall temperature, and
vibrations [52]. Among all factors, it has been found that transition can become
completely dominated by roughness or ablation effects [53]. This allows the current
work to take an approach that assumes, somewhere along the surface of the vehicle,
an isolated roughness is present and the effects are strong enough to dominate
transition.
3.1 Roughness-Induced Transition
There appears to be at least three fundamentally different ways by which
roughness affects transition [53]:
1. Roughness generates a wake with streamwise vorticity and a possibly unstable
shear layer. This wake may transition to turbulence immediately behind an
36
effective roughness.
2. Streamwise vorticity behind a small roughness element can grow via instability
mechanisms such as stationary crossflow, Go?rtler, or transient growth. These
instabilities then lead to transition.
3. Roughness can interact with acoustic waves or freestream disturbances gener-
ating instability waves.
The first of these three modes is the most understood and was the focus for devel-
oping the model in the present study. The goal of the current work was to model
the roughness-induced transition as computationally inexpensively as possible while
retaining as much physical accuracy as possible. This was done by analyzing the
physical effects of surface roughness on transition, identifying any parametric trends
that exist, then establishing simplifying assumptions. The following sections out-
line this approach and give a summary of the physical effects, parametric trends,
and the model which was developed. For a more detailed analysis, a comprehensive
literature review of the effects of roughness on hypersonic boundary-layer transition
was written by Schneider [53].
3.1.1 Physical Effects
Surface roughness is usually categorized into two types: isolated or distributed.
Isolated roughness includes steps, gaps, joints, or discrete machining flaws. Dis-
tributed roughness may include machining or finishing patterns on a surface or
screw threads [53]. The present study focused on the effects of an isolated rough-
37
ness since more experimental and computational work has been carried out with
this type.
An isolated surface roughness promotes boundary-layer transition producing a
turbulence wedge downstream of the disturbance. The establishment of a turbulent
wedge starts with a single disturbance source in a laminar boundary-layer. A stand-
ing horseshoe vortex system develops around the roughness element and trails some
distance downstream until a breakdown into turbulence occurs. The breakdown
produces a near perfect symmetrical pattern of turbulent flow. Unsteady structures
within the turbulent region break down into smaller substructures which grow in
the spanwise direction as they propagate downstream [59, 66]. Turbulence wedges
generated by the presence of an isolated surface roughness element exist in subsonic,
supersonic, and hypersonic flows. Examples of typical turbulence wedges produced
by an isolated roughness in a hypersonic flow field are shown in Figure 3.1.
3.1.2 Parametric Trends
3.1.2.1 Roughness Height and Effective Roughness
In general, roughness affects boundary-layer transition by moving the transi-
tion point upstream of the transition location for a smooth wall. The amount of
upstream movement increases with roughness height. At a particular height, the
transition location will be located close to the roughness element, and any further
increase in height will cause a larger disturbance but will not move the location any
further upstream. This roughness condition is termed effective [53].
38
(a) Instantaneous velocity contours at a wall normal
plane (M? = 3.37). [66]
(b) Sublimation pattern and turbulent heat flux of
wake past isolated roughness (Me = 2.9). [91]
Figure 3.1: Example turbulence wedge generated from isolated surface roughness
[66,91].
The effective roughness height is dependent on several factors including boundary-
layer edge Mach number, position of the roughness with respect to the leading edge,
and wall temperature. The present work assumes that any disturbance modeled on
a vehicle will be effective. Shocks are regularly seen in schleiren images for roughness
elements of heights larger than the boundary-layer thickness [53]. A key assump-
tion in the current work is that the roughness height will be on the order of the
boundary-layer thickness or smaller. Therefore, no additional shocks as a result of
the roughness were modeled .
39
3.1.2.2 Roughness Shape
Effective isolated roughnesses promote boundary-layer transition efficiently,
yet transition may still occur at a certain distance downstream of the trip location.
The distance between the trip location and transition location is found to be affected
by factors such as boundary-layer edge Mach number and the roughness shape.
Tirtey et al. [91] experimentally produced turbulence wedges using a variety
of isolated roughness shapes shown in Figure 3.2(a). Infrared thermography mea-
surements were taken to characterize the heat load distribution in the wake of the
roughness. A modified Stanton number distribution along the centerline of the tur-
bulence wedge was calculated for each roughness shape. Resulting distributions
were compared to laminar and turbulent theoretical results (calculated using Eck-
ert?s reference temperature method). The comparison plot is given in Figure 3.2(b).
The dashed lines show that for all cases, areas outside of the turbulent wedge wake
are well approximated by laminar flat plate theory. The solid lines show that up-
stream of the disturbance, the flow is also approximately laminar. At the centerline
of the turbulence wedges downstream of the roughness, the solid lines show that the
flow will eventually fully transition to turbulence. When fully turbulent, the flow is
well approximated by theoretical calculations. For the same freestream conditions,
the shape of the roughness strongly influences how quickly or effectively the flow is
tripped from laminar to turbulent.
The model constructed in the present work did not assume any particular
roughness geometry. Rather, an assumption was made that at a particular location
40
on the surface of the vehicle, an isolated roughness disturbance exists such that flow
will transition to a fully turbulent boundary-layer immediately downstream.
(a) Roughness element geometries
(b) Modified Stanton number distribution along wake axis compared to laminar and
turbulent heating theoretical predictions
Figure 3.2: Effect of roughness shape on transition from Tirtey et al. [91]
3.1.2.3 Turbulence Wedge Spreading Angle
The spanwise growth of the turbulent region as it moves downstream of a sur-
face roughness is measured by a lateral spreading angle. The lateral spreading angle
is defined as the half-angle of the turbulence wedge assuming symmetrical growth in
the lateral direction. Fischer [60] shows a parametric trend existing between Mach
number and spreading angle. Although not pursued in the current work, the prob-
41
lem still remains to be studied to rule out any other parametric trends that may
exist (e.g., possibly with wall temperature [63] or tunnel noise [53]).
Fischer compiled lateral spreading angle data from a wide variety of turbulence
wedge experiments [60]. Variation of the lateral spreading angle with boundary-layer
edge Mach number is presented in Figure 3.3. The figure shows that with increasing
boundary-layer edge Mach number, the spreading angle decreases.
Figure 3.3: Variation of turbulent spreading with boundary-layer edge Mach number
adapted from Fischer [60].
The turbulence wedge model constructed for the current work (described in
Section 3.1.3) requires a lateral spreading angle to define the area on the waverider
where the turbulence wedge exists. This spreading angle was approximated by the
value of the data point that vertically bisects the shaded region of Figure 3.3 at a
given Mach number. For example, at a Mach number of 4, the shaded region spans
from 3.0? to 5.0?. So, the resulting spreading angle for a design Mach number of 4
would be approximated as 4.0?.
42
3.1.3 Turbulence Wedge Model
One of the goals in this study is to inexpensively model the aerodynamics,
in efforts to employ optimization methods. To accomplish this, roughness-induced
boundary-layer transition was not modeled using any CFD or numerical boundary-
layer transition theory (e.g. Johnson and Candler [67]). To model the influence of
surface deformations or roughness on a hypersonic vehicle, several simplifying as-
sumptions are made. The disturbance is assumed to be severe enough to fully trip
the boundary-layer immediately downstream of the isolated roughness. Therefore,
the vertex of the turbulence wedge coincides with the disturbance point. The tur-
bulence wedge centerline is oriented parallel to the vehicle body x-axis. The Mach
number dependence of the turbulence wedge is included in the model using a tur-
bulence wedge half-angle parameter. This angle is defined by the curve presented
by Fischer [60]. The boundary-layer outside of the turbulence wedge is modeled as
completely laminar, while the flow inside the turbulence wedge is modeled as fully
turbulent; no transition region is defined. The location and geometry of the turbu-
lence wedge can now be described by four parameters: x position, y position, upper
or lower surface, and turbulence wedge half-angle. Figure 3.4 shows an example of
the parametrically described asymmetry as applied to an example vehicle geometry
from Section 2.1.
43
?x y
z
Isolated roughness
Turbulence wedge
half-angle (?)
Fully turbulent
Fully laminar
(a) Isometric view
?
y
x
Isolated roughness
Turbulence wedge
half-angle, (?)
Fully turbulent
Fully laminar
(b) Planform view
Figure 3.4: Example of an isolated roughness element applied to upper surface of a
generic hypersonic vehicle.
44
Chapter 4: Optimization Methods, Algorithms, and Definitions
This chapter presents the optimization methods, algorithms, and definitions
that were used in this work. All optimization problems formed in this study were
solved using genetic algorithms to help ensure global optimality and to remove the
requirement of numerically calculating derivative information. The use of genetic
algorithms were possible due to the relatively inexpensive nature of the objective
functions.
A brief overview of traditional genetic algorithm optimization methods is pro-
vided in Section 4.1. This section gives the basics of genetic algorithms as well as
the current most established methods and concepts within a multi-objective genetic
algorithm context. One of the main contributions in this work is the formulation of a
novel robust optimization method. This work provides comparisons between the per-
formance of the proposed method to common robust optimization techniques found
in the literature. A detailed discussion on these robust optimization methods is
provided in Section 4.2. To illustrate the use of the proposed robust multi-objective
method, an illustrated example is provided in Section 4.3. In particular, Section 4.4
focused on analyzing the changes in convergence and accuracy of the method with
changes in the genetic algorithm initial conditions and data handling settings.
45
4.1 Traditional Optimization Methods
4.1.1 Genetic Algorithms
This section serves as a cursory introduction into genetic algorithms. It will
provide the basic definitions of the terms used in this study as well as a top level
explanation into the process of solving an optimization problem using a standard
genetic algorithm. An in-depth discussion on the history, theory, and use of genetic
algorithms is available in a book by Mitchell [92].
The theory behind genetic algorithms stems from the biological process of
evolution. In nature, as a population of a particular species reproduces and survives
from generation to generation, the species evolves. Evolution is the mechanism for
which changes and improvements are made based on which particular traits generate
?fit? individuals who can survive. In other words, it is a game of ?survival of the
fittest?. This basic idea is applied to optimization problems by starting with a set of
possible solutions. Each solution is encoded into a sequence of genes (hence the name
genetic algorithms), where they compete by being measured against the objective
function of the optimization problem. These individuals are then ranked by their
level of optimality, reproduce and mutate to form a second generation of solutions.
The theory is that as the generations progress, the new solutions will become better
and better over time. To help understand the genetic algorithm process and the
discussions presented in this study, some of the basic terms are listed below. These
terms were adapted from Starkey [87]:
46
? Population: A group of potential solutions.
? Generation: A particular time or iteration in the process of the algorithm.
? Encoding: The method in which real numbers are converted into binary.
? Chromosome: The compilation of parameters for a potential solution (i.e. the
encoded version of the design vector).
? Gene: The building blocks of a chromosome (i.e. the encoded version of one
design variable in the design vector).
? Objective Function (or Fitness Function): A performance indicator, taken as
the function to be maximized or minimized in the optimization problem. The
value is often called its fitness.
? Reproduction: The process in which a new generation is created.
? Selection: The procedure by which survivors are selected to reproduce and
populate the next generation.
? Parents: Survivors of the selection process who may reproduce.
? Offspring/Children: The resulting designs from mating pairs of parents.
? Mating: The mechanism through which two parent chromosomes can produce
offspring by the exchange of information.
? Crossover: The actual exchange of information between survivors.
47
? Mutation: A random process whereby a chromosome undergoes a small change
(such as a flipping of a binary bit from 0 to 1 or vice versa).
Now with a better understanding of the terminology used in the field of genetic
algorithms, the actual algorithm or process of solving an optimization problem using
genetic algorithms can be discussed. Figure 4.1 presents a flowchart illustrating
the basic procedures found in genetic algorithms. This process is generally the
steps taken and form a base for such algorithms. Many algorithms have made
tweaks, updates, and changes to this procedure to make it more accurate, quick,
and effective.
The general genetic algorithm optimization process starts with the creation
of an initial population of designs. This population can be generated randomly or
?intelligently? using statistical methods (e.g., Latin hypercube sampling). Next,
each individual in the (initial) population is evaluated using the objective function
to quantify its fitness. Then, the convergence criteria is checked. Usually this checks
to see if a desired optimum value is met, or if the progression has ceased to improve
from one generation to the next. If the population is still the initial population, it is
likely that the convergence criteria hasn?t been satisfied. If not, then the population
follows the mating process.
The mating process includes a variety of steps. First, the individuals in the
population are ranked according to their fitness. Then, the numerous pairs of de-
signs, known as parents, are selected. The selection process varies from algorithm
to algorithm. Once the parent pairs are decided upon, they follow the crossover
48
Begin
Generate initial
population of size N
Calculate fitness of each
individual in population
Converged?
End
Mate:
create offspring
Generation +1
Start
Rank
population
Selection
of parents
Crossover
Mutation
Stop
yes
no
Figure 4.1: Generalized genetic algorithm process.
49
and mutation processes. The parent design vectors are encoded into binary chro-
mosomes, where the information is split and spliced forming new designs, called
offspring. A fraction of these offspring are subject to random binary bit flips in
a process called mutation. This completes the reproduction process through mat-
ing; the next generation of designs has been formed. This new generation then
repeats the process, starting with calculating the fitness of the individuals in the
new population. This process is repeated until the convergence criteria are met.
The overarching premise is that the solutions with the best fitness are more
likely to survive the selection process and mate to produce the next generation
of designs. Over time, the new generations build off the success of the survival
of the previous generations, thus ?evolving? into the best fit (or optimum) designs.
Genetic algorithms are beneficial in many optimization problems due to their ability
to solve the problem without the need of gradient information. It is a zero-th order
method, meaning only the objective function value is needed. Also, the mutation
and crossover processes give more assurance of reaching a global optimum, where
gradient-based methods can only guarantee finding a local optimum.
4.1.2 Multi-Objective Genetic Algorithms (MOGA)
This section serves as a brief introduction into the additions necessary to solve
multi-objective optimization problems (MOOP) using genetic algorithms. It will
provide the basic definitions of the terms used in this study as well as a top level
explanation into the process of solving an optimization problem using a standard
50
genetic algorithm discussed in Section 4.1.1. There are many different schemes,
additions, methods, and algorithms in the literature aimed at solving multi-objective
problems using genetic algorithms. Yet, for the most part, they stick to the overall
premise of genetic algorithms. This section will base its discussion off of the NSGA-II
algorithm developed by Deb et al. [93]. A detailed discussion on MOGA algorithms
and nuances to the processes involved are available in a text written by Deb [94].
To help understand the multi-objective genetic algorithm process and the dis-
cussions presented in this study, some of the basic terms used in the field of multi-
objective genetic algorithms are listed below. These definitions were adapted from
Deb [94]:
? Domination: A design x1 is said to dominate another design x2 if both condi-
tions are true:
1. The design x1 is no worse than x2 in all objectives.
2. The design x1 is strictly better than x2 in at least one objective.
Note, that if x1 dominates x2 it is also customary to say that x1 is non-
dominated by x2.
? Non-dominated set: Among a set of designs P , the non-dominated set of
designs P ? are those that are not dominated by any member of the set P .
? Global Pareto-optimal set/front: The non-dominated set of the entire feasible
search space is the global Pareto-optimal set. Each design within the set is
said to be globally Pareto-optimal.
51
? Local Pareto-optimal set/front: The non-dominated set of a local feasible
search space is the local Pareto-optimal set. Each design within the set is said
to be locally Pareto-optimal.
? Extreme Pareto-optimal points: The points in the Pareto-optimal set which
have the best value for each individual objective function (See Figure 6.1 for
an example).
? Rank: An attribute of each design, given to quantify its level of Pareto-
optimality. Among set of solutions P , the designs in the local Pareto-optimal
set are said to have a rank of 1. The designs in the local Pareto-optimal set
that exists if all rank 1 designs are removed from P , are said to have a rank
of 2. This ranking process is repeated until all designs in P are ranked.
? Diversity: A quality metric for a set of Pareto-optimal designs. The goal is
to have a set of designs with maximum diversity. Two designs are diverse in
objective space if their Euclidean distance in the objective space is large. For
a set of results, this means the desired effect is to have large, even spacing be-
tween individuals in the Pareto-optimal set. Maximum diversity helps ensure
an equal and comprehensive sampling of the trade-offs in the design space.
? Spread: A quality metric for a set of Pareto-optimal designs. The goal is to
maximize spread. The spread is the total Euclidean distance in objective space
between the extreme Pareto-optimal points. A maximum spread helps ensure
a wide sampling of the trade-offs in the design space.
52
Now with a better understanding of the terminology used in multi-objective
genetic algorithms, the process of solving a multi-objective optimization problem
using genetic algorithms can be discussed. The process is illustrated by discussing
changes or additions made to the general genetic algorithm procedure given in Figure
4.1. Once again, the process discussed here is generalized and forms a base for
MOGAs. Many different algorithms exist which may make tweaks, updates, and
changes to this procedure to make it more accurate, quick, and effective.
The first difference in following the generalized genetic algorithm procedure for
MOOPs is in the convergence decision phase. Here, different criteria are necessary
to see if a set of solutions has sufficiently progressed towards an optimum, rather
than just the single most optimum design out of a given population. Many different
criteria exist which tie in the metrics such as diversity and spread, among others.
The second difference is in the Rank population step in the flowchart. Note, do
not confuse the word Rank in Rank population with the definition of rank in the
above list of terms. In this step, since there are multiple objectives, ranking each
of individuals from best to worst is not straight forward. Most often than not, an
aggregation-based scheme is taken to add up the fitness between the objectives,
while scaling the values by taking into account a design?s rank and diversity. Most
of the other procedures within the flowchart remain the same for multi-objective
problems.
53
4.2 Robust Optimization Methods
Consider the following optimization problem:
min
~x
f (~x, p)
~x ? X = {~x|~xlb ? ~x ? ~xub} ,
(4.1)
where ~x is a vector of design variables and p is an operating condition or an envi-
ronmental parameter. The problem defined by Equation 4.1 will be referred to as
the original optimization problem. The original optimization problem can contain
different kinds of uncertainties. Uncertainties can exist in the design vectors (~x),
the environmental parameter (p), or in the objective function output (f (~x, p)) [4].
This work deals with solving the original problem while accounting for changing
environmental parameters or operating conditions due to uncertainty. For a more
detailed review on robust optimization and the current state-of-the-art methods, a
comprehensive survey is provided by Beyer and Sendhoff [4].
Generally, robust optimization methods incorporate uncertainty by construct-
ing robust counterpart objective functions F to the original objective function f [4].
For example, worst-case scenario methods use a minimax-type formulation in or-
der to construct a robust counterpart function. Aggregation methods form robust
counterparts based on the mathematical moments of f over the uncertainty. For
the methods described in Section 4.2.1, set f = FI . For Section 4.2.2, set f = FII
or FIII . For the method in Section 4.2.3, the robust counterpart formulation is
multi-objective, so the original objective function f is replaced by a set of M robust
counterpart objective functions, f = [FIV1 , ..., FIVM ].
54
In the present work, a novel multi-objective robust counterpart formulation
is proposed, where post-optimality data handling techniques can identify, among a
set of Pareto-optimal solutions, the optimal worst-case scenario and expected value
solutions for any uncertainty PDF. Solving the robust multi-objective formulation
once produces a set of designs, thus giving the designer more information about the
design space. The foundations behind a robust regularization-based method, two
aggregation-based methods, and the novel multi-objective based robust optimiza-
tion method are presented and compared in this Sections 4.2.1, 4.2.2, and 4.2.3,
respectively. The optimization problem formulation, convergence criteria, and data
handling techniques behind the robust multi-objective optimization method are dis-
cussed in detail in Section 4.2.3.
The robust multi-objective optimization method is applied to an example op-
timization problem in Section 4.3. The robust optimum results found using the
data handling techniques are compared to the true robust optimum results found
analytically. Accuracy and convergence studies were conducted by running multi-
ple simulations with varying method initial conditions and data handling technique
settings. The method accuracy was measured by calculating the average error and
root-mean-square error between the predicted and analytical robust optimum re-
sults.
The robust multi-objective formulation and three commonly used robust method
formulations were applied to the ABLT design problem (Section 5.3.2) and the re-
sults from each method were compared (Section 6.2.2). It was found that the robust
multi-objective optimization method was able to find the same results as the other
55
robust optimization methods, while providing more insight into the physics behind
the design problem by illustrating trade-offs that exist in the design space.
4.2.1 Robust Regularization
This form of accounting for uncertainty is directly analogous to the minimax
approximation in the field of applied mathematics [69]. This method appears in
research spanning multiple disciplines and may be referred to as minimax, robust
regularization [7], or simply as a robust counterpart approach [70].
A robust regularization-based method accounts for uncertainty by formulat-
ing a robust counterpart function that represents the worst-case scenario. For a
minimization problem, as in Equation 4.1, a robust regularization-based method
minimizes the maximum f value over all possible uncertainty in parameter p. The
robust counterpart function is defined as the following:
FI (~x) = sup
p?P
f (~x, p) , (4.2)
where P is the full range of possible operating conditions due to the uncertainty
in parameter p. The set P is defined by the upper and lower bounds, pub and plb
respectively:
p ? P = {p|plb ? p ? pub} . (4.3)
Equation 4.3 defines the parameter domain for nominal value of p plus any possible
associated uncertainty. Given this definition, the robust regularization technique
is nominally used when uncertainty is modeled deterministically. Optimization
problems where the uncertainty is modeled probabilistically, probability measures
56
describing the likelihood for a particular value of p in that domain to occur are pro-
vided [4]. Probabilistically modeled uncertainty problems can also be solved with
this method, as long as an upper and lower bound on p can be defined. Inserting
Equation 4.2 and 4.3 into the original optimization problem (Equation 4.1) forms the
robust optimization problem defined by using a robust regularization-based method:
min
~x
sup
p
f (~x, p)
~x ? X = {~x|~xlb ? ~x ? ~xub}
p ? P = {p|plb ? p ? pub} .
(4.4)
4.2.2 Aggregation Approach
The robust regularization approach is considered as a worst-case scenario and
may be too conservative for some applications. In such applications, a probabilistic
aggregation approach to forming a robust counterpart function is often used. The
most common technique is to use an approximation of the conditional expectation
of f [5, 6, 8, 9, 71?75]. The simplest way of calculating the expected value of the
original objective function for a particular distribution of uncertainty is through
explicit averaging [6,8,9]. Given a nominal value p and some associated uncertainty,
the robust counterpart function using explicit averaging is defined by,
FII (~x, p) =
n?
i=1
1
n
f (~x, p+ ?pi) , (4.5)
where ?pi is a small perturbation generated from the realization of the probabil-
ity density function of the uncertainty in the nominal value of the environmental
parameter p.
57
In the context of applying evolutionary or genetic algorithms, Tsutsui et al.
[72,73] developed a more efficient variation of this method. The method was reduced
from taking an average fitness over n samples to adding a stochastic perturbation
to each design in the population for a given generation. Much of the previous work
using this perturbation method covered uncertainty in the design variables, while
a few others [5, 74] have found applications using other types of uncertainties. The
current work will apply the perturbation method to uncertainty in the environmental
parameter p, thus forming the following robust counterpart function:
FIII (~x, p (t)) = f (~x, p+ ?p (t)) , (4.6)
where t is the generation index and ?p (t) is the random perturbation added to the
nominal value of the parameter at generation t. The distribution of the random
perturbation should be consistent with the probability density function of the un-
certainty associated with the environmental parameter. For an infinite population
size, it has been shown that the perturbation method is equivalent to the explicit
averaging method [4, 8, 74, 75]. Also note that if n = 1 and ?p is chosen randomly
per generation, then the explicit averaging method is equivalent to the perturbation
method (i.e., Equation 4.5 reduces to Equation 4.6).
4.2.3 Novel Multi-Objective Approach
The explicit averaging and perturbation methods described in Section 4.2.2
form an optimization problem aimed at minimizing the expected value of the original
objective function over the uncertainty in the environmental parameter. Optimizing
58
the expected value may lead to a design with a large variance; this is not always
desirable [4]. The current state-of-the-art robust optimization methods begin to
address this problem by forming a multi-objective problem using the variance and
expected value as conflicting objective functions [4, 9, 74]. This is a beginning step
to providing the designer with more insight into the nature of the robust design
space. Yet, these methods are still limited to one particular uncertainty PDF for
one nominal value.
The current work developed a novel robust multi-objective optimization for-
mulation and method which allows for post-optimality data handling techniques to
find robust optimum designs that follow the worst-case scenario and best mean fit-
ness approaches of robust regularization and explicit averaging. The method has
the ability to post-optimally pick designs that are most robust for any uncertainty
PDF and for various nominal values or the uncertain parameter. Therefore, many
types of results can be identified and analyzed by solving only one multi-objective
optimization problem.
The robust multi-objective formulation should be solved using a multi-objective
optimization solver, algorithm, or method which is most appropriate for the applied
application. For example, if the objective function in the original optimization
problem (Equation 4.1) is expensive, a gradient-based method using meta-models
or surrogates may be most appropriate. If the objective function contains many lo-
cal minima or maxima, and the objective function is relatively inexpensive, a genetic
algorithm-based solver would be more appropriate. If the expected value integral is
inexpensive and accurate to approximate using Monte Carlo integration, it may be
59
most efficient to couple the aggregation technique with a gradient-based method.
Keep in mind the pros/cons of the methods reviewed in Section 1.4.3.1 should be
consulted to evaluate which optimization formulation and solver (gradient-based
with explicit averaging, perturbation-based with genetic algorithms, meta-models,
etc.) would be best.
4.2.3.1 Robust MOOP Formulation
For the original optimization problem in Equation 4.1, the robust multi-
objective optimization method forms a set of robust counterpart functions given
by,
min
~x
Fm (~x) = E [f |~x, pm] m = 1, 2, . . . ,M
~x ? X
X = {~x|~xlb ? ~x ? ~xub}
(pm + ?m) ? Pm
?m ? U (?lbm , ?ubm)
P =
M?
m=1
Pm ,
(4.7)
where P is the full range of possible operating conditions due to the uncertainty
in parameter p. The set P is defined by the upper and lower bounds, pub and plb
respectively:
p ? P = {p|plb ? p ? pub} . (4.8)
The method assigns each of the M objective functions with the task of minimizing
the expected value of f over an uncertainty subset Pm. Each subset has a defined
nominal value pm and associated uncertainty ?m, which is sampled uniformly, such
that (pm + ?m) always exists within Pm. Note that this removes the functional
60
dependence on the nominal value pm, thus Fm is solely a function of ~x. The number
of parameter uncertainty subsets, M , and the size of each subset Pm is chosen by the
designer. Once this multi-objective optimization problem is solved, post-optimality
data handling equations can be used to predict different types of robust optimum
designs.
As mentioned in the introduction to this chapter, Equation 4.7 can be solved
using any multi-objective optimization method (direct search methods, surrogate
assisted methods, genetic algorithms, gradient-based methods, etc.). From this
point on, all investigations and applications of the novel robust MOOP formulation
and data handling method in this work are carried out from a genetic algorithm
perspective. In particular, the perturbation approach by Tsutsui [72, 73] is used
to handle the expected value operators, E [f |~x, pm], in Equation 4.7. If another
method is taken to solve the MOOP, use the data handling analysis in Section
4.2.3.2 starting with Equation 4.11, replacing ?FIVm (~x)?h with Fm (~x) when using
the post-optimality equations (Equations 4.11, 4.12, 4.13, 4.14, and 4.15).
The robust MOGA equivalent of Equation 4.7, using the perturbation method,
is given by,
min
~x
FIVm (~x, p?m (t)) = f (~x, p?m (t)) m = 1, 2, . . . ,M
~x ? X
X = {~x|~xlb ? ~x ? ~xub}
p?m (t) ? U (plbm , pubm)
Pm = {p|plbm ? p ? pubm}
P =
M?
m=1
Pm .
(4.9)
61
Parameter t is the genetic algorithm generation index, and p?m (t) is a random pa-
rameter value chosen over a uniform distribution of points belonging to Pm. The
number of parameter uncertainty subsets, M , and the size of each subset Pm is
chosen by the designer. The versatility in predicting results for different proba-
bility distributions increases with M . This will become apparent when discussing
the post-optimality data handling techniques. Yet, as M increases, multi-objective
genetic algorithms have a more difficult time advancing the Pareto-optimal front
and producing diverse results [94, 95]. Once again this formulation can viewed as a
multi-objective implementation of the perturbation method by Tsutsui et al. [72,73],
as applied to parameter uncertainties.
4.2.3.2 Convergence Criteria and Data Handling
The robust multi-objective optimization method presented here utilizes a per-
turbation-based approach. The robust counterpart environmental parameter p?m (t)
changes from generation to generation, rendering generation history plots based on
objective function values noisy and difficult to interpret. Therefore, this method can-
not use convergence criteria based on information from a single generation. Conver-
gence criteria of this type are typical in non-perturbation-based optimization [94,96].
Instead, the Pareto-optimal solution set and the convergence status of the method
at a given generation needs to be an average over previous generations.
At any generation index t, there exists a local set of Pareto-optimal solutions.
These points are the non-dominated designs in the particular objective space defined
62
by the parameter set [p?1 (t) , ..., p?M (t)]. Qualitatively, the set of true Pareto-optimal
solutions, x?, are the ones that are most common for all possible sets of p?. An
approximated set of true Pareto-optimal solutions, x?, at generation t is calculated
by the following steps:
1. For each generation from 1 to t, identify the local Pareto-optimal solutions
and add them to a running archive.
2. Generate a v-dimensional histogram with bins containing counts of the design
vectors in the archive, where v is the length of ~x.
3. Isolate the bins with bin counts greater than the mean bin count over all
bins. The mean values of x1, ..., xv, in each of these v-dimensional bins, form
a design.
4. Collectively, the designs from Step 3 form the approximation of the set of true
Pareto-optimal solutions, x?t.
As the generations progress, changes in the set of designs x?t should become less fre-
quent and the set should eventually converge. No rigorous proof of convergence has
been performed, but convergence towards a particular x? was demonstrated for the
example optimization problem in Section 4.3 and the ABLT optimization problem
results in Section 6.2.2.4. Accuracy and convergence trends corresponding to this
method are investigated in Section 4.4. Since x?t is an approximation, the design
vectors in this set are not guaranteed to be non-dominated. Therefore, the final step
after convergence is to filter out all the dominated points to ensure that the final
63
set is completely Pareto-optimal.
Changes in x?t from one generation to the next were measured as the total
number of unique designs that were added to or subtracted from x? from one gen-
eration to the next. This metric, represented as a percent change with respect to
the average size of x? over generation t and t ? 1, is referred to as ?PM (t). For
example, given a single design variable (v = 1), if the approximated set of true
Pareto-optimal solutions at generation 100 is x?100 = {1, 2, 3, 4} and at generation
101 is x?101 = {1, 2, 3, 4, 5, 6}, then the value of ?PM (101) = 40%. Written in set
notation, ?PM (t) is defined by,
?PM (t) = 100 ?
|x?t4 x?t?1|
1
2 (|x?t|+ |x?t?1|)
. (4.10)
The results from the perturbation-based robust multi-objective formulation are in
the form of a generation-dependent Pareto-optimal hypersurface in objective func-
tion space. The Pareto-optimal solutions most common in a large generation by
generation archival set are the results of interest.
Once an approximation to the true set of Pareto-optimal solutions is identified,
post-optimality data handling techniques can be employed to identify different types
of ?preferred? robust optimal designs [94]. An optimal worst-case scenario solution is
found from the set of Pareto-optimal solutions using the following weighted minimax
formula:
~xwcs = min
~x?x?t
(
M
max
m=1
wm?FIVm (~x)?h
)
, (4.11)
where, ?FIVm (~x)?h is the mean m-th objective function value over h generations (1
to t) for the Pareto-optimal design ~x ? x?t. The weights w1 through wm are defined
64
by the following:
wm =
?
????
????
0, if Pm ? Pwcs
1, if Pm 6? Pwcs
,
Pwcs ? P ,
(4.12)
where Pwcs is the uncertainty subset over which the worst-case scenario is being
analyzed. In other words, the designer may not always want to find the worst-case
scenario over the entire uncertainty set P, but yet a subset Pwcs.
A minimized expected value solution is found by using a weighted L1-norm
formula defined by,
~xev = min
~x?x?t
(
M?
m=1
wm|?FIVm (~x)?h ? z
?
m|
)
, (4.13)
where z?m is the ideal m-th objective function value, equal to the best possible value
found over the entire Pareto-optimal set defined by:
z?m = min
~x?x?t
?FIVm (~x)?h . (4.14)
The weights w1 through wM are set by the nominal parameter value, p, and accom-
panying perturbation ? (a realization from the uncertainty PDF). Typically, these
weights are defined such that they sum to one [94]. The weights can be calculated
as follows:
wm =
Pr (p+ ? ? Pm)
Pr (p+ ? ? P)
. (4.15)
For example, consider an optimization problem aimed to minimize the ex-
pected value of f (~x, p) where p = 0.5 and the uncertainty in p is assumed to be uni-
formly distributed about p with a lower and upper bound of plb = 0.0 and pub = 1.0.
65
Assume the robust MOGA method problem is using five objective functions (M = 5)
where the uncertainty is split into equally spaced subsets Pm :
P1 = {p|0.0 ? p < 0.2} ,
P2 = {p|0.2 ? p < 0.4} ,
P3 = {p|0.4 ? p < 0.6} ,
P4 = {p|0.6 ? p < 0.8} ,
P5 = {p|0.8 ? p ? 1.0} .
(4.16)
Following Equation 4.15 the calculation of wm would produce
w1 = w2 = w3 = w4 = w5 =
1
5
. (4.17)
Now consider the same problem for uncertainty that follows a truncated normal
distribution rather than a uniform distribution. Here p values are sampled from a
truncated normally-distributed probability density function with a mean of 0.5,
standard deviation of 0.2, and bounded by plb = 0.0 and pub = 1.0. Equation 4.15
would then produce the following weights:
w1 = 0.0614 ,
w2 = 0.2449 ,
w3 = 0.3879 ,
w4 = 0.2449 ,
w5 = 0.0614 .
(4.18)
66
4.3 Illustrative Example
An example optimization problem will be used in this section to illustrate the
use of the robust MOGA formulation proposed in Section 4.2.3. An optimization
problem with two design variables and one environmental parameter will be used as
the original problem. The original optimization problem formulation is presented
in Section 4.3.1. Worst-case scenario robust optimum solutions and minimum ex-
pected value robust optimum solutions were derived analytically in Section 4.3.2.
The multi-objective robust counterpart optimization problem was then formed and
solved using genetic algorithms in Section 4.3.3. The optimization results and the
overall effectiveness of the proposed method was analyzed by comparisons with the
analytical solutions. The optimization problem formulation and set of results in
Section 4.3.3 will be referred to as the baseline simulation.
4.3.1 Original Optimization Problem
Consider the following original optimization problem:
min
~x
fi.e. (~x, p) = 1?
(
0.4 +
1
2 + x1
)
exp
(
?
(p? x2)
2
2x12
)
~x ? X
X = {(x1, x2) |x1 ? [0.1, 1] , x2 ? [0, 1]} ,
(4.19)
where the environmental parameter p ranges from zero to one.
Figures 4.2 and 4.3 present two-dimensional slices of the design space defined
by Equation 4.19. The objective function value fi.e. as a function of p for increasing
values of x1 and a constant x2 forms inverted bell shaped curves with increasing
67
0 0.2 0.4 0.6 0.8 10.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
p
f i.e.
(x,p
)
x1 = 0.1
x1 = 0.2
x1 = 0.4
x1 = 0.6
x1 = 0.8
x1 = 1
~x
(a)
0 0.2 0.4 0.6 0.8 10.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
p
f i.e.
(x,p
)
 
 
x2 = 0.0x2 = 0.2x2 = 0.4x2 = 0.6
~x
(b)
Figure 4.2: fi.e. curves as a function of p and ~x.
68
0 0.2 0.4 0.6 0.8 10.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x1
f i.e.
(x,p
)
p = 0
p = 0.05
p = 0.1
p = 0.3
p = 0.5
p = 1
~x
(a)
0 0.2 0.4 0.6 0.8 10.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x1
f i.e.
(x,p
)
p = 0.5
p = 0.55
p = 0.6
p = 0.8
p = 1
~x
(b)
Figure 4.3: fi.e. curves as a function of p and ~x.
69
width (Figure 4.2(a), where x2 = 0.0). As x1 increases from 0.1 to 1.0, the width of
the curve increases and the minimum value of fi.e. increases slightly, thus flattening
and raising the curve. The values of fi.e. as a function of p and x2 for a constant x1
form an inverted bell shape of constant width, with x2 defining the center of the bell.
Contours of fi.e. vs. p for varying values of x2 with x1 = 0.1 are plotted in Figure
4.2(b). The objective function when plotted as a function of x1 for given values of
p and a constant x2 is presented in Figures 4.3(a) and 4.3(b), where x2 = 0.0 and
x2 = 0.5, respectively. This particular objective function was chosen due its smooth
monotonic nature and its ability to render a range of robust optimum solutions
spanning the design space.
4.3.2 Analytical Robust Optimum Solutions
4.3.2.1 Optimal Expected Value
A robust counterpart approach to Equation 4.19, where the goal is to minimize
the expected value of fi.e., forms the following robust optimization problem:
min
~x?X
E [fi.e.|~x, p] . (4.20)
Inserting fi.e. from Equation 4.19, E [fi.e.|~x, p] becomes
??
??
[
?
(
0.4 +
1
x1 + 2
)
exp
(
?
(p+ ? ? x2)
2
2x12
)
+ 1
]
? (?) d? , (4.21)
where ? is a perturbation away from a nominal value, p, due to uncertainty and ? (?)
is the accompanying probability density function of the uncertainty. After applying
the assumption that for any value of p the perturbation ? is a realization from a
70
zero mean normally distributed PDF with a variance of ?2,
? (?) = N
(
0, ?2
)
, (4.22)
the expression in Equation 4.21 simplifies to
?
0.4 + 1x1+2?
?2
x21
+ 1
exp
(
? (p? x2)
2
2 (?2 + x21)
)
+ 1 . (4.23)
The global robust optimum design ~x? that minimizes Equation 4.23 with respect to
~x is
~x? = (x?1, x
?
2) ,
?2 =
5x?1
3
2x?1
2 + 8x?1 + 18
,
x?2 = p .
(4.24)
4.3.2.2 Optimal Worst-Case Scenario
Using a robust regularization approach, where the uncertainty in p is defined
over some uniformly distributed neighborhood  ? [a, b], the robust optimum coun-
terpart to Equation 4.19 would be,
min
~x?X
sup
?[a,b]
fi.e. (~x, ) , (4.25)
where 0.0 ? a ? b and b ? 1.0. Solving the inner optimization problem, i.e. finding
the supremum, gives the following:
sup
?[a,b]
fi.e. (~x, ) =
?
????
????
A, if B or C
D, if E or F
, (4.26)
71
where,
A = ?
(
0.4 +
1
x1 + 2
)
exp
(
?
(b? x2)
2
2x12
)
+ 1 ,
B = x2 ? a ,
C = a ? x2 ? b and x2 ?
(
b? a
a
+ a
)
,
D = ?
(
0.4 +
1
x1 + 2
)
exp
(
?
(a? x2)
2
2x12
)
+ 1 ,
E = x2 ? b ,
F = a ? x2 ? b and x2 ?
(
b? a
a
+ a
)
.
(4.27)
The global robust optimum design ~x? that minimizes Equation 4.26 with respect to
~x (i.e. solving the outer optimization problem in Equation 4.25) is
~x? = (x?1, x
?
2) ,
(
b?a
2
)2
=
5x?1
3
2x?1
2 + 13x?1 + 18
, (4.28)
x?2 =
(
b
a
+ a? 1
)
.
4.3.3 Robust MOGA Formulation & Data Handling Solutions: Base-
line Simulation
The first step in forming the robust optimization problem using the proposed
multi-objective method is to choose the subsets Pm for the problem. For this ex-
ample, lower and upper bounds on the possible range of values of p will be 0.0 and
1.0, respectively. Therefore, using Equation 4.8, P = {p|0.0 ? p ? 1.0}. The full
72
uncertainty set was split into the following five subsets:
P1 = {p|0.0 ? p ? 0.2} ,
P2 = {p|0.2 ? p ? 0.4} ,
P3 = {p|0.4 ? p ? 0.6} ,
P4 = {p|0.6 ? p ? 0.8} ,
P5 = {p|0.8 ? p ? 1.0} .
(4.29)
These subsets were chosen in order to equally sample the uncertainty space and
balance the trade-off between accuracy of the post-optimality handling predictions
and the MOGA problem complexity (as the number of subsets increases the problem
becomes more difficult to solve [94, 95]). Following the method outlined by Equa-
tion 4.9 and given the original optimization problem in Equation 4.19, the robust
multi-objective optimization method forms the following MOGA-based optimization
problem:
min
~x
FIVm (~x, p?m (t)) = fi.e. (~x, p?m (t)) m = 1, 2, 3, 4, 5
~x ? X
X = {(x1, x2) |x1 ? [0.1, 1] , x2 ? [0, 1]}
p?1 (t) ? U (0.0, 0.2)
p?2 (t) ? U (0.2, 0.4)
p?3 (t) ? U (0.4, 0.6)
p?4 (t) ? U (0.6, 0.8)
p?5 (t) ? U (0.8, 1.0) .
(4.30)
A qualitative way of reading the optimization problem outlined in Equation 4.30 is
as follows:
? Objective 1 (min
~x
FIV1) will find the value of ~x such that the mean value of fi.e.
over p = [0.0, 0.2] is minimized.
73
? Objective 2 (min
~x
FIV2) will find the value of ~x such that the mean value of fi.e.
over p = [0.2, 0.4] is minimized.
? Objective 3 (min
~x
FIV3) will find the value of ~x such that the mean value of fi.e.
over p = [0.4, 0.6] is minimized.
? Objective 4 (min
~x
FIV4) will find the value of ~x such that the mean value of fi.e.
over p = [0.6, 0.8] is minimized.
? Objective 5 (min
~x
FIV5) will find the value of ~x such that the mean value of fi.e.
over p = [0.8, 1.0] is minimized.
With the true robust optimum designs known analytically (Section 4.3.2),
comparisons can be made to the results predicted by the robust MOGA formulation
when solving Equation 4.30 using a MOGA solver and applying the post-optimality
data handling techniques. The multi-objective optimization problem was solved
using a modified version of the MOGA solver in the MATLAB Genetic Algorithm
and Direct Search toolbox (GADS). The multi-objective genetic algorithm solver
in the GADS toolbox employs a version of an elitist non-dominated sorting genetic
algorithm method, NSGA-II [93]. The genetic algorithm options and settings used
to solve Equation 4.30 are listed in Table 4.1.
The simulation was run for 500 generations in order to illustrate the genera-
tion history trends of the convergence metric ?PM . Figure 4.4 plots the moving
average of ?PM over 10 generations as a function of current generation number.
The value of the moving average of ?PM at generation 1 is about 27. As gener-
ations progress, the moving average value of ?PM decreases rapidly until about
74
Table 4.1: MOGA options used for solving Equation 4.30.
Genetic algorithm option Setting
Population size 600
Generations 500
Initial population Random within ~xlb and ~xub
Fitness scaling Based on the individual?s rank
in the sorted scores
Selection Tournament selection
Crossover Randomized
Mutation Gaussian
Crossover fraction 0.8
generation 25. After generation 25, the overall trend continues to decrease slowly
until about generation 125. From generation 125 onward, no appreciable decrease
is seen. Similar convergence trends were found when the robust multi-objective
method was applied to the vehicle design optimization problem in Section 6.2.2.4.
This simulation was said to have converged when the moving average of ?PM over
10 generations is equal to zero. This first occurs at generation 135. Therefore, the
approximation to the true Pareto-optimal solution at generation 135 (x?135) will be
used in the analysis.
The set of Pareto-optimal solutions at generation 135 contains 337 individuals
in the population. The predicted optimal worst-case scenario and expected value
solutions can be picked out of the set of 337 solutions using Equations 4.11 and 4.13,
respectively. Figure 4.5 and Figure 4.6 compare the analytical optimal expected
value designs using Equation 4.24 to the predicted optimal expected value designs
using Equation 4.13 for p values spanning from 0.2 to 0.8 in increments of 0.1
75
0 100 200 300 400 500
0
5
10
15
20
25
30
Generation (t)
Mov
ing?
PM
aver
age
over
10g
ens.
(a) Full generation history.
0 100 200 300 400 500
0
0.5
1
1.5
Generation (t)
Mov
ing?
PM
aver
age
over
10g
ens.
(b) Generation history zoomed in.
Figure 4.4: Mean ?PM generation history.
and ? values from 0.02 to 0.2 in increments of 0.02. Errors in predicted x?1 values
over the p-? combinations reveal particular areas of accuracy and inaccuracy when
compared to the analytical values. When p = 0.3, 0.5, and 0.7, predicted values of
x?1 were consistently within ?0.075 of the analytical values for all values of ?. The
predictions for p = 0.2, 0.4, 0.6, and 0.8, diverged from the analytical values when ?
dropped below a value of about 0.15, 0.08, 0.08, and 0.08, respectively. These areas
of accuracy and inaccuracy are governed by the calculated weights in Equation 4.15.
When the environmental parameter p is 0.2, 0.4, 0.6, or 0.8, for the subsets
given by Equation 4.29, the uncertainty PDF for the lower values of ? renders
wm values that mostly span two adjacent Pm subsets. Therefore, for ? values less
than about 0.1, the weights for the two adjacent subsets Pm are about equal and
Equation 4.13 selects the same solution from the Pareto-optimal set. The solutions
thus diverge from the analytical solutions for decreasing ?. As ? becomes greater
than about 0.1, the weights wm become distributed more evenly among the other
76
0 0.05 0.1 0.15 0.200.1
0.20.3
0.40.5
0.60.7
0.80.9
?
x? 1,x? 2
 
 Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(a)
0 0.05 0.1 0.15 0.200.1
0.20.3
0.40.5
0.60.7
0.80.9
?
x? 1,x? 2
 
 Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(b)
0 0.05 0.1 0.15 0.200.1
0.20.3
0.40.5
0.60.7
0.80.9
?
x? 1,x? 2
 
 Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(c)
0 0.05 0.1 0.15 0.20
0.10.2
0.30.4
0.50.6
0.70.8
0.9
?
x? 1,x? 2
 
 
Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(d)
Figure 4.5: Baseline analytical and predicted optimum expected value design vectors
for p values of 0.2, 0.3, 0.4, and 0.5, respectively.
three subsets, which in turn enables Equation 4.13 to select different solutions within
the Pareto-optimal set.
When the environmental parameter p is 0.3, 0.5, or 0.7, the uncertainty PDF
is centered within subset P2, P3, and P4, respectively. This creates wm values that
are more distributed and dynamically changing with ?, even for values lower than
about 0.1. For ? less than about 0.04, the uncertainty PDF is mostly contained
within a single subset Pm (wm is close to unity) and Equation 4.13 selects the same
77
0 0.05 0.1 0.15 0.20
0.10.2
0.30.4
0.50.6
0.70.8
0.9
?
x? 1,x? 2
 
 
Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(a)
0 0.05 0.1 0.15 0.20
0.10.2
0.30.4
0.50.6
0.70.8
0.9
?
x? 1,x? 2
 
 
Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(b)
0 0.05 0.1 0.15 0.20
0.10.2
0.30.4
0.50.6
0.70.8
0.9
?
x? 1,x? 2
 
 
Analytical x?1Predicted x?1Analytical x?2Predicted x?2
(c)
Figure 4.6: Baseline analytical and predicted optimum expected value design vectors
for p values of 0.6, 0.7, and 0.8, respectively.
solution from the Pareto-optimal set. Errors in the value of the predicted x?2 for 66
out of the 70 p-? combinations were all within ?0.03 of their respective analytical
counterparts. The four other predicted values were within ?0.08 of the analytical
values.
To produce more accurate results, the designer would have to decrease the
size of the parameter uncertainty subsets and increase the number subsets. This
would produce a more dynamic trade-off in the weighted L1-norm selection of the
78
Pareto-optimal designs. This effect is illustrated by the results in Section 4.4.3.
Figure 4.7 compares the analytical optimal worst-case scenario designs using
Equation 4.28 to the predicted optimal worst-case scenario designs using Equation
4.11 for fifteen uncertainty neighborhoods defined by,
Pab = {p|a ? p ? b} . (4.31)
The values of a and b for each of the 15 cases are listed in Table 4.2. Figure 4.7
shows that the proposed robust MOGA method with post-optimality data handling
is successful in identifying robust optimum designs based on a worst-case scenario for
various different neighborhoods of uncertainty. Predicted values of x?1 and x
?
2 follow
the same trends as their analytical counterparts. Errors in predicted x?1 values were
all within ?0.1 except for Cases 12 and 14, where the values were overpredicted by
0.112 and 0.164, respectively. Errors in predicted x?2 values were all within ?0.03
except for Cases 1 and 5, where the values were overpredicted by 0.0773 and 0.0747,
respectively. Accuracy in the predicted solutions is dependent on the number of
subsets M used in the analysis, as well as the number and diversity of the individuals
in the Pareto-optimal set.
79
1 2 3 4 5 6 7 8 9 10 11 12 13 14 150
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
[a, b] case number
x? 1,x
? 2
 
 
Analytical x?1Predicted x?1
(a)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 150
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
[a, b] case number
x? 1,x
? 2
 
 
Analytical x?2Predicted x?2
(b)
Figure 4.7: Analytical and predicted optimum worst-case scenario design vectors for
various [a, b] values.
80
Table 4.2: Worst-case scenario comparison settings in Figure 4.7.
Case number [a, b] values
1 [0.0, 0.2]
2 [0.2, 0.4]
3 [0.4, 0.6]
4 [0.6, 0.8]
5 [0.8, 1.0]
6 [0.0, 0.4]
7 [0.2, 0.6]
8 [0.4, 0.8]
9 [0.6, 1.0]
10 [0.0, 0.6]
11 [0.2, 0.8]
12 [0.4, 1.0]
13 [0.0, 0.8]
14 [0.2, 1.0]
15 [0.0, 1.0]
4.4 Accuracy and Convergence Tests
Multiple simulations were run with varying method initial conditions and data
handling technique settings in order to test the accuracy and convergence charac-
teristics of the proposed robust MOGA method. Each simulation perturbs one of
three initial conditions and settings from the baseline simulation. The three initial
conditions and settings which were subject to change are the following:
? The initial population size used to solve the robust optimization problem with
the NSGA-II genetic algorithm
81
? The number of histogram bins used in Step 2 of the data handling process
outlined in Section 4.2.3.2
? The number of objective functions and uncertainty subsets (M in Equation
4.9) used in the multi-objective optimization problem formulation
The optimum solutions predicted by each simulation were judged using three
criteria: accuracy, quality, and speed. These criteria quantitatively describe the
results and help compare the effects of changing the simulation settings. The defi-
nitions of accuracy, quality, and speed are listed below.
? Accuracy: The accuracy of a predicted solution set is measured by the mean
absolute error and root-mean-square error of the predicted optimum solutions
when compared to the analytical optimum solutions. All errors are given in
percentage of the design variable range.
? Quality: The quality of the predicted solution set is measured by the number
of Pareto-optimal points in the final Pareto-optimal front.
? Speed: The speed of the simulation is measured by how many generations
were needed to progress the simulation to a particular level of convergence.
All the simulations were run with the same convergence criterion as the baseline
case: the simulation has converged when the moving average value of ?PM over
10 generations is equal to zero.
Table 4.3 is the test matrix which contains a list of the simulations that were
run and all of the associated input settings. Simulations P1 through P8 vary the
82
size of the initial population used in the genetic algorithm, while keeping the num-
ber of histogram bins and uncertainty subsets equal to that used in the baseline
simulation. Simulations B1 through B4 vary the number of histogram bins used in
the data handling procedure, while keeping the initial population size and number
of uncertainty subsets equal to that of the baseline. Simulations S1 and S2 vary
the number of uncertainty subsets in the optimization problem formulation while
keeping the other two parameters the same as the baseline. Accuracy, quality, and
speed comparisons for these simulations are given in Sections 4.4.1, 4.4.2, and 4.4.3.
Table 4.3: Accuracy and convergence input settings test matrix.
Number of
Number of
Simulation Initial population size
histogram bins
uncertainty subsets
(M in Equation 4.9)
Baseline 600 20 5
P1 1200 20 5
P2 300 20 5
P3 150 20 5
P4 75 20 5
P5 37 20 5
P6 18 20 5
P7 9 20 5
P8 5 20 5
B1 600 60 5
B2 600 40 5
B3 600 10 5
B4 600 5 5
S1 600 20 8
S2 600 20 3
83
4.4.1 Effects of Changes in Initial Population Size
Nine simulations (P1-P8 and Baseline) were run to compare how changes in the
size of the initial population affect the speed, quality, and accuracy of the predicted
robust optimum design vectors. Graphs illustrating these effects are provided in
Figure 4.8.
Figure 4.8(a) shows that, for a population size between 37 and 1200, a simula-
tion with a larger population size converges after a smaller number of generations.
The generation number at convergence is reduced from 406 to 156 over this range.
For the three data points with a population size less than 37, the generation number
at convergence falls off quickly from 406 to 84. This trend exists due to the influence
of the population size on the NSGA-II algorithm effectiveness. As the population
increases, a smaller fraction of points in the population are non-dominated. This
gives the algorithm a more diverse set of designs to rank and sort using an elitism
scheme. More diverse designs result in a more effective genetic algorithm and quicker
convergence. If the population is below a certain threshold, most if not all the points
in the population are non-dominated. When this happens, there is no fitness advan-
tage to any of the solutions. Therefore, the mutation and crossover operators are
solely responsible for creating solutions in a better Pareto-optimal front, allowing
the algorithm to continue to make progress [94]. Consequently, changes in the set
of Pareto-optimal points are sparse. This results in premature convergence when
using the ?PM convergence criterion. This concept is discussed in more detail in
Section 8.8.2 in Deb [94].
84
Figure 4.8(b) shows that increasing the initial population size from 5 to about
37 will increase the number of Pareto-optimal points in the final set of optimum de-
signs from 53 to 283. This drastic change in quality is likely to be a consequence of
the premature convergence found in the analysis of the population effects on speed.
This trend quickly levels off for a population size between 37 and 1200 individuals
(note the x-axis is shown in log scale), where the number of Pareto-optimal points
grows from 283 to 337. This relatively mild increase is attributed to the fact that
a larger population aids in the identification of Pareto-optimal points by providing
a better probability that any particular Pareto-optimal design will be involved in
the optimization process. A smaller population has a higher likelihood that a non-
dominated candidate will be skipped over due to limited sampling. This mild but
steady trend from 283 to 337 Pareto-optimal points looks to be asymptotically pro-
gressing towards a theoretical maximum number of Pareto-optimal points possible
for that optimization problem. This theoretical maximum number is found to be
determined by the number of histogram bins.
Figure 4.8(c) shows that with increasing population size, the mean error and
RMS error of the predicted designs decrease. Note that the x-axis is shown in log
scale, highlighting that the accuracy metrics quickly improve from a population size
of 5 to about 37. This again is most likely attributed to premature convergence.
For a population size of 37 to 1200, improvement in accuracy is gradual, where the
mean error reduces from 9.7% to 8% and the RMS error reduces from 7.2% to 6%.
This reduction is caused by an increasing number of samples present in each bin
when performing Step 3 in Section 4.2.3.2 and the averaging operators in Equations
85
4.11 and 4.13. Averages over a larger sample size produce more accurate results
with less dispersion.
4.4.2 Effects of Changes in Number of Histogram Bins
Five total simulations (B1-B4 and Baseline) were run to compare how changes
in the number of histogram bins affect the speed, quality, and accuracy of the
predicted robust optimum design vectors. Graphs illustrating these changes are
provided in Figure 4.9.
Figure 4.9(a) shows that as the number of bins used in the data handling
process increases, the simulation takes more generations to converge. Increasing the
number of bins from 5 to 60 forces the simulation to need 939 generations to converge
rather than 13. This is a direct result of the effect seen in Figure 4.9(b). As the
number of histogram bins increases, the number of possible Pareto-optimal points
increases. When a high number of possible Pareto-optimal points exist, a generation
to generation analysis is more likely to find a large number of new Pareto-optimal
points, thus a large value of ?PM . The value of ?PM will need to progress through
more generations to satisfy the convergence criterion.
Figure 4.9(b) shows that increasing the number of bins used in the data han-
dling process will increase the number of Pareto-optimal points in the final set of
optimum designs. Increasing the number of bins from 5 to 60 increases the result-
ing number of Pareto-optimal points from 26 to 2214. This drastic change is a
result of the data handling technique outlined in Section 4.2.3.2. The number of
86
histogram bins in this process defines a theoretical upper limit on the number of
possible Pareto-optimal points for a given simulation. The number of histogram
bins discretizes the predicted Pareto-optimal hypersurface by setting a resolution or
fineness to each design vector in Step 3 of the data handling process.
Figure 4.9(c) shows that a particular optimum number of bins exist which
minimizes the mean error and RMS error in the predictions. If there are not enough
or too many bins, the accuracy suffers. The optimal number of histogram bins in
this study is 20. If there are too few bins (10 or 5 bins in this example), the number
of Pareto-optimal points decreases significantly (Figure 4.9(b)). As the number of
Pareto-optimal points decreases, the post-optimality data handling equations have
significantly less designs to choose from. This would result in a loss of accuracy for
each selected design. If there are too many bins (40 or 60 bins in this example), the
number of samples present in each bin when performing Step 3 in Section 4.2.3.2
and the averaging operators in Equations 4.11 and 4.13 decreases. This means the
averaging operators would have a smaller sample size and become less accurate.
This effect was also seen and discussed in Section 4.4.1.
4.4.3 Effects of Changes in Number of Uncertainty Subsets
Three total simulations (S1, S2, and Baseline) were run to compare how
changes in the number of uncertainty subsets affect the speed, quality, and accuracy
of the predicted robust optimum design vectors. Graphs illustrating these changes
are provided in Figure 4.10.
87
Figure 4.10(a) shows that as the number of uncertainty subsets used in the
optimization formulation increases, the simulation takes longer to converge. In-
creasing the number of uncertainty subsets from 3 to 8 increased the generation
number at convergence from 108 to 187. This effect is due to the same phenomena
discussed in Section 4.4.1. By the optimization problem definition in Equation 4.9,
more uncertainty subsets means the multi-objective optimization problem involves
more objective functions. The increase in the number of more objective functions
increases the fraction of points in the population which are non-dominated, thus de-
creasing the effectiveness of the NSGA-II algorithm (see Section 8.8.2 in Deb [94]).
Figure 4.10(b) shows that increasing the number of uncertainty subsets used
in the optimization formulation from 3 to 5 will substantially increase the number
of Pareto-optimal points in the final set of optimum designs from 177 points to
337 points. Adding objective functions inherently increases the number of possible
Pareto-optimal points, due to the increase in the dimensions of the Pareto-optimal
front. Increasing the number of subsets from 5 to 8 shows a minimal increase from
337 to 348 points. This is due to either the lack of individuals in the population or
a lack of histogram bins for an optimization problem with eight objective functions.
This effect is coupled with the effects seen in Sections 4.4.1 and 4.4.2.
Figure 4.10(c) shows that the accuracy (both mean and RMS errors) improves
with an increasing number of uncertainty subsets. This effect on accuracy proved
to be the strongest out of all changes seen in Section 4.4. Increasing the number
of uncertainty subsets from 3 to 8 decreased the mean error from 21.5% to 5.9%
and the RMS error from 14.6% to 3.8%. This is largely due to the fact that the
88
approximations of the uncertainty PDFs through the calculations of the weights in
Equations 4.12 and 4.15 are more accurate. This is especially true for the smaller
values of ?. More uncertainty subsets are utilized when making the predictions
and the changes in the weights are more dynamic. This concept is discussed and
illustrated in the full analysis of the baseline results in Section 4.3.3.
89
0 200 400 600 800 1000 120050
100
150
200
250
300
350
400
450
Gen
erat
ion
#a
tco
nve
rgen
ce
Initial population
 
 Baseline
(a) Speed vs. initial population size
101 102 1030
50
100
150
200
250
300
350
400
Num
ber
ofP
aret
opo
ints
Initial population
 
 Baseline
(b) Quality vs. initial population size
101 102 103
4
6
8
10
12
14
16
18
20
22
Mea
n/R
MS
erro
rs(i
n%
)
Initial population
 
 Mean errorRMS errorBaseline
(c) Accuracy vs. initial population size
Figure 4.8: Speed, quality, and accuracy changes with change in initial population
size.
90
0 10 20 30 40 50 60
0
200
400
600
800
1000
Gen
erat
ion
#a
tco
nve
rgen
ce
Number of bins 
 Baseline
(a) Speed vs. number of histogram bins
0 10 20 30 40 50 60
0
500
1000
1500
2000
2500
Num
ber
ofP
aret
opo
ints
Number of bins 
 Baseline
(b) Quality vs. number of histogram bins
0 10 20 30 40 50 60
4
6
8
10
12
14
16
18
20
22
Mea
n/R
MS
erro
rs(i
n%
)
Number of bins 
 Mean errorRMS errorBaseline
(c) Accuracy vs. number of histogram bins
Figure 4.9: Speed, quality, and accuracy changes with change in number of his-
togram bins used in the data handling procedure.
91
2 3 4 5 6 7 8 980
100
120
140
160
180
200
220
Gen
erat
ion
#a
tco
nve
rgen
ce
Number of uncertainty subsets
 
 Baseline
(a) Speed vs. number of uncertainty subsets
2 3 4 5 6 7 8 9
50
100
150
200
250
300
350
400
450
500
Num
ber
ofP
aret
opo
ints
Number of uncertainty subsets 
 Baseline
(b) Quality vs. number of uncertainty subsets
2 3 4 5 6 7 8 9
4
6
8
10
12
14
16
18
20
22
Mea
n/R
MS
erro
rs(i
n%
)
Number of uncertainty subsets 
 Mean errorRMS errorBaseline
(c) Accuracy vs. number of uncertainty subsets
Figure 4.10: Speed, quality, and accuracy changes with change in number of uncer-
tainty subsets in the multi-objective optimization problem formulation.
92
Chapter 5: Vehicle Design Optimization Methodology
Both multi-objective and robust optimization studies were performed to find
optimum vehicle geometries and to analyze trade-offs between vehicle performance
and controllability under asymmetric boundary-layer conditions. Objective func-
tions of lift-to-drag ratio, internal volumetric efficiency, maximum heat flux, and
vehicle controllability were included in the studies. The multi-objective studies
produce sets of Pareto-optimal solutions. These solutions illustrate the trade-offs
between the competing objective functions. The L2-norm optimum Pareto-optimal
points were found, identifying single optimum designs of minimum compromise [97].
The robust optimization studies focused on the objective function pertaining to ve-
hicle controllability. Several commonly used robust optimization techniques were
applied in order to identify vehicle designs which minimize the worst-case scenario
and expected value of the influence of roughness-induced ABLT. This vehicle opti-
mization problem was also used as a test-bed for the the novel robust optimization
method presented in Section 4.2.3. All optimization problems were solved using a
multi-objective genetic algorithm to predict global optimality.
Chapter 5 is organized as follows. Section 5.1 will discuss the design point
specifications for all the optimization problems investigated in this study. This
93
discussion includes design vector formulation, the side constraints applied to the
design variables, and pertinent center of gravity (CG) assumptions. Within the
multi-objective context, two sets of optimization problems were solved. One was
a case where the objective functions were defined and solved without the presence
of an isolated roughness; no asymmetric boundary-layer was induced. The second
problem involved finding optimized geometries which are most robust to a single
disturbance which could occur on the top surface of the vehicle at any spanwise
location. In other words, the resulting optimum shapes are vehicles which are the
least affected by a random asymmetric boundary-layer condition.
The two multi-objective based optimization problems, with and without an
ABLT event, will be formed in Section 5.2 and Section 5.3.1, respectively. In ad-
dition, each scenario was solved for both infinitely sharp leading edge and blunted
leading edge waveriders. The concept of a roughness-induced ABLT event occurring
on the vehicle with some probability is extended in the robust optimization studies
in Section 5.3.2. The objective functions and constraint functions for each class of
optimization problems are presented in their respective sections. Finally, Section
5.4 concludes Chapter 5 with information regarding the optimization method and
parameters used to solve each problem.
5.1 Optimization Design Point, Constraints, and Assumptions
All the optimization problems were solved at a particular set of design point
specifications. These specifications define the freestream conditions and geometric
94
constraints of the problem. The design point parameters are listed in Table 5.1.
Most of these specifications match those of Starkey [30] to aid in the validation of
the models and comparisons of the results. As a side note, the simulations were also
run at a higher Mach number (M? = 18) where it was found that the geometric
trends of the results stayed the same. The differences existed only in the magnitudes
of the objective functions. However, the geometric constraints did prove to make
significant changes to the shape of the optimum waverider designs. Studies on this
effect is suggested for future work.
Note that CG assumption I was used for all multi-objective formulations in
Sections 5.2 and 5.3.1. This was initially used as a first order assumption. Once
simulation results were collected and analyses were conducted, the decision was made
to update the CG assumption to a center of volume assumption (CG assumption
II). CG assumption II was applied to the formulations for the robust optimization
problems in Section 5.3.2.
Table 5.1: Design point specifications.
Parameter Value
Freestream Mach number 8
Altitude 60,000 ft
Freestream pressure 150 lb/ft2
Freestream temperature 389.99 ?R
Minimum vehicle internal volume 88,286 ft3
Maximum vehicle internal volume 105,944 ft3
CG (x,y,z) assumption I
(
L
2 , 0,
L
4 tan (?local)
)
CG assumption II Center of volume
95
Using the methodology outlined in Section 2.1, the design vector contains
variables fully describing a vehicle geometry. A design vector is defined by
~x =
[
n m ? ? L w
]
, (5.1)
for sharp leading edge waveriders and by
~x =
[
n m ? ? L w d
]
, (5.2)
for blunted leading edge waveriders. These vector definitions are used for all opti-
mization formulations in Chapter 5. The design vectors defined by Equations 5.1
and 5.2 belong to the design vector set space X, bounded by user defined lower
and upper bounds. For the multi-objective problems in Sections 5.2 and 5.3.1, the
bounds of the design vector were set to be
~xlb ? ~x ? ~xub =
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
0.001 ? n ? 1.0
0.001 ? m ? 1.0
1? ? ? ? 10?
1? ? ? ? 20?
1 ft ? L ? 220 ft
1 ft ? w ? 120 ft
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
, (5.3)
96
for sharp waveriders and
~xlb ? ~x ? ~xub =
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
0.001 ? n ? 1.0
0.001 ? m ? 1.0
1? ? ? ? 10?
1? ? ? ? 20?
1 ft ? L ? 220 ft
1 ft ? w ? 120 ft
0.0164 ft ? d ? 1.64 ft
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
, (5.4)
for blunt waveriders. Due to floating-point precision limits on the calculations of the
CG location for CG assumption II, the lower limit of variables n and m for Section
5.3.2 was changed from 0.001 to 0.05.
5.2 Baseline Multi-Objective Optimization
The baseline multi-objective optimization problem was formulated as follows:
min
~x
fq,scaled (~x) q = 1, 2, ..., Q
s.t.
gj (~x) ? 0 j = 1, 2, ..., J
~x ? X
X = {~x|~xlb ? ~x ? ~xub} ,
(5.5)
97
where Q is the total number of objective functions and J is the total number of
constraint functions. The objective functions for the baseline case are defined by
the following:
f1 (~x) = ?
Lift
Drag
, (5.6)
f2 (~x) = ?? , (5.7)
f3 (~x) = ?~CMCG,tot? , (5.8)
f4 (~x) = max (qw) . (5.9)
Functions f1 through f3 are used for the sharp leading edge waverider case (Q = 3)
and f4 is added for the blunt waverider case (Q = 4). The constraints added to the
formulated optimization problem are the following:
g1 (~x) = (V ? 105, 944) ? 0 , (5.10)
g2 (~x) = (88, 286? V ) ? 0 , (5.11)
g3 (~x) =
(
1?
m
n
tan ?
tan ? ? tan ?
)
< 0 , (5.12)
g4 (~x) = (max (?local)? ?max) ? 0 . (5.13)
Functions g1 through g3 are used for the blunt leading edge waverider case (J = 3)
and g4 is added for sharp leading edge case (J = 4).
Objective function f1 is the lift-to-drag ratio of the vehicle preceded with a
negative sign. It is calculated using the aerodynamic methods described in Chapter
2. The goal here is to maximize the lift-to-drag ratio. To keep with a minimization
convention, the optimizer attempts to minimize f1. Objective function f2 is the
volumetric efficiency, ?, of the vehicle preceded with a negative sign. The volumetric
98
efficiency is defined by comparing the vehicle?s wetted surface area to the surface
area of a sphere of equal volume. This was done due to the fact that a sphere is the
most volumetrically efficient shape possible. The volumetric efficiency is defined by,
? ?
(36pi)
1
3 V
2
3
Sw
, (5.14)
where V is the vehicle internal volume and Sw is the wetted surface area. The goal
of this vehicle design problem is to maximize the volumetric efficiency by minimizing
f2. Objective function f3 is a measure of the vehicle?s trimmability and is defined as
the magnitude of the moment coefficient given by Equation 2.14. The smaller the
magnitude of f3, the smaller the control moments need to be to trim the vehicle at
the design attitude. Therefore, f3 is minimized by the optimizer. Objective function
f4 is the maximum heat flux that occurs for any panel on the blunted leading edge of
the vehicle. This is an important parameter for heat shield design and TPS sizing.
This value is minimized by the optimizer.
Objective functions are scaled so that they are all on the same order of magni-
tude and so they also convert the problem to a minimization type [98]. To accomplish
this, an estimate of the best and worst possible values of each objective function are
used. These points are called ?good? points and ?bad? points, respectively. Scaling
an objective function using these points is accomplished by the following equation:
fq,scaled =
fq ? fq,good
fq,bad ? fq,good
. (5.15)
The scaled values of the objective functions are in the range of [0,1], where a smaller
the value is more optimum. The values of the corresponding good and bad points
for each objective function are given in Table 5.2 in Section 5.4. Note that the good
99
and bad points were chosen by observing test results of preliminary optimization
runs and may not always successfully scale the objective functions to lie exactly in
the range of [0,1].
The g1 and g2 constraints are required to keep the internal volume within the
requirements set in Table 5.1. Constraint g3 is necessary for a geometrically feasible
concave vehicle (when ? < ?) to be generated by the method outlined in Section
2.1. This constraint keeps the bottom surface of the waverider from crossing the top
surface [30]. The last constraint, g4, is necessary for an attached shock condition
at the leading edge as discussed in Chapter 2. Constraint g4 is only applied for the
sharp leading edge cases. The most common way of handling constraints in genetic
algorithms is by the use of a penalty method [98]. An exterior penalty method was
used to handle the constraints in this work. The idea is to alter the fitness value
of an individual by a severe penalty if it violates a constraint. The exterior penalty
method applied to the scaled objective functions is given by
fq,con = fq,scaled +K
J?
j=1
{gj} , (5.16)
where K is a penalty factor of 105 and the bracket function {gj} is defined as
{gj} =
?
????
????
gj, if gj ? 0
0, if gj< 0
. (5.17)
100
5.3 Influence of Asymmetric Boundary-Layer
5.3.1 Multi-Objective Optimization
The multi-objective optimization problem under the influence of ABLT is for-
mulated to be the same as the baseline :
min
~x
fq,scaled (~x) q = 1, 2, ..., Q
s.t.
gj (~x) ? 0 j = 1, 2, ..., J
~x ? X
X = {~x|~xlb ? ~x ? ~xub} ,
(5.18)
where now a single isolated disturbance is introduced at the leading edge of the
top surface of the vehicle. The lateral position of the disturbance is added to the
optimization problem as an integer valued environmental parameter p. The value
of p represents a spanwise position of the isolated roughness in increments of one-
eighth the total span of the waverider. Four lateral positions were used spanning
outward from the centerline of the vehicle as illustrated in Figure 5.1.
The idea behind incorporating an ABLT event to the multi-objective opti-
mization problem was to identify vehicles that are the most capable of withstanding
the effect of a turbulence wedge and remaining controllable. A vehicle state is said
to be controllable if there exists a control signal, over some finite time interval, that
can bring the vehicle to a future state equal to the desired equilibrium state [99].
In other words, any vehicle state which the control system has enough authority to
bring it to a trimmed condition is considered a controllable state. The set of all
101
??
?
?
w
8
w
8
w
8
w
8
p = 1
p = 2
p = 3
p = 4
Figure 5.1: Parametric locations of leading edge isolated roughness (planform view).
states which are controllable define the vehicle controllability region. The size of
the controllability region can be considered as a measure of the allowable level of
external disturbances [100].
With these definitions, two things can be taken into account when predicting
if an external event will cause the loss of a vehicle: the size of the controllability
region and the size of the external disturbance. Minimizing the chance of losing
a vehicle due to an external disturbance can be accomplished by maximizing the
controllability region of the vehicle and/or minimizing the size of the external dis-
turbance. To accomplish the former, in the current work, the dynamics of the range
of waveriders modeled would need to be calculated. This, as discussed in Section
1.4.1, can be difficult and potentially inaccurate. In addition, assumptions pertain-
ing to the control system (flaps, thrusters, etc.) would need to be made. This is
a very difficult task to do a priori without any knowledge of vehicle shape and in-
102
ternal volume distribution. Instead, the best first step is to accomplish the latter
and minimize the external disturbance on the vehicle. This requires no restrictive
assumptions to be made about the type and strength of the control system on the
vehicle. This, in turn will minimize the chances of the external disturbance moving
the system out of the controllability region, where no controller is able to keep the
system stable [100].
To capture this, a controllability objective function was defined by modifying
the definition of f3 in Equation 5.8. The resulting objective function f3 specifically
for Section 5.3.1 is defined by
f3 (~x) =
1
4
(f3,span (~x, 1) + f3,span (~x, 2) + f3,span (~x, 3) + f3,span (~x, 4)) , (5.19)
where
f3,span (~x, p) = abs
(
?~CMCG,tot?with??~CMCG,tot?without
)
. (5.20)
Minimizing f3,span at a particular spanwise position minimizes the absolute value
of the difference between the length of the moment vector with a turbulence wedge
present and the length of the moment vector without a turbulence wedge present.
When Equation 5.19 is minimized, the effect due to an isolated disturbance which
may occur at any of the four spanwise positions is minimized. Therefore, the vehicle
design which minimizes Equation 5.19 best maintains its controllability character-
istics if an isolated roughness event would happen to occur at any of these four
spanwise locations. The the rest of the objective functions, constraint functions,
scaling method, and constraint handling method described in Section 5.2 is the
same for Section 5.3.1.
103
Calculating the moment coefficient characteristics at four different spanwise
locations in Equation 5.19 increased the computational expense of the optimization
problem. To more efficiently use the calculations that are being performed, the
objective function f3 in Equation 5.18 was changed to a robust counterpart function
using the Tsutsui perturbation approach [72,73]:
min
~x
F3,scaled (~x, p (t)) = f3,span,scaled (~x, p+ ?p (t))
s.t.
gj (~x, p) ? 0 j = 1, 2, ..., J
~x ? X
X = {~x|~xlb ? ~x ? ~xub}
(p+ ?p (t)) ? UD (p?)
p? = {1, 2, 3, 4} ,
(5.21)
where t is a generation index and p + ?p (t) is a random value (1, 2, 3, or 4) at
generation t.
The goal of the original optimization problem and f3 of Equations 5.18 and
5.19 was to identify waverider geometries that are robust to an unintended roughness
element that may appear anywhere on the top surface of the waverider. To mirror
this effect in Equation 5.21, the value of p+?p (t) is set randomly at the start of each
generation within the MOGA solver. All calculations performed for each individual
in the population for a given generation uses that particular value of p + ?p (t).
Therefore, for a large number of generations, each spanwise location (1, 2, 3, or 4)
would tend to be equally represented and the final optimum design should be robust
to a turbulence wedge which may occur at any of these spanwise locations. The last
step when evaluating Equation 5.21 within a MOGA solver was to take the solution
104
set of Pareto-optimal points from the evaluate them for the original f3 in Equation
5.19.
5.3.2 Robust Optimization Formulations
The robust optimization formulations aid in further exploring and analyzing
the design space by focusing on uncertainty in the position of the surface roughness
parameter. The waverider geometries used in this section are all assumed to have an
infinitely sharp leading edge. To quantify the effect of an asymmetric flow field on
the controllability of the vehicle, the sole objective function used in the optimization
studies in this section was the following:
f (~x, p) = abs
(
?~CM?with??~CM?without
)
, (5.22)
where ?~CM? is the magnitude of the moment coefficient vector taken about the
center of gravity of the vehicle. Note, ?~CM? is shorthand and is equivalent to
?~CMCG,tot?. Subscripts with and without denote conditions with and without the
asymmetric flow field (turbulence wedge) present, respectively. This formulation
isolates the effect of the roughness and defines the moment increment, ?~CM , which
needs to be attenuated using vehicle stability and control mechanisms. Minimizing
this moment increment minimizes the amount of control needed to keep the system
at trimmed conditions and gives the best likelihood that the vehicle will not be
perturbed outside its controllability region.
Parameter p denotes a lateral distance of the isolated roughness on the leading
edge of the top surface of the vehicle. Parameter p spans from the center line of
105
the vehicle (p = 0.5) to the full half-span of the vehicle (p = 0.0). Therefore, the
bounded environmental parameter uncertainty set for the hypersonic vehicle design
problem is defined by
p ? P = {p|0.0 ? p ? 0.5} . (5.23)
Note, the parameter uncertainty space here is continuous and connected, unlike the
discrete uncertainty space in Section 5.3.1. Figure 5.2 provides example locations
of the isolated roughness for various values of p.
?
?
?
?
p = 0.125
p = 0.25
p = 0.375
p = 0.5
p0 0.5
Figure 5.2: Example locations of leading edge isolated roughnesses as function of p
(planform view).
Using the terminology in Chapter 4, the original optimization problem is for-
mulated as follows:
min
~x
f (~x, p) = abs
(
?~CM?with??~CM?without
)
~x ? X = {~x|~xlb ? ~x ? ~xub} ,
(5.24)
106
where ~x is the design vector defined by Equation 5.1. The design vector belongs to
the design vector set space X defined by Equation 5.3, where the lower bound on
variables n and m are 0.05.
Each of the four methods described in Section 4.2 forms a robust counterpart
optimization problem to the original problem defined by Equation 5.24. The formu-
lated robust counterpart optimization problem for the robust regularization-based
method is given by,
min
~x
FI (~x) = max
p?p?
abs
(
?~CM (~x, p) ?with??~CM (~x, p) ?without
)
~x ? X
p? = {p|0.0, 0.05, 0.1, . . . , 0.45, 0.5} ? P .
(5.25)
Note here, the supremum over the uncertainty set, sup
p?P
, was approximated by a
maximization over a discrete subset p?. This was done as an attempt to solve the inner
optimization problem with less computational expense than a genetic algorithm or
gradient-based method, while maintaining accuracy in the solution.
The explicit averaging method forms the following robust counterpart opti-
mization problem:
min
~x
FII (~x, p) =
10?
i=1
1
10
abs
(
?~CM (~x, p+ ?pi) ?with??~CM (~x, p+ ?pi) ?without
)
~x ? X
(p+ ?pi) ? U (0.0, 0.5) .
(5.26)
Note here, the Monte Carlo integration was performed using 10 samples. This
number of samples was found to perform adequately while maintaining a feasible
computational expense. The value of the environmental parameter was chosen to
107
be sampled from a uniform distribution over the uncertainty set. This distribution
was chosen as a baseline and is suggested to be used if the designer has no educated
guess on the probability of a roughness occurring on the surface. If, for example,
the designer knows that the vehicle is made up of different ablative materials, an
educated guess of where a roughness is more likely occur is possible. Or if the
vehicle has machining marks, screw heads, or rivets, areas with a higher probability
of transition can be identified. These example scenarios would lend towards a non-
uniform uncertainty PDF.
The perturbation method forms the following robust counterpart optimization
problem:
min
~x
FIII (~x, p (t)) = abs
(
?~CM (~x, p+ ?p (t)) ?with??~CM (~x, p+ ?p (t)) ?without
)
~x ? X
(p+ ?p (t)) ? U (0.0, 0.5) .
(5.27)
The robust multi-objective approach forms the following robust MOGA coun-
terpart optimization problem:
min
~x
FIVm (~x, p?m (t)) = abs
(
?~CM (~x, p?m (t)) ?with??~CM (~x, p?m (t)) ?without
)
m = 1, 2, 3, 4, 5
~x ? X
p?1 (t) ? U (0.0, 0.1)
p?2 (t) ? U (0.1, 0.2)
p?3 (t) ? U (0.2, 0.3)
p?4 (t) ? U (0.3, 0.4)
p?5 (t) ? U (0.4, 0.5) .
(5.28)
108
Note the robust MOGA counterpart formulation set M = 5, creating five objec-
tive functions. This value was chosen as a trade-off between accuracy of the post-
optimality handling prediction (Equations 4.11 and 4.13) and the MOGA problem
complexity (as M gets larger, the problem becomes more difficult to solve [94,95]).
5.4 MOGA Method and Parameters
5.4.1 Multi-Objective Studies MOGA Parameters
The multi-objective optimization problems outlined in Sections 5.2 and 5.3.1
were solved using the MATLAB GADS toolbox. The multi-objective genetic al-
gorithm solver in the GADS toolbox employs a modified version of an elitist non-
dominated sorting genetic algorithm method, NSGA-II [94]. Table 5.2 lists all the
parameters which were necessary to run the MOGA solver.
Deb suggests that for a multi-objective optimization problem with four objec-
tive functions a population size of at least 100 is to necessary to maintain reasonable
effectiveness using the elitism scheme (see Figure 276 in Deb [94]). Therefore, the
population size was set to 100 for all optimization runs. The MOGA solver was
set to terminate after a particular number of generations was computed. The num-
ber of generations was determined by running several test cases and observing the
change in the Pareto spread from generation to generation. For the baseline multi-
objective studies in Section 5.2, it was decided that after 500 generations the average
change in the Pareto spread from generation to generation was not significant. In
other words, the Pareto-optimal front was not making significant progress towards
109
more optimum designs. For the multi-objective optimization studies under the in-
fluence of an asymmetric boundary-layer transition (Section 5.3.1), this was found
to be the case after 2000 generations. These numbers were set subjectively and
may not always guarantee that the final Pareto-optimal points adequately cover the
entire design space or accurately represent the true Pareto-optimal front. Genetic
algorithms have no guarantee of global optimality and/or sufficient coverage of the
design space. In an attempt to minimize the chances of obtaining an inadequate or
inaccurate set of Pareto-optimal points, the current study ran each of the four cases
25 times.
110
T
ab
le
5.
2:
M
O
G
A
pa
ra
m
et
er
s
fo
r
th
e
m
ul
ti
-o
b
je
ct
iv
e
op
ti
m
iz
at
io
n
st
ud
ie
s.
P
ar
am
et
er
C
as
e
B
as
el
in
e
ca
se
:
sh
ar
p
B
as
el
in
e
ca
se
:
bl
un
t
A
B
LT
ca
se
:
sh
ar
p
A
B
LT
ca
se
:
bl
un
t
P
op
ul
at
io
n
si
ze
10
0
10
0
10
0
10
0
N
um
b
er
of
ge
ne
ra
ti
on
s
50
0
50
0
20
00
20
00
N
um
b
er
of
M
O
G
A
ru
ns
25
25
25
25
F
1,
go
od
10
10
10
10
F
2,
go
od
1
1
1
1
F
3,
go
od
0
0
0
0
F
4,
go
od
N
/A
30
W
/c
m
2
N
/A
30
W
/c
m
2
F
1,
ba
d
0
0
0
0
F
2,
ba
d
0
0
0
0
F
3,
ba
d
0.
05
0.
05
5
?
10
?
6
5
?
10
?
6
F
4,
ba
d
N
/A
40
0
W
/c
m
2
N
/A
40
0
W
/c
m
2
111
5.4.2 Robust Optimization Studies MOGA Parameters
The four robust counterpart optimization problems described by Equation 5.25
through Equation 5.28 were solved using the MATLAB GADS toolbox. The MOGA
solver in the GADS toolbox employs a modified version of the elitist non-dominated
sorting genetic algorithm NSGA-II [94]. Table 5.3 lists all the genetic algorithm
settings which were necessary to run the MOGA solver for each robust formulation.
Table 5.3: Genetic algorithm options for the robust optimization studies.
Genetic algorithm option Setting
Population size 50 (50?M for MOGA)
Initial population Random within xlb and xub
Fitness scaling Normalized to range between 0 and 1
Selection Stochastic uniform selection
Crossover Randomized scattered
Crossover fraction 0.8
Mutation Gaussian
Convergence criterion Run for 1000 generations (2000 for MOGA)
A population size of 50 was used for the single-objective methods (robust
regularization and the two aggregation methods). The methods, when using a pop-
ulation of this size, were observed to make significant progress from generation to
generation in the initial stages of the simulation while maintaining a reasonable
computational expense. A population size of 50M (or 250 for five objectives) was
used for the robust MOGA method. This ensured that the population size exceeded
the minimum values necessary to maintain reasonable effectiveness when using the
112
NSGA-II elitism scheme (see Figure 276 in Deb [94]). Objective function values
were all scaled using the same method outlined in Section 5.2. The simulation
was configured to run for a set number of generations. The number of generations
was set such that the generation histories show a significant level of convergence.
Generation histories are provided in the following sections illustrating the extent of
convergence for each method as a function of generation number.
113
Chapter 6: Results
This chapter presents the results of the optimization problems formulated in
Chapter 5. The results presented here were found by solving the baseline optimiza-
tion problems and the ABLT influenced optimization problems (Sections 5.2 and
5.3, respectively), under the assumptions outlined in Section 5.1 using the method
and parameters given in Section 5.4.
This chapter parallels Chapter 5, where the baseline optimization results are
given in Section 6.1 and the ABLT influenced optimization results are given in Sec-
tion 6.2. Within this chapter, two categories of results are presented: multi-objective
optimization results and robust optimization results. Multi-objective optimization
results are found in Sections 6.1 and 6.2.1, while the robust optimization results are
found in Section 6.2.2. These two categories of results are presented with slightly
different nomenclature, and were found using slightly different convergence criteria.
The following paragraphs will provide a brief introduction into how the results were
found and will be presented in the subsequent sections in this chapter.
The multi-objective optimization results in Sections 6.1 and 6.2.1 are split into
four cases: a baseline case with infinitely sharp leading edges, a baseline case with
blunted leading edges, an ABLT case with sharp leading edges, and an ABLT case
114
with blunted leading edges. Each of the four optimization cases described in Table
5.2 were run 25 times using the MOGA solver, producing 25 sets of Pareto-optimal
points per case. The 25 sets were compared for dominance and combined to create
an overall set of non-dominated Pareto-optimal points. The results presented in
Sections 6.1 and 6.2.1 analyze the two following types of points that exist in the
overall set Pareto-optimal points for each case:
? extreme Pareto-optimal points
? L2-norm optimum Pareto-optimal point
The extreme Pareto-optimal points are the points in the set which have the min-
imum value for each individual objective function. Designs which represent the
extreme Pareto-optimal points that minimize the scaled objective functions will be
termed ?best? designs for the corresponding objective functions. These points help
illustrate the trade-offs that exist between the objectives, showing the full range
of compromise. The L2-norm optimum Pareto-optimal point represents the design
that minimizes the compromise between all the objectives by minimizing
?
?
?
?
Q?
q=1
(
abs (fq ? fq,min)
2) , (6.1)
where fq,min is the extreme Pareto-optimal point for objective q. A simplified two
objective example is given in Figure 6.1 to help illustrate where these types of
Pareto-optimal points exist on the Pareto-optimal front.
The robust optimization results in Sections 6.2.2.1, 6.2.2.2, and 6.2.2.3 are
found from solving a single-objective optimization problem. These sections therefore
115
f1
f2
L2-norm optimum
f2,min, extreme Pareto-optimal point
f1,min, extreme Pareto-optimal point
Pareto-optimal points
Figure 6.1: Two objective function example showing extreme Pareto-optimal points
and L2-norm optimum Pareto-optimal point.
provide a single robust optimum design. The results in Section 6.2.2.4 are found
by solving a multi-objective optimization problem and therefore exist as a set of
optimum designs. The nuances to this particular Pareto-optimal set of designs is
discussed in detail in Section 4.2.3. In addition, an illustrative example problem
formulation and corresponding results are presented in Section 4.3.
6.1 Baseline Multi-Objective Optimal Geometries
6.1.1 Infinitely Sharp Case
Two classes of shapes were identified from the extreme Pareto-optimal designs
for the infinitely sharp baseline case: wedge-like and cone-like. Wedge-like shapes
are produced when the variable n approaches zero and ? ? ?. Cone-like shapes are
produced when n approaches unity for a moderately convex vehicle (? > 2?). The
116
best f1 and best f3 designs were wedge-like in shape and had a small wedge angle.
The best f2 design was cone-like. Isometric views of the best designs are shown
in Figure 6.2. The corresponding objective function values and design vectors are
presented in Table 6.1 and Table 6.2, respectively.
(a) Best f1 (b) Best f2 (c) Best f3
Figure 6.2: Extreme Pareto-optimal point vehicle designs with best f1, f2, f3 values
for the infinitely sharp leading edge baseline case (isometric views).
The wedge-like shape is optimal for maximizing the lift-to-drag ratio because
it best aligns the resultant pressure force vector with the z-axis. The lower surface
has little curvature and has a small wedge angle ?, therefore the force components
in the drag and side force directions are minimal . The cone-like shape is the most
volumetrically efficient shape in the design space, due to its rounded lower surface.
As a cone-like shape becomes more wedge-like or flat, the volume-to-surface area
ratio decreases.
The total moment force coefficient is comprised of a both a pressure and viscous
component. For the wedge-like and cone-like shapes, it was found that the moment
coefficient is dominated by the pressure force component. Therefore, the best design
for f3 was driven by the moments created by pressure forces. The wedge-like design
was best because it minimizes the pressure induced pitching moment and balances
117
any pressure induced roll or yaw moments due to symmetry. The wedge shape has
the minimal pitching moment because is has a small wedge angle at all spanwise
locations, thus minimizing the pressure force on the lower surface of the vehicle.
The extreme Pareto-optimal designs for the infinitely sharp baseline case show
that a trade-off exists between wedge-like and cone-like shapes. If a design aims
solely for a high lift-to-drag ratio and a low moment coefficient, it will have poor
a volumetric efficiency. For example, the wedge-like optimum shape for the design
point in this study had a high lift-to-drag ratio and low moment coefficient mag-
nitude of about 8 and 2.8? 10?5, respectively, yet suffered in volumetric efficiency
with a value of about 0.2. Any gain in volumetric efficiency will come at the cost
of the performance (lift-drag-ratio) and trimmability (moment coefficient). In this
study, the cone-like shape had approximately three times the volumetric efficiency
of the wedge-like shape, yet saw about a 75% reduction in lift-to-drag ratio and a
jump in moment coefficient magnitude of four orders of magnitude compared to the
wedge-like counterpart.
The L2-norm optimum design is the design which equally balances the three
objectives. The L2-norm design for the infinitely sharp baseline case is shown in
Figure 6.3. The objective function values and design vector are listed in Tables
6.1 and 6.2. Geometrically, the L2-norm design is a hybrid shape balancing both
wedge-like and cone-like shape characteristics. The design vector values reinforce
this balance. The value of n is 0.127. This balance between zero and one keeps
the long slender shape of the cone while still flattening the shape out. The ? and
? values are about 5.956?. and 4.9? respectively. This produces a slightly convex
118
shape of a moderate wedge angle. The L2-norm shape exhibits a lift-to-drag ratio of
6.185, ? of 0.458, and moment coefficient magnitude of 1.451? 10?3. This amounts
to approximately twice the volumetric efficiency of the wedge-like shape with a
25% reduction in lift-to-drag ratio and a two orders of magnitude jump in moment
coefficient magnitude compared to the wedge-like shape. These performance and
trimmability characteristics are the values of minimum compromise between the
three objectives.
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.3: L2-norm Pareto-optimal point vehicle design for the infinitely sharp
leading edge baseline case.
The current work found that a wedge-like shape produces a maximum lift-to-
drag ratio vehicle. This was verified by comparing the results of the current study to
the results found by Starkey [30]. Starkey found a wedge-like optimum shape when
using a single-objective optimization study to find the maximum value of lift-to-drag
ratio multiplied by internal volume. Starkey also optimized for maximum lift-to-drag
ratio multiplied by volumetric efficiency. This objective function is analogous to a
L2-norm optimum between lift-to-drag ratio and volumetric efficiency. Therefore, it
is no surprise that the L2-norm optimum shape in the current work exhibits similar
characteristics to the single-objective optimum in Starkey?s work.
119
T
ab
le
6.
1:
P
ar
et
o-
op
ti
m
al
p
oi
nt
ob
je
ct
iv
e
fu
nc
ti
on
va
lu
es
fo
r
th
e
in
fin
it
el
y
sh
ar
p
le
ad
in
g
ba
se
lin
e
ed
ge
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
L
if
t
D
ra
g
?
?
~ C
M
C
G
,t
o
t
?
F
ig
ur
e
B
es
t
f 1
8.
13
2
0.
22
9
2.
93
2
?
10
?
5
F
ig
ur
e
6.
2(
a)
B
es
t
f 2
1.
96
6
0.
63
3
4.
03
1
?
10
?
2
F
ig
ur
e
6.
2(
b)
B
es
t
f 3
8.
13
1
0.
22
9
2.
76
7
?
10
?
5
F
ig
ur
e
6.
2(
c)
L
2
-n
or
m
6.
18
5
0.
45
8
1.
45
1
?
10
?
3
F
ig
ur
e
6.
3
T
ab
le
6.
2:
P
ar
et
o-
op
ti
m
al
p
oi
nt
de
si
gn
ve
ct
or
va
lu
es
fo
r
th
e
in
fin
it
el
y
sh
ar
p
le
ad
in
g
ed
ge
ba
se
lin
e
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
n
m
?
(d
eg
)
?
(d
eg
)
L
(f
t)
w
(f
t)
F
ig
ur
e
B
es
t
f 1
0.
00
3
0.
66
8
2.
70
2
2.
95
3
19
0.
40
5
11
2.
56
7
F
ig
ur
e
6.
2(
a)
B
es
t
f 2
1.
00
0
0.
14
6
9.
99
7
2.
43
0
20
5.
75
7
52
.9
75
F
ig
ur
e
6.
2(
b)
B
es
t
f 3
0.
00
1
0.
66
4
2.
70
1
2.
95
4
19
0.
40
5
11
2.
56
8
F
ig
ur
e
6.
2(
c)
L
2
-n
or
m
0.
12
7
0.
12
1
5.
95
6
4.
90
0
20
6.
04
0
51
.7
71
F
ig
ur
e
6.
3
120
6.1.2 Blunted Case
The extreme Pareto-optimal design results for the blunted baseline case, shown
in Figure 6.4, follow the same patterns as the sharp case. The designs are either
wedge-like or cone-like for the same corresponding best f1, f2, and f3 objectives. The
best f1 design for the blunted case is not as wedge-like as what was found in the
sharp case. However, this not a product of adding a blunted leading edge. A fully
wedge-like extreme Pareto-optimal shape is possible with a blunted leading edge
because it was found for the best f1 blunted shape in Section 6.2.1 (Figure 6.13(a)).
The optimization simulations in Section 6.2.1 were run for more generations and the
Pareto-optimal designs had more time to progress towards the wedge-like extreme
design. If given more generations, the best f1 blunted baseline design would match
the wedge shapes in Section 6.2.1 and the previous section (Figures 6.13(a) and
6.2(a)). The maximum possible lift-to-drag ratio for the best blunted f1 shape does
however differ from the sharp case. Blunting the leading edge adds more pressure
drag and decreases the lift-to-drag ratio by about 9%. The optimum design therefore
tries to mitigate the loss of lift-to-drag ratio by decreasing the separation distance
d to its minimum possible value.
The best f2 result is the same as the sharp case, yet now the optimizer gains
some volumetric efficiency (an increase of about 10%) by maximizing the separation
distance to its maximum possible value. This shape is also the best f4 design because
the maximum value of d gives a minimum peak heat flux on the leading edge. The
resulting design for f4 may show that minimizing peak heat flux does not drive the
121
design, and will always follow the best volumetrically efficient design. Future studies
may want to also factor in a method of minimizing the total heat load on the leading
edge.
The characteristics of the best f3 and L2-norm designs agree with the analysis
presented in the infinitely sharp baseline results section. The blunted L2-norm
shape, shown in Figure 6.5, favors a more cone-like shape compared to the sharp
L2-norm shape. This is due to the inclusion of the f4 objective function. The best f4
shape is cone-like, which weights the L2-norm more towards a cone-like shape. The
corresponding objective function values and design vectors for the extreme Pareto-
optimal designs and L2-norm design are presented in Table 6.3 and Table 6.4.
(a) Best f1 (b) Best f2 (c) Best f3 (d) Best f4
Figure 6.4: Extreme Pareto-optimal point vehicle designs with best f1, f2, f3, f4
values for the blunted leading edge baseline case (isometric views).
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.5: L2-norm Pareto-optimal point vehicle design for the blunted leading
edge baseline case.
122
T
ab
le
6.
3:
P
ar
et
o-
op
ti
m
al
p
oi
nt
ob
je
ct
iv
e
fu
nc
ti
on
va
lu
es
fo
r
th
e
bl
un
te
d
le
ad
in
g
ed
ge
ba
se
lin
e
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
L
if
t
D
ra
g
?
?
~ C
M
C
G
,t
o
t
?
m
ax
(q
w
)
(W
/c
m
2
)
F
ig
ur
e
B
es
t
f 1
7.
41
6
0.
28
0
5.
26
6
?
10
?
4
34
6.
06
7
F
ig
ur
e
6.
4(
a)
B
es
t
f 2
1.
71
8
0.
70
2
3.
83
9
?
10
?
2
39
.3
85
F
ig
ur
e
6.
4(
b)
B
es
t
f 3
3.
33
2
0.
24
1
1.
29
5
?
10
?
5
85
.4
38
F
ig
ur
e
6.
4(
c)
B
es
t
f 4
1.
71
8
0.
70
2
3.
83
9
?
10
?
2
39
.3
85
F
ig
ur
e
6.
4(
d)
L
2
-n
or
m
4.
45
8
0.
53
5
5.
16
2
?
10
?
3
95
.4
05
F
ig
ur
e
6.
5
T
ab
le
6.
4:
P
ar
et
o-
op
ti
m
al
p
oi
nt
de
si
gn
ve
ct
or
va
lu
es
fo
r
th
e
bl
un
te
d
le
ad
in
g
ed
ge
ba
se
lin
e
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
n
m
?
(d
eg
)
?
(d
eg
)
L
(f
t)
w
(f
t)
d
(f
t)
F
ig
ur
e
B
es
t
f 1
0.
11
5
0.
86
9
2.
90
1
3.
53
1
21
2.
03
6
88
.9
84
0.
02
1
F
ig
ur
e
6.
4(
a)
B
es
t
f 2
0.
98
4
0.
17
3
9.
98
9
2.
65
3
16
5.
17
0
53
.1
58
1.
64
0
F
ig
ur
e
6.
4(
b)
B
es
t
f 3
0.
00
3
0.
84
9
2.
36
8
1.
50
1
21
2.
57
2
90
.5
55
0.
35
2
F
ig
ur
e
6.
4(
c)
B
es
t
f 4
0.
98
4
0.
17
3
9.
98
9
2.
65
3
16
5.
17
0
53
.1
58
1.
64
0
F
ig
ur
e
6.
4(
d)
L
2
-n
or
m
0.
36
4
0.
60
8
8.
41
7
3.
54
5
19
7.
91
6
59
.2
52
0.
28
2
F
ig
ur
e
6.
5
123
6.2 Optimal Geometries under Asymmetric Boundary-Layer Condi-
tions
6.2.1 Multi-Objective Results
The extreme Pareto-optimal designs for f1, f2, and f4 for all of the ABLT
influenced optimization cases are very similar to the designs found for the baseline
cases. Both the wedge-like and cone-like shapes were found as extreme Pareto-
optimal designs and the analysis and discussion for this section parallels that of
Section 6.1. The analyses for the L2-norm designs in Section 6.1 are also applicable
in this section. In lieu of these similarities, the following section will only focus on
the f3 extreme Pareto-optimal designs and the influence of an isolated disturbance
and asymmetric boundary-layer.
6.2.1.1 Infinitely Sharp Case
Minimizing f3 under the influence of an ABLT event aims to minimize the
change in the magnitude of the moment vector due to an asymmetric boundary-layer
at multiple possible spanwise locations. Subtracting the moment vector magnitude
for a design without a turbulence wedge from the moment vector magnitude for a
design with a turbulence wedge isolates the viscous effect on vehicle controllability.
In other words, the moments due to pressure are subtracted out.
To aid in the analysis of the results, the best and worst f3 designs are given
in Figures 6.6(c) and 6.6(d). The design that is most robust to spanwise turbulence
124
wedge disturbances is the small wedge angle wedge-like geometry. The wedge-like
geometry had a f3 value of 1.453? 10?7. The design that is most sensitive to the
presence of a turbulence wedge is a cone-like shape. The cone-like and wedge-like
shapes illustrate the trade-off that exists between volumetric efficiency, lift-to-drag
ratio, and controllability. The cone-like shape had 3.5 times the volumetric effi-
ciency of the wedge-like shape, yet compromised with an approximate 75% decrease
in lift-to-drag ratio and one full order of magnitude jump in the value of f3 (the
controllability factor) when compared to the wedge-like counterpart. Note that the
base plane cross section of the cone-like shape differs from the extreme cone-like
shape baseline results. The worst f3 shape for the ABLT case has a more rounded,
circular top surface (n ? 0.5), whereas the baseline design has a flatter upper surface
(n ? 1).
(a) Best f1 (b) Best f2 (c) Best f3 (d) Worst f3
Figure 6.6: Extreme Pareto-optimal point vehicle designs with best f1, f2, f3, values
and the design with the worst f3 value for the ABLT infinitely sharp leading edge
case (isometric views).
As expected, the L2-norm design, seen in Figure 6.7, is a balance of these two
extreme designs. The L2-norm design had approximately 2.3 times the volumetric
efficiency of the wedge-like shape with a 25% reduction in lift-to-drag ratio and
125
about five times (half an order of magnitude) the moment coefficient magnitude of
the wedge-like shape.
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.7: L2-norm Pareto-optimal point vehicle design for the ABLT infinitely
sharp leading edge case.
An understanding of why these shapes are the best, worst, and L2-norm designs
with respect to f3 requires both a qualitative and quantitative analysis. A qualitative
look at how the turbulence wedge geometrically covers the vehicles is necessary to
visualize the effect on the moments. Views of the turbulence wedge as applied to
these the best, worst, and L2-norm Pareto-optimal designs as a function of p are
shown in Figures 6.8, 6.9, and 6.10. The corresponding quantitative look of the
spanwise effects of the turbulence wedge is given as the spanwise distribution of
F3,scaled, shown in Figure 6.11. Note, the figure presents the robust counterpart
function F not the original function f (see Equation 5.21 in Section 5.3.1).
126
(a) p = 1 (b) p = 2 (c) p = 3 (d) p = 4
Figure 6.8: Isometric views of the best f3 Pareto-optimal design with a spanwise
disturbance at p = 1, 2, 3, and 4 for the ABLT infinitely sharp case.
(a) p = 1 (b) p = 2 (c) p = 3 (d) p = 4
Figure 6.9: Isometric views of the worst f3 Pareto-optimal design with a spanwise
disturbance at p = 1, 2, 3, and 4 for the ABLT infinitely sharp case.
(a) p = 1 (b) p = 2 (c) p = 3 (d) p = 4
Figure 6.10: Isometric views of L2-norm Pareto-optimal design with a spanwise
disturbance at p = 1, 2, 3, and 4 for the ABLT infinitely sharp case.
127
1 1.5 2 2.5 3 3.5 40
0.1
0.2
0.3
0.4
0.5
0.6
p
F 3,r
efo
rm,
sca
led
 
 Worst F3 Pareto extreme pointBest F3 Pareto extreme pointL2-norm Pareto point
F 3
,s
ca
le
d
p
f
f3
Figure 6.11: Spanwise distribution of F3,scaled for best, worst, and L2-norm vehicle
designs for the ABLT infinitely sharp leading edge case.
128
The wedge-like shape keeps the spanwise value of F3,scaled at a low value, fluc-
tuating only slightly as p goes from 1 to 4. The cone-like shape starts at a low value
similar to the wedge-like shape at p = 1. But as the turbulence wedge approaches
the center of the vehicle, the value of F3,scaled increases monotonically. The roles of
the three components to F3 (roll, pitch, and yaw) are important in understanding
why the wedge and cone shapes produce these curves. The x component, or roll
component, is zero due to the assumption that the direction of the shear force vector
is purely in the x direction (see Chapter 2). The y and z components to F3 are the
pitch and yaw components, respectively. The spanwise distributions of the y and z
components of F3 are shown in Figure 6.12.
At p = 4, the moment created by the turbulence wedge is purely a pitching
moment. Since the wedge-like (best) design has a very small moment arm in the
z direction, the pitching moment value is low. Conversely, the cone-like (worst)
design has a much larger moment arm, so the pitching moment value is large. As
the turbulence wedge progresses away from the centerline, the moment arm for the
wedge-like vehicle stays constant, so the pitching moment stays the same. The
moment arm for the cone-like vehicle decreases as the wedge moves away from the
centerline and the area covered by the turbulence wedge becomes smaller. Therefore,
the cone-like design shows a decrease in pitching moment as p goes from 4 to 1. This
effect is shown by the curves in Figure 6.12(a).
The yaw moment (Figure 6.12(b)) for both the best and worst designs start
at zero and becomes more negative as the turbulence wedge moves away from the
centerline. This is due to the increase in the moment arm in the y direction. The area
129
covered by the turbulence wedge for the wedge-like (best) design remains constant,
so the yaw moment is only affected by the linear increase in moment arm length.
Consequently, the yaw moment becomes more negative nearly linearly as p goes
from 4 to 1. The area covered by the turbulence wedge for the cone-like (worst)
design decreases as p goes from 4 to 1 and eventually counteracts the increase in
moment due to an increasing moment arm length. This is shown by the decrease
from p = 4 to p = 2 then increase from p = 2 to p = 1 in the worst design curve.
For the Pareto-optimal shapes, the moment vector without a turbulence wedge
points solely in the y direction (As described in the baseline results section). Con-
sequently, relatively small changes to the moment vector?s z component will change
the direction of the vector, but only weakly affect the change in magnitude of the
vector. In comparison, relatively small changes in the y component will change the
direction of the vector, but strongly affect the change in magnitude. F3 was defined
as the change in the magnitude of the moment vector due to a turbulence wedge.
Changes in direction do not affect the objective function value. Therefore, the F3 y
component curves dictate which design is optimum for objective function f3. That
is why the trends in Figure 6.11 closely match those in Figure 6.12(a).
The best designs for f1 and f2 match the analysis given for the baseline case.
For reference, the tables containing the Pareto-optimal point information is given
below. Table 6.5 gives the objective function values and Table 6.6 gives the design
vector values.
130
1 1.5 2 2.5 3 3.5 40
0.5
1
1.5
2
2.5
3
3.5 x 10?6
p
(F 3,
ref
orm
) y
 
 Worst F3 Pareto extreme pointBest F3 Pareto extreme pointL2-norm Pareto point
(F
3)
y
p
f3
f3
(a) F3 y component (pitching moment) vs. p
1 1.5 2 2.5 3 3.5 4?4.5
?4
?3.5
?3
?2.5
?2
?1.5
?1
?0.5
0 x 10?6
p
(F 3,
ref
orm
) z
 
 Worst F3 Pareto extreme pointBest F3 Pareto extreme pointL2-norm Pareto point
(F
3)
z
p
f3
f3
(b) F3 z component (yaw moment) vs. p
Figure 6.12: Spanwise distribution of the y and z components of F3 for best, worst,
and L2-norm vehicle designs for the ABLT infinitely sharp leading edge case.
131
T
ab
le
6.
5:
P
ar
et
o-
op
ti
m
al
p
oi
nt
ob
je
ct
iv
e
fu
nc
ti
on
va
lu
es
fo
r
th
e
A
B
LT
in
fin
it
el
y
sh
ar
p
le
ad
in
g
ed
ge
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
L
if
t
D
ra
g
?
f 3
F
ig
ur
e
B
es
t
f 1
8.
63
3
0.
17
8
1.
45
3
?
10
?
7
F
ig
ur
e
6.
6(
a)
B
es
t
f 2
1.
95
8
0.
63
3
1.
59
1
?
10
?
6
F
ig
ur
e
6.
6(
b)
B
es
t
f 3
8.
47
4
0.
19
0
1.
12
1
?
10
?
7
F
ig
ur
e
6.
6(
c)
W
or
st
f 3
4.
27
6
0.
60
9
1.
96
7
?
10
?
6
F
ig
ur
e
6.
6(
d)
L
2
-n
or
m
6.
47
2
0.
41
2
6.
85
0
?
10
?
7
F
ig
ur
e
6.
7
T
ab
le
6.
6:
P
ar
et
o-
op
ti
m
al
p
oi
nt
de
si
gn
ve
ct
or
va
lu
es
fo
r
th
e
A
B
LT
in
fin
it
el
y
sh
ar
p
le
ad
in
g
ed
ge
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
n
m
?
(d
eg
)
?
(d
eg
)
L
(f
t)
w
(f
t)
F
ig
ur
e
B
es
t
f 1
0.
00
7
0.
03
8
1.
78
9
1.
98
7
21
9.
86
7
11
9.
84
7
F
ig
ur
e
6.
6(
a)
B
es
t
f 2
0.
99
9
0.
14
8
10
.0
00
2.
29
1
19
7.
41
8
50
.9
40
F
ig
ur
e
6.
6(
b)
B
es
t
f 3
0.
00
8
0.
23
8
1.
81
3
2.
96
8
21
6.
27
7
11
9.
84
6
F
ig
ur
e
6.
6(
c)
W
or
st
f 3
0.
43
1
0.
24
9
9.
75
4
3.
77
2
20
1.
14
2
43
.6
20
F
ig
ur
e
6.
6(
d)
L
2
-n
or
m
0.
23
6
0.
28
4
5.
39
8
3.
69
5
20
6.
55
8
69
.6
84
F
ig
ur
e
6.
7
132
6.2.1.2 Blunted Case
Effects of the blunted leading edge for the ABLT case parallel that of the
baseline case. The analysis of the blunt leading edge results in the baseline case
should be consulted to understand the results in this section. The analysis of f3 for
the blunted ABLT influenced case is analogous to the discussion for the infinitely
sharp case. The ABLT infinitely sharp case section can be consulted to understand
the f3 results in this section. The tables and figures below are provided for reference.
133
(a) Best f1 (b) Best f2 (c) Best f3 (d) Worst f3
(e) Best f4
Figure 6.13: Extreme Pareto-optimal point vehicle designs with best f1, f2, f3, f4
values and the design with the worst f3 value for the ABLT blunted leading edge
case (isometric views).
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.14: L2-norm Pareto-optimal point vehicle design for the ABLT blunted
leading edge case.
134
T
ab
le
6.
7:
P
ar
et
o-
op
ti
m
al
p
oi
nt
ob
je
ct
iv
e
fu
nc
ti
on
va
lu
es
fo
r
th
e
A
B
LT
bl
un
te
d
le
ad
in
g
ed
ge
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
L
if
t
D
ra
g
?
f 3
m
ax
(q
w
)
(W
/c
m
2
)
F
ig
ur
e
B
es
t
f 1
7.
78
3
0.
20
2
1.
25
7
?
10
?
7
39
3.
37
7
F
ig
ur
e
6.
13
(a
)
B
es
t
f 2
1.
64
3
0.
70
3
1.
19
9
?
10
?
6
39
.3
69
F
ig
ur
e
6.
13
(b
)
B
es
t
f 3
5.
99
1
0.
20
1
9.
74
8
?
10
?
8
16
9.
59
7
F
ig
ur
e
6.
13
(c
)
W
or
st
f 3
3.
40
7
0.
62
6
1.
78
7
?
10
?
6
80
.1
23
F
ig
ur
e
6.
13
(d
)
B
es
t
f 4
1.
29
0
0.
69
1
9.
23
5
?
10
?
7
39
.3
40
F
ig
ur
e
6.
13
(e
)
L
2
-n
or
m
4.
22
4
0.
48
8
6.
84
7
?
10
?
7
88
.0
30
F
ig
ur
e
6.
14
T
ab
le
6.
8:
P
ar
et
o-
op
ti
m
al
p
oi
nt
de
si
gn
ve
ct
or
va
lu
es
fo
r
th
e
A
B
LT
bl
un
te
d
le
ad
in
g
ed
ge
ca
se
.
P
ar
et
o-
op
ti
m
al
p
oi
nt
n
m
?
(d
eg
)
?
(d
eg
)
L
(f
t)
w
(f
t)
d
(f
t)
F
ig
ur
e
B
es
t
f 1
0.
00
5
0.
11
0
2.
24
9
3.
59
9
19
0.
75
9
11
9.
68
0
0.
01
7
F
ig
ur
e
6.
13
(a
)
B
es
t
f 2
1.
00
0
0.
16
4
10
.0
00
2.
58
1
16
0.
48
4
55
.1
60
1.
64
0
F
ig
ur
e
6.
13
(b
)
B
es
t
f 3
0.
01
0
0.
18
6
1.
88
5
4.
37
6
18
9.
78
5
11
8.
11
9
0.
08
9
F
ig
ur
e
6.
13
(c
)
W
or
st
f 3
0.
56
3
0.
16
4
9.
11
8
4.
00
4
20
7.
78
8
42
.2
60
0.
39
9
F
ig
ur
e
6.
13
(d
)
B
es
t
f 4
1.
00
0
0.
08
3
9.
99
9
2.
28
4
15
0.
76
1
69
.7
20
1.
64
0
F
ig
ur
e
6.
13
(e
)
L
2
-n
or
m
0.
50
4
0.
44
5
7.
07
4
3.
46
8
18
8.
60
1
77
.4
53
0.
33
1
F
ig
ur
e
6.
14
135
6.2.2 Robust Optimization Results
Section 6.2.2 provides all the results from the four robust optimization prob-
lems. Generation histories and optimum design vectors are provided and discussed
from an optimization perspective. Section 6.2.2.5 draws comparisons between the
robust optimum designs, resulting hypersonic vehicle shapes, and the corresponding
method effectiveness and expense. Section 6.2.2.6 provides an example of how the
robust MOGA method has the unique capability to provide the designer with more
information about the design space.
6.2.2.1 Robust Regularization Results
Figure 6.15 provides a generation history plot for the robust regularization-
based method. The figure plots the minimum FI value found in the population
at a given generation number. This value decreases from generation to generation,
eventually leveling off. After about generation 250, little to no improvement in the
objective function is made.
It is reasonable to say (without any strict convergence criterion) that the
optimization problem converges at generation 250. The optimum design vector ~x?
is the value of ~x which gives the minimum value of FI at generation 250. The value
of ~x? for the robust regularization-based method was the following:
~x? =
[
0.05 0.999 4.414 4.015 146.424 120.0
]
. (6.2)
136
For the implementation defined in Equation 5.25, the optimum result was found
after 137,500 function evaluations.
0 100 200 300 400 500 600 700 800 900 10000
0.05
0.1
0.15
0.2
0.25
Generation
mini
mum
F I(~x
)
Figure 6.15: Robust regularization-based method generation history.
6.2.2.2 Explicit Averaging Results
Figure 6.16 provides a generation history plot for the explicit averaging method
using Monte Carlo integration. The figure plots the minimum FII value found in
the population at a given generation number. Due to the nature of the sampling
technique in this method, the plot is noisy; FII doesn?t monotonically decrease with
increasing generation number. Nonetheless, Figure 6.16 shows a rapid improvement
in the objective function value which quickly levels off at about generation 100 . The
value continues to decrease until somewhere between generation 200 and 300. After
this, improvement in the objective function value is indistinguishable through the
noise. It will be assumed that the optimum design vector was found at generation
250.
137
The optimum design vector ~x? is the value of ~x which gives the minimum value
of FII at generation 250. The value of ~x? for the explicit averaging method was
~x? =
[
0.999 0.903 8.559 3.16 176.631 120.0
]
. (6.3)
For the implementation defined in Equation 5.26, the optimum result was found
after 125,000 function evaluations.
0 100 200 300 400 500 600 700 800 900 1000
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Generation
mini
mum
F II(
~x,p)
(a) Explicit averaging generation history.
0 100 200 300 400 500 600 700 800 900 10000
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Generation
mini
mum
F II(
~x,p)
(b) Explicit averaging generation history
(zoomed in).
Figure 6.16: Explicit averaging generation history.
6.2.2.3 Perturbation Results
The nature of the perturbation method makes generation histories based on
FIII difficult to interpret. The value of p is generation dependent, meaning the plot
would be very noisy and may potentially mask any convergence trends. Tsutsui
et al. [72, 73] determined convergence by plotting the mean value of the design
variable(s) over the population at each generation. Figure 6.17 shows these trends
for the optimization problem in this study. Each design variable fluctuates with
138
different magnitudes and frequencies as the generations progress. Eventually, each
of the design variables finds a nearly constant value and no longer changes with
subsequent generations. Once all of the variables reach this condition, the problem
has converged to a solution. From Figure 6.17, convergence was found to occur at
about generation 750.
The optimum design vector ~x? was taken to be the vector ~x in the population
which has the smallest departure from the mean value at generation 750. The value
of ~x? for the perturbation method was
~x? =
[
0.999 0.999 5.578 1.671 214.01 119.934
]
. (6.4)
For the implementation defined in Equation 5.27, the optimum result was found
after 37,500 function evaluations.
139
0 100 200 300 400 500 600 700 800 900 10000.5
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
Generation
mea
nn
(a) Mean value of n
0 100 200 300 400 500 600 700 800 900 10000.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Generation
mea
nm
(b) Mean value of m
0 100 200 300 400 500 600 700 800 900 10004.5
5
5.5
6
Generation
mea
n?(
deg)
(c) Mean value of ?
0 100 200 300 400 500 600 700 800 900 10000
2
4
6
8
10
12
Generation
mea
n?(
deg)
(d) Mean value of ?
0 100 200 300 400 500 600 700 800 900 1000100
120
140
160
180
200
220
Generation
mean
L(ft
)
(e) Mean value of L
0 100 200 300 400 500 600 700 800 900 100060
70
80
90
100
110
120
Generation
mean
w(ft
)
(f) Mean value of w
Figure 6.17: Generation history of the mean ~x values over the population for the
perturbation method.
140
6.2.2.4 Robust MOGA Results
As discussed in Section 4.2.3.2, convergence of the robust MOGA optimization
method is measured using Equation 4.10. Figure 6.18 plots the average ?PM over
the previous 10 generations as a function of current generation number. Starting at
generation 10, the average value of ?PM is about 52. As generations progress, the
average value of ?PM decreases rapidly. With increasing generation number the
average value of ?PM levels off, between generation 100 and 500, to a value which
fluctuates between 0.0 and 0.5. Over the next 1500 generations, these fluctuations
slowly die out. Typical values of ?PM which correspond to a particular level of
accuracy in the true Pareto-optimal set approximation is not fully understood, yet
a preliminary analysis on the convergence trends and its link to accuracy, quality,
and speed of the optimization process is given in Section 4.4. Overall, a converging
trend is observed in Figure 6.18.
0 500 1000 1500 2000
0
10
20
30
40
50
60
Generation (t)
Mov
ing?
PM
aver
age
over
10g
ens.
(a) Full generation history.
0 500 1000 1500 2000
0
0.5
1
1.5
Generation (t)
Mov
ing?
PM
aver
age
over
10g
ens.
(b) Generation history zoomed in.
Figure 6.18: Robust MOGA mean ?PM generation history.
141
In attempt to deliver results from a simulation which is as converged as pos-
sible, the true Pareto-optimal solution set will be taken as the approximation at
generation 2000 (x?2000). The set of Pareto-optimal solutions at generation 2000
contains 96 individuals in the population. For each individual in the Pareto-optimal
set, the mean value of FIV1 through FIV5 was calculated over the 2000 generations
(i.e., ?FIVm (~x)?t seen in Equation 4.11). The mean objective function values, when
plotted in objective function space, form the Pareto-optimal front.
Figure 6.19 shows the value of ?FIV2 (~x)?t vs. the value of ?FIV1 (~x)?t for each
of the 96 Pareto-optimal solutions. From this plot, it is seen that these two objective
functions are non-competing. Therefore, only one of these two functions needs to
be plotted in order to form the Pareto-optimal front.
0 0.002 0.004 0.006 0.008 0.01 0.0120
0.005
0.01
0.015
0.02
0.025
0.03
mean FIV1 over all gens.
mea
nF
IV 2
ove
ral
lge
ns.
Figure 6.19: Average FIV2 vs. average FIV1 over all generations for each of the 96
Pareto-optimal solutions.
142
The Pareto-optimal front in four dimensions is shown in Figure 6.20. The x-,
y-, and z-axis in the plot are represented by ?FIV1 (~x)?t, ?FIV3 (~x)?t, and ?FIV4 (~x)?t,
respectively. The value of ?FIV5 (~x)?t is represented by color. To gain a better
feel and understanding for the four-dimensional Pareto-optimal hypersurface, three-
dimensional projections of Figure 6.20 is provided in Figure 6.21.
0.0050
0.0100
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.025
0.03
0.035
0.04
0.045
0.05
0.055
 
mean FIV1over all gens.
mean FIV3over all gens.
 
me
an
F I
V 4
ov
er
all
ge
ns
.
me
an
F I
V 5
ov
er
all
ge
ns
.
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0.11
Figure 6.20: Four-dimensional Pareto-optimal front.
Robust regularization and explicit averaging equivalent solutions can be picked
out of the set of Pareto-optimal solutions using Equation 4.11 and Equation 4.13,
respectively. The resulting optimum design vectors using these equations were the
following:
~x?rr =
[
0.113 0.793 3.688 1.656 184.8 119.314
]
, (6.5)
~x?eap =
[
0.989 0.351 5.596 2.421 188.763 119.875
]
. (6.6)
143
The robust MOGA formulation using the prescribed MOGA population settings
required 500,000 function evaluations after 2000 generations and resulted in 96 dif-
ferent optimum design vectors.
Note, this is not necessarily an ?apples to apple? comparison, when comparing
the number of total function evaluations for each method. The implementation and
convergence criterion greatly affects this value. The main purpose of Section 6.2.2
was to compare the resulting robust optimum vehicle designs from each method.
The given total number of function evaluations for each method here is given in
an attempt to be concise. The efficiency and convergence aspects of the proposed
robust MOGA method is given in Section 4.4.
144
0
0.0
05
0.0
1
0.0
15
0.0
2
0.0
250.0
3
0.0
350.0
4
0.0
450.0
5
0.0
550.0
6
me
an
F IV
1
ove
ra
llg
ens
.
meanFIV4 overallgens.
0
0.0
1
0.0
2
0.0
3
0.0
4
0.0
2
0.0
250.0
3
0.0
350.0
4
0.0
450.0
5
0.0
550.0
6
me
an
F IV
3
ove
ra
llg
ens
.
meanFIV4 overallgens.
0
0.0
1
0.0
2
0.0
3
0.0
4
0
0.0
02
0.0
04
0.0
06
0.0
080.0
1
0.0
12
me
an
F IV
3
ove
ra
llg
ens
.
meanFIV1 overallgens.
F
ig
ur
e
6.
21
:
T
hr
ee
-d
im
en
si
on
al
pr
oj
ec
ti
on
s
of
th
e
P
ar
et
o-
op
ti
m
al
fr
on
t.
145
6.2.2.5 Result Comparisons and Vehicle Analysis
This section will aim at comparing the results from the different robust op-
timization methods. A focus will be set on analyzing the effectiveness of using a
robust MOGA method coupled with post-optimality data handling techniques to
pick preferred solutions from a set of Pareto-optimal solutions that match the phi-
losophy of the other robust methods.
Table 6.9 provides a list of the robust optimum design vectors found by each
method. The first two entries (above the dashed line) represent the results based
on the worst-case scenario philosophy. The next three (below the dashed line)
represent the results based off of the philosophy of minimizing the expected value
of the original objective function over the uncertainty set. To aid in the analysis of
the robust optimum design vectors, the robust optimum vehicle shapes are provided
in Figure 6.22 through Figure 6.26. A few trends were common to all the robust
optimum designs. The width parameter, w, of the robust optimum geometry for
all cases is very close to the upper design limit of 120 feet. Although not shown in
the table, the internal volume of each robust optimum design is near the lower limit
constraint of 88,286 cubic feet. Also, all the shapes are convex, meaning ? > ?.
146
T
ab
le
6.
9:
R
ob
us
t
op
ti
m
um
de
si
gn
s.
A
lg
or
it
hm
~x
?
n
m
?(
de
g)
?(
de
g)
L
(f
t)
w
(f
t)
R
ob
us
t
re
gu
la
ri
za
ti
on
(E
qu
at
io
n
5.
25
)
0.
05
0.
99
9
4.
41
4
4.
01
5
14
6.
42
4
12
0.
0
M
O
G
A
:
ro
bu
st
re
g.
eq
ui
va
le
nt
(E
qu
at
io
ns
4.
30
,4
.1
3)
0.
11
3
0.
79
3
3.
68
8
1.
65
6
18
4.
8
11
9.
31
4
A
gg
re
ga
ti
on
:
ex
pl
ic
it
av
er
ag
in
g
(E
qu
at
io
n
5.
26
)
0.
99
9
0.
90
3
8.
55
9
3.
16
17
6.
63
1
12
0.
0
A
gg
re
ga
ti
on
:
p
er
tu
rb
at
io
n
(E
qu
at
io
n
5.
27
)
0.
99
9
0.
99
9
5.
57
8
1.
67
1
21
4.
01
11
9.
93
4
M
O
G
A
:
ag
gr
eg
at
io
n
eq
ui
va
le
nt
(E
qu
at
io
ns
4.
30
,4
.1
1)
0.
98
9
0.
35
1
5.
59
6
2.
42
1
18
8.
76
3
11
9.
87
5
147
The two solutions based on a worst-case scenario are characterized by small
n values near or at the lower limit of 0.05 and large m values near or at the upper
limit of 1.0. The two shapes are convex, but the values of ? are quite close to the
value of ?. The optimal shape for the robust regularization-based method has a ?
value that is only 0.4 degrees greater than ?. For the MOGA robust regularization
equivalent solution, ? is about two degrees greater than ?. Figure 6.22 and Figure
6.23 show that these design vector qualities produce a flat, spatula or wedge-like
shape.
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.22: Optimal vehicle design for the robust regularization-based method.
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.23: Optimal vehicle design for the robust MOGA method with robust
regularization equivalent data handling.
The three minimized expected value robust optimum solutions are character-
ized by large values of n very close to the design variable upper limit of 1.0. This
148
produces a triangular planform shape. The explicit averaging and perturbation
method optimum designs have large m values near the limit of 1.0. This produces
a diamond shape cross section as seen in Figure 6.24 and Figure 6.25. The MOGA
explicit averaging and perturbation equivalent solution has an optimum m value of
0.351. This value isn?t quite are large as the other two. Yet judging from Figure
6.26, qualitatively the overall shape remains the same (namely a pyramidal shape).
To gain a better understanding of the implications of this difference, a more quan-
titative look at the corresponding f values is needed.
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.24: Optimal vehicle design for the explicit averaging method.
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.25: Optimal vehicle design for the perturbation method.
In order to quantitatively analyze each optimal design, values of f as a func-
tion of p are provided in Figure 6.27. The curve marked with square data points
corresponds to the robust regularization-based method optimum design. The curve
149
(a) Isometric view (b) Planform view (c) Rear view (d) Side view
Figure 6.26: Optimal vehicle design for the robust MOGA method with explicit
averaging and perturbation equivalent data handling.
marked with downward pointing triangles is the MOGA robust regularization equiv-
alent solution. These solutions are optimum based on a worst-case scenario philos-
ophy. Therefore, they are the designs for which the max value of f (~x, p) over all p
values is minimized. This can be seen in the plot, where both maximum values of
f on the curves are around 0.04 to 0.05. The trends of both curves exhibit similar
characteristics. They are generally flat from p = 0.5 down to some value below
p = 0.25. The curves then dip down to a value close to zero at about p = 0.05.
Between p = 0.0 and p = 0.05 both curves exhibit a bump. The bump in the curve
for the robust regularization-based method is significantly larger than the bump in
the robust MOGA curve. This difference is due to the formulation of the robust
regularization optimization problem in Equation 5.25, where p is sampled from a
discrete set p? that only includes values 0.0, 0.05, 0.1, and so on. Therefore, p values
between 0.0 and 0.05 do not get evaluated and are not taken into account when
searching for an optimum design. The MOGA formulation however, accounts for p
over the entire range 0.0 ? p ? 0.5. The two curves show us that the MOGA-based
answer is therefore more successful in minimizing the maximum value of f .
150
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
p
f(~x
,p)
 
 ~x = ~x?I~x = ~x?I I~x = ~x?I I I~x = ~x?eap~x = ~x?rr
Figure 6.27: Original objective function values of each optimum design vector as a
function of environmental parameter p.
The other three designs are represented by the curves which are marked with
the star, upward pointing triangle, and circle data points. These methods are based
off the philosophy of minimizing the expected value of f over uncertainty in p, which
for this study was set to include all values in the range 0.0 ? p ? 0.5. The optimum
designs aim to minimize the area under the curves presented in Figure 6.27. The
three curves for the methods based on an optimum expected value exhibit very
similar characteristics. The overall trends are all the same. The area under the
curve between p = 0.0 and p = 0.25 is almost reduced to zero. From p = 0.25 to
p = 0.5 the curve increases monotonically from f = 0.0 to a value of f ranging
between 0.12 and 0.17. The curve corresponding to the robust MOGA method after
applying the data handling falls somewhere between the other two methods. From
151
this, it can be concluded that the slightly different value of m in the optimum design
vector does not affect the overall goal of minimizing the expected value.
The vehicle designs that are wedge-like in shape are the designs that minimize
the maximum change in moment coefficient due to an asymmetric boundary-layer
transition occurring anywhere on the leading edge. The vehicle accomplishes this
by distributing its internal volume equally along the spanwise direction, namely
in a wedge shape. This means, an asymmetric viscous drag, occurring on the top
surface the vehicle would contribute mostly to a change in yaw moment. The change
in pitching moment is minimized by having a small wedge angle, bringing the top
surface close to the center of gravity. From the analysis in Section 6.2.1, it was
shown that changes in pitching moment are the main contributors to changes in f .
Therefore, a wedge-like shape produces relatively constant value of f over different
values of p.
The vehicle designs that are pyramidal in shape are the designs that minimize
the expected change in moment coefficient due to an asymmetric boundary-layer
transition, if the event is equally probable to happen at any spanwise location along
the leading edge. The vehicle has a strongly swept leading edge. This renders
a minimal area covered by the turbulence wedge for small p values. Also, the
diamond shape cross section acts to minimize the moment arm with respect to the
center of gravity for small p values. The internal volume of the vehicle is focused
more towards the centerline of the body creating more distance between the CG and
the top surface for high values of p. This means, if the isolated roughness event takes
place near the centerline, a large change in pitching moment will occur rendering a
152
large value of f . However, in this study, the robust optimization formulations based
on conditional expectation state that any spanwise location is equally probable,
therefore, the expected value of f over all p is still minimized.
An inherent trade-off exists between these two shapes. The wedge-like shape
has a worst-case scenario f value of about 0.04 and a half-span equally probable
expected f value of about 0.0375. The pyramidal shape can improve on this half-
span equally probable expected value by about 12%, but at the cost of increasing
the worst-case scenario f value by 350%.
6.2.2.6 Design Space Analysis
The optimization results reveal that the robust MOGA method has the ca-
pability to identify the optimum designs which follow the same philosophies as the
robust regularization- and aggregation-based methods. The robust MOGA method
accomplished this task by using post-optimality data handling to pick preferred so-
lutions from the set of 96 designs in the Pareto-optimal front. The Pareto-optimal
front exists as an organized set of designs which inherently illustrates design trends
or trade-offs within the design space. The robust MOGA method therefore en-
ables the designer to investigate the trade-offs within the design space by analyzing
the Pareto-optimal front. Single-objective methods like the robust regularization-
and aggregation-based methods cannot provide such information without repeat-
edly reforming and solving the optimization problem. This tends to get tedious
and computationally expensive. This section will provide an example of how the
153
designer could utilize the information provided by the robust MOGA method results
to investigate design trends or trade-offs.
The comparisons of the results in Section 6.2.2.5, specifically discussions on
Figure 6.27, show that the results based on robust regularization and aggregation
trade off between minimizing the mean values of objective functions FIV5 and FIV1 ,
respectively. The previous discussion linked this trade-off to the optimum design?s
allocation of internal volume and the corresponding effect on f . To further in-
vestigate this trade-off, the Pareto-optimal front is presented as a projection in
?FIV5?t-?FIV1?t space. This plot is given in Figure 6.28(a). To isolate the trade-off,
only the non-dominated points in this space are examined. Figure 6.28(b) removes
all of the dominated points, leaving 23 non-dominated points which illustrate the
trade-off between ?FIV5?t and ?FIV1?t.
Eight designs are picked within this set to aid in the discussion and are labeled
A through H in the plot. These points were chosen in attempt to equally sample
the Pareto-optimal curve. The vehicle geometries of the eight designs are given in
Figure 6.29. From a preliminary qualitative look at the designs and a quantitative
look at the corresponding mean FIV5 and FIV1 values, the trade-off between robust
regularization optimum designs and aggregation optimum designs is apparent. This
information provides the designer with a better understanding of how the optimum
designs presented in Section 6.2.2.5 are coupled. Designs A through H also provide
the designer with options of deigns which balance worst-case scenario characteristics
and expected value characteristics (balancing the robust regularization philosophy
and the aggregation philosophy, respectively).
154
0 0.002 0.004 0.006 0.008 0.01 0.0120.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0.11
0.12
mean FIV1 over all gens.
mea
nF
IV 5
over
allg
ens.
(a) FIV5 vs FIV1 Pareto-optimal projection
0 0.002 0.004 0.006 0.008 0.01 0.0120.03
0.04
0.05
0.06
0.07
0.08
0.09
mean FIV1 over all gens.
mea
nF
IV 5
over
allg
ens.
A B
C
D
E
F
G
H
(b) Non-dominated designs in FIV5 -FIV1 space
Figure 6.28: Two-dimensional projections of Pareto-optimal front.
155
(a) Design A (b) Design B (c) Design C
(d) Design D (e) Design E (f) Design F
(g) Design G (h) Design H
Figure 6.29: Vehicle designs A through F from Figure 6.28(b).
156
Chapter 7: Conclusions
This research provides a simplified methodology for modeling boundary-layer
effects of isolated surface roughnesses for hypersonic vehicle applications without
the use of expensive boundary-layer or CFD solvers. The effect of induced moments
on vehicle controllability for a particular vehicle can be quantified, determining the
amount of control authority necessary for vehicle survival. Inexpensive calculations
of vehicle asymmetric aerodynamics made it feasible to apply multi-objective and
robust optimization schemes, thus providing an analysis of the design trade-offs
between the performance, controllability, and usability of hypersonic systems. This
design space was further explored through the use of a novel robust optimization
method. It was found that this method is more versatile and has proved to be more
informative than traditional robust optimization methods.
The results from the multi-objective optimization studies identified two classes
of shapes on the boundaries of the Pareto-optimal front: wedge-like and cone-like.
Objective function values of the optimum designs were found to compete between
these two classes of shapes. When analyzing the effect of ABLT, it was found that
the cone-like shape had 3.5 times the volumetric efficiency of the wedge-like shape,
yet compromised with an approximate 75% decrease in lift-to-drag ratio and one
157
full order of magnitude jump in the controllability factor when compared to the
wedge-like counterpart. A design of minimum compromise (the L2-norm optimum)
was found to be a hybrid shape between a wedge of small deflection angle and a
cone. The L2-norm design had approximately 2.3 times the volumetric efficiency of
the wedge-like shape with a 25% reduction in lift-to-drag ratio and about five times
(half an order of magnitude) the moment coefficient magnitude of the wedge-like
shape.
The results from this study identified important design considerations when
looking at the effects of an asymmetric boundary-layer. It was found that the pitch-
ing moment produced by the spanwise turbulence wedge is the design driver for
controllability. The results from this study show that once the pitching moments
are trimmed, the turbulence wedge induced yawing moments are no longer domi-
nated by the induced pitching moments. Therefore both moment components need
to be quantified for adequate design of a vehicle?s guidance and control system.
Strong adverse pitching moment characteristics are often augmented by longitudi-
nal movement of the CG. The design point specification of the longitudinal CG
position is important and plays a significant role in the resulting optimum shapes.
Studies on how the CG placement affects the Pareto-optimal designs is suggested
as future work.
The aerodynamic model and optimization methodology for an infinitely sharp
waverider with a fully laminar boundary-layer was verified with the results given by
Starkey [30]. Adding the capability of a blunted leading edge to the design space
did not drastically change the resulting optimum vehicle designs. Incorporating an
158
objective of minimizing total heat load in conjunction with minimizing total heat
flux would be worth investigating as future work.
A novel multi-objective robust counterpart formulation to solve single-objective
optimization problems with environmental parameter uncertainty is proposed. The
method uses post-optimality data handling techniques to identify, among a set of
Pareto-optimal solutions, the optimal worst-case scenario and expected value solu-
tions for any uncertainty PDF. Solving the robust multi-objective formulation once
produces a set of designs, thus giving the designer more information about the design
space. This robust multi-objective method was applied to an example optimization
problem, where the predicted robust optimum results were compared to the true
robust optimum results derived analytically. These results were used as a baseline
to study how changes in the method initial conditions and data handling settings
affect the accuracy and convergence qualities of the method.
The baseline results show that the post-optimality data handling techniques
were successful in identifying optimum worst-case scenario and expected value solu-
tions for a wide range of uncertainty conditions. The worst-case scenario optimum
design variable values predicted by the robust multi-objective method deviated from
the analytical values by no more than 20% of the design variable range. Most of
the predicted design variables deviated from the analytical values by 3%-7% of the
design variable range. The robust multi-objective method predicted optimum ex-
pected value designs with similar accuracy for most uncertainty PDF mean and
variance values. Error in the predictions became significant if the uncertainty PDF
variance was not large enough to adequately cover multiple uncertainty subsets.
159
The accuracy and convergence tests revealed that a dynamic relationship ex-
ists between the initial conditions and data handling settings of the method and the
accuracy, quality, and speed of the results. The number of histogram bins had the
strongest effect on the speed and quality of the results. Increasing the number of
bins from 5 to 60 increased the generation number at convergence from 13 to 939
and increased the number of Pareto-optimal points from 26 to 2214. The number of
uncertainty subsets had the strongest effect on the accuracy of the results. Increas-
ing the number of uncertainty subsets from 3 to 8 decreased the mean error from
21.5% to 5.9% of the design variable range and the RMS error from 14.6% to 3.8%
of the design variable range. The effects of changing the population size, number of
bins, and number of subsets were found not to be completely independent. There-
fore, a further investigation of the coupled relationships through a DOE analysis is
suggested as future work.
The robust multi-objective optimization formulation was also applied the hy-
personic ABLT vehicle design problem. The results from the robust multi-objective
method and three traditional robust optimization methods, when applied to the
design problem, were compared and analyzed. The analysis showed that the novel
robust multi-objective optimization formulation coupled with post-optimality data
handling was capable of finding results that would be expected from typical robust
optimization methods. For changes in the environmental parameter uncertainty
probability density function, the traditional robust optimization methods would
need to be solved again. This is because these methods only produce one optimal
design for each PDF. However, when the robust multi-objective method optimiza-
160
tion problem is solved, a set of solutions is produced. Therefore, for any uncertainty
PDF, the optimum solution can be identified using a post-optimality data handling
equation. This not only reduced computational expense, but allows for the designer
to investigate design trends which may exist for the problem at hand. This was
exemplified in the design space analysis in Section 6.2.2.6.
The hypersonic vehicle design problem provided a valuable test-bed for the
robust optimization methods investigated in this study. It was found that differ-
ent robust design philosophies produce significantly different shaped vehicles that
exhibit different responses to isolated roughness-induced boundary-layer transition.
Vehicles that are spatula or wedge shaped are ideal for a worst-case scenario. Highly
swept pyramidal shaped vehicles minimize the aerodynamic effects if it is equally
probable that the roughness is present anywhere along the leading edge. An inherent
trade-off was found to exist between these two shapes. The wedge-like shape had a
worst-case scenario f value of about 0.04 and a half-span equally probable expected
f value of about 0.0375. The pyramidal shape improved upon this half-span equally
probable expected value by about 12%, but at the cost of increasing the worst-case
scenario f value by 350%.
The conclusions from this work can be summarized by the following:
? Multi-objective optimization studies identified two categories of shapes which
bound the Pareto-optimal front for the trade-off between vehicle performance,
controllability, and usability. These shape categories are wedge-like and cone-
like.
161
? Robust optimization studies identified two categories of shapes (one additional
to the categories found in the multi-objective studies) which are optimal for
combating the unwanted effects on controllability caused by roughness-induced
ABLT. Spatula or wedge shaped vehicles are ideal when taking a worst-case
scenario approach. Highly swept pyramidal shaped vehicles minimize the aero-
dynamic effects if it is equally probable that the roughness is present anywhere
along the leading edge.
? This research revealed a missing capability in state-of-the-art robust optimiza-
tion methods, especially when applied to engineering problems. This work
identified the need for an all encompassing robust optimization method that
could provide multiple types of robust optima for many uncertainty distribu-
tions without having to repeatedly resolve the problem. The current study
fulfilled this need by providing a novel robust multi-objective formulation.
? The proposed robust multi-objective optimization method to solve single-
objective optimization problems with parameter uncertainty via the use of
post-optimality data handling techniques was successful in finding the com-
mon results found by traditional methods. Moreover, the method solves for
a coupled set of robust optimum designs. This provides the designer with an
opportunity to examine the trade-offs in the design space at no extra cost. The
ability to investigate the design space and compare multiple robust optimum
designs is a particularly important contribution to the field of engineering op-
timization, where providing context with an optimum design is often valuable.
162
? The proposed novel robust optimization method provided an organized Pareto-
optimal front within the ABLT vehicle design space. The results showed that
a wedge shape and pyramidal shape vehicle bound a trade-off between induced
moments in the pitch plane and yaw plane.
163
Bibliography
[1] E. Reshotko. Transition issues for atmospheric entry. Journal of Spacecraft
and Rockets, 45(2):161?164, 2008.
[2] T. J. Horvath, J. N. Zalameda, W. A. Wood, S. A. Berry, R. J. Schwartz,
R. F. Dantowitz, T. S. Spisz, and J. C. Taylor. Global infrared observations
of roughness induced transition on the space shuttle orbiter. RTO AVT-200
RSM-030 Specialists? Meeting on Hypersonic Laminar-Turbulent Transition,
April 2012.
[3] Press release: Engineering review board concludes re-
view of htv-2 second test flight. DARPA, April 2012.
http://www.darpa.mil/NewsEvents/Releases/2012/04/20.aspx.
[4] H. Beyer and B. Sendhoff. Robust optimization - a comprehensive survey.
Computer Methods in Applied Mechanics and Engineering, 196(33-34):3190?
3218, 2007.
[5] Y. Ong, P. B. Nair, and K. Y. Lum. Max-min surrogate-assisted evolutionary
algorithm for robust design. Evolutionary Computation, IEEE Transactions
on, 10(4):392 ?404, August 2006.
[6] G. Ridolfi, E. Mooij, D. Dirkx, and S. Corpino. Robust multi-disciplinary
optimization of unmanned entry capsules. AIAA Modeling and Simulation
Technologies Conference, August 2012.
[7] A. Lewis. Robust regularization. Technical report, Simon Fraser University,
September 2002.
[8] Y. Jin and J. Branke. Evolutionary optimization in uncertain environments-a
survey. Evolutionary Computation, IEEE Transactions on, 9(3):303?317, June
2005.
[9] J. W. Kruisselbrink. Evolution Strategies for Robust Optimization. PhD thesis,
Leiden University, May 2012.
164
[10] K. M. Ryan and M. J. Lewis. Multi-objective optimization of hypersonic
vehicles under asymmetric boundary-layer conditions. 18th AIAA/3AF Inter-
national Space Planes and Hypersonic Systems and Technologies Conference,
September 2012.
[11] K. M. Ryan, M. J. Lewis, and K. H. Yu. Comparison of robust optimization
methods applied to hypersonic vehicle design. AIAA Modeling and Simulation
Technologies (MST) Conference, August 2013.
[12] K. M. Ryan, M. J. Lewis, and K. H. Yu. A multi-objective post-optimality
data handling approach to robust optimization. 10th AIAA Multidisciplinary
Design Optimization Conference, January 2014.
[13] M. Adler, M. Wright, C. Campbell, I. Clark, W. Engelund, and T. Rivellini.
Entry, descent, and landing roadmap: Technology area 09. Technical report,
National Aeronautics and Space Administration, April 2012.
[14] J. Knittel, M. J. Lewis, and K. H. Yu. Optimization of a martian aero-gravity
assist. AIAA Atmospheric Flight Mechanics Conference, 2014.
[15] T. R. F. Nonweiler. Aerodynamic problems of manned space vehicles. Journal
of the Royal Aeronautical Society, 63:521?528, September 1959.
[16] J. G. Jones, K. C. Moore, J. Pike, and P. L. Roe. A method for designing lift-
ing configurations for high supersonic speeds, using axisymmetric flow fields.
Ingenieur-Archiv, 37(1):56?72, 1968.
[17] K. B. Center, H. Sobieczky, and F. C. Dougherty. Interactive design of hyper-
sonic waverider geometries. AIAA Paper, June 1991.
[18] R. P. Starkey and M. J. Lewis. A simple analytical model for parametric stud-
ies of hypersonic waveriders. AIAA International Space Planes and Hypersonic
Systems and Technologies Conference, April 1998.
[19] J. D. Anderson. Hypersonic and High-Temperature Gas Dynamics. American
Institute of Aeronautics and Astronautics, second edition, 2006.
[20] L. V. Rudd, D. J. Pines, and C. Tarpley. Dynamic stability of mission-oriented
hypersonic waveriders. AIAA, Aerospace Sciences Meeting and Exhibit, 37th,
Reno, NV, 1999.
[21] L. Rudd and D. J. Pines. Dynamic control of mission oriented hypersonic
waveriders. AIAA Paper, 1999.
[22] C. Tarpley and M. J. Lewis. Stability derivatives for engine integrated wa-
veriders with viscous and pitch effects. AIAA, Aerospace Sciences Meeting
and Exhibit, 32nd, Reno, NV, 1994.
165
[23] D. B. Johnson, R. Thomas, and D. Manor. Stability and control analysis of
a wave-rider tsto second stage. AIAA 10th International Space Planes and
Hypersonic Systems and Technologies Conference, 2001.
[24] M. W. Oppenheimer and D. B. Doman. A hypersonic vehicle model developed
with piston theory. AIAA Atmospheric Flight Mechanics Conference, 2006.
[25] D. D. Liu, Z. X. Yao, D. Sarhaddi, and F. Chavez. From piston theory to
a unified hypersonic-supersonic lifting surface method. Journal of Aircraft,
34(3):304?312, 1997.
[26] M. L. Rasmussen. Stability derivatives for hypersonic waveriders according to
newtonian theory. AIAA Paper, 1995.
[27] M. L. Rasmussen. Effects of angle-of-attack and sideslip on the forces and
moments of hypersonic waveriders. AIAA Atmospheric Flight Mechanics Con-
ference, 1996.
[28] M. L. Rasmussen. Effects of anhedral and finlets on lateral stability of hyper-
sonic waveriders. AIAA Aerospace Sciences Meeting & Exhibit, 1997.
[29] M. L. Rasmussen and B. Brandes-Duncan. Some general features of waveriders
derived from power-law shocks. AIAA Paper, 1996.
[30] R. P. Starkey. A parametric study of l/d and volumetric efficiency trade-offs
for waverider based vehicles. Master?s thesis, University of Maryland, 1998.
[31] S. Corda and J. D. Anderson. Viscous optimized hypersonic waveriders de-
signed from axisymmetric flow fields. AIAA Aerospace Sciences Meeting, 1988.
[32] T. F. O?Brien and M. J. Lewis. Power law leading edges for waveriders designed
with shock attachment. AIAA Paper 98-600, 1998.
[33] C. E. Cockrell, L. D. Huebner, and D. B. Finley. Aerodynamic characteristics
of two waverider-derived hypersonic cruise configurations. NASA Technical
Paper 3559, 1996.
[34] D. B. Finley and C. Cockrell. Control effectiveness and lateral-directional
stability for two wave-rider derived hypersonic cruise configurations. AIAA
Applied Aerodynamics Conference, 1995.
[35] M. J. Gillum and M. J. Lewis. Experimental results on a mach 14 waverider
with blunt leading edges. Journal of Aircraft, 34(3):296?303, 1997.
[36] M. L. Blosser, I. M. Blankson, S. Schoerke, D. Mrunson, and P. Hagseth. Wing
leading-edge design concepts for airbreathing hypersonic waveriders. 32nd
Aerospace Sciences Meeting and Exhibit, January 1994.
[37] R. P. Starkey. Design of waverider based re-entry vehicles. AIAA/CIRA 13th
International Space Planes and Hypersonics Systems and Technologies, 2005.
166
[38] X. Chen, Z. Hou, J. Liu, and X. Gao. Bluntness impact on performance of
waverider. Computers and Fluids, 48(1):30?43, 2011.
[39] M. W. Oppenheimer, D. B. Doman, J. J. McNamara, and A. J. Culler. Viscous
effects for a hypersonic vehicle model. AIAA Atmospheric Flight Mechanics
Conference, 2008.
[40] M. A. Bolender, M. W. Oppenheimer, and D. B. Doman. Effects of unsteady
and viscous aerodynamics on the dynamics of a flexible air-breathing hyper-
sonic vehicle. AIAA Atmospheric Flight Mechanics Conf, & Exhibit, 2007.
[41] L. Rudd, J. Hodgkinson, R. Parker, and C. Tarpley. Hypersonic stability
derivative modeling issues. AIAA Atmospheric Flight Mechanics Conference,
2010.
[42] N. Takashima and M. J. Lewis. Navier-stokes computation of a viscous opti-
mized waverider. Journal of Spacecraft and Rockets, 31(3):383?391, 1994.
[43] J. M. Burt, E. Josyula, F. Ferguson, and I. M. Blankson. Automated aerody-
namic optimization for lifting hypersonic vehicles at high altitude. 50th AIAA
Aerospace Sciences Meeting and Exhibit, 2012.
[44] C. Tarpley and M. J. Lewis. Stability and control of hypersonic waveriders.
31st AIAA Aerospace Sciences Meeting, 1993.
[45] N. Takashima and M. J. Lewis. Optimization of waverider-based hypersonic
cruise vehicles with off-design considerations. Journal of Aircraft, 36(1):235?
245, 1999.
[46] I. I. Mazhul and R. D. Rakhimov. Hypersonic power-law shaped waveriders
in off-design regimes. Journal of Aircraft, 41(4):839?845, 2004.
[47] Y. Shi, B. J. Tsai, J. B. Miles, and K. M. Isaac. Cone-derived waverider
flowfield simulation including turbulence and off-design conditions. Journal of
Spacecraft and Rockets, 33(2):185?190, 1996.
[48] Y. Shi, J. B. Miles, and K. Isaac. Computational fluid dynamics simulation of
turbulent waverider flowfield with sideslip. Journal of Spacecraft and Rockets,
34(1):76?82, 1997.
[49] S. A. Berry, R. Kimmel, and E. Reshotko. Recommendations for hypersonic
boundary layer transition flight testing. AIAA Fluid Dynamics Conference,
2011.
[50] S. A. Berry, T. J. Horvath, B. R. Hollis, R. A. Thompson, and H. H. Hamilton.
X-33 hypersonic boundary-layer transition. Journal of Spacecraft and Rockets,
38(5):646?657, 2001.
167
[51] S. A. Berry, R. A. King, M. A. Kegerise, W. A. Wood, C. B. McGinley, K. T.
Berger, and B. Anderson. Orbiter boundary layer transition prediction tool
enhancements. AIAA Aerospace Sciences Meeting, 2010.
[52] K. F. Stetson. Hypersonic boundary-layer transition. In Advances in Hyper-
sonics Volume 1. 1992.
[53] S. P. Schneider. Effects of roughness on hypersonic boundary-layer transition.
AIAA Aerospace Sciences Meeting, 2007.
[54] X. Zhong and X. Wang. Direct numerical simulation on the receptivity, insta-
bility, and transition of hypersonic boundary layers. Annual Review of Fluid
Mechanics, 44:527?561, 2012.
[55] I. Tani, R. Hama, and S. Mituisi. On the permissible roughness in the laminar
boundary layer. Aeronautical Research Institute, Tokyo Imperial University,
15(13):417?428, 1940. Report Number 99.
[56] N. Gregory and W. S. Walker. The effect on transition of isolated surface
exrescences in the boundary layer. Technical Report R&M 2779, Aeronautical
Research Council, 1956.
[57] T. N. Canning, M. E. Tauber, M. E. Wilkins, and G. T. Chapman. Orderly
three-dimensional processes in turbulent boundary layers on ablating bodies.
Technical Report AGARD CP-30, NASA, May 1968.
[58] M. S. Acarlar and C. R. Smith. A study of hairpin vortices in a laminar
boundary layer. part 1. hairpin vortices generated by a hemisphere protur-
bance. Journal of Fluid Mechanics, 175:1?41, 1987.
[59] F. G. Ergin. Roughness Induced Transition. PhD thesis, Case Western Reserve
University, 2005.
[60] M. C. Fischer. Spreading of turbulent disturbance. AIAA Journal, 10(7):957?
959, 1972.
[61] A. Fiala and R. Hillier. Aerothermodynamics of turbulent spots and wedges
at hypersonic speeds. Fifth European Symposium on Aerothermodynamics for
Space Vehicles, 563:615, 2005.
[62] L. Krishnan and N. D. Sandham. Effect of mach number on the structure of
turbulent spots. Journal of Fluid Mechanics, 566:225?234, 2006.
[63] J. A. Redford, N. D. Sandham, and G. T. Roberts. Numerical simulations of
turbulent spots in supersonic boundary layers: Effects of mach number and
wall temperature. Progress in Aerospace Sciences, 52:67?79, 2012.
168
[64] P. M. Danehy, A. P. Garcia, S. Borg, A. A. Dyakonov, S. A. Berry, J. A.
Inman, and D. W. Alderfer. Fluorescence visualization of hypersonic flow
past trianglar and rectangular boundary-layer trips. 45th AIAA Aerospace
Sciences Meeting and Exhibit, January 2007.
[65] P. M. Danehy, B. Bathel, C. Ivey, J. A. Inman, and S. B. Jones. No plif study
of hypersonic transition over a discrete hemispherical roughness element. 47th
AIAA Aerospace Sciences Meeting Including The New Horizons Forum and
Aerospace Exposition, January 2009.
[66] P. S. Iyer, S. Muppidi, and K. Mahesh. Transition of hypersonic flow past flat
plate with roughness elements. 40th Fluid Dynamics Conference and Exhibit,
July 2010.
[67] H. B. Johnson and G. V. Candler. Three-dimensional hypersonic boundary
layer stability analysis with stabl-3d. AIAA Fluid Dynamics Conference, 2010.
[68] G. Taguchi. Introduction to quality engineering: designing quality into prod-
ucts and processes. 1986.
[69] A. Gil, J. Segura, and N. M. Temme. Numerical methods for special functions.
Society for Industrial and Applied Mathematics, 2007.
[70] A. Ben-Tal and A. Nemirovski. Robust convex optimization. Mathematics of
Operations Research, 23(4):769?805, 1998.
[71] D. Wiesmann, U. Hammel, and T. Back. Robust design of multilayer optical
coatings by means of evolutionary algorithms. Evolutionary Computation,
IEEE Transactions on, 2(4):162?167, 1998.
[72] S. Tsutsui and A. Ghosh. Genetic algorithms with a robust solution search-
ing scheme. Evolutionary Computation, IEEE Transactions on, 1(3):201?208,
1997.
[73] S. Tsutsui, A. Ghosh, and Y. Fujimoto. A robust solution searching scheme in
genetic search. Parallel Problem Solving from NaturePPSN IV, pages 543?552,
1996.
[74] Y. Jin and B. Sendhoff. Trade-off between performance and robustness: an
evolutionary multiobjective approach. In Evolutionary Multi-Criterion Opti-
mization, pages 237?252. Springer, 2003.
[75] J. Branke. Creating robust solutions by means of evolutionary algorithms. In
Parallel Problem Solving from NaturePPSN V, pages 119?128. Springer, 1998.
[76] T. W. Simpson, J. Peplinski, P. N. Koch, and J. K. Allen. On the use of
statistics in design and the implications for deterministic computer experi-
ments. Design Theory and Methodology-DTM?97, pages 14?17, 1997.
169
[77] Y. Jin. A comprehensive survey of fitness approximation in evolutionary com-
putation. Soft Computing, 9(1):3?12, 2005.
[78] C. Elster and A. Neumaier. A grid algorithm for bound constrained optimiza-
tion of noisy functions. IMA Journal of Numerical Analysis, 15(4):585?608,
1995.
[79] N. D. Lagaros, V. Plevris, and M. Papadrakakis. Multi-objective design op-
timization using cascade evolutionary computations. Computer Methods in
Applied Mechanics and Engineering, 194(30):3496?3515, 2005.
[80] J. Zhang, M. M. Wiecek, and W. Chen. Local approximation of the efficient
frontier in robust design. Journal of Mechanical Design, 122:232, 2000.
[81] J. S. Kang, M. H. Suh, and T. Y. Lee. Robust economic optimization of process
design under uncertainty. Engineering Optimization, 36(1):51?75, 2004.
[82] I. Das. Robustness optimization for constrained nonlinear programming prob-
lems. Engineering Optimization, 32(5):585?618, 2000.
[83] W. Chen. Quality utility?a compromise programming approach to robust de-
sign. PhD thesis, University of Illinois, 1998.
[84] N. Rolander, J. Rambo, Y. Joshi, J. K. Allen, and F. Mistree. An approach to
robust design of turbulent convective systems. Journal of Mechanical Design,
128:844, 2006.
[85] T. Ray. Constrained robust optimal design using a multiobjective evolutionary
algorithm. In Evolutionary Computation, 2002. CEC?02. Proceedings of the
2002 Congress on, volume 1, pages 419?424. IEEE, 2002.
[86] K. So?rensen. Tabu searching for robust solutions. In Proceedings of the 4th
metaheuristics international conference, pages 707?712. Citeseer, 2001.
[87] R. P. Starkey. Investigation of Air-Breathing Hypersonic Missile Configura-
tions Within External Box Constraints. PhD thesis, University of Maryland,
2002.
[88] F. M. White. Viscous Fluid Flow. McGraw-Hill, third edition, 2006.
[89] G. Emanuel. Analytical Fluid Dynamics. CRC Press, first edition, 1994.
[90] G. Emanuel and Y. Y. Bae. Boundary-layer tables for similar compressible
flow. AIAA Journal, 27(9):1163?1164, 1989.
[91] S. C. Tirtey, O. Chazot, and L. Walpot. Characterization of hyper-
sonic roughness-induced boundary-layer transition. Experiments in Fluids,
50(2):407?418, 2011.
[92] M. Mitchell. An Introduction to Genetic Algorithms. The MIT Press, 1998.
170
[93] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan. A fast and elitist mul-
tiobjective genetic algorithm: NSGA-II. Evolutionary Computation, IEEE
Transactions on, 6(2):182?197, 2002.
[94] K. Deb. Multi-Objective Optimization using Evolutionary Algorithms. John
Wiley and Sons, 2001.
[95] H. Ishibuchi, N. Tsukamoto, and Y. Nojima. Evolutionary many-objective
optimization: A short review. In Evolutionary Computation, 2008. CEC
2008.(IEEE World Congress on Computational Intelligence). IEEE Congress
on, pages 2419?2426. IEEE, 2008.
[96] J. Wu and S. Azarm. Metrics for quality assessment of a multiobjective design
optimization solution set. Journal of Mechanical Design, 123:18?25, 2002.
[97] N. Palli, P. McCluskey, S. Azarm, and R. Sundararajan. An interactive mul-
tistage epsilon-inequality constraint method for multiple objectives decision
making. ASME Journal of Mechanical Design, 120:678?686, 1998.
[98] A. Kurpati, S. Azarm, and J. Wu. Constraint handling improvements for
multi-objective genetic algorithms. Structural and Multidisciplinary Optimiza-
tion, 23(3):204?213, 2002.
[99] R. Kalman. On the general theory of control systems. IRE Transactions on
Automatic Control, 4(3):110?110, 1959.
[100] M. G. Goman and M. N. Demenkov. Computation of controllability regions
for unstable aircraft. AIAA Journal of Guidance, Control, and Dynamics,
27(4):647?656, 2004.
171