Structural Optimization Techniques Guide
Structural Optimization Techniques Guide
Structural Optimization
3 Design Evaluation 15
3.1 Local Optimality Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Global Objective Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Constraining Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4 Several Design Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4.1 Pareto Optimality . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4.2 Substitute Problem and Preference Function or Scalarization . . . . 17
3.5 Transformation Methods and Pseudo Objectives . . . . . . . . . . . . . . . 18
3.5.1 Penalty Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.5.2 Method of Multipliers . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5.3 Unconstrained Lagrange Problem Formulation . . . . . . . . . . . . 21
3.6 Fitness Function for Evolutionary Algorithms . . . . . . . . . . . . . . . . . 22
3.6.1 Mapping Functions for Objectives . . . . . . . . . . . . . . . . . . . 23
3.6.2 Constraint Mapping Functions . . . . . . . . . . . . . . . . . . . . . 24
3.7 Design Evaluation Exemplified on Selected Problems . . . . . . . . . . . . . 27
3.7.1 Weight Minimization of a Motorcycle Tubular Frame . . . . . . . . . 27
3.7.2 Racing Car Rim Design Evaluation . . . . . . . . . . . . . . . . . . . 28
3.7.3 Maximum-Strength Flywheel Design . . . . . . . . . . . . . . . . . . 30
3.7.4 Maximum Bond-Strength Design . . . . . . . . . . . . . . . . . . . . 32
3.7.5 Composite Boat Hull . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
optimizers.
Generally, real optimization problems yield objective functions that can not be processed
with mathematical programming. Stochastic search methods do not suffer limitations like
mathematical programming and the interfacing problems are often less severe. They are
thus better suited for many practical problems but they require a much higher number of
function evaluations unless the optimum is hit early by chance. Consequently, the latest
research [1] focuses on rendering stochastic search methods more efficient by borrowing
concepts from mathematical programming. The more efficient stochastic search methods
are based on concepts inspired by the evolution of life or bacterial search mechanisms,
for instance, and involve the evaluation of whole populations of trial design solutions (or
their genotypes) within one generation of the evolution process. Since the individuals of
one population can easily be evaluated in-parallel, the use of modern massively parallel
computer architecture, such as Beowulf clusters with hundreds of processors, mitigates the
efficiency deficiency of the stochastic search methods.
2
G
G
G !
G "
tions in terms of load-carrying frames, often used for bridges, motorcycle frames, or other
structures, have a fixed connectivity of the various members with each other, or topology.
If the topology is already specified, further possibility for improving the performance of
such structures lies in the adjustment of the various member’s size. The sizing changes
the area values, or moments of inertia, of members such as trusses or beams. The illus-
tration shown on Fig. 1.1(a) indicates not only the sizing of member properties but also
the change of the position of the nodes where members connect.
The position change of connectors indicated in Fig. 1.1(a) is very similar to the shape
optimization indicated in Fig. 1.1(b). The example shows an optimized design with a
shape that resulted, through a shape-optimization process, from a rectangular-shaped ini-
tial design.
Sizing and shape optimization always require some initial design where the topology is
already fixed. Less information, or pre-existing knowledge, requires the topology opti-
mization. It requires only a definition of the physical or geometric design space and the
specification of geometric boundary conditions and loads. By redistributing the material,
which is initially evenly distributed in the design space, some topology as indicated in Fig.
1.1(c) is automatically created. The method, along with a specific solution technique [12],
is explained in section 8.1. It was used to create the title page illustration. The shown
result is the stiffest structure with respect to the sketched boundary conditions and the
fineness of the finite-element mesh. The simulated best design hints at the optimal Michell
structures [13].
Advanced composite materials consist of fibers of high stiffness and strength embedded
in some matrix the material of which may be rather weak and compliant. The resulting
composite, for instance with unidirectional reinforcement, has mechanical properties that
are highly direction-dependent, or anisotropic. Such materials are most often used for thin
shell structures, where the walls consist of laminates of several layers of the composite ma-
terial. The orientation of the reinforcement of the individual layers, see Fig. 1.1(d), can be
chosen to obtain some desired global structural behavior. Fiber orientation is one of the
internal parameters characterizing a laminate, and algorithms for the optimum adjustment
of these are called internal material parameter optimization.
Of course, different structural optimization types may be combined to solve one opti-
mization task. Topology optimization may be followed by shape optimization or shape
optimization may be coupled with internal material parameter optimization.
The motor cycle manufacturer Ducati, active in racing, uses a tubular steel trellis frames
such as the one of the Ducati 996R Superbike shown in Fig. 1.2(a). A co-operation
with Ducati lead to a student graduate project and inspired the development of some
solution technique to be published soon [3]. The company is interested in making the frame
(a) (b)
structure as lightweight as possible whilst at the same time preserving prescribed structural
stiffness properties. The structural stiffness influences the racing performance of the bike
because it contributes to the suspension characteristics. The springs and shocks provide
the suspension characteristics of the bike in the straight position. The forced displacements
or dynamic loads, exerted by the road onto the bike at certain speeds, act then parallel
to the vertical axis of the bike. In curves, surrounded at high speeds, the forces act more
sideways and the suspension system takes only a small component in the vertical direction
of the bike. The lateral component must be damped by the structural compliance of the
bike in the lateral direction and with respect to the contact points between wheels and
road. A significant factor determining the lateral suspension properties to be provided by
the frame is its torsional stiffness. Therefore, the torsional stiffness of the frame, measured
in terms of twisting moment related to relative twist between its front and rear ends, is a
fixed value that must be kept constant when reducing the weight. Such a requirement is
called a constraint, specifically an equality constraint. Another constraint is derived from
the fact that the frame must withstand the loads to be transferred by it. This implies
that the stresses induced by external loading must not exceed the material strength at
any location of the frame. In order to obtain the necessary information on the structural
stiffness and the stressing of the frame, its behavior and performance under loads must
be simulated. The simulations are provided by the finite-element method FEM and Fig.
1.2(b) shows the frame geometry that was used for the FEM modelling.
Figure 1.3: CAD-Model of the racing car rim (courtesy O. König [5])
not only influence the performance as a part of the cars overall weight. Since the wheels
belong to the so called unsprung mass, a low rim weight improves the mechanical grip of
the car especially on bumpy road surfaces. The moment of inertia along the wheels
rotational axis should be minimal for several reasons. Low moments of inertia allow faster
acceleration and deceleration of the wheels and therefore of the whole car. Furthermore,
the moment of inertia leads through the gyro effect to higher steering forces as well as a
higher inertia of the car with respect to direction changes. Finally, the stiffness of the
rim is of high importance in turns at high speed. Vertical loads of 5700N , resulting from
the car’s weight and aerodynamical descending forces, as well as maximum lateral forces
of 7000N , as a result of the centripetal forces, build up a bending moment on the rim.
Additionally, strength requirements must be fulfilled. Plastic yield must not occur in use,
whereas some parts reach temperatures well above 200◦ C. A special magnesium alloy
is used where maximum yield stress does hardly decrease with higher temperatures and
the mass-specific stiffness ratio is high. Manufacturing starts with a forging blank, the
rim’s bed is shaped by CNC-lathe, and the spokes form the interspace of CNC-milled
pockets. The forging blank’s shape is not to be changed, as this would exceed costs. FIA1
regulations affect the bead diameter as well as dimensions of the lower rim-bed.
For the optimization presented in here, maximum bending stiffness is defined as main
design objective from the rim manufacturer. Nevertheless, the desired properties
• Low mass
of the existing design should also be matched or even surpassed by the optimized design.
Furthermore, the optimization must also take into account the functional, regulatory, and
manufacturing requirements discussed earlier in this section. Altogether, this constitutes
1
Fédération International d’Automobiles ([Link]
a highly constrained optimization problem, which can not be tackled with classical math-
ematical optimization techniques. Since all original data is confidential, arbitrary new
geometries for the rim and the forging blank, as well as modified load cases are created
and optimized for the presentation in here.
(a) (b)
Figure 1.4: Flywheel with central bore, 2-dimensional FEM model and stresses [2]
r s a n d w ic h
r o n s e rt
r b o re
O n s e rt h o n s e rt
B o n d L a y e r h g lu e
F a c e S h e e t
C o re h c o re
F a c e S h e e t h fa c e
(a) (b)
(c) (d)
Figure 1.5: Onsert design demonstrator (a) and geometry model (b)
The hull of a sail boat, made of composite materials, should be as stiff as possible under
typical service loads. The particular model shown in Fig. 1.6 should be low priced for
marketing reasons. The sample problem is thus useful to demonstrate a problem where a
Figure 1.6: ANSYS model of the sail boat hull with composite material patches (Courtesy
N. Zehnder [14])
2
[Link]
3
[Link]
Treatment of a Structural
Optimization Problem
Optimizing a structure by some automated numerical procedure may seem very complex
and difficult to organize. The concept presented in the following sections was worked out
by Eschenauer [18] and decomposes the task into manageable subtasks so that it can be
solved in a straightforward manner. Although Eschenauer’s concept seems to have been
developed with regard to the optimization algorithms labelled by the term mathematical
programming, it is also valid when other solution techniques such as genetic algorithms
are used.
d e c is io n d a ta
m a k e r
in p u t
o p tim a l
d e s ig n
x * in itia l d e s ig n s tru c tu ra l
x 0 p a ra m e te rs
tra n s fo rm e d a n a ly s is
tr a n s fo r m a tio n d e s ig n m o d e l
v a r ia b le z ® x d e s ig n x ® y v a r ia b le
v a r ia b le
z x y
o p tim iz a tio n s tru c tu ra l m o d e l
a lg o r ith m o p tim iz a tio n m o d e l s tr u c tu r a l a n a ly s is
u = u (y )
F ,g ,h B ,f ,g ,h
o p tim iz a tio n e v a lu a tio n
s tr a te g ie s m o d e l u
¶ F ¶ g ¶ h
, ,
¶ x ¶ x ¶ x
s e n s itiv ity
a n a ly s is
Design Evaluation
Design evaluation and the setting up of an evaluation model is discussed in chapter 3
on the basis of the sample problems presented in chapter 1. In structural optimization,
one uses generally FEM to obtain the response of a structure to loads under specified
geometric boundary conditions. The solution of the numerical system of equations obtains
the primary solution in terms of the nodal-point degrees of freedom. In structural analysis,
the degrees of freedom are displacements. From the primary solution other results such
as stresses can be obtained. The stresses may be used to formulate an objective when
the strength of a part is to be maximized. Stresses may also be used to formulate a
constraint when, for instance, the weight of a load-carrying part is to be minimized but a
required strength must be preserved. The weight is then calculated from the integral of
the material densities over the volume of the considered part and does not depend on the
load response. Generally, however, the structural response is needed for the evaluation of
both, the objective and the constraints.
Design Parameterization
Transformed Variables
The design variables may be transformed to meet certain requirements of the optimization
algorithm. For instance, genetic algorithms operate upon genotype variables that may be
defined in terms of bit strings whose encoded information must be transformed into the
phenotype design variables. The topic of transformed variables is therefore discussed in
context with genetic algorithms in section 6.3.
Input files for the simulation programs. The actual simulations to be carried out are
established with the input files, as well as the effective results to be calculated and
stored in output files.
Mapping from genotype to input files. Every optimization parameter in the input
files must be linked to the appropriate gene of the eoUniGene genotype. This is
done by storing the exact row and column, where every parameter must be inserted
in the input files.
Simulation
Objective & Constraint Program
Value Reading Manager File File File
File File
Simulation Simulation
Program B File Program A
Figure 2.2: DynOPS evaluator to calculate fitness values for a population of individuals
using external simulation software.
Objectives, constraints, and fitness function. How to read the objectives and con-
straints from the appropriate output files, as well as how the actual fitness should
be calculated from these values must be defined.
The simulation program manager starts the first simulation program together with the
appropriate input files, and waits for job completion. The results from the evaluation
stage stored in output files are either used as input for the next simulation program (e.g.
a geometry file), or are directly used for fitness calculation (e.g. a mass evaluated in the
CAD system). The next evaluation stage in the sequence is then started, and so forth.
After completion of the sequence of simulation programs, objective and constraint values
stored in result files are transferred back to the DynOPS evaluator. The evaluation of an
individual is completed by computing its fitness value. This evaluation loop is repeated
for the whole population.
Design Evaluation
and if the objective is actually to maximize some property such as volume V , the objective
function can be set up so that minimizing it is equivalent to maximizing that property.
This is easily achieved by multiplying the property with a negative number, f (x) = −V (x)
1
or dividing a positive number by it, f (x) = V (x) . In context with evolutionary-algorithm
solution techniques, the objective function is also called a fitness function.
At some stage during an multi-objective optimization process the situation appears that
a further minimization of one objective function goes on account of increasing some other
objective function value. Such a situation is called an objective conflict because none of
the possible solutions allows for simultaneous optimum fulfillment of all objectives.
Fig. 3.1 shows a projection from the two-dimensional design space X into the objective
function space Y . The Pareto optimal solutions lie on the lines AB. The designer may
choose from these solutions from assessment of the relative values of the two objective
functions.
Figure 3.1: Mapping of a feasible design space into the criteria space [18]
such that
p [f (x̃)] = min p [f (x)] . (3.9)
X∈R\
Eschenauer [18] cites various formulations of the preference functions from which here only
the sum of weighted objectives is mentioned:
m
X
p [f (x)] := [wj fj (x)] , x ∈ Rn . (3.10)
j=1
m
X
0 ≤ wj ≤ 1, wj = 1. (3.11)
j=1
If all objectives are convex, a full set of Pareto-optimal solutions can be generated by
running a sequence of substitute scalar problems where the preference functions cover an
appropriate range of values for the weighting factors wj . If one or several of the objectives
are not convex, the Pareto-optimal set is not so easily generated and Eckhart Zitzler
explains the topic in his lecture class Bio-Inspired Computation and Optimization at ETH
Zurich.
The function Ω can be defined so that either the exterior point method or the interior
point method results [20]. A disadvantage of both penalty methods shown in the follow-
ing subsections is that the respective pseudo-objective functions topologies make it more
difficult for deterministic optimization methods to perform well.
Therefore, of the inequality constraining functions g(x), only the active ones may be consid-
ered in the penalty formulation (3.13). Inserting it into (3.12) results in an unconstrained
function the minimum point of which lies outside the feasible region, hence the name of
this method. As Fig. 3.2(a) shows, with increasing values of the penalty parameter R the
minimum moves closer to the feasible region but it can never quite reach it.
Note that the interior-point penalty method does not work with equality constraining
functions h as violations of these may lead to either positive or negative values. Fig.
3.2(b) illustrates that with decreasing value of the penalty parameter R0 the minimum
point of the transformed objective function moves closer to the infeasible region or the
constrained minimum point.
4 4
3 3
2 2
1 1
0 0
-2 -1 0 1 2 3 -2 -1 0 1 2 3
Figure 3.2: Examples for the exterior (a) and interior (b) penalty functions
Thus, the feasible region of the optimization variable is 1 < x < 2. The two side constraints
define the inequality constraining functions
It can be seen that the minimum point of the constrained original problem is at x = 1.
The substitute problem resulting from the transformation with the outer penalty method
has its minimum point between that of the unconstrained and the constrained original
objective functions, −2 < x∗ < 1. In Fig. 3.2(a) the unconstrained function f and three
transformed functions p are plotted for the penalty parameter values 1, 10, and 100. It
can be seen how the minimum of the transformed function moves closer to x = 1 as the
penalty parameter increases from 1 to 100.
The inner penalty method obtains an unconstrained substitute function whose minimum
point lies within the feasible region. As the penalty parameter R0 decreases from 100 to 1
the minimum point moves closer to the edge of the infeasible region. It can also be seen
from Fig. 3.2(b) that it will take very small values of R0 to move the minimum of p close
to x∗ = 1.
One might conclude that solutions close to the constrained minimum point can be obtained
simply by using very high or small values for the penalty parameters R or R0 , respectively,
and just minimizing the pseudo objective for these values. However, it can be seen from
Fig. 3.2 that the unconstrained pseudo objectives are more distorted if compared to the
original objectives. The search methods of mathematical programming (section 6.3) are
designed to work best on objective functions that behave almost like quadratic functions
and might fail at the highly distorted pseudo-objectives resulting from choosing the penalty
parameters so that the minimum point of the unconstrained pseudo objective is very close
to the true constrained minimum point.
It is then inevitable that the constrained optimization problem is solved by a sequence
of unconstrained subproblems, where the penalty parameters are updated at each step.
Considering the exterior point method, the parameter R is chosen small, for instance zero,
at the first stage, and gradually increased with the subsequent stages. For the interior
point method, one starts with a high value of R0 and decreases it from stage to stage. The
minimum point of each subproblem is then used as a starting point for solving the next
subproblem. So the considered regions in search space become smaller as the pseudo objec-
tives become more distorted. Nevertheless it is inevitable that the generated subproblems
become progressively ill-conditioned so that, at one point, the sequence terminates not
because of finding a very close approximation to the true constrained minimum point but
because of failure of the search algorithms.
where R is a constant scale factor (R may vary from constraint to constraint but remains
constant from stage to stage), and the bracket operator is defined as
α if α > 0
hαi = . (3.19)
0 if α ≤ 0
The σj and τk parameters are constant during each unconstrained minimization but are
updated from stage to stage. It is not necessary for the starting vector x0 to be feasible,
and the parameters can be conveniently chosen for the first stage as σ = τ = 0. Thus the
first minimization stage is identical to the first unconstrained minimization using standard
exterior point method penalty terms.
Multiplier estimates for the (t + 1)st stage are formed according to the following rules:
D E
(t+1) (t)
σj = gj (x(t) ) + σj j = 1, 2, 3, ...J
(3.21)
(t+1) (t) (t)
τj = hk (x ) + τk k = 1, 2, 3, ...K
Because of the bracket operator, σ has no negative elements, whereas the elements of τ
can take either sign.
Suppose that improvement is sought by moving from a reference point, which must be
feasible, along a search direction s. The search direction is the linear combination of
N
Ñ D
I ×Ñ D =
I
Ñ B
D = I = - Ñ B - l Ñ D
N
a usable direction, along which smaller values of f are found and an obvious choice of
which is the direction of steepest descent, and the gradients of the equality constraining
functions:
s = −∇f − λT ∇h (3.23)
The search direction is feasible, or will not violate the linear constraint, if it is orthogonal
to the constraining function gradient:
sT ∇h = 0 (3.24)
Inserting the definition (3.23) of the search direction into the orthogonality condition (3.24)
yields the Lagrange factors λi :
∇f T ∇hi
λi = − (3.25)
∇hTi ∇hi
The Lagrange function, or Lagrangian, is the function whose negative gradient provides a
feasible search direction obeying all constraints.
m
X l
X
L(x, λ) = f (x) + λj gj (x) + λm+k hk (x) (3.26)
j=1 k=1
More information on the method of feasible directions, the Lagrangian, and the Kuhn-
Tucker conditions for constrained optimization problems is given in Sections 5.4 through
5.6. Section 5.5 on page 55 explains why the stationary point of the Lagrangian is a
saddle point and Section 5.7 on page 59 introduces the concept of duality, where the dual
problem is that of solving the constrained optimization problem in terms of the Lagrange
multipliers. Section 6.9 considers numerical solution of the Lagrangian by searching for a
minimum. If the constraining function g is not linear, the found search direction is feasible
only at the reference point and moving along it will eventually violate the inequality
constraints. It will then be necessary to remove the violations, or find a feasible point
close to the infeasible one just obtained. Such an algorithm is explained in Section 6.9.2.
The demands represent ratings for one or several objectives and constraints so that the
fitness function F appears on first sight analogous to the pseudo-objective function P
explained in section 3.5.1.
Recent research [6, 5] has elaborated schemes for defining the ratings in such a way that
the problem of finding appropriate weight values dissolves. The following is direct citation,
or uses material, from Oliver König’s Ph.D. thesis [5].
In order to avoid that one of these terms becomes much larger than the other
ones and therefore dominant, only bounded functions scaled to the interval
[0, 1] are used. Moreover, this facilitates the adjustment of the weight coeffi-
cients wi . Further, to enhance general usability of the fitness formulations, the
mapping functions Di (~ p) are defined range-independent. This means that a
certain mapping function does not change its behavior if applied to objectives
or constraints operating in different number ranges.
Based on these requirements, functions Di for the different possible types of
optimization objectives and constraints are presented. The focus for the formu-
lation of these mapping functions is put on good practical usability. The user
of an evolutionary design optimization program shall be able to define good
mapping functions by only bringing in know-how about the problem he wants
to solve. Thus three types of general mapping functions are defined with their
defining parameters as listed in Table 3.1. The defining parameters are chosen
so that they relate directly to engineering practice. For a problem at hand,
the parameters Oinit and Oestim refer to the initial value and the estimated
best-possible value of the design objective respectively. The design objective
function can optionally be modified using an amplification factor αpost . A
limit-value constraint is defined through Climit for the given limit value, and
through a parameter Cfeas tol specifying a tolerance range for constraint values
still considered feasible for the problem at hand. Finally there are target-value
constraints defined through three parameters. Ctarget defines the target value
to be achieved. It is useful to specify an admissible tolerance Cadm tol , defining
an interval for the target value to be reached at the end of the optimization.
Additionally a feasible value tolerance Cfeas tol should be specified, defining
which constraint values should still be considered during optimization.
In the following, fitness functions for the different types and parameters are
presented.
1. The resulting fitness values must fit into the interval [0, 1].
2. Relevant design improvements should be reflected in distinct decreases of
Di (O0 ).
3. Selection pressure is initially strong and slows down close at the optimum.
Di (O = Oinit ) = 1
(3.30)
Di (O = Oestim ) = 0.1
Oinit represents an initial value of the design objective, which shall result in
the maximum fitness value 1. Oestim is the estimated goal value that can be
achieved in the optimization. The fitness value 0.1 for Oestim was adjusted to-
gether with an exponential factor α = 5 in order to fulfill the third requirement
defined above. Furthermore, a small fitness value for the estimated objective
also ensures conformance with the second requirement defined. The scaling
factors a and b can be computed as:
√
1 − α 0.1
a = (3.31)
Oinit − Oestim
b = 1 − aOinit
The user has only to specify Oinit and Oestim to define the fitness function for
a design objective of a problem at hand. The bold line in Fig. 3.4 pictures an
O b je c tiv e V a lu e
5 0 6 0 7 0 8 0 9 0 1 0 0 1 1 0
1 .2
1 .0
0 .8
0 .6
F itn e s s
0 .4 3
4
0 .2 5
6
0 .0 8
1 0
-0 .2
Figure 3.4: Fitness function for a design objective defined through Oinit and Oestim .
example of the fitness function computed for Oinit = 100 and Oestim = 60. The
graph demonstrates that the three given requirements for the fitness function
are met for arbitrary ranges of objective values. As a last tuning parameter
αpost is introduced. Leaving the scaling parameters a and b constant, the
exponential factor α can be varied subsequently to adapt for the problem at
hand as presented in Fig. 3.4. Usually these variations will be of negligible
influence on the overall performance of the algorithm.
With this penalty function a strong selection pressure is exerted towards so-
lutions that meet the constraint. On the other hand, the constraint function
makes no distinction between values that hurt the constraint only marginally
and values that clearly violate the restriction. However the former values can
contain valuable information for the problem and should therefore not be ex-
cluded strictly. In order to take this fact into account, a smoothed step function
is defined as
1
Di (C) = −λ(C(~p)−Climit −∆)
(3.33)
1+e
allowing the algorithm to also extract valuable information from solutions with
marginally violated constraints. The step function is controlled with two pa-
rameters: ∆ adjusts the horizontal positioning of the function, and λ deter-
mines the steepness of the step. However, the parameters ∆ and λ are difficult
to adjust correctly, since they depend on the range of occurring constraint val-
ues. Furthermore it is difficult to estimate their quantitative effect on the step
function.
Based on these findings, another definition of the step function has been de-
veloped. For practical purpose it would be much more comfortable to define
the step function only with its limit value Climit , and an additional tolerance
Cfeas tol giving an upper limit of feasible constraint values to be taken into ac-
count by the EA. A method to define the step functions that way is to specify
two conditions
where Dlimit is the penalty value typically reached at the end of an optimiza-
tion, where the EA has found an equilibrium between the different demands
of the fitness function. Dfeas corresponds to the penalty value which is typi-
cally still taken into account by the EA. For the normalized fitness formulation
used within this thesis, Dlimit = 0.01 and Dfeas = 0.5 proved to be reasonable.
With these two conditions given, the original parameters of Equation 3.33 can
be computed as
1 1 1
λ = ln − 1 − ln −1 (3.35)
Cfeas tol Dlimit Dfeas
1 1
∆ = ln −1
λ Dlimit
Fig. 3.5 presents the mapping functions for a critical value given as Climit = 60,
and for tolerance values Cfeas tol = 0.6...6. This corresponds to a 1 − 10%
tolerance of the limit value Climit . For this formulation, it has to be Cfeas tol ∈
R+ for upper limits and Cfeas tol ∈ R− for lower limits, respectively.
1 − e− 2σ 2 : |C(~
p) − Ctarget | ≥ Cadm tol
(3.36)
where C(~ p) refers to the actual constraint value. The parameter σ determines
how fast the penalty value increases when the acceptance interval is left. In this
C o n s tr a in t V a lu e
4 0 4 5 5 0 5 5 6 0 6 5 7 0 7 5 8 0
1 .2
1 .0
0 .8
0 .6
F itn e s s
0 .4 0 .6
1 .0
0 .2 2 .0
c to l 3 .0
0 .0 4 .0
6 .0
-0 .2
Figure 3.5: Upper limit constraint penalty functions defined through Climit and Cfeas tol .
form, the penalty function can be used to solve practical problems. However,
the parameter σ depends on the number range of the constraint considered,
and is therefore very difficult to adjust for a problem at hand. To paraphrase
a constraint function defined through the parameters given in Table 3.1, an
additional condition is introduced.
As introduced before, the feasible value tolerance Cfeas tol determines which
constraint values should still be taken into consideration during optimization.
Dfeas corresponds to the penalty value which is typically still taken into account
by the EA. For practical purpose, Dfeas = 0.5 proved to be reasonable. The
initial parameter σ can now be determined as
c fe a s _ to l
c a d m in _ to l = 5
F B
M T
F B
Figure 3.7: FEM model of the frame and load cases, courtesy =. König[5]
and the tubes are suffering mainly loads in the axial directions. However, since they
are welded together at the joints, bending will also contribute slightly to the structural
stiffness. Thus, the fitness of the design is more sensible to changes of the cross-sectional
areas of the tubes but depends also slightly on their second moments of inertia. The area
A, depending on the geometrical parameters diameter D and wall thickness t, influences
both the weight and the extensional stiffness linearly. Given a certain area value, the
bending stiffness increases quadratically with increasing diameter:
1
I = AD2 . (3.40)
8
Due to the influence of bending stiffness on the fitness tubes will tend to grow to large
diameters and small wall thicknesses. It is more difficult to weld tubes with small than
those with large wall thickness values. Tubes with smaller diameter and larger wall thick-
ness are therefore preferable. A further term is introduced to the fitness function, giving
tubes with thin walls a slight penalty,
PNtubes 1
i=1 (αi (ti (g)−tmin )+1)βi
i
Dthick (g) = , (3.41)
Ntubes
B2
B1
• strength constraints
• geometric constraints
The mass and the moment of inertia around the wheel’s axis of the existing design are low
and the same values should be attained by any new design. This is achieved by using the
target-value constraint mapping function 3.36. The complex geometry makes it necessary
to use CAD as well as FEM to simulate and evaluate the design. Therefore, the mass and
moment of inertia can be evaluated by the CAD system but the bending-stiffness objective
and the strength constraints require costly FEM analysis.
For evaluating the strength constraints, the load case combined from the loads indicated in
Fig. 3.9 is considered. The loads are applied to the rim’s shoulders – since this is the only
area of contact with the tire – and include horizontal (car’s weight plus aerodynamic force),
vertical (centripetal), and rotational forces (braking). Essential boundary conditions of
the analysis model simulate clamping at the contact area with the car’s suspension. The
margin of safety for mechanical stresses is defined as
σ
msaf ety = − 1, (3.42)
σyield
where σ is the acting stress and σyield is the yield stress of the material at a given tem-
perature. Values less than zero indicate plastic yield to occur. The margin of safety is
calculated for every node, a counter S(p) records the number of nodes with msaf ety < 0.
• manufacturing constraints
• assembly constraints
and are kept by the parameterization of the design model, which will be explained in
Section 4.2.2.
The manufacturing constraints include the given shape of the forging blank. A 1mm
distance to its outer contour must be kept to allow for properly machined surfaces. In
addition, manufacturing techniques must not be changed. Finishing is restricted to a
CNC-lathe and a CNC-mill. Also, the rim’s bed wall thickness must not be thinner than
2mm.
The assembly constraints include contact areas to the car’s suspension as well as to the
nut holding the wheel and the tire. An additional constraint is a 3mm distance to the
brake assembly positioned on the inside of the rim.
The FIA regulations applying to the rim include that the maximum bed diameter is limited
to 330mm. The minimum depth of the lower bed is 13.57mm and the maximum distance
from the outside surface is 43.3mm.
f = trω 2 πρ (3.44)
t
(σr t) ,r + (σr − σθ ) + trω 2 ρ = 0. (3.45)
r
The equilibrium equation 3.45 shows that the radial and the circumferential stresses, σr
and σθ , increase quadratically with increasing rotational speed ω. Thus, the rotational
speed, and with it the kinetic energy, of the flywheel are limited by the maximum stress
that the material can bear.
As outlined in section 1.3.3 and shown in Fig. 1.4, the distribution of stresses in a flywheel
of constant thickness and with a central bore is uneven. Particularly, the circumferential
stress increases steeply towards the edge of the bore. Failure will initiate at the bore
although other regions are not critically stressed.
1
Fédération International d’Automobiles ([Link]
Let the symbol σ stand for a failure criterion such as maximum principal stress σI or von
Mises equivalent stress σeqv . The problem with evaluating the objective (3.46) is that not
only the stress values at all points along r must be calculated but that the location at
which the maximum stress appears may change when the design changes. This causes the
first and second derivatives of the objective function to be discontinuous. Consequently,
one would be forced to use genetic algorithms instead of a mathematical programming
technique which would increase the numerical solution effort considerably.
A continuous objective function with continuous derivatives can be constructed by the
following argument. The maximum stress at some point in the flywheel is minimized when
its value equals the minimum stress at some other point. Then, the stress distribution
must be a constant. This idea was expanded by Stodola [24] who derived analytically
the shape of a turbine disk of constant strength. The constant stress distribution implies
that the local stress σ(r) is everywhere equal to the mean stress σ̄ and that therefore its
variance (average quadratic deviation) is zero:
Z
σ(r) = σ̄ → s = r (σ(r) − σ̄)2 dr = 0. (3.47)
r
In reality a perfectly spatially constant stress state can not be reached for physical reasons:
the surfaces perpendicular to the radial direction at the bore and at the outside are stress-
free and the radial direct stress must drop off to zero. This causes the stress distribution
inevitably not be a constant. However, a global objective function f can be based on
(3.47) so that the objective is reached by minimizing
Z
f = r (σ(r) − σ̄)2 dr. (3.48)
r
Now let us consider again the two assumptions of the mass and the rotational moment of
inertia remaining constant during the optimization process. They constitute two equality
constraints which can be written in terms of constraining functions h1 and h2 ,
Z
h1 = 2πρ r (tr − t0 ) dr = 0 (3.49)
r
and Z
2
h2 = ω πρ r3 (tr − t0 ) dr = 0. (3.50)
r
H E
Figure 3.10: Region out of the Flywheel Analysis Model with Nodal and Optimum Stress
Points of Quadratic Serendipity Type Finite Elements
from the mean stress σ eqv . For p = 1, the objective function becomes the variance (square
of the standard deviation), of the stress distribution. Higher values of p further penalize
stress peaks so that the absolute values of minimum and maximum stress deviate less from
the mean stress. A value of p = 4 has been used for the sample calculations. The stresses
are evaluated at the optimum stress points of each finite element in the bonding layer.
The parameterization of the onsert problem will be described in section 4.2.4. However,
the shape optimization, to maximize the bond strength, regards the thickness distribution
of the onsert. There it is assumed that the interface between the onsert and the bonding
layer remains straight, preserving the initial constant thickness distribution of the bonding
layer. Thus, the thickness t(x) = p2 (x) − p1 depends only on the vertical position of the
points p2 (x) on the upper onsert surface. The thickness must be positive and greater than
or at least equal to some predefined small minimum-thickness value t which gives the
constraining function
g(x) = t − t(x) = t + p1 − p2 (x) ≤ 0 (3.52)
Figure 3.11: ANSYS model of the sail boat hull with composite material patches
• sufficient strength
The parameterization of the hull is explained in Chapter 4 and allows the laminate con-
struction to change over the hull’s shell area. It is foreseen that the hull is made from
glass-fiber as well as carbon-fiber reinforced prepreg material. Carbon fibers are much
more expensive than glass fibers and so the material cost responds to the total amount
as well as the amount ratio of the two types of material. However, under the given mass
constraint the carbon fiber may be more effective, if used in certain regions, because of its
low mass and high stiffness properties. The upper limits on mass and cost are expressed
as demands D by using the mapping function 3.32. Their evaluations do not require FEM
analysis.
The finite element model is needed for the stiffness and strength evaluations. It is assem-
bled from shell elements. Their extensional and bending stiffness depends on the local
laminate construction which is given by the number of layers, their respective thicknesses
and the materials and orientation of principal material axes used for each layer, and is
evaluated by the theory of laminates plates [25]. Once the primary unknowns of the finite-
element model are known, the stresses and failure probabilities of the laminate can be
evaluated.
sym
sym
Fb
ps
Figure 3.12: Quarter end-plate model in the initial design (courtesy [5])
The other objective of having an even pressure distribution on the bipolar plate cross-
section areas is reached independently of, or after completing, the optimization procedure.
In the middle of the stack, the cross section will be plane because of symmetry. When
the other cross sections also remain plane, the interface between the fuel-cell stack and
the end plates should be so as well. Then, the stack is under a state of plane strain
but the pressure will not necessarily by constant everywhere because of Poisson’s ratio
effects. However, the animus behind the constant pressure idea is that the bipolar plates
be gas tight and one may expect that that goal is equally well reached by demanding the
longitudinal strain to be constant. A structural analysis modelling the optimized design
will reveal the out-of-plane displacement field in the interface between the end plate and
the stack. The displacement field is interpreted as a shape and the end plate’s bottom
surface will then be given the negative of that shape. The displacements under load will
then cancel out that shape so that the bottom becomes plane, consistent with a plane
state of the stack. The surface shape variations are expected to be small enough so that
there effect on the bending stiffness is negligible, decoupling the objective of gas tightness
from the minimum weight objective and stress constraints.
a ) C o n s tr u c tio n
b ) T o p o lo g y
M a te r ia l
c ) S te e l A lu m in u m C o m p o s ite
P r o p e r tie s
G e o m e try
d )
o r S h a p e
e ) S u p p o rt
o r lo a d in g
f) C r o s s s e c tio n
o r s iz in g
Figure 4.1: Classification of design optimization problems for truss-like structures in terms
of different types of design variables, after Eschenauer [18]
a) Constructive layout
The determination of the best suited layout requires optimizing each layout coming
into consideration and comparing the calculated optimum solutions.
b) Topology
The topology or arrangement of the elements in a structure is often described by
parameters that can be modified in discrete steps only. Different topologies can
also be obtained by eliminating nodes and linking elements. Note also the topology
optimization method introduced by Bendsøe and Kikuchi [12] explained in section
9.1.
c) Material properties
The material properties of isotropic building materials such as steel or aluminum
describe the stiffness in terms of Young’s modulus or Poisson’s ratio, the strength
in terms of the yield stress or other strength limit, or the weight in terms of specific
weight or mass density. The designer can often select the most suitable material
from a selection of alloys and sometimes he may decide whether steel, aluminum
or any other type of metal alloy would best meet the requirements of the design
objective. All of these choices are discrete in nature, leading to a discrete design
variable set in terms of the materials contained in a data base. Only laminates made
from anisotropic composite materials have some continuous design variables in terms
of fiber orientation.
f) Sizing
Structures in terms of members such as bars, trusses, beams, plates, or shells and also
their FEM models offer properties such as thickness, cross-sectional area, moment
of inertia, or thickness to be used as design variables. It is important to distinguish
between independent and dependent design variable. When a cross-section geometry
is determined by some variable, the geometrical properties listed above are depen-
dent on that variable. Sizing usually leads to a discrete optimization problem as
commercially available members with I-sections or channel sections come in different
discrete sizes.
The foregoing classification is exemplified on truss-like structures but remains valid for
shells or solid body structures. A number of design parameterization models and trans-
formations between design and analysis variables are discussed on the basis of the sample
structural optimization problems used for the lecture class.
ra= 1 5 m m
ra= 2 0 m m
ra= 2 5 m m
one can choose from a catalogue or a set of existing discrete design variable values.
Considering custom made tubes, one can freely adjust the inner and outer diameters as
continuous variables but the number of the customized tube sizes should be small, say
three. All tubes of the frame should then be chosen from the small number of custom
made tubes and when it is not predefined which of the tubes are to be of the same size,
a choice must be made. That choice again renders the optimization problem discrete.
The structural model must be so that all finite elements along one of the tubes must
1 1
2 3
2 15 9
1 1 2
2 1 2 14
2 1 2 1
2 2 8
1 2 12 13
2 4 7
3 1
1 2 2
1 3 6
2 11
1 10
2 1 1 3
2 5
1
2 2
Tube types: 1 : ri= 12.5mm, t=1.5mm
1
2 : ri= 9.5mm, t=1.5mm
3 : ri= 6.75mm, t=1.5mm
Figure 4.3: FEM-Model and Parameterization of the Motorcycle Frame, after [5, 6]
have the same properties determined by the respective design variable. Consequently, the
independent design variables must be transformed into the dependent geometric element
cross-section property value variables, and these must be assigned to all the elements
constituting the respective tube. Also, the frame must be symmetric with respect to the
midplane spanned by vertical and longitudinal directions so that each tube on one side of
the midplane must have the same dimensions as the respective tube on the other side of
it.
Half of the rim bed. Spoke body with the two pocket features subtracted
All these features are built on fully parametric two-dimensional sketches. In addition to
these basic features, various chamfers are applied, chosen completely parametric as well.
The basic design is taken from an existing racing car rim.
Parameterization
For performance reasons as outlined before, a parameterization including implicitly as
many of the mentioned constraints as possible has to be found, without excluding relevant
feasible solutions. The bed is parameterized by 9 wall thicknesses including the inner and
the outer bead as shown in Figure 4.5. With this approach the manufacturing constraint
of a minimum wall thickness of 2 mm can directly be included. The bed’s outer contour
can not be altered in some sections. The outer bead diameter and some shoulder diameters
are restricted by FIA regulations. The inner bed’s maximum diameter is limited to the
shoulder’s diameter to allow the assembly of the tire. The lower bed’s minimum depth
and the maximum distance from the rim’s outer face are also subject to FIA regulations.
Therefore the coordinates of the point in question also form parameters allowing to comply
with those restrictions. For the remaining degrees of freedom of the bed, contour lines of
the forging blank and the brake assembly limit the geometric design space. To make sure
wall thicknesses
distances
point
coordinates
radii
the minimum distances of 1 mm to the forging blank and 3 mm to the braking contour
can be approached as closely as possible without crossing them, the contour lines form
construction elements in the sketch and the bed’s contour is directly dimensioned to those
contours, setting the distances as parameters. The rotational body for the spokes is
parameterized in a way similar to the rim’s bed. Front and rear contour are dimensioned
to the blank’s contours. The parts of the contour-forming interfaces to the suspension
and the nut are non-parametric. The interface between spokes and bed depends on the
shape of both features. The sketch for the bed contains the spokes’ rear and front contour
lines as construction elements. The second sketch for the rotational body of the spokes is
referenced to those construction elements and the interface line with reference dimensions.
This way, the parameterization for both features is done in only one sketch and update
loops are avoided. At the same time, two separate features are needed, because the base
part between the spokes is not defined by prismatic pockets but the contour of the bed.
The pockets are removed from the spoke body going beyond the outer diameter. The
remaining spokes are then added to the rim bed. The sharp-angled base of the spokes is
then smoothed out by two parametric fillet features. The two pockets are parameterized in
a third sketch (see Fig. 4.6). A base line with a constant radius is defined. Parameters for
Poc
ket2
Pocket 1
the larger pocket include height and width at the base as well as at a continuous transition
point, and the radius of the upper rounding. The gap between the two pockets, forming
the spoke, is parameterized by two widths, one at the base and one at the same transition
point. This leaves only the upper rounding as a free parameter of the small pocket.
t0+ , t t0
Figure 4.7: Relation between nodal point coordinates and shape parameter t
[2]
very little pre-existing knowledge, allowing a great design freedom: assuming an analysis
model with 100 element columns in the radial direction, there are 101 element interface
lines associated with independent thickness values to define a design with individual fea-
tures. Thus, it became possible that the numerical shape optimization results suggested
a simplified mechanical model for finding the optimum flywheel design by exact formulae
[2].
p 1, p 2 i
F o p tim u m - s h a p e p o s itio n s (p 2 i - p 1) m in
the values of p2 changes the thickness distribution of the onsert. This allows to study the
influence of the onsert structural stiffness properties on the bonding strength. In the initial
configuration the onsert has a rectangular cross section with a regular mesh consisting of
rectangular elements evenly arranged in rows and columns. Mesh distortion is minimum
when the positions of nodes between the unmovable bottom node and the top node are
adjusted so that a constant spacing in the axial direction is maintained as can be seen
from the mesh plot shown in Fig.
The number of entries of p increases with increasing fineness of the finite-element mesh in
the radial direction and may be quite high when an accurate stress analysis is required.
The number of optimization variables becomes decoupled from mesh size when the nodal
point positions p are made to depend on a different set of design parameters x. By
choosing that the onsert surface shape is composed of simple polynomial functions P, the
optimization variables x control the original shape parameters p as weight coefficients of
the polynomials P:
p = P n xn , 0 ≤ n ≤ N. (4.1)
Einstein’s summation convention is implied in (4.1) and the set of polynomials P is com-
plete, ranging in degree from zero to a specified maximum value n.
The shape depicted in Fig. 4.8 corresponds with a polynomial including the constant,
linear, and quadratic terms.
t(ξ) = t0 + t1 ξ + t2 ξ 2 + t3 ξ 3 (4.5)
A first step is to require that the minimum thickness is kept at the points ξ = −1, ξ = 0,
and ξ = 1:
g1 = g(−1) = − t0 + t1 − t2 + t3 ≤ 0
g2 = g( 0) = − t0 ≤ 0 (4.6)
g3 = g( 1) = − t0 − t1 − t2 − t3 ≤ 0
If these conditions are satisfied, violations in the interior can only exist together with real
polynomial extrema. These are found with the condition:
1
q
2 2
g(ξ),xi = 0 → t1 + 2t2 ξ + 3t3 ξ = 0 → ξ1,2 = − t2 ∓ t2 − 3t1 t3 (4.7)
3t3
Real extrema do not exist if the radikant of the root is negative, or if the discriminat D is
positive:
D = 3t1 t3 − t22 ≥ 0 . (4.8)
This could be used complementary with (4.6) to form a set of conditions for keeping the
minimum thickness everwere within the considered domain. However, the search space
would be unnecessarily reduced. The maximum possible constrained search space is main-
tained by first checking with the discriminat whether real extremas exist. Then, the
position of the maximum
1
q
ξmax = − 2
t2 − t2 − 3t1 t3 (4.9)
3t3
is substituted in (4.5). This result is then used as complementary condition (4.6)
g4 = − tξmax ≤ 0 . (4.10)
assembly of prepreg sheets each of which is characterized by its instance in the lay-up
sequence, position in space, shape, and size. These parameters are geometrical patch
parameters. Each prepreg sheet constitutes a patch and is in addition assigned a choice of
material and orientation of the principal material axes, which are internal material patch
parameters. For each point on the shell-type structure, the local laminate construction is
determined by the geometrical and internal material patch parameters. Of the complete
set of parameters, the lay-up sequence and the choice of material are discrete, rendering
the parameterization suitable for Evolutionary Algorithms rather than for Mathematical
Programming optimization engines.
The composite boat hull problem is multi-objective where stiffness and prize are in conflict
with each other. This assigns a crucial role to the choice and placement of materials which
can be chosen from a range of inexpensive low-stiffness glass-fiber and expensive high-
stiffness carbon-fiber reinforced plastics. A patch pattern on the boat hull is indicated in
Fig. 1.6.
xui
Figure 4.11: CAD features: lower plate, upper plate and a rib
plate is bounded through a planar functional face at the bottom and through assembled
face segments at the top. These segments are defined through equidistant sampling points
defining six optimization variables tl1 ...tl6 of the lower plate. Additionally, the edges of the
top faces are chamfered. The upper plate is defined through two sets of sampling points,
i.e. six equidistant height parameters hu1 ...hu6 defining the bottom face and six thickness
variables tu1 ...tu6 defining the top face of the upper plate. Again, the edges of these faces
are chamfered. Finally, each rib is defined through three optimization variables: a lower
position xli , an upper position xui , and a thickness tri . For all these optimization parame-
ters a range and a step size σ for Gaussian mutation or initialization are assigned as listed
in Table 4.1. The genotype for the Evolutionary Algorithm is defined as follows. The ribs
are chosen as CAD features to be optimized; one rib is defined by three parameters that
are represented in one gene. For the lower and upper plate sampling positions are chosen
as genes. That means for the lower plate each thickness represents a gene, and for the
upper plate the height and the thickness at a sampling position form a gene. This leads
to the following genotype, where each gene is marked through accolades:
{tl1 } {hu1 , tu1 } {tl2 } {hu2 , tu2 } {tl3 } {hu3 , tu3 } {tl4 } {hu4 , tu4 } {tl5 } {hu5 , tu5 } {tl6 }
{hu6 , tu6 } {xl1 , xu1 , tr1 } {xl2 , xu2 , tr2 } {xl3 , xu3 , tr3 } {xl4 , xu4 , tr4 }.
This representation implicitly fulfills the manufacturing requirement, i.e. extrusion mold-
Table 4.1: Ranges and step sizes for the optimization variables.
Parameters Range Step
[mm] [mm]
tl1 ...tl6 [1, 7] 1
hu1 ...hu6 [1, 33] 4
tu1 ...tu6 [1, 7] 1
xl1 ...xl4 [0.1, 73.3] 8
xu1 ...xu4 [0.1, 73.3] 8
tr1 ...tr4 [1, 7] 1
ing. For the end plate, the holes for medium flow and the tension bolts are made in a
subsequent machining process, and the global vertical edges are rounded. Since material
can only be removed from the extruded base block, planar horizontal faces have to be
machined into the upper plate for a well defined load introduction from the bolts to the
end plate. The position of these contact faces must be adapted to the varying slope and
curvature of the top face of the upper plate, whereas for some solutions this face even can
be split in two subregions.
K n o w -H o w / C o n s tr a in t D e n s ity
variables. The horizontal axis indicates the measure of constraint density, invested pre-
existing know-how, or parameterization sophistication. The indicated example problems
are placed in the scheme according to the summarizing description below.
Motorcycle Frame The problem was solved by assigning tube dimensions to each of the
15 independent tubes. An additionally posed, and interesting, problem was that only
three different dimensions where to be used but these had to be adjusted optimally.
On the other hand, the design freedom is reduced by prescribing fixed values for the
node positions where the tubes are connected with each other. This can be regarded
as pre-existing knowledge.
Race Car Rim The problem was solved by adjusting about 30 parameter values of CAD
entities. The problem is highly constrained as explained in Section 3.7.2.
Fuel-Cell-Stack End Plate Since four ribs where foreseen for the extruded-profile de-
sign, there are 27 free CAD parameters used in the optimization. Thickness values
must be larger than the minimum value that can be realized by the extrusion mold-
ing technique. The upper plate must be in a position above the lower plate. The four
ribs must be positioned within the geometric design space of the structure. Apart
from these obvious constraints that must be imposed to guarantee meaningful de-
sign solutions, there is much freedom to arrive at an optimum design solution whose
shape was not intuitively expected.
to use such different parameterizations? For both problems seem to be very similar re-
garding the structural models, the objectives, the formulation of the objective functions,
and the side constraints.
Both mechanical problems are rotationally symmetric and use in fact the same differen-
tial equations and the same quadratic eight-node Serendipity finite-elements to create the
structural stiffness matrix.
The objectives are to maximize strength and the basic idea behind reaching these ob-
jectives is to make the distribution of stresses, within the regions of interest, as even as
possible. So both objective functions globally penalize the sum of the quadratic (if n = 0)
deviations of the local stresses from the mean stress.
The two problems have also in common that thickness values must be positive, yielding
basically identical side constraints.
Only the loading sets the two problems apart. The flywheel is loaded by inertia body forces
in the radial direction due to the rotation while the onsert is loaded by some external force
in the axial direction.
Anticipating the study of the behavior of the two optimization programs, the flywheel
optimization procedure yields, without any problems, shape results such as shown in Fig.
1.1(b). On the other hand, an early version of the onsert program, using the same parame-
terization as the flywheel program, failed to produce a smooth shape such as shown in Fig.
4.9 [28, 29]. Rather, jagged shapes such as shown in Fig. 4.13 were obtained [30]. The
Figure 4.13: Jagged onsert shape obtained with mesh-dependent analysis variables [30]
parameterization based on global polynomial shape representation makes sure that the ob-
tained shapes are nice in appearance and easy to manufacture although the jagged shape
results are quite correct regarding the successful minimization of the objective function.
Another positive effect of the global shape representation is that the reduced number of
independent optimization variables tends to speed up convergence and reduce the number
of necessary design evaluations.
The question remains as to why the flywheel problem is better-natured than the onsert
problem. The answer lies in its loading situation and how it affects the stress distribution,
forming an objective function topology favoring smooth shape results. The loading of the
flywheel is in the radial direction and, assuming a mesh with only one element row, the
finite elements are loaded in-series. Consequently, when one of these elements becomes
thinner, the radial stress must increase and as it becomes thicker, the stress decreases. But
these stress deviations are immediately penalized by the objective function. Particularly,
the objective function tends to prevent exceedingly small thickness values, having the same
effect as if the side constraint were implemented via a penalty method transformation. On
the other hand, the loading of the onsert with respect to the stressing of the bonding layer
is partially in-parallel and does not produce the same stabilizing effect occurring in the
flywheel problem.
The different design variable models place the two problems at different positions on the
parameterization spectrum. Also, one could say that the program for finding maximum-
strength onsert shapes is useful as a preliminary design tool while the program for finding
maximum-strength flywheel shapes gave some deeper understanding of the mechanical
problem and inspired a new simplified mechanical model [2].
The vector x is called vector of design variables. The objective function as well as the
constraining functions may be linear or nonlinear functions of x. They may be explicit or
implicit in x and may be evaluated by any analytical or numerical techniques we have at
our disposal. However, if mathematical programming is used, it is important that these
functions be continuous and have continuous first derivatives in x. If these conditions
are not satisfied, for instance when discrete-valued variables appear, one must either in-
vent homogenization techniques such as [12] or resort to other methods such as genetic
algorithms.
= > ?
Figure 5.1: Convex (a), concave (b), and neither convex nor concave function (c)
f (x) is convex if for any two points x1 and x2 contained in the set it holds that
4
1 6
3
8
2 4
1 2
1
1
2
0 0 .5
0
0 1 2 3 4
is convex and its feasible region is bounded by a convex constraining function. Then, the
design space is a convex set.
x 2 g ra d f(x (0 ) )
fe a s ib le s e c to r
x (0 )
s
g ra d g (x (0 )
) F (x )= c o n s t.
u s a b le a n d
fe a s ib le s e c to r
u s a b le s e c to r
g (x (0 )
)= 0
x 1
in the figure. The line passing through the reference point x( 0) and being vertical to the
object-function gradient ∇f divides the two-dimensional search space into the usable and
unusable sectors. The line passing through the reference point x( 0) and being vertical to
the constraining-function gradient ∇g divides the feasible and the unfeasible sectors. For
a usable search direction it must be that
∇f (x(0) )s ≤ 0 (5.5)
and for a feasible search direction
∇g(x(0) )s ≤ 0 (5.6)
must hold. The feasible search direction with steepest descent is tangential to the contour
line of the constraining function, g = const. However, if the constraining function is
nonlinear, the feasible search direction of steepest descent may stray into the infeasible
region which is unwanted. Therefore, one can introduce into (5.6) a positive parameter
θ (push-off factor), pushing the search direction away from the tangent at the infeasible
region:
∇g(x(0) )s + θ ≤ 0 (5.7)
The concept of including push-off factors and the resulting algorithms are well described
in Vanderplaat’s textbook [20].
The weight factors λi must now be adjusted so that the search direction s is orthogonal
to the gradients of the constraining functions:
After executing the multiplications and re-ordering the terms one obtains the system of
equations for the weight factors λi
The system of equations can be resolved when the determinant D of the coefficient matrix
is not equal to zero. The determinant D
1.5 0.5
f g L
0.4
1
Lagrangian
0.3
0.5
0.2
0 0.1
‐2 ‐1 0 1 2 3
0
‐0.5 x 0 0.2 0.4 0.6
a b
The second derivative is a matrix called Hessian H and, at the minimum point, it must
be positive definite:
2
∂ f (x) ∂ 2 f (x) ∂ 2 f (x)
. . .
∂x21 ∂x1 ∂x2 ∂x1 ∂xn
∂ 2 f (x) ∂ 2 f (x) ∂ 2 f (x)
∂x ∂x ∂x22
. . . ∂x2 ∂xn
1 2
|H|> 0, . (5.28)
.. .. .. ..
. . . .
2
∂ f (x) ∂ 2 f (x) ∂ 2 f (x)
∂x1 ∂xn ∂x2 ∂xn . . . ∂x2 n
Positive definiteness means that this matrix has all positive eigenvalues. If the gradient is
zero and the Hessian matrix is positive definite for a given x, this insures that the design
is at least a relative minimum, but does not insure that the design is a global minimum.
Only when the objective is known to be convex will (5.27) and (5.28) suffice to identify
the global minimum.
1. x is feasible
gj (x∗ ) ≤ 0 j = 1, m
∗
hk (x ) = 0 k = 1, l
2. λ∗j gj (x∗ ) =0 j = 1, m λ∗j ≥0
m
X l
X
3. ∇x f (x∗ ) + λ∗j ∇x gj (x) + λ∗m+k ∇x hk (x) = 0
j=1 k=1
λ∗j ≥ 0
λ∗m+k unrestricted in sign
The first of the Kuhn-Tucker conditions requires that the minimum of the constrained
function lie within the feasible region. The second condition implies that inactive con-
straints, as the one inactive constraint shown in Fig. 5.5, must not be considered. The
third condition can be understood by the method of feasible directions: At the constrained
minimum point, the length of the feasible search direction becomes zero. This appears
with one constraint, see Fig. 3.3, or several active constraints, see Fig. 5.5.
Ñ B x *
N B x = c o n s t l Ñ C x *
l Ñ C x *
Ñ C x * Ñ B x *
Ñ C x *
C x =
C x =
N
5.7 Duality
Section 5.5 explains that the critical, or stationary, points of the Lagrangian are saddle
points. The saddle-point nature is connected with the max-min problem considered on
Section 5.7.1 and the duality concept introduced with 5.7.2; the presentation of both of
these topics rely on textbook material [20, 32, 33, 34].
For finding a solution to the optimization problem the Lagrangian must be maximized
with respect to λ, or
max L(λ) = max min L(λ, x) (5.30)
λ λ x
- Separable objective with linear and quadratic terms and linear constraints,
- Separable objective and constraints both with linear and quadratic terms, and
Separable Objective with Linear and Quadratic Terms and Linear Constraints
The primal problem is
n
X 1
f (x) = ai xi + bi x2i (5.45)
2
i=1
n
X
g(x) = cji xi j = 1, m (5.46)
i=1
Referring to the result (5.43) the following function is minimized with respect to x:
m
1 X
min ai xi + bi x2i + λj cji xi (5.47)
xi 2
j=1
At the minimum the derivative with respect xi must vanish which condition yields
m
X
ai + bi xi + λj cji = 0 (5.48)
j=1
In case of bi being equal to zero, xi is set to its lower bound if the numerator in (5.49)
is positive and vice versa. The result is to be inserted into (5.43) which is then to be
maximized with respect to the dual variables λj .
Separable Objective and Constraints Both with Linear and Quadratic Terms
n
X 1 2
f (x) = ai xi + bi xi (5.51)
2
i=1
n
X 1
gj (x) = cji xi + dji x2i j = 1, m (5.52)
2
i=1
Inserting the objective and constraints into (5.43) gives the sub problem
m
1 X 1
min ai xi + bi x2i + λj cji xi + dji x2i (5.53)
xi 2 2
j=1
There are many different methods from which optimization algorithms are derived. Müller
provides in her recent dissertation [1] a concise scheme for dividing the methods. The direct
methods require only the values of the objective function itself and the indirect methods
need additional information in terms of first or higher derivatives. Complete evaluation is
an extremely simple and inefficient direct method devoid of any strategy to utilize topol-
ogy information from previous function evaluations for improving convergence. The more
sophisticated direct methods are further subdivided into stochastic and deterministic ones.
Stochastic methods use random numbers to identify new design solutions. The stochastic
methods include, for instance, evolutionary algorithms some of which are explained in
chapter 7.
The deterministic methods rely on some information on the local objective function topol-
ogy to identify more systematically promising points in variable space. They include the
simplex search method or Powell’s conjugate direction method, explained in sections 6.1
and 6.6, respectively.
The indirect methods are all deterministic. One distinguishes there between first-order
methods, needing first-derivative information, and second-order methods, needing in ad-
dition second-derivative information. The first-order methods include the methods of
steepest descent (section 6.2) or the nonlinear conjugated gradient method of Fletcher
and Reeves (section 6.5). The second-order methods include the Newton’s method (sec-
tion 6.4) or surface-response methods (section 6.7).
In the literature [18], the deterministic methods, or the algorithms based on them, are cus-
tomarily called mathematical programming and the direct deterministic methods go under
the label zeroth-order methods of mathematical programming. The illustration in Fig. 5.7
adopts parts from [1] but also indicates the realm of mathematical programming. Apart
from the method order indicated in Fig. 5.7, search methods are characterized by the
model order which is is the order of the Taylor series expansion that approximates the ob-
jective function. Fig. 5.8 arranges the methods in a matrix whose rows and columns stand
for the model and method orders, respectively. Search methods based on quadratic models
of the objective function converge very quickly on quadratic functionals. For instance, the
Newton and the surface-response methods converge in one step on quadratic functionals
P a r a m e te r O p tim iz a tio n T e c h n iq u e s
D ir e c t In d ir e c t
S to c h a s tic D e te r m in is tic U s e f, f U s e f, f, 2
f
- S im p le x M e th o d
- E v o lu tio n a r y - C a u c h y 's M e t h o d
- P o w e ll's M e t h o d - N e w t o n 's M e t h o d
C o m p u ta tio n - F le tc h e r R e e v e s M .
- R e s p o n s e -S u rfa c e
0 th
o rd e r 1 s t o rd e r 2 n d
o rd e r
M a th e m a tic a l P r o g r a m m in g
- R e s p o n s e -S u r fa c e M . - F le tc h e r -R e e v e s C G M . - N e w to n M e th o d
2
- P o w e ll's C D M e t h o d - V a r . M e tr ic M e th o d - Q u a s i-N e w to n M .
r d e r
1 - C a u c h y M e th o d
O
o d e l
M
0 - S im p le x S e a r c h M e th o d
0 1 2
M e t h o d O r d e r
Figure 5.8: Search methods arranged after model and method orders
although the methods are of order two and zero, respectively. The lower triangle of that
matrix is empty because model order limits method order.
x * P o in t M e s h
x 2
x 2
o p t
f ( x ,y ) = c o n s t.
x 1
o p t x 1
He uses this to underline the importance of the much more efficient mathematical pro-
gramming algorithms.
Mathematical Programming (MP) is a rather generic term and the algorithms embraced by
it assume that in general the objective and the constraining functions are continuous and
at least twice differentiable [18]. Moreover, MP assumes convex objective functions. Given
a non-convex objective and a starting point, MP will generally obtain a local minimum
and an additional effort, for instance by employing a multi-start technique, is needed to
find other local minima with smaller function values, and, hopefully, the global minimum
too. Some mathematical programming methods are regarded as highly efficient but it is
also important to note that the efficiency of these methods depends very much on the
actual objective function. Therefore, general statements, claiming one certain method to
be always better than another certain method, are usually untenable.
The terms Sequential Linear Programming and Sequential Quadratic Programming are
used to characterize solution methods for nonlinear constrained optimization problems
which use linear and quadratic approximations, respectively, to the problems. The follow-
ing definitions have been taken from Vanderplaats’ textbook [20].
problem.
Because of the approximation both methods require iterations and each ap-
proximation is called a quadratic sub problem. If all constraints are equality
constraints, the search method remains identical to that of the Newton method.
If the problem includes inequality conditions, the search becomes similar to the
modified Newton method but a nonlinear line search along may be replaced
with a direct step to activate the closest and as of yet inactive constraining
equation.
Fig. 5.7 illustrates how the methods of mathematical programming are divided into the
categories direct methods, gradient-based methods, and second-order methods. Zeroth-
order, or direct, methods require only objective-function values. First-order, or gradient-
based, methods require in addition the gradient (first derivative). Second-order methods
require also the Hesse matrix (second derivative).
Some methods are derived upon the model of a quadratic function, and thus have a the-
oretical basis. The textbook by Reklaitis, Ravindran, and Ragsdell [23] gives two reasons
for choosing a quadratic model:
1. It is the simplest type of nonlinear function to minimize, and hence any general
technique must work well on a quadratic if it is to have any success with a general
function
The direct method due to Powell [35] is also based on a quadratic model. Cauchy’s method
of steepest descent [36] is a gradient method which is based on a linear model: improved
points are suspected in directions along the most local descent. Since no assumption is
made how that descent property might change at some distance away from the reference
point, the method is often less efficient than those based on quadratic models. Its one
advantage over more sophisticated strategies lies in its descent property.
C e n tr o id o f x (2 )
a n d x (3 )
x (1 )
x (1 ) x (4 )
( a ) o ld s im p le x (a ) n e w s im p le x
x (3 )
x (3 )
x (1 )
, x (2 )
, x (3 )
x (2 )
, x (3 )
, x (4 )
The method begins by setting up a regular simplex in the space of the independent vari-
ables and evaluating the function at each vertex. The vertex with highest functional value
is located. This vertex is then reflected through the centroid to generate a new point,
which is used to complete the next simplex. As long as the function values obtained at
each new point decrease monotonously, the iterations move along until either the min-
imum point is straddled or the iterations begin to cycle between two or more vertices.
The minimum point is straddled when always the same vertex is reflected back and forth
through the iterations. These situations are resolved using the following three rules:
1. Minimum ”Straddled”
If the vertex with the highest function value was generated in the previous iteration,
then choose instead the vertex with the next highest value for reflection.
2. Cycling
If a given vertex remains unchanged for a M iterations, reduce the size of the simplex
by some factor. Set up a new simplex with the currently lowest point as the base
point. Spendley et al. suggest that M be predicted via M = 1.65N + 0.05N 2 , where
N is the problem dimension and M is rounded to the nearest integer. This rule
requires the specification of a reduction factor.
3. Termination Criterion
The search is terminated when the simplex gets small enough or else if the standard
deviation of the function values at the vertices gets small enough. This rule requires
the specification of a termination parameter.
Apart from the objective evaluations, there are only two types of calculation required for
the simplex search algorithm: (1) generation of simplex at a given base point x(0) and
scale factor α in N -dimensional space and (2) calculation of the reflected point. For a
given base point the other vertices of the simplex are calculated from
(0)
xj + δ 1 if j = i
(i)
xj = (6.3)
(0)
xj + δ 2 if j 6= i
for i and j = 1, 2, ..., N . The increments δ1 and δ2 , which depend only on N and the
selected scale factor α, are calculated from
h√ i
δ1 = √ −1 α
N +1+N
N 2
h√ i . (6.4)
N +1−1
δ2 = √
N 2
α
The design of a simplex in two dimensions is illustrated in Fig. 6.2(a). Suppose x(j) is the
d 2 d 1- d 2 N
x 4
x 3
a d 1- d 2
d 1 a x c
x c-x 1
x 2
d x
a 2 1
d 1 N
(a) (b)
Figure 6.2: Design Principle (a) and Reflection (b) of a Simplex in Two Dimensions
point to be reflected as illustrated in Fig. 6.2(b). Then the centroid of the remaining N
points is
N
1 X (i)
xc = x , i 6= j. (6.5)
N
i=0
Choosing λ = 2 will yield the new vertex point so that the regularity of the simplex is
retained. Thus,
(j)
x(j)
new = 2xc − xold . (6.7)
s = −∇f (6.8)
Thus, the direction of steepest descent is a useful search direction along which to locate
points with lower function values. Moving along the search direction will initially obtain
ever smaller function values but after reaching a minimum the values will, from there on,
again increase. That minimum point along the search direction can be used as a new
reference point. This generates the sequence
xk+1 = xk + αk sk (6.9)
where the αk are obtained by line searches (see section 6.8). The slope f 0 is the component
of the gradient in the normalized search direction,
sk
f 0 = (∇f (xk+1 ))T (6.10)
|sk |
It vanishes at the minimum point along the search direction because there the search
direction and the gradient are perpendicular to each other. Therefore, any two consecutive
line searches are perpendicular to each other as Fig. 6.3 illustrates. However, exact
values of the step-length αk can be calculated only when the objective is a quadratic
functional. Even then the sequence 6.9 is not guaranteed to yield the absolute minimum
of a convex function in a finite number of iteration steps - the zigzagging course indicated
in Fig. 6.3 could be continued forever. Moreover, the step-length can only be numerically
estimated and very accurate estimations may cost many function evaluations, further
reducing the efficiency of the sequence 6.9. Therefore, using in-exact line searches with
reduced numerical effort have the effect that two consecutive search directions are no
longer perpendicular to each other but the overall efficiency may increase as long as any
new point generated by the line search has a lower function value than the reference point.
f (x) = a0 + a1 x + a2 x2 , (6.11)
f 0 (x) = a1 + 2a2 x
. (6.12)
f 00 (x) = + 2a2
Setting the first derivative equal to zero determines the point xe where the function value
is an extremum:
a1
xe = − . (6.13)
2a2
If the value of the second derivative, or a2 , is positive, the extremum is a minimum. When
the minimum point is considered an optimum, the two conditions are called optimality
conditions.
Next we generalize the concept to quadratic functions depending on several variables xi
that are arranged in the variables vector x,
where the function f and the constant coefficient p are scalars, the variables vector x and
the coefficient p of the linear term are vectors, and the coefficient P is a matrix. The
extreme point xe in variables space is identified by setting the first derivative of (6.14)
equal to zero:
∂f
∇f (x) = = p + 2Px = 0. (6.15)
∂x
The first derivative of a function of several variables, or variables vector, is a vector called
gradient ∇f . From (6.15) follows the extremum point xe
1
xe = − P−1 p. (6.16)
2
The extreme point is a minimum if for any arbitrary vector v with non-zero length it
holds:
vT Pv > 0. (6.17)
Setting the gradient of the Tailor-series approximation fˆ(x0 + 4x) equal to zero yields the
distance 4x from the reference point x0 to the extremum point xapp e :
xapp
e = x0 + 4x = x0 − H−1 ∇f (x0 ). (6.22)
The extremum is a minimum if, for any arbitrary vector v, it holds that vT Hv > 0.
The extremum point of the Tailor-series approximation is generally not identical with the
extremum of the objective function itself. Under certain conditions it will be closer to
the true extremum than the reference point. It can then replace the reference point and
applying (6.22) again yields an improved estimate so that the continued sequence,
converges to the true extremum. It is understood that the gradient and the Hesse matrix
are evaluated at the point xk . This procedure is called the Newton Method.
The original objective function may not be well approximated by its quadratic Taylor-
Series expansion. This is likely to occur when the reference point is too far away from
the minimum point of the original objective. Then, the calculated extremum point of the
approximating function may be farther away from the true minimum than the reference
point so that the sequence does not converge. Introducing a variable step length αk leads
to the so-called Modified Newton Method,
The step length is not beforehand known and must be determined by a so-called line
search. Line search methods are explained in Section 6.7.
In cases the Newton and Modified Newton Methods may be very efficient. However, in
structural optimization they enjoy only limited popularity. This is because the derivatives
can only be extracted numerically by using the methods explained in Section 6.8. The
number of function evaluations to calculate the Hesse matrix scales quadratically with the
number of design variables so that for a large number of design variables the numerical
effort becomes too high. The methods presented in the following sections avoid the need
for explicit calculation of the Hesse matrix.
s0 T Cs1 = 0, (6.25)
s (0 )
r (1 ) = s (1 )
x (1 )
s (0 )
x x
s (1 )
x (1 )
r (1 )
usually leads to nonlinear objective functions so that the algorithm that is well-developed
for quadratic functionals cannot be applied. Then, the Method of Conjugated Directions
developed by Fletcher and Reeves [39] can be used. It generates new search directions
that are approximately C-conjugated to the previous one from the information of previous
line searches after the formula
|∇f k |2 k−1
sk = −∇f k + s . (6.27)
|∇f k−1 |2
The sequence requires memorizing the previous search-direction and gradient vectors. At
the first step of the sequence, these quantities are not known so that the first step is the
same as for the Method of Steepest Descent.
In deriving the sequence (6.27) one assumes a quadratic functional (6.26) so that the
gradient at a given point x is given by
∇f = Cx. (6.28)
The search starts at a point x0 and the initially chosen search directions is that of steepest
descent at this point,
The vector of optimization variables along the search direction depends on a factor λ,
That factor is chosen to minimize the value of f along the search direction. At that
minimum point the derivative of f with respect to λ must vanish. This condition is used
along with (6.29) and (6.30) to derive the minimizing value of λ:
f,λ = xT Cx,λ = 0
We obtain the minimum point along the search direction s(0) by inserting the just derived
minimizing value of λ in (6.30) and call it x(1) . The gradient at this point,
T
(1) (0) (0) (0) ∇f (0) ∇f (0)
∇f = C(x + λs ) = ∇f − T
C∇f (0) , (6.32)
∇f (0) C∇f (0)
is naturally orthogonal to ∇f (0) . Guided by the pre-existing knowledge, that the mini-
mum of quadratic functionals is obtained after a finite number of iteration steps if two
consecutive search directions fulfill C-conjugacy, we require that the form
s(0)T Cs(1) = 0
is derived by analyzing the scalar product of the gradient ∇f (1) with itself:
T
∇f (0) ∇f (0)
∇f (1)T ∇f (1) = ∇f (1)T ∇f (0) − T
∇f (1)T C∇f (0)
∇f (0) C∇f (0)
T
∇f (0) ∇f (0)
= − T
∇f (1)T C∇f (0) . (6.35)
∇f (0) C∇f (0)
T
∇f (1) C∇f (0) ∇f (1)T ∇f (1)
→ T
= −
∇f (0) C∇f (0) ∇f (0)T ∇f (0)
∇f (1)T ∇f (1)
β= (6.36)
∇f (0)T ∇f (0)
Generalizing the indices 0 and 1 to k − 1 and k, respectively, gives the formula (6.27)
found by Fletcher and Reeves [39]. The derivation is based on a quadratic functional and
the formula will result in a better choice of search directions than Cauchy’s Method if the
objective function is dominated by the quadratic terms of its Taylor series representation.
The quadratic terms tend to overwhelm the other terms with decreasing distance to the
minimum point.
If many line searches must be performed within a region of the objective where it deviates
much from a quadratic functional, the search directions may lose their descent property
so that it becomes useful to occasionally restart the method by setting β = 0. The Polak-
Ribière method [40],
∇f (i)T ∇f (i) − ∇f (i−1)
PR
βi = , (6.38)
∇f (i−1)T ∇f (i−1)
allows for re-starting the search through periodically setting βi to zero, thereby offering
potential performance enhancement for non-quadratic functions compared to Fletcher-
Reeves. The underlying mechanism can be explained as follows [41]. Given a quadratic
objective function, both strategies coincide (∇f (i)T ∇f (i−1) = 0). The presence of non-
quadratic terms, potentially combined with insufficiently accurate line searches, is likely
to induce conjugacy loss as the search progresses. This manifests itself in the determined
search direction being nearly orthogonal to ∇f (i) , which is accompanied by:
Hence βi becomes almost zero, directing the search towards the steepest descent:
Q(x) = xT Cx (6.41)
is generally not a sum of perfect squares unless all off-diagonal elements are zero. The
process of transforming it into a sum of perfect squares is equivalent to finding a transfor-
mation matrix T,
x = Tz, (6.42)
so that the functional becomes a sum of perfect squares in terms of the transformed
variables z and the diagonal matrix D,
We realize that the columns tj of the transformation matrix T give a new set of coordinates
that, because they diagonalize the quadratic C, correspond to its principal axes:
x = Tz = t1 z1 + t2 z2 + ... + tN zN . (6.44)
The new coordinates tj are called conjugate directions and the remaining problem is how
to calculate such a set of conjugate directions.
The basic approach to finding the conjugate directions is based on the parallel subspace
property of a quadratic function. The parallel subspace property is explained by using
the two-dimensional example illustrated in Fig. 6.5. Consider two arbitrary but distinct
points x(1) and x(2) and a direction vector d. Let y(1) be the point corresponding to the
minimum of Q(x(1) + λd) and (y)2 be the solution to Q(x(2) + λd). Then the direction
(y(2) − y(1) ) is C conjugate to d. From Fig. 6.5 it is apparent that two line line searches
determine the points y(1) and y(2) , establishing the set of C conjugate directions d and
(y(2) −y(1) ), and a third line search with reference point y(1) or y(2) and along the direction
(y(2) −y(1) ) finds the minimum point of the quadratic. All this is achieved without gradient
calculation.
The foregoing can be extended to give a elucidating definition of conjugacy that is taken
from the textbook [23]. Given an N × N symmetric matrix C, the directions s(1) , s(2) , ...
, s(r) , r ≤ N are called C conjugate if the directions are linearly independent, and
The points along the direction d from x(1) depend on the single variable λ,
x = x1 + λd. (6.47)
N
O
N
O
The minimum of q along the line defined by (6.47) is obtained by finding λ∗ such that the
derivative of q with respect to λ is zero:
∂q ∂q ∂x
= = (bT + xT C)d. (6.48)
∂λ ∂x ∂λ
Accordingly, the directions (y(2) − y(1) ) and d are C conjugate, and the parallel subspace
property of quadratic functions has been verified.
It remains to develop the actual minimization algorithm. The parallel subspace property
has been explained by using two starting points and a direction. Having to generate a
number of starting points is, however, not elegant from a computational point of view.
Therefore we construct a way to find C conjugate directions and the minimum by using
only one starting point. To this end we employ the coordinate unit vectors e(1) , e(2) ,
..., e(N ) . We consider, as before, the two-dimensional case and let e(1) = [1, 0]T and
e(2) = [0, 1]T . We use a starting point x(0) and calculate λ(0) so that f (x(0) + λ(0) e(1) ) is
minimized. This gives us the new point x(1)
x1 = x0 + λ0 e1 . (6.52)
From the new point we start a line search in the direction of e(2) so that f (x(1) + λ(1) e(2) )
is minimized which obtains a second new point
x2 = x1 + λ1 e2 . (6.53)
(0 ) ( 2 ) (1 ) (1 ) ( 2 ) ( 2 ) (3 ) (3 ) (1 )
I = A I = A I = A I = N - N
( 4 ) *
N (3 ) N = N
(1 ) ( 2 )
N N
( 0 )
N
Next, we use the first coordinate unit vector again for a line search to calculate λ(2) so
that f (x(2) + λ(2) e(1) ) is minimized and let
x3 = x2 + λ2 e1 . (6.54)
Then, the directions (x(3) − x(1) ) and e(1) will be conjugate as Fig. 6.6 illustrates. The
construction can be extended from 2 to N dimensions and yields Powell’s Conjugate
Direction Method :
1. Define the starting point x(0) and a set of N linearly independent directions, possibly
s(i) = e(i) , i = 1, 2, 3, ..., N .
2. Minimize along the N + 1 directions, using the previous minimum to begin the next
search and letting s(N ) be the first and last searched
3. Form the new conjugate direction using the extended parallel subspace property
4. Delete s(1) and replace it with s(2) , and so on. Place the new conjugate direction in
s(N ) . Goto step 2.
Figure 6.7 illustrates the path taken by Powell’s method through the variables space of a
non-quadratic function.
Figure 6.7: Path Taken by Powell’s Method Algorithm Through Non-Quadratic Function
Space. Source: [20]
The lower and upper bounds xL,i and xU,i are initially imposed by the designer. To reduce
the number of costly function evaluations, a response surface model, also called a surrogate,
is built for the objective and possibly also for the constraining functions. The search for
the optimum design can then be partly make use of the surrogate functions and then the
original optimization problem is replaced by
n o
minn f˜(x)|g̃i ≤ 0, h̃k = 0 , j = 1, ..., J, k = 1, ..., K, xL,i ≤ xi ≤ xU,i (6.56)
x∈R
where the tilde symbol indicates that the respective surrogates are meant instead of the
original functions.
The response surface methodology was originally developed for constructing empirical
response functions based on physical experiments. When the physical experiments are
replaced by numerical function evaluations, the methodology can be used to find minimum
points of the function. For that purpose quadratic polynomials are obviously suited.
They are the simplest functions with a minimum point and can easily be constructed and
evaluated. The methodology is based on the elements:
1. Construct the response surface model for the given number of variables
3. Find response surface model parameters from the supporting point evaluations
4. Find minimum point of the response surface model, or surrogate objective function
Response surface approximations shift the computational burden from the optimization
problem to the problem of constructing the approximations, and accommodate the use of
detailed analysis techniques without the need of derivative information, [44]. Additionally,
response surface approximations filter out numerical noise inherent to most numerical
analysis procedures, by providing a smooth approximate response function, and simplify
the integration of the analysis and the optimization codes.
where c contains the polynomial powers and pk is the k th supporting point. The parameters
ai of Π2 can be adjusted to fit the objective function values fk at the supporting points
pk . The number ns of supporting points spanning the search space must at least equal the
number nc of parameters or coefficients in Π2 . That number depends on the number nx
of optimization variables:
1
nc = 1 + nx + [nx (nx + 1)] . (6.59)
2
If the number ns of supporting points equals the number nc of the parameters of the
quadratic approximation, the parameters a are obtained by inverting (6.58):
a = C−T f . (6.60)
It is also possible to use a higher number of supporting points for constructing the quadratic
model. In general the approximation cannot match the function values for all of the
supporting points if their number is higher than the number of terms in the polynomial,
ns > nc , so that a residual r remains:
CT a − f = r (6.61)
A best-fit approximation follows from minimizing the inner norm rT r of the residual r
with respect to the polynomial parameters a:
∂ rT r
−1
=0 ⇒ a = CCT C f. (6.62)
∂a
If there are more supporting points ns than model parameters nc , the coefficient matrix
C has dimensions ns × nc .
Π2 = p + xT p + xT Px, (6.63)
where the parameters a are contained in the scalar p, the vector p, and the matrix P. As
a necessary condition, the gradient must vanish at an extremum,
vT Pv > 0, |v| =
6 0. (6.66)
& '
# $
!
"
Figure 6.8: Placement of supporting points x(1) through x(10) in 3-dimensional variables
space
is the reference for constructing the other supporting points. Assuming the placement
of the points on a regular lattice, indicated in Fig. 6.8, with a spacing D along the
respective coordinate directions, the coordinate variations of the individual supporting
points with respect to the reference point x(1) are given in the matrix presentation of
Table 6.1. The rows correspond to the optimization variables and the columns correspond
to the supporting points x(2) through x(10) . The coordinates of the points can generally
2 3 4 5 6 7 8 9 10
1 −D D 0 0 D 0 0 D 0
2 0 0 −D D D 0 0 0 D
3 0 0 0 0 0 −D D D D
be identified for nx variables by the scheme developed in the following. The points are
identified successively from variable 1 through nx . For the ith variable, a number of i + 1
supporting points must be added to the set of nc = nc (i − 1) (6.59) existing points. The
first two of the new points span the direction of the ith variable:
(nc +1) 1 0, k 6= nx + 1
xk = xk − (6.67)
D, k = nx + 1
and
(n +2) 0, k 6= nc + 1
xk c = x1k + . (6.68)
D, k = nc + 1
The remaining i − 1 points are variations of the new point x(nc +2) (6.68):
(nc +2+l) (nc +2) 0, k 6= l
xk = xk + , l = 1, i − 1 (6.69)
D, k = l
The supporting points span some appropriate region of the search space. Along with the
function values at each point, they are used to construct the response-surface approxima-
tion.
It is interesting to examine the presented supporting-point arrangement for the NM. The
NM uses a set of supporting points in the vicinity of a reference point to calculate the
gradient and Hessian matrix at that point. Point x(1) then becomes the reference point
at which the gradient and Hessian matrix must be obtained by the method of differences.
The other supporting points must then move very close to the reference point which is
achieved by choosing a very small value for the distance D:
D = << 1. (6.70)
It is interesting to note that the point arrangement allows exact calculation of both first
and second derivatives of quadratic polynomials independent of the value of . This is
illustrated on the three dimensional case, where the points 1 through 10 are arranged as
shown in Fig. 6.8. Let the gradient and Hesse matrix be calculated from
f3 − f2
1
∇f = f5 − f4 (6.71)
2
f8 − f7
and
1 a b c
∇∇f = Hij + H ij + Hij , (6.72)
2
and Hc is formed in terms of the points shifted in away from the coordinate axes in other
directions:
0 f6 f9
Hijc =
f6 0 f10
. (6.75)
f9 f10 0
Obviously, the gradient of quadratic functions is calculated exactly because (6.71) is the
central differences method. It can be seen from (6.73), (6.74), and (6.75) that (6.72)
produces a symmetric Hessian matrix. Moreover, it can be shown that the evaluations
reproduce the symmetric version of the matrix P, multiplied by two:
The presented point arrangement and evaluation scheme calculate for quadratic functions
f := Rn → R and regardless of the value of the epsilon parameter not only the Hessian
but also the gradient exactly.
The RSM evaluation scheme obtaining the approximation parameters p, p, and P, ex-
plained in Subsections 6.1 and 6.2, can therefore be replaced by p = f (x(1) , p = ∇f , and
P = 12 ∇∇f where ∇f , and ∇∇f are obtained by (6.71) and (6.72), respectively for arbi-
trary values of D. For quadratic functions, both methods obtain the same approximation
parameters. Therefore, the evaluation scheme of the Newton method becomes identical
with that of the RSM method at the limit of D approaching very small values D → .
The entries of the gradient vector, where (6.71) represents the example for three variables,
are calculated by the general scheme
fnc (k−1)+2 − fnc (k−1)+1
∇fk = . (6.77)
2
The entries of the matrices Ha , Hb , and Hc follow from
(
2f1 , j = i
Hija = − , (6.78)
−f1 , j 6= i
(
(fnc (i−1) + fnc (j−1) ) , j = i
Hijb = , (6.79)
−(fnc (i−1)+2 + fnc (j−1)+2 ) , j =
6 i
and (
0 , j=i
Hijc = . (6.80)
fnc (j)−i+1 , j 6= i
Figure 6.9: Updated supporting point sets around successful minimum point estimates
by (6.65) does not yield a smaller function value than any of the supporting points, it is
discarded and the strategy considers two cases.
In the one case, the current reference point has a smaller function value than any other
point of the set, thus remaining the best estimate of the objective minimum point. Then,
the reference point is kept and convergence is facilitated by shrinking the√set of supporting
points about the reference point. The shrinking factor can be chosen 1/ 2, for instance.
In the other case, a point of the set other than the referenced point appears as best
estimate of the objective minimum point. The whole point set is then shifted by moving
the reference point to the point with smallest function value. Instead of using this mental
picture, one can also say that say that the point with smallest function value becomes the
new reference point and the distance value defining the extension of the region covered by
the point set is kept constant. This keeping constant of the distance value is motivated
by the opportunity to reduce the number of function evaluations of the new point set
as some points of the new set will then have the same coordinates as other points of the
Figure
√ 6.10: Shrinking of supporting point region around the reference point by a factor
of 1/ 2
previous set whose function values are already known. Fig. 6.11 illustrates the shifts to the
points 2 through 6, respectively. The figure obtains that in each case three points of the
shifted sets coincide with previous points so that three function evaluations can be saved,
making the optimization process more time-efficient. In general, however, when moving
5 6
2 1 3
c u r r e n t p o in t s e t s h ifte d to p o in t 6
4
s h ifte d to p o in t 5 s h ifte d to p o in t 3
s h ifte d to p o in t 2 s h ifte d to p o in t 4
Figure 6.11: Supporting point set shifted to center around the respective points with
smallest function value
the supporting point set so that the reference point of the shifted set coincides with one of
the other points in the original position, the total number of coinciding points depends on
the point to which the reference point is shifted. Table 6.2 shows the number of coinciding
2 3 4 5 6 7 8 9 10 11 12 13 14 15
2 2
3 3 3 3 3
4 4 4 4 3 4 4 3 4
5 5 5 5 3 5 5 3 3 5 5 3 3 3
Table 6.2: Number of coinciding points for the different reference point positions of the
shifted sets corresponding to 1-, 2-, 3-, and 4-dimensional variable spaces.
points when the reference point of the new supporting point set coincides with the other
points of the previous set. In the 1-dimensional variables space there are two other points
where the reference point can be shifted and in each case there are two coinciding points.
The situation in the 2-dimensional variables space is illustrated in Fig. 6.11. In the 3-
dimensional variables space with ten supporting points the number of coincidences depends
on which point of the previous point set the reference point of the new point set is located.
The number of coincidences is either 3 or 4. In the 4-dimensional space the number of
coincidences is either 3 or 5. For a set of supporting points corresponding to nx variables,
the average number nc of coinciding points can be shown to be
nx
nc = [2 (nx + 1) + 3 (nx − 1)] (6.82)
np − 1
The minimum, maximum, and average numbers of coinciding points, depending on the
number of variables, are shown in Fig. 6.12. The maximum number increases linearly
Figure 6.12: Maximum, average, and minimum number of coinciding points versus number
of variables
and the minimum number remains constant at three for variables spaces higher than one-
dimensional. The average number shown in the figure is calculated with (6.82). The points
of each shifted set coinciding with previously evaluated points need not be evaluated again.
Fig. 6.13 shows the point set size and its reduction due to coinciding points depending on
the number of variables. The potential computing time savings are rather large when the
number of variables is small but tend to become insignificant at large numbers of variables
as Fig. 6.14 illustrates. The figure plots the point set size over the point set size reduced
Figure 6.13: Points set size and reduction due to coinciding points versus number of
variables
which is a measure of the potential efficiency gain. As the figure shows, a program may
run up to three times faster when searching for the minimum point of only one variable, or
twice as fast when there are two variables, but only insignificantly faster when the numbers
of variables is large.
Figure 6.14: Relative computing time savings due to coinciding points versus number of
variables
$
#
"
!
! " # $
ARSM because the points are easy to construct and they fill the complete reducing design
space.
The quadratic response surface model, or surrogate, is fitted to the data using the least
square method. When constraints are considered, a global optimization algorithm is used
to find the optimum. Following this step, the value of the actual objective function at the
optimum of the surrogate is calculated through an evaluation of the computation-intensive
objective function. If the value of the actual objective function at the surrogate optimum
is better than the values at all other experimental designs, the point is added to the set of
experimental designs for the following iteration because the point represents a promising
search direction.
4x = 4xs̄. (6.86)
At the minimum point the slope f 0 , or the derivative of the objective function f with
respect to the distance variable x, along the search direction s̄ must vanish. The slope can
be calculated as scalar product of the gradient vector and the normalized search direction:
y m + 1
y m
y m + 1
x m -1 x m + 3 x m + 4 x m x m + 5 x m + 2 x m + 1
• The relative change between two consecutive function values is smaller than a preset
value
The golden-section method uses an optimized division of the intervals to obtain the desired
containment of the minimum point with the least possible number of iteration steps. For
deriving the method, consider a function with one single minimum within the normalized
interval xL = 0 and xU = 1, see Fig. Next, two more points x1 and x2 are chosen inside
of the initial interval. The example illustrated in Fig. 6.17 indicates a higher function
value at x1 so that this point yields an improved lower bound. The coordinates of the new
points are chosen such that the interval is reduced by the same fraction. Since generally
y o
y u
y 1
y 2
x u x 2x 1 x o
it is not known which of the points x1 and x2 will become the new bound, the new points
must be chosen symmetrical to the midpoint of the interval:
xU − x2 = x1 − xL . (6.88)
The ratio between the old and the new intervals remains constant at each iterate, if the
points abide to
x1 − xL x2 − x1
= . (6.89)
xU − xL xL − x1
From Fig. 6.17 we have that
x2 = 1 − x1 (6.90)
1 − 2x1
x1 = → x21 − 3x1 + 1 = 0 → x1 = 0.38197. (6.91)
1 − x1
x2
= 1.61803 (6.92)
x1
This ratio has some significance in fine arts, architecture, and philosophy: proportions in
the golden-section ratio are sensed esthetic or harmonious. In the more general case, when
xL 6= 0 and xU 6= 1, new bounds are obtained by the formulae
x1 = xL + τ (xU − xL )
. (6.93)
x2 = xU − τ (xU − xL )
f = a0 + a1 x + a2 x2 . (6.94)
One version uses three supporting points (x1 , f1 ), (x2 , f2 ), and (x3 , f3 ) to determine the
coefficients ai of (6.94) fitting through these points. An illustrative example of this is given
in Fig. 6.18 The coefficients a0 , a1 , and a2 follow from the system of equations
1 x1 x21
a0
f1
2
1 x2 x2 a = f (6.95)
1 2
2
1 x3 x3 a2 f3
1 .8
1 .6
1 .4
1 .2
1
f(x )
0 .8
0 .6
0 .4
0 .2
0
-1 -0 .8 -0 .6 -0 .4 -0 .2 0 0 .2 0 .4 0 .6 0 .8 1
x
If one places value upon a translucent symbolic resolution, the form [23]
f = b0 + b1 (x − x1 ) + b2 (x − x2 )(x − x3 )
, (6.97)
f0 = b1 + b2 (2x − x2 − x3 )
after substituting the supporting points for the variable x, leads to the less strongly coupled
system of equations
a0 = b0 − b1 x1 + b2 x2 x3
a1 = b1 − b2 (x2 + x3 ) . (6.100)
a2 = b2
The extreme xe is given by (6.13). The iteration requires replacement of one of the three
supporting points by the newly found point xe . The function value fe at the new point
follows from evaluating the original objective function. For instance, one can discard the
old supporting point with the highest function value. The algorithm is relatively easy to
implement and works without additional gradient information.
The cubic approximation method uses the model
f = a0 + a1 x + a2 x2 + a3 x3
. (6.101)
f0 = a1 + 2a2 x + 3a3 x2
It needs only two supporting points but at each of these the function value f as well as the
slope f 0 must be calculated. An illustrative example is given in Fig. 6.19. The coefficients
f = b0 + b1 (x − x1 ) + b2 (x − x1 )(x − x2 ) + b3 (x − x1 )2 (x − x2 )
. (6.103)
f 0 = b1 + b2 [(x − x1 ) + (x − x2 )] + b3 [(x − x1 )2 + 2(x − x1 )(x − x2 )]
f1 = f (x1 ) = b0
f2 = f (x2 ) = b0 + b1 4x
. (6.104)
f10 = f 0 (x 1) = b1 − b2 4x
f20 = f 0 (x2 ) = b1 + b2 4x + b3 4x2
The result is
4f 4f − 4x4f 0 4xf10 + 4xf20 − 24f
b0 = f1 , b1 = , b2 = , b 3 = . (6.105)
4x 4x2 4x3
The coefficients ai of the polynomial (6.101) are related with the bi by
a0 = b0 − b1 x1 + b2 x1 x2 − b3 x21 x2
a1 = b1 − b2 (x1 + x2 ) + b3 (2x1 x2 + x21 )
. (6.106)
a2 = b2 − b3 (2x1 + x2 )
a3 = b3
The extreme points then follow from setting the derivative (6.101) equal to zero:
a2 a1
x2e + 2pxe + q = 0, p= , q= (6.107)
3a3 3a3
If the original objective is a quadratic, the approximation can only be quadratic as well so
that the coefficient a3 vanishes. Then xe follows from (6.13). Else the solution of (6.107)
is given by p
−a2 ± a22 − 3a1 a3
xe = (6.108)
3a3
The sign of the root must be chosen correctly. The cubic approximation offers a maximum
of two extreme points the one of which is the desired minimum and the other one is a
maximum. At the minimum point the second derivative is positive:
Combining (6.109) with (6.108) yields after some simplifications the result
q
± a22 − 3a1 a3 ≥ 0 (6.110)
and reveals that by choosing the positive sign of the root in (6.108) we always obtain
the minimum point xe = xmin . Extreme values exist only if the argument of the root
is positive. An unfortunate choice of the supporting points may yield an approximating
polynomial without extreme points. This can be avoided by checking the condition
x ≤ xe → f0 ≤ 0 ∧ x ≥ xe → f 0 ≥ 0. (6.111)
is fulfilled or that the extreme point is bounded between the two supporting points.
Often, the numerical effort, or the programming difficulty, of obtaining derivative infor-
mation lets the cubic approximation method based on four supporting points appear more
attractive. Here, the parameters ai can easily be obtained analytically by using the form
where 4xij = xi − xj . The equations (6.113) are resolved to give the parameters bi :
4f32
b1 = ;
4x32
b0 = f3 − b1 4x31 ,
. (6.114)
f4 − b0 − b1 4x41
b2 = ,
4x42 4x43
f1 − b0 − b2 4x12 4x13
b3 =
4x12 4x13 4x14
Ordering the terms of (6.112) after powers of the variable x yields the form
N N N ! N "
Figure 6.20: Golden section method supporting point on a range bracketing the minimum
point
with the smallest objective function value. A set of three points is now defined that
consequently contains x2 , x3 , and X4 of the four-point set. The lower limit xL of the
narrowed range is set equal to x1 of the three-point set and the upper limit is set to x3
of the three-point set as shown in Fig. 6.21. In the Situation depicted in Fig. 6.21 the
X 4 1 X 3 1= X 4 2 X 3 2= X 4 3 X 3 3= X 4 4
Figure 6.21: Narrowing the bracket by reduction of four supporting points to three
minimum point, in terms of the original four supporting points, is point x4 . Using the
quadratic approximation method for estimating a new minimum point implies a risk that
the new minimum point is estimated outside of the current bracketed range as indicated
in Fig. 6.22. This is not wanted because it has been established before that the minimum
must be inside the bracketed range. Therefore, the integer switch iQ is set equal to zero
so that the golden section method instead of the quadratic approximation methods is
selected. The golden section method uses the reduced range defined by the set of three
X 4 1 X 3 1= X 4 2 X 3 2= X 4 3 X 3 3= X 4 4
supporting points to replace the previous set of four points by a new set of four points.
Based on the situation depicted so far, the result is shown in Fig. 6.23. Characteristic for
the golden section method is that the new four points now divide the new and narrower
range with the same ratios as the previous set divided the previous range. This is achieved
by using the selected three points of the previous set and adding just one new point. The
new point is point x3 of the new set in Figure 6.23. The objective function is evaluated
at the new point and it is determined which point of the set has the minimum function
value found so far. As Fig. 6.23 indicates, point x4 remains the minimum point for the
considered example. The minimum point and function value found by the golden section
method are called xGS and fGS . One could continue iterations on the golden section
method by restarting the narrowing of the bracketed range as explained above. However,
the routine CBS EARCH performs a quadratic approximation as well because the case
of the estimated minimum point lying outside of the bracketed range is only a possibility
but, generally, not a certainty.
The quadratic approximation method with three supporting points is explained in section
6.8.3. Basically, a quadratic polynomial is fitted to the three supporting points of the
current three-point set. As a condition for an extreme point, the spatial derivative of the
polynomial is set equal to zero. The resulting equation is resolved for x. The second
spatial derivative must be positive if x is a minimum point. If it is not positive or if the
minimum point x is calculated outside the bracketed range, the result is disregarded and
c is placed in the middle between the two brackets. In the program, x is named XQA and
the objective function is evaluated and the value is called F QA. Even if the quadratic
approximation as such failed, it is possible that the newly found point in the middle of the
range might be a valid minimum point. Therefore, the results of the golden section and
the quadratic approximation codes are compared with each other.
The step of comparing the golden-section and the quadratic-approximation methods is
skipped when the golden section method has not been used. A criterion is introduced
that characterizes the shape of the objective function within the bracketed range. If
the function is extremely steep, numerical round-off errors may impair the quality of the
quadratic-approximation results. If the function is not terribly steep and the function
value found by the quadratic approximation method is smaller than that found by the
golden section method, the result of the quadratic approximation is chosen for the new
minimum point and the switch IQ is set equal to 1 so the golden section method is not
used anymore in order to save computing time. Referring to the situation depicted in Figs.
6.20 through 6.23 one can expect one or two more iterations on the golden section method
before the program switches to the then more effective quadratic approximation method.
After it has been decided upon which one of the points XGS and XQA determines the
new minimum-point step length XN EW , the new point x0 in variable space is calculated.
If F N EW is smaller than the previous found minimum value F M IN , the minimum point
is update by replacing XM IN with XN EW and F M IN with F N EW . The newly found
step-length value XN EW must now replace one of the points of the existing four-point
sets according to its value with respect to the other points.
The newly found step-length value XN EW and the current three-point produce the new
four-point set. In case the new point is found by quadratic approximation the golden-
section ratio is generally lost. Fig. 6.24 sketches that a newly found point is located
between points two and three of the current three-point set so that it becomes point three
of the new four-point set. Its position deviates from that corresponding with the golden
section ratio, which does not really matter because when the quadratic search starts being
preferred over the golden section method, the latter one is not used anymore during the
remainder of the current line search. The point with the smallest function value has the
: ! : ! : ! !
step length XM IN and the function value F M IN . The point with the second smallest
function value has the step length XSEC and the function value F SEC. The points are
needed for the termination or convergence criteria, respectively. Referring to Fig. 6.24,
the third point would be XM IN and the fourth point would be XSEC.
There are two criteria that must both be fulfilled simultaneously in order to terminate the
line search.
1. Close at the true minimum point, the actual function value FUN must be very close
to the function value F N EW A as estimated by the quadratic approximation curve.
The variable DF carries a relative error that is calculated as the absolute value of
the difference divided by the sum of the two values
2. Close at the true minimum point, the relative change of the step length with respect
to the absolute value is small. The criterion is set up so that smaller values result
as the minimum point moves closer to the center of the bracketed range.
A third criterion terminates the search even if the two described above do not. A perfectly
solved DoD problem corresponds to a zero objective function value. It is therefore safe
to terminate the search when the value of the objective function itself is smaller than the
chosen epsilon value.
If the criteria do not terminate the search, the iterations continue with extracting a new
three-point set from the four-point set shown in Fig. 6.24.
∇x L = ∇x f + λg ∇x g + λh ∇x h ; ∇λ L = g + h (6.116)
If no other means are available, the partial derivatives can be approximated with:
Figure 6.25: Two-dimensional variable space with two linear constraining functions
nomial, and that the method of conjugate gradients by Fletcher and Reeves finds the
constrained minimum exactly within two line searches.
N
D = D =
N
? N
N L
L
u @
u @
N L
Ñ D
Ñ D
Figure 6.26: Two-dimensional variable space with two linear constraining functions
Let xv be a violating point in an infeasible region of the search space. It is then desired
to find a new point x0 where all the violated constraints remain active with zero values.
The vector pointing from xv to x0 is called c and is a linear combination of the gradients
∇gi of the constraining functions at the point xv :
c = µi ∇gi . (6.120)
It is assumed that xv is close enough to the feasible region so that the constraining func-
tions gi can be linearly approximated along c:
∇giT ∇gi
gi = βi gi0 = βi = βi |∇gi |. (6.121)
|∇gi |
In the case of side constraints, the assumption made above is fully justified because the
constraining functions are linear anyway and (6.121) holds exactly. The distances of the
functions (or hyper planes) gi = 0 to the point xv are
gi
βi = . (6.122)
|∇gi |
The constraining function gradients ∇gi will generally not be linearly independent. There-
fore one may not confuse the linear-combination weight factors µi with the respective dis-
tances βi . The following method for constructing vectors di that have components in the
direction of the gradients ∇gi only, i.e. that are completely decoupled from all the other
gradients ∇gj
di T ∇gj = 0, i 6= j, (6.123)
requires that not two or several of the gradients ∇gi are equal to each other. The vectors
di are then constructed from all remaining violated constraining function gradients dj :
After substituting (6.124) in the condition (6.123) and re-ordering one obtains the factors
γij :
i=j → γij = 0
T (6.125)
i 6= j → γij = ∇giT ∇gj
∇g ∇gj j
Of course the correction vector c (6.120) can as well be written in terms of the vectors di :
c = νi di . (6.126)
In contrast to the gradients ∇gi the vectors di are linearly independent from each other
and so the weight factors νi can be calculated from the respective single constraining
functions. The component νi di of the vector di in the direction of the gradient ∇gi must
be equal to the distance βi :
∇gi
νi di = βi . (6.127)
|∇gi |
Solving (6.127) for νi and substituting (6.122) for βi yields
gi
νi = T
. (6.128)
di ∇gi
In a program the first step is to identify the violated constraints and to obtain the respective
gradients. The factors γij are calculated after (6.125) and the linearly independent set of
vectors di is computed after (6.124). Next the factors νi are calculated after (6.128) and
the correction vector c follows from (6.126). In cases where the assumption (6.121), the
constraining functions being linear along c, is only an approximation, the new point must
be checked for remaining constraint violations. If any of the constraining functions have
values greater than zero the procedure must be repeated until all constraints are satisfied.
An illustration to the here presented concept is given in Fig. 6.26. If multiple vectors
appear, they must first be removed from the set. From a set of multiple gradient vectors
∇gi , the one associated with the constraining function with the greatest value should be
kept.
The differences method implies that the structural model must be evaluated N times for
each gradient calculation and is therefore also referred to as a brute-force method. The
expense for one gradient calculation equals roughly N function evaluations. For problems
with a high number of optimization variables the expense for one gradient calculation may
become much greater than one function evaluation, abolishing the efficiency of gradient-
based minimization methods.
Kũ = r (6.130)
that must be solved for the unknown nodal-point displacements, or degrees of freedom, ũ.
The stiffness matrix K is a numerical representation of the behavior of the structure and
is assembled from the stiffness matrices of the individual finite elements into which the
structure domain is divided. The right-hand-side vector r contains the nodal forces derived
from the loads acting upon the structure. Once (6.130) is resolved, or the primary solution
ũ is obtained, any other quantity is calculated as a function of the primary solution, and
possibly model parameters, by the postprocessing step. Having gradient calculation in
mind, we assume that the objective function depends explicitly and implicitly, through
the primary solution, on the optimization variables x, f = f (x, ũ(x)). Therefore, the
sensitivity of the objective with respect to a change of the ith optimization variable is
obtained by the rule of chains:
∂f ∂f ∂ũ
∇fi = + (6.131)
∂xi ∂ũ ∂xi
The postprocessing is most often numerically insignificant if compared with the effort to
∂f ∂f
solve (6.130). Therefore, the quantities ∂x i
and ∂ũ are obtained at low numerical cost.
The cost of calculating the sensitivity of the primary solution ũ upon the optimization
variables x is much reduced by the formula derived in the following.
The first step in deriving the sensitivity formula is to form the partial derivative of (6.130)
with respect to the ith optimization variable xi ,
∂K ∂ũ ∂r
ũ + K = . (6.132)
∂xi ∂xi ∂xi
Resolving (6.132) for the sensitivity of the primary solution upon xi yields the desired
sensitivity formula
∂ũ ∂r ∂K
= K−1 − ũ . (6.133)
∂xi ∂xi ∂xi
On first sight one could conclude that nothing has been gained since, after all, the inversion
of the stiffness matrix is at least numerically equivalent to resolving (6.130). In fact,
however, (6.132) allows to calculate the gradient very efficiently, especially for large models.
Let us realize that the inverse at the reference point, at least in terms of a triangulated
matrix, exists at the point (x0 , f0 ) because the objective has just been evaluated there.
It remains to calculate the terms within the braces. The sensitivity of the nodal forces
upon the optimization variables, which are design parameters in structural optimization,
is easily calculated at practically no cost. Often, if the design is subject to fixed and
concentrated loads, the nodal forces are independent of the design parameters so that the
first term in braces becomes zero. The second term in braces is a product of the sensitivity
of the stiffness matrix with the existing primary solution vector. Often, the change of a
design parameter affects only part of the structural model, or only a small subset of all
the finite elements making up the whole model. It is therefore not necessary to assembly
the complete stiffness matrix. Instead, one may only consider the elements whose shape or
other properties are affected by a change of xi . Then, the numerical effort may be much
less than for a complete stiffness matrix assembly, which is again much less expensive than
a solution of the system equations. Fig. 6.27 illustrates the point. The example refers
Figure 6.27: Shape optimization and sensitivity formula for gradient calculation
to shape optimization where the positions of boundary nodes is variable in the direction
perpendicular to the boundary. The variation of the position of one node affects the shape
of the two adjacent finite elements while all other elements remain unchanged. The stiffness
matrices of the two elements with the reference shapes are subtracted from those with the
modified shapes. The difference is multiplied with the existing displacement vector, where
only the displacements on the nodes of the affected elements need to be considered.
Although the sensitivity formula allows gradient calculation at much lower cost than the
differences method, it may be difficult to implement in existing general-purpose FEM
programs.
• Multi-start method
• Tunneling method
x * p o in t m e s h
x 2
x o p t
2
f(x 1,x 2)= c o n s t.
x 1
o p t x 1
Figure 6.28: Function with minima in considered region and starting-point mesh.
these points are used as starting points for a minimum search by using mathematical
programming. If the search from each starting points leads to the same minimum point,
one may be dealing with a convex function. Operating upon a non-convex function, the
multi-start method will generally obtain several minimum points (if there are several in
the feasible region). The one with the lowest function value may be taken for the global
minimum although there is no proof that the absolute minimum may not have escaped the
search. The probability of finding the global minimum increases with increasing density
of the starting-point mesh. However, a high mesh implies an immense number of starting
points, or optimization processes, when the number of dimensions is high. Covering the
range of each of the N variables with M points yields M N starting points altogether. It
goes without saying that the method is forbidding when it comes to large optimization
problems.
1. Minimization phase: starting from the point x1 , a new local minimum point x∗1
is obtained by mathematical programming
2. Tunneling phase: starting from the local minimum point x∗1 , a different point x2 ,
having the same function value than x∗1 or fulfilling f (x2 ) = f (x∗1 ), is searched
After the tunneling phase follows the return to the minimization phase. With given min-
imum point x∗i a tunneling function t,
F (x) − f (x∗i )
t(x, x∗i , si ) = , (6.134)
[(x − x∗i )T (x − x∗i )]si
is defined. The parameter si must be chosen so that the denominator tends stronger
against zero than the numerator when x → x∗i so that the tunneling function value is
different from zero:
t(x → x∗i , x∗i , si ) 6= 0. (6.135)
On the other hand, the value si should also be such that the tunneling function decreases
with increasing distance from the latest minimum point x∗1 :
If the conditions (6.135) and (6.136) are satisfied, any null of the tunneling function t is a
point where the objective function f has the same value as the reference minimum point.
This new point is a suitable starting point for a new minimum search and guarantees that
the next obtained minimum point has a lower function value than the previous one. The
tunneling function corresponding to the global minimum has no null.
The tunneling method is demonstrated on a polynomial of degree six that depends on only
one variable. The objective function shown in Fig. 6.29(a),
has a minimum point at x∗1 = 1 with a value f (x∗1 ) = 20. For purpose of demonstrating
an unusable tunneling function we choose a value si = 0.5 for the root parameter. Then,
the tunneling function shown in Fig. 6.29(b) is obtained by subtracting 20 from the
polynomial (6.137) and dividing the result by x − 1. Obviously, that tunneling function
does not satisfy the condition (6.135) because it has a null at x∗1 or t(x∗1 ) = 0. However,
choosing si = 0.5 obtains the tunneling function shown in Fig. 6.29(c) by subtracting 20
Figure 6.29: Function with objective (a) and unusable (b) and usable (c) tunneling func-
tions
from the polynomial (6.137) and dividing the result by (x − 1)2 . This tunneling function
is usable because its value at x∗1 is different from zero and it decreases monotonously
from there to the next null at x2 = 2. The null is quickly approximated by employing
the Raphson method. From there, mathematical programming will obtain the other local
minimum point x∗2 ≈ 2.517.
Stochastic Search
A general optimization problem can be characterized by the structure of its search space, its
objective space and its objective function. The optimization variables are parameterized
within the search space (e.g. Rn in a real-valued optimization problem). Candidate
solutions are rated in the objective space (e.g. R for single-objective problems or Rn for
multi-objective problems). The methods from the field of Mathematical Programming are
well suited for continuous, homogeneous search spaces, continuous, scalar objective spaces
and smooth, convex objective functions. These qualities can not always be expected to be
present in an engineering problem at hand, indicating that deterministic algorithms are
inappropriate. Therefore, a lot of research has focused on the development and applications
of stochastic search methods to overcome this limitation. This chapter gives an overview
over the field of stochastic optimization. Some of the concepts presented earlier in this
lecture are directly transferable to stochastic search. However, it should be pointed out
that the separation of the parameterization concept (or in general the representation) and
the search method does no longer hold for all search methods presented in the following.
Some of the methods evolved closely affiliated representation concepts dedicated at a
specific problem type.
• There is a random choice made in the search direction as the algorithm iterates
toward a solution.
These properties contrast with one of the basic hypothesis of Mathematical Programming,
where it is assumed that one has perfect knowledge about the objective and depending
on the method on its derivatives and that this information is used to determine a search
direction through to a deterministic rule. In many applications, this information is not
available. This is usually referred to as black-box optimization.
A large group of stochastic search methods is inspired by biology or physics, e.g.
Simulated Annealing, Evolutionary Algorithms, Ant Colony Optimization and Swarm Op-
timization. Others imitate learning mechanisms, e.g. Tabu Search or Neural Networks.
The methods presented in the following sections rely on two model assumptions con-
cerning the objective function:
1. Pointwise sampling of the search space allows to get a kind of a problem landscape
at least locally.
These assumptions are less restrictive than the typical polynomial models employed in
Mathematical Programming. However, there are problem types which are excluded by
this model. A typical optimization problem which is not a model of the above concept is :
0 if x = 0
minimize F (x) = (7.1)
1 otherwise
neighborhood(x, r) = x0 ∈ X||x0 − x| ≤ r ⊆ X
(7.2)
In a binary search space the Hamming distance H(x, x0 ) may serve as distance measure.
It denotes the number of bits which need to be flipped to change x into x0 .
5: while continue(M, t) = 1 do
6: t←t+1
7: Create a set of alternative solutions {x̂i } ∈ X
8: Calculate F (x̂i )
9: Update the memory M ← update (M, {(x̂i , F (x̂i ))})
10: end while
11: Return the best solution x∗ in M
The steps 1, 7 and 9 may contain a random component. The function continue(M, t)
stands for some kind of termination criterion. It returns 1 as long as the current state
of the search (i.e. the iteration counter t and the solutions in the memory M ) does not
indicate a termination.
U
9: Add the new solution to the memory M ← {(x, F (xi ))}
10: end while
11: Return the best solution x∗ in M
Random search applies a pure explorative search by just randomly sampling the search
space. The creation of new solutions is not based on solutions stored in the memory.
Figure 7.1: Metropolis Algorithm: Probability p for the acceptance of a new solution in
dependence of the energy difference ∆E and the temperature T .
The probability for uphill steps increases with increasing temperature, thus exploration is
emphasized. For low temperatures the algorithm becomes similar to the stochastic descent
algorithm.
There are several alternatives introduced concerning the temperature update. A common
choice is a linear rule:
Tt+1 = update(Tt , t) = αTt (7.4)
The constant α (0 < α < 1) is called annealing constant. Simulated annealing is a global
search algorithm. The probability Pt that simulated annealing finds a global optimal
solution in t iterations converges to one:
lim Pt = 1 (7.5)
t→∞
This is a special property of the Simulated Annealing algorithm within the field of stochas-
tic search. However, it is usually not relevant for application since the value of t has to
take huge values in order to achieve a high probability Pt .
E v o l u t i o n a r y C o m p u t a t i o n
E v o lu tio n a r y G e n e tic E v o lu tio n
P r o g r a m m in g A lg o r ith m s S tr a te g ie s
E v o l u t i o n a r y A l g o r i t h m s
Evolutionary algorithms are inspired by and based upon evaluation in nature. They typi-
cally use an analogy to natural evolution to perform search by evolving solutions to prob-
lems, usually working with large collection or populations of solutions at a time. From the
point of view of classical optimization, EAs represent a powerful stochastic zeroth order
method which can find the global optimum of very rough functions.
In Structural Optimization, evolutionary algorithms are used in many different forms.
Commonly they are divided into four categories with respect to their historical background.
They are:
• Genetic Algorithms (GAs) developed by John Holland [55]. The basic terminology
of genetic search and its principal components are discussed by Goldberg [56]. An
introduction to the application of Genetic Algorithms to Structural Optimization
using traditional binary string coding is given by Hajela [57].
• Evolutionary Programming (EP) created by Lawrence Fogel [58] and developed fur-
ther by his son David Fogel [59].
• Evolution Strategies (ES) created by Ingo Rechenberg [60] and today strongly pro-
moted by Thomas Bäck [61].
• Genetic Programming (GP) is the most recent development in the field by John
Koza [62]
However, they are all based on the same evolutionary principles. Therefore we will use
a more modern terminology also used by Bentley [63] and Schoenauer [64] and generally
speak just of Evolutionary Algorithms (EAs). All the above listed strategies can be seen
as specializations of general EAs which we will describe below.
• Reproduction
• Inheritance
• Variation
• Selection
They perform the reproduction of individuals, either directly cloning parents or by using
recombination and mutation operators to allow inheritance with variation. These opera-
tors may perform many different tasks, from a simple random mutation to a complete local
search algorithm. All EAs also use some form of selection to determine which solutions will
have the opportunity to reproduce, and which will not. The key thing to remember about
selection is that it exerts selection pressure, or evolutionary pressure to guide the evolu-
tionary process towards specific areas of the search space. To do this, certain individuals
must be allocated a greater probability of having offspring compared to other individuals.
As will be shown below, selection does not only mean parent selection - it can also be
performed using fertility, replacement, or even death operators. It is also quite common
for multiple evolutionary pressures to be exerted towards more than one objective in a
single EA.
But unlike natural evolution, EAs further require three other important features:
• Initialization
• Evaluation
• Termination
Because we are not prepared to wait for the computer to evolve for several million of
generations, EAs are typically given a head start by initializing (or seeding) them with
solutions that have fixed structures and meanings, but random values. This means we
feed a certain amount of knowledge already at the beginning of the evolution into the
algorithm. Evaluation in EAs is responsible for guiding evolution towards better solutions.
Unlike natural evolution, evolutionary algorithms do not have a real environment in which
the survivability or goodness of its solutions can be tested, they must instead rely on
simulation, analysis and calculation to evaluate solutions.
Extinction is the only guaranteed way to terminate natural evolution. This is obviously
a highly unsuitable way to halt EAs, for all the evolved solutions will be lost. Instead,
explicit termination criteria are employed to halt evolution, typically when a good solution
has evolved or when a predefined number of generations have passed. In practice even more
often it happens that just the programmers patience is exceeded and he manually stops
the evolution.
There is one more important processes, which, although not necessary to trigger or
control evolution, will improve the capabilities of evolution: Mapping. Even though not
always necessary we typically separate between the search space (genotype space) and
the phenotype space. The search space is a space of coded solutions (genotypes) to the
problem, and the solution space is the space of actual solutions (phenotype). The genotype
can be understood as a recipe of how to build an actual solution. E.g., it encodes attributes
of a geometry model. Coded solutions must be mapped onto actual solutions, before the
fitness of each solution can be evaluated. The process of decoding the recipe and building
the actual solution is called mapping.
Figure 7.3 shows the general architecture of EAs. This architecture should be regarded
as a general framework for EAs, not an algorithm itself. Indeed, most EAs use only a subset
of the stages listed. In the following we will now briefly examine each of these possible
stages of EAs.
Initialization EAs typically seed the initial population with entirely random values.
Evolution is then used to discover which of the randomly sampled areas of the search
space contain better solutions, and then to converge upon that area. Sometimes the entire
population is constructed from random mutants of a single user-supplied solution. Often
random values are generated inside specified ranges (a form of constraint handling). It is
not uncommon for explicit constraint handling to be performed during initialization, by
deleting any solutions which do not satisfy the constraints and subsequently creating new
ones. More complex problems often demand alternative methods of initialization. Some
researches provide the EA with embryos - simplified non-random solutions which are then
used as starting points for the evolution of more complex solutions. Some algorithms
actually attempt to evolve representations or low-level building blocks first, then use the
results to initialize another EA which will evolve complex designs using these representa-
tions or building blocks. Although most algorithms do use solutions with fixed structures
(i.e. a fixed number of decision variables), some allow the evolution of the number and
organization of parameters in addition to parameter values. In other words, some evolve
structure as well as details. For such algorithms, initialization will typically involve the
seeding of solutions with both random values and random structures. paragraphMap Since
typically EAs distinguish between search space (genotype) and solution space (phenotype)
they require a mapping stage to convert genotypes into phenotypes. This process, known
by biologists as embryogeny, is highly important for the following reasons:
• Improved constraint handling. Mapping can ensure that all phenotypes always
satisfy all constraints, without reducing the effectiveness of the search process in any
way, by mapping every genotype onto a legal phenotype.
Especially in the field of Structural Optimization this mapping is a very crucial point for
the performance of the whole optimization. Therefore we will have a closer look at it in
Section 7.3.
Evaluation Every new phenotype must be evaluated to provide a level of goodness for
each solution. Often a single run of an EA will involve thousands of evaluations, which
means that almost all computation time is spent performing the evaluation process. In
Structural Optimization, evaluation is often performed by dedicated analysis software
(CAD-software, FEM-tools etc.) which can take minutes or even hours to evaluate a
single solution. Therefore often a strong emphasis exists toward reducing the number
of evaluations during evolution. Sometimes one even knows at the beginning how many
evaluations can be afforded and therefore the EA should be designed to make the best out
of it.
Evaluation involves the use of fitness functions to assign fitness scores to solutions.
These fitness functions can have single or multiple objectives, they can be unimodal or
multi-modal, continuous or discontinuous, smooth or noisy, static or continuously chang-
ing. EAs are known to be proficient at finding good solutions for all these types of fitness
functions. Nevertheless, the implementation of the fitness function often has tremendous
influence on the performance of the EA.
Mating Selection Parent solutions are always required in an EA, otherwise no child
solutions can be generated. However, their preferential selection of some parents instead
of others is not essential for the EA to work. Evolution will still occur without it, as long
as evolutionary pressure is exerted by one of the the other selection methods: fertility,
replacement and death. Nevertheless, most forms of EAs perform parent selection.
Choosing the fitter solutions to be parents of the next generation is the most common
and direct way of inducing a selective pressure towards the evolution of fitter solutions.
Typically, one of three selection methods are utilized: fitness ranking, tournament selection
or fitness proportionate selection. Fitness ranking sorts the population in order of the
fitness values and bases the probability of a solution being selected for parenthood on its
position in the ranking. Tournament selection bases the probability of a solution being
selected on how many other randomly picked individuals it can beat. Fitness proportionate
selection (or roulette wheel selection) bases the probability of a parent being selected on
the relative fitness score of each individual, e.g. a solution ten times as fit as another is ten
times more likely to be picked as parent (Goldberg [56]). This method also incorporates
fertility selection (see below).
Although normally fitter parents are selected, this does not have to be the case. It is
possible to select parents based on how many constraints they satisfy, or how well they
fulfill other criteria, as long as a fitness-based selecting pressure is introduced elsewhere in
the algorithm. In algorithms that record the age of individuals, parent selection may be
limited to individuals that are mature or individuals which are below their maximum life
spans.
parent solution in order to generate the child, some mutate the solution during the appli-
cation of the recombination operators, others use recombination to generate children and
then mutate these children. In addition, the probability of mutation varies depending on
the EA.
There are huge numbers of different mutation operators in use today. Examples include:
bit-mutation, translocation, segregation, inversion, structure mutation, permutation, edit-
ing, encapsulation, mutation directed by strategy parameters, and even mutation using
local search algorithms (see [56, 62, 65, 61]).
An important feature of both recombination and mutation is non-disruption. Although
variation between parent and child is essential, this variation should not produce excessive
changes to phenotypes. In other words, child solutions should always be near to their
parent solutions in the solution space. If this is not the case, i.e. if huge changes are per-
mitted, then the semblance of inheritance from parent to child solutions will be reduced,
and their position in the solution space will become excessively randomized. Evolution
relies on inheritance to ensure the preservation of useful characteristics of parent solution
in child solutions. When disruption is too high, evolution becomes no more than a random
search algorithm.
Environmental Selection Once offspring have been created, they must be inserted
into the population. EAs usually maintain populations of fixed sizes, hence for every new
individual that is inserted into the population, an existing individual must be deleted.
Therefore, this stage is also called replacement. The simpler EAs just delete every individ-
ual and replace them with new offspring. However, some EAs use an explicit replacement
operator to determine which solutions are replaced by the offspring. Replacement is often
fitness-based, i.e. children always replace solutions less fit than themselves, or the weakest
in the population are replaced by fitter offspring.
Replacement is clearly a third method of introducing evolutionary pressure to EAs, but
instead of being a selection method, it is a negative selection method. In other words,
instead of choosing which individuals should reproduce or how many offspring they should
have, replacement chooses which individuals will die.
Replacement needs not be fitness-based, it can be based on constraint satisfaction, the
similarity of genotypes, the age of solutions, or any other criterion, as long as a fitness-
based evolutionary pressure is exerted elsewhere in the EA. Replacement is also limited by
speciation within EAs: a child from two parents of one species/population/island should
not replace an individual in a different species/population/island.
Float-list-gene is a list of arbitrary floating point values. The parameter always has to
represent one of the values specified in the definition of the gene.
Const-float-list-gene are equally distributed floating-point values, i.e. they have a con-
stant distance between two neighboring values.
This gene is quite similar to the integer gene. Optionally, limits or a mutation step
size can be provided.
String-list-gene is a list of arbitrary discrete values upon which no norm or ordering can
be applied.
This is in contrast to the integer-gene, the float-list-gene, and the const-float-list-
gene.
For the definition of the listed gene types, the standard deviation σ used for Gaussian
mutation must be specified additionally for every gene.
Float-gene, integer-gene, float-list-gene, and const-float-list-gene can further be pro-
vided with so called cyclic properties.
This is suitable when between two possible gene values a distance but no absolute
order can be defined, as e.g. for angle values.
• [float gene 1.2 true 0.5 false 0.0 3.0 5.0 false] denotes a float gene with a default
value of 1.2, with a lower limit 0.5, and with no upper limit defined (the booleans
before the limit values specify, whether a limit exists or not).
Further the mutation parameters = 3.0 and σ = 5.0 are specified.
denotes the range for uniform mutation of an unbounded gene, and σ defines the
standard deviation to be used for Gaussian mutation.
The last boolean value indicates that this example gene does not have cyclic prop-
erties.
• [int gene 3 true 0 true 360 0 10.0 true] denotes an int-gene with a default value
of 3, with a lower limit 0, and an upper limit of 360.
Further, mutation parameters = 0 (useless for this fully bounded gene), and σ =
10.0 are specified.
In addition, this gene has cyclic properties, indicated by the last boolean value.
• [bool gene true 2.0] denotes a bool-gene with default value true.
A standard deviation σ = 2.0 for some sort of Gaussian mutation is also defined.
• [float list gene 4 3 0.0 1.0 false 1.2 2.5 4.5 4.9] denotes a float-list-gene with 4
values, where the value with index 3 (starting from 0) is the actual default value.
Mutation parameters = 0.0 (useless for this gene) and σ = 1.0 follow.
Then it is defined that this gene has no cyclic properties, and the four possible float
values are appended.
• [const float list gene 1.75 1.0 true 3.0 false 8 0.5 1.2 false] denotes a const-
float-list-gene with a default value of 1.75.
The gene has an active lower limit of 1 and an inactive (indicated through the boolean
after the limit value) upper limit of 3.
By specifying the number of intervals to be 8, the possible values for this genes are
[1.0, 1.25, 1.5, ...].
Finally, = 0.5, σ = 1.2 and no cyclic properties are given.
• [string gene 3 1 1.0 blue green red ] denotes a string-gene with 3 possible values
and where the value with index 1 (again starting from 0) is chosen to be the actual
default value. Further a σ = 1.0 value is specified, and the different possible values
are appended.
During the initialization of the Evolutionary Algorithm, a list containing such genes is
read from file.
This list defines the genotype structure for a problem at hand.
Further, the given default values can also be used to incorporate existing solutions for
knowledge-based initialization.
Composite Structures
This chapter discusses the structural optimization of laminated composite materials with
focus on fiber reinforced materials. The thickness of these structures is usually thin com-
pared to the other dimensions which is why they are typically represented in FEM simula-
tions with layered shell finite elements. The first section gives a short introduction on the
design of laminated composites including a selection of applications. Consequently, the
Classical Lamination Theory (CLT) is introduced which is the standard calculus for lam-
inated composites. Afterward, the basics of a shell finite element are explained. Together
with the CLT, the element is enhanced to a layered shell element which can be used for
the computational representation of laminates. The majority of these fundamentals are
based on the doctoral thesis of B. Schläpfer [66]. The second part of this chapter discusses
different disciplines in laminate optimization including fiber orientation, laminate thick-
ness, stacking sequence and material optimization. Finally, the theories are specialized to
optimization techniques for finding locally varying laminates which is also called laminate
tailoring. Additionally, a selection of laminate tailoring methods is presented.
sented within the next sections. A selection of applications with shell characteristics are
(a) Fuselage section of the Air- (b) Payload fairing of the (c) CFRP-monocoque of a
bus A350 Ariane 5 launcher Lamborghini Aventador
presented in Figure 8.1. Figure 8.1(a) shows a fuselage section of the Airbus 350 which
is entirely made of Carbon Fiber Reinforced Plastics (CFRP). Considering light-weight
structures or especially aircraft structures, the thickness of almost any components is
usually much smaller compared to the other dimensions wherefore they are prevailingly
modeled with shells. Another typical shell application is the payload fairing of the Ariane 5
launcher shown in Figure 8.1(b). Its mission is to protect the payload during the launch
and it should of course be as light as possible in order to maximize the payload capacity of
the launcher. Additionally, it must fulfill acoustic and dynamic requirements. The third
example shown in Figure 8.1(c) is a CFRP-monocoque of a Lamborghini. With increasing
fuel costs, weight saving has become an important issue in automotive engineering and
components are increasingly often made of fiber reinforced plastics.
are highly anisotropic. While the mass-specific properties in fiber direction may be one
magnitude higher than for metallic materials, the material properties perpendicular to the
fiber orientation, which are mainly governed by the properties of the matrix material, are
comparatively low. Therefore, structures are built by stacking several plies with different
orientation angles which are then called laminated composites (see Figure 8.2(b)). Basi-
cally, the orientation angles can be chosen in a way so that the overall mechanical behavior
becomes nearly isotropic (also called quasi-isotropic). This is achieved by distributing the
fiber orientations regularly in all directions. A big advantage of laminated composites
is the possibility to specifically design their structural response to different loadings and
requirements. In contrast to isotropic materials for which the amount of material, in
particular the thickness, is the only design parameter, the behavior of composites can be
designed by varying the orientation angles, the number of plies, the stacking sequence or
in special cases the fiber volume content. Material can be added easily to highly loaded
regions and the fibers can be oriented dependent on the principal directions of the loads.
Furthermore, the layered design enables to build laminates of different ply materials. Due
to the behavior of anisotropic materials, the layered design method and the resulting large
number of design variables, the design process of laminated composites becomes complex
and time consuming. The experience and the intuition of structural engineers may be
pushed to the limit and, excepting structures with a low degree of complexity, e.g. plates
or cylinders, finding solutions with good light-weight properties is hard to be accomplished
manually. The application of computer-based design methods is indispensable aiming for
high quality designs in a reasonable time frame.
• the structure be thin compared to the other dimensions and a constant thickness.
• the Bernoulli-theory be valid, namely plane sections and no transverse shear strains.
• deformations be small.
Supplementary to the assumptions for homogeneous materials, it is assumed that the single
layers are bonded perfectly. Due to the stacking of several layers, the material stiffness
properties of a laminate through the thickness are inhomogeneous. In order to provide
a linear relation between plate deformations to plate line loads of the global laminate,
the CLT performs a stiffness homogenization. This requires that the layer stiffnesses,
which are usually given in material principal coordinates 1,2, be transformed to the global
coordinates x, y of the problem.
In order to express the stiffness entries as functions of the engineering constants such as
the Young’s moduli, the Poisson ratios and the shear moduli, it is easier to formulate
the so-called compliance matrix S which is the inverse of the stiffness matrix C. The
compliance matrix of an orthotropic material is defined as
1
− Eν21 − Eν31
E11 22 33
0 0 0
− ν12 1
− Eν32 0 0 0
E11 E22 33
− ν13 − Eν23 1
0 0 0
E11 E33
Sortho = 22 (8.2)
1
0 0 0 0 0
G23
0 1
0 0 0 G31 0
1
0 0 0 0 0 G12
Measured engineering constants can be inserted directly and the stiffness matrix results
of the inversion. Considering a unidirectional reinforced composite material, the material
properties are transversely-isotropic. Such a material is defined with 5 elastic coefficients
which are usually the Young’s modulus in fiber direction E11 , the Young’s modulus per-
pendicular to the fibers E22 , the in-plane shear modulus G12 , the in-plane Poissons ratio
ν12 and the transverse shear modulus G23 . The transversal shear Poisson ratio is related
to the variables above with
E22
ν23 = −1 (8.3)
2G23
Due to the assumption of a plane stress state, the stresses with out-of-plane components
are zero and the compliance matrix S can be reduced to a 3×3-matrix. Consequently, the
stiffness matrix C is reduced to the so-called reduced stiffness matrix Q and the Hooke’s
law is simplified to
σ1 Q11 Q12 0 ε1
σ2 = Q12 Q22 0 ε2 (8.4)
τ12 0 0 Q66 γ12
Keep in mind that the strains having out-of-plane components, namely ε3 , γ13 and γ23 ,
are not zero, they are just not considered. Strains and stresses are expressed in material
principal coordinates 1,2. However, the orientation angle of a unidirectional laminate
layer can be chosen arbitrarily wherefore the material principal directions do not coincide
with the global coordinates of the problem formulation. Consider a laminate illustrated
in Figure 8.3 where the material coordinates 1,2 are rotated with respect to the global
coordinates x,y by a given angle ϕ. A relation which transforms the stresses and the
2 y
φ
x
strains from the material coordinate system to the global coordinate system is needed.
In order to be able to work directly with the engineering strains without a pre-factor,
Reuter [68] introduced the simple matrix R.
1 0 0
R = 0 1 0 (8.9)
0 0 2
An application of this Reuter matrix reduces the potential for mistakes to a minimum
since it is clear that only engineering stains are utilized. The engineering strains are then
transformed with equation (8.10).
Contrariwise, the global strains and stresses can be mapped to the material principal
strains and stresses performing a multiplication with the inverse rotation matrix T−1 .
The connection between the strains and stresses in global coordinates is based on the
transformed reduced stiffness matrix Q.
σx εx
σy = Q εy (8.11)
τxy γxy
The transformed reduced stiffness matrix Q can be derived by using equations (8.4), (8.6)
and (8.10) and is therefore a function of Q, R and T only.
The transformed reduced stiffness matrix Q has to be evaluated for every laminate layer
in order to be able to homogenize the material data of the entire laminate stack. However,
the evaluation is straight forward and the computational costs are low.
Stiffness Homogenization
In contrast to a homogeneous material, the material properties of a laminate are alter-
nating through the thickness. This complicates the analysis but also the formulation of
a finite element, since the integration through the thickness becomes more complex. The
stiffness discontinuity due to the layer-wise design technique results in a discontinuous
stress distribution. To quantify a stress state in a laminate, the components have to be
evaluated for each layer particularly. In order to have a load unit that includes all layers,
the line loads are introduced, namely the force per unit length N and the moment per
unit length M. The force per unit length N is the integration of the stress components
over the laminate thickness t.
t
zk
zk-1 hk
hk-1 t
zj 2
hj
z1 h2
z0 h1
Z t n Z zj n Z zj
2 X X
ε0 + zκ zdz
M= σzdz = Qj εzdz = Qj (8.18)
− 2t j=1 zj−1 j=1 zj−1
n
X 1 2 2
0 1 3 3
= Qj zj − zj−1 ε + zj − zj−1 κ (8.19)
2 3
j=1
Considering equations (8.17) and (8.19), the matrices A, B and D can be extracted which
are all of dimension 3×3. The A-matrix connects the membrane strains ε0 with the force
per unit length N.
Xn
A= Qj (zj − zj−1 ) (8.20)
j=1
and the D-matrix connects the plate curvatures κ with the moments per unit length M.
n
1X
Qj zj3 − zj−1
3
D= (8.21)
3
j=1
The B-matrix is responsible for the coupling of membrane and bending components.
n
1X
Qj zj2 − zj−1
2
B= (8.22)
2
j=1
Dependent on the laminate layup, some entries may become zero. In case of a symmet-
ric laminate, there is no coupling between bending and membrane effects wherefore the
B-matrix vanishes. These matrices build the so-called ABD-matrix which is the main
achievement of the homogenization process.
A B ε0
N
= (8.23)
M B D κ
ε0x
Nx A11 A12 A16 B11 B12 B16
Ny A12 A22 A26 B12 B22 B26 ε0y
ε0xy
Nxy
= A16 A26 A66 B16 B26 B66
(8.24)
Mx
B11 B12 B16 D11 D12 D16
κx
My B12 B22 B26 D12 D22 D26 κy
Mxy B16 B26 B66 D16 D26 D66 κxy
The application of shell elements requires that the thickness of the represented structures
is much smaller compared to its other dimensions. Thin structures are usually not modeled
with 3D elements, especially if analyzing plate bending. If 3D elements are modeled thin
in only the thickness direction according to Figure 8.5(a), there may be the problem
of shear locking and ill-conditioning [69]. This problem can be avoided by using many
elements according to Figure 8.5(b). However, the number of degrees-of-freedom (d.o.f.) is
increased drastically and the models become large wherefore solving the problems becomes
computationally expensive. Using shell elements and the underlying plate theory, the
problem of having a large number of d.o.f. is mitigated (see Figure 8.5(c)).
(a) 3D-modeling with large aspect ratios may lead to shear lock-
ing and ill-conditioning
Stiffness Matrix
k = km + k c + kb + ks (8.25)
with
n Z
X
km = BTm Qj Bm dA (zj − zj−1 ) (8.26)
j=1 A
n Z
X 1 2
BTm Qj Bb BTb Qj Bm 2
kc = + dA zj − zj−1 (8.27)
A 2
j=1
n Z
X 1 3
BTb Qj Bb dA 3
kb = z − zj−1 (8.28)
A 3 j
j=1
n Z
X S
ks = BTs κQj Bs dA (zj − zj−1 ) (8.29)
j=1 A
Consider that the matrices Bm , Bb and Bs are strain-displacement matrices (see Appendix
A.2) and should not be confused with the B-matrix of the lamination theory. In contrast to
homogeneous materials, an additional term from the coupling stiffness matrix kc appears
which becomes zero again for symmetric laminates. However, no additional components
are need since it is a combination of membrane and bending parts.
Mass Matrix
The mass matrix formulation for a layered shell element is simple when using a lumped
mass model where rotary inertia is neglected. It is feasible for small deformations which is
usually true for harmonic vibration problems. The lumped mass matrix is expressed with
equation
Xn Z
T
m= ρj tj φ φ dA (8.30)
j=1 A
df f (x + ∆x) − f (x)
∇f = = lim (8.31)
dx ∆x→0 ∆x
A simple method for the numerical determination of the sensitivities is the finite difference
approximation. The forward finite difference approximation is defined as
of the sensitivities may become too expensive. Additionally, the numerical determination
includes a numerical error which is caused by the finite value of ∆xi . However, a numerical
determination of the sensitivities is the only possibility if the objective function is not
known explicitly, for example in black box optimizations.
Considering for instance the common sensitivity equations for minimal compliance
∂W
= −uT ∂K T ∂r W = uT r
∂x ∂x u + 2u ∂x (8.33)
for displacements which has already been derived in Section 6.10.2 (see equation (6.133))
∂u
= K−1
∂r ∂K
∂x ∂x − ∂x {u = Kr} (8.34)
or eigenfrequencies
∂K ∂M
∂λn ΦT
n ( ∂x −λn ∂x )Φn
∂x = T
Φn MΦn
{(K − λn M) Φn = 0}
(8.35)
it can be noticed that they both require the sensitivities of the stiffness matrices. Since the
formulation of the finite elements, and therefore the connections between design variables
and the stiffness matrix are known, the sensitivities can be derived analytically which is
done here for the thickness of the layers and the orientation angles.
j
X
zj = hk (8.36)
k=1
The finite element stiffness matrices for membrane, bending and transverse shear conse-
quently yield
n/2 Z j j−1
" !#
X X X
km = 2 BTm Qj Bm dA hk − hk (8.37)
j=1 A k=1 k=1
n/2 Z
j
!3 j−1
!3
X 1
BTb Qj Bb dA
X X
kb = 2 hk − hk (8.38)
A 3
j=1 k=1 k=1
n/2 j j−1
"Z !#
X S X X
ks = 2 BTs κQj Bs dA hk − hk (8.39)
j=1 A k=1 k=1
Due to the symmetry of the laminate, the coupling stiffness matrix kc is zero. Taking
advantage of relation
Xj j−1
X
hk − hk = hj (8.40)
k=1 k=1
the stiffness for the membrane and shear parts can be simplified to
n/2 Z
X
km = 2 BTm Qj Bm dA hj (8.41)
j=1 A
n/2 Z
X S
ks = 2 BTs κQj Bs dA hj (8.42)
j=1 A
Consequently, the sensitivities of these parts are only dependent on the material properties
of the respective layer and yield
Z
dkm
=2 BTm Ql Bm dA (8.43)
dhl A
Z
dks S
=2 BTs κQl Bs dA (8.44)
dhl A
For simplicity, it is assumed that the stiffness correction factor κ is independent on the
thickness hl . The derivative of bending part is more complex because the thickness appears
with a higher order. Considering equation (8.38), it must be distinguished whether the
layer with respect to which the derivative is taken is part of the summation or not. This
can be expressed with
j
!
∂ X 1 for j ≥ l
hk = (8.45)
∂hl 0 for j < l
k=1
Z l
!2
dkb X
=2 BTb Ql Bb dA hk
dhl A k=1
n/2
j
!2 j−1
!2
X Z X X
+2 BTb Qj Bb dA hk − hk (8.46)
j=l+1 A k=1 k=1
Here it becomes obvious that a change of the bending stiffness is caused by two different
effects. The first term expresses the stiffness change due to the thickness change of the
layer itself. The terms within the summation consider the change of the stiffness caused by
pushing outward the overlaying layers which results in a higher area moment of inertia (see
Figure 8.8). This relation can be written alternatively with expressions (8.47) and (8.48).
The top layer (l = n/2) is independent on the subjacent layers wherefore its sensitivity is
expressed with
2
Z n/2
dkb X
=2 BTb Ql Bb dA hk (8.47)
dhl A k=1
The lower layers (l = 1, ..., n/2 − 1) take into account the stiffness contribution of the
overlaying layers wherefore
Z Z X l
!2
dkb dkb
=2 BTb Ql Bb dA − BTb Ql+1 Bb dA hk + (8.48)
dhl A A dhj+1
k=1
The total sensitivities of the stiffness matrix are built with an addressed summation of the
several parts corresponding to equation (8.25).
dk dkm dkb dks
= + + (8.49)
dhl dhl dhl dhl
All the derivations above have been performed on the element level for the element stiffness
k. However, the sensitivity equations (8.33) through (8.35) require the sensitivities of the
global stiffness matrix K. They can however be derived directly from the element stiffness
parts. The global stiffness matrix is an addressed summation of all the element matrices
which is denoted schematically in equation (8.50).
k1 0
k2
K= (8.50)
k3
..
0 .
The element stiffness matrix derivatives are dependent only on the thickness of the respec-
tive element. The derivatives with respect to the thicknesses of other elements are all zero.
dK
Thus, the global stiffness matrix derivative dh l
contains only zeros except of the entries
corresponding to the considered element l which is shown schematically in equation (8.51).
0 0
dK dkl
= (8.51)
dhl hl
0 0
The sensitivity of the mass matrix with respect to a layer thickness change is simply
Z
dm
= ρl φT φ dA (8.52)
dhl A
z,w
A
y,v
t
x,u
zn/2 t
zj 2
z1
∆t z3
k3(z2,z3) 3
z z2 k3(z2+∆h,z3+∆h)
z
k2(z1,z2) 2 z1 k2(z1+∆h,z2+∆h)
z k1(z0+∆h,z1)
k1(z0,z1) z1 z0
0
∂Qj
Z
∂km
= BTm Bm dA (zj − zj−1 ) (8.53)
∂ϕj A ∂ϕj
since only the transformed reduced material stiffness matrix Qj is dependent on the fiber
orientations This is the same for all other stiffness parts given in equations (8.27) through
(8.29). The derivative if the reduced material stiffness matrix starts from equation (8.12)
and yields
∂Q ∂T−1 ∂T −1
= QRTR−1 + T−1 QR R (8.54)
∂ϕj ∂ϕj ∂ϕj
The derivative of the transformation matrix introduced in equation (8.5) is
2 cos2 ϕ − 2 sin2 ϕ
−2 sin ϕ cos ϕ 2 sin ϕ cos ϕ
∂T
= 2 sin ϕ cos ϕ −2 sin ϕ cos ϕ −2 cos2 ϕ + 2 sin2 ϕ (8.55)
∂ϕj
− cos2 ϕ + sin2 ϕ cos2 ϕ − sin2 ϕ −4 cos ϕ sin ϕ
The mass matrix is not sensitive to a change of the fiber orientations wherefore their
sensitivities become zero.
with
Φ = {cos 2ϕ, sin 2ϕ, cos 4ϕ, sin 4ϕ} (8.63)
Assuming that the thickness of all layers is unique and taking into account the fact that
the orientation angles ϕ are constant through the thickness of one layer, the integrals can
be replaced with summations according to equations (8.64), (8.65) and (8.66).
n
A A A A 1X
V1 , V2 , V3 , V4 = [Φi ] (8.64)
n
i=1
n 2
B B B B 2 X n 2 n
V1 , V2 , V3 , V4 = 2 −i+1 − −i Φi (8.65)
n 2 2
i=1
n 3
D D D D 4 X n 3 n
V1 , V2 , V3 , V4 = 3 −i+1 − −i Φi (8.66)
n 2 2
i=1
with
Φi = {cos 2ϕi , sin 2ϕi , cos 4ϕi , sin 4ϕi } (8.67)
Laminate optimization by means of the lamination parameter theory operates directly on
the lamination parameters ViA,B,D and not on the physical design variables such as for
example the orientation angle or layer thickness. Consequently, the number of linearly in-
dependent parameters is 12, for laminates which are symmetric and balanced even only 4.
Naturally, the lamination parameters cannot be chosen arbitrarily since they are bounded
to a feasible domain through the underlying trigonometric functions. A few publications
[75, 76, 77] focus on the description of the feasible lamination parameter domains, but so
far, there exists no generally valid analytical approach. However, if assuming a balanced
and symmetric laminate, the feasible domain of the in-plane lamination parameters can
be found rather easily. The contours of the feasible domain for the in-plane lamination
parameters can be determined, by evaluating equation (8.64) for the border case. An
illustration of the feasible domain for the lamination parameters V1A and V3A is given
in Figure 8.9. The lower boundary is defined by (±ϕ)S -laminates. The upper straight
VA3
(0/903)S 1 (03/90)S
(904)S (04)S
(02/902)S
(0/45/902)S (02/45/90)S
(±752)S (±152)S
(45/903)S (03/45)S
(0/452/90)S V A1
−1 (452/902)S (02/452)S 1
(±302)S
(±602)S (453/90)S (0/453)S
−1 (±452)S
Figure 8.9: Feasible domain of the in-plane lamination parameters V1A, , V3A,
boundary is given by (0j /90n−j )S -laminates. For symmetric and balanced laminates, the
in-plane lamination parameters V2A and V4A are always zero due to the integration over
the thickness. Additionally, all coupling lamination parameters ViB are equal to zero for
the symmetric case. For simplicity, only the +45◦ -layers are labelled in Figure 8.9. Since
the lamination parameters V1A and V3A are based on cosine-functions, each +45◦ -layer can
directly be replaced with a −45◦ -layer
It can be shown that the feasible domain is convex wherefore an optimal solution of the
ABD-matrix can be found efficiently by using algorithms of mathematical programming.
However, the major drawback of using lamination parameters is the fact that the infor-
mation of the ABD-matrix does not explicitly include the information of the physical
laminate. The back-substitution from the ABD-matrix or from the lamination parame-
ters, respectively, to the physical laminate with a known stacking sequence and orientation
angles is not unique and includes a second optimization problem. Considering again the
feasible domain in Figure 8.9, the lamination parameters can take every value in the do-
main. However, the feasible physical laminates are distributed discretely in the design
domain. Having a given number of layers, which is a basic assumption of the theory, the
optimal lamination parameters can only be represented approximately by a physical lam-
inate. In case of a symmetric and balanced laminate as shown above, the complexity of
the back-substitution is acceptable. The problem becomes much more complex for general
laminates.
The same calculations can be done for a laminate with maximal in-plane shear stiffness:
Alternatively, we may seek for a laminate which has equal stiffness in all directions, which
is also known as quasi-isotropic laminate. The equations characterizing a quasi-isotropic
laminate are
A11 = A22 → U1 + V1A · U2 + V3A · U3 = U1 − V1A · U2 + V3A · U3
1
A66 = 2 (A11 − A12 ) → U5 − V3A · U3 = 12 U1 + V1A · U2 + V3A · U3 − U4 + V3A · U3
V3A = 0
Since we assume the lamination parameters V2A and V4A to be zero, the laminate must
be balanced wherefore the solution corresponds to (0/45/-45/90)S -laminate. However,
this is not the only solution for a quasi-isotropic laminate. Also a laminate with layup
(0/60/-60)S is quasi-isotropic, symmetric and balanced and results in the same lamina-
tion parameters (V1A = V3A = 0). This demonstrates that the back-substitution for the
lamination parameters to the physical laminate is not is not unique. Of course, the two
laminates (0/45/-45/90)S and (0/60/-60)S produce different entries ViD . However, the
difference becomes small if having a large number of layers.
x0
Optimizer
x Optimization
Environment
FEM FEM
FEM FEM
Output FEM-Solver Input
Model Environment
File File
is predefined and fixed during the optimization. The orientation angles are taken as
design variables wherefore the optimization problem becomes continuous. Figure 8.11(a)
schematically illustrates a laminate consisting of four layers whose orientation angles are
optimized. The design variables can simply be arranged in an array which is illustrated
φ1 = x1,1 φ1 = x1,n
φ2 = x2,1 φ2 = x2,n
φ3 = x3,1 φ3 = x3,n
φ4 = x4,1 φ4 = x4,n
(a) Fiber orientation optimization scheme
in Figure 8.11(b). In general, the fiber orientation values are real numbers (ϕj ∈ R).
Due to the periodicity of the orientation angles, their are often bounded to an interval of
180◦ , such as [0◦ , 180◦ ] or [−90◦ , 90◦ ]. Consequently, the search space is restricted without
reducing the potential solution quality. However, some optimization algorithms cannot
handle constrained problems. Of course, the search space is not restricted using this type
of algorithms. The continuity of the design variables allows to solve fiber orientation
optimization with algorithms of the class of mathematical programming. However, it
must be considered that the search space is multi-modal which may be problematic for
this class of algorithms since they risk to get stuck in a local optima. The non-convexity
of a search space with two orientation angles has been visualized by Keller [78]. Figure
0.17 0.2
60 0.18 0.18 60
0.21
0.2 0.18 0.21 0.19
45 0.18 45
30 30
15 15
θ2
y -15 -15
-30 -30
0.18
-45 -45
0.21 0.18 0.21 0.19
θ
-60 0.18 -60
x 0. 0.17 .2
-75 0.19 0.16 -75
0.18 0.18 0.18
-90 0.18 0.18 -90
-90 -75 -60 -45 -30 -15 0 15 30 45 60 75
θ1
(a) The structure is clamped at the (b) Contour plot of the objective for the tensile
left side and a uniform line-load of specimen experiment as a function of the ply an-
0.01 N/mm is applied to right side. gles φ1 and φ2 (in degrees).
Figure 8.12: Geometry and objective of the tensile specimen experiment. Source: Keller
[78]
8.12(a) illustrates a simple plate which is clamped at the left side and loaded with a
uniform line load. The plate consists of two layers with orientation angles x1 and x2 .
Figure 8.12(b) shows the corresponding contour plot of the objective function which is
obviously multi-modal. There will even be more local optima when using more than two
layers.
Allowing the orientation angles to be real may lead to solutions with many positions
after decimal points which is inappropriate for manufacturing since processes are often re-
stricted to discrete steps, for instance to 5◦ -steps or higher. This renders the problem to be
discrete wherefore the application of algorithms of the category of mathematical program-
ming becomes impossible. Stochastic algorithms may mitigate that problem. Alterna-
tively, continuous orientation angles may be rounded to the next feasible angle. However,
if the discrete steps are too large, the solution may change its behavior significantly.
Another optimization discipline, which has already been performed in the early stages of
computer-aided laminate design, is the optimization of the laminate stacking sequence. It
basically modifies the order of the laminate stack by exchanging the layers or its position
in the laminate, respectively. Stacking sequence optimization is important if the load cases
include bending. For such a case, it is absolutely essential which layers are located in the
outer positions since these layers have a greater impact on the area moment of inertia.
Moreover, the stacking sequence is important for the mechanical coupling between bend-
ing and twisting which is mathematically represented by the B-matrix. If the layers are
distributed appropriately, an in-plane load may induce out-of-plane deformations and vice
versa which may be interesting for specific applications. Additionally, the interlaminar
stress distribution is strongly dependent on the distribution of the single layers. Figure
8.13(a) schematically illustrates a stacking sequence optimization. The parametrization of
A B
B D
C A
D C
(a) Stacking sequence optimization scheme
A B C D B D A C
(b) Stacking sequence optimization parametrization
a stacking sequence is usually solved by arranging the numbered layers in an array accord-
ing to Figure 8.13(b). An exchange of the array sequence is equivalent to an exchange of
the corresponding layers. When exchanging two layers in the laminate stack, the objec-
tive function changes by leaps and bounds wherefore the optimization problem becomes
discrete. Consequently, these category of problems can only be solved with stochastic
algorithms. A type of algorithms which has been used very often to optimize the stacking
sequence are Genetic Algorithms (GA) [79, 80, 81, 82, 83, 84, 85, 86, 87]. Basically, a
stacking sequence problem can be transformed in to a problem of fiber orientation op-
timization. Instead of exchanging the layers, their fiber orientations may be varied and
potentially, the same solution can be found.
The weighting factors wi , which are bounded by 0 and 1, are used as design variables where-
fore the problem becomes continuous. According to the topology optimization method of
Bendsøe and Kikuchi (see Section 9.4), the design variables are pushed towards their
boundaries and the final weighting factor array of each layer must have one value of 1
and the rest must be zero. If defining several constitutive matrices made of one material
but different orientation angles, the method can be used for fiber orientation or stacking
sequence optimizations.
t1 = x1,1 t1 = x1,n
t2 = x2,1 t2 = x2,n
t3 = x3,1 t3 = x3,n
t4 = x4,1 t4 = x4,n
φ1 = x1,n
φ1 = x1,1
φ2 = x2,n
φ2 = x2,1
φ3 = x3,n
φ3 = x3,1
φ4 = x4,n
φ4 = x4,1
φ5 = x5,n
(a) Scheme of a problem with variable number of layers
optimization problem with a variable number of layers. For each new appearing layer, the
fiber orientation, the thickness, the material and the position in the laminate stack has to
be determined (even if in Figure 8.15 only the orientation angles are indicated), which is
usually done randomly or by choosing them randomly from a predefined database, respec-
tively. Due to the variable length of the design variable array, these kinds of problems can
not be solved with algorithms from Mathematical Programming. Usually, Genetic Algo-
rithms with variable genotype length are applied. It is important to implement appropriate
reproduction and mutation routines in order to come close to optimal solutions.
Homogeneous
optimal laminate which can be load
found rather easily.
distributions occurHowever, in general,
only in academic the internal loads
examples
are inhomogeneous and the stress state at each point of the laminate is different. Figure
Motivation
8.17 shows the strain field of a plate with a centered hole which is loaded uniaxially. Even
18. December 2012 Doctoral Examination – Benjamin Schläpfer 7
Find global
Figure laminate layup strain
8.17: Inhomogeneous which field
is theofbest solution
a plate with for the hole
centered
structural behavior → local loadings are not respected
Laminatoptimierung
if this problem (Laminate-Tailoring)
is quite simple, the internal load distribution is highly inhomogeneous.
Tailor
In order to the local
fully exploit thelaminate properties
potential to the the
of composites, locallaminate
needs by taking have to be
properties
advantage of locally varying thickness or fiber orientations
tailored to the local load conditions. This can be achieved by splitting the design domain
locally
Ω into sub-domains Ω varying
Die grundlegenden anisotropy
each havingParametrisierungen
a unique layup, und Optimierungsmethoden
which is schematically illustrated in
bleibeni dieselben
Figure 8.18. Consequently, the objective here is to find an optimal
Die Design-Domain wird in sogenannte Sub-Domains aufgeteilt laminate for each sub-
ȍ3
ȍ1 ȍ2
a relatively new discipline since the computer system requirements are rather high. The
number of design variables increases linearly with the number of sub-domains which makes
the optimization more computational expensive.
Performing laminate tailoring, there might be the risk of finding solutions which can-
not be physically realized since they are not manufacturing friendly. Consequently, man-
ufacturing aspects should be regarded in the parametrization concept. For instance, the
Elemente
Sub-domain
Region 11 Sub-domain
Region 2 2 Sub-domain
Region 33
sub-domains. Layers A, C and E are covering the entire domain and a global connectiv-
ity is given. A simple method to guarantee connectivity is to define a number of plies
which cover the entire domain and may not be adapted during the optimization. This
may restrict again the search space.
Basically, the geometric representation of the sub-domains can be predefined and as-
sumed to be fixed during the optimization. This results in a constant number of design
variables which simplifies the optimization, but also restricts the search space and the so-
lution quality, respectively. More freedom in terms of design is achieved if the parametriza-
tion additionally includes the geometrical representation of the sub-domains. A variable
number of sub-domains leads to a variable number of design variables which can only
be solved by means of stochastic algorithms. Two parametrization approaches for the
sub-domains are presented in the next two sections.
Ω1 Ω4 Ω1
Ω2 Ω…
Ω3 Ωn Ω2 Ωn
(a) Parametrization with single elements (b) Parametrization with multiple ele-
representing a sub-domain ments pooled into a sub-domain
dependent on the finite element mesh which may restrict the design freedom, but only if a
rough mesh is chosen. However, no additional parametrization of the geometry is needed.
Alternatively, the number of design variables can be reduced by pooling a given number
of finite elements into groups which then represent the design domains. A schematic
sketch is shown in Figure 8.20(b). One possibility is to pre-define the elements belonging
to one sub-domain and keep it fixed during the optimization. This of course restricts
the freedom of design but it might be helpful to make the solutions more manufacturing
friendly. Alternatively, the elements which belonging to a sub-domain may be chosen to
be variable. The realization and the optimization become more complex but the potential
space for improved solutions becomes bigger. Consider that for both parameterizations,
the solution quality, but also the computational costs are dependent on the finite element
mesh.
Giger et al.[94] proposed a graph-based parametrization scheme which allows to make
the elements which are pooled in a sub-domain to be variable. This parametrization
scheme allows maximal freedom in terms of design since regions having a unique laminate
are morphing. The material data and the orientation angle of a so called patch, which
can be understood as a layer that partially covers the design domain, is stored in a patch
vertex as shown in Figure 8.21. The elements belonging to a patch, but also the material
and the orientation angle can be modified by means of a genetic algorithm. In order to
guarantee the connectivity of a patch, only boundary elements may be added or removed
which are located automatically by a computer routine. The total laminate layup in each
finite element is set up by summarizing the patch layers the respective element belongs
to in the right order (according to Figure 8.21 for instance from left to right). This
parametrization scheme is inspired by the manufacturing process and the solutions are
manufacturing-friendly.
An alternative approach for laminate tailoring based on a finite element based
parametrization is proposed by Schläpfer and Kress [93, 95, 66]. The method basically
locally reinforces a predefined laminate structure with laminate patches regarding a spe-
cific design criterion such as stiffness, strength or dynamic behavior by simultaneously
aiming for minimal mass. Layers with given material and fiber orientation are predefined
in the design domain. Some of these layers have a thickness of zero wherefore they do not
really contribute to the finite element model. These zero-thickness layers act as potential
reinforcements. The local reinforcements are generated by partially setting the thickness
of these layers to values of the semi-finished material. The decision, whether a layer in
one finite element is set to a real value or not is made on the sensitivities which have been
0° 0°
45° 45°
-45° -45°
90° 90°
low high
(a) Red areas indicate high sensitivities which (b) The resulting local reinforcements are rep-
implies that the thickness of the corresponding resented by the black areas
layer has to be increased in order to improve the
design value
Figure 8.22: Sensitivities and local reinforcements of a vibration plate for which the second
and the third eigenfrequency are separated. Source: Schläpfer [66]
Figure 8.23: Connection between global layers and laminate regions. Source: Zehnder [96]
rather computational expensive. Also this method works with a genetic algorithm.
r e g io n to b e o p tim iz e d
c o n s ta n t- th ic k n e s s s u r fa c e la y e r o r c a m b iu m
Figure 9.1: Part with stress raising geometry and modelled cambium
two phases. The first phase determines the stress distribution within the cambium, due to
the specified geometric boundary conditions and applied forces, by a structural analysis.
The stresses are reduced to some equivalent stress distribution. During the second phase,
a temperature distribution is derived from the equivalent stress distribution. A tempera-
ture expansion coefficient is appropriately chosen and Young’s modulus of the cambium is
reduced by some factor, for instance 400. The structural analysis of the second phase then
simulates a growth of the cambium effected by the temperature strains. The strains tend
to be higher where higher equivalent stress values were calculated by the preceding first
phase. The strains imply a growth of the structural model. The growth is numerically
recorded by adding the nodal-point displacements to the nodal point reference positions.
The so obtained geometry changes can have the effect that the stress concentrations due
to the mechanical loading are reduced. The process is illustrated in Fig. 9.2 During the
iterations the geometry changes automatically so that the stress concentrations are sys-
tematically reduced.
The advantage of the method with its ingenious simplicity lies in the relatively low
numerical effort to obtain a significantly improved design. The method was used with
much success to reduce stress concentrations, and dramatically increase lifetime, in load-
carrying automotive parts that could be re-designed to achieve lower weight.
The disadvantage of it lies in the very nature of is its local approach: more sophis-
ticated solutions other than increasing the thickness of highly stressed parts can not be
found by it. For example, the maximum-strength flywheel design, see section 9.3, could
not be solved by CAO.
Figure 9.4: Initial flywheel shape and radial and circumferential stress distributions
ary condition at the bore and at the rim. The circumferential stress σθ is higher than the
radial stress everywhere and has a stress peak at the bore. Because of the inhomogeneous
distribution of both the normal and the circumferential stresses, the material of the disk
of constant thickness is not economically used.
the central disk are taken to be equivalent to an average radial stress which enters the
mathematical model as a natural boundary condition. By choosing the more general
formulation - allowing for a changing thickness t - of the equilibrium of forces in the radial
direction,
t
(σrt) ,r + (σr − σθ ) + trω 2 ρ = 0 , (9.1)
r
σr = σθ = σ , (9.2)
the general equilibrium equation is reduced to a differential equation where only the thick-
ness remains as a dependent variable
ρω 2
t,r +t0 r=0, (9.3)
σ
and the shape of the evenly stressed flywheel is defined by the solution
ρω 2 2
t + t0 e− 2σ
r
(9.4)
It must, however, be kept in mind that the model is based on radial equilibrium (9.1)
only, wherefore it is accurate only if the reference thickness t0 is small compared to the
diameter of the disk. Then, the normal stress in the axial direction σz and the shear stress
τrz can be neglected. As many flywheels possess a central bore for the purpose of fixing
them on a rotating shaft, the study of Stodola’s problem provokes the question whether
a shape can be found for which an almost even stress distribution exists in the presence
of a central bore and in the absence of external radial stress. An answer to this question
is found by using a numerical shape optimization scheme where a FEM-based structural
model provides the system equations.
t0+ , t t0
linearly on the number of elements in the radial direction but is independent of the number
of elements in the axial direction. The position of the nodes on the element sides is always
on the straight line connecting the respective corner nodes. The objective function is min-
imized using the method of feasible directions according to the textbook of Vanderplaats
[20] as well as exact gradient information from the objective and the constraint functions.
Additional corrections to ensure that the equality constraints remain satisfied are made
after each line search. Quadratic approximations of the objective function are used for
quickly finding the minimum value along the current search directions as described in the
textbooks of Reklaitis et al. [33] and and Baier et al. [98].
Figure 9.7: Shape optimized after maximum stress criterion and stress distributions
evenly distributed circumferential stress. Since the circumferential stress has higher val-
ues everywhere in the initial configuration, it is identical to the first principal stress and
minimizing the objective represents the engineering aspect of seeking a spatially constant
probability of brittle failure in the flywheel. Fig. 9.8 shows the shape corresponding to
the most evenly distributed von Mises stress. The stress distribution obtained here comes
closest to a spatially constant probability of yield failure in the flywheel. Thus, it may
be concluded that local details of the shape depend on the choice of stress for the objec-
tive function while global features such as the existence of regions I, II, and III remain
Figure 9.8: Shape optimized after yield stress criterion and stress distributions
the same. Both the shape of region II and the almost even distribution of the stresses
within it bring to mind Stodola’s evenly stressed turbine disk. The deviations from per-
fect constancy of the stresses are explained by the effects of steep thickness changes on the
two-dimensional stress equilibrium. The hypothesis is proposed that region I alleviates
the stress perturbations of the central hole and that region III introduces radial stresses
on the outer edge of region II which are used in Stodola’s model to simulate the forces
exerted by the blades on the rotating turbine disk.
d is c r e te r in g w ith d is c r e te r in g w ith
a r e a A 1 fo r r e g io n I a r e a A 2 fo r r e g io n II
a x is o f r o ta tio n
S t o d o la 's d is k o f v a r y in g
th ic k n e s s fo r r e g io n II
Figure 9.9: Mechanical model consisting of Stodola’s disk and two discrete rings repre-
senting regions I and III
a thickness distribution as given by Stodola [24]. The inner and outer rings at positions
r1 and r2 have the cross-sectional areas A1 and A2 , respectively.
The line load Nr resulting from the integration of the radial stress σr over the thickness
is assumed to be positive, and the circumferential stresses σθ in the rings are
r1
σθ1 = ρω 2 r12 + Nr1
A1
r2
σθ2 = ρω 2 r22 − Nr2
A2
where the first terms on the right-hand side give the circumferential stress in the inde-
pendently rotating ring and the second terms make corrections for the tensile radial line
load which tends to pull the inner ring apart and the outer ring together. The rotational
symmetry of the problem automatically satisfies the condition of force equilibrium in the
circumferential direction. The condition of force equilibrium in the radial direction re-
quires that the radial line loads be continuous at the interfaces between the disk and the
rings. If the material of the flywheel is homogeneous, the kinematics are satisfied if the
circumferential stress is also continuous. The particular shape of the disk as defined by
Stodola’s equation guarantees the constant stress distribution indicated in (9.2) as long as
the conditions at the discrete rings
Nri = σt|ri , i = 1, 2 (9.8)
apply. This determines the cross-sectional areas of the rings, A1 and A2 ,
σr1 ρω 2 r 2
− 2σ 1
A1 = t 0 e
σ − ρω 2 r12
σr1 ρω 2 r 2
− 2σ 2
A1 = t 0 e
σ − ρω 2 r22
The cross-sectional areas of the rings are positive as long as σ is higher than the circum-
ferential stress due to the inertia body forces that would arise in the rotating inner ring
as an isolated system, and also lower than the circumferential stress that would arise in
the rotating outer ring:
ρω 2 r12 ≤ σ ≤ ρω 2 r22 . (9.9)
Then the stress σ, acting in the disk, tends to pull the inner ring apart in addition to
its own inertia effect and to restrain the outer ring against the inertia effect acting on it.
Thus, the stress σ can be chosen freely as long as it lies between the bounds defined by
(9.9). In accordance with the finite-element shape optimization described above it was
decided to fix the mass M as well as the rotational energy U of the disk to the respective
values of the initial shape (IS),
Z r2 2
2 2 − ρω r 2
M = tIS πρ r2 − r1 = 2πρ r1 A1 + r2 A2 + t0 re 2σ , (9.10)
r1
Z r2
1 3 − ρω
2
r2
U = tIS πρω 2 4 4 2 3 3
r2 − r1 = 2πρω r1 A1 + r2 A2 + t0 r e 2σ . (9.11)
4 r1
This can be achieved by adjusting the reference thickness t0 and the constant stress value
σ. This stress value will always be lower than the maximum circumferential stress at the
edge of the hole in the initial design. The practical significance of the optimized shape
lies in the wider choice of less expensive materials that can be used for the flywheel. The
integrals in (9.10) and (9.11) can not be solved algebraically. They can, however, be
eliminated by using the identity
Z r2 Z r2
ar2 r2 ar2 r2 2 ρω 2
re dr ≡ e −a r3 ear dr , a=− . (9.12)
r1 2 r1 r1 2σ
Combining (9.10), (9.11), and (9.12) leads to the stress σ,
U 1
= ρω 2 r22 + r22 .
σ=ρ (9.13)
M 4
The maximum circumferential stress at the edge of the hole of the initial configuration
depends on the material’s Poisson’s ratio ν,
1
σθIS = ρω 2 r22 (3 + ν) + r12 (1 − ν) .
(9.14)
4
Its value, even when neglecting the influence of Poisson’s ratio, can be up to three times
higher than that of the constant stress of the optimized flywheel.
Substituting the result given in (9.13) into the exponent of the thickness function of the
disk,
2r 2
ρω 2 2 M ω2 2 −
tr = t0 e− = t0 e−
r r 2 +r 2
r2
2σ 2U = t0 e 1 , (9.15)
reveals that the diameters of the rim and the central hole uniquely determine the shape
of the optimized disk as well as the ring areas,
2r 2
r2 + r12 − 2 12
A1 = t0 22 r1 e r2 +r1
r2 − 3r12
2r 2
r2 + r2 − 2 22
A2 = t0 22 12 r2 e r2 +r1 . (9.16)
3r2 − r1
Under the constraints of having the same mass and rotational energy as a reference flywheel
of constant thickness, a shape for constant stressing can be found only when the outer
diameter is at least 1.73205 times larger than the diameter of a central hole. Introducing
the ratio α of the two radii
r1 1
α= , 0≤α≤ √ . (9.18)
r2 3
yields more transparent equations for the ring areas and the thickness function of the disk,
1 + α2 2
− 2α 2
A1 = t 0 r1 e 1+α ,
1 − 3α2
1 + α2 − 2
A2 = t 0 2
r2 e 1+α2 ,
3−α
2r/r2
−
t(r) = t0 e 1+α2 . (9.19)
The ring areas depend on the diameter ratio as well as linearly on the absolute size of the
flywheel. The thickness function is independent of the diameter of the flywheel. Thus,
the ratio α uniquely determines the overall shape of the optimized flywheel. In order to
determine the remaining model parameter t0 it is necessary to evaluate one of the integrals
in (9.10) or (9.11), for instance, by using a numerical integration scheme.
Figure 9.10 gives the areas of the two rings and the disk over the full permitted range
of α as (9.18) indicates. The considered initial-shape flywheels have a unit radius and a
thickness of one hundredth.
In the case of very small values of α, the disk and the outer ring A2 dominate the area
of the optimized flywheel. With increasing α, the areas of the inner ring A1 and the disk
quickly increase and decrease, respectively. In the case of values of α higher than 0.3 the
inner ring tends to take up the major portion of the total area. The area of the outer ring
decreases slowly with increasing α. As α approaches the upper limit, the whole design
degenerates in the sense that all the available area moves into the inner ring.
The shapes according to the simplified model are visualized on the left-hand side of Fig.
Figure 9.10: Relative Values of the areas of the inner and outer rings and the disk
Figure 9.11: Optimum shapes for different values of the diameter ratio. Source: [2]
9.11 where ring areas are represented by solid circles. On the right-hand side of Fig.
9.11 the results of the numerical two-dimensional FEM analysis, based on the same input
data, are shown. It appears that the shape results agree well for lower values of α up
to 0.3, which seems remarkable in view of the crudeness of the simplified model. The
main difference in shape is that the extra areas at the inner and the outer rim are more
smoothly distributed according to the FEM model than according to the simplified model.
The two models agree in the extra area at the outer rim remaining fairly constant, the
extra area at the inner rim increasing quickly, and the thickness of the disk decreasing, all
with increasing values of α.
The reason for the degeneration of the agreement between the two models with increasing
values of α becomes obvious when examining the stress plot in Fig. 9.7. The radial stress
has to vanish at the boundaries causing a perturbation to the even stress distribution. The
perturbation is modeled well by the FE method but is not considered by the simplified one-
dimensional model. Thus, the shape plotted in Fig. 9.7 is optimal although the stresses
are not even throughout the domain. Equation (9.13) also implies that the capacity of
an evenly stressed flywheel for storing kinematic energy per unit mass equals the specific
strength of the material,
U σ
= σ∗ = . (9.20)
M ρ
This provides an immediate estimate of the energy storage capacity per unit mass at
optimum shape for any given material.
nT = 2n , (9.21)
and the sketch in Fig. 9.12 illustrates the equation with n = 4. The moderate-sized sample
H = 0 .0 0 H = 0 .2 5 H = 0 .5 0 H = 0 .7 5 H = 1 .0 0
problem given in Fig. 9.14 uses n = 300 finite elements and the number of unconstrained
topologies is nT = 2.037 × 1090 . Of course, a mass, or volume, constraint reduces the
number of feasible topologies drastically. Let n continue to give the mesh size and introduce
m to give the number of elements that must be filled with material due to a specified
average density in the design space. Then the number of combinations is only
n!
nT = (9.22)
m! (n − m)!
and Fig. 9.12 illustrates how, for instance, a mean density ρ = 0.25 reduces the number of
solutions from 16 to 6. The number of combinations for the sample problem with n = 300
and specified material density , giving m = 75, is nT = 9.796 × 1072 . A further constraint
arises from the consideration that the topology must at least connect the boundaries
where displacements are prescribed with those where non-vanishing stresses or forces are
applied. Design solutions violating that requirement are called illegal. A simple formula
for obtaining the legal solutions for any mesh does of course not exist but the illustrative
sample shown in Fig. 9.13 obtains that, under the mass constraint ρ̄ = 0.5, only four of the
924 design solutions after (9.22) are legal. Here it would appear attractive to investigate
Figure 9.13: The legal (top) and some illegal (bottom) topologies with 4 by 3 elements
the possibilities of stochastic or even exhaustive search provided that illegal topologies
could be systematically identified and kept from being numerically evaluated.
However, the minimum compliance problem has been proven to be convex and the well-
known homogenization approach of Bendsøe and Kikuchi [12] uses a distributed-function
parameterization rendering the objective function continuous. Each element is equipped
with a variable density function whose minimum value is zero, representing a void, and
whose maximum value represents the density of the of the massive material. Between those
bounds are intermediate states of more or less thinned out material. These intermediate
states are not wanted in the resulting topology, or layout, but accepted in order to enable a
continuous optimization process. The final design should have only material-filled or void
elements and the filled elements present the generated shape of the part to be designed.
Then the best design solution can be found quickly by using some method provided by
mathematical programming.
Figure 9.14: Topology optimization example: initial (a), intermediate (b through e), and
final (f) states for an average density ρ̄ = 0.252
Fig. 9.14(a), is high enough to reveal the deformation of the finite-element mesh with the
highest displacement at the node where the force attacks. The displacement scale factor is
the same for all plots and it can seen that the mesh deformation reduces quickly through
the optimization process. The mass density is presented by gray level. Unfortunately, there
are not enough different gray levels to portray the density differences between elements
in more detail. Nevertheless, it can be seen that already after the first iteration step a
differentiated mass density distribution with concentrations at the corners of the clamped
boundary emerges, Fig. 9.14(b). The state after 20 iterations, Fig. 9.14(b), reminds of a
sandwich design where less strained regions are already thinned out. After 40 iterations
emerges some differentiation in the interior in terms of diagonally arranged members,
Fig. 9.14(d). The intermediate crossing-members design, emerging after 60 iterations and
shown in Fig. 9.14(e), dissolves again in favor of the simple diagonal member pattern of
the final design obtained after 83 iterations and shown in Fig. 9.14(f).
The load vector r remains constant. The displacements ũ depend on the state of the
design parameters and require the solution of system of equations
Kũ = r (9.24)
assembled by the FEM model. The tilde reminds of the fact that the displacements are
discrete values defined on the nodal points of the finite-element mesh. Then, the objective
function calculation itself requires only to perform the scalar product (9.23).
9.4.3 Parameterization
The structural stiffness depends on the stiffness contributions of the finite elements. The
finite elements do not change their size or shape but their stiffness is controlled by Young’s
modulus E. It depends, see Fig. 9.15 in turn on the density distribution function ρ so
1 .0
N o r m a liz e d Y o u n g 's m o d u lu s
0 .8
0 .6
F
- = - r
0 .4
0 .2
0 .0
0 .0 0 .2 0 .4 0 .6 0 .8 1 .0
D e n s ity
that
E = E0 ρp , → 0 ≤ ≤ 1, (9.25)
where E0 is a reference stiffness value and is a small but finite number so that ρ is
prevented to become exactly zero for numerical reasons. The exponent in (9.25), usually
chosen from the range 3 ≤ p ≤ 4, is justified by some material modelling considerations
but also effects Young’s modulus moving against its bounds more quickly. In the original
approach, the values of the density functions ρ are directly related to the each finite
element and used as variable design variables. The lower bound signifies a void and the
upper bound signifies a region of massive material extending over the sub-domain of just
one finite element. Thus, there are as many design variables as finite elements. If the
vector x denotes the sought topology result, the density distribution function ρ(x) is
if x ∈ Ωv
ρ(x) = (9.26)
1 if x ∈ Ωs
where Ωv and Ωs denote the void and solid regions in the design space Ω, respectively. So
we have a typical case of a mesh-dependent parameterization and the number of variable
design parameters, or optimization variables, may become quite high: the relatively simple
example with its course mesh presented in section 9.4.1 already features 300 independent
optimization variables. An evenly spaced mesh, where all finite elements have the same
shape and size, reduces the numerical effort because all element matrices can be traced
back to a reference matrix K0 so that
Kk = Ek K0 (9.27)
and the reference element stiffness matrix needs to be calculated only once. The stiffness
matrix K of the whole structure, appearing in the system of equations (9.24), is thus
assembled from the element matrices by
NEL
ρpk .
X
K = E0 K0 (9.28)
k=1
Fig. illustrates how the nodal point addresses are used to connect the element stiffness
matrix entries with those of the system matrix. Finally it is important to note that the
5 ,6 1 1 ,1 2 1 7 ,1 8 2 3 ,2 4
3 ,4 9 ,1 0 1 5 ,1 6 2 1 ,2 2
1 ,2 7 ,8 1 3 ,1 4 1 9 ,2 0
F in ite E le m e n t M o d e l S y s te m S tiffn e s s M a tr ix
total mass must be constrained to some initially specified value M0 , or mean density ρ̄ in
the volume V ,
N
X EL
ρi Vi = M0 = ρ̄V , (9.29)
i=1
where Vi is the volume of one finite element, because otherwise the objective of maximum
structural stiffness would be reached by the trivial solution that the design space Ω is
completely filled with material.
The side constraints on the design variables ρi are cast into the inequality constraining
functions
gi+ = ρmin − ρi ≤ 0 (9.32)
and
gi− = ρi − 1 ≤ 0 (9.33)
Optimal topology is obtained at the saddle point of the Lagrangian and the derivative of
it with respect to the optimization variables must vanish:
∂L(ρ, Λ, λ+ , λ− )
=0 (9.35)
∂ρi
The next step is to derive the gradients of all parts of the Lagrange function.
∂ũ (T )
∂W ∂r
= r + ũ(T ) (9.36)
∂ρi ∂ρi ∂ρi
The sensitivity of the displacement solution ũ to the optimization variable ρi is described
by the sensitivity formula (6.133). Substituting it into (9.36) yields
∂K (T ) (−T )
∂W ∂r ∂r
= − ũ K r + ũ(T ) (9.37)
∂ρi ∂ρi ∂ρi ∂ρi
Because of the symmetry of the global stiffness matrix K we have that K(−T ) r = ũ and
therefore (9.36) simplifies to
∂W ∂r ∂K
= 2ũ(T ) − ũ(T ) ũ (9.38)
∂ρi ∂ρi ∂ρi
A change of the density value ρi of the ith finite element influences only the stiffness
of it and not the stiffness of any other element. If it influences, in the case of density-
dependent weight, only the load of the ith finite element and not the load of other elements,
the sensitive of the objective with respect to the variable ρi can be calculated on element
level,
∂W (T ) ∂ri (T ) ∂Ki
= 2ũi − ũi ũi , (9.39)
∂ρi ∂ρi ∂ρi
where ũi and Ki are the displacement vector and stiffness matrix of the ith finite element.
Fixed loads do not depend on the state of ρi and vanish from (9.39). With a reference
stiffness matrix K0 for a unit Young’s modulus E = 1 the stiffness matrix of the ith finite
element is obtained by
Ki = ρ4i K0 . (9.40)
The reference stiffness matrix must be calculated only once if the mesh is a rectangular
array or otherwise regular in a way so that all elements have the same geometry. This
reduces the numerical effort for the assembly of the global stiffness matrix drastically.
Returning to our gradient-calculation problem, we discover that the derivative of the
stiffness matrix can now be obtained analytically,
∂Ki
= 4ρ3i K0 . (9.41)
∂ρi
Last we assume fixed loads and insert (9.41) in the element-level sensitivity equation (9.39):
∂W (T )
= −4ρ3i ũi K0 ũi . (9.42)
∂ρi
∂hM (ρ)
= Vi (9.43)
∂ρi
∂gi+ (ρ)
= −1 (9.44)
∂ρi
∂gi− (ρ)
=1 (9.45)
∂ρi
Substituting the result into (9.46) and taking the negative of the gradient gives a search
direction along which the mass remains the same:
N EL
p
pE0 ρp−1 ρp−1
X
si = i ũTi K0 ũi + E0 T
k ũk K0 ũk (9.51)
V
k=1
Since it has been assumed that the inequality constraints are inactive, it must be discussed
what happens if they become active. Then, since they apply individually to the density
functions, the respective entries in the gradient vectors can be considered constants and
the respective entries in all gradients can be set equal to zero:
where σ̄ and ū are prescribed natural and geometric boundary conditions, respectively,
and b is distributed body force. The first term on the right-hand-side carries the same
meaning as (9.23), with the only difference being that the latter is expressed with the
discrete forces and displacements of a finite element model. The second term gives the
compliance due to prescribed displacement rather than traction and one wishes, for a stiff
structure, the reaction stress as high as possible which explains the negative sign. The
third term gives the work done by distributed body force such as self-weight (sometimes
also called dead weight). The total potential energy Π of an elastic system is the difference
between the stored elastic deformation energy U and the external work W or
Π=U −W . (9.55)
The principle of the minimum of the total potential energy requires that Π be a minimum
with respect to the displacement. Then, at the equilibrium point, Clapeyron’s theorem
states that Π = −W/2 or
Z Z Z
T
Π(u, ρ) = udΩ − σ̄ udΓ − bT udΩ (9.56)
Ω Γσ Ω
where u denotes the strain energy density function (please distinguish this from the bold
printed displacement) which for a linear elastic material is given by
u = T C (9.57)
where denotes the linearized strain tensor in vector notation and C the materials law in
matrix notation.
Jog defines the dual problem by the task of finding the Lagrange multiplier Λ associated
with the volume constraint that solves
where Z
L(Λ) = max min Π(u, ρ) − Λ ρdΩ − M0 (9.59)
ρ u Ω
From Section 5.7.4 we recall that dual formulations require the primal problem consist of
separable functions. Jog proceeds with establishing a separable approximation for L(Λ)
and approximates the total potential energy by using the reciprocal variables κ to denote
the field 1/ρ. Inventing the symbols u0 and ρ0 to denote the equilibrium displacement field
corresponding to the design ρ0 , respectively, a first-order approximation for Π by Taylor
series is given by Z
0 0 0 ∂u 1 1
Π(u , ρ) ≈ Π(u , ρ ) + − dΩ (9.60)
Ω ∂κ κ
0 ρ ρ0
The total potential energy can now be maximized without considering the constant term.
Jog next introduces the discrete representation corresponding to the finite element method.
Then, the Lagrangian can be written as a sum over the separate contributions of each finite
element, or
!
X Z ∂u 1 1
XZ
L(Λ) = max − 0 dΩ − Λ ρdΩ − M 0 (9.61)
ρ Ωi ∂κ κ0 ρ ρ Ωi
i i
where Omegai is the domain occupied by the ith finite element. Since the density values
are assumed to be constant over each element, the densities can be removed from within
the integrals, or integrals can be replaced by summations. Thus we get
X 2 1 1
Z
∂u
0
L(Λ) = max ρi − dΩ − Λρi Vi + ΛM 0 (9.62)
ρ
i
ρ0i ρi Ωi ∂ρ ρ
i i
0
The maximization is carried out by maximizing each term under the summation sign, or
Z
X
02 1 1 ∂u
L(Λ) = max ρi − dΩ − Λρi Vi + ΛM 0 (9.63)
i
ρ ρ0i ρi Ωi ∂ρi ρ 0
i
The density ρi can take on either of the two values (ρmin , ρmax ). If the maximum of the
term in square brackets occurs at the value of (ρmax ), it holds that
Z Z
02 1 1 ∂u 02 1 1 ∂u
ρi 0 − dΩ − Λρmax Vi > ρi 0 − dΩ − Λρmin Vi
ρi ρmax Ωi ∂ρi ρi
0 ρi ρmin Ωi ∂ρi ρi
0
(9.64)
and simplification leads to the condition
2
ρ0i
Z
∂u
ρi = ρmax if dΩ > Λρmin Vi (9.65)
ρmax ρmin Ωi ∂ρi ρ0i
In fact the dual problem has only the one single dual variable Λ which needs to be found for
minimizing the Lagrangian, all other variables follow from the above conditions. The dual
optimization problem requires an iterative procedure. In alternations the displacement
field is solved for given density distribution and the Lagrangian function is then minimized
with respect to the Lagrangian multiplier Λ. A line search is conducted for the Lagrange
multiplier only. Its value is increased if the current mass is higher than the specified mass,
and vice versa. Jog reports difficulty with obtaining topologies close to the optimum
arising from the fact that the linear approximations of the Lagrangian close to a reference
point permit only small changes in the design variables. As a remedy he uses a filtering
technique.
1 6 c m
5 c m
P = 4 8 0 0 0 N
Figure 9.17: Topology design problem considered by Jog. Source: Jog [99]
modulus and Poisson’s ratio E = 2.1 × 107 and 0.25, the force attacks at the domain
center, and the average mass density is 30 per cent. He utilizes the existing symmetry and
uses for the half-model a mesh of 40 × 25 finite elements with quadratic shape functions
(Q9). His solution, see Fig. 9.18, is obtained with the dual problem and here compared
Figure 9.18: Topology design solution after Jog. Source: Jog [99]
with a solution, see Fig. 9.19 obtained with the demonstration program TOP which is
based on the primal problem. The performance data, as reported by Jog and obtained
Figure 9.19: Topology design solutions with half (left) and full (right) models
with our program TOP, are summarized in Table 9.1. Using a half model with TOP yields
a design with almost the same compliance as the one obtained by Jog; it exhibits more of
the expected symmetries and the number of iterations is even smaller.
The ground structure approach is well established for the optimization of the geometry
and topology of trusses. A number of mathematical formulations for solving the problem
of finding minimum compliance is explained in the review paper [100] and the material
presented in here is selected from it. The layout of a truss structure is found by allowing a
certain set of connections between a fixed set of nodal points as active structural members
or vanishing members. Fig. 9.20 exemplifies a ground structure in two dimensions with
= >
? @
fifteen points and four pre-defined sets of connections with increasing complexity. In case
of the so-called complete ground structure, where each nodal points is connected with all
others and which is illustrated by Fig. 9.20(d), the number of bars, m, increases with the
square of the number of points, n, by
1
m = n (n − 1) . (9.67)
2
Then the number of bars is much higher than the number of degrees-of-freedom of the truss-
structure finite-element model. Please note that with the topology optimization method,
commented in Section 9.4, the numbers of density variables and degrees-of-freedom grow
only linearly with the mesh density. Also, the full ground structure produces a structural
stiffness matrix lacking any sparseness and bandedness. It should therefore not surprise
that, for the two approaches, the respective most efficient problem formulations and solu-
tion methods are not the same.
Obviously, the ground-structure approach is only a sizing problem, where the optimum
combination of the member cross-sectional areas must be found, if the areas must not
be smaller than a lower limit. Permitting zero area values combines the sizing with the
topology optimization problem. If, moreover, the positions of the nodal points of the
structure are allowed to move, shape optimization joins in and the combination of all
three disciplines is called a layout problem.
m
X
T
min r u, subject to: ti Ki u = r,
u,t
i=1
m
X
ti = V, ti ≥ 0, i = 1, . . . , m (9.68)
i=1
Here, ti symbolizes the volume of the ith bar, and ti = ai li is introduced to achieve a more
compact notation. In local coordinates, the stiffness matrix of a bar reads
" #
EA 1 −1
K= (9.69)
L −1 1
and consequently the Ki appearing in (9.68) must be normalized with respect to bar
volume: " #
E 1 −1
Ki = 2 (9.70)
li −1 1
If a non-negative lower bound is imposed on the volumes ti , the stiffness matrix remains
positive definite for all ti > 0 and the remaining problem of just adjusting the bar volumes
for maximum stiffness with respect to the defined forces has been shown by Svanberg [101]
to be convex with assured existing solutions.
The zero lower bound on the variables ti implies that bars of the ground structure can be
removed and the problem statement thus covers topology design. It also implies that the
stiffness matrix is not necessarily positive definite and that the displacement solution u can
not simply be removed from the problem by solving the finite-element-method equations.
M m
k kT
X X
k
min w r u , subject to: ti Ki uk = rk ,
u,t
k=1 i=1
m
X
k = 1, . . . , M, ti = V, ti ≥ 0, i = 1, . . . , m . (9.71)
i=1
w2 Ki
K̂i = (9.72)
..
.
w M Ki
are introduced. This allows writing problem (9.71) as
m
T
X
min r̂k ûk , subject to: ti K̂i ûk = r̂k ,
u,t
i=1
m
X
ti = V, ti ≥ 0, i = 1, . . . , m . (9.73)
i=1
m
X m
X m
X
subject to: ti Ki u = r + ti gi , ti = V,
i=1 i=1 i=1
ti ≥ 0, i = 1, . . . , m . (9.74)
m m
!
X X
+Λ ti − V + λi (−ti ) . (9.75)
i=1 i=1
∂L ∂L
= 0; → =0, (9.76)
∂u ∂t
obtain the necessary conditions
m
X m
X
ti Ki ǔ = r + ti gi
i=1 i=1
m
!
X
ǔT Ki u − 2gi = Λ − λi
i=1
λi ≤ 0; λi ti = 0; i = 1, . . . , m; Λ≥0. (9.77)
Let µ(u) denote the maximal mutual energy with self-weight uT (Ki u − 2gi ) of the indi-
vidual bars, i.e.
µ(u) = max uT (Ki u − 2gi ) |i = 1, . . . , m ,
(9.78)
and let J(u) denote the set of bars for which the mutual energy attains this maximum
level,
J(u) = i|uT (Ki u − 2gi ) = µ(u) .
(9.79)
After defining the non-dimensional element volumes t̃i = ti /V the necessary conditions
are satisfied with
ǔ = u; ti = t̃i V, i ∈ Ju; ti = 0 i ∈
/ J(u); Λ = U (u)
provided that there exists a displacement field u with corresponding set J(u) and non-
dimensional element volumes t̃i , i ∈ J(u), such that
X X X
V t̃i Ki u = r + V t̃i gi ; t̃i = 1 (9.81)
i∈J(u) i∈J(u) i∈J(u)
The reduced optimality conditions (9.81) state that a convex combination of the gradients
of the quadratic functions
1 T T
V u Ki u − gi u , i ∈ J(u) , (9.82)
2
Π=U −W , (9.83)
7 2 9
K~ 0 K~
Figure 9.21: Total potential energy visualization for only one dependent variable
W = −2Π (9.84)
which is visualized in Fig. 9.21. Here, we identify the external work W and the total
potential energy Π with
m
!T m
!
X 1 T X
W = r+ ti gi u, Π= u ti Ki u − 2r − 2 ti gi (9.85)
2
i=1 i=1
Substituting the expressions (9.85) into the equilibrium condition (9.84) gives
m
!T m
!T m
X X X
r+ ti gi u=2 r+ ti gi u− ti uT Ki u (9.86)
i=1 i=1 i=1
Since the maximal mutual energy µ(u) is a constant, the sum over the volumes of all bars
can be replaced with the total volume V of the truss. Changing the volumes (which is
traced back to changing the cross-sectional areas) gives a different truss design si but let
the maximum mutual energy µ(u) continue to correspond to the fully-stressed design ti .
These thoughts give
Xm
2rT u − V µ(u) = 2rT u − si µ(u) (9.88)
i=1
However, the other design si may have different bar volumes only for the active bars but
generally it must be assumed that volumes are assigned to bars which are inactive in the
fully stressed design ti . The mutual energy of the inactive bars is smaller than µ(u) and
therefore it holds that
m
X m
X
2rT u − si µ(u) ≤ 2rT u − si uT (Ki u − 2gi ) (9.89)
i=1 i=1
Next observe that the design si can not be at equilibrium with the here assumed dis-
placements corresponding with the fully stressed design ti . The extremum principle for
equilibrium requires that Π be a minimum or, by (9.84), W be maximum. Let the maxi-
mum of W for design si be found by adjusting a variable displacement vector w:
!T
m m m
X X 1 X
2rT u − si uT (Ki u − 2gi ) ≤ 2 max r + si gi w− si wT Ki w (9.90)
w 2
i=1 i=1 i=1
Practically the equilibrium adjustment of w is found by simply solving the elasticity prob-
lem for given design si for the displacements which are here here labelled v. Using (9.84)
again identifies
m
!T m m
!T
X X X
2 r+ s i gi v− si vT Ki v = r+ s i gi v (9.91)
i=1 i=1 i=1
The sequence (9.87) through (9.91) of equalities and inequalities proves that
m
!T m
!T
X X
r+ ti gi u≤ r+ si gi v (9.92)
i=1 i=1
which says that the best design solution is the fully-stressed design.
The optimality criterion (9.81) with its implicated fully stressed design solution allows
devising a very simple and effective search method which is known in the literature [104,
105] as optimality criterion method. The iterative method assign in each step volumes to
the bars proportionally to the respective mutual energies in order to the state of constant
mutual energy in the active bars. At iteration step k the following computations are
performed:
Pm
- find actual volumes satisfying volume constraint V k = k tki = ζik V /V k
i=1 ζi ,
Demonstration Programs
Two demonstration programs are used for exercises and are available to the students. The
data input files of both programs are explained in the following sections.
1 : Q U A D R A T I C P O L Y N O M I A L
2 : H I M M E L B L A U F U N C T I O N
3 : R O S E N B R O C K F U N C T I O N
4 : F E N T O N E A S O N F U N C T I O N
5 : C A N T I L E V E R B E A M P R O B L E M
N O I S E : F R E Q U E N C Y A N D A M P L I T U D E 0 . 0 0 E + 0 0 0 . 0 0 E + 0 0
D E F I N E V I E W I N G R E G I O N : X M I N X M A X Y M I N Y M A X
2 . 0 0 3 5 . 0 0 2 . 0 0 2 5 . 0 0
S T A R T V A L U E S : X 0 Y 0 S I Z E
5 . 0 0 0 1 0 . 0 0 0 0 . 1 0 E + 0 0
S E L E C T O N E O F T H E M E T H O D S L I S T E D B E L O W : 4
1 S I M P L E X S E A R C H
2 R E S P O N S E S U R F A C E : S P E C I F I Y V E R S I O N 0 - 1 1
3 P O W E L L S M E T H O D
4 S T E E P E S T D E S C E N T ( E X A C T L I N E S E A R C H E S ) 1
5 F L E T C H E R - R E E V E S
6 N E W T O N M E T H O D
7 M O D I F I E D N E W T O N M E T H O D
F O R W A R D ( 0 ) O R C E N T R A L ( 1 ) D I F F E R E N C E S M E T H O D 0
N U M B E R A N D D I S T R I B U T I O N O F C O N T O U R L I N E S : 1 0 0 8
D E M O N S T R A T I O N T I M E S P A N : 1
ferent test functions, where the functions 2 through 4 often appear in the literature. The
user defines a viewing region and choices of the corner coordinates are suggested in table
10.1. The table suggests also starting point coordinates and the number and distribution
quadratic polynomial
Himmelblau −6.0 6.0 −6.0 6.0 100 4
Rosenbrock −3.0 3.0 −4.0 2.0 100 4 −1.2 1.0
Fenton-Eason 0.1 3.0 0.1 3.0 100 16 0.5 0.5
Cantilver Beam 1.0 40.0 1.0 30.0 100 16 2.0 2.0
of contour lines. According to the suggested values DEMO OPT produces a series of plots
showing each iteration of the optimization process and Fig. 10.2 shows the respective last
#
! "
images. The parameter size has an effect on the initial size of the simplex, the initial lat-
tice spacing of the supporting point set of the response-surface method, and takes on the
meaning of a fixed step length when these are combined with choosing steepest-descend
search directions. Choosing the central differences method over the forward differences
method tends to increase computation time. The parameter demonstration time span can
be adjusted to view the optimization process as quickly or slowly as desired.
L E N G T H A N D H E I G H T O F D E S I G N S P A C E : X L E N G T H Y H E I G H T Y H E IG H T
1 0 . 0 1 0 . 0
X L E N G T H
N O . O F E L E M E N T S I N X A N D Y D I R E C T I O N S : N E L X N E L Y
1 5 0 1 5 0
S U P P . E D G E S ( 0 = F R E E , 1 = X , 2 = Y , 3 = B O T H ) : 0 0 0 0 N E L Y
S U P P . C O R N E R S ( 0 = F R E E , 1 = X , 2 = Y , 3 = B O T H ) : 3 3 3 3
E D G E T R A C T I O N S : 0 . 0 0 0 0 . 0 0 0 1
N E L X
E D G E T R A C T I O N S : 0 . 0 0 0 0 . 0 0 0 1
E D G E T R A C T I O N S : 0 . 0 0 0 0 . 0 0 0 1 E D G E 3
E D G E T R A C T I O N S : 0 . 0 0 0 0 . 0 0 0 1
E D G E 4 E D G E 2
C O R N E R F O R C E S : 1 . 0 0 0 - 1 . 0 0 0 7 2
C O R N E R F O R C E S : 1 . 0 0 0 1 . 0 0 0 7 2
C O R N E R F O R C E S : - 1 . 0 0 0 1 . 0 0 0 7 2 E D G E 1
C O R N E R F O R C E S : - 1 . 0 0 0 - 1 . 0 0 0 7 2
M I D S I D E F O R C E S : 0 . 0 0 0 0 . 0 0 0 1 C O R N E R 4 C O R N E R 3
M I D S I D E F O R C E S : 0 . 0 0 0 0 . 0 0 0 1
M I D S I D E F O R C E S : 0 . 0 0 0 0 . 0 0 0 1
M I D S I D E F O R C E S : 0 . 0 0 0 0 . 0 0 0 1
G R A V I T Y 0 . 0 0 0 C O R N E R 1 C O R N E R 2
A V E R A G E D E N S I T Y 0 . 3 2 0
E X P O N E N T P 4 . M ID S ID E 3
T O L E R A N C E R A N G E [ % ] 2 . 0
F I L T E R O F F O R O N ? 0 O R 1 1
N U M B E R O F I T E R A T I O N S 1 0 0 0 M ID S ID E 4 M ID S ID E 2
R E A D M O V I E F I L E ? 0 O R 1 0
M ID S ID E 1
by which the respective load introduction points are moved away from the boundary into
the domain. Rather time consuming optimization processes, for which the presented data
set is an example, can be viewed at much higher speed after the end of the computations.
Then, the encircled number must be set to 1 and the intermediate topology information
is read from a data file and plotted.
7.1 Metropolis Algorithm: Probability for the acceptance of a new solution . . 125
7.2 Selected mathematical programming methods and ordering scheme . . . . . 126
7.3 General architecture of an evolutionary algorithm. . . . . . . . . . . . . . . 129
9.1 Part with stress raising geometry and modelled cambium . . . . . . . . . . 167
9.2 Two-phase CAO process after Mattheck . . . . . . . . . . . . . . . . . . . . 168
9.3 Soft-Kill Option process after Mattheck . . . . . . . . . . . . . . . . . . . . 169
9.4 Initial flywheel shape and radial and circumferential stress distributions . . 171
9.5 Blueprint of a turbine disk design solution after Stodola . . . . . . . . . . . 171
9.6 Parameterization of thickness and mesh generator . . . . . . . . . . . . . . 172
9.7 Shape optimized after maximum stress criterion and stress distributions . . 173
9.8 Shape optimized after yield stress criterion and stress distributions . . . . . 174
9.9 Mechanical model consisting of Stodola’s disk and two discrete rings . . . . 174
9.10 Relative Values of the areas of the inner and outer rings and the disk . . . . 177
9.11 Optimum shapes for different values of the diameter ratio . . . . . . . . . . 177
9.12 Illustration of topologies for a mesh of 1 by 4 elements . . . . . . . . . . . . 179
9.13 The legal and some illegal topologies with 4 by 3 elements . . . . . . . . . . 180
9.14 Topology optimization example . . . . . . . . . . . . . . . . . . . . . . . . . 180
9.15 Young’s modulus and density . . . . . . . . . . . . . . . . . . . . . . . . . . 181
9.16 Illustration to the assembly of the global stiffness matrix . . . . . . . . . . . 182
9.17 Topology design problem considered by Jog . . . . . . . . . . . . . . . . . . 188
9.18 Topology design solution after Jog . . . . . . . . . . . . . . . . . . . . . . . 188
9.19 Topology design solutions with half and full models . . . . . . . . . . . . . . 188
9.20 Various ground structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
9.21 Total potential energy visualization for only one dependent variable . . . . 193
[1] Müller, S.D. Bio-inspired optimization algorithms for engineering applications. Dis-
sertation ETH No. 14719, Swiss Federal Institute of Technology Zurich, 2002.
[2] Kress, G.R. Shape optimization of a flyweel. Struct. Multidisc. Optim., 19:74–81,
2000.
[7] Kress, G., P. Naeff, M. Niedermeier, and P. Ermanni. Onsert strength design. Int.
J. of Adhesion & Adhesives, 24:201–209, 2004.
[8] Kress, G. and P. Ermanni. The onsert: A new joining technology for sandwich
structures. In Proc. International Conference on Buckling and Postbuckling Behavior
of Composite Laminated Shell Structures, Eilat, Israel, 2004.
[9] Zehnder, B. and P. Ermanni. A methodology for the global optimization of laminated
composite structures. Composite Structures, 72(3):311–320, 2006.
[10] Zehnder, B. and P. Ermanni. Optimizing the shape and placement of patches of
reinforcement fibers. Composite Structures, 77(1):1–9, 2007.
[11] Gätzi, R., M. Uebersax, O. König. Structural optimization tool using genetic al-
gorithms and ansys. Proc. 18. CAD-FEM User’s Meeting, Internationale FEM-
Technologietage, Graf-Zeppelin-Haus, Friedrichshafen, 2000.
[12] Bendsœ, M.P., N. Kikuchi. Generating optimal topologies in structural design using a
homogenization method. Computer Methods in Applied Mechanics and Engineering,
71:197–224, 1988.
[13] Michell, A.G.M. The limits of economy of material in frame-structures. Phil. Mag.
(Series 6), 8:589–597, 1904.
[14] N. Zehnder. Global Optimization of Laminated Structures. PhD thesis, ETH Zurich,
2008. DISS. ETH. NO. 17573.
[15] Ruge, M. Entwicklung eines flüssigkeitsgekühlten Polymer-Elektrolyt-Membran-
Brennstoffzellenstapels mit einer Leistung von 6.5 kW. Ph.D thesis, Swiss Federal
Institute of Technology, VDI Verlag GmbH, Reihe 6, Nr. 494, Fortschrittberichte
VDI, 2003.
[16] Schmid, D. Entwicklung eines Brennstoffzellenstapels für portable Aggregate unter-
schiedlicher Leistungsbereiche. Ph.D thesis, Swiss Federal Institute of Technology,
VDI Verlag GmbH, Reihe 6, Nr. 500, Fortschrittberichte VDI, 2003.
[17] Evertz, J. and M. Günthart. Structural concepts for lightweight and cost-effective
end plates for fuel cell stacks. In In 2nd European PEFC Forum, Lucerne, Switzer-
land, 2003.
[18] H. Eschenauer, H., N. Olhoff, W. Schnell. Applied Structural Mechanics. Springer,
1997.
[19] Courant, R. Variational methods for the solution of problems of equilibrium and
vibrations. Bull. Am. Math. Soc., pages 1–23, 1943.
[20] Vanderplaats, G.N. Numerical Optimization Techniques for Engineering Design:
with Applications. McGraw-Hill: Series in Mechanical Engineering, 1984.
[21] Schuldt, S.B. A method of multipliers for mathematical programming methods with
equality and inequality constraints. JOTA, 17:155–162, 1975.
[22] Schuldt, S.B., G.A. Gabriele, R.R. Root, E. Sandgren, and K.M. Ragsdell. Applica-
tion of a new penalty function method to design optimization. ASME J. Eng. Ind.,
99:31–36, 1977.
[23] Reklaitis, G.V., A. Ravindran, K.M. Ragsdell. Engineering Optimization - Methods
and Applications. John Woley and Sons, 1983.
[24] Stodola, A. Dampf- und Gasturbinen. Springer, 1924.
[25] Jones, R.M. Mechanics of Composite Materials. Hemisphere Publishing Corporation,
1975.
[26] Schmit, L.A., R.H. Mallet. Structural synthesis and design parameters. Hierarchy
Journal of the Struct. Division, Proceedings of the ASCE, 89:269–299, 1963.
[27] Olhoff, N., J.E. Taylor. On structural optimization. J. Applied Mechanics, 50:1139–
1151, 1983.
[28] Fritsche, D. Lasteinleitungselemente maximaler verbindungsfestigkeit f/”ur sand-
wichbauteile. Student project thesis no. winter semester 2000/2001, Structure Tech-
nologies, Institute of Mechanical Systems, ETH Zürich.
[29] Naeff, P. Experimentelle verifikation strukturell geklebter lasteinleitungselemente.
Student project thesis no. 02-115 summer semester 2002, Structure Technologies,
Institute of Mechanical Systems, ETH Zürich.
[31] Kuhn, H.W., A.W. Tucker. Nonlinear programming. Proc. Second Berkeley Symp.
on Math. Statist. and Probability, pages 481–492, 1951. J. Neyman, Ed.
[35] Powell, M.J.D. An efficient method for finding the minimum of a function of several
variables without calculating derivatives. Computer J., 7:155–162, 1964.
[36] Cauchy, A. Method generale pour la resolution des systemes d’equations simultanees.
Compt. Rend. Acad. Sci., 25:536–538, 1847.
[37] Spendley, W., G.R. Hext, F.R. Himsworth. Sequential application of of simplex
designs in optimization and evolutionary operation. Technometrics, 4:441–461, 1962.
[38] Shewchuck, J.R. An introduction to the conjugate gradient method without the
agonizing pain. School of Computer Science, Carnegie Mellon University, Pittsburgh,
PA, 15213, Aug. 4, 1994.
[39] Fletcher, R., C.M. Reeves. Function minimization by conjugate gradients. Computer
J., 7(5):149–154, 1964.
[40] Polak, E., G. Ribière. Note sur la convergence de méthodes de directions conjuguées.
Revue française d’informatique opérationelle, série rouge, 3(1):35–43, 1969.
[41] Navon, I.M., D.M. Legler. Conjugate-gradient methods for large-scale minimization
in meteorology. Mon. Wea. Rev., 115:1479–1502, 1987.
[42] Zangwill, W.I. Minimizing a function without calculating derivatives. Computer J.,
10:293–296, 1967.
[43] Brent, R.P. Algorithms for Minimization without Derivatives. Prentice-Hall, Engle-
wood Cliffs, NJ, 1973.
[44] Venter, G. Non-dimensional response surfaces for structural optimization with un-
certainty. Dissertation, University of Florida, 1998.
[45] Kress, G. and P. Ermanni. Comparison between newton and response surface meth-
ods. Struct. Multidisc. Optim., accepted for publication, 2004.
[46] Wang, G.G. Adaptive response surface method using inherited latin hypercube
design points. Journal of Mechanical Design, 125:210–220, 2004.
[47] McKay, M.D., R.J. Bechmann, and W.J. Conover. A comparison of three methods
for selecting values of input variables in the analysis of output from a computer code.
Technometrics, 21(2):239–245, 1997.
[48] Levy, A.V., S. Gomez. The tunneling method applied to global optimization. Proc.
SIAM Conf. on Num. Optimization, pages 213–244, 1984.
[50] D.H. Wolpert and W.G. Macready. No free lunch theorems for optimization,. IEEE
Transactions on Evolutionary Computation, 1:67–82, 1997.
[52] Fogel, L.J., A.J. Owens, A.J. Walsh. Artificial Intelligence Through Simulated Evo-
lution. Wiley, New York, 1966.
[54] Holland, J.H. Adaptation in Natural and Artificial Systems. Michigan Press, Ann
Arbor, MI, 1975.
[55] J.H. Holland. Adaptation in Natural and Artificial Systems. University Michigen
Press, Ann Arbor, MI, 1975.
[56] D.E. Goldberg. Genetic Algorithms in Search, Optimization and Machine Learning.
Addison-Wesley Publishing Company, Inc., 1989.
[59] D. Fogel. Evolving Artificial Intelligence. PhD thesis, University of California, San
Diego, CA, 1992.
[61] T. Bäck. Evolutionary Algorithms in Theory and Practice. Oxford Univerity Press,
New York, 1996.
[65] N.J. Radcliffe and P.D. Surry. Formal memetic algorithms. Technical report, Edin-
burgh Parallel Computing Centre, 1994.
[67] J. E. Gordon. The new science of strong materials: or, Why you don’t fall through
the floor. Pelican books. Penguin, 1991.
[75] H. Fukunaga and H. Sekine. Stiffness design method of symmetric laminates using
lamination parameters. AIAA Journal, 30:2791–2793, 1992.
[76] C.G. Diaconu, M. Sato, and H. Sekine. Layup optimization of symmetrically lami-
nated thick plates for fundamental frequencies using lamination parameters. Struc-
tural and Multidisciplinary Optimization, 24(4):302–311, 2002.
[81] R. Le Riche and R.T. Haftka. Improved genetic algorithm for minimum thickness
composite laminate design. Composites Engineering, 5(2):143 – 161, 1995.
[83] Z. Gürdal, R. Haftka, and P. Hajela. Design and Optimization of Laminated Com-
posite Materials. John, 1999.
[86] A. Rama Mohan Rao and N. Arvind. A scatter search algorithm for stacking sequence
optimisation of laminate composites. Composite Structures, 70(4):383–402, 2005.
[87] J.-S. Kim. Development of a user-friendly expert system for composite laminate
design. Composite Structures, 79(1):76 – 83, 2007.
[89] J. Stegmann and E. Lund. Discrete material optimization of general composite shell
structures. International Journal for Numerical Methods in Engineering, 62(14):2009
– 2027, 2005.
[92] Shahriar Setoodeh, Mostafa M. Abdalla, and Zafer Grdal. Design of variablestiffness
laminates using lamination parameters. Composites Part B: Engineering, 37(45):301
– 309, 2006.
[96] N. Zehnder and P. Ermanni. A methodology for the global optimization of laminated
composite structures. Composite Structures, 72(3):311–320, 2006.
[97] Mattheck, C. Design in der Natur - der Baum als Lehrmeister Bd. 1. Rombacher
Ökologie, 1993.
[98] Baier, H., Ch. Seeßelberg, B. Specht. Optimierung in der Struktrumechanik. Vieweg,
1994.
[99] Jog, C.S. A robust dual algorithm for topology design of structures in discrete
variables. Int. J. Numer. Meth. Engng, 50:1607–1618, 2001.
[100] Bendsøe, M.P., A. Ben-Tal, J. Zowe. Optimization methods for truss geometry and
topology design. Structural Optimization, 7:141–159, 1994.
[102] Diaz, A., Bendsøe, M.P. Shape optimization of structures for multiple loading con-
ditions using a homogenization method. Structural Optimization, 4:17–22, 1992.
[103] Taylor, J.E. Maximum strength elastic structural design. Proc. ASCE, 95:653–663,
1969.
[104] Olhoff, N., Taylor, J.E. On structural optimization. J. Appl. Mech., 50:1134–1151,
1983.
[105] Rozvani,. Structural Design via Optimality Criteria. Dordrecht: Kluwer, 1989.
[107] H.R. Schwarz. Methode der finiten Elemente: eine Einführung unter besonderer
Berücksichtigung der Rechenpraxis. Leitfäden der angewandten Mathematik und
Mechanik. Teubner, 1991.
[109] A. E. H. Love. The small free vibrations and deformation of a thin elastic shell.
Philosophical Transactions of the Royal Society of London. A, 179:pp. 491–546, 1888.
[110] M. Iura and S. N. Atluri. Formulation of a membrane finite element with drilling
degrees of freedom. Computational Mechanics, 9(6):417–428, 1992.
[111] R. D. Cook. Four-node flat shell element: Drilling degrees of freedom, membrane-
bending coupling, warped geometry, and behavior. Computers & Structures,
50(4):549–555, 1994.
Within this section, a short overview of the Finite Element Method (FEM) is given. The
content is mainly based on the textbooks of Cook et al. [69], Reddy [106] and Schwarz
[107]. The FEM is a numerical method for solving partial differential equations. It can
be applied to a wide field of physical problems such as structural analysis, heat transfer,
magnetic fields, flow processes and many more. The geometry of the physical problem is
discretized by dividing the domain into smaller parts called the finite elements. A finite
element has a domain, a boundary and so-called nodes on which the degrees of freedom
are defined. While the partial differential equations can only be solved analytically for
simple geometries, there is no geometric restriction using the FEM. Moreover, there is
no restriction to boundary conditions or material properties which makes the method ap-
plicable to any continuum mechanical problems. There exist numerous commercial finite
element codes which are still enhanced and refined today. While a deeper understanding of
the method requires some effort, the usage of FEM-codes is possible with little knowledge
of the method and the underlying problem. However, the consequences of an incorrect
application ”may range from embarrassing to disastrous” [69].
Within this thesis, the further explanations are focused on the FEM for structural anal-
ysis problems. The solution of the elasticity problem demands for the fulfillment of the
fundamental equation
LT CLu + f = ρü (A.1)
within a given domain Ω taking into consideration the prescribed surface stresses σ̂ and
displacements û on the surface Γ as well as the internally acting body forces f . As men-
tioned, an exact solution of the problem can only be found for simple geometries and
simple boundary conditions, e.g. a bar under uniaxial load. Generally, the FEM provides
only approximate solutions. Accepting a sufficient numerical effort, the obtained solution
may come close to the true solution.
is defined as
d ∂L ∂L
− =0 (A.2)
dt ∂ ẋ ∂x
where L denotes the Lagrangian and x the generalized coordinates. The Lagrangian L is
defined as
L = T − Π = T − (U − W ) (A.3)
where T denotes the kinetic energy and Π the potential energy. The Hamilton’s principle
for an elastic body is represented with equation
Z t2 Z t2
δ [T − (U − W )] dt = δ Ldt = 0 (A.4)
t1 t1
which calls the time integral of the Lagrangian L to be stationary. The fundamental lemma
of the calculus of variation shows that solving the Lagrange equations is equivalent to find-
ing a solution of the Hamilton’s principle. However, taking advantage of the Lagrangian
mechanics transforms the equation of motion to its weak form which is needed for finite
element formulation.
The potential energy Π is the sum of the deformation energy U and the negative of the
potential of the external forces W . If the displacements u are small, which is required for
the linear elasticity problem, the velocities can be approximated with displacement time
derivatives u̇. The kinetic energy can thus be formulated as the integral over the domain Ω
of the density ρ and the scalar product or the velocities u̇.
Z
1
T = ρu̇T u̇dΩ (A.5)
2 Ω
The deformation energy U is defined as the domain integral of the scalar product of
stresses σ and strains ε. Z
1
U= εT σdΩ (A.6)
2 Ω
and the potential of the external forces W is dependent on the body forces f and the
surface stresses σ̂. Consider that the surface stresses are integrated over the surface Γ.
Z Z
W = f udΩ + σ̂ T udΓ
T
(A.7)
Ω Γ
The Lagrangian L is formulated for the entire domain Ω. In order to get the FEM for-
mulation, the domain must be discretized into smaller sub-domains Ωe , namely the finite
elements. The integration over the domain Ω is replaced with a summation of the integrals
over the sub-domains Ωe . Additionally, local approximation functions φ, which are also
called shape functions, are defined. They map the finite element nodal displacements ũ,
which represent also the degrees of freedom, to the continuous displacements u. The nodal
velocities ũ˙ and the accelerations ũ
¨ are mapped analogously.
u ≈ φT ũ , u̇ ≈ φT ũ˙ , ü ≈ φT ũ
¨ (A.8)
A strain-displacement matrix B, which contains the spacial derivatives of the local shape
functions, is built by applying the differential matrix operator L to the shape functions φ.
The total Lagrangian for the discretized system then takes the discrete form
X 1 Z X Z
1
L= ˙ T T ˙
ρũ φφ ũdΩe − T T
ũ B CBũdΩe
nelm
2 Ωe nelm
2 Ωe
X Z X Z
+ f T φT ũdΩe + σ T φT ũdΓe (A.11)
nelm Ωe nelm Γe
Finally, the unknown nodal displacements ũ are taken as generalized coordinates wherefore
the Lagrange-Euler equation (A.2) is modified to
d ∂L ∂L
− =0 (A.12)
˙
dt ∂ ũ ∂ ũ
The evaluation of the Lagrangian formalism leads to the equation of motion for the dis-
cretized linear system.
X Z X Z
T¨ T
ρφφ ũdΩe + B CBũdΩe
nelm Ωe nelm Ωe
X Z X Z
− φf dΩe − φσ̂dΓe = 0 (A.13)
nelm Ωe nelm Γe
The equation is then rearranged so that terms with no dependence on the nodal
displacements ũ are brought to the right hand side.
X Z X Z
T
B CBdΩe ũ + T ¨
ρφφ dΩe ũ
nelm Ωe nelm Ωe
X Z X Z
= φf dΩe + φσ̂dΓe (A.14)
nelm Ωe nelm Γe
The sums on the left side represent the global stiffness matrix K and the global mass
matrix M while the right hand side represents the load vector r. With these symbols,the
basic problem of the FEM for the linear elastic case is written as
¨=r
Kũ + Mũ (A.15)
¨ vanish so that the equation simplifies to
Assuming a static problem, the accelerations ũ
Kũ = r (A.16)
which is a simple linear equation system for the unknown displacements ũ. The equation
of motion can also be used for the determination of the harmonic eigenfrequencies of the
structural system. There, it is assumed that the load vector r is zero. Taking advantage
of the harmonic approach, the equation of motion can be transferred into an eigenvalue
problem with the unknown eigenvalues λ and eigenvectors Φ, whereas λ is equal to the
square of the angular frequency ω.
The combination of the equation of motion (A.15) and the harmonic approach (A.17, A.18)
yields to the eigenvalue problem for the harmonic vibration.
K − ω2M Φ = 0
(A.19)
Considering equation (A.14), the global system matrices can be extracted directly. The
global stiffness matrix K is defined as the sum of the element stiffness matrices
X Z
T
K= B CBdΩe (A.20)
nelm Ωe
Analogously, the global mass matrix M and the element mass matrices m are defined as
X Z
T
M= ρφφ dΩe (A.22)
nelm Ωe
and Z
m= ρφφT dΩe (A.23)
Ωe
The load vector r contains the internal body forces and the prescribed forces on the surface
of the structure. The shape functions φ distribute the loads to the nodes.
X Z X Z
r= φf dΩe + φσ̂dΓe (A.24)
nelm Ωe nelm Γe
The summations must consider the connectivity of the respective element with the globally
numbered mesh nodes.
u = u0 + z θy (A.25)
v = v0 + z (−θx ) (A.26)
The chosen conventions for the derivation of the finite shell element are shown in
Figure A.1. However, an alternative convention is feasible as well. Taking advantage
z,w
z,w
zθy
θx y,v
θy Ny x,u
Mxy θy P z w
My Nxy
x,u Mx
Nxy
Nx Mxy
where ε0 denotes the mid-plane or membrane strains and κ the plate curvatures. Thin
plate theory assumes the strain distribution to be linear through the thickness. The out-
of-plane deflections w are connected to the in-plane displacement u, v with the transverse
shear which is
∂w ∂u
γxz +
∂x ∂z
= ∂v ∂w (A.28)
γyz
+
∂z ∂y
Alternatively it can be expressed in terms of θx and θy taking advantage of the derivatives
of the kinematic relations (A.25) and (A.26). The derivatives yield
∂u
= θy (A.29)
∂z
∂v
= −θx (A.30)
∂z
which leads to
∂w
γxz + θy
∂x
=
∂w
(A.31)
γyz
− θx
∂y
Taking advantage of the Finite Element Method (FEM) (see Appendix A), the strains
are composed of the matrix multiplication of the differential operator L, the shape func-
tions φ and the nodal displacements ũ (see equation (A.10)).
Since membrane (m), bending (b) and transverse shear (s) parts are differently dependent
on the out-of-plane-coordinates z (see equation (A.27)) and on different material laws,
respectively, the integration over the shell thickness is separated. The membrane stiffness
matrix assuming a homogeneous material yields
Z Z t Z
2
km = BTm Cf Bm dzdA = t BTm Cf Bm dA (A.34)
A − 2t A
with
u v
∂φ
∂x 0
∂φ (A.35)
Bm =
0
∂y
∂φ ∂φ
∂y ∂x
and
1 ν 0
E ν 1
Cf = 0 (A.36)
(1 − ν)2 1−ν
0 0 2
with
θx θy
∂φ
0
∂x
∂φ (A.38)
Bb = − 0
∂y
∂φ ∂φ
−
∂x ∂y
w θ x θy
∂φ
0 φ (A.41)
∂x
Bs =
∂φ
−φ 0
∂y
The total element stiffness matrix k is the addressed summation of the single parts which
means, that the degrees-of-freedom indicated at the top of equations (A.35), (A.38) and
(A.41) must be respected.
k = k m + kb + ks (A.42)
A general shell element must have 6 degrees-of-freedom in order to be feasible for spatial
3D-modeling. However, the formulation above covers only the 5 d.o.f. u, v, w, θx and θy
but not θz The so called drilling degree-of-freedom θz must be introduced artificially (see
[110, 111]) in order to make the equation system solvable. Considering equations (A.34),
(A.37) and (A.39), it becomes obvious that no re-evaluation of the integrals is needed if the
thickness t is changed since it is contained explicitly. The thickness can be adapted with
very low additional computational cost which makes the application of the shell elements
very efficient and adequate for preliminarily design.
The theory of lamination parameters introduced in Section 8.4 require the matrices Γi
containing the material invariants of orthotropic materials. Tsai and Pagano [70] accom-
plished to formulate the reduced stiffness transformation equations by a combination of
trigonometric identities and 5 material invariants Ui (which are invariant under rotations
about the z-axis). After the textbook of Jones [25] the relations are
in which
1
U1 = (3Q11 + 3Q22 + 2Q12 + 4Q66 ) (B.7)
8
1
U2 = (Q11 − Q22 ) (B.8)
2
1
U3 = (Q11 + Q22 − 2Q12 − 4Q66 ) (B.9)
8
1
U4 = (Q11 + Q22 + 6Q12 − 4Q66 ) (B.10)
8
1
U5 = (Q11 + Q22 − 2Q12 + 4Q66 ) (B.11)
8
whereas Qii are the entries of the material stiffness matrix.
with
U1 U4 0
Γ0 = U4 U1 0 (B.13)
0 0 U5
U2 0 0
Γ1 = 0 −U2 0 (B.14)
0 0 0
0 0 U2
1
Γ2 = 0 0 U2 (B.15)
2
U2 U2 0
U3 −U3 0
Γ3 = −U3 U3 0 (B.16)
0 0 −U3
0 0 U3
Γ4 =0 0 −U3 (B.17)
U3 −U3 0
These are also the matrices Γi which are employed in the lamination parameter theory in
Section 8.4.