HEURISTIC OPTIMIZATION METHODS FOR THE CHARACTERIZATION OF
DYNAMIC MECHANICAL PROPERTIES OF COMPOSITE MATERIALS
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee. This
dissertation does not include proprietary or classified information.
Klaus H. Hornig
Certificate of Approval:
Malcolm J. Crocker
Distinguished University Professor
Mechanical Engineering
George T. Flowers, Chair
Professor
Mechanical Engineering
Gerry V. Dozier
Associate Professor
Computer Science and Software
Engineering
Subhash C. Sinha
Professor
Mechanical Engineering
Joe F. Pittman
Interim Dean
Graduate School
HEURISTIC OPTIMIZATION METHODS FOR THE CHARACTERIZATION OF
DYNAMIC MECHANICAL PROPERTIES OF COMPOSITE MATERIALS
Klaus H. Hornig
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
May 10, 2007
iii
HEURISTIC OPTIMIZATION METHODS FOR THE CHARACTERIZATION OF
DYNAMIC MECHANICAL PROPERTIES OF COMPOSITE MATERIALS
Klaus H. Hornig
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon request of individuals or institutions and at their expense. The author
reserves all publication rights.
Signature of Author
05/10/2007
Date of Graduation
iv
VITA
Klaus Harald Hornig Hollstein, son of Carlos E. Hornig and Inge A. Hollstein,
was born on June 1, 1976, in Puerto Montt, Chile. He obtained his Bachelor of Science
degree (Summa cum Laude) in Acoustical Engineering from Universidad Austral de
Chile in July 2001. Working initially in acoustic room design, and projects in
environmental noise control, he joined the Mechanical Engineering Department at
Auburn University as a Research Assistant in August 2001, and became a Doctor of
Philosophy candidate in November 2005. His scientific interests cover vibration analysis,
nonlinear systems, optimization, and acoustics, in general. He is a member of the
International Institute of Acoustics and Vibration (IIAV), and the American Society of
Mechanical Engineers (ASME). He married Carolina D. Rojas on December 17, 2005.
v
DISSERTATION ABSTRACT
HEURISTIC OPTIMIZATION METHODS FOR THE CHARACTERIZATION OF
DYNAMIC MECHANICAL PROPERTIES OF COMPOSITE MATERIALS
Klaus H. Hornig
Doctor of Philosophy, May 10, 2007
(B.Sc. Universidad Austral de Chile, July, 2001)
138 Typed Pages
Directed by George T. Flowers
Generally speaking, of the fundamental dynamic mechanical properties - mass,
damping, and stiffness, damping is usually the most difficult to quantify. This is perhaps
particularly true for composite materials which tend to have substantially higher damping
than comparable isotropic materials and therefore having an accurate representation is
correspondingly more important. Accordingly, some heuristic optimization techniques for
the identification of the dynamic characteristics of honeycomb-core sandwich composite
materials have been suggested. More specifically, Genetic Algorithms (GA) and Particle
Swarm Optimization (PSO) have been used in the present work as the optimization
techniques for the system identification process. Experimental measurements of the
dynamic responses (in the form of hysteresis loops) of simply-supported composite beam
samples have been carried out, and a simplified semi-empirical mathematical model has
vi
been developed for such a system, tailored from individual experimental observations of
the dynamic behavior of the samples when they are excited at their mid-points by
sinusoidal displacement waves. The hysteresis loops that were obtained are for several
frequencies and excitation amplitudes around the first mode of vibration. The analytical
model contains four unknown system parameters, which must be identified by both
optimization techniques utilized. The performance of these optimization methods are
compared with computer-generated and experimental hysteresis loops. In addition, the
effect of noise contamination in the signals has been studied in order to assess the search
accuracy of the optimization algorithms under such conditions.
vii
ACKNOWLEDGMENTS
I would like to thank my advisor Dr. George T. Flowers for his wise guidance,
encouragement, and support in this research work, while I was pursuing my Doctoral
studies in Mechanical Engineering at Auburn University. I also would like to thank Dr.
Malcolm J. Crocker, Dr. Subhash C. Sinha, and Dr. Gerry V. Dozier, mainly for their
directions in the areas of vibration analysis, and artificial intelligence (focused on
optimization techniques), which have been very important for the fulfillment of my
degree. I also thank Dr. Winfred A. Foster for his valuable comments and corrections,
and for his great help for being my outside reader. I would like to give thanks in a very
special way to my parents, family, and friends, that without their constant support this
?journey? would have been very difficult to accomplish. And finally, I would like to
express my greatest gratitude and gratefulness to my wife, Carolina Rojas, whose
company, love, and perpetual support have been invaluable in the final steps of my
degree.
viii
Style manual or journal used: International Journal of Acoustics and Vibration
Computer software used: Microsoft Office 2003, Matlab R2006a, LabView 8.0,
Adobe Photoshop CS2
ix
TABLE OF CONTENTS
LIST OF FIGURES ........................................................................................................... xi
LIST OF TABLES.............................................................................................................xv
CHAPTER 1. INTRODUCTION ........................................................................................1
1.1. Composite Materials and Damping ..................................................................1
1.2. Scope of Study ..................................................................................................5
CHAPTER 2. LITERATURE REVIEW .............................................................................6
2.1. Hysteresis Models and Applications.................................................................6
2.2. Damping, and Measurement Techniques........................................................12
2.3. Experimental Considerations..........................................................................20
2.4. Heuristic Optimization Algorithms in the Study of Composites....................22
CHAPTER 3. HEURISTIC OPTIMIZATION METHODS .............................................26
3.1. Genetic Algorithms (GA) ...............................................................................26
3.2. Particle Swarm Optimization (PSO)...............................................................37
3.2.1. PSO: The Basic Algorithm ..............................................................40
CHAPTER 4. MATHEMATICAL MODEL DEVELOPMENT ......................................42
4.1. Experimental Setup.........................................................................................42
4.2. Mathematical Model Baseline ........................................................................47
4.3. Preliminary Experimental Observations.........................................................49
4.4. Final Model.....................................................................................................52
x
CHAPTER 5. IMPLEMENTATION AND RESULTS.....................................................54
5.1. Optimization Algorithms Implementation......................................................54
5.2. Benchmark Test of GA and PSO....................................................................58
5.2.1. Benchmark Test Results ..................................................................60
5.3. Material Samples Description and Experimental Results...............................75
5.3.1. Optimization Results from Experimental Data................................77
CHAPTER 6. DISCUSSION AND CONCLUSIONS ......................................................89
6.1. Suggestions for Future Research ....................................................................95
REFERENCES ..................................................................................................................96
APPENDIX A..................................................................................................................101
APPENDIX B ..................................................................................................................120
xi
LIST OF FIGURES
Figure 1.1. Honeycomb sandwich composite ...................................................................2
Figure 1.2. Mechanical hysteresis loops examples ...........................................................4
Figure 2.1. Plot of a hysteron of the second kind..............................................................8
Figure 2.2. Chua-Stromsmoe hysteresis loop....................................................................8
Figure 2.3. Input-output behavior of a hysteretic relay.....................................................9
Figure 2.4. Weighted summation of 3 hysteretic relays....................................................9
Figure 2.5. Geometric interpretation of the Preisach model ...........................................10
Figure 2.6. Hysteresis loops from Bouc-Wen model ......................................................12
Figure 2.7. Loss factor (?) as a function of frequency, of a 2024-T4 aluminum alloy
cantilever beam, obtained with the Zener?s damping model ........................15
Figure 2.8. Compliance FRF, and half power bandwidth method ..................................16
Figure 2.9. General behavior of hysteresis in a system, when crossing resonance (f
0
)...19
Figure 2.10. Categories of optimization techniques..........................................................23
Figure 3.1. Schematic of a parameter estimation method ...............................................26
Figure 3.2. The two versions of GA chromosome representation...................................28
Figure 3.3. Main steps for a typical Genetic Algorithm (GA) ........................................31
Figure 3.4. Pseudo-code of a typical Genetic Algorithm (GA).......................................31
Figure 3.5. Roulette wheel analogy, as a parent selection strategy.................................33
Figure 3.6. Two types of crossover points ......................................................................34
xii
Figure 3.7. Mating process, when the crossover point is between genes........................35
Figure 3.8. Mating process, when the crossover point lies over a gene..........................35
Figure 3.9. PSO method of particle position correction at j-th iteration.........................38
Figure 3.10. PSO neighborhood topology types ...............................................................39
Figure 3.11. PSO pseudo-code for synchronous implementation.....................................41
Figure 3.12. PSO pseudo-code for asynchronous implementation ...................................41
Figure 4.1. Concept drawing of the test rig for hysteresis loops generation...................42
Figure 4.2. Pictures of the test rig for hysteresis loops generation .................................43
Figure 4.3. Experimental setup schematic, for hysteresis loops generation....................44
Figure 4.4. Pictures of the experimental setup for hysteresis loops generation ..............45
Figure 4.5. Compliance FRF comparison, between a mid-point excited beam, and
an equivalent 1-DOF system.........................................................................48
Figure 4.6. Observed behavior of a sandwich composite beam; plot of f
n
vs. X
0
...........49
Figure 4.7. Observed behavior of a sandwich composite beam; plot of ? vs. X
0
............50
Figure 4.8. Observed behavior of a sandwich composite beam; plot of ? vs. f at
different excitation amplitudes .....................................................................50
Figure 4.9. Equivalent experimental mechanical system, considering the stiffness
of the supports...............................................................................................51
Figure 4.10. Overlay plot example of the compliance frequency response of a
sandwich composite beam sample around its first resonance,
at different excitation amplitudes..................................................................52
Figure 5.1. Comparison of hysteresis loops for error calculation
(displacement loading)..................................................................................56
Figure 5.2. Scatter plots from the results of the Synth set, for ?
0
....................................63
Figure 5.3. Scatter plots from the results of the Synth set, for S
?
...................................64
Figure 5.4. Scatter plots from the results of the Synth set, for ?
0
...................................65
xiii
Figure 5.5. Scatter plots from the results of the Synth set, for S
?
...................................66
Figure 5.6. Convergence plots from the results of the Synth set, for 0% noise ..............67
Figure 5.7. Convergence plots from the results of the Synth set, for 10% noise ............68
Figure 5.8. Convergence plots from the results of the Synth set, for 20% noise ............69
Figure 5.9. Convergence plots from the results of the Synth set, for 30% noise ............70
Figure 5.10. Comparison plots for some of the hysteresis loops of the Synth set;
0% noise........................................................................................................71
Figure 5.11. Comparison plots for some of the hysteresis loops of the Synth set;
10% noise......................................................................................................72
Figure 5.12. Comparison plots for some of the hysteresis loops of the Synth set;
20% noise......................................................................................................73
Figure 5.13. Comparison plots for some of the hysteresis loops of the Synth set;
30% noise......................................................................................................74
Figure 5.14. Sketch of a sandwich composite beam sample.............................................75
Figure 5.15. Parameters scatter plots from the optimization results of beam sample 1....80
Figure 5.16. Parameters scatter plots from the optimization results of beam sample 2....81
Figure 5.17. Parameters scatter plots from the optimization results of beam sample 3....82
Figure 5.18. Optimization convergence plots for beam sample 1.....................................83
Figure 5.19. Optimization convergence plots for beam sample 2.....................................84
Figure 5.20. Optimization convergence plots for sample beam 3.....................................85
Figure 5.21. Optimization comparison plots for some of the hysteresis loops
of sample 1....................................................................................................86
Figure 5.22. Optimization comparison plots for some of the hysteresis loops
of sample 2....................................................................................................87
Figure 5.23. Optimization comparison plots for some of the hysteresis loops
of sample 3....................................................................................................88
xiv
Figure 6.1. Overlay compliance FRF plots of all 3 beam samples utilized,
for low, medium, and high excitation amplitudes.........................................91
xv
LIST OF TABLES
Table 3.1. Parameters properties for decoding the binary chromosome of Fig. 3.2........29
Table 5.1. Mutated parameters of the best partial solution (of 4 parameters),
for the creation of an additional set of candidate solutions (L=3) for
local search.....................................................................................................55
Table 5.2. Excitation frequencies of the computer-generated Synth hysteresis set ........59
Table 5.3. System parameters used to generate the Synth hysteresis set ........................59
Table 5.4. Search space for the Synth hysteresis set .......................................................60
Table 5.5. Optimization results for the Synth data set; GA, 0% noise............................61
Table 5.6. Optimization results for the Synth data set; GA, 10% noise..........................61
Table 5.7. Optimization results for the Synth data set; GA, 20% noise..........................61
Table 5.8. Optimization results for the Synth data set; GA, 30% noise..........................61
Table 5.9. Optimization results for the Synth data set; PSO, 0% noise ..........................62
Table 5.10. Optimization results for the Synth data set; PSO, 10% noise ........................62
Table 5.11. Optimization results for the Synth data set; PSO, 20% noise ........................62
Table 5.12. Optimization results for the Synth data set; PSO, 30% noise ........................62
Table 5.13. Differences between the resultant averages and the actual
global optimum, for GA and PSO, at different noise levels...........................75
Table 5.14. Properties of the sandwich composite beams studied ....................................75
Table 5.15. Excitation frequencies used with sample 1 ....................................................76
Table 5.16. Excitation frequencies used with sample 2 ....................................................76
xvi
Table 5.17. Excitation frequencies used with sample 3 ....................................................77
Table 5.18. Excitation displacement amplitudes used with the beam samples studied ....77
Table 5.19. Search spaces for all 3 beam samples studied................................................77
Table 5.20. GA optimization results for beam sample 1...................................................78
Table 5.21. GA optimization results for beam sample 2...................................................78
Table 5.22. GA optimization results for beam sample 3...................................................78
Table 5.23. PSO optimization results for beam sample 1 .................................................79
Table 5.24. PSO optimization results for beam sample 2 .................................................79
Table 5.25. PSO optimization results for beam sample 3 .................................................79
1
CHAPTER 1
INTRODUCTION
1.1. Composite Materials and Damping
Composite materials are being increasingly used as alternatives to conventional
materials for highly demanding structural applications, such as the construction of
marine, automobile, and spacecraft structures, which include the equipment panels,
payload platforms, solar panels, and antenna reflectors
[1]
. Composites are mainly used
for these applications because they have the desired properties of high specific strength,
specific stiffness and tailorable design, along with a low mass. However, they can be
quite different from isotropic materials, because the former may have (for instance)
certain failure characteristics, such as core or matrix crazing, delamination, and interface
debonding, which typically do not occur in more conventional materials.
Damping is an important parameter related to the study of the dynamic behavior
of composite structures in general. There are numerous types of composites currently
available, but the present study deals with the use of a specific kind of composite
construction called ?sandwich composite?. A sandwich composite (or ?sandwich panel?)
is a structure consisting of two thin faces bonded to a thick lightweight core
[2-3]
. The
materials used to construct the faces are usually aluminum or a composite laminate, and
the core can be a lightweight foam or a honeycomb structure of some other material, such
as shown in Fig. 1.1.
2
Figure 1.1. Honeycomb sandwich composite; (a) detail of the core, (b) composite beam
(a)
(b)
3
Chandra et al
[4]
make an extensive review of damping mechanisms in composites, some
of which are detailed below in order to enhance the understanding of this important
physical phenomenon:
? Viscoelastic nature of matrix and/or fiber materials: the matrix shows the most
significant contribution to the damping. In the case of high damping fibers, such
as in carbon and kevlar composites, the effect of the fibers on damping should be
included in the analysis.
? Damping due to interfaces: the interfaces are the regions where adjacent surfaces
exist. In fiber-reinforced composites, the interfaces are the areas where the fibers
are in direct contact with the matrix. In sandwich composites, on the other hand,
the interfaces are the areas where the different layers are in contact with one
other.
? Damping due to damage: this may be caused by frictional damping due to slip in
the unbound regions between the fiber/matrix interface, and may also be caused
by delaminations, and damping due to energy dissipation in the area of matrix
cracks, broken fibers, etc.
There are also other mechanisms of energy dissipation in composites, including
viscoplastic and thermoelastic damping, which depend on the nonlinearities seen for high
amplitude vibration levels, and on the cyclic heat flow from regions of compressive to
tensile stresses, respectively.
4
Displacement / Strain
Forc
e / Stres
s
Displacement / Strain
Forc
e / Stres
s
Several damping measurement techniques have been developed, and the
applicability of each of these methods is normally dependent on the frequency range of
interest, and whether or not the damping in a resonance frequency is to be measured.
Among the several approaches available for damping measurement, the half power
method, logarithmic decrement, and hysteresis loops analysis are generally the most
popular. This topic will be discussed in more detail in section 2.2.
A mechanical hysteresis loop consists of a phase plot of the force vs. the
displacement usually at a specific point in a structure. Equivalently, it may be also
considered as a plot of stress vs. strain. For dynamical systems that include damping, the
displacement/force relationship is typically out of phase near resonance, generating an
elliptically-shaped hysteresis loop, for a linear system, or a sharp-edged loop, for a non-
linear system, as shown in Fig. 1.2. The area enclosed by the hysteresis loop is
proportional to the damping in the system.
Figure 1.2. Mechanical hysteresis loops examples; (a) linear system, (b) non-linear
system
(a) (b)
5
1.2. Scope of Study
The present work examines the hysteresis loop technique for the characterization
of the vibrational dynamics of three beam samples of sandwich composite materials,
consisting of a paper honeycomb core filled with foam and external layers of interlaced
strips of carbon fibers, such as was shown in Fig. 1.1b. This combination provides great
bending stiffness and endurance to a very light structure. These characteristics, along
with the great number of interfaces in a sandwich composite of this nature, typically
results in higher damping than standard homogeneous materials.
The main goal is to fit a parametric mathematical model to several experimental
hysteresis loops obtained for specific sandwich composite materials, simultaneously for
different displacement amplitudes and frequencies around the first resonance of the beam
sample. Experimental observations of the behavior of the materials are made, in order to
tailor a simple baseline mathematical model to another that includes such effects. Due to
the large scale of the system identification process, and the number of parameters to
identify, two popular heuristic optimization methods have been used in this problem:
Genetic Algorithms (GA), and Particle Swarm Optimization (PSO). Before using these
optimization techniques with actual experimental data obtained from the composite
samples with an unknown global optimum solution, computer-generated hysteresis loops,
with a known global optimum solution, and with similar vibrational behavior of actual
composites, are created. This is done to test the search performance of both optimization
techniques, with added Gaussian random noise of different levels.
6
CHAPTER 2
LITERATURE REVIEW
2.1. Hysteresis Models and Applications
The hysteresis phenomenon appears in several physical systems, such as
mechanical and magnetic materials and devices. In some cases, it can cause some
undesirable effects, such as loss of stability robustness; but it is mostly a way in how
energy is dissipated in a system, whose behavior typically characterizes damping. There
are several mathematical models of hysteresis, which are useful for different types of
systems. Their usefulness depends on the field of study to which they are applied, and
also depends on the type of materials to be characterized and the amount of nonlinearities
present. Sain
[5]
presents a survey on hysteresis models, focusing on four popular
representations, which have been used with success on general control designs. They are:
hysterons, the Chua-Stromsmoe, Preisach, and the Bouc-Wen model.
Hysterons were developed by Krasnosel?skii and Pokorovskii
[6]
, with the intent
of establishing a basic and general model of hysteresis. In essence, a hysteron represents
a transducer, which maps input functions into output functions. These output functions
are characterized by a family of curves, used to determine the output value of a specific
input, depending on its previous coordinate position and direction. In this way, the model
allows the state of the system to oscillate in a closed-loop manner, as pictured in Fig. 2.1.
7
The Chua-Stromsmoe hysteresis model
[7]
was created primarily with the intent of
modeling the behavior of ferromagnetic hysteretic inductors, with some suitable
characteristics, such as frequency-dependence control of the width of the loops. It is
given by the expression: [ ])()()()( xftugxhuwx ?= , with x(t
0
)=x
0
, where u is the input,
and x is the output. The functions named f and g represent dissipative and restorative
phenomena, respectively; h controls the width of the loop, and w controls the frequency-
dependent widening behavior. A plot describing this frequency widening effect is shown
in Fig. 2.2, for w( u )=h(x)=1, f(x)=100x
3
, g(u)=u
5
, and u(t)=sin(?t).
The Preisach model of hysteresis
[8]
is based on a weighted summation of the
outputs of a set of hysteretic relays, each of them integrated over a set of individual
switching thresholds, a and b. This means that, when a threshold is attained for a single
relay, the output changes between ?
1
and ?
0
, (which are the high and low levels,
respectively), and remains at a specific level until the next threshold is attained, where it
switches back to the initial level. The basic behavior of such relay is shown in Fig. 2.3. A
first discrete approximation of the behavior of the Preisach model can be seen on Fig. 2.4,
where the weighted summation of the outputs of three parallel relays generates a very
basic hysteresis loop. In the limit, where an infinite number of relays are working in
parallel, a continuous response can be attained, and it can be interpreted into an intuitive
geometrical representation. If we consider an a-b plane, and for a monotonically
increasing input value, the output of the Preisach model corresponds to the area obtained
from the bottom left corner of a minimum value over an a=b line up to the current value
of the input (Fig. 2.5a). Then, for a monotonically decreasing input, it vertically wipes
this area , from the maximum value reached up to the current input value (Fig. 2.5b).
8
Output
Increasing
frequency
Input
Figure 2.1. Plot of a hysteron of the second kind
Figure 2.2. Chua-Stromsmoe hysteresis loop: Effect of increasing the loop area,
for increasing excitation frequencies
Input
Output
9
Figure 2.3. Input-output behavior of a hysteretic relay
Figure 2.4. Weighted summation of 3 hysteretic relays
Input
Output
?
1
?
0
b a
b
1
a
1
b
2
a
2
b
3
a
3
?
Input
Output
b
1
b
3
b
2
a
3 a2 a1
10
Figure 2.5. Geometric interpretation of the Preisach model;
(a) monotonically increasing input, (b) monotonically decreasing input
Finally, the Bouc/Wen model
[9-10]
is one of the most popular mathematical
representations of several hysteresis phenomena, but it is mostly used in vibrations
problems, as it is based on the differential equation of motion of a spring-mass-viscous
damper mechanical system. The difference from the standard linear model, is that it adds
an extra term in the equation, referred as the hysteretic nonlinear term. This model can be
expressed as a system of two coupled nonlinear differential equations, defined by
Eq. (2.1).
a
b
a=b line
monotonically
increasing
This shaded area is the
value of the output
?
Input
Output
(a)
a
b
a=b line
monotonically
decreasing
?
Input
Output
This shaded
area is the
value of the
output
(b)
11
( )
()
22
1
a) 2 1 ( ) ,
b) ,
nn n
nn
x xx zft
zxzzxzax
?? ?? ? ?
??
?
+++?=
=? ? +
(2.1)
where the parameters of the system are:
? : rigidity ratio (0 ? ? ? 1),
? : linear elastic viscous damping ratio (0 ? ? ? 1),
?
n
: pseudo-natural frequency of the system (rad/s),
a : parameter controlling hysteresis amplitude,
?,?,n : parameters controlling hysteresis shape (?>0, n?1) .
The f(t) term in Eq. (2.1a) is a mass normalized forcing function (N/kg). The
hysteretic (or memory) variable z is a ?fictitious? displacement related to the actual
displacement, x. This variable accounts for the dependence of the nonlinear restoring
force (last term on the left-hand side of Eq. (2.1a)) on not only the current displacement
value, but on the previous time history of the motion, as well. Different combinations of
the parameters yield different shapes for the hysteresis loops. To illustrate this let n=1
and vary only the parameters ? and ?. Different possible stable hysteresis loops
[11]
are
shown in Fig. 2.6. The versatility of the system of differential equations described in
Eq. (2.1) can be observed. Numerous researchers have successfully used this model. For
example, Ni
[12]
and Constantinou
[13]
used it for the study of nonlinear hysteretic
isolators; Smyth
[14]
for an adaptive application and Heine
[15]
for an optimization
approach to degrading hysteretic joints with slack behavior.
12
Figure 2.6. Hysteresis loops from Bouc-Wen model, obtained by varying ? and ? ; n=1;
(a) ?+?>0 , ???<0 , (b) ?+?>0 , ???=0 , (c) ?+?>???>0 ,
(d) ?+?=0 , ???<0 , (e) ???< ?+?<0
2.2. Damping, and Measurement Techniques
Extensive research has been performed on the topic of energy dissipation of
materials, whose damping mechanisms and modeling depend on specific types of
materials. For metals, viscoelastic materials, etc. Macioce
[16]
gives a brief but precise
explanation of damping and some of its measurement techniques, as follows.
A material or structure subjected to vibrations has a combination of potential and
kinetic energy, along with a dissipative element, which converts some of this mechanical
energy into heat during each cycle of motion, which is then lost when it is transferred to
the surrounding environment. The amount of dissipated energy is related to the damping
of the structure, which is a measure of the ability of the system to dissipate mechanical
energy into thermal energy. Several forms of damping have been identified, like
F
x
F
x
F
x
F
x
F
x
(a) (b) (c)
(d)
(e)
13
Coulomb friction, fluid viscosity, air damping, particle damping, magnetic hysteresis, and
acoustic radiation damping. As it can be seen, the environment can also play a role in
providing a mechanism for added external damping, but it is insignificant in most cases.
The most common models of damping used in solving vibrations problems, are the
viscoelastic damping (also called ?structural?, ?material?, or ?hysteretic? damping),
viscous damping, and Coulomb damping.
Viscoelastic damping is found in a wide range of materials (known as
?viscoelastic materials?), such as rubbers, epoxies, and foams. One main characteristic is
that it is rate-independent, which means that the energy stored per cycle remains constant
when changing the excitation frequency. And another important characteristic is that its
modulus of elasticity (a function of the temperature T and frequency f, in general) is
expressed as a complex term, which is a sum of a conservative part (E
1
) and a dissipative
part (E
2
), as shown in Eq. (2.2).
( )
121
(, ) 1 ,ET f E iE E i?=+ = + (2.2)
where ?=loss factor (a descriptor of damping). So, if the behavior of this type of material
is included in a 1-DOF (Degree of Freedom) system, its equation of motion is represented
in Eq. (2.3).
( )() 1 () () .mx t k i x t F t?+?+ = (2.3)
On the other hand, the viscous damping model is rate-dependent, and it may
generally be used for some non-viscoelastic materials. A 1-DOF equation of motion
including this type of damping is shown in Eq. (2.4), where ? =?/2, and ?
n
=(k/m)
1/2
.
2
() () () () () 2 () () () .
nn
mxtcxtkxtFt xt xt xtFt?? ?++= ? + + = (2.4)
14
Finally, Coulomb damping (also called ?dry damping?)
[17-18]
is produced from the
sliding friction of two adjacent dry surfaces. The dry friction force that appears is parallel
to the surface, and its magnitude is proportional to the normal force on the body (such as
its own weight), but in a direction opposite to that of the velocity. A 1-DOF equation of
motion with this type of damping is shown in Eq. (2.5), with F
d
=??N=magnitude of the
damping force, ?=coefficient of friction, N=force normal to the sliding surface, and
sgn( )x xx==sign of the velocity.
() sgn( ) () () .
d
mx t F x kx t F t+ += (2.5)
Some other models of damping have been developed, such as Zener?s damping
model for metals
[19]
. As explained by Crandall
[20]
, this model is frequency dependent,
and is based on the assumption of the existence of heat flow from warmed compression
fibers in a material to cooled tension fibers. The mathematical expression for the loss
factor, ?, for this model is shown in Eq. (2.6)
[20-21]
. A plot of the modeled loss factor for
a 2024-T4 aluminum alloy cantilever beam is shown in Fig. 2.7.
22
,
1
s
r
s
??
?
??
=?
+
(2.6)
where:
r
? =
2
TT
p
ET
C
?
?
= relaxation strength,
?
s
=
2
p
C
d
K
?
?
??
??
??
= thermal relaxation time,
? = frequency of vibration,
?
T
= temperature expansion coefficient,
E
T
= Young?s modulus at constant temperature,
15
T = absolute temperature,
?C
p
= specific heat per unit volume,
d = material sample thickness,
K = thermal conductivity.
0.1 1 10
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
x 10
-3
Normalized Frequency (f / f
0
)
L
o
ss F
a
ct
o
r
(
?
)
Figure 2.7. Loss factor (?) as a function of frequency, of a 2024-T4 aluminum alloy
cantilever beam, obtained with the Zener?s damping model
In several applications, the damping is usually very light, and its effect is
observable only around resonances. Also, in many practical situations, it is common to
represent the damping of a system with an equivalent viscous damper (dashpot), which
simplifies the mathematical calculations, and provides a sufficiently accurate
approximation to the total damping. In order to experimentally quantify the amount of
damping present in a material or structure, there are several measurement techniques
developed by researchers in the area of vibration analysis. Three of the most popular
16
techniques used are
[16,17,21]
: the half power bandwidth method, logarithmic decrement,
and hysteresis loops analysis.
The half power bandwidth method consists of determining the modal damping
from a classic frequency response function (FRF), such as the compliance response
(displacement/force), also known as admittance, or receptance. The amount of damping
for a specific mode of vibration will depend on the sharpness of the resonance peak, i.e.
the rounder the peak, the larger the damping. From the compliance function shown in
Fig. 2.8, the damping ratio ? can be determined with Eq. (2.7).
21
0
,
22
? ??
?
?
?
==
?
(2.7)
where ?
0
=frequency of resonance peak, ?
1
,?
2
=lower and upper frequencies of -3dB from
the resonance peak. This method is also often referred as the ?3dB method?.
Figure 2.8. Compliance FRF, and half power bandwidth method
Com
p
li
ance
[
dB
]
Frequency
H
0[dB]
H
-3[dB]
?
0
?
1
?
2
()
X
H
F
? =
17
The logarithmic decrement method consists of computing the damping of an
unforced underdamped system from the rate of decay of its response in the time domain.
The amount of amplitude that has decreased in the system between two consecutive
cycles (at times t
1
and t
2
, respectively) is shown in Eq. (2.8)
[17]
, and since
t
2
=t
1
+T=t
1
+2?/?
d
, the expression reduces to Eq. (2.9).
( )
()
1
2
1
1
22
cos
,
cos
n
n
t
d
t
d
Ae t
x
xAe t
??
??
??
??
?
?
?
=
?
(2.8)
1
1
1
()
2
.
n
n
n
t
T
tT
x e
e
xAe
??
??
??
?
?+
== (2.9)
Because of the exponential nature of Eq. (2.9), it is proper to introduce the notation:
1
2
2
2
ln ,
1
n
x
T
x
??
???
?
===
?
(2.10)
where ? is the logarithmic decrement. Now, solving for ?, the damping ratio (in terms of
? ) is given by Eq. (2.11).
1
22
2
1
ln .
22
(2 )
x
x
??
?
??
??
??
=?=
??
+ ??
(2.11)
The approximation given in Eq. (2.11) is true for small values of damping, which is often
encountered in most systems. It is also possible to determine the damping ratio by
measuring the logarithmic decrement of two displacements separated by N cycles,
utilizing Eq. (2.12).
11
ln ln .
2
NN
xx
Nx N x
??
?
++
?? ??
=??
?? ??
?
?? ??
(2.12)
18
The hysteresis loop method of damping analysis consists in calculating the energy
loss per cycle at steady-state harmonic loading. This may be achieved by creating a plot
of displacement vs. force (or equivalently, strain vs. stress) of the system at a fixed
excitation frequency, and determining the area enclosed by this hysteresis loop, which is
precisely the energy dissipated per cycle, and is proportional to the damping in the
system. This type of plot is equivalent to the Lissajous figures of two signals being (in
general) out of phase. Below resonance (and far from it), the hysteresis plot usually
shows a straight line with positive slope. This means that both, the displacement and the
force waves, are in phase. When the excitation frequency is close to resonance and below
it, both signals are slightly out of phase (displacement lags the force), and the hysteresis
loop is a tilted ellipse with positive slope (for linear systems, as in the viscous and
viscoelastic damping models), as shown in Fig. 1.2a, or a closed but not-elliptical loop
for the nonlinear case, as shown in and example in Fig. 1.2b. At resonance, both signals
are exactly 90? out of phase, and the hysteresis is a non-tilted closed loop, or in some
cases, a circle, for the linear case. Above resonance and close to it, the hysteresis loop is a
closed loop with negative slope. Finally, above resonance and far from it, the hysteresis
turns again into a straight line, but with negative slope, which means that the
displacement lags the force in 180?. A summary of all this behavior is shown in Fig. 2.9.
An equivalent methodology of damping analysis using this technique consists of
performing a curve-fitting of the mathematical model of motion of the system to the
actual loops obtained experimentally, this way finding a set of parameters from the model
that is able to represent the actual dynamics of the system in reality.
19
Figure 2.9. General behavior of hysteresis in a system, when crossing resonance (f
0
);
hysteresis loops, along with time trace overlays of displacement (x) and force (F)
x
F
x
F
x
F
x
F
x
F
0
f f
0
f f<
0
f f=
0
f f>
0
f f
F
x
t
F
x
t
F
x
t
F
x
t
F
x
t
20
A variety of novel ideas for damping measurement have been introduced by
researchers, who have also described the critical factors for specific problems.
Hooker
[22]
, for example, explains how damping measurements for biaxial hysteresis
loops yield apparent loss factors, which may be related to the true loss factors of the
system. On the other hand, So, et al
[23]
shows a method for frequency-response-based
damping measurement using a periodic chirp excitation, similar to pseudo-random noise
excitation, but with a much higher peak-to-rms ratio. Also, McDaniel, Dupont, and
Salvino
[24]
introduce an approach for the determination of a frequency-dependent
damping from transient loadings.
2.3. Experimental Considerations
The common methodology used as a first approach in the development of
mathematical models is assuming that the system is ?perfect?, in the sense that its
behavior is supposed to be completely represented by the proposed equations. In reality,
this is almost impossible to achieve, as several additional factors related to the system,
such as effects from the environment, material damage, or attachments, may have not
been considered. This may alter the accuracy of the model predictions. Hence,
experimental considerations must be made in order to incorporate these effects into the
equations, with the purpose of getting as close as possible to the exact solutions. Several
researchers have addressed these practical issues as a way to improving the predictions of
their models. Some of them are presented next, which are relevant to the present work.
Renji and Shankar Narayan
[1]
, along with Brown and Norton
[25]
, explained the
importance of considering the effect of mass of the force transducer (or impedance head)
21
in their experiments. They argue that the measured force by the transducer will be
actually different from the actual applied force. Its mass causes a variation in the
impedance at the attachment point. Thus, the driving point admittance Y, as well as the
transducer admittance Y
m
, must be considered in order to obtain the actual force f
a
from
the measured force f
m
, as shown in Eq. (2.13).
()
.
1/
m
a
m
f
f
YY
=
+
(2.13)
The admittance introduced by the mass of the transducer is simply 1/( )
m
YjM?= , and the
driving point admittance is Y=?
fv
/?
ff
, where ?
ff
is the auto-spectral density of the force,
and ?
fv
is the cross-spectral density between the force and the velocity. The direct use of
the measured force (without corrections) in models may lead to errors in the modal
damping values calculated, although it has been seen that this effect is more significant at
high frequencies rather than low frequencies, where in the latter the effect may be
neglected without significant loss of accuracy.
On the other hand, Boyaci
[26]
explains the effects of non-ideal boundary
conditions as follows. Beam samples are commonly used for the study of physical
properties of materials, and their supports may have a direct effect on the experimental
results. Different types of boundary conditions are defined, such as clamped, free, simply
supported, etc. However, there may be small deviations from these ideal conditions, such
as small gaps and friction in the attachments, which may introduce small deflections and
moments at the ends. These non-ideal boundary conditions are modeled using
perturbation theory. More specifically, the method of multiple scales is used. It was found
22
that while the stretching effect increases the frequencies, the non-ideal boundary
conditions may increase or decrease them.
Finally, Qiao, Pilpchuk, and Ibrahim
[27]
made a similar analysis, focusing on
parameter uncertainties and relaxation of joints. They suggest that relaxation effects in
the attachments at the boundaries will cause time-dependent boundary conditions, and
will depend on the level of structural vibration. The main problems that are encountered
in those cases are the existence of random eigenvalues, and probability of failure. The
analysis of stiffness uncertainties at the joints under dynamic excitation were performed,
with and without the effect of the relaxation process. It has been found that the frequency
of resonance of a material is reduced when the ?level of clamping? is reduced. Also, the
relaxation of the joints produces a gradual decrease in the resonance frequency.
2.4. Heuristic Optimization Algorithms in the Study of Composites
The design of materials and system identification of model parameters provide an
optimization problem which happens to be difficult in most cases. First because of the
functional forms of the optimization descriptors often encountered in engineering
problems, and secondly, because of large dimensionalities of the search spaces to be
analyzed. Several computational tools are available for optimization, which in turn may
be categorized as shown in Fig. 2.10
[15]
. As can be seen, heuristic optimization
corresponds to the continuous and global search methods, which makes this category
very powerful and versatile as compared to other techniques. Their most important
characteristics are: (1) they don?t need a starting search point (initial guess), as they
generate populations of candidate solutions, usually spanning a great area of the search
23
space simultaneously, and (2) they work in a random fashion (i.e. not using calculus
operations in the search), looking for areas of interest in the search space, and gradually
converging to optimum solutions within these areas. Because of this random nature of the
methodology, different runs of the optimization computer routines yield also different
results, but generally close to each other.
Figure 2.10. Categories of optimization techniques
Heuristic optimization methods have been extensively used for solving problems
in different areas, such as in mechanics, electronics, economy, and several others.
Focusing on the topics relevant to this work, we have found publications on the
utilization of Genetic Algorithms (GA) and Particle Swarm Optimization (PSO) in
Optimization
Methods
Discrete Continuous
Local
Global
Calculus-based
Indirect Direct
Heuristic Enumerative Exact
Possible
combination
Genetic Algorithms Particle Swarm
Optimization
24
composite materials, for design optimization, extraction of physical properties, and
system identification of model parameters.
Liu, Han, and Lam
[28]
used a combination of GA with a nonlinear least squares
method (LSM) in a problem of material characterization, where the elastic constants,
(Young?s moduli, shear modulus, and Poisson?s ratios) of some composite plates were
required to be determined. They used GA to find a solution close to the optimum, which
is then used as an initial guess for LSM to finally converge to that point. Matsuoka,
Yamamoto, and Takahara
[29]
used a combination of GA and finite element stress analysis
on the search of material structures that give way to desired elastic properties. Muc and
Gurba
[30]
used a combination of GA and finite element analysis for layout optimization
of composite structures, which consists in finding information on stacking sequences,
shape, and sizes of the material layers, in order to achieve optimum mechanical
properties. Gantovnik, G?rdal, and Watson
[31]
also used GA for the optimal design of
laminated sandwich composite panels, in which memory is included for continuous
variables in order to preserve data from previous analyzed designs, thus avoiding the
repetition of an analysis already performed, and then reducing the computational costs.
Regarding the use of PSO in composites design, Kovacs, Groenwold, Jarmai, and
Farkas
[32]
have used this technique for the optimal design of a carbon-fiber-reinforced
plastic (CFRP) sandwich-like structure with aluminum webs, in order to achieve a
minimum cost and maximum stiffness. The CFRP plate is optimized with respect to ply
arrangement, while the complete structure is optimized with respect to a combination of
manufacturing and material cost. Sekishiro, and Todoroki
[33]
also used PSO for the
optimal design of laminated composite structures, focused in achieving a minimum
25
weight of a hat-stiffened composite panel subjected to a buckling constraint. In order to
reach their goal, they used the fractal branch and bound method (a stacking sequence
optimization method) integrated into PSO.
In the topic of heuristic optimization methods applied to mechanical hysteresis,
the publications available are focused more on GA, rather than PSO, but in neither case
focused specifically on composites. Heine
[15]
, however, studied the response of
degrading hysteretic joints with slack behavior in connected timber members, by using
GA for the identification of the parameters of a modified Bouc-Wen hysteresis model.
Kyprianou and Worden
[34]
carried out a generalized study of the performance of GA on
system identification of the Bouc-Wen hysteresis model. Zhang, et al
[35]
, also presented a
system identification problem on the parameters of the Bouc-Wen hysteresis model by
utilizing GA, in which the crossover and mutation operators are adopted in between of
the gene code strings of the same parameters, improving the convergence of the
algorithm. On the other hand, Peng, Li, and Suzuki
[36]
successfully used PSO in the
identification of the parameters of a nonlinear differential hysteresis model with pinching,
based on the Duhem hysteresis operator.
In general, both GA and PSO have been proven to perform well in most
engineering problems, but several researchers
[37-38]
have reported that PSO has been
much more efficient than GA in terms of both, quality of the solutions, and simplicity of
implementation.
26
CHAPTER 3
HEURISTIC OPTIMIZATION METHODS
3.1. Genetic Algorithms (GA)
As explained in section 2.4, characterizing the parametric values of a
mathematical model is some times a challenging problem due to the relative complexity
of the functional forms of the models often used, especially in nonlinear problems. A
possible solution is the application of an optimization strategy to determine the specific
parametric values that best fit the observed behavior. In general, a parameter estimation
method is an iterative process that takes a known input u(t) of a system, evaluates it into a
mathematical model with a starting set of parameter values, compares its result with the
trial data (thus obtaining a value of error between them), and then modifies the values of
the model parameters, expecting the new set to get closer to the optimum solution. This
process is summarized in Fig. 3.1
[15]
.
Figure 3.1. Schematic of a parameter estimation method
Model Experiment
Parameter
adjustment
?
()f t
Error
()ut
()F t
27
The complexity and high dimensionality of some models lead to the use of a
heuristic search method. For those problems, Genetic Algorithms (GA) has proven to be
a suitable optimization tool for a wide selection of problems.
This method was developed by John Holland
[39]
over the course of the 1960?s and
1970?s and finally popularized by one of his students, David Goldberg, who was able to
solve a difficult problem involving the control of a gas-pipeline transmission for his
dissertation
[40]
. Some of the advantages of a GA include that it
[41]
:
? Optimizes with continuous or discrete parameters,
? Doesn?t require derivative information,
? Simultaneously searches from a wide sampling of the cost surface,
? Deals with a large number of parameters,
? Is well suited for parallel computers,
? Optimizes parameters with extremely difficult cost surfaces; it can jump out of a local
minimum,
? Provides a list of optimum parameters, not just a single solution,
? May encode the parameters so that the optimization is done with the encoded
parameters, and
? Works with numerically generated data, experimental data, or analytical functions.
GAs may be also used to identify an initial guess for calculus-based optimization
techniques, in order to locally refine the search.
28
GAs have been a breakthrough to optimization theory because of their versatility,
performance, and ease of use. The GA method is based on the ideas associated with the
natural evolution of species, in which the ?fittest? organisms in a population survive,
passing their genes to the next generation, and so forth. Basically, the numerical version
of this mechanism is a population composed of a large set of chromosomes (set of
parameters), which generate results from a parametric mathematical model. These results
are then compared to an objective function, which is usually the Mean Square Error
(MSE) function, in order to determine the fitness of a particular chromosome. This yields
a set of fitness values, which are used to sort the population, eliminate the badly-fitted
individuals, mate the best-fitted chromosomes, and then propagate the ?good? genes
(parameter combinations) to the upcoming generation. This process may be repeated a
fixed number of generations or indefinitely, until a certain error threshold is attained.
There are two versions of GA: one with binary encoding (all parameters are encoded by
long strings of ones and zeros), and another one with floating-point-number parameters
(?Differential Evolution Algorithm?, or ?Real-Coded GA?). An example of these two
versions is shown in Fig. 3.2, for a search space with 4 parameters.
Figure 3.2. The two versions of GA chromosome representation; upper, binary
representation, with 5 bits for each gene (parameter); lower, real-coded representation
0 0 1 1 0 0 0 0 1 1 1 1 0 1 0 0 1 1 0 0
Parameter #1 Parameter #2 Parameter #3 Parameter #4
0.1935 -26.667 8387.097 7.125
CHROMOSOME
GENE
ALLELE
29
In the binary-coded representation, a population of several random candidate
solutions (chromosomes) is first generated, where all parameters are encoded into strings
of ones and zeros. The algorithm works in this ?encoded environment? (genotype space)
while modifying the population. Then, it must decode each parameter in each
chromosome into real numbers (phenotype space) in order to evaluate them with a
suitable fitness function. This decoding may be carried out by using Eq. (3.1)
[42]
.
()()
(,,,) ,
(2 1)
i
ii i
iiiii iL
ub lb d P
decode ub lb L P lb
? ?
=+
?
(3.1)
where:
ub
i
= upper bound of parameter i,
lb
i
= lower bound of parameter i,
P
i
= binary string of parameter i,
d(P
i
) = decimal conversion of binary string P
i
,
L
i
= length of binary string of parameter i.
In Fig. 3.2, the values presented in the example of a real-coded GA, correspond to
the decoded values of each binary string in the binary GA. The decoding was performed
assuming the values given by Table 3.1.
Property Parameter #1 Parameter #2 Parameter #3 Parameter #4
Upper bound (ub) 1 40 10000 20
Lower bound (lb) 0 -40 0 1
String length (L) 5 5 5 5
Table 3.1. Parameters properties for decoding the binary chromosome of Fig. 3.2
30
The real-coded representation works in a similar way as in the binary version, but
no conversions are done. Instead, it works just in the phenotype space, so this method is
usually much easier to implement.
Both versions perform the same genetic operations, but with different
mechanisms. The operations are, in order of occurrence: (1) initial population generation
(generates randomly a set of chromosomes), (2) fitness evaluation (determines how fit is
each of the chromosomes to the objective function), (3) selection (determines pairs of
chromosomes for mating/combination, either depending statistically on their fitness
values, ranking, or tournament results), (4) crossover (for each pair of parents, it sets a
parameter to be combined), (5) mating (combines the parameters determined by the
crossover stage and swaps a certain number of parameters between each of the selected
parents), and (6) mutation (randomly changes the value of a certain amount of parameters
in the new offspring; this mutation version is called ?uniform mutation?). Usually the
best chromosome is allowed to propagate to the next generation unchanged by mutation.
Such an approach is called an ?elitist strategy?. Fig. 3.3 summarizes the main operations
of the algorithm, where Fig. 3.4 shows it in pseudo-code. The aforementioned steps are
explained next in more detail.
The initial population is generated first by defining the search space (i.e. the
bounds ub and lb for each parameter j), and then using Eq. (3.2).
( )
() rand(0,1) ,
ij j j j
p tublb lb=?? + (3.2)
where:
p
ij
(t) = parameter j of chromosome i at itereation t,
rand(0,1) = uniform random number between 0 and 1.
31
Mutate
Define: parameters,
bounds, cost function
Generate Initial
Population
Evaluate costs
(fitness)
Select
mates
Assign crossover
points
Mate
parents
STOP
YES
START
Evaluate costs
(fitness)
MAIN LOOP
?converged?
NO
Figure 3.3. Main steps for a typical Genetic Algorithm (GA)
Figure 3.4. Pseudo-code of a typical Genetic Algorithm (GA)
Procedure GA {
t=0;
Create_Population(P(t));
Evaluate_Fitness(P(t));
while (not done) {
t=t+1;
M(t)=Select_Mates(P(t));
Assign_Xover(M(t));
Offspring(t)=Mate(M(t));
P(t)=Add_to_Population(Offspring(t));
Mutate(P(t));
Evaluate_Fitness(P(t));
Select_Survivors(P(t));}
}
32
The fitness of a candidate solution (chromosome) is evaluated by a suitable error
mathematical expression, such as the Ordinary Least Squares (OLS) error function, given
by Eq. (3.3), with the goal of finding a chromosome with the lowest possible value of
OLS.
()()
2
12
1
?
OLS ( ) , ,..., ; ( ) ( ) ,
N
ik nk
k
tfppputFut
=
? ?
=?
? ?
?
(3.3)
where:
OLS
i
(t) = ordinary least square error of chromosome i at iteration t,
?
k
f = model output data point number k,
F
k
= experimental output data point number k,
N = total number of data points (samples) to be evaluated.
The parents selection may be carried out by several ways. The most popular
techniques are: fitness-based, ranking-based, and tournament selection. In the fitness-
based selection, a probability P
i
of being chosen as a parent for recombination is
calculated for each chromosome (at iteration t) depending on their fitness value, by using
Eq. (3.4).
1
()
() ,
()
i
i
M
q
q
Ct
Pt
Ct
=
=
?
(3.4)
where C
i
(t)=[Fitness
i
(t)-Fitness
M+1
(t)], Fitness
i
(t) is the fitness of chromosome i at
iteration t, and M is the total number of chromosomes to keep from the population, for the
general case of discarding a specific amount of ?badly-fitted? chromosomes at each
iteration. In the ranking-based selection, however, the probability of being chosen is
33
calculated from the ranking position, where the chromosome with the best fitness (lowest
error) is ranked #1, and so forth. The probability in this case is given by Eq. (3.5).
( )
()
1
21
1
() .
1
i M
r
Mi
Mi
Pt
MM
r
=
??+
?+
==
+
?
(3.5)
In both, fitness-based, and ranking-based selection, a technique called ?roulette wheel? is
used in order to pick pairs of parents for recombination, which consists in an analogy of
spinning a roulette wheel which has M sections, each of them with a size corresponding
to the probability determined by either Eq. (3.4) or (3.5), as shown in Fig. 3.5 for a
population where 4 individuals are to be kept. After spinning the wheel a first time, an
individual will be randomly picked as the first parent. Similarly, the second parent is then
picked after spinning the wheel a second time.
Figure 3.5. Roulette wheel analogy, as a parent selection strategy
P
1
P
2
P
3
P
4
34
The tournament selection is the easiest way of obtaining couples from the population. It
consists of simply selecting a random subset (of q individuals, where q is known as
?tournament size?) from the mating pool, and the best one (with the lowest value of error)
is selected as the first parent. Then, the process is repeated in order to obtain the second
parent.
The crossover operator is utilized to determine swapping points within the
chromosomes of the selected couples, with the purpose of using this information for the
mating process. Particularly, in the real-coded chromosome representation, the crossover
point may be selected between genes, or it may be over a gene, as shown in Fig. 3.6.
Notice that multiple crossover points may be selected for each couple, although single
point crossover is often employed.
Figure 3.6. Two types of crossover points; (a) between genes, (b) over a gene
The mating process is used to recombine each couple in order to generate
offsprings. In this step, the genes at one side of the crossover point selected for a couple
are swapped, thus creating a new pair of chromosomes. This process is shown in Fig. 3.7.
D
1
D
2
D
3
D
4
D
5
D
6
M
1
M
2
M
3
M
4
M
5
M
6
DAD:
MOM:
CROSSOVER POINT
D
1
D
2
D
3
D
4
D
5
D
6
M
1
M
2
M
3
M
4
M
5
M
6
DAD:
MOM:
CROSSOVER POINT
(a)
(b)
35
Figure 3.7. Mating process, when the crossover point is between genes
In the case where the crossover points are implemented to lie over the i-th gene in a
couple, the process is similar as in Fig. 3.7 for all other genes. The crossover genes are
then blended using the generalized Michalewicz
[43]
crossover method (Eq. (3.6)), which
performs a weighted average on them, having the capacity of extending the new values
out of the range defined by both genes at the crossover point in a proportion defined by
the optimization parameter ?. This type of crossover method is known as BLX-?. The
procedure is shown in Fig. 3.8.
( )
()
1
2
,
,
new i i i
new i i i
pD DM
pM DM
??
??
=?? ?
=+? ?
(3.6)
where ?=rand(0,1), a uniform random number between 0 and 1, in which the same value
must be used for both expressions.
Figure 3.8. Mating process, when the crossover point lies over a gene
D
1
D
2
D
3
D
4
D
5
D
6
M
1
M
2
M
3
M
4
M
5
M
6
OFFSPRING 1:
OFFSPRING 2:
CROSSOVER POINT
D
1
D
2
p
new2
D
4
D
5
D
6
M
1
M
2
p
new1
M
4
M
5
M
6
OFFSPRING 1:
OFFSPRING 2:
CROSSOVER POINT
36
After the mating process, mutation is applied to some genes in the whole new
population, in order to avoid letting the algorithm converge to suboptimal solutions, thus
being able to jump out of these local minima and keep searching for the global optimum.
It may be either applied just to genes of the newly created offspring, or to both, the
population of parent chromosomes and offspring. There can be mutation of a specific
percentage of genes of the population, or a mutation probability may be set in order to
determine, for each gene, if mutation is going to be performed. The value to be set for a
specific mutated gene can be either a uniform random number within the bounds defined
for it (uniform mutation), or a number around the current value, at a distance defined by a
mutation amount factor, the parameter bounds, and a Gaussian normally distributed
random number (Gaussian mutation).
In the final steps of GA, the fitnesses of the new chromosomes in the population
(offspring and chromosomes with mutated genes) are calculated, followed by a process of
survival, in order to keep a certain amount of the best performing chromosomes. The
survival can be either by using a (?,?) strategy (where ? newly-generated chromosomes
substitute the ? worst performing ones from the parent population of size ?), or a (?+?)
strategy (where the survivors are the ? best performing chromosomes from a pool
composed by the parent population of size ?, and the ? newly-generated ones). The
process is then repeated (from the selection step) until a certain amount of iterations or a
fitness threshold is reached.
37
3.2. Particle Swarm Optimization (PSO)
The Particle Swarm Optimization (PSO) is an optimization technique for
continuous nonlinear functions. It was originally created for the simulation of a
simplified social model, such as bird flocking or fish schooling, where a set of individuals
(particles) interact in order to find an optimal solution for either maximization or
minimization problems. It was first introduced by James Kennedy and Russ Everhart in
1995
[44]
. However, there are many variations of this approach, each created by different
researchers in order to improve the search performance for specific problems. Shi and
Eberhart
[45]
have shown that PSO searches wide areas effectively, but usually lacks
precision for local search. Their suggestion for solving this issue is to introduce an inertia
factor (?) in the standard equation, that adjusts the velocity correction over time,
gradually concentrating the algorithm into a local search, as shown in Eq. (3.7).
( ) ( )
12
,a) rand(0,1) rand(0,1)
b) ' ,
id id id id gd id
id id id
vv px px
xxv
??+=?? + ? ? ? ? ? ?
=+
(3.7)
where i is the particle index, d is the dimension index (parameter), g stands for ?global?,
v is the velocity (correction), ?
1
and ?
2
determine the relative influence of the ?cognition?
and ?social? components of the velocity, p is the location of the previous best position, x
is the current particle position, and rand(0,1) is a uniform random number between 0
and 1. Maurice Clerc
[46]
introduced a so-called ?constriction factor? (K), that improves
the control of the velocities, given by Eq. (3.8).
,
2
2
24
K
?? ?
=
?? ?
(3.8)
where ?=?
1
+?
2
; ?>4. The velocity correction equation for the PSO is given by Eq. (3.9).
38
()
()()
12
.rand(0,1) rand(0,1)
id id id gd idid
vKv px px??=? +? ? ? +? ? ? (3.9)
From Eqs. (3.7) to (3.9), it can be seen that PSO moves each particle of the swarm
(population) in a direction determined by the current particle?s best and the global best,
each weighted by the cognition and social constants (?
1
and ?
2
, respectively), as can be
seen in Fig. 3.9.
Figure 3.9. PSO method of particle position correction at j-th iteration
The main advantage of PSO is that it is simple to implement, as compared to other
evolutionary computation optimization techniques and generally has an outstanding
performance. For example, for some popular testing functions
[47]
, like the Sphere (De
Jong?s F1), Rosenbrock, generalized Rastrigin, generalized Griewank, and Schaffer?s F6
functions, success rates of 100% and low numbers of function evaluations may be easily
obtained.
x
i
(j)
?
2
?
1
SEARCH SPACE
SOCIAL COMP
ONEN
T
CO
GNIT
I
O
N COMPONENT
v
i
(j)
x
i
(j+1)
P
i
(j)
P
g
(j)
39
PSO?s can have different neighborhood topologies. A particle?s neighborhood is
defined as a set of particles that influence an individual particle movement through the
search space. It defines what is considered as P
g
(global best position) for a specific
particle being evaluated. The larger the neighborhood size, the more particles included in
a particular particle?s neighborhood
[47]
. There are two basic topologies for PSO: the ring
topology and the star (or global) topology. In the ring topology, the particles are thought
to be held in a wrap-around array, with each particle?s neighborhood consisting of
themselves, the particle to their left, and the particle to their right. The number of
particles to include in the neighborhood set can expand or contract, adding the particles
that are two positions away, three positions away, etc., thus increasing the neighborhood
size. If the neighborhood set includes every single particle in the swarm, then the
neighborhood is ?fully-connected?. A fully-connected neighborhood is said to be in a star
(or global) topology. Fig. 3.10 shows these two types of topology for a swarm of N
particles, where p
1
to p
N
are the individual particles, and p
i
is the current particle to be
analyzed.
Figure 3.10. PSO neighborhood topology types;
(a) Ring topology, (b) Star topology (global neighborhood)
p
1
p
i-1
p
i
p
i+1
p
N
p
1
p
i-1
p
i
p
i+1
p
N
(a)
(b)
40
There are two ways of determining the global best (P
g
) position. The first is to
only update P
g
after every iteration. The particles are all updated in parallel, and only
then a P
g
is determined. This method is called a ?synchronous? update. The second
method is to determine the P
g
after each particle has been updated. This form is called an
?asynchronous? update. Asynchronous updates tend to find solutions quicker (in terms of
number of iterations), due to their ability to use a newly found P
g
in subsequent particle
updates immediately in the same iteration, instead of having to wait for the next iteration,
as is done in synchronous updates. This advantage usually makes asynchronous updates a
better choice for a standard PSO.
3.2.1. PSO: The Basic Algorithm
As it was mentioned in section 3.2, there are two methods of velocity update in
PSO: synchronous, and asynchronous. The pseudo-codes describing the baseline
algorithms for these PSO?s are shown in Fig. 3.11 and 3.12, respectively. Notice that the
synchronous method is easier to implement, since the algorithm works with blocks of
data (when all particle positions have been just updated), after which the global best (P
g
)
is determined. The asynchronous method is somewhat more difficult to implement, as the
algorithm must evaluate the particles one by one, where P
g
is updated as soon as a
particle is determined to be better than the current global best. As was mentioned
previously, the asynchronous method generally provides a better performance, and so it is
the usually recommended update method.
41
Figure 3.11. PSO pseudo-code for synchronous implementation
Figure 3.12. PSO pseudo-code for asynchronous implementation
Procedure PSO_Async {
Initialize x; /* swarm of size S */
Evaluate(x_fitness);
P = x;
Pg = best_P;
t = 1;
while (t <= number_of_iterations){
for (i=1, i<=S, i++){
/* evaluate every particle */
Calculate v
i
;
x
i
= x
i
+ v
i
;
Evaluate(x
i
_fitness);
Pi = best x
i
so far;
Pg = best P so far; }
t = t + 1; }
}
Procedure PSO_Sync {
Initialize x; /* create swarm */
while (t <= number_of_iterations){
Evaluate(x_fitness);
P = best x?s so far;
Pg = best P so far;
Calculate v;
x = x + v;
t = t + 1; }
}
42
CHAPTER 4
MATHEMATICAL MODEL DEVELOPMENT
4.1. Experimental Setup
In order to develop a suitable mathematical model, able to predict the dynamic
response of honeycomb sandwich composites, observations of individual experimental
hysteresis loops were first performed for numerous beam samples of these materials.
These experimental hysteresis loops were obtained by constructing a customized test rig
that allows the mounting of the beam samples with simply supported boundary
conditions, and mid-point excitation. A concept drawing and pictures of the test rig are
shown in Fig. 4.1 and 4.2, respectively.
Figure 4.1. Concept drawing of the test rig for hysteresis loops generation
FORCE TRANSDUCER
DYNAMIC SHAKER HEAD
ATTACHMENT
BEAM SAMPLE
TENSIONED STEEL WIRE
FRAME
BEAM/WIRE COUPLER
43
Figure 4.2. Pictures of the test rig for hysteresis loops generation;
(a) test rig, (b) wire attachment detail
(a)
(b)
44
Ni, Ko, and Wong
[12]
have performed system identification experiments with
hysteretic systems, where the displacement-controlled tests were shown to have a much
better accuracy than the force-controlled ones. Accordingly, the experimental tests
performed in this work consist of a simply-supported (by steel wires at both ends) beam
sample of the composite material, excited by a pure sinusoidal displacement wave with a
dynamic shaker. This applied displacement is measured at the mid-point of the beam with
a laser displacement vibrometer, and the resulting force is measured by a force transducer
positioned in-line between the shaker contact and the test sample. The coupling is
attained by using a rigid point connector, bonded to the beam with Cyanoacrylate. The
signals are then collected using a dynamic signal analyzer, and also alternatively using a
PC with a DAQ card and National Instrument?s LabView v8.0. A schematic and a picture
of the complete experimental setup are shown in Fig. 4.3 and 4.4, respectively.
Figure 4.3. Experimental setup schematic, for hysteresis loops generation
HP 35665A DYNAMIC SIGNAL ANALYZER
POLYTEC OFV 2610
VIBROMETER CONTROLLER
POLYTEC OFV 300
LASER VIBROMETER HEAD
KISTLER 5010B
CHARGE AMPLIFIER
KISTLER 9212
FORCE TRANSDUCER
LDS PA 500L POWER
AMPLIFIER
SHAKER
SAMPLE
45
Figure 4.4. Pictures of the experimental setup for hysteresis loops generation;
(a) complete setup, (b) detail of the excitation/measurement contact point
SHAKER
FORCE
TRANSDUCER
COMPOSITE SAMPLE
LASER BEAM
(b)
(a)
46
A very important experimental consideration is explained as follows. A charge-
mode signal conditioning amplifier has a time constant (TC) determined by TC=C
r
?R
t
,
where C
r
=range capacitor capacitance, and R
t
=time constant resistor resistance. This
time constant has a direct effect on the output of the amplifier, in relation to the
expressions given by Eqs. (4.1) and (4.2), for its voltage transfer and phase responses,
respectively.
( )
()()
2
2TC
,
12 TC
out
in
f
V
V
f
?
?
?
=
+?
(4.1)
()
1
1
tan .
2TCf
?
?
?
??
=
??
?
??
(4.2)
The time constant has the same effect as a single-pole high-pass filter, whose cutoff
frequency is defined by Eq. (4.3).
()
1
.
2TC
c
f
?
=
?
(4.3)
According to Eqs. (4.1) to (4.3), the higher the time constant, the better response the
charge amplifier will have, considering that the frequency range of operation increases
(as it improves the low frequency response), the voltage transfer approaches to unity, and
the phase approaches to zero. This last result (regarding the phase output) is critical for
hysteresis analysis, as the width of the loops are directly related to the phase between
both signals. So, any external source of signal phase shift may cause the experimental
data to be incorrect, hence leading to erroneous system identification results. If the
frequency of the signal is sufficiently high, these issues barely affect the response of the
amplifier, but considering that the experimental analysis of hysteresis is usually
47
performed at relatively low frequencies, the time constant must be always set to the
largest possible value.
4.2. Mathematical Model Baseline
There have been a number of earlier attempts at representing the dynamics of
sandwich composites with nonlinear models, such as the Bouc/Wen hysteresis model
[2]
.
However, the almost perfectly elliptical hysteresis loops that were obtained
experimentally, along with a lack of consistent convergence in many of the system
parameters for the nonlinear model, led to the conclusion that a simpler semi-empirical
equivalent linear model would yield a much better representation of these materials
behavior. A baseline model consisting of a linear 1-DOF mechanical system with an
equivalent viscous damping has been used as a starting point. The differential equation
for such system is shown in Eq. (4.4), where ?=damping ratio, ?
n
=natural frequency of
vibration, and M
eq
=equivalent mass.
2
()
2 .
nn
eq
Ft
xxx
M
?? ?++= (4.4)
Since a sinusoidal displacement excitation is to be used in this study, the time functions
for x, x , and x are given by Eq. (4.5), where X
0
=amplitude, and ?=excitation frequency.
( )
()
()
0
0
2
0
sin ,
cos ,
sin .
xX t
xX t
x Xt
?
??
??
=?
=?
=? ?
(4.5)
Substituting Eq. (4.5) into Eq. (4.4), a closed-form solution for F(t) is given by Eq. (4.6).
()()
2
2
22 1
0 22
2
( ) 2 sin tan .
n
eq n n
n
Ft M X t
?? ?
?? ??? ?
??
?
??
???
=?? ? + ?? +
??
??
?
??
??
(4.6)
48
Additional tailoring of this expression to the observed experimental behavior has
been performed, as explained in the following. It should be pointed out first that, since
the experimental setup consists of a point excitation at the beam center, and the sample is
simply-supported at both ends, the Euler?s beam equation has been approximated into
this equivalent 1-DOF equation. The mid-point excitation eliminates the presence of the
second natural mode of vibration of the beam, thus having the first and third modes well
separated, by a factor of nine. This allows the infinite-DOF system that characterizes a
beam to be approximated into a 1-DOF system (around the first mode of vibration) with
sufficient accuracy, as can be observed in the compliance FRF for this case, in Fig. 4.5.
10
1
10
2
10
3
-160
-150
-140
-130
-120
-110
-100
-90
-80
-70
-60
Frequency [Hz]
C
o
m
p
l
i
a
nc
e M
a
gn
i
t
u
d
e
[
d
B
]
Beam
1-DOF
f
1
f
3
=9?f
1
Figure 4.5. Compliance FRF comparison, between a mid-point excited beam, and an
equivalent 1-DOF system
49
4.3. Preliminary Experimental Observations
As a first step in identifying appropriate modifications to the baseline model,
several hysteresis loops have been obtained experimentally for different material
samples, at different excitation frequencies and displacement amplitudes. Two unusual
phenomena have been noted after the preliminary identification of the system parameters
for individual loops, using an implementation of Genetic Algorithms, as it will be
described in more detail in section 5.1. Both, the resonance frequency (?
n
) and the
damping ratio (?) are linear functions of the excitation amplitude, as can be seen in the
graphs of Figs. 4.6 and 4.7, respectively, as an example of the results obtained for one of
the composite beam samples to be studied. Similar behavior has been observed for
several other samples as well.
86.78
86.8
86.82
86.84
86.86
86.88
86.9
86.92
0.04 0.06 0.08 0.1 0.12 0.14 0.16
Excitation amplitude [mm]
f
n
[Hz]
Figure 4.6. Observed behavior of a sandwich composite beam; plot of f
n
=?
n
/2? vs. X
0
50
0.0055
0.006
0.0065
0.007
0.0075
0.008
0.0085
0.009
0.04 0.06 0.08 0.1 0.12 0.14 0.16
Excitation amplitude [mm]
Da
mp
in
g r
a
ti
o
Figure 4.7. Observed behavior of a sandwich composite beam; plot of ? vs. X
0
The damping ratios at fixed excitation amplitudes, however, remained relatively constant
for frequencies around the first resonance, as in the example shown in Fig. 4.8.
0.005
0.0055
0.006
0.0065
0.007
0.0075
0.008
0.0085
0.009
80 82 84 86 88 90 92
Frequency [Hz]
Damping Ratio
0.05 mm
0.10 mm
0.15 mm
Resonance Frequencies
Figure 4.8. Observed behavior of a sandwich composite beam; plot of ? vs. f at different
excitation amplitudes
51
The probable cause of the amplitude dependence of the resonance frequency are
the nonlinear effects, which arise from the flexing of the wire supports at both ends of the
beam. This produces a decrease in the effective stiffness of the system. This effect is
illustrated in Fig. 4.9, with both (a) the beam sample, and (b) an equivalent system
consisting of a mass attached to two springs in series, and a viscous damper.
Figure 4.9. Equivalent experimental mechanical system, considering the stiffness of the
supports
Then the total effective stiffness of the system is given by Eq. (4.7), where K
b
is the beam
equivalent stiffness, K
s
is the supports nominal stiffness, and A is a stiffness decay factor.
1
0
11
.
bs
K
KKAX
?
??
=+
??
??
??
(4.7)
This expression for the stiffness leads to a behavior similar to that shown in Fig. 4.6.
As for the amplitude dependence of the damping ratio, this effect was further
verified (by generating compliance frequency response plots of some beam samples) that,
for increasing vibration amplitudes, the peaks of the first resonance of the system also
increased in level, as in the example shown in Fig. 4.10, for one of the beam samples.
Hence, according to the half-power method damping equation given by Eq. (2.7), the
-3dB point frequencies are located progressively in narrower bandwidths, reducing the
value of the damping ratio as the excitation level is increased.
?
K
b
, C
b
, M
eq
K
s
/2 K
s
/2
M
eq
K
s
(x
pk
) Kb
C
b
(a) (b)
52
80 85 90 95 100 105
-10
-5
0
5
10
15
20
Frequency [Hz]
C
o
mp
lia
n
c
e
[
d
B
]
HIGH Amp. (damping=0.00719)
MEDIUM Amp. (damping=0.00774)
LOW Amp. (damping=0.00806)
Figure 4.10. Overlay plot example of the compliance frequency response of a sandwich
composite beam sample around its first resonance, at different excitation amplitudes
From the behavior observed in Figs. 4.6 to 4.8, functional expressions for the
damping ratio (?) and the resonance frequency (?
n
=2?f
n
) can be developed.
4.4. Final Model
For the values of the damping ratio ?, and for the resonance frequency ?
n
, linear
functions in X
0
will be used (Eq. (4.8) and (4.9), respectively, where S
?
and S
?
are
referred as the slopes of these linear functions, respectively). As a result, the force
response can be represented by the steady-state expression shown in Eq. (4.10), for a
sinusoidal displacement excitation, as given by Eq. (4.5).
53
00 0
() ,XSX
?
? ?= +? (4.8)
00 0
() ,
n
XSX
?
? ?= +? (4.9)
()()
2
2
22
000
1 00
22
0
() ( ) 2 ( ) ( )
2( ) ( )
sin tan .
()
eq pk n n
n
n
Ft M x X X X
XX
t
X
?????
?? ?
?
??
?
=?? ? + ? ??
??
????
+
????
?
??
??
(4.10)
Therefore, the system parameters to be identified is the set composed by: (?
0
, S
?
, ?
0
, S
?
).
54
CHAPTER 5
IMPLEMENTATION AND RESULTS
5.1. Optimization Algorithms Implementation
In order to identify the system parameters from the model given by Eqs. (4.8) to
(4.10) for specific beam samples of honeycomb sandwich composite materials, computer
programs have been written in Matlab R2006a, implementing customized codes for GA
and PSO (as shown in Appendix A). It is expected that, in general, the search algorithms
will find parametric solutions relatively close to each other for different runs of the
computer codes, simultaneously for different excitation conditions. It should be pointed
out that the equivalent mass of the system (M
eq
) is a 1-DOF approximation of the
continuum of the simply-supported beams, which happens to be equal to the total mass of
the samples being analyzed.
Some additional features were also implemented from the baseline characteristics
of both algorithms, which greatly improved their search performance. First, a large initial
population (at iteration 0) is initially created, in order to find good candidates for the
upcoming iterations. For PSO, another feature is that when a ?particle position update?
attempts to go beyond the search bounds set by the user, the particle remains motionless
only for those parameters (not set to x
min
or x
max
, as is more typically done), and V
id
is set
to zero for the next velocity calculation. In addition to the regular offsprings that are
created in both algorithms, an additional set of candidate solutions (CS) is generated,
55
which consists of several different combinations of Gaussian-mutated versions of the
current best partial solution (BPS), with the intention of performing an exhaustive local
search around the BPS for an even better solution. The number of members of this group
depends on an ?add-on level? (L) set by the user, which is the maximum number of
parameters to be mutated in the BPS in order to create additional candidate solutions, as
in the example shown in Table 5.1, for a search space of 4 parameters, and L=3.
Table 5.1. Mutated parameters of the best partial solution (of 4 parameters, p
i
, where
i=1,2,3,4), for the creation of an additional set of candidate solutions (L=3)
for local search
X
X
X
X
X X
X X
X X
X X
X X
X X
X X X
X X X
X X X
X X X
Mutated parameters from the
best partial solution
p
1
p
2
p
3
p
4
LEVEL 1
LEVEL 2
LEVEL 3
56
For a search space of size N, the number of added individuals is given by Eq. (5.1).
11
!
AddPopNumber .
!( )!
LL
ii
N
N
i iN i
==
??
??
?
??
??
(5.1)
In the extreme case, where L=N (as used in the present study), the result may be
simplified as: AddPopNumber=2
N
-1 (i.e. for N=4 ? AddPopNumber=15). Finally, the
method of error calculation between the trial and model was performed by using the
standard Ordinary Least Squares (OLS) error function given by Eq. (3.3), comparing the
difference between the force data for every displacement sample value (Fig. 5.1) in both,
the top and bottom parts of the hysteresis loop. Interpolations are performed in order to
enforce the requirement that the model must match the trial data sample points.
Figure 5.1. Comparison of hysteresis loops for error calculation (displacement loading);
trial, model, sample points, error
X
min
X
max
Force
Displacement
57
The characteristics and optimization parameters in the PSO routine created were
taken mainly from the recommendations of Carlisle and Dozier
[47]
, typical values used
by other researchers, and some specific values considered appropriate for this problem.
These characteristics are:
? Asynchronous PSO
? Global topology
? ?
1
=?
2
=2.05
As for the GA, several versions of it have been tested, but the best results were
obtained for the following definitive implementation:
? Real-coded GA
? Rank-based parent selection
? Single-point BLX-? crossover (with ?=0.3)
?
()()1??++ survival selection strategy
? Elitist strategy
? Mutation rate=5%
? Gaussian mutation amount=30% (amount allowed for the mutated parameters to
move).
Both algorithms shared the following optimization parameters, which provided
with sufficiently accurate results:
? Initial population/swarm=6000 (at iteration 0)
? Population/swarm size=500 (after iteration 0)
58
? Iterations=18
? Additional-population Gaussian mutation amount=20% (amount allowed for the
mutated parameters of the additional population to move).
? Additional-population ?add-on? level: L=4
It should be noted that the main steps the GA takes at each iteration, are as
follows:
1. From the parent population of size ?, generate an offspring also of size ?.
2. Generate a set of ?AddPopNumber? chromosomes for the additional population,
based on the best current solution available.
3. Identify the best chromosome from the additional population generated.
4. Select the best chromosomes from a pool of size (?+?+1), composed of the parent
population, the offspring, and the best performing chromosome from the additional
population.
5.2. Benchmark Test of GA and PSO
Before making use of the computer optimization programs with actual
experimental data, computer-generated hysteresis loops (Synth set) have been created, at
40 different frequencies, and at 3 different excitation amplitudes each, making a total of
120 loops. The purpose of this is to carry out a benchmark test of the search performance
of both PSO and GA with a known global solution. Also, the influence of Gaussian
random noise present in the signals has been studied, in order to verify its effect on the
optimization results.
59
The equivalent mass used to generate the Synth data set is M
eq
=40e-3 kg, with
excitation amplitudes (X
0
) of 0.5, 1.0, and 1.5mm, and the excitation frequencies given in
Table 5.2.
Excitation Frequencies [Hz]
25 26 26.5 27 27.5 28 28.2 28.4 28.6 28.8
29 29.1 29.2 29.3 29.4 29.5 29.6 29.7 29.8 29.9
30 30.1 30.2 30.3 30.4 30.5 30.6 30.7 30.8 30.9
31 31.2 31.4 31.6 31.8 32 32.5 33 33.5 34
Table 5.2. Excitation frequencies of the computer-generated Synth hysteresis set
The system parameters used (global optimum) are shown in Table 5.3, where
?
0
=2??30.1=189.1239, and S
?
=2??(-200)=-1256.64.
?
0
S
?
?
0
S
?
0.0042 -0.4 189.1239 -1256.64
Table 5.3. System parameters used for generating the Synth hysteresis set
The data given in Table 5.3 correspond to the following constant values of
resonance frequency (?
n
) and damping ratio (?) of the general baseline model solution
(Eq. (4.6)) at specific excitation amplitudes:
? At X
0
=0.5mm : ?
n
=2??30.0 [rad/s] ; ?=0.0040
? At X
0
=1.0mm : ?
n
=2??29.9 [rad/s] ; ?=0.0038
? At X
0
=1.5mm : ?
n
=2??29.8 [rad/s] ; ?=0.0036
The noise levels used for this test (with respect of the peak values of force of each
hysteresis loop) were: 0% (clean signals), 10%, 20% and 30%. Finally, the search space
used (parameter bounds) is shown in Table 5.4.
60
Parameter Min Max
?
0
0 0.1
S
?
-10 10
?
0
100 500
S
?
-10000 10000
Table 5.4. Search space for the Synth hysteresis set
5.2.1. Benchmark Test Results
A total of 4 runs of each optimization routine were performed for each noise level
case of the Synth hysteresis set. The parametric results for each case, along with their
average values ( x ), standard deviations (?), and normalized standard deviations ( x? ,
for a relative comparison among the results of each parameter), are shown in Tables 5.5
to 5.8 for GA, and Tables 5.9 to 5.12 for PSO.
In order to determine the effect of the added Gaussian random noise on the
optimization results, scatter plots describing the spread ranges of the results of each
parameter, and for each optimization algorithm, are shown in Figs. 5.2 to 5.5, where the
dot over each line represents the average solution, and the red dotted line represents the
global optimum. Also, convergence plots (i.e. the fitness histories over all iterations) for
both algorithms and at each noise level (overlaying all 4 runs in each case, and with
matching scales, for comparison purposes), are shown in Figs. 5.6 to 5.9. Some plots
showing comparisons between trial hysteresis loops and model results (typical for both
optimization techniques), are shown in Figs. 5.10 to 5.13, for all 3 excitation amplitudes,
and all 4 noise levels. Finally, as a way of measuring the quality of the solutions of both
optimization techniques utilized, the differences of their averages from the actual global
optimum are compared in Table 5.13, for all 4 noise levels considered.
61
GA
0% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.00413 0.004444 0.004306 0.003081 0.00399 0.00062 0.155336
S
?
-0.4 -0.33347 -0.58421 -0.42632 0.372055 -0.24299 0.42289 1.740382
?
0
189.1239 189.1282 189.1831 189.222 189.1089 189.1605 0.051628 0.000273
S
?
-1256.64 -1256.03 -1311.81 -1309.71 -1298.78 -1294.08 26.00358 0.020094
Table 5.5. Optimization results for the Synth data set; GA, 0% noise
GA
10% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.00379 0.003934 0.004319 0.004104 0.004037 0.000227 0.056337
S
?
-0.4 -0.00494 -0.22136 -0.52046 -0.32666 -0.26836 0.214907 0.800830
?
0
189.1239 189.0764 189.1619 189.2238 189.1275 189.1474 0.061878 0.000327
S
?
-1256.64 -1199.93 -1303.86 -1327.24 -1256.3 -1271.83 56.29433 0.044262
Table 5.6. Optimization results for the Synth data set; GA, 10% noise
GA
20% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.004398 0.004363 0.004103 0.004045 0.004227 0.000179 0.042303
S
?
-0.4 -0.61225 -0.50622 -0.33247 -0.25812 -0.42726 0.161295 0.377506
?
0
189.1239 188.9366 189.1497 189.0965 189.1741 189.0892 0.106813 0.000565
S
?
-1256.64 -1092.83 -1243.7 -1226.6 -1303.53 -1216.66 88.90028 0.073069
Table 5.7. Optimization results for the Synth data set; GA, 20% noise
GA
30% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.004771 0.003594 0.003973 0.003385 0.003931 0.000611 0.155381
S
?
-0.4 -0.71675 0.208465 -0.25031 0.246317 -0.12807 0.452737 3.535067
?
0
189.1239 189.0589 189.052 189.1316 189.0822 189.0812 0.035996 0.000190
S
?
-1256.64 -1200.27 -1201.69 -1228.52 -1233.67 -1216.04 17.52276 0.014410
Table 5.8. Optimization results for the Synth data set; GA, 30% noise
62
PSO
0% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.004437 0.004611 0.004557 0.004614 0.004555 8.27E-05 0.018151
S
?
-0.4 -0.56995 -0.81059 -0.70964 -0.67819 -0.69209 0.099097 0.143185
?
0
189.1239 189.1421 189.0642 189.0937 189.1477 189.1119 0.039992 0.000211
S
?
-1256.64 -1302.46 -1192.85 -1245.96 -1288.61 -1257.47 49.3357 0.039234
Table 5.9. Optimization results for the Synth data set; PSO, 0% noise
PSO
10% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.004785 0.004713 0.003923 0.003881 0.004326 0.00049 0.113324
S
?
-0.4 -0.79555 -0.81552 -0.23406 -0.14349 -0.49716 0.358097 0.720290
?
0
189.1239 188.9828 189.0352 189.1041 189.088 189.0525 0.054989 0.000291
S
?
-1256.64 -1131.53 -1192.03 -1253.12 -1236.37 -1203.26 54.32264 0.045146
Table 5.10. Optimization results for the Synth data set; PSO, 10% noise
PSO
20% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.003963 0.004014 0.003709 0.00389 0.003894 0.000134 0.034335
S
?
-0.4 -0.14437 -0.26865 -0.02538 -0.24365 -0.17052 0.110643 0.648877
?
0
189.1239 189.0428 189.1218 189.1098 189.1164 189.0977 0.036953 0.000195
S
?
-1256.64 -1194.78 -1245.47 -1245.32 -1240.31 -1231.47 24.57532 0.019956
Table 5.11. Optimization results for the Synth data set; PSO, 20% noise
PSO
30% noise
ACTUAL RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.0042 0.004226 0.004108 0.004937 0.003767 0.004259 0.000492 0.115522
S
?
-0.4 -0.4518 -0.4938 -1.07537 -0.1138 -0.53369 0.399173 0.747946
?
0
189.1239 189.1569 189.1508 189.3972 189.1963 189.2253 0.116378 0.000615
S
?
-1256.64 -1264.21 -1262.62 -1442.37 -1347.87 -1329.27 85.26707 0.064146
Table 5.12. Optimization results for the Synth data set; PSO, 30% noise
63
3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5
x 10
-3
0
10
20
30
?
0
Range
N
o
i
s
e Le
v
e
l
(
%
)
?
0
Ranges vs. Noise Level (GA)
3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5
x 10
-3
0
10
20
30
?
0
Range
N
o
i
s
e Lev
el
(
%
)
?
0
Ranges vs. Noise Level (PSO)
Figure 5.2. Scatter plots from the results of the Synth set, for ?
0
; (a) GA, (b) PSO
(a)
(b)
64
-1 -0.8 -0.6 -0.4 -0.2 0 0.2
0
10
20
30
S
?
Range
No
i
s
e
L
e
v
e
l
(%
)
S
?
Ranges vs. Noise Level (GA)
-1 -0.8 -0.6 -0.4 -0.2 0 0.2
0
10
20
30
S
?
Range
No
i
s
e
L
e
v
e
l
(%
)
S
?
Ranges vs. Noise Level (PSO)
Figure 5.3. Scatter plots from the results of the Synth set, for S
?
; (a) GA, (b) PSO
(a)
(b)
65
188.9 188.95 189 189.05 189.1 189.15 189.2 189.25 189.3 189.35 189.4
0
10
20
30
?
0
Range
No
i
s
e
L
e
v
e
l
(%
)
?
0
Ranges vs. Noise Level (GA)
188.9 188.95 189 189.05 189.1 189.15 189.2 189.25 189.3 189.35 189.4
0
10
20
30
?
0
Range
No
i
s
e
L
e
v
e
l
(%
)
?
0
Ranges vs. Noise Level (PSO)
Figure 5.4. Scatter plots from the results of the Synth set, for ?
0
; (a) GA, (b) PSO
(a)
(b)
66
-1450 -1400 -1350 -1300 -1250 -1200 -1150 -1100 -1050
0
10
20
30
S
?
Range
No
i
s
e
L
e
v
e
l
(%
)
S
?
Ranges vs. Noise Level (GA)
-1450 -1400 -1350 -1300 -1250 -1200 -1150 -1100 -1050
0
10
20
30
S
?
Range
No
i
s
e
L
e
v
e
l
(%
)
S
?
Ranges vs. Noise Level (PSO)
Figure 5.5. Scatter plots from the results of the Synth set, for S
?
; (a) GA, (b) PSO
(a)
(b)
67
0 2 4 6 8 10 12 14 16 18
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Synth (0% noise) - GA
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Synth (0% noise) - PSO
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.6. Convergence plots from the results of the Synth set, for 0% noise;
(a) GA, (b) PSO
Iteration number
F
i
tness value
Iteration number
F
i
tness value
(a)
(b)
68
0 2 4 6 8 10 12 14 16 18
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Synth (10% noise) - GA
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Synth (10% noise) - PSO
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.7. Convergence plots from the results of the Synth set, for 10% noise;
(a) GA, (b) PSO
Iteration number
F
i
tness value
Iteration number
F
i
tness value
(a)
(b)
69
0 2 4 6 8 10 12 14 16 18
0.4
0.5
0.6
0.7
0.8
0.9
1
Synth (20% noise) - GA
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0.4
0.5
0.6
0.7
0.8
0.9
1
Synth (20% noise) - PSO
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.8. Convergence plots from the results of the Synth set, for 20% noise;
(a) GA, (b) PSO
Iteration number
F
i
tness value
Iteration number
F
i
tness value
(a)
(b)
70
0 2 4 6 8 10 12 14 16 18
0.9
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
Synth (30% noise) - GA
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0.9
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
Synth (30% noise) - PSO
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.9. Convergence plots from the results of the Synth set, for 30% noise;
(a) GA, (b) PSO
Iteration number
F
i
tness value
Iteration number
F
i
tness value
(a)
(b)
71
-5 0 5
x 10
-4
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
29.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-6
-4
-2
0
2
4
6
x 10
-3
30 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
30.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.05
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
0.05
29.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
29.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.05
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
0.05
30.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 -1 0 1 2
x 10
-3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
28 [Hz]
Disp. [m]
For
c
e [
N
]
-2 0 2
x 10
-3
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
29.3 [Hz]
Disp. [m]
For
c
e [
N
]
-2 0 2
x 10
-3
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
29.8 [Hz]
Disp. [m]
For
c
e [
N
]
-2 0 2
x 10
-3
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
30.3 [Hz]
Disp. [m]
For
c
e [
N
]
-2 0 2
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
32 [Hz]
Disp. [m]
For
c
e [
N
]
Figure 5.10. Comparison plots for some of the hysteresis loops of the Synth set; 0%
noise; 1st.row: 0.5mm; 2nd.row: 1.0mm; 3rd.row: 1.5mm; trial, model
72
-5 0 5
x 10
-4
-0.15
-0.1
-0.05
0
0.05
0.1
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
29.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-8
-6
-4
-2
0
2
4
6
8
x 10
-3
30 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
30.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
29.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
29.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
30.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 -1 0 1 2
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
29.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
29.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
30.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
Figure 5.11. Comparison plots for some of the hysteresis loops of the Synth set; 10%
noise; 1st.row: 0.5mm; 2nd.row: 1.0mm; 3rd.row: 1.5mm; trial, model
73
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
29.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-8
-6
-4
-2
0
2
4
6
8
10
x 10
-3
30 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
30.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
29.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
29.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
30.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 -1 0 1 2
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.15
-0.1
-0.05
0
0.05
0.1
29.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
29.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
30.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
Figure 5.12. Comparison plots for some of the hysteresis loops of the Synth set; 20%
noise; 1st.row: 0.5mm; 2nd.row: 1.0mm; 3rd.row: 1.5mm; trial, model
74
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
29.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.01
-0.008
-0.006
-0.004
-0.002
0
0.002
0.004
0.006
0.008
0.01
30 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
30.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
29.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
29.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
30.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 -1 0 1 2
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
28 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
29.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
29.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
30.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
32 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
Figure 5.13. Comparison plots for some of the hysteresis loops of the Synth set; 30%
noise; 1st.row: 0.5mm; 2nd.row: 1.0mm; 3rd.row: 1.5mm; trial, model
75
0% noise 10% noise 20% noise 30% noise
Parameter
GA PSO GA PSO GA PSO GA PSO
?
0
2.10E-4 3.55E-4 1.63E-4 1.26E-4 0.27E-4 3.06E-4 2.69E-4 0.59E-4
S
?
0.1570 0.2920 0.1320 0.0972 0.0273 0.2290 0.2720 0.1340
?
0
0.0366 0.0120 0.0235 0.0714 0.0347 0.0262 0.0427 0.101
S
?
37.4 0.834 15.2 53.4 40.0 25.2 40.6 72.6
Table 5.13. Differences between the resultant averages and the actual global optimum,
for GA and PSO, at different noise levels
5.3. Material Samples Description and Experimental Results
As explained in section 1.2, three sandwich composite material beams have been
studied. They consist of a paper honeycomb core filled with foam and external layers of
interlaced strips of carbon fibers, as shown in Fig. 1.1b. The samples used are described
in Table 5.14, based on the sketch shown in Fig. 5.14.
Figure 5.14. Sketch of a sandwich composite beam sample
PROPERTY SAMPLE 1 SAMPLE 2 SAMPLE 3
Length (L) [mm] 590 590 590
Width (W) [mm] 25 25 25
Core Thickness (H
1
) [mm] 6.5 6.5 13
Face Thickness (H
2
) [mm] 0.2 0.2 0.2
Faces per side 1 2 1
Mass [g] 28.7 36.8 48.0
Table 5.14. Properties of the sandwich composite beams studied
LENGTH (L)
WIDTH (W)
CORE THICKNESS (H
1
)
FACE THICKNESS (H
2
)
76
In a way similar to that described in section 5.2 for the Synth set, several
hysteresis loops have been experimentally obtained for all 3 beam samples. For samples 1
and 2, a total of 40 frequencies have been used, at 3 different displacement excitation
amplitudes each, making a total of 120 hysteresis loops for each of these samples. For
sample 3, however, 31 frequencies have been used, at 3 different displacement excitation
amplitudes each, making a total of 93 hysteresis loops for that sample. The excitation
frequencies used for each sample sweep a range that crosses their first resonance
frequency, and they are shown in Tables 5.15 to 5.17. The displacement excitation
amplitudes used are indicated in Table 5.18, and the search spaces used for the beam
samples are shown in Table 5.19.
Excitation Frequencies [Hz] - Sample 1
28 29 30 30.5 31 31.2 31.4 31.6 31.8 32
32.2 32.4 32.6 32.7 32.8 32.9 33 33.1 33.2 33.3
33.4 33.5 33.6 33.7 33.8 33.9 34 34.2 34.4 34.6
34.8 35 35.2 35.4 35.6 36 36.5 37 38 39
Table 5.15. Excitation frequencies used with sample 1
Excitation Frequencies [Hz] - Sample 2
37 38 38.5 39 39.5 40 40.2 40.4 40.6 40.8
41 41.2 41.4 41.5 41.6 41.7 41.8 41.9 42 42.1
42.2 42.3 42.4 42.5 42.6 42.7 42.8 42.9 43 43.2
43.4 43.6 43.8 44 44.2 44.4 45 45.5 46 47
Table 5.16. Excitation frequencies used with sample 2
77
Excitation Frequencies [Hz] - Sample 3
81 82 83 84 84.5 85 85.6 85.8
86 86.1 86.2 86.3 86.4 86.5 86.6 86.7
86.8 86.9 87 87.1 87.2 87.3 87.4 87.6
87.8 88 88.2 88.4 89 89.5 90
Table 5.17. Excitation frequencies used with sample 3
Amplitude [mm] Sample 1 Sample 2 Sample 3
Low 0.5 0.2 0.05
Medium 0.8 0.5 0.10
High 1.0 0.8 0.15
Table 5.18. Excitation displacement amplitudes used with the beam samples studied
Sample 1 Sample 2 Sample 3
Parameters
Min Max Min Max Min Max
?
0
0 0.05 0 0.1 0 0.05
S
?
-20 20 -10 10 -100 100
?
0
100 400 100 500 400 700
S
?
-4000 4000 -10000 10000 -10000 10000
Table 5.19. Search spaces for all 3 beam samples studied
5.3.1. Optimization Results from Experimental Data
A total of 4 runs of each optimization routine were performed, for the
identification of the system parameters of all 3 beam samples studied. As it has been done
in section 5.2.1, the parametric results for each case, along with their average values ( x ),
standard deviations (?), and normalized standard deviations ( x? ), are shown in
78
Tables 5.20 to 5.22 for GA, and Tables 5.23 to 5.25 for PSO. Scatter plots, comparing the
GA and PSO spread ranges in the results for each parameter, are shown in Figs. 5.15 to
5.17. Also, convergence plots for both algorithms (overlaying all 4 runs in each case, and
with matching scales, for comparison purposes), and for each beam sample, are shown in
Figs. 5.18 to 5.20. Some plots showing comparisons between trial hysteresis loops and
model results (typical for both optimization techniques), are shown in Figs. 5.21 to 5.23,
for all 3 excitation amplitudes.
GA
SAMPLE 1
RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.00554 0.005378 0.005608 0.004806 0.005333 0.000364 0.068322
S
?
-2.41365 -2.22411 -2.46943 -1.77114 -2.21958 0.316861 0.142757
?
0
209.9339 209.9262 209.8417 210.3053 210.0018 0.206592 0.000984
S
?
-1257.99 -1249.85 -1160.66 -1702.94 -1342.86 244.0698 0.181754
Table 5.20. GA optimization results for beam sample 1
GA
SAMPLE 2
RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.004618 0.004864 0.005017 0.004822 0.00483 0.000164 0.034017
S
?
-0.08142 -0.31777 -0.61774 -0.29375 -0.32767 0.220629 0.673327
?
0
265.8243 265.8762 265.7912 265.7695 265.8153 0.046433 0.000175
S
?
-1479.36 -1512.6 -1409.21 -1378.99 -1445.04 61.61231 0.042637
Table 5.21. GA optimization results for beam sample 2
GA
SAMPLE 3
RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.009875 0.009658 0.00995 0.009857 0.009835 0.000125 0.012664
S
?
-25.7341 -24.0335 -26.4077 -25.7074 -25.4707 1.01142 0.039709
?
0
546.2193 546.3223 546.3425 546.389 546.3183 0.071626 0.000131
S
?
-5584.3 -6259.44 -6498.74 -6853.79 -6299.07 535.4282 0.085001
Table 5.22. GA optimization results for beam sample 3
79
PSO
SAMPLE 1
RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.004985 0.005251 0.005896 0.00591 0.00551 0.000466 0.084637
S
?
-1.81823 -2.04107 -2.87999 -3.05194 -2.44781 0.609253 0.248897
?
0
210.035 209.8965 209.9815 209.9329 209.9615 0.060139 0.000286
S
?
-1376.72 -1206.84 -1314.3 -1258.59 -1289.11 73.04978 0.056667
Table 5.23. PSO optimization results for beam sample 1
PSO
SAMPLE 2
RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.004922 0.004993 0.004957 0.004672 0.004886 0.000146 0.029812
S
?
-0.31858 -0.6121 -0.3394 -0.10522 -0.34382 0.207813 0.604417
?
0
265.8805 265.8835 265.8225 265.8239 265.8526 0.034008 0.000128
S
?
-1559.42 -1508.21 -1481.76 -1462.2 -1502.9 42.13699 0.028037
Table 5.24. PSO optimization results for beam sample 2
PSO
SAMPLE 3
RUN 1 RUN 2 RUN 3 RUN 4 x
?
? x
?
0
0.009784 0.009884 0.009892 0.009974 0.009884 7.79E-05 0.007879
S
?
-25.0854 -25.9808 -26.1065 -26.3107 -25.8709 0.540977 0.020911
?
0
546.2392 546.2741 546.2434 546.4486 546.3013 0.099383 0.000182
S
?
-5740.92 -5792.44 -5850.91 -7493.71 -6219.49 850.6676 0.136774
Table 5.25. PSO optimization results for beam sample 3
80
4.8 5 5.2 5.4 5.6 5.8 6
x 10
-3
PSO
?
0
Range (Sample 1)
GA
-3.2 -3 -2.8 -2.6 -2.4 -2.2 -2 -1.8 -1.6
PSO
S
?
Range (Sample 1)
GA
209.8 209.9 210 210.1 210.2 210.3 210.4
PSO
?
0
Range (Sample 1)
GA
-1800 -1700 -1600 -1500 -1400 -1300 -1200 -1100
PSO
S
?
Range (Sample 1)
GA
Figure 5.15. Parameters scatter plots from the optimization results of beam sample 1
81
4.6 4.65 4.7 4.75 4.8 4.85 4.9 4.95 5 5.05
x 10
-3
PSO
?
0
Range (Sample 2)
GA
-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0
PSO
S
?
Range (Sample 2)
GA
265.76 265.78 265.8 265.82 265.84 265.86 265.88 265.9
PSO
?
0
Range (Sample 2)
GA
-1560 -1540 -1520 -1500 -1480 -1460 -1440 -1420 -1400 -1380 -1360
PSO
S
?
Range (Sample 2)
GA
Figure 5.16. Parameters scatter plots from the optimization results of beam sample 2
82
9.65 9.7 9.75 9.8 9.85 9.9 9.95 10
x 10
-3
PSO
?
0
Range (Sample 3)
GA
-26.5 -26 -25.5 -25 -24.5 -24
PSO
S
?
Range (Sample 3)
GA
546.2 546.25 546.3 546.35 546.4 546.45
PSO
?
0
Range (Sample 3)
GA
-7500 -7000 -6500 -6000 -5500
PSO
S
?
Range (Sample 3)
GA
Figure 5.17. Parameters scatter plots from the optimization results of beam sample 3
83
0 2 4 6 8 10 12 14 16 18
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Sample 1 - GA
Iteration number
F
i
tn
e
s
s va
l
u
e
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Sample 1 - PSO
Iteration number
F
i
tn
e
s
s va
l
u
e
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.18. Optimization convergence plots for beam sample 1
84
0 2 4 6 8 10 12 14 16 18
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Sample 2 - GA
Iteration number
F
i
tn
e
s
s va
l
u
e
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Sample 2 - PSO
Iteration number
F
i
tn
e
s
s va
l
u
e
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.19. Optimization convergence plots for beam sample 2
85
0 2 4 6 8 10 12 14 16 18
0
0.001
0.002
0.003
0.004
0.005
0.006
0.007
0.008
0.009
0.01
Sample 3 - GA
Iteration number
F
i
tn
e
s
s va
l
u
e
RUN 1
RUN 2
RUN 3
RUN 4
0 2 4 6 8 10 12 14 16 18
0
0.001
0.002
0.003
0.004
0.005
0.006
0.007
0.008
0.009
0.01
Sample 3 - PSO
Iteration number
F
i
tn
e
s
s va
l
u
e
RUN 1
RUN 2
RUN 3
RUN 4
Figure 5.20. Optimization convergence plots for sample beam 3
86
-5 0 5
x 10
-4
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
31 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
32.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-6
-4
-2
0
2
4
6
x 10
-3
33.3 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
33.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
36 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
31 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
32.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-8
-6
-4
-2
0
2
4
6
8
x 10
-3
33.2 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
33.7 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
36 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
31 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
32.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 -1 0 1 2
x 10
-3
-8
-6
-4
-2
0
2
4
6
8
x 10
-3
33.2 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-3
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
33.7 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 -1 0 1
x 10
-3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
0.25
36 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
Figure 5.21. Optimization comparison plots for some of the hysteresis loops of sample 1;
1st.row: 0.5mm; 2nd.row: 0.8mm; 3rd.row: 1.0mm; trial, model
87
-2 -1 0 1 2
x 10
-4
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
39 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2 4
x 10
-4
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
41.7 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-4 -2 0 2 4
x 10
-4
-6
-4
-2
0
2
4
6
x 10
-3
42.2 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
42.7 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-2 0 2
x 10
-4
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
45 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
39 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
41.8 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
42.2 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
42.6 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-4
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
45 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
39 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.05
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
0.05
41.7 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
42.1 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.05
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
0.05
42.5 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-3
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
45 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
Figure 5.22. Optimization comparison plots for some of the hysteresis loops of sample 2;
1st.row: 0.2mm; 2nd.row: 0.5mm; 3rd.row: 0.8mm; trial, model
88
-5 0 5
x 10
-5
-0.05
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
0.05
84 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-5
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
86.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-5
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
86.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-5
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
87.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-5 0 5
x 10
-5
-0.04
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
0.04
88.9 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-4
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
84.15 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
86.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
86.85 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
87.35 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
88.75 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 -0.5 0 0.5 1
x 10
-4
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
84.15 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
86.4 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
86.85 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.03
-0.02
-0.01
0
0.01
0.02
0.03
87.35 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
-1 0 1
x 10
-4
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
88.75 [Hz]
Disp. [m]
Fo
r
c
e
[
N
]
Figure 5.23. Optimization comparison plots for some of the hysteresis loops of sample 3;
1st.row: 0.05mm; 2nd.row: 0.10mm; 3rd.row: 0.15mm; trial, model
89
CHAPTER 6
DISCUSSION AND CONCLUSIONS
In this dissertation, a study to identify model parameters related to beams of
sandwich composite materials has been presented. The experimental testing apparatus
and some representative hysteresis loops that were obtained have been also presented. A
simplified mathematical model for the representation of the dynamic response of a
honeycomb-core sandwich composite material beam sample under a displacement
excitation at its mid-point has been developed and tailored based upon experimental
observations of the system behavior.
It can be seen from the simulation results obtained from both the PSO and GA,
over a set of computer-generated hysteresis loops (Synth set), that the parametric
solutions that were found are usually quite close (among different runs) to the actual
values of the parameters when the global optimum was known. In addition, it has been
observed that in spite of the addition of high levels of Gaussian random noise to the
simulated signals, the solutions in these cases still remain close to the global optimum.
However, the noise in many cases in the simulation study was deliberately set at very
high levels, and the solutions obtained under these conditions are still quite consistent and
acceptable for most applications.
The experimentally identified parameters are also quite consistent and lend
confidence to both the proposed mathematical model and to the identification algorithms.
90
From the experimental results obtained with the optimization techniques in Tables 5.20 to
5.25, it can be observed that nominal values of the damping ratio for each sample (i.e. for
samples 1, 2, and 3, at ?medium? excitation amplitudes) are 0.0035, 0.0047, and 0.0073,
respectively, when the results are evaluated using Eq. (4.8). These damping ratios are
relatively high even for composite materials and this phenomenon is probably caused by
the relative motion at the various interfaces within the composite and their interaction
with the foam filling in the honeycomb core structure. This makes this type of material a
good choice when low weight and good vibration energy dissipation performance are of
importance for the construction of a structure, such as what is needed in aircrafts.
Considering sample 1 as a standard for comparison with the other samples
utilized, this result shows that the addition of an extra face layer (represented by sample
2) increased the damping ratio by about 34%, whereas doubling the core thickness
(represented by sample 3) led to an increase of the damping ratio by about 108%.
The resulting damping ratios from the hysteresis loop approach are consistent
with the results obtained by using the half-power bandwidth method, described in section
2.2. Compliance frequency response plots for each beam sample, near their first
resonance frequency, are shown in Fig. 6.1. Each of these plots are shown as an overlay
of the frequency response of a beam at different excitation levels, generically named
?low?, ?medium?, and ?high?, respectively. The resonance frequency (in magenta color)
and the -3dB frequency points (in cyan color), for a medium excitation amplitude, are
also shown in each case, whose values (when evaluated into Eq. (2.7)) lead to similar
damping ratio results as with the hysteresis loops technique.
91
20 25 30 35 40 45
-75
-70
-65
-60
-55
-50
-45
-40
-35
-30
Compliance Frequency Responses - Sample 1
Frequency [Hz]
C
o
m
p
l
i
anc
e [
d
B
]
LOW
MEDIUM
HIGH
f
0
=33.25[Hz]
f
1
=33.15[Hz] f
2
=33.38[Hz]
25 30 35 40 45 50
5
10
15
20
25
30
35
40
45
Frequency [Hz]
C
o
m
p
l
i
anc
e [
d
B
]
Compliance Frequency Responses - Sample 2
LOW
MEDIUM
HIGH
f
1
=42.01[Hz] f
2
=42.41[Hz]
f
0
=42.13[Hz]
80 85 90 95 100 105
-15
-10
-5
0
5
10
15
20
Frequency [Hz]
C
o
mp
lia
n
c
e
[
d
B
]
Compliance Frequency Responses - Sample 3
LOW
MEDIUM
HIGH
f
1
=86.18[Hz] f
2
=87.45[Hz]
f
0
=86.85[Hz]
Figure 6.1. Overlay compliance FRF plots of all 3 beam samples utilized, for low,
medium, and high excitation amplitudes
92
It can be also verified analytically from Fig. 6.1 that the higher damping ratios
correspond to the cases with the smoothest resonance peaks. This result was expected,
bearing in mind the explanation in section 2.2, concerning the relationship between
damping and the sharpness of the resonance peaks.
With regard to the effect of Gaussian random noise in the system response signal
(force, in this case), it was expected that, in general, both optimization algorithms would
perform well in the search for an ?average? hysteresis loop. The explanation for this lies
in the nature of Gaussian random noise, as explained next. The total mean square error of
a specific candidate solution (CS) to a set of noisy sample points (from a trial set of size
N) is given by:
()
2
1
1
?
,
N
CS i i
i
MSE F f
N
=
=?
?
(6.1)
where ()
?
iii
FFn=+ is a noise-contaminated trial sample point, F
i
=clean trial point i,
n
i
=Gaussian random noise of point i, and f
i
=corresponding sample point i from the model
solution. We have then, omitting the i subscripts, for convenience:
() []()
2
2
2
?
2 .Ff Fn f Ff Ffn
??
???=+?=?+??
??
??
(6.2)
Since Gaussian random noise has an expected value (i.e. the estimated average for an
infinitely long data set) of E(n)=0, and a standard deviation of 1, we have:
()()() ()
22
11
2 .
NN
CS ii ii ii
ii
MSE F f F f E n F f
==
=?+??=?
??
(6.3)
This result shows that, in theory, Gaussian random noise present in a system response
signal should not interfere in the quality of the system identification results. In other
words, the optimization algorithms should be able to find solutions similar to those of the
93
corresponding noiseless signals. This conclusion was verified from the results of the
Synth computer-generated hysteresis set, shown in Tables 5.5 to 5.12, and Figs. 5.10 to
5.13.
An interesting optimization performance result can be observed from all the
convergence plots shown in this work, in Figs. 5.6 to 5.9 (for the Synth hysteresis set),
and 5.18 to 5.20 (for the experimental data). In spite of the fact that both optimization
algorithms used provided sufficiently accurate results, and relatively low standard
deviations, PSO produced a much better search performance than GA. This may be
verified by observing the slopes of the convergence plots in the initial iterations, which
are much steeper in PSO than in GA, in all cases. This means that PSO is able to find
lower-error candidate solutions much faster than GA. Besides, PSO is much easier to
implement than GA, as it requires fewer steps to update the candidate solutions during
iterations, therefore requiring several fewer lines of computer code. This makes PSO a
definite choice between both algorithms, which is in agreement with what has been
reported by several researchers
[37-38]
, as explained in section 2.4.
The large scale of the present study (considering the amount of hysteresis loops
simultaneously used in each test), and the relatively high dimensionality of the model,
makes it a difficult problem to solve for any optimization algorithm, in general. It seems
fair then to consider that both, GA and PSO, have been shown to be good choices for
solving the problems presented in this study, with low error values and few iterations.
Some fundamental contributions of the presented research work can be identified
as follows:
94
? Developed techniques for the determination of the damping characteristics of
sandwich composite materials from hysteresis curves.
? Created computer programs for the implementation of different customized heuristic
optimization methods, such as GA and PSO.
? Tested the performance and reliability of these optimization methods when used for
system characterization of the semi-empirical model that has been developed.
? Validated the experimental results utilizing these two different optimization methods.
? Analyzed the related theory in order to develop all the necessary mathematical
expressions needed for the present research.
Finally, some potentially important applications of this research work are:
? The extraction of important information of the physical/mechanical properties of the
vibratory dynamics of composite materials (or any type of material, in general) that
may be used for design and for further research in vibration analysis.
? The ability to more accurately predict the dynamical behavior in time of such
materials through computer simulation once the system has already been
characterized (parameters of the system are finally known).
? A better understanding of the application of optimization methods in solving
engineering problems, focused on two powerful techniques, GA and PSO, which have
shown to be of popular usage, and proven to be successful in a wide array of
problems. This knowledge may be utilized for future research as well, in order to
substantially reduce the effort required to implement such techniques, this way being
95
able to create optimized computer programs, thus reducing the CPU time required to
converge to good solutions.
6.1. Suggestions for Future Research
1. Although the present work deals with the study of a generalized set of honeycomb-
core sandwich composite beams, several more samples should be studied (following
the presented approach), having different combinations of core materials and
structures, number of face layers, and dimensions. Composite materials of different
nature should also be studied, outside of the honeycomb-core category.
2. A study on the variation of the optimization parameters for both algorithms should be
carried out, in order to analyze their effect on the optimization performances. More
specifically, for GA, different population sizes, number of iterations, parent selection
techniques, mutation rates, and survivor selection strategies should be tested, in order
to find combinations that lead to fewer function evaluations and less CPU time
required. As for PSO, different swarm sizes, in conjunction with different cognition
and social component factor values should also be tested.
3. In order to improve the performance of the time consuming task of finding optimal
solutions for large sets of trial hysteresis loops, a networked collaboration of
computer systems may be implemented. In other words, several networked computers
could run the same optimization routines simultaneously, sharing all the best partial
solutions found, then incorporating them into their own populations and continue
their search. Thus, the time needed to get optimal solutions may be dramatically
reduced.
96
REFERENCES
1
Renji, K., and Shankar Narayan, S. Loss Factors of Composite Honeycomb Sandwich
Panels, J. of Sound and Vib., 250(4), 745-761, (2002).
2
Hornig, K.H., and Flowers, G.T. Parameter Characterisation of the Bouc/Wen
Mechanical Hysteresis Model for Sandwich Composite Materials using Real Coded
Genetic Algorithms, Intl. J. of Acoustics and Vib. (IJAV), 10(2), 73-81, (2005).
3
Nilsson, E. Some Dynamic Properties of Honeycomb Panels, Report, Royal Inst. of
Tech., Stockholm, Sweden, (2000).
4
Chandra, R., Singh, S.P., and Gupta, K. Damping Studies in Fiber-Reinforced
Composites - a Review, Composite Structures, 46(1), 41-51, (1999).
5
Sain, P.M., Sain, M.K., and Spencer, B.F. Models of Hysteresis and Application to
Structural Control, Proc. of the Amer. Control Conf., Albuquerque, NM, 16-20,
(1997).
6
Krasnosel?skii, M.A., and Pokorovskii, A.V. Systems with Hysteresis, Springer-
Verlag, New York, (1989).
7
Chua, L.O., and Bass, S.C. A Generalized Hysteresis Model, IEEE Trans. Circuit
Theory, 19, 36-48, (1972).
8
Preisach, P. ?ber die Magnetische Nachwirkung, Zeitschrift f?r Physik, 94, 277-302,
(1938).
9
Bouc, R. Forced Vibration of Mechanical Systems with Hysteresis, Proc. of the 4th
Conf. on Nonlinear Oscillation, Prague, Czechoslovakia, (1967).
10
Wen, Y.K. Method for Random Vibration of Hysteretic Systems, J. Engr. Mech. 102,
249-263, (1976).
97
11
Spencer, B.F. Reliability of Randomly Excited Hysteretic Structures, Springer-
Heidelberg, New York, (1986).
12
Ni, Y.Q., Ko, J.M., and Wong, C.W. Identification of Non-Linear Hysteretic Isolators
form Periodic Vibration Tests, J. of Sound and Vib., 217 (4), 737-756 (1998).
13
Constantinou, M.C., and Tadjbakhsh, I.G. Hysteretic Dampers in Base Isolation:
Random Approach, J. of Struct. Eng., ASCE 111 (ST4), 705-721, (1998).
14
Smyth, A.W., Masri, S.F., Kosmatopoulos, E.B., Chassiakos, A.G., and Caughey,
T.K. Development of Adaptive Modeling Techniques for Non-linear Hysteretic
Systems, Intl. J. of Non-Linear Mech., 37(8), 1435-1451, (2002).
15
Heine, C.P. Simulated Response of Degrading Hysteretic Joints With Slack Behavior,
Ph.D. Dissertation, Virginia Polytech. Inst. and State Univ., (2001).
16
Macioce, P. Viscoelastic Damping 101, Sound and Vibration Magazine, 4, 4-5,
(2003).
17
Meirovitch, L. Elements of Vibration Analysis, McGraw-Hill, New York, (1986).
18
Sun, C.T., and Lu, Y.P. Vibration Damping of Structural Elements, Prentice Hall,
New Jersey, (1995).
19
Zener, C. Elasticity and Inelasticity of Metals, University of Chicago Press, Chicago,
(1948).
20
Crandall, S.H. The Role of Damping in Vibration Theory, J. of Sound and Vib., 11(1),
3-18, (1970).
21
Lazan, B.J. Damping of Materials and Members in Structural Mechanics, Pergamon
Press, UK, (1968).
22
Hooker, R.J. On the Interpretation of Biaxial Hysteresis Loops, J. of Sound and Vib.,
76(3), 463-466, (1981).
23
So, C.K., Lai, T.C., and Tse, P.C. The Measurement of Material Damping by Free-
Vibration Technique with Periodic Excitation, Experimental Techniques, 14(3), 41-
42, (1990).
98
24
McDaniel, J.G., Dupont, P., and Salvino, L.A. Wave Approach to Estimating
Frequency-Dependent Damping under Transient Loading, J. of Sound and Vib.,
231(2), 433-449, (2000).
25
Brown, K.T., and Norton, M.P. Some Comments on the Experimental Determination
of Modal Densities and Loss Factors for Statistical Energy Analysis Applications, J.
of Sound and Vib., 102, 588-594, (1985).
26
Boyaci, H. Vibrations of Stretched Damped Beams under Non-Ideal Boundary
Conditions, S?dhan?, 31(1), 1-8, (2006).
27
Qiao, S.L., Pilipchuk, V.N., and Ibrahim, R.A. Modeling and Simulation of Elastic
Structures with Parameter Uncertainties and Relaxation of Joints, J. of Vib. and
Acoustics, 123, 45-52, (2001).
28
Liu, G.R., Han, X., and Lam, K.Y. A Combined Genetic Algorithm and Nonlinear
Least Squares Method for Material Characterization using Elastic Waves, Comput.
Methods in Appl. Mech. and Eng., 191(17), 1909-1921, (2002).
29
Matsuoka, T., Yamamoto, S., and Takahara, M. Prediction of Structures and
Mechanical Properties of Composites using a Genetic Algorithm and Finite Element
Method, J. of Materials Science, 36, 27-33, (2001).
30
Mac, A., and Gurba, W. Genetic Algorithms and Finite Element Analysis in
Optimization of Composite Structures, Composite Structures, 54, 275-281, (2001).
31
Gantovnik, V.B., G?rdal, Z., and Watson, L.T. A Genetic Algorithm with Memory
for Optimal Design of Laminated Sandwich Composite Panels, Composite Structures,
58, 513-520, (2002).
32
Kovacs, G., Groenwold, A.A., Jarmai, K., and Farkas, J. Analysis and Optimum
Design of Fibre-Reinforced Composite Structures, Struct. and Multidisc.
Optimization, 28(2), 170-179, (2004).
33
Sekishiro, M., and Todoroki, A. Optimizations of Stiffened Composite Panel Using
Fractal Branch and Bound Method, Proc. of the 46th Struct. Dynamics and Materials
Conf., Austin, TX, 1525-1534, (2005).
99
34
Kyprianou, A., and Worden, K. Identification of Hysteretic Systems using the
Differential Evolution Algorithm, Proc. of the 23rd Intl. Conf. on Noise and
Vibration Engineering (ISMA), Leuven, Belgium, 557-568, (1998).
35
Zhang, X., Huang, Y., Liu, Y., Wang, X., and Gao, F. Method Identifying the
Parameters of Bouc-Wen Hysteretic Nonlinear Model Based on Genetic-Algorithm,
Proc. of the IEEE Intl. Conf. on Intelligent Proc. Systems (ICIPS), Beijing, China,
602-605, (1997).
36
Peng, J.Y., Li, H., and Suzuki, Y. Modeling of Nonlinear Hysteresis with Pinching,
Journal of Shenyang Jianzhu University, 21(4), 325-328, (2005).
37
Fu, J., Hou, C., and Dou, L. Application of Particle Swarm Optimization Algorithm
in the Design of Multilayered Planar Shielding Body, Chinese Journal of Electronics,
14(2), 357-360, (2005).
38
Chen, J., Ge, R., and Wei, J. Optimization of the Reliability of Laminated Plates
Based on the PSO, J. of Huazhong Univ. of Sc. and Technology, 34(4), 96-98, (2006).
39
Holland, J.H. Adaptation in Natural and Artificial Systems, Ann Arbor, Michigan,
The University of Michigan Press, (1975).
40
Goldberg, D.E. Genetic Algorithms in Search, Optimization, and Machine Learning,
Reading, MA, Addison-Wesley, (1989).
41
Haupt, R.L., and Haupt, S.E. Practical Genetic Algorithms, 2nd Ed., Wiley, New
York, (2004).
42
Zilouchian, A., and Jamshidi, M. (Editors). Intelligent Control Systems Using Soft
Computing Methodologies, CRC Press, USA, (2001).
43
Michalewicz, Z. Genetic Algorithms + Data Structures = Evolution Programs, 3rd
Ed., Springer-Verlag, New York, (1996).
44
Kennedy, J., and Eberhart, R. Particle Swarm Optimization, Proc. of the IEEE Intl.
Conf. on Neural Networks, Perth, Australia: IEEE Service Center, Piscataway, NJ,
1942-1948, (1995).
100
45
Shi, Y.H., and Eberhart, R.C. Parameter Selection in Particle Swarm Optimization,
Proc. of the 7th Annual Conf. on Evolutionary Programming, New York, NY, 591-
600, (1998).
46
Clerc, M. The Swarm and the Queen: Towards a Deterministic and Adaptive Particle
Swarm Optimization, Proc. of the 1999 ICEC, Washington, DC, 1951-1957, (1999).
47
Carlisle, A., and Dozier, G. An Off-The-Shelf PSO, Proc. of the 2001 Workshop on
Particle Swarm Optimization, Indianapolis, IN, 1-6, (2001).
101
APPENDIX A
PROGRAMS CREATED IN MATLAB R2006a FOR THE IMPLEMENTATION
OF GENETIC ALGORITHMS (GA) AND PARTICLE SWARM
OPTIMIZATION (PSO) FOR SYSTEM IDENTIFICATION OF
SANDWICH COMPOSITE MATERIALS
102
%Optimization Program for the Hysteresis model using GA
%(c)2006 - Klaus H. Hornig
%____________________________________________________________
clear all; clc;
%Parameters: [xi0, S_xi, w0, Sw]
pb=[0,.1;-10,10;100,500;-10000,10000]; %parameter bounds
np=size(pb,1); %np=number of parameters
tspan=.4*ones(1,120); %time span
Meq=36.8e-3; %equivalent mass (sample #2)
%GA operators
ni=18; %Number of iterations of the whole GA code (main loop)
nipop=6000; %Number of initial population
rempop=500; %Number of remaining population
alphaBLX=0.3; %alpha parameter for X-over
mutrate=.05; %Mutation rate
mutamount=.3; %Gaussian mutation amount
addpopmut=.2; %Additional population mutation amount
addpoplevel=4; %Additional population "add-on" level (L)
interppts=100; %Resampling resolution (initial interpolation)
randn('state',0);
rand('state',0);
%Generate the initial population
ipop=genipop(pb,nipop);
ipop=[ipop,zeros(nipop,1)];
%load experimental loops (files that contain both variables, HexpT and HexpB,
%for TOP and BOTTOM parts, respectively)
hysloopsnames=['SAMPLE3_02mm_370Hz';... %0.2mm
'SAMPLE2_02mm_380Hz';'SAMPLE2_02mm_385Hz';'SAMPLE2_02mm_390Hz';...
'SAMPLE2_02mm_395Hz';'SAMPLE2_02mm_400Hz';'SAMPLE2_02mm_402Hz';...
'SAMPLE2_02mm_404Hz';'SAMPLE2_02mm_406Hz';'SAMPLE2_02mm_408Hz';...
'SAMPLE2_02mm_410Hz';'SAMPLE2_02mm_412Hz';'SAMPLE2_02mm_414Hz';...
'SAMPLE2_02mm_415Hz';'SAMPLE2_02mm_416Hz';'SAMPLE2_02mm_417Hz';...
'SAMPLE2_02mm_418Hz';'SAMPLE2_02mm_419Hz';'SAMPLE2_02mm_420Hz';...
'SAMPLE2_02mm_421Hz';'SAMPLE2_02mm_422Hz';'SAMPLE2_02mm_423Hz';...
'SAMPLE2_02mm_424Hz';'SAMPLE2_02mm_425Hz';'SAMPLE2_02mm_426Hz';...
'SAMPLE2_02mm_427Hz';'SAMPLE2_02mm_428Hz';'SAMPLE2_02mm_429Hz';...
'SAMPLE2_02mm_430Hz';'SAMPLE2_02mm_432Hz';'SAMPLE2_02mm_434Hz';...
'SAMPLE2_02mm_436Hz';'SAMPLE2_02mm_438Hz';'SAMPLE2_02mm_440Hz';...
'SAMPLE2_02mm_442Hz';'SAMPLE2_02mm_444Hz';'SAMPLE2_02mm_450Hz';...
'SAMPLE2_02mm_455Hz';'SAMPLE2_02mm_460Hz';'SAMPLE2_02mm_470Hz';...
'SAMPLE2_05mm_370Hz';... %0.5mm
'SAMPLE2_05mm_380Hz';'SAMPLE2_05mm_385Hz';'SAMPLE2_05mm_390Hz';...
'SAMPLE2_05mm_395Hz';'SAMPLE2_05mm_400Hz';'SAMPLE2_05mm_402Hz';...
'SAMPLE2_05mm_404Hz';'SAMPLE2_05mm_406Hz';'SAMPLE2_05mm_408Hz';...
'SAMPLE2_05mm_410Hz';'SAMPLE2_05mm_412Hz';'SAMPLE2_05mm_414Hz';...
'SAMPLE2_05mm_415Hz';'SAMPLE2_05mm_416Hz';'SAMPLE2_05mm_417Hz';...
'SAMPLE2_05mm_418Hz';'SAMPLE2_05mm_419Hz';'SAMPLE2_05mm_420Hz';...
'SAMPLE2_05mm_421Hz';'SAMPLE2_05mm_422Hz';'SAMPLE2_05mm_423Hz';...
'SAMPLE2_05mm_424Hz';'SAMPLE2_05mm_425Hz';'SAMPLE2_05mm_426Hz';...
'SAMPLE2_05mm_427Hz';'SAMPLE2_05mm_428Hz';'SAMPLE2_05mm_429Hz';...
'SAMPLE2_05mm_430Hz';'SAMPLE2_05mm_432Hz';'SAMPLE2_05mm_434Hz';...
'SAMPLE2_05mm_436Hz';'SAMPLE2_05mm_438Hz';'SAMPLE2_05mm_440Hz';...
103
'SAMPLE2_05mm_442Hz';'SAMPLE2_05mm_444Hz';'SAMPLE2_05mm_450Hz';...
'SAMPLE2_05mm_455Hz';'SAMPLE2_05mm_460Hz';'SAMPLE2_05mm_470Hz';...
'SAMPLE2_08mm_370Hz';... %0.8mm
'SAMPLE2_08mm_380Hz';'SAMPLE2_08mm_385Hz';'SAMPLE2_08mm_390Hz';...
'SAMPLE2_08mm_395Hz';'SAMPLE2_08mm_400Hz';'SAMPLE2_08mm_402Hz';...
'SAMPLE2_08mm_404Hz';'SAMPLE2_08mm_406Hz';'SAMPLE2_08mm_408Hz';...
'SAMPLE2_08mm_410Hz';'SAMPLE2_08mm_412Hz';'SAMPLE2_08mm_414Hz';...
'SAMPLE2_08mm_415Hz';'SAMPLE2_08mm_416Hz';'SAMPLE2_08mm_417Hz';...
'SAMPLE2_08mm_418Hz';'SAMPLE2_08mm_419Hz';'SAMPLE2_08mm_420Hz';...
'SAMPLE2_08mm_421Hz';'SAMPLE2_08mm_422Hz';'SAMPLE2_08mm_423Hz';...
'SAMPLE2_08mm_424Hz';'SAMPLE2_08mm_425Hz';'SAMPLE2_08mm_426Hz';...
'SAMPLE2_08mm_427Hz';'SAMPLE2_08mm_428Hz';'SAMPLE2_08mm_429Hz';...
'SAMPLE2_08mm_430Hz';'SAMPLE2_08mm_432Hz';'SAMPLE2_08mm_434Hz';...
'SAMPLE2_08mm_436Hz';'SAMPLE2_08mm_438Hz';'SAMPLE2_08mm_440Hz';...
'SAMPLE2_08mm_442Hz';'SAMPLE2_08mm_444Hz';'SAMPLE2_08mm_450Hz';...
'SAMPLE2_08mm_455Hz';'SAMPLE2_08mm_460Hz';'SAMPLE2_08mm_470Hz'];
load efSample2; %load variable "ef", with the excitation frequency of each loop
numloops=length(ef);
userweights=ones(1,numloops);
for loopsload=1:numloops %create cell array, storing all loops {TOP,BOTTOM}
load(hysloopsnames(loopsload,:));
[HexpT,HexpB]=reinterpolate(HexpT,HexpB,interppts); %reinterpolate loops
maxF(loopsload)=max(abs([HexpT(:,2);HexpB(:,2)]));
rnT=(maxF(loopsload)*noiselevel/100)*randn(size(HexpT,1),1);
rnB=(maxF(loopsload)*noiselevel/100)*randn(size(HexpB,1),1);
HexpT=[HexpT(:,1),HexpT(:,2)+rnT];
HexpB=[HexpB(:,1),HexpB(:,2)+rnB];
hysteresis_loops{loopsload,1}=HexpT;
hysteresis_loops{loopsload,2}=HexpB;
end
maxFall=max(maxF);
for i=1:numloops %get peak displacement for each loop
Xpk(i)=max(abs([hysteresis_loops{i,1}(:,1);hysteresis_loops{i,2}(:,1)]));
amplitudescaling(i)=maxFall/maxF(i);
end
loopsweights=userweights.*amplitudescaling;
nap=0;
for lvl=1:addpoplevel %calculate number of additional population (AddPopNumber)
nap=nap+nchoosek(np,lvl);
end
%CALCULATE FITNESS FOR EACH CHROMOSOME
for i=1:nipop
fprintf('Calculating fitness for chromosome N?: %g\n',i);
p=ipop(i,1:np);
lasterr='a';
while lasterr~=' '
lasterr=' ';
try
ipop(i,np+1)=0;
for j=1:numloops
HexpT=hysteresis_loops{j,1};
104
HexpB=hysteresis_loops{j,2};
[NOS,dummy]=size([HexpT;HexpB]); %NOS=Number Of Samples
ipop(i,np+1)=ipop(i,np+1)+loopsweights(j)*...
erreval(HexpT,HexpB,p,ef(j),Xpk(j),Meq,tspan(j))/NOS;
ipoperror=ipop(i,np+1);
end %for j
fprintf('%g\n\n',ipop(i,np+1));
catch
fprintf('\nCouldn''t calculate. Generating new chromosome\n\n');
ipop(i,1:np)=genipop(pb,1);
p=ipop(i,1:np);
lasterr='a';
end %try-catch
end %while
end %for i
fprintf('\nKeeping the best %g chromosomes...\n\n\n',rempop);
iipop=sortrows(ipop,np+1);
ipop=iipop(1:rempop,:); %keep just the best 'rempop' ones
BC=ipop(1,1:np); %Best Chromosome
BCf=ipop(1,np+1); %Best Chromosome fitness
fprintf('The best chromosome so far is:\n---------------------------\n');
bp_string_values=sprintf('[ %g ',BC(1));
if np>1
for i=2:np
bp_string_values=[bp_string_values,sprintf(' , %g',BC(i))];
end
end
fprintf([bp_string_values,' ]\n\n\n']);
%Fitness history value
fhv=zeros(1,ni+1);
fhv(1)=BCf;
save allvariables;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% MAIN LOOP %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for loop=1:ni
%Calculate mating probability distribution, according to RANKING
cumprob=zeros(rempop,1);
prob=((rempop-(1:rempop)+1)/rempop)*2/(rempop+1);
for i=1:rempop
cumprob(i,1)=sum(prob(1:i));
end
%CHOOSE THE MATES
randmates=rand(rempop,1);
mates=zeros(rempop,1); %initialize variable "mates"
for i=1:rempop
j=1; %initialize variable 'j'
while randmates(i)>cumprob(j)
j=j+1;
end
mates(i)=j; %it's the entry number of the parent
105
end
%check to avoid a parent to mate with itself
for i=1:rempop/2 %because there are rempop/2 couples
while mates(2*i-1)==mates(2*i)
temprandmate=rand;
for ii=1:rempop
jj=1; %initialize variable 'jj'
while temprandmate>cumprob(jj)
jj=jj+1;
end
mates(2*i)=jj;
end %for ii
end %while mates
end %for i
%MATE THE PARENTS
offspring=zeros(rempop,np); %initialize variable "offspring"
for i=1:(rempop/2) %couple number
crosspnt=ceil(np*rand); %CROSSOVER POINT!!!
for parameter=1:np
if parameterpb(parameter,2)
offspring(2*i-1,parameter)=ipop(mates(2*i-1),parameter);
end
offspring(2*i,parameter)=ipop(mates(2*i),parameter)...
+temprand*(ipop(mates(2*i-1),parameter)...
-ipop(mates(2*i),parameter));
%check if out of bounds
if offspring(2*i,parameter)pb(parameter,2)
offspring(2*i,parameter)=ipop(mates(2*i),parameter);
end
end
if parameter>crosspnt
offspring(2*i-1,parameter)=ipop(mates(2*i),parameter);
offspring(2*i,parameter)=ipop(mates(2*i-1),parameter);
end
end %for parameter
end %for i
%Create a new population of (Gaussian) mutated versions of BC
APcount=0; %additional population count
addpop=ones(nap,1)*BC;
for j=1:addpoplevel
C=nchoosek(1:np,j)';
for ap=1:size(C,2)
APcount=APcount+1;
for i=C(:,ap)'
newvalue=addpop(APcount,i)+randn*addpopmut*(pb(i,2)-pb(i,1));
106
if newvalue>pb(i,1)&&newvaluepb(rpc,1)&&newvalue1
for i=2:np
bp_string_values=[bp_string_values,sprintf(' , %g',BC(i))];
end
end
fprintf([bp_string_values,' ]\n\n\n']);
%Fitness history value
fhv(loop+1)=BCf;
save allvariables;
%**************************************************************
end % END OF MAIN LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!*
%**************************************************************
108
%Optimization Program for the Hysteresis model using PSO
%(c)2006 - Klaus H. Hornig
%____________________________________________________________
clear all; clc;
%Parameters: [xi0, S_xi, w0, Sw]
pb=[0,.1;-10,10;100,500;-10000,10000]; %parameter bounds (for constant damping)
np=size(pb,1); %np=number of parameters
tspan=.4*ones(1,120);
Meq=36.8e-3; %equivalent mass (sample #2)
%GA operators constants
ni=18; %Number of iterations of the whole GA code (main loop)
nipop=6000; %Number of initial population
rempop=500; %Number of remaining population
interppts=100; %resampling resolution (initial interpolation)
addpopmut=.2; %Additional population mutation amount
addpoplevel=4; %max. number of parameters to be changed in addpop
randn('state',0);
rand('state',0);
%Initialize PSO variables
Pi=zeros(nipop,np);
Pif=zeros(nipop,1);
Vi=zeros(rempop,np);
phi1=2.05; %PSO cognition factor
phi2=2.05; %PSO social factor
phi=phi1+phi2;
K=2/abs(2-phi-sqrt(phi^2-4*phi)); %PSO constriction factor
%Generate the initial population
ipop=genipop(pb,nipop);
ipop=[ipop,zeros(nipop,1)];
%load experimental loops (files that contain both variables, HexpT and HexpB,
%for TOP and BOTTOM parts, respectively)
hysloopsnames=['SAMPLE3_02mm_370Hz';... %0.2mm
'SAMPLE2_02mm_380Hz';'SAMPLE2_02mm_385Hz';'SAMPLE2_02mm_390Hz';...
'SAMPLE2_02mm_395Hz';'SAMPLE2_02mm_400Hz';'SAMPLE2_02mm_402Hz';...
'SAMPLE2_02mm_404Hz';'SAMPLE2_02mm_406Hz';'SAMPLE2_02mm_408Hz';...
'SAMPLE2_02mm_410Hz';'SAMPLE2_02mm_412Hz';'SAMPLE2_02mm_414Hz';...
'SAMPLE2_02mm_415Hz';'SAMPLE2_02mm_416Hz';'SAMPLE2_02mm_417Hz';...
'SAMPLE2_02mm_418Hz';'SAMPLE2_02mm_419Hz';'SAMPLE2_02mm_420Hz';...
'SAMPLE2_02mm_421Hz';'SAMPLE2_02mm_422Hz';'SAMPLE2_02mm_423Hz';...
'SAMPLE2_02mm_424Hz';'SAMPLE2_02mm_425Hz';'SAMPLE2_02mm_426Hz';...
'SAMPLE2_02mm_427Hz';'SAMPLE2_02mm_428Hz';'SAMPLE2_02mm_429Hz';...
'SAMPLE2_02mm_430Hz';'SAMPLE2_02mm_432Hz';'SAMPLE2_02mm_434Hz';...
'SAMPLE2_02mm_436Hz';'SAMPLE2_02mm_438Hz';'SAMPLE2_02mm_440Hz';...
'SAMPLE2_02mm_442Hz';'SAMPLE2_02mm_444Hz';'SAMPLE2_02mm_450Hz';...
'SAMPLE2_02mm_455Hz';'SAMPLE2_02mm_460Hz';'SAMPLE2_02mm_470Hz';...
'SAMPLE2_05mm_370Hz';... %0.5mm
'SAMPLE2_05mm_380Hz';'SAMPLE2_05mm_385Hz';'SAMPLE2_05mm_390Hz';...
'SAMPLE2_05mm_395Hz';'SAMPLE2_05mm_400Hz';'SAMPLE2_05mm_402Hz';...
'SAMPLE2_05mm_404Hz';'SAMPLE2_05mm_406Hz';'SAMPLE2_05mm_408Hz';...
'SAMPLE2_05mm_410Hz';'SAMPLE2_05mm_412Hz';'SAMPLE2_05mm_414Hz';...
'SAMPLE2_05mm_415Hz';'SAMPLE2_05mm_416Hz';'SAMPLE2_05mm_417Hz';...
109
'SAMPLE2_05mm_418Hz';'SAMPLE2_05mm_419Hz';'SAMPLE2_05mm_420Hz';...
'SAMPLE2_05mm_421Hz';'SAMPLE2_05mm_422Hz';'SAMPLE2_05mm_423Hz';...
'SAMPLE2_05mm_424Hz';'SAMPLE2_05mm_425Hz';'SAMPLE2_05mm_426Hz';...
'SAMPLE2_05mm_427Hz';'SAMPLE2_05mm_428Hz';'SAMPLE2_05mm_429Hz';...
'SAMPLE2_05mm_430Hz';'SAMPLE2_05mm_432Hz';'SAMPLE2_05mm_434Hz';...
'SAMPLE2_05mm_436Hz';'SAMPLE2_05mm_438Hz';'SAMPLE2_05mm_440Hz';...
'SAMPLE2_05mm_442Hz';'SAMPLE2_05mm_444Hz';'SAMPLE2_05mm_450Hz';...
'SAMPLE2_05mm_455Hz';'SAMPLE2_05mm_460Hz';'SAMPLE2_05mm_470Hz';...
'SAMPLE2_08mm_370Hz';... %0.8mm
'SAMPLE2_08mm_380Hz';'SAMPLE2_08mm_385Hz';'SAMPLE2_08mm_390Hz';...
'SAMPLE2_08mm_395Hz';'SAMPLE2_08mm_400Hz';'SAMPLE2_08mm_402Hz';...
'SAMPLE2_08mm_404Hz';'SAMPLE2_08mm_406Hz';'SAMPLE2_08mm_408Hz';...
'SAMPLE2_08mm_410Hz';'SAMPLE2_08mm_412Hz';'SAMPLE2_08mm_414Hz';...
'SAMPLE2_08mm_415Hz';'SAMPLE2_08mm_416Hz';'SAMPLE2_08mm_417Hz';...
'SAMPLE2_08mm_418Hz';'SAMPLE2_08mm_419Hz';'SAMPLE2_08mm_420Hz';...
'SAMPLE2_08mm_421Hz';'SAMPLE2_08mm_422Hz';'SAMPLE2_08mm_423Hz';...
'SAMPLE2_08mm_424Hz';'SAMPLE2_08mm_425Hz';'SAMPLE2_08mm_426Hz';...
'SAMPLE2_08mm_427Hz';'SAMPLE2_08mm_428Hz';'SAMPLE2_08mm_429Hz';...
'SAMPLE2_08mm_430Hz';'SAMPLE2_08mm_432Hz';'SAMPLE2_08mm_434Hz';...
'SAMPLE2_08mm_436Hz';'SAMPLE2_08mm_438Hz';'SAMPLE2_08mm_440Hz';...
'SAMPLE2_08mm_442Hz';'SAMPLE2_08mm_444Hz';'SAMPLE2_08mm_450Hz';...
'SAMPLE2_08mm_455Hz';'SAMPLE2_08mm_460Hz';'SAMPLE2_08mm_470Hz'];
load efSample2; %load variable "ef", with the excitation frequency of each loop
numloops=length(ef);
userweights=ones(1,numloops);
for loopsload=1:numloops %create cell array, storing all loops {TOP,BOTTOM}
load(hysloopsnames(loopsload,:));
[HexpT,HexpB]=reinterpolate(HexpT,HexpB,interppts); %reinterpolate loops
maxF(loopsload)=max(abs([HexpT(:,2);HexpB(:,2)]));
rnT=(maxF(loopsload)*noiselevel/100)*randn(size(HexpT,1),1);
rnB=(maxF(loopsload)*noiselevel/100)*randn(size(HexpB,1),1);
HexpT=[HexpT(:,1),HexpT(:,2)+rnT];
HexpB=[HexpB(:,1),HexpB(:,2)+rnB];
hysteresis_loops{loopsload,1}=HexpT;
hysteresis_loops{loopsload,2}=HexpB;
end
maxFall=max(maxF);
for i=1:numloops %get peak displacement for each loop
Xpk(i)=max(abs([hysteresis_loops{i,1}(:,1);hysteresis_loops{i,2}(:,1)]));
amplitudescaling(i)=maxFall/maxF(i);
end
loopsweights=userweights.*amplitudescaling;
nap=0;
for lvl=1:addpoplevel %calculate number of additional population (AddPopNumber)
nap=nap+nchoosek(np,lvl);
end
%CALCULATE FITNESS FOR EACH CHROMOSOME
for i=1:nipop
fprintf('Calculating fitness for particle N?: %g\n',i);
p=ipop(i,1:np);
lasterr='a';
while lasterr~=' '
110
lasterr=' ';
try
ipop(i,np+1)=0;
for j=1:numloops
HexpT=hysteresis_loops{j,1};
HexpB=hysteresis_loops{j,2};
NOS=size([HexpT;HexpB],1); %NOS=Number Of Samples
ipop(i,np+1)=ipop(i,np+1)+loopsweights(j)*...
erreval(HexpT,HexpB,p,ef(j),Xpk(j),Meq,tspan(j))/NOS;
ipoperror=ipop(i,np+1);
end %for j
fprintf('%g\n\n',ipop(i,np+1));
catch
fprintf('\nCouldn''t calculate. Generating new particle\n\n');
ipop(i,1:np)=genipop(pb,1);
p=ipop(i,1:np);
lasterr='a';
end %try-catch
end %while
end %for i
fprintf('\nKeeping the best %g particles...\n\n\n',rempop);
iipop=sortrows(ipop,np+1);
ipop=iipop(1:rempop,:); %keep just the best 'rempop' ones
Pi=ipop(:,1:np);
Pif=ipop(:,np+1);
%Find BEST Pi (least error)
idx=min(find(Pif==min(Pif)));
Pg=Pi(idx,:);
Pgf=Pif(idx);
fprintf('The best particle so far is:\n---------------------------\n');
bp_string_values=sprintf('[ %g ',Pg(1));
if np>1
for i=2:np
bp_string_values=[bp_string_values,sprintf(' , %g',Pg(i))];
end
end
fprintf([bp_string_values,' ]\n\n\n']);
%Fitness history value
fhv=zeros(1,ni+1);
fhv(1)=Pgf;
save allvariables;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% MAIN LOOP %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for loop=1:ni
for par=1:rempop %particle index
rnd1=rand(1,np);
rnd2=rand(1,np);
Vi(par,:)=K*(Vi(par,:)+phi1*rnd1.*(Pi(par,:)-ipop(par,1:np))+...
111
phi2*rnd2.*(Pg-ipop(par,1:np))); %velocity
xt=ipop(par,1:np)+Vi(par,:); % temporal x-positions
for i=1:np %update particle position
if (xt(i)>pb(i,1)) && (xt(i)pb(i,1)&&newvalue1
for i=2:np
bp_string_values=[bp_string_values,sprintf(' , %g',Pg(i))];
end
end
fprintf([bp_string_values,' ]\n\n\n']);
%Fitness history value
fhv(loop+1)=Pgf;
save allvariables;
end % END OF MAIN LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!
113
function ipop=genipop(pb,nipop)
%genipop(pb,nipop):
%-----------------
%
% Generates the initial population.
%
% where: pb = parameters bounds vector
% nipop = number of individuals to be created
np=size(pb,1); %np=number of parameters
ipop=zeros(nipop,np); %initialize ipop
for parameter=1:np
ipop(:,parameter)=(pb(parameter,2)-pb(parameter,1))...
*rand(nipop,1)+pb(parameter,1);
end
function fitness=erreval(HexpT,HexpB,p,ef,Xpk,Meq,tspan)
% ERREVAL function to evaluate the error of a hysteresis loop
% from the hysteresis model compared to experimental data.
%
% Syntax: erreval(HexpT,HexpB,p,ef,Xpk,Meq,tspan)
%
% where: HexpT = experimental input data (TOP part of curve)
% HexpB = experimental input data (BOTTOM part of curve)
% p = chromosome (set of 'np' parameters)
% ef = excitation frequency [Hz]
% Xpk = peak displacement amplitude [m]
% Meq = 1-DOF equivalent mass of the beam (kg)
% tspan = time span of numerical integration [sec]
%
%(c)2006 - Klaus H. Hornig
w=2*pi*ef;
t=0:1e-4:tspan;
t=t';
xi=p(1)+p(2)*Xpk;
wn=p(3)+p(4)*Xpk;
X=Xpk*sin(w*t);
Fcf=Xpk*Meq*((wn^2-w^2)*sin(w*t)+2*xi*wn*w*cos(w*t));
hys1=[X,Fcf];
[modelT,modelB]=getloopdisp(hys1);
YmodT=interp1(modelT(:,1),modelT(:,2),HexpT(:,1),'linear'); %interp. TOP data
%Remove Nan's from YmodT
YmodTl=length(YmodT);
if isnan(YmodT(1))
YmodT=YmodT(2:YmodTl);
HexpT=HexpT(2:YmodTl,:);
end
YmodTl=length(YmodT);
if isnan(YmodT(YmodTl))
YmodT=YmodT(1:YmodTl-1);
114
HexpT=HexpT(1:YmodTl-1,:);
end
YmodB=interp1(modelB(:,1),modelB(:,2),HexpB(:,1),'linear'); %interp. BOTTOM data
%Remove Nan's from YmodB
YmodBl=length(YmodB);
if isnan(YmodB(1))
YmodB=YmodB(2:YmodBl);
HexpB=HexpB(2:YmodBl,:);
end
YmodBl=length(YmodB);
if isnan(YmodB(YmodBl))
YmodB=YmodB(1:YmodBl-1);
HexpB=HexpB(1:YmodBl-1,:);
end
fitness=leastsqdisp([HexpT;HexpB],[YmodT;YmodB]); %get the OLS error
function value=leastsqdisp(A,B)
[NRA,NCA]=size(A);
[NRB,NCB]=size(B);
if NRA==NRB
valuevector=(A(:,2)-B).^2;
value=sum(valuevector);
else
disp('Least-squares error: both matrices must have the same dimensions');
end
function [top,bot]=getloopdisp(hys1)
[NR,NC]=size(hys1);
i=0;
while hys1(NR-i,1)==hys1(NR-i-1,1)
i=i+1;
end
%_______________________________________
while hys1(NR-i,1)hys1(NR-i-1,1)
115
i=i+1;
end
F2=NR-i; %2nd. end
%advance if next is the same
F2f=0;
while hys1(NR-i,1)==hys1(NR-i-1,1)
i=i+1;
F2f=F2f+1;
end
%_______________________________________
while hys1(NR-i,1)hys1(NR-i-1,1)
i=i+1;
end
F4=NR-i; %4th. end
if (hys1(NR,1) positive slope loop / end bottom
TOPmax=F1-F1f;
BOTmin=F3;
BOTmax=F2-F2f;
elseif (hys1(NR,1)>hys1(NR-1,1)) & (hys1(NR,2) negative slope loop / end top
TOPmax=F3-F3f;
BOTmin=F3;
BOTmax=F2-F2f;
elseif (hys1(NR,1)>hys1(NR-1,1)) & (hys1(NR,2)>hys1(NR-1,2))
TOPmin=F4; %=> positive slope loop / end top
TOPmax=F3-F3f;
BOTmin=F3;
BOTmax=F2-F2f;
elseif (hys1(NR,1)hys1(NR-1,2))
TOPmin=F2; %=> negative slope loop / end bottom
TOPmax=F1-F1f;
BOTmin=F3;
BOTmax=F2-F2f;
end
top=hys1(TOPmin:TOPmax,:);
bot=hys1(BOTmin:BOTmax,:);
116
function [HexpT,HexpB]=reinterpolate(HexpTo,HexpBo,ints,plotflag,inttype)
%
%Program to reinterpolate the trial hysteresis curves, adding
%a specified discretization of the displacement axis.
%
%HexpTo : original trial TOP curve
%HexpBo : original trial BOTTOM curve
%ints : interpolation size
%plotflag : comparison plot flag (1=yes)
%inttype : interpolation type
% (nearest,linear,spline,pchip,cubic,v5cubic)
if nargin==2
ints=100;
plotflag=0;
inttype='linear';
elseif nargin==3
plotflag=0;
inttype='linear';
elseif nargin==4
inttype='linear';
end
%TOP curve
Xmin=min(HexpTo(:,1));
Xmax=max(HexpTo(:,1));
X=linspace(Xmin,Xmax,ints);
X=X';
F=interp1(HexpTo(:,1),HexpTo(:,2),X,inttype,'extrap');
HexpT=[X,F];
%BOTTOM curve
Xmin=min(HexpBo(:,1));
Xmax=max(HexpBo(:,1));
X=linspace(Xmin,Xmax,ints);
X=X';
F=interp1(HexpBo(:,1),HexpBo(:,2),X,inttype,'extrap');
HexpB=[X,F];
if plotflag
%PLOT comparison
figure %TOP
plot(HexpTo(:,1),HexpTo(:,2))
hold on
plot(HexpTo(:,1),HexpTo(:,2),'kx','MarkerSize',14)
plot(HexpT(:,1),HexpT(:,2),'r')
plot(HexpT(:,1),HexpT(:,2),'+g')
figure %BOTTOM
plot(HexpBo(:,1),HexpBo(:,2))
hold on
plot(HexpBo(:,1),HexpBo(:,2),'kx','MarkerSize',14)
plot(HexpB(:,1),HexpB(:,2),'r')
plot(HexpB(:,1),HexpB(:,2),'+g')
end
117
%PrepareExpData.m
%Program to create files with HexpT and HexpB variables (which are the TOP and
%BOTTOM parts of each hysteresis loop) from experimental results
clear all; clc;
%Sample #2
freqs=[370,380,385,390,395,400,402,404,406,408,410,412,414,415,416,417,418,...
419,420,421,422,423,424,425,426,427,428,429,430,432,434,436,438,440,...
442,444,450,455,460,470];
j=0;
for i=1:40
j=j+1;
file=sprintf('TRAC%g',i); %import TRAC***.txt files (experimental data)
eval(['load ',file,'.lvm;']);
eval(['trac=',file,';']);
trac(:,1)=trac(:,1)*1e-3; %transform Displacement from MM to METERS
overagain=0;
while overagain==0
Xs=smooth(trac(:,1)); %smooth displacement (moving average - 5 pts)
F=smooth(trac(:,2));
loopchoice=0;
while loopchoice==0
try
[HexpBo,HexpTo,TOPmin,TOPmax,BOTmin,BOTmax]=getloopdisp([Xs,F]);
catch
fprintf('\nThis is the last loop! What do you want to do?\n')
overagain=input('OVER AGAIN(0), SELECT THIS ONE(1): ');
fprintf('\n');
break
end
[HexpT,HexpB]=centercurve1loop(HexpTo,HexpBo);
figure(1); clf;
plot(trac(:,1),trac(:,2),'k');
hold on;
plot(HexpT(:,1),HexpT(:,2),'b','LineWidth',2);
plot(HexpB(:,1),HexpB(:,2),'r','LineWidth',2);
title(sprintf('SAMPLE3_02mm_trac%g_%gHz.mat;',...
i,freqs(j)),'Interpreter','none');
loopchoice=input('Like it? NO(0) / YES(1) / EDIT(2): ');
if loopchoice==1
overagain=1;
elseif loopchoice==2
break
end
Hmin=min([TOPmin,TOPmax,BOTmin,BOTmax]);
Xs=Xs(1:Hmin); F=F(1:Hmin);
end %while loopchoice
if loopchoice==2 %Input points with mouse
epoa=0; %EditPoints OverAgain
while epoa==0
figure(1); clf; plot(trac(:,1),trac(:,2),'k');
title(sprintf('SAMPLE2_02mm_%gHz.mat;',freqs(j)),...
118
'Interpreter','none');
fprintf('\nPick points with LEFT mouse button\n');
fprintf('When ready, click RIGHT mouse button\n');
fprintf('(Create OVER a whole loop, go CLOCKWISE!)\n');
button=1; k=1; clear xi yi;
while button==1
[x,y,button]=ginput(1);
if button==1
xi(k)=x; yi(k)=y;
if k==1
hold on;
plot(xi,yi,'xm','MarkerSize',8);
else
hold on;
plot(xi(k),yi(k),'xm','MarkerSize',8);
plot(xi(k-1:k),yi(k-1:k),'r','LineWidth',2);
end
k=k+1;
end
end
[HexpTo,HexpBo,TOPmin,TOPmax,BOTmin,BOTmax]=...
getloopdisp([xi',yi']);
[HexpT,HexpB]=centercurve1loop(HexpTo,HexpBo);
figure(1); clf;
plot(trac(:,1),trac(:,2),'k');
hold on;
plot(HexpT(:,1),HexpT(:,2),'b','LineWidth',2);
plot(HexpB(:,1),HexpB(:,2),'r','LineWidth',2);
title(sprintf('SAMPLE2_02mm_%gHz.mat;',freqs(j)),...
'Interpreter','none');
epoa=input('Like it? NO(0) / YES(1): ');
if epoa==1
overagain=1;
end
end
end
end %while overagain
clear HexpTo HexpBo;
eval(sprintf('save SAMPLE2_02mm_%gHz HexpT HexpB;',freqs(j)));
end
function [HexpT,HexpB]=centercurve1loop(HexpTo,HexpBo)
%This function centers a hysteresis loop
Tlength=length(HexpTo);
Blength=length(HexpBo);
x=[HexpTo(:,1);HexpBo(:,1)];
F=[HexpTo(:,2);HexpBo(:,2)];
xmax=max(x);
xmin=min(x);
Fmax=max(F);
119
Fmin=min(F);
xpk=(xmax-xmin)/2;
Fpk=(Fmax-Fmin)/2;
dx=xmax-xpk;
dF=Fmax-Fpk;
Hc=[HexpTo(:,1)-dx , HexpTo(:,2)-dF;...
HexpBo(:,1)-dx , HexpBo(:,2)-dF];
HexpT=Hc(1:Tlength,:);
HexpB=Hc(Tlength+1:Tlength+Blength,:);
120
APPENDIX B
PROGRAM CREATED IN NATIONAL INSTRUMENTS LABVIEW 8.0 FOR
THE ACQUISITION OF EXPERIMENTAL HYSTERESIS LOOPS
121
FRONT PANEL:
122
DIAGRAM: