University of Southern California Dissertations and Theses

Studies into computational intelligence approaches for the identification of complex nonlinear systems

(USC Thesis Other)

Studies into computational intelligence approaches for the identification of complex nonlinear systems

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Studies into Computational Intelligence Approaches for the Identification
of Complex Nonlinear Systems
by
Seyed Ali Bolourchi Yazdi
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements of the Degree
DOCTOR OF PHILOSOPHY
(Engineering)
May 2014
Copyright 2014 Seyed Ali Bolourchi Yazdi

Dedicated to my father, my mother and my brother.
ii
Acknowledgements
First and foremost, I would like to thank my advisor, Professor Sami Masri, for his inspiration,
motivation, guidance and support throughout my post-graduate study. His breadth of knowledge,
excitement for new discoveries, appreciation for innovation and enthusiasm for crossing disci-
plinary borders to explore profound concepts encouraged me to conduct an unconventional and
exceptional multidisciplinary research crossed between Structural Engineering, Electrical Engi-
neering, Computer Science, Mechanical Engineering, and Biology. I am so honored and delighted
to have him as my mentor, and have learnt everlasting lessons for my life.
I would like to thank the members of my defense and qualifying committee, Professors Roger
Ghanem, Carter Wellford, Jerry Mendel, and James Anderson for their valuable recommenda-
tions and encouragement.
I would like to thank my friends in our research group that have helped me throughout my jour-
ney, Drs. Reza Jafarkhani, Mohammadreza Jahanshahi, Hae-Bum (Andrew) Yun, and (soon-to-
be Drs.) Miguel Hernandez-Garcia, Armen Derkevorkian, Vahid Keshavarzzadeh and Preetham
Aghalaya Manjunatha. Thank you for making our experience so memorable. I would also like to
extend my gratitude to Professor Najmedin Meshkati and Dr. Farzad Naeim for their support.
I would like to thank Viterbi School of Engineering, the Department of Civil and Environmen-
tal Engineering and the Department of Physics and Astronomy for awarding me fellowships,
teachings assistantships and other awards. Special thanks to Professor Erik Johnson, Professor
Jean-Pierre Bardet, Professor Lucio Soibelman, Dr. Gokhan Esirgen, Margeri Berti, Gennifer
Gerson and Jennifer Craig.
iii
My education at USC could not be more exciting without countless friends and colleagues I have
made throughout this journey. I also need to thank all my students in classrooms and labs that
made teaching one of the most phenomenal and enjoyable experiences of my entire life.
Finally, I extend my biggest thank you to my dad Hossein, my mom Aghdas and my brother Mo-
hammad for their unconditional love, support and motivational words during tough times. By all
means, I owe my success to them.
iv
Table of Contents
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 The Philosophical Motivation behind Pursuing PhD . . . . . . . . . . . . 1
1.1.2 The Scientific Motivation behind Conducting this Research . . . . . . . . 4
1.2 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Computational Intelligence, Stochastic Optimization and Evolutionary Algorithms 9
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Genetic Programming: An Overview . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Population . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Evolutionary Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.3 Fitness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Bloat Control using Linear Algebra . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Parameter Optimization using Genetic Algorithms . . . . . . . . . . . . . . . . . 19
2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3 Identification of SDOF Systems with Polynomial Nonlinearities 21
v
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.2 Review of Nonparametric Approaches for the Identification of Nonlinear
Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.3 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2 Identification Scheme Composed of GP, LA and GA . . . . . . . . . . . . . . . 25
3.2.1 Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.2 Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3 Identification of Linear Oscillator . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4 Duffing Oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5 Identification of van der Pol Oscillator . . . . . . . . . . . . . . . . . . . . . . . 44
3.6 Identification of Noisy Duffing-van der Pol Oscillator . . . . . . . . . . . . . . . 50
3.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4 Identification of SDOF Systems with Non-Polynomial Nonlinearities 60
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.1.1 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2 Identification Scheme Composed of GP, LA and GA . . . . . . . . . . . . . . . 62
4.2.1 Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.3 Identification of Coulomb Friction Oscillator . . . . . . . . . . . . . . . . . . . 65
4.4 Identification of Tri-Linear Oscillator . . . . . . . . . . . . . . . . . . . . . . . 70
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.5.1 Computational Advantages of the Introduced Identification Technique . . 76
vi
4.6 Advantages and Limitations of GP Compared to other Identification Methods . . 79
4.6.1 Advantages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.6.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.7 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5 Identification of SDOF Systems with Hysteretic Behavior 82
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.1.2 Review of Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.1.3 Developments in Computational Intelligence Approaches . . . . . . . . . 84
5.1.4 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.2 Identification Scheme Using Genetic programming . . . . . . . . . . . . . . . . 86
5.2.1 General Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.2.2 Investigated Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2.3 External Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2.4 Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2.5 Identification Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.3 Identification of SDOF Hysteretic Systems Governed by Bouc-Wen Model . . . . 93
5.4 Identification of SDOF Hysteretic Systems Exhibiting Bilinear Behavior . . . . . 94
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
6 Modeling of Human Spine 108
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.1.1 Review of Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.1.2 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
vii
6.2 The Method of Experimentation . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.2.1 Experimental Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.2.2 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.2.3 Data Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.3 Non-Parametric Modeling of Human Spine Using Genetic programming . . . . . 115
6.3.1 Modeling Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.4 Identification Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.6 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
7 Identification of MDOF Systems with Hysteretic Behavior 125
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
7.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
7.1.2 Review of Related Studies . . . . . . . . . . . . . . . . . . . . . . . . . 126
7.1.3 Advances in Computational Intelligence . . . . . . . . . . . . . . . . . . 127
7.1.4 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
7.2 Introduction to Genetic Programming . . . . . . . . . . . . . . . . . . . . . . . 128
7.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7.2.2 Population . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
7.2.3 Evolutionary Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.2.4 Fitness Criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.3 Identification of MDOF Systems: General Procedure . . . . . . . . . . . . . . . 131
7.3.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.3.2 Identification Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . 136
viii
7.4 Identification of MDOF Hysteretic Systems Characterized by Bouc-Wen Model . 141
7.4.1 Training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
7.5 Validation and Discussion: Statistical Analysis . . . . . . . . . . . . . . . . . . 146
7.6 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
8 Summary and Conclusion 163
Bibliography 167
ix
List of Figures
2.1 Crossover operation on chromosomes. . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Crossover of trees in GP: (a) Parent 1 before crossover; (b) Parent 2 before
crossover; (c) Offspring 1 after crossover; (d) Offspring 2 after crossover. . . . 14
2.3 Mutation operation on a chromosome. . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Mutation operation on a tree in GP: (a) A tree before mutation with a randomly
selected branch; (b) A newly generated tree to replace a branch of another tree;
(c) The result of the mutation operation. . . . . . . . . . . . . . . . . . . . . . 16
3.1 Flowchart of the proposed technique for the identification of systems with polynomial-
type nonlinearities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 The trees of the linear oscillator’s governing equation of motion. . . . . . . . . 33
3.3 The external force and statistics of the linear oscillator response used for training. 34
3.4 Comparison of the actual and identified response of the linear oscillator to the
exciting force used for training. . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.5 Comparison of the actual and identified response of the linear oscillator to a new
exciting force for validation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.6 The trees of the Duffing oscillator’s governing equation of motion. . . . . . . . 40
3.7 The external force and statistics of the Duffing oscillator response used for train-
ing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.8 Comparison of reference and identified response of the Duffing oscillator to the
exciting force used for training. . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.9 Comparison of reference and identified response of the Duffing oscillator to a
new exciting force for validation. . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.10 The tree structures of the van der Pol oscillator’s governing equation of motion. 46
3.11 The external force and statistics of the van der Pol oscillator response used for
training. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.12 Comparison of the reference and identified response of the van der Pol oscillator
to the exciting force used for training. . . . . . . . . . . . . . . . . . . . . . . 48
3.13 Comparison of the reference and identified response of the van der Pol oscillator
to a new exciting force for validation. . . . . . . . . . . . . . . . . . . . . . . 49
3.14 The tree of the noisy Duffing-van der Pol oscillator’s governing equation of
motion obtained with 5% noise. . . . . . . . . . . . . . . . . . . . . . . . . . 54
x
3.15 External force and statistics of the noisy Duffing-van der Pol oscillator response
used for training: (a) excitation; (b) standard deviation of displacement. . . . . 55
3.16 Comparison of the actual and identified response of the noisy Duffing-van der
Pol oscillator to the exciting force used for training: (a) time history of dis-
placement; (b) time history of velocity; (c) time history of acceleration; (d) time
history of restoring force; (e) velocity vs. displacement; (f) restoring force vs.
displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.17 Comparison of the actual and identified response of the noisy Duffing-van der
Pol oscillator to a new exciting force for validation: (a) validation excitation
f(t); (b) external disturbance p(t); (c) time history of acceleration; (d) time
history of restoring force; (e) velocity vs. displacement; (f) restoring force vs.
displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1 Flow chart of the proposed identification scheme. . . . . . . . . . . . . . . . . 64
4.2 SDOF oscillator with Coulomb friction: (a) behavior of the friction force of
the Coulomb friction oscillator; (b) the tree of governing equation of motion
obtained with 5% noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.3 Exact vs. rational polynomials, both obtained from GP. . . . . . . . . . . . . . 69
4.4 Comparison of the actual and identified response of the oscillator with Coulomb
friction to a new exciting force for validation: (a) validation excitation; (b)
time history of displacement; (c) time history of velocity; (d) restoring force vs.
velocity; (e) Coulomb friction vs. velocity; (f) restoring force vs. displacement.
Note that three superposed curves are shown in each figure. . . . . . . . . . . . 69
4.5 Behavior of the stiffness force in the tri-linear oscillator. . . . . . . . . . . . . . 73
4.6 Tree representation of GP-obtained actual differential operator of the tri-linear
oscillator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.7 Comparison of the actual and identified response of the tri-linear oscillator to a
new exciting force for validation. (a) time history of validation excitation; (b)
time history of displacement; (c) time history of velocity; (d) time history of
acceleration; (e) time history of restoring force; (f) velocity vs. displacement;
(g) restoring force vs. displacement;(h) restoring force vs. velocity; (i) stiffness
force vs. displacement. Note that four superposed curves are shown in each figure. 75
4.8 The comparison of different combinations of the utilized techniques to verify the
superiority of the proposed approach in discovering the restoring force equation
of the oscillator with Coulomb friction . . . . . . . . . . . . . . . . . . . . . . 78
5.1 Flowchart of the proposed hybrid evolutionary technique for the identification
of hysteretic systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2 Comparison of the reference and estimated response of Bouc-Wen oscillator to
the exciting force used for training. Note 4 superposed curves are shown in each
panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
xi
5.3 Comparison of the reference and estimated response of Bouc-Wen oscillator to
the exciting force used for validation. Note that 4 superposed curves are shown
in each panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.4 Comparison of the phase plots of the reference and estimated restoring force
of Bouc-Wen oscillator to the exciting force used for validation. Note that 4
superposed curves are shown in each panel. . . . . . . . . . . . . . . . . . . . 102
5.5 Behavior of the stiffness force in the bilinear hysteretic oscillator. . . . . . . . . 103
5.6 Comparison of the reference and estimated response of the bilinear hysteretic
oscillator to the exciting force used for training. Note that 4 superposed curves
are shown in each panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.7 Comparison of the reference and estimated response of the bilinear hysteretic
oscillator subjected to the exciting forces used for validation. Note that 4 super-
posed curves are shown in each panel. . . . . . . . . . . . . . . . . . . . . . . 104
5.8 Comparison of the phase plots of the reference and estimated restoring force of
the bilinear hysteretic oscillator subjected to the exciting force used for valida-
tion with
v1
= 0:33. Note that 4 curves are superposed for different noises. . . 105
5.9 Comparison of the phase plots of the reference and estimated restoring force of
the bilinear hysteretic oscillator subjected to the exciting force used for valida-
tion with
v2
= 1:00. Note that 4 curves are superposed for different noises. . . 106
5.10 Comparison of the phase plots of the reference and estimated restoring force of
the bilinear hysteretic oscillator subjected to the exciting force used for valida-
tion with
v3
= 3:00. Note that 4 curves are superposed for different noises. . . 107
6.1 The test setup [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.2 The motion segment is shown schematically. The top picture is intact, and the
bottom picture is with X-stop implant [2]. . . . . . . . . . . . . . . . . . . . . 113
6.3 The time history of the data used for the modeling of spine. . . . . . . . . . . . 114
6.4 The time history of the normalized data used for the modeling of spine. . . . . 115
6.5 The phase plots of the data used for the modeling of spine. . . . . . . . . . . . 116
6.6 The phase plots of the normalized data used for the modeling of spine. . . . . . 117
6.7 The phase plots of torsional torque vs. torsional rotation. . . . . . . . . . . . . 118
6.8 (a) The bar graph of the coefficients of the differential equations for spinal mo-
tion segment intact (blue), injured (green) and X-stop (red)). . . . . . . . . . . 121
6.9 The solution of the short and long ODE 1 against the reference response for the
intact spine segment motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.10 The solution of the short and long ODE 1 against the reference response for the
injured spine segment motion. . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.11 The solution of the short and long ODE 1 against the reference response for the
spine segment motion with X-stop. . . . . . . . . . . . . . . . . . . . . . . . . 124
7.1 The benchmark 3 degree-of-freedom system with interconnecting links and sup-
ports. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
xii
7.2 The flowchart of the methodology presented for the identification of hysteretic
MDOF systems. The unidirectional arrows show the main stream of the training
process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
7.3 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 1. The system is stimulated by 3 distinct external training ex-
citations at all degrees-of-freedom. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 143
7.4 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 2. The system is stimulated by 3 distinct external training ex-
citations at all degrees-of-freedom. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 143
7.5 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 3. The system is stimulated by 3 distinct external training ex-
citations at all degrees-of-freedom. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 144
7.6 The error of the estimated response of the 3-DOF hysteretic system at all 3
degrees-of-freedom denoted as D1, D2, D3. The system is stimulated by 3
distinct external training excitations at all degrees-of-freedom. . . . . . . . . . 144
7.7 Phase plots of the force at interconnecting elements ofg
1
tog
6
vs. the relative
displacement of the corresponding ends when the 3-DoF hysteretic system is
subjected to distinct training excitations with
v
= 1:0 at all degrees of freedom. 145
7.8 Distribution of the percentage of the model response error at every degree of
freedom when subjected to 20 validation excitations with
v
= 0:5. Having

v
= 0:5, these excitations are less intense than the identification excitation hav-
ing
t
= 1:0. The vertical axes show the number of samples and the horizontal
axes show the corresponding error value of the same number of excitations in
the bin. To obtain an estimate about the average error, the normal distribution
of the data are also provided in each figure. . . . . . . . . . . . . . . . . . . . 148
7.9 Distribution of the percentage of the model response error at every degree of
freedom when subjected to 20 validation excitations with
v
= 1:0. Having

v
= 1:0, these excitations are as intense as the identification excitation having

t
= 1:0. The vertical axes show the number of samples and the horizontal
axes show the corresponding error value of the same number of excitations in
the bin. To obtain an estimate about the average error, the normal distribution
of the data are also provided in each figure. . . . . . . . . . . . . . . . . . . . 149
7.10 Distribution of the percentage of the model response error at every degree of
freedom when subjected to 20 validation excitations with
v
= 1:5. Having

v
= 1:5, these excitations are 1.5 times more intense than the identification
excitation having
t
= 1:0. The vertical axes show the number of samples and
the horizontal axes show the corresponding error value of the same number of
excitations in the bin. To obtain an estimate about the average error, the normal
distribution of the data are also provided in each figure. . . . . . . . . . . . . . 150
xiii
7.11 The average of the validation error associated with the model response subject to
the validation excitations with respect to the reference response. The stimulation
level of the validation excitations are substantially different from the training
excitation: Excitation 1 (
v
= 0:5) is less intense than the training excitation
(
t
= 1:0), Excitation 2 (
v
= 1:0) is as intense as the training excitation (
t
=
1:0) and Excitation 3 (
v
= 1:5) is more intense than the training excitation
(
t
= 1:0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
7.12 Phase plots of the force at interconnecting elements ofg
1
tog
6
vs. the relative
displacement of the corresponding ends when the 3-DOF hysteretic system is
subjected to distinct validation excitations with
v
= 0:5 at all degrees of free-
dom. These plots are obtained from reference responses for providing insight to
the hysteretic behavior of interconnecting elementsg
1
tog
6
. . . . . . . . . . . 153
7.13 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 1. The system is stimulated by 3 distinct external validation
excitations with
v
= 0:5 at all degrees-of-freedom which is less intense than
the training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 154
7.14 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 2. The system is stimulated by 3 distinct external validation
excitations with
v
= 0:5 at all degrees-of-freedom which is less intense than
the training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 154
7.15 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 3. The system is stimulated by 3 distinct external validation
excitations with
v
= 0:5 at all degrees-of-freedom which is less intense than
the training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 155
7.16 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 1. The system is stimulated by 3 distinct external validation
excitations with
v
= 1:0 at all degrees-of-freedom which is as intense as the
training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 156
7.17 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 2. The system is stimulated by 3 distinct external validation
excitations with
v
= 1:0 at all degrees-of-freedom which is as intense as the
training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 157
7.18 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 3. The system is stimulated by 3 distinct external validation
excitations with
v
= 1:0 at all degrees-of-freedom which is as intense as the
training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 157
xiv
7.19 Phase plots of the force at interconnecting elements ofg
1
tog
6
vs. the relative
displacement of the corresponding ends when the 3-DoF hysteretic system is
subjected to distinct validation excitations with
v
= 1:0 at all degrees of free-
dom. These plots are obtained from reference responses for providing insight to
the hysteretic behavior of interconnecting elementsg
1
tog
6
. . . . . . . . . . . 158
7.20 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 1. The system is stimulated by 3 distinct external validation
excitations with
v
= 1:0 at all degrees-of-freedom which is as intense as the
training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 159
7.21 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 2. The system is stimulated by 3 distinct external validation
excitations with
v
= 1:0 at all degrees-of-freedom which is as intense as the
training excitation with
t
= 1:0. The graphs of the reference and estimated
response are plotted using solid and dotted lines, respectively. . . . . . . . . . . 159
7.22 Comparison of the reference and estimated response of the 3-DOF hysteretic
system at DOF 3. The system is stimulated by 3 distinct external validation
excitations with
v
= 1:5 at all degrees-of-freedom which is 1.5 times more
intense than the training excitation with
t
= 1:0. The graphs of the reference
and estimated response are plotted using solid and dotted lines, respectively. . . 160
7.23 Phase plots of the force at interconnecting elements ofg
1
tog
6
vs. the relative
displacement of the corresponding ends when the 3-DoF hysteretic system is
subjected to distinct validation excitations with
v
= 1:5 at all degrees of free-
dom. These plots are obtained from reference responses for providing insight to
the hysteretic behavior of interconnecting elementsg
1
tog
6
. . . . . . . . . . . 161
xv
List of Tables
2.1 The properties of different elements in trees. . . . . . . . . . . . . . . . . . . . . 19
3.1 Properties of the exciting force for both training and validation, used in the linear,
Duffing, van der Pol and noisy Duffing-van der Pol oscillators. . . . . . . . . . . 29
3.2 Parameters of the linear oscillator: actual vs. identified. . . . . . . . . . . . . . . 31
3.3 The actual and discovered equations of motion of the linear oscillator and their
associated fitness error for different noise levels. . . . . . . . . . . . . . . . . . . 32
3.4 The normalized error associated with states in the training and validation phase
due to solving the discovered ODEs corresponding to the linear oscillator. . . . . 33
3.5 Parameters of the Duffing oscillator: actual vs. identified by GP. . . . . . . . . . 38
3.6 The actual and discovered equations of motion of the Duffing oscillator and their
associated fitness for different noise levels. . . . . . . . . . . . . . . . . . . . . . 38
3.7 The normalized error associated with states in training and validation phase due
to solving the GP-found ODEs. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.8 Parameters of the van der Pol oscillator: actual vs. identified. . . . . . . . . . . . 45
3.9 The actual and discovered equations of motion of the van der Pol oscillator and
their associated fitness for different noise levels. . . . . . . . . . . . . . . . . . . 45
3.10 The normalized error associated with states in the training and validation phase
due to solving the GP-found ODEs corresponding to the van der Pol oscillator. . . 45
3.11 Parameters of the exciting force for both training and validation and the evolu-
tionary parameters of the noisy Duffing-van der Pol. . . . . . . . . . . . . . . . . 51
3.12 Parameters of the noisy Duffing-van der Pol oscillator, actual vs. identified. . . . 53
3.13 The actual and discovered differential equations of the noisy Duffing-van der Pol
oscillator and their associated RMS error for different noise levels. . . . . . . . . 53
3.14 The normalized error associated with estimated states in training and validation
phase due to solving the GP-found ODEs corresponding to the noisy Duffing-van
der Pol oscillator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1 Properties of the applied force for both training and validation used in a SDOF
system with Coulomb friction and the corresponding evolutionary parameters. . . 67
4.2 Properties of the oscillator with Coulomb friction, actual vs. identified. Mass
assumed known. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
xvi
4.3 The actual and discovered equations of motion of the oscillator with Coulomb
friction and their associated RMS error for different noise levels. . . . . . . . . . 68
4.4 The RMS error associated with the estimated states in the training and validation
phase due to solving the GP-found ODEs corresponding to the oscillator with
Coulomb friction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.5 Properties of the tri-linear oscillator, the applied excitations for both training and
validation, and the evolutionary parameters. . . . . . . . . . . . . . . . . . . . . 72
4.6 The actual and discovered equations of the restoring force of the tri-linear oscil-
lator and their associated RMs error for different noise levels. . . . . . . . . . . . 72
4.7 The RMS error associated with the estimated states in training and validation
phase due to solving the GP-found ODEs corresponding to the tri-linear oscillator. 73
5.1 Properties of GP and GA for discovering the optimum structure and optimum
parameters, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.2 Properties of the Bouc-Wen model, the applied excitations for both training and
validation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.3 The actual and discovered equations of motion of the Bouc-Wen oscillator and
their associated RMS error when polluted and unpolluted data are used. . . . . . 95
5.4 The normalized error associated with the estimated states in the training and vali-
dation phase due to solving the GP-found ODEs corresponding to the Bouc-Wen
oscillator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.5 The normalized error associated with the estimated states in the training and val-
idation phase due to solving the GP-found ODEs corresponding to the bilinear
hysteretic oscillator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.6 Properties of the bilinear hysteretic oscillator and the applied excitations for both
training and validation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.7 The actual and discovered equations of motion of the bilinear oscillator and their
associated RMS error, when polluted and unpolluted data are used. . . . . . . . . 100
6.1 Properties of evolutionary based identification approach for discovering the dif-
ferential equation that governs the motion of the spine. . . . . . . . . . . . . . . 119
6.2 The optimized coefficients of the ODEs and the RMS error associated with the
rotation in the denormalized domain. . . . . . . . . . . . . . . . . . . . . . . . . 120
7.1 Properties of the Bouc-Wen model and the applied excitationsf
i
(t) for training
and validation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
7.2 Properties of GA for parameter optimization. . . . . . . . . . . . . . . . . . . . 140
xvii
Abstract
This study builds on major advances in the field of Computational Intelligence to develop a state-
of-the-art data-driven methodology that provides parsimonious optimized computational models
in the form of systems of differential equations that characterize the behavior of complex nonlin-
ear phenomena observed in mechanical and biological systems. The proposed hybrid identifica-
tion scheme integrates various stochastic optimization methods and computer algebra techniques,
such as Genetic Programming and Genetic Algorithms, to evolve structures of differential equa-
tions, to optimize their parameters, and to reduce their complexity for controlling bloat. The in-
vestigated scenarios include systems that exhibit polynomial-type nonlinearities in their response,
systems that show discontinuity in their nonlinear behavior, systems with memory-dependent and
dissipative characteristics, as well as the human spine. The investigations are conducted by pro-
cessing input and output data obtained from synthetic simulations as well as experiments. It is
shown that the proposed technique yields reduced-order, reduced-complexity, optimized differen-
tial equations, that accurately characterize the behavior of the investigated systems, and provide
accurate estimates. The generalization extent of the discovered models is scrutinized by assess-
ing their performance in new dynamical environments through applying validation excitations
that are substantially different from the excitations employed for training. Findings reveal that
the resulting models provide reasonably accurate estimates, even when models are subjected to
new stimulations with various intensities. Thus, the proposed approach of this study presents
a robust data-driven methodology based on evolutionary computation techniques that provides
elegant computational models to represent variety of complex nonlinear systems.
xviii
Chapter 1
Introduction
1.1 Motivation
T
HE motivation behind this research is described from two different perspectives: the
philosophical inspiration and the technical aspect.
1.1.1 The Philosophical Motivation behind Pursuing PhD
How does this world function? Human beings have repeatedly asked this question for thousands
of years. This question is rooted in the curiosity trait of humans, and their essential need to em-
ploy and shape the nature to fulfil their goals.
The answer to this question has always been limited to the knowledge of the people of the same
era, as well as the sophistication of the available tools to perceive their environment. In the ab-
sence of modern science and currently available advanced instruments, humans’ understanding
of natural phenomena was drastically influenced by their inaccurate observations. It was only a
few centuries ago when the emergence of modern science opened new horizons and resulted in
unprecedented discoveries and developments that, according to the history of science, had never
1
existed with this degree in ancient times.
The scientific revolution helped humans ignite a growing and evolving machine that is self-
constructing, self-correcting and is fueled by the thirst for knowledge of human beings and our
increasing need of resources. This machine is the most important and the most valuable legacy
of the human population on Earth throughout the history. I name this machine HESPHE because
it is made of the knowledge from all major existing fields: Humanities, Engineering, Science,
Professional fields, Health and Education. While we, as humans, improve the Hesphe machine,
this machine is intertwined with our society and culture, and the mutual impact of our beliefs and
science is not hidden to the eye.
The Hesphe machine has significantly changed our view of the world. New discoveries are
made everyday and advanced technologies are introduced to the society on a daily basis. In early
stages of these advancements, new findings were helping us better understand the world, the envi-
ronment, living beings, and every thing around us. Bringing more convenience to our lives, new
discoveries facilitated the formation of societies and social interactions. They were very useful
tools in our hands. However, during the past few decades, the blossom of science and technol-
ogy started to change human’s mind and lifestyle. New explorations questioned many religious
believes and challenged personal integrity. Soon we realized that our societies are significantly
affected by new innovations that are developed just a few years ago in the laboratory of high-tech
companies. They are widely used worldwide, and bringing extremely fast and inevitable changes
to society, our lives and our thoughts. The danger is that we are not using technology as a tool any
more, but it seems sometimes we are some tools in the hands of technology. We force ourselves
to start using state-of-the-art technology to stay up to date, but soon we realize that we cannot
2
live without them. For example the rise of social media is undoubtedly one of the most influential
innovations of our era. They affect our relations with our government, society, family, friends,
and our interactions with ourselves. Publicly available information from our family, friends and
colleagues on social media make us think twice about ourselves, how we behave and how we
should behave to achieve more and to feel accepted in the society.
The destination of human race in such a fast transforming society is by no means certain nor
predictable. Perhaps we need more time to fully digest different aspects of modern discoveries
and innovations, and learn how to use them in our favor rather than using them recklessly as
slaves. Still new scientific discoveries will challenge our beliefs more, and we should revisit our
thoughts and figure out what the roots of our beliefs are, the reason we find them valuable, and
how we should preserve them against unwillingly bombarding transformations around us. Then,
we are able to pave the route that is made by our will, founded on our deep knowledge and wis-
dom, and equipped with assisting technologies to take us to the goal we have defined in our life
for ourselves. This exclusive path neither is similar to the path of others nor should be influenced
by the advertisement of big corporations whose aim is to sell their products having short lifes-
pan. Our dream goal can be the happiness of ourselves, our family and friends, our people, the
environment, and the world as a whole. Therefore, the opportunity for every one to be able to
pursue happiness forms the foundation of all values and beliefs. It is the birthplace of life, justice,
equality, truth, freedom, self-responsibility, diversity, integrity, donation, caring of others, sacri-
fice, self-respect and others. We spend our entire life in this path to accomplish these goals. But
life is a unique opportunity that is given to us, and time is the most valuable asset we have. Thus,
we should be very careful about how we use our life, what we use it for and what we will achieve
down the road after we are old. We would like to be proud of ourselves when we look back at our
3
accomplishments, relations and whatever we leave behind. Various paths can be defined, however
the route should not distract us from our final goals. It’s our dream goals that keep us hopeful and
alive. Because in the storms of our life, it is hope that keeps our inner candle lit.
The research presented in this dissertation tries to add a piece to the Hasphe machine to
advance science. That’s one of the major steps in returning to society after studying and learning
for many years. Pursuing PhD as a whole, and conducting research as an element of the PhD,
equip us with fantastic skills, make us ready for our future stages of life, and eventually help us
follow our dreams, to make world a better place.
1.1.2 The Scientific Motivation behind Conducting this Research
Modeling plays an important role in our interaction with nature and human-made devices. Mod-
eling provides information about underlying systems, and enables us to understand them, predict
their behavior, identify their problems and suggest approaches to overcome their existing issues.
The complexity of many systems and their dependency on many parameters cause the model-
ing to be challenging. Various techniques and approaches are introduced in the literature for the
identification and modeling complex systems. These approaches can be classified into two ma-
jor categories: Parametric Modeling and Non-Parametric Modeling. In parametric modeling, a
model is assigned to the system to evaluate the performance of the system or monitor its changes.
This approach is more applicable to well-studied systems for which a clear model is available.
However, finding a properly-defined model for many complex systems is not feasible. Lack of
information about the system, dependency on many parameters and the complex nature of many
phenomena are among the limitations that direct researchers toward Non-Parametric Modeling.
Non-Parametric approaches mainly use the available data to infer about the unknown systems.
4
They may not result in physically interpretable models in many scenarios, but still yield surrogate
models that can represent the underlying systems fairly accurately. Yet, the scarcity of Non-
Parametric approaches that provide robust high-fidelity parsimonious models that can accurately
characterize complex nonlinear systems exist in the literature.
The advent of computers and the emergence of fast computing capabilities and artificial in-
telligence in the past few decades, sparked a new era in human history: human-made machines
could “think” instead of their creator with the degree of elegance that was not achievable before.
Scientists have used computers to better understand very big phenomena in our universe such as
the clash of galaxies, the evolution of human’s race, social relations in societies, climate change,
nuclear explosions as well as very small phenomena such as the behavior of virus in a victims’s
body, the interaction of subatomic particles, the effect of certain drugs on neurons, the influence
of a specific gene on cancer development, etc. Hence, growing computational abilities of com-
puters has widened the perspective of researchers, and equipped them with a powerful tool based
on which new approaches for solving complicated problems can be founded.
Evolutionary approaches are a class of computationally intensive algorithms extensively used
to solve complex problems. These approaches are rooted in Darwin’s theory of natural selection
and population genetics, that explain how living beings have evolved. Evolutionary algorithms
are based on the following: A population of suitable individuals is generated at random. Each
individual is formed by certain rules to fit the environment, and is potentially capable of perform-
ing the required task. However, how well an individual perceives the environment, determines its
desirability, its fitness level and its competence to transfer its genes to the next generation. Two
major evolutionary operations, crossover and mutation, combine the genes of the population of
5
the former generation to construct a new generation. Theory of natural selection suggests that as
the population evolve, individuals adapt better to the environment, and consequently the fitness
of the population improves. The evolution stops when an acceptable individual is found to fulfill
the desired goal.
Genetic Programming (GP) is an advanced variant of Evolutionary Algorithms for evolving
computer programs. GP utilizes the same evolutionary principals to advance the evolution of
the population. Much research has been conducted on the efficiency and applicability of GP for
modeling various phenomena. Hence, it has the potential to be considered for the non-parametric
identification of different classes of problems. Thus, the effectiveness of GP together with other
improving techniques are investigated in this research for identification purposes. The eventual
goal is to provide an elegant methodology based on Computational Intelligence and Evolutionary
approaches to fill the existing gap in the literature for the identification of various classes of
complex nonlinear dynamical systems.
1.2 Scope
The presented research builds on the advancements in the field of Computational Intelligence and
Evolutionary Algorithms to establish the foundations of a reliable methodology for the identifica-
tion of various complex nonlinear dynamical phenomena, ranging from mechanical to biological
systems. The proposed technique employs Genetic Programming as the heart of the identification
procedure, and combines it with stochastic optimization techniques for parametric optimization,
as well as Linear Algebra rules for symbolic simplifications. The proposed identification tech-
nique yields robust representative models in the form of differential equations to fully characterize
6
the behavior of the underlying systems. The material in this dissertation is organized in the fol-
lowing order:
Chapter 2 provides an overview of Computational Intelligence, Stochastic Optimization and
Evolutionary Algorithms to base the required foundation for presenting the proposed identifica-
tion procedure for discovering the governing differential equations of nonlinear systems. The
foundations of Genetic Programming for structure optimization, Genetic Algorithms for param-
eter optimization, and Linear Algebra for symbolic simplification as well as bloat control, are
presented in this chapter.
Chapter 3 presents the proposed technique for the identification of SDOF systems with poly-
nomial nonlinearities. The investigated systems in this chapter include linear, Duffing, van der
Pol and noisy Duffing-van der Pol oscillators. The identification results are presented and the
effectiveness of the approach is discussed.
Chapter 4 presents the suggested technique for the identification of SDOF Systems with non-
polynomial nonlinearities. The studied systems include Coulomb friction and the tri-linear oscil-
lator. These dynamical systems exhibit discontinuity and sharp changes in their response, respec-
tively. The suggested identification method is expanded for this class of nonlinear systems, and
the identification results are presented and discussed.
Chapter 5 presents the suggested approach for the identification of SDOF systems with hys-
teretic behavior that includes the Bouc-Wen and bilinear hysteretic oscillators. These memory-
dependent phenomena are modeled by a system of differential equations, which are obtained
7
using the proposed identification scheme. The identification results are presented, and important
points are discussed.
Chapter 6 presents the modeling of a human spine using experimental data. The method of
the experimentation, the properties of the data set, the modeling procedure, and the results are
presented. The generalizability of the discovered models and the effectiveness of the proposed
approach are discussed.
Chapter 7 provides the proposed technique for the identification of MDOF systems exhibiting
hysteretic behavior. First, distinct differential equations are obtained for each degree-of-freedom
using Genetic Programming. Then, the un-optimized system of differential operators is formed,
and Genetic Algorithms are applied to optimize the parameters of the model using the error of
the response estimate as the cost function. The identification procedure is explained in detail, and
the identification results are presented. A statistical analysis is performed on the accuracy of the
response estimates when the MDOF system is stimulated by various validation excitations with
different intensities.
Chapter 8 provides the summary as well as the concluding remarks of this dissertation.
8
Chapter 2
Computational Intelligence, Stochastic Optimization and
Evolutionary Algorithms
2.1 Introduction
O
VER the past decades, the advent of computers and the blossoming of machine intelli-
gence sparked a new era in computational abilities, and led to the emergence and flourish
of new techniques that overcame the limitations that have been challenging researchers for many
years. Having benefited significantly from fast computing capabilities, Computational Intelli-
gence (CI) is an area in which design of intelligent agents is studied and researched. Intelligent
agents are systems that perceive their environment and consequently perform proper actions to
accomplish a desired task, e.g., provide the solution to a problem, or make decisions. In CI, solv-
ing problems is carried out by searching the domain of possible solutions intelligently. Various
optimization techniques are utilized to conduct the search for a solution more efficiently and effec-
tively. Inspired by Darwin’s theory of evolution [3], Evolutionary Algorithms, such as Evolution
Strategies and Genetic Algorithms, use the model of natural selection as the main core of their
9
optimization machine. The main idea behind all evolutionary algorithms is to generate a popula-
tion whose elements can potentially be a solution, and then evolve them by applying appropriate
evolutionary operations in consecutive generations with the consideration of individuals’ fitness
to obtain the desired solution.
Genetic Algorithms (GA) are known to be the most popular type of Evolutionary Algorithms.
They are widely used for solving optimization problems whose solutions are usually in the form
of numbers. GAs have many application in various fields such as chemistry [4], mathematics
[5], computer-aided design [6], control engineering [7], electronic circuit design [8], computa-
tional finance [9], music [10], parametric optimization [11], etc. The complexity of many prob-
lems prevents most of the non-heuristic optimization algorithms, such as direct search technique,
gradient-based approaches and the downhill simplex method, to perform effectively. Therefore,
due to the complexity of the investigated models in this research, GAs are employed for the pa-
rameter optimization of the problems under discussion in this research. More information on GAs
for parametric optimization and its application in various disciplines can be found in [12].
Genetic Programming (GP), first introduced by Koza in 1992 [13], is a special variant of
Genetic Algorithms, and a unique optimization scheme through which a population of computer
programs, represented in expressions trees, evolve to satisfy a fitness landscape determined based
upon the considered task. This technique has been utilized to tackle various problems in different
disciplines such as decision making [13], pattern recognition [14], robotic networks [15], music
[16], bioinformatics [17], data mining [18], traffic analysis [19], finance [20], antenna design [21],
helicopter hovering control [22], physics [23], cancer prediction [24], etc.
10
In the proposed approach, Genetic Programming is used as the heart of the optimization tech-
nique for finding suitable structures, and Genetic Algorithms for optimizing the parameters of
structures. This approach is also integrated with Linear Algebra techniques for carrying out sym-
bolic simplifications. An overview of GP and the way it is adapted in this research for identifica-
tion purposes is described next.
2.2 Genetic Programming: An Overview
Evolutionary Algorithms (EAs) are global optimization techniques, rooted in the principles of
natural selection [3] and population genetics [25], and are incorporated in all four revolutionary
fields of Evolutionary Computation, i.e., Evolutionary Programming [26], Evolution Strategies
[27, 28], Genetic Algorithms [29] and Genetic Programming [13]. The primary idea behind EAs
is that first, a population, whose individuals can potentially perform the task in hand, is generated.
A fitness measure, based on a pre-established criterion and the desired goal is established and
assessed for the evaluation of every individual, indicating how successful each individual is in
achieving the desired goal, and consequently in transferring its genes to the next generation.
After individuals are ranked based on the assigned fitness, genetic operators, that are mutation,
recombination, and representation, manipulate the genetic materials, to form a new population.
Exercising evolutionary operations enhances the quality of the population over generations. The
entire process repeats until the termination criteria is reached.
GP is a landmark in evolutionary approaches that was first popularized by Koza in 1992
[13]. It utilizes the concept of EAs to evolve a population of computer programs presented in
expression trees. The most significant advantage of GP is its capability of embracing any form of
11
structure whose functionality can be evaluated quantitatively with respect to a favorable perfor-
mance. Hence, GP can be adapted to accommodate a wide range of problems in the science and
engineering fields to obtain accurate solutions with a degree of “elegance” that were not attainable
before. The current study presents the adaptation of GP for identifying different nonlinear com-
plex systems via constructing the corresponding governing ODEs and incorporating parameter
optimization methods using GAs and bloat control techniques using linear algebra.
2.2.1 Population
A suitable computer program maintains the population in data-structures. The population size is
specified based on the computational capacities. The tree structures are constructed by a com-
position of variables, numbers, and mathematical operations, and then their behavior in the en-
vironment is estimated quantitatively by means of the pre-established fitness measure, and with
the consideration of the input data. GP selects the structural components from a pool of operators
to construct trees, without any human interference. Terminal nodes are filled with variables and
numbers, while non-terminal nodes are occupied by operators including basic algebraic opera-
tions (+,,,=), functions such assin,cos,abs,sgn andstep or even user-defined functions,
e.g.,f(x) =cos(1:25x) +x
2
. The evaluation of trees is accomplished by knowing the input and
considering the inherent trait of the operators appearing in tree nodes.
Helpful but not necessary, any prior knowledge about the physics of the system being investi-
gated narrows down the number of operators in the library and boosts the convergence rate. Yet,
the library should be sufficiently comprehensive to contain a wide range of functions that may
contribute to a potential solution. Examples of typical individuals are presented in Fig. (2.2) and
Fig. (2.4). Terminal nodes in circles, correspond to numbers or input variables, and non-terminal
nodes in squares, correspond to mathematical operations.
12
2.2.2 Evolutionary Operators
In order to simulate mating, mutation and ageing of the population, evolutionary operators are
implemented as follows:
Recombination
The recombination, also called corssover, is an evolutionary operator that comprises the inter-
change of genetic material between parents (chromosomes), schematically described in Fig. (2.1).
In GP this operation is executed by exchanging a randomly selected branch of a parent tree with
another randomly selected branch from another parent tree to produce two new offsprings, illus-
trated in Fig. (2.2). This operation preserves the genetic materials of parents and transfers them
to the next generation but combines them differently.
Figure 2.1: Crossover operation on chromosomes.
Mutation
Mutations are generally meant to alter the genetic materials of chromosomes slightly, and in a
random manner, as described schematically in Fig. (2.3). Yet, in GP this operation may account
for a significant change in the tree behavior because it is carried out by randomly selecting and
13
(a) (b)
(c) (d)
Figure 2.2: Crossover of trees in GP: (a) Parent 1 before crossover; (b) Parent 2 before crossover;
(c) Offspring 1 after crossover; (d) Offspring 2 after crossover.
replacing a tree branch by a randomly generated new tree as demonstrated in Fig. (2.4). The
mutation operation is performed at a low rate, and is less likely to find an optimum solution
but intended to provide the chance of recovering good genetic materials if they are lost during
recombination and selection operations. Mutations guarantee that the probability of exploring
14
any point on the entire search domain in not zero, and also assist the optimization algorithm to
escape from local minima.
Figure 2.3: Mutation operation on a chromosome.
Representation
In order to imitate aging and take advantage of highly-fit individuals in a population for more
than a generation, the fittest individuals are transferred to the next generation while still intact.
Representation also preserves the fittest individuals for the next generations against destructive
manipulations of the mutation and recombination.
Selection
Performing genetic operations requires selecting individuals from the population. A variety of
selection schemes can be found in the literature. Truncation selection is ranking the candidate so-
lutions based on their fitness and the survival of the fittest for reinsertion into the population in the
next generations, or for consideration in other evolutionary operations [30]. Fitness-proportionate
selection, also known as roulette-wheel selection, is based on assigning a probability of survival
to every individual on the basis of the relative fitness of the individual and overall fitness of the
population. Yet, it produces bias and leads to genetic drift [31]. Stochastic universal sampling
(SUS) is an alternative approach to fitness-proportionate selection that eliminates the statistical
bias [32]. Another method, extensively implemented for selection and suitable for parallelism, is
15
(a) (b)
(c)
Figure 2.4: Mutation operation on a tree in GP: (a) A tree before mutation with a randomly
selected branch; (b) A newly generated tree to replace a branch of another tree; (c) The result of
the mutation operation.
tournament selection [33]. It is performed by randomly choosing a tour of individuals from the
population and selecting the fittest. The two latter methods provide the opportunity for weaker
individuals to survive because they may contain high-quality genetic material that may aid the
16
evolution of generations in future. Tournament selection is the approach considered in this inves-
tigation [13].
2.2.3 Fitness
The primary goal of EAs is to direct the optimization process to increase the level of adaptation
of individuals to the environment. Due to the high dependency of the evolution process on how an
individual perceives its environment and the construction of the fitness measure, establishing the
fitness criterion requires close investigation and scrutiny in order to carefully reflect the eventual
aim. Since, for identification purposes, the agreement between a candidate model response and
the actual response of the system under investigation represents an admissible behavior, herein,
the normalized error between the target signal and its estimate is considered for the fitness error
determination. The fitness error associated with the response estimate ^ s with respect to the
reference response s is calculated as the ratio of the root-mean-square (RMS) value of s and the
RMS value of the deviation e = s ^ s, as follows:
=
ks ^ sk
ksk
=
kek
ksk
(2.1)
A higher fitness value increases the chance of survival of an individual, and consequently the
presence of its genes in future generations.
2.3 Bloat Control using Linear Algebra
Bloat is the indefinite growth of trees that results in consuming available of computer memory
and increases the optimization time to reach the target solution [34–36]. Also, it is understood
17
that more parsimonious individuals exhibit better generalization abilities in future (not-seen-in-
training) dynamical environments [37, 38]. Since genetic operators reorganize the genetic mate-
rial blindly, the formation of redundant and dependent terms is inevitable, and causes bloat. There
is an extensive research literature on bloat control through manipulating the evolution process,
regardless of the problem under discussion. These techniques include improving the resource
allocation, dynamically changing the population size, varying the maximum depth of the trees
based on the performance of the population, etc. We leave the devolvement of such techniques to
researchers in the field of computer science who investigate evolutionary algorithms. However,
the nature of the problems of our interest is such that Linear Algebra rules can dramatically re-
duce the size of trees, while at the same time, this technique can be coupled with the available
methods for bloat control based on the evolution process manipulations. In order to control bloat
and increase the efficiency of the optimization technique, in this study, the mathematical form of
each tree is simplified utilizing Linear Algebra (LA) rules, through the MATLAB symbolic tool-
box, to group dependent terms and biases. It consequently results in a smaller tree when the tree
representation is reconstructed from the simplified expression. In order to perform this task, tree
structures are converted to mathematical expressions, in Matlab. These mathematical expression
are simplified symbolically, and are then re-converted to the expression trees. The conversion
from a mathematical expression to a tree structure requires extensive programming because all
of the algebraic operations, functions, variables and numbers should be identified separately, and
be treated differently, depending on their characteristic traits. Hence, the members of the mathe-
matical expressions are identified according to the following categories to be placed in the right
position in the applicable trees:
Numbers and variables are positioned in the terminal nodes while functions and algebraic op-
erators are put in the non-terminal nodes. Later in this dissertation, we will see that a maximum
18
Tree Members Characteristic
+,, and= Non-terminal nodes with two operands
sin,cos,abs,sgn andstep Non-terminal nodes with one operand
X1,X2,X3 andX4 Terminal nodes
Real numbers2 R Terminal nodes
Table 2.1: The properties of different elements in trees.
of four variables is considered in single degree-of-freedom systems that represent displacementx,
velocity _ x, acceleration x and external disturbancep(t), denoted asX1,X2,X3 andX4, respec-
tively. On the other hand, investigated multi degree-of-freedom systems in Chapter 7 involve more
than four variables. The method of investigation of these systems is discussed more extensively
in Chapter 7.
2.4 Parameter Optimization using Genetic Algorithms
Ordinary GP uses randomly generated numbers in tree structures. However, inappropriate num-
bers influence the associated differential equations drastically, and can cause them to blow up.
Moreover, unsuitable numbers lead to poor fitness, though the structure of the tree can be very
suitable. Therefore, optimized parameters improve the fitness and save the favorable structures
that might be driven to elimination otherwise. Due to the inclusion of discontinuous basis func-
tions, such as the Heaviside step function in equations, many optimization schemes such as direct
search technique, gradient based approaches and downhill simplex method are not suitable. Ge-
netic Algorithms, that are founded on the basis of evolving a population by means of evolutionary
operators to satisfy a certain fitness criterion, don’t suffer from the issues associated with many
conventional methods. Neither non-differentiability and discontinuity of the objective function,
19
nor the nature of the function results in the false detection of local extrema rather than global ex-
trema. GAs have been used extensively for parameter optimization. Therefore, in our approach,
the Genetic Algorithm (GA) is employed to optimize the parameters of the model that was shrunk
before using Linear Algebra techniques. The population size, the fitness criterion, and the termi-
nation point are specified in each problem later in the dissertation. It will be shown that exercising
GA for parameter optimization enhances the discovery of suitable differential equations that rep-
resent the underlying systems.
2.5 Summary
This chapter presented a brief introduction to Computational Intelligence approaches that incor-
porate stochastic optimization methods based on Evolutionary approaches. The eventual goal is
to provide the foundation needed for presenting the optimization technique that will be introduced
later for the identification of complex nonlinear system through discovering a suitable differential
equation. Within this concept, an overview of Genetic Programming as the heart of the pro-
posed identification technique, was provided. Important elements of GP, that are the population
type, evolutionary operators, and the fitness criterion are defined. Moreover, the bloat control
strategy, that uses Linear Algebra rules, is elaborated. The Genetic Algorithm was also intro-
duced briefly because it is extensively utilized for the parametric optimization of the discovered
structures by GP. Therefore, the required background for the introduction of the proposed iden-
tification scheme in the following chapters is provided. In the following chapters, the suggested
identification scheme and its applications to various classes of problems will be presented.
20
Chapter 3
Identification of SDOF Systems with Polynomial
Nonlinearities
3.1 Introduction
3.1.1 Background
O
VER the past several decades, significant progress has been achieved in the development
of data processing tools and algorithms employing a variety of approaches, such as data
mining, machine learning and pattern recognition, to extract useful information from data sets
corresponding to complex processes. However, in spite of these achievements, there is still a
paucity of tools that are suitable for utilization in the Applied Mechanics field to identify complex
nonlinear systems in a format that is suitable for computational simulation purposes. There is
an extensive literature in the Applied Mechanics field dealing with the broad classes of nonlinear
(i.e., not-necessarily-linear) system identification approaches (parametric as well as nonparamet-
ric methods), where the advantages and limitations of various techniques have been documented.
Among the more promising techniques, non-parametric approaches (such as artificial neu-
ral networks) have been shown to be robust for providing reduced-order, reduced-complexity,
21
computational models. The following section provides an overview of developments in the non-
parametric identification of nonlinear systems typically encountered in the Applied Mechanics
field, and it also motivates the need for capitalizing on some powerful tools that have been devel-
oped in the area of evolutionary optimization to “discover” the governing equation that reflects
the underlying nonlinear phenomena in an observed system, based on its “signature” analysis.
3.1.2 Review of Nonparametric Approaches for the Identification of Nonlinear
Systems
The identification of nonlinear structural systems based on measured input and output data has
been receiving considerable attention in the past decades. Early studies in non-parametric identifi-
cation of nonlinear structural systems from experimental data in the time domain were conducted
by Ibanez [39], and Masri and Caughey [40] back in 1970’s. In the work of Masri and Caughey,
orthogonal Chebyshev polynomials were utilized to identify nonlinear single-degree-of-freedom
(SDOF) systems. In this non-parametric approach, the restoring force of the system is approx-
imated by two-dimensional orthogonal Chebyshev polynomials. Hence, the obtained restoring
force surface (RFS) is an indicator of the nonlinearity-type of the system. This approach has
been applied for identifying a variety of systems including multi-degree-of-freedom (MDOF)
structures [41], buildings [42], nonlinear dampers [43, 44], automotive structures [45] and Mi-
croElectroMechanical Systems (MEMS) [46]. Suggested by Crawley [47] and others [48–50],
utilizing force-state mapping and ordinary polynomials instead of Chebyshev polynomials [51]
were other variations to Masri and Caughey’s RFS method. However, choosing the appropriate
order of polynomial models, in both Chebyshev polynomials and ordinary polynomials, is not
backed-up by rigorous mathematical support. Often, the appropriate order is chosen based on
trial and error that involves experimental investigation and prior knowledge of the system being
22
investigated. Yet, it is not guaranteed that these models do not suffer from over-fitting or under-
fitting. In addition, the discovered models, in general, may not provide information about the
physical nature of the system being investigated.
Another popular approach for system identification is employing Auto Regressive Moving
Average (ARMA) models, first proposed by Box and Jenkins [52]. Typically, the input and output
signals are time histories of the excitation and displacement response, respectively. Suggested
by Leontaritis and Billings [53, 54], this methodology was generalized to be implemented on
nonlinear systems, called Nonlinear Auto Regressive Moving Average with eXogeneous Input
(NARMAX). Auto Regressive models are widely used as damage indicators in civil and mechan-
ical structures, first for being under the class of non-parametric identification, which facilitates the
identification of complex and large structures, and secondly, the identification can be performed
based on the output data without the need for knowing the excitation [55]. Despite the versatility
of the approach and the advantage of exploiting a noise model, not providing physical insight
directly about the phenomenological nature of the system being studied [56] and the non exis-
tence of a rigorous criterion with mathematical support for choosing the right order of models,
are significant disadvantages attributed to this method [57].
Being under the class of black-box modeling, the use of Neural Networks is another approach
for identifying complex nonlinear systems under the category of non-parametric methods. Some-
times, in real world problems, selecting a physical model for complex systems that results in
accurate estimations is not feasible. In that sense, Neural Networks are a very powerful tool for
understanding and describing the behavior of complex systems on the basis of the input and output
data, without possessing any prior knowledge about the nature of the system. Neural Networks
were employed for describing various systems in different disciplines ranging from biology to
engineering. Masri et al. [58, 59] and Worden et al. [60, 61] were amongst the first researchers
23
that utilized the technique for the identification of nonlinear systems in the field of Applied Me-
chanics. Neural Networks were employed for the identification of MDOF systems [62, 63], the
detection of changes in system parameters [64, 65] and modeling nonlinear forces [66, 67]. How-
ever, the negative aspect of describing systems on the basis of data merely by means of Neural
Networks is the lack of a physical model, and consequently not obtaining insight into the nature
of the system being studied.
3.1.3 Scope
Considering the limitations of available methods for the identification of nonlinear complex sys-
tems, the investigation in this chapter explores the potential of evolutionary approaches and adapts
them for the identification of complex systems with polynomial nonlinearities generally encoun-
tered in the field of Applied Mechanics. This chapter introduces a proposed hybrid technique
based on GP with a focus on problems of interest to the Applied Mechanics community, to es-
tablish the foundations of the employment of evolutionary techniques for the identification of
different dynamical systems. The proposed approach combines ordinary Genetic Programming,
Genetic Algorithms (GAs), and Linear Algebra to intelligently find the differential equations that
govern highly nonlinear systems by:
Utilizing GP for finding suitable structures.
Optimizing the parameters of the evolving structures using GAs.
Incorporating Linear Algebra techniques to simplify tree structures, and consequently con-
trol bloat.
It will be shown that applying the approach of this chapter overcomes the aforementioned limita-
tions of conventional methods by: (1) Eliminating the requirement of prior knowledge, or intuitive
24
assumptions about the system up-front, and (2) Providing insight into the nature of the problem,
and unveiling the physics of the system.
This chapter introduces a hybrid Genetic Programming-based technique to “discover” a par-
simonious differential operator that represents an optimum match to the governing differential
equation of the target underlying complex nonlinear system. The material in this chapter is orga-
nized in the following sequence: Section 3.2 presents the evolutionary based optimization scheme
for discovering the governing differential equations of nonlinear systems; Sections 3.3, 3.4, 3.5,
3.6 present the results of the applied technique for the identification of Linear, Duffing, van der
Pol and noisy Duffing-van der Pol oscillators, respectively; Section 3.7 provides the discussion,
and Section 3.8 summarizes this chapters.
3.2 Identification Scheme Composed of GP, LA and GA
The generic form of the governing ODE of the class of nonlinear SDOF systems being studied in
this chapter is:
f(t) =m x +r(x; _ x) +p(t)x (3.1)
wheref(t) is the external excitation,m is the mass of the system,r(x; _ x) is the restoring force,x,
_ x and x represent the displacement, velocity, and acceleration of the system, respectively, andp(t)
is the external disturbance that appears in the noisy Duffing-van der Pol oscillator ODE only, and
otherwise is 0 for other classes of nonlinear systems. The external excitation and the acceleration
of the system are assumed to be available from measurements. One may obtain the velocity and
displacement by careful integration of the acceleration, where strategies for obtaining independent
states and noise-accumulation prevention should be placed. Therefore, the identification is carried
out with the knowledge of the exciting forcef(t) and all three states, displacementx, velocity _ x
25
and acceleration x, for all the oscillators with the addition of known disturbancep(t) in the case
of the noisy Duffing-van der Pol oscillator. Systems that are investigated in this chapter include
linear, van der Pol and noisy Duffing-van der Pol oscillators.
The proposed identification scheme, whose flowchart is presented in Fig. (3.1), requires two
fundamental steps. First, the training phase, where the differential operator of the system is re-
vealed after the termination of the evolution process. Second, the validation phase, where the
generalizability of the discovered ODE is evaluated with unseen data. The first step in the train-
ing phase is to generate the reference excitation and the corresponding response of the system.
In order to provide a more realistic simulation,x, _ x and x are polluted by 1% and 5 % additive
noise, i.e., zero-mean Gaussian stationary white noise. Next, the training data set is provided to
the GP module.
In the GP module, the first population is constructed by selecting randomly-generated num-
bers, available state variables and mathematical operators at random. There are two different
approaches for constructing the initial trees of the population: the “full method” and the “grow
method”. The full method produces full trees for the first population, while in the grow method,
unequal length of branches is allowed [13]. Applied in this investigation, the ramped half-and-
half technique produces half of the population using the full method and the rest using the grow
method [13]. Each tree represents a mathematical expression whose intent is to truly map states to
the target signal, i.e., either the external excitation or the restoring force. In other words, the target
signal acts similar to a function whose variables and parameters are to be determined in order to
minimize the error between the target signal and the GP estimate. Therefore every individual in
the population is a candidate solution to represent the entire governing equation of motion, or the
restoring force of the system. Three introduced statesx, _ x and x, and the external disturbance
26
p(t), denoted asX1,X2,X3 andX4, respectively, construct four-variable tree structures.
After the creation of the new generation, individuals are evaluated and a fitness value is as-
signed to every tree. Then, the genetic materials of the population are manipulated by the genetic
operators. The mathematical form of the new trees are simplified by Linear Algebra rules using
the Matlab symbolic toolbox. Then, the parameters of these reduced complexity models are op-
timized by means of Genetic Algorithms, using a population of 30. The parameter optimization
algorithm stops when change in the best fitness value, over 50 consecutive generations, doesn’t
exceedTolFun = 0:01. This value is obtained by experimental study of different scenarios to
ensure accurate results, while the computational cost is affordable. The optimized individuals
are transformed to the expressions trees, their fitness values are updated, and a new population is
formed.
The evolution process terminates when the fitness of the best individual doesn’t improve more
than 0.1% in 10
5
consecutive generations. If unsure about the termination point, one may exam-
ine the fittest individual throughout the course of the optimization process several times through
comparing its solution to the actual response of the system. If an acceptable solution is found,
the evolution process can be stopped. Upon reaching the ultimate answer, the solution of the
discovered ODE is compared to the reference response employed for training initially.
For validation, the discovered ODE is subjected to a new excitation f(t) and disturbance
p(t), that possess similar statistical properties as the training phase, in order to investigate the
performance of the discovered ODE when a not-previously-experienced excitation is applied.
It is noteworthy that any successful identification scheme requires the above methodology to
discover the characteristics of the system in the training phase, and verifies the findings in the
27
validation phase, regardless of the technique that is employed as the core of the identification
procedure (i.e., instead of the genetic programming module used here).
Figure 3.1: Flowchart of the proposed technique for the identification of systems with polynomial-
type nonlinearities.
3.2.1 Excitation
All systems discussed in this chapter are subjected to zero-mean stationary Gaussian white noise
excitation, with known standard deviation

. The duration of the excitation should be sufficiently
long so that the response of the system reaches the stationary phase, where the standard deviation
28
of the response stays almost unchanged. In this investigation, the duration of the excitation record
is such that the response will be approximately 10% in the transition region, and 90% in the
stationary region. It is noteworthy that in system identification, both the transition and stationary
regions have valuable information about the characteristics of the system, and are essential for
accurate inference about the system being investigated. The properties of the excitation used for
the identification of all oscillators are provided in Table 3.1.
Excitation Parameters Values
Mean 0
Gaussian stationary white noise Standard Deviation 1.0
Duration (Periods) 200T
Sampling Frequency t=T 0.05
Table 3.1: Properties of the exciting force for both training and validation, used in the linear,
Duffing, van der Pol and noisy Duffing-van der Pol oscillators.
3.2.2 Noise
In the real world, data are acquired from sensors mounted on systems and transferred to the
processing units. Yet, data acquisition and transmission devices may introduce pollution to the
data during the operation process. Thus, in order to determine a model, obtained as a consequence
of a realistic manner, incorporating noise in the data is essential. The additive noise associated
with the acquired data is assumed to be zero-mean stationary Gaussian white, with a certain
standard deviation. The standard deviation (SD) of the noise

is assumed to be a fraction of the
SD of the signal
X
as follows:
29

X
=X
RMS
=
v
u
u
t
1
N
N
X
i=1:n
x
2
i
(3.2)

=
'
100

X
(3.3)
whereN is the number of data pointsx
i
, and' is the percentage of the introduced noise with
respect to the signal amplitudes.
3.3 Identification of Linear Oscillator
For “calibration” purposes, GP is employed for the identification of the governing equation of
motion of a linear oscillator, denoted as follows:
f(t) =m x +c _ x +kx (3.4)
wherem is the mass,c is the damping factor,k is the stiffness,f(t) is the exciting force, and x, _ x
andx are the acceleration, velocity and displacement of the system, respectively. The properties
of the exciting forcef(t) and the parameters of the system are listed in Table (3.1) and Table (3.2),
respectively. The exciting force and the polluted states obtained from the numerically solved ODE
of the system with 1% and 5% additive noise are employed for training. As the standard deviation
of the response in Fig. (3.3) suggests, the excitation is lasting a sufficient time to ensure that the
oscillator has reached the steady state.
Actual and discovered equations of motion, and the corresponding tree structures are pre-
sented in Table (3.3) and Fig. (3.2), respectively. The associated fitness value is based upon the
30
Actual Identified
1% noise 5% noise
Values Values error Values error
m 2.000 1.998 0.1% 1.964 2.0%
c 0.080 0.080 0.0% 0.078 2.5%
k 0.320 0.319 0.3% 0.313 2.2%
Frequency!
n
(Hz) 0.400 0.400 0.0% 0.399 0.2%
Damping Ratio 0.050 0.050 0.0% 0.049 0.5%
Table 3.2: Parameters of the linear oscillator: actual vs. identified.
normalized error between the actual and the resultant exciting force. Actual and extracted pa-
rameters from the discovered ODE that are displayed alongside in Table (3.2) suggest that GP is
able to discover both the form of the governing equation, and the parameters of the system, fairly
accurately. Still, as the noise increases, the accuracy of the identified parameters declines, but no
extra term appears in the ODE, as a positive consequence of using a long duration data set.
In furtherance of the investigation, to demonstrate the validity of the identified ODEs, uti-
lizing the excitation of the training phase, the attained differential equations in Table (3.3) are
solved numerically and the solution for the finer range of 70T - 80T is compared to the reference
response used for training in Fig. (3.4) . The associated normalized error for the entire range
of 200T is reported in Table (3.4). In addition, plots of the restoring force and velocity against
displacement are presented in Fig. (3.4(f)) and Fig. (3.4(e)) to enhance the data interpretation by
bringing insight into the problem.
31
To broaden the study and examine the generalization extent of the equations to unseen data,
possessing the same properties as the exciting force used for training, a new white Gaussian
random excitation is produced by utilizing a different seed number. Only the first 80 periods of
the new excitation are used for the validation. The obtained differential equations in Table (3.3)
are subjected to the validation excitation, and solved numerically. The response for periods of
70T - 80T is presented in Fig. (3.5) and the corresponding error for the entire range of 0T - 80T
is tabulated in Table (3.4). The results of the investigation suggest that GP not only is capable of
finding the form of the differential equation and the key parameters of the system, irrespective of
the data provided for training, but also will lead to the accurate estimation of the states, even when
the system is subjected to a not-previously-experienced excitation in a new dynamical condition.
Type Differential Equations Fitness Error
Actual F (t) = 2:000 x + 0:080 _ x + 0:320x -
1% noise F (t) = 1:998 x + 0:080 _ x + 0:319x 2.37%
5% noise F (t) = 1:964 x + 0:078 _ x + 0:313x 11.65%
Table 3.3: The actual and discovered equations of motion of the linear oscillator and their associ-
ated fitness error for different noise levels.
32
States Training Data Validation Data
1% noise 5% noise 1% noise 5% noise
Displacementx 3.17% 3.40% 3.21% 3.47%
Velocity _ x 3.83% 3.89% 3.82% 3.92%
Acceleration x 2.72% 2.61% 2.73% 2.64%
Restoring Forceg(x; _ x) 3.19% 2.19% 2.23% 2.22%
Table 3.4: The normalized error associated with states in the training and validation phase due to
solving the discovered ODEs corresponding to the linear oscillator.
0.320 X1
times
0.080 X2
times
plus
2.000 X3
times
plus
(a) Actual equation of motion
0.313 X1
times
0.078 X2
times
plus
1.964 X3
times
plus
(b) Equation of motion with 5% noise
Figure 3.2: The trees of the linear oscillator’s governing equation of motion.
33
0 100 200
−4
0
4
Excitation
t/T
(a) Excitation
0 50 100 150 200
0
5
Standard Deviation of Displacement
t/T
(b) Standard deviation of displacement
0 50 100 150 200
0
2
Standard Deviation of Velocity
t/T
(c) Standard deviation of velocity
0 50 100 150 200
0
1
Standard Deviation of Acceleration
t/T
(d) Standard deviation of acceleration
Figure 3.3: The external force and statistics of the linear oscillator response used for training.
34
70 72 74 76 78 80
−10
0
10
Displacement
t/T

0% noise
1% noise
5% noise
(a) Time history of displacement
70 72 74 76 78 80
−3
0
3
Velocity
t/T

0% noise
1% noise
5% noise
(b) Time history of velocity
70 72 74 76 78 80
−2.5
2.5
Acceleration
t/T

0% noise
1% noise
5% noise
(c) Time history of acceleration
70 72 74 76 78 80
−3
0
3
Restoring Force
t/T

0% noise
1% noise
5% noise
(d) Time history of restoring force
−10 0 10
−3
0
3
Restoring Force
Displacement

0% noise
1% noise
5% noise
(e) Restoring force vs. Displacement
−10 0 10
−3
0
3
Velocity
Displacement

0% noise
1% noise
5% noise
(f) Velocity vs. Displacement
Figure 3.4: Comparison of the actual and identified response of the linear oscillator to the exciting
force used for training.
35
0 80
−4
0
4
Excitation
t/T
(a) Validation excitation
70 72 74 76 78 80
−10
0
10
Displacement
t/T

0% noise
1% noise
5% noise
(b) Time history of displacement
70 72 74 76 78 80
−3
0
3
Velocity
t/T

0% noise
1% noise
5% noise
(c) Time history of velocity
70 72 74 76 78 80
−2.5
2.5
Acceleration
t/T

0% noise
1% noise
5% noise
(d) Time history of acceleration
70 72 74 76 78 80
−3
0
3
Restoring Force
t/T

0% noise
1% noise
5% noise
(e) Time history of restoring force
−10 0 10
−3
0
3
Restoring Force
Displacement

0% noise
1% noise
5% noise
(f) Restoring force vs. Displacement
−10 0 10
−3
0
3
Velocity
Displacement

0% noise
1% noise
5% noise
(g) Velocity vs. Displacement
Figure 3.5: Comparison of the actual and identified response of the linear oscillator to a new
exciting force for validation.
36
3.4 Duffing Oscillator
The proposed method is applied for the identification of the Duffing oscillator, defined as follows:
f(t) =m x +c _ x +k(x +x
3
) (3.5)
The parameters of the Duffing oscillator and the properties of the exciting force are the same as
the linear oscillator with the addition of the cubic stiffness nonlinearity term, and the correspond-
ing parameter. The same identification procedure discussed above is applied and the ODEs in
Table (3.6), the corresponding trees in Fig. (3.6) and the system parameters in Table (3.5) are
obtained. The reference and identified responses are also compared in Fig. (3.7) and Fig. (3.8),
respectively, and the normalized error is reported in Table (3.7).
Since the importance of any identification scheme is to provide reliable tools to infer about
the behavior of the system to a previously unexperienced external stimulation, validation results
more effectively express the robustness of the approach. Hence, a new exciting force is used for
the validation. The reference response and the response estimates of the system stimulated by
the validation excitation are shown in Fig. (3.9), and the associated error is reported in Table
(3.7). Acceptable match of the reference response and its estimate implies the effectiveness of the
proposed approach for the identification of the Duffing oscillator.
37
Parameters Actual Identified
1% noise 5% noise
Values Values error Values error
m 2.000 1.996 0.2% 1.918 4.1%
c 0.080 0.080 0.0% 0.075 6.3%
k 0.320 0.320 0.0% 0.309 3.4%
0.060 0.060 0.0% 0.058 3.3%
Initial Frequency!
n
(Hz) 0.400 0.400 0.0% 0.401 0.3%
Damping Ratio 5.0% 5.0% 0.0% 4.9% 2.49%
Table 3.5: Parameters of the Duffing oscillator: actual vs. identified by GP.
Type Differential Equations Fitness Error
Actual F (t) = 2:000 x + 0:080 _ x + 0:320x + 0:019x
3
-
1% noise F (t) = 1:996 x + 0:080 _ x + 0:320x + 0:019x
3
3.76%
5% noise F (t) = 1:918 x + 0:075 _ x + 0:309x + 0:018x
3
18.23%
Table 3.6: The actual and discovered equations of motion of the Duffing oscillator and their
associated fitness for different noise levels.
38
States Training Data Validation Data
1% noise 5% noise 1% noise 5% noise
Displacementx 0.7% 4.7% 0.8% 4.9%
Velocity _ x 0.8% 6.3% 0.8% 6.5%
Acceleration x 0.9% 2.8% 0.8% 3.3%
Restoring Forcer(x; _ x) 2.4% 5.3% 2.7% 5.3%
Table 3.7: The normalized error associated with states in training and validation phase due to
solving the GP-found ODEs.
39
0.019
X1
X1 X1
times
times
times
0.320 X1
times
plus
0.080 X2
times
plus
2.000 X3
times
plus
(a) Actual equation of motion
0.018
X1
X1 X1
times
times
times
0.309 X1
times
plus
0.075 X2
times
plus
1.918 X3
times
plus
(b) Equation of motion with 5% noise
Figure 3.6: The trees of the Duffing oscillator’s governing equation of motion.
40
0 100 200
−4
0
4
Excitation
t/T
(a) Excitation
0 50 100 150 200
0
3
Standard Deviation of Displacement
t/T
(b) Standard deviation of displacement
0 50 100 150 200
0
2
Standard Deviation of Velocity
t/T
(c) Standard deviation of velocity
0 50 100 150 200
0
1
Standard Deviation of Acceleration
t/T
(d) Standard deviation of acceleration
Figure 3.7: The external force and statistics of the Duffing oscillator response used for training.
41
70 72 74 76 78 80
−6
0
6
Displacement
t/T

0% noise
1% noise
5% noise
(a) Time history of displacement
70 72 74 76 78 80
−4
0
4
Velocity
t/T

0% noise
1% noise
5% noise
(b) Time history of velocity
70 72 74 76 78 80
−3
3
Acceleration
t/T

0% noise
1% noise
5% noise
(c) Time history of acceleration
70 72 74 76 78 80
−6
0
6
Restoring Force
t/T

0% noise
1% noise
5% noise
(d) Time history of restoring force
−6 0 6
−6
0
6
Restoring Force
Displacement

0% noise
1% noise
5% noise
(e) Restoring force vs. Displacement
−6 0 6
−4
0
4
Velocity
Displacement

0% noise
1% noise
5% noise
(f) Velocity vs. Displacement
Figure 3.8: Comparison of reference and identified response of the Duffing oscillator to the excit-
ing force used for training.
42
0 80
−4
0
4
Excitation
t/T
(a) Validation Excitation
70 72 74 76 78 80
−6
0
6
Displacement
t/T

0% noise
1% noise
5% noise
(b) Time history of displacement
70 72 74 76 78 80
−4
0
4
Velocity
t/T

0% noise
1% noise
5% noise
(c) Time history of velocity
70 72 74 76 78 80
−3
3
Acceleration
t/T

0% noise
1% noise
5% noise
(d) Time history of acceleration
70 72 74 76 78 80
−6
0
6
Restoring Force
t/T

0% noise
1% noise
5% noise
(e) Time history of restoring force
−6 0 6
−6
0
6
Restoring Force
Displacement

0% noise
1% noise
5% noise
(f) Restoring force vs. Displacement
−6 0 6
−4
0
4
Velocity
Displacement

0% noise
1% noise
5% noise
(g) Velocity vs. Displacement
Figure 3.9: Comparison of reference and identified response of the Duffing oscillator to a new
exciting force for validation.
43
3.5 Identification of van der Pol Oscillator
The aforementioned identification technique is utilized again for the identification of the van der
Pol oscillator, characterized by:
f(t) =m x(1x
2
) _ x +kx (3.6)
The parameters of the van der Pol oscillator, and the properties of the exciting force are similar
to the linear oscillator with the exception of the presence of the velocity-displacement nonlinear-
ity term, and the corresponding coefficient. The characteristics of the exciting force and the key
parameters of the system are presented in Table (3.1) and Table (3.8), respectively. The identi-
fication procedure presented earlier in this chapter is applied to the synthetic data obtained from
the van der Pol oscillator for training. The discovered ODEs are presented in Table (3.9), the
corresponding tree structures are shown in Fig. (3.10), and the extracted parameters of the system
form the discovered equations are tabulated in Table (3.8). The actual response of the system is
also compared to the model estimated in Fig. (3.11) and Fig. (3.12).
The validation excitation is applied to assess the admissibility of the obtained models using
the proposed evolutionary technique. The reference response is compared to the modal solution
in Fig. (3.13), and the error is reported in Table (3.10). The agreement of the discovered ODEs
and the actual state equation as well as the perfect match of the subsequent numerical solutions to
the reference response of the system verify the strength of the proposed approach for the identifi-
cation of the van der Pol oscillator.
44
Actual Identified
1% noise 5% noise
Values Values error Values error
m 2.000 1.996 0.2% 1.916 4.2%
k 1.000 0.998 0.2% 0.953 4.7%
0.2 0.199 0.5% 0.188 6.0%
Table 3.8: Parameters of the van der Pol oscillator: actual vs. identified.
Type Differential Equations Fitness Error
Actual F (t) = 2:000 x + 0:200 _ xx
2
0:200 _ x + 1:000x -
1% noise F (t) = 1:996 x + 0:199 _ xx
2
0:200 _ x + 0:998x 3.77%
5% noise F (t) = 1:916 x + 0:188 _ xx
2
0:189 _ x + 0:953x 18.29%
Table 3.9: The actual and discovered equations of motion of the van der Pol oscillator and their
associated fitness for different noise levels.
States Training Data Validation Data
1% noise 5% noise 1% noise 5% noise
Displacementx 0.3% 8.3% 0.3% 8.3%
Velocity _ x 0.3% 8.3% 0.3% 8.3%
Acceleration x 16.0% 18.0% 16.0% 18.0%
Restoring Forceg(x; _ x) 0.4% 9.1% 0.4% 9.0%
Table 3.10: The normalized error associated with states in the training and validation phase due
to solving the GP-found ODEs corresponding to the van der Pol oscillator.
45
1.000 X1
times
0.200 X2
times
minus
2.000 X3
times
plus
0.200
X1 X1
times X2
times
times
plus
(a) Actual equation of motion
0.953 X1
times
0.188 X2
times
minus
1.916 X3
times
plus
0.199
X1 X1
times X2
times
times
plus
(b) Equation of motion with 5% noise
Figure 3.10: The tree structures of the van der Pol oscillator’s governing equation of motion.
46
0 100 200
−4
0
4
Excitation
t/T
(a) Excitation
0 50 100 150 200
0
2
Standard Deviation of Displacement
t/T
(b) Standard deviation of displacement
0 50 100 150 200
0
1.5
Standard Deviation of Velocity
t/T
(c) Standard deviation of velocity
0 50 100 150 200
0
1
Standard Deviation of Acceleration
t/T
(d) Standard deviation of acceleration
Figure 3.11: The external force and statistics of the van der Pol oscillator response used for
training.
47
70 72 74 76 78 80
−4
0
4
Displacement
t/T

0% noise
1% noise
5% noise
(a) Time history of displacement
70 72 74 76 78 80
−3
0
3
Velocity
t/T

0% noise
1% noise
5% noise
(b) Time history of velocity
70 72 74 76 78 80
−3
3
Acceleration
t/T

0% noise
1% noise
5% noise
(c) Time history of acceleration
70 72 74 76 78 80
−6
0
6
Restoring Force
t/T

0% noise
1% noise
5% noise
(d) Time history of restoring force
−4 0 4
−6
0
6
Restoring Force
Displacement

0% noise
1% noise
5% noise
(e) Restoring force vs. Displacement
−4 0 4
−3
0
3
Velocity
Displacement

0% noise
1% noise
5% noise
(f) Velocity vs. Displacement
Figure 3.12: Comparison of the reference and identified response of the van der Pol oscillator to
the exciting force used for training.
48
0 80
−4
0
4
Excitation
t/T
(a) Validation excitation
70 72 74 76 78 80
−4
0
4
Displacement
t/T

0% noise
1% noise
5% noise
(b) Time history of displacement
70 72 74 76 78 80
−3
0
3
Velocity
t/T

0% noise
1% noise
5% noise
(c) Time history of velocity
70 72 74 76 78 80
−3
3
Acceleration
t/T

0% noise
1% noise
5% noise
(d) Time history of acceleration
70 72 74 76 78 80
−6
0
6
Restoring Force
t/T

0% noise
1% noise
5% noise
(e) Time history of restoring force
−4 0 4
−6
0
6
Restoring Force
Displacement

0% noise
1% noise
5% noise
(f) Restoring force vs. Displacement
−4 0 4
−3
0
3
Velocity
Displacement

0% noise
1% noise
5% noise
(g) Velocity vs. Displacement
Figure 3.13: Comparison of the reference and identified response of the van der Pol oscillator to
a new exciting force for validation.
49
3.6 Identification of Noisy Duffing-van der Pol Oscillator
The Noisy Duffing-van der Pol oscillator is composed of various linear and nonlinear terms:
velocity-dependent viscous damping, linear stiffness, cubic stiffness nonlinearity, displacement-
velocity dependent van der Pol nonlinearity, and a displacement dependent term that incorporates
an additional disturbancep(t) [68]. Therefore, accomplishing the identification of this complex
nonlinear system guarantees the success of the proposed identification scheme if simpler systems,
such as the linear, Duffing and van der Pol oscillators, are to be identified. GP is employed for
the identification of the governing equation of motion of the noisy Duffing-van der Pol oscillator,
defined as:
f(t) =m x +c _ x +kx +x
3
+x
2
_ x +p(t)x (3.7)
wheref(t) is the primary exciting force, m is the mass, c is the viscous damping factor, k
is the linear stiffness, is the cubic stiffness nonlinearity, is velocity-displacement dependent
nonlinearity coefficient, p(t) is the external disturbance, and x, _ x and x are the acceleration,
velocity and displacement of the system, respectively. In addition to the states, introduction of the
additional disturbancep(t) for training, adds one more variable to the equations. The properties
of the exciting forcef(t) and the parameters of the system are listed in Table (3.11) and Table
(3.12), respectively. The primary excitationf(t) and the disturbancep(t) have similar statistical
properties, but the standard deviation of the uncorrelated disturbancep(t) is assumed to be 30%
of the primary excitationf(t).
The exciting force and the polluted states, obtained from the numerically solved ODE of the
system with 1% and 5% additive noise, are employed for training. As the standard deviation of
50
the displacement in Fig. (3.15) shows, the excitation is lasting enough to ensure that the oscillator
has reached the stationary phase of the motion.
Excitation Parameters Values Evolutionary Parameters Values
Mean 0 Population Size 100
Standard Deviation 1.00 Crossover Rate 0.80
Duration (Periods) 200T Mutation Rate 0.20
Sampling Freq. t=T 0.05 Representation 10%
Basis Functions: Frequency of Simplification
step; abs; times; plus; minus and Parameter Optimization 50 gen.
Table 3.11: Parameters of the exciting force for both training and validation and the evolutionary
parameters of the noisy Duffing-van der Pol.
Actual and identified equations of motion, the extracted parameters, and the tree structure of
the identified ODE, when polluted data with 5% noise are used, are presented in Table (3.13),
Table (3.12) and Fig. (3.14), respectively. The numbers are reported with three decimal places
to show the strength of the approach in optimizing the parameters. The outcomes suggest that
the proposed hybrid technique, is able to reveal both the form of the governing equation and the
parameters of the system, reasonably accurately. Still, as the noise increases, the accuracy of the
identified parameters decreases, but no excessive terms appear in the ODE, as the result of the
intelligent selection of basis functions, and the use of a long-duration data set.
In furtherance of the investigation, and to prove the validity of the identified ODEs, the at-
tained differential equations in Table (3.13) are solved, and the solution for the shorter range of
51
70T – 80T is compared to the reference response used for training in Fig. (3.16) . The associated
normalized error for the entire range of 200T is reported in Table (3.14). In addition, plots of the
restoring force and velocity against displacement are presented in Fig. (3.16(e)) and Fig. (3.16(f)).
To examine the generalization extent of the equations to new excitations in future situations,
possessing the same statistical properties as that of the training, the validation excitation is pro-
duced by a different random seed number. Only the first 80 periods of the new excitation are
utilized for the validation and the actual response is compared to the solution of the GP-obtained
ODEs. Thus the obtained differential equations in Table (3.13) are subjected to the new excita-
tion. The response for periods of 70T – 80T is presented in Fig. (3.17) and the corresponding
error for the entire range of 0T – 80T is listed in Table (3.14).
The results of the investigation in this section suggest that, in addition to the accurate esti-
mation of the states even when the training is conducted by noise-polluted data, the proposed
method results in the identification of the correct form of the governing differential equation and
the key parameters of the system. The exclusion of extraneous arguments in the identified state
equations, when noise is present, is another significant achievement acquired through GP, while
other approaches such as the conventional least squares method might not perform as well [68].
52
Actual Identified
1% noise 5% noise
Values Values error Values error
m 2.000 1.999 0.1% 1.889 5.6%
c 0.080 0.080 0.0% 0.090 12.5%
k 0.320 0.320 0.0% 0.321 0.3%
0.060 0.060 0.0% 0.049 18.3%
0.200 0.200 0.0% 0.174 20.0%
Table 3.12: Parameters of the noisy Duffing-van der Pol oscillator, actual vs. identified.
Data Differential Equations Fitness Error
Reference f(t) = 2:000 x + 0:200 _ xx
2
+ 0:080 _ x + 0:320x + 0:060x
3
+ 1:000p(t)x -
1% noise f(t) = 1:999 x + 0:200 _ xx
2
+ 0:080 _ x + 0:320x + 0:060x
3
+ 0:980p(t)x 2.2%
5% noise f(t) = 1:889 x + 0:174 _ xx
2
+ 0:090 _ x + 0:321x + 0:049x
3
+ 0:927p(t)x 10.8%
Table 3.13: The actual and discovered differential equations of the noisy Duffing-van der Pol
oscillator and their associated RMS error for different noise levels.
53
States Training Data Validation Data
1% noise 5% noise 1% noise 5% noise
Displacementx 0.1% 5.2% 0.1% 5.0%
Velocity _ x 0.1% 5.6% 0.1% 5.3%
Acceleration x 0.1% 6.0% 0.1% 5.9%
Table 3.14: The normalized error associated with estimated states in training and validation phase
due to solving the GP-found ODEs corresponding to the noisy Duffing-van der Pol oscillator.
1.889 X3
times
0.090 X2
times
plus
0.321 X1
times
plus
0.927
X1 X4
times
times
plus
0.174
X1 X1
times X2
times
times
plus
0.049
X1
X1 X1
times
times
times
plus
Figure 3.14: The tree of the noisy Duffing-van der Pol oscillator’s governing equation of motion
obtained with 5% noise.
54
0 100 200
−4
0
4
Excitation
t/T
(a)
0 50 100 150 200
0
1.5
STD of Displacement
t/T
(b)
Figure 3.15: External force and statistics of the noisy Duffing-van der Pol oscillator response used
for training: (a) excitation; (b) standard deviation of displacement.
55
70 72 74 76 78 80
−4
0
4
Displacement
t/T
(a)
70 72 74 76 78 80
−2.5
0
2.5
Velocity
t/T
(b)
70 72 74 76 78 80
−3
3
Acceleration
t/T
(c)
70 72 74 76 78 80
−4
0
4
Restoring Force
t/T
(d)
−4 0 4
−2.5
0
2.5
Velocity
Displacement
(e)
−4 0 4
−4
0
4
Restoring Force
Displacement

Exact
1% noise
5% noise
(f)
Figure 3.16: Comparison of the actual and identified response of the noisy Duffing-van der Pol
oscillator to the exciting force used for training: (a) time history of displacement; (b) time history
of velocity; (c) time history of acceleration; (d) time history of restoring force; (e) velocity vs.
displacement; (f) restoring force vs. displacement.
56
0 80
−4
0
4
Excitation
t/T
(a)
0 80
−4
0
4
Disturbance
t/T
(b)
70 72 74 76 78 80
−3
3
Acceleration
t/T
(c)
70 72 74 76 78 80
−4
0
4
Restoring Force
t/T
(d)
−4 0 4
−2.5
0
2.5
Velocity
Displacement
(e)
−4 0 4
−4
0
4
Restoring Force
Displacement

Exact
1% noise
5% noise
(f)
Figure 3.17: Comparison of the actual and identified response of the noisy Duffing-van der Pol
oscillator to a new exciting force for validation: (a) validation excitationf(t); (b) external dis-
turbancep(t); (c) time history of acceleration; (d) time history of restoring force; (e) velocity vs.
displacement; (f) restoring force vs. displacement.
3.7 Discussion
The proposed identification scheme is able to accurately discover the differential equations that
govern the behavior of the oscillators under the discussion. Using a persistent excitation is the
57
key to the discovery of correct differential equations, even when noisy data are used. Other-
wise, a short excitation may not be able to effectively stimulate the oscillators to reveal their
nonlinearities. Another drawback of using a short length excitation is over parametrization. The
identification scheme may mistakenly fit a higher order polynomial to the data to accurately fit
the entire span.
Utilizing evolutionary techniques for parameters optimization, resulted in accurate estimates
of the parameters, that was not easily achievable if GP was not coupled with a parameter optimiza-
tion routine. Moreover, the employment of Linear Algebra rules for the symbolic simplification
resulted in parsimonious models, and subsequently enhanced the evolution. The suggested ap-
proach is also capable of disclosing the physical nature of the problem. Careful investigation of
the obtained differential equations results in the discovery of the parameters of the system and
their corresponding values. These parameters may even disclose more information about the sys-
tem. For example the coefficient adjacent to the displacementx is known to be the stiffness of the
systemk. If the identification is performed in multiple stages at different times, the deterioration
of the system can be reflected in the stiffness changes. Therefore, not only is the differential equa-
tion capable of predicting the performance of the system when subjected to a certain stimulation,
but also it can raise a warning flag when extreme changes in the parameters of the system occur.
Hence, the presented methodology is an effective approach for structural health monitoring pur-
poses.
Many of the identification approaches that are based on the investigation of the restoring
force of the system, require an estimate of the mass up-front. However, the proposed technique,
58
presented in this chapter, has the advantage that the mass doesn’t need to be known a priori. More-
over, once the differential equation of the system is discovered, the investigation of the restoring
force using the available approaches in the literature is easily achievable.
3.8 Conclusion
Being an idea rooted in computational intelligence for the evolution of structures aimed to ac-
complish a defined task, this investigation combines Genetic Programming, for discovering suit-
able structures, together with Linear Algebra techniques for symbolic manipulations as well as
Genetic Algorithms for the parameter optimization, to provide a robust methodology for the non-
parametric identification of complex nonlinear systems in the field of Applied Mechanics. The
procedure is based on direct identification of the governing differential equation of the system by
optimizing the choice of various basis functions, including polynomials and discontinuous func-
tions, through the evolution process. No model is postulated up-front, and the program is only
provided with input and output data. Different classes of nonlinear single-degree-of-freedom sys-
tems with polynomial-type nonlinearities are studied to investigate the utility and generalizability
of the proposed approach. It is shown that the proposed method yields accurate estimates of the
response, and is very effective in unveiling the physical nature of the underlying problems.
59
Chapter 4
Identification of SDOF Systems with Non-Polynomial
Nonlinearities
4.1 Introduction
L
AST chapter presented an evolutionary based identification technique and its application
to complex systems exhibiting polynomial-type nonlinearities. It was shown that this
methodology yields accurate differential equations that represent dynamical systems exhibiting
polynomial-type nonlinearities. These models can accurately estimate the response of the sys-
tem in future dynamical situations. However, nonlinearities associated with many phenomena in
nature are not definable by polynomials.
The systems studied in this chapter involve discontinuous and non-differentiable functions in
their governing equation of motion. Consequently, they exhibit sharp jumps in their response,
or have unsmooth corners, yet without hysteretic behavior (i.e., they are single-valued nonlinear
functions). Thus, non-parametric identification of their state equation using polynomials may
not result in accurate estimates. Much research has been conducted on the approximation of
60
discontinuous functions such as thesgn and the Heavisidestep function by means of Pad´ e ap-
proximations [69], and its extensions [70–72]. The Pad´ e approximation approach proposes the
employment of the ratio of two polynomials to estimate the behavior of systems that exhibit dis-
continuity in their response. This approximation can be considered as a generalized Taylor’s
series expansion. However, one of the well-studied issues of Pad´ e approximation is the Gibbs
phenomenon that causes unfavorable oscillation and a significant mismatch of the estimate at un-
smooth points [73]. Much effort has been done to decrease the oscillations through expanding
Pad´ e approximation and filtering [74–76]. However, there still exist limitations tied up with the
approach in general, such as not being able to find the polynomial that has a root at the disconti-
nuity point and no where else, and not being able to satisfy the boundary conditions [76]. Unlike
other approaches that provide a surrogate model to the restoring force surface using GP [77], the
application of the introduced hybrid approach, results in the discovery of the exact model that
represents the system.
The proposed approach of this study is exploited for the identification of the restoring force
r(x; _ x) of systems that exhibit non-polynomial nonlinearities. It is assumed that the states and
the exciting force are available from measurements. Since each state accounts for a dimension in
the search domain, eliminating one state, increases the convergence rate significantly. Therefore,
for complex nonlinear systems in this section, GP is only supplied by x and _ x, denoted as X1
andX2, respectively to reduce the search domain for the identification ofr(x; _ x). The only draw
back of this strategy is the necessity of an estimate on the massm. Hence, beside the availability
of displacement and velocity from measurements, the restoring force should be calculated by
r(x; _ x) = f(t)m x. Thus, GP tries to map the displacement and velocity to the restoring
force through a mathematical expression by minimizing the normalized error between the actual
restoring force and its estimate, and eventually, to discover the differential operator of the system.
61
4.1.1 Scope
Considering the limitations of available methods for the identification of nonlinear complex sys-
tems with discontinuity in their response, this chapter explores the potential of evolutionary ap-
proaches for the identification of complex systems with non-polynomial nonlinearities which ex-
hibit discontinuity in their response. The proposed approach combines ordinary Genetic Program-
ming for evolving structures, Genetic Algorithms (GAs) for parameter optimization, and Linear
Algebra for complexity reduction and bloat control, to intelligently find the differential equations
that characterize the dynamics of the systems. The material in this chapter is organized as fol-
lows: Section 4.2 introduces the proposed identification technique. The flowchart of the proposed
algorithm is also presented in this section; Section 4.3 presents the results of the identification
of the Coulomb friction oscillator; Section 4.4 presents the results of the identification of a Tri-
Linear oscillator; Section 4.5 discusses the computational advantages of the proposed technique
as well as the advantages and disadvantages of this technique over other available identification
techniques, and Section 4.7 presents the conclusion.
4.2 Identification Scheme Composed of GP, LA and GA
The generic form of the governing ODE of the class of nonlinear SDOF systems being studied is:
f(t) =m x +r(x; _ x) (4.1)
wheref(t) is the external excitation,m is the mass of the system,r(x; _ x) is the restoring force,
x, _ x and x represent the displacement, velocity, and acceleration of the system, respectively. The
62
external excitation and the acceleration of the system are assumed to be available from mea-
surements. One may obtain the velocity and displacement by integrating the acceleration, and
avoiding the accumulation of error. Therefore, the identification is performed by knowing the
exciting forcef(t) and all three states, displacementx, velocity _ x and acceleration x. Systems
that are investigated in this chapter include Coulomb friction oscillator and Tri-Linear oscillator
exhibiting non-hysteretic yielding behavior.
The flowchart of the proposed identification scheme, is shown Fig. (4.1) which is identical to
the process elaborated in chapter (3). The parameters of the Genetic Algorithms for the parameter
optimization are also identical.
63
Figure 4.1: Flow chart of the proposed identification scheme.
4.2.1 Excitation
All systems are subjected to zero-mean stationary Gaussian white noise excitation, with known
standard deviation. The duration of the excitation should be sufficiently long so that the systems
yield. The properties of the external excitation will be presented next.
64
4.3 Identification of Coulomb Friction Oscillator
The behavior of the Coulomb friction oscillator is characterized by:
f(t) =m x +f
c
( _ x) +kx; f
c
( _ x) =csgn( _ x) =
8
>
>
>
<
>
>
>
:
c if _ x< 0:
c if _ x> 0:
(4.2)
wheref(t) is the applied force,m is the mass,f
c
( _ x) is the friction force,c is the friction factor,
k is the linear stiffness, and x, _ x and x are the acceleration, velocity, and displacement of the
system, respectively. Shown schematically in Fig. (4.2(a)), the system experiences a sharp jump
in the response, when the direction of motion alters. This system is excited by a white Gaussian
excitation whose properties are specified in Table (4.1).
2.000 X3
times
0.154
X2
sign
times
plus
0.318 X1
times
plus
(a) (b)
Figure 4.2: SDOF oscillator with Coulomb friction: (a) behavior of the friction force of the
Coulomb friction oscillator; (b) the tree of governing equation of motion obtained with 5% noise.
65
The fittest individuals of the population, after the termination of the optimization procedure,
are presented in Table (4.3). Since bothsgn and Heavisidestep functions were placed in the pool
of potential basis functions, they were both found in the tree structure of the identified ODEs,
with estimates that are the same and virtually identical to the exact system, when data pollution is
not significant.
The direct identification of discontinuous functions by GP eliminates any need for the poly-
nomial approximation, and as a consequence, removes any unwanted oscillations about the un-
smooth points that a polynomial approximation normally introduces. The extracted parameters
from the ODEs that involve thesgn function are compared with the corresponding real values in
Table (4.2). Ifsgn and Heavisidestep functions are taken out of the function pool, when noise
free data are used, GP converges to the rational polynomial approximation in Table (4.3), with
inevitable oscillations around the jump as shown in Fig. (4.3) for the interval of _ x: [-12 12].
Since the effectiveness of any identification scheme is essential in predicting the behavior of
the system subjected to not-formerly-experienced external stimulants, validation results are more
valuable and are better reflections of the efficiency of the approach. Therefore, only validation
results will be reported next (i.e., the training results are skipped).
The solution of the ODEs under the validation excitation is compared with the actual response
in Fig. (4.4), and the corresponding error is reported in Table (4.4). The results demonstrate the
effectiveness of the proposed approach in discovering the differential operator for the system that
exhibit discontinuity in its response, without causing the Gibbs phenomenon.
66
Excitation Parameters Values Evolutionary Parameters Values
Mean 0 Population Size 100
Standard Deviation 1.00 Crossover Rate 0.80
Duration (Periods) 18T Mutation Rate 0.20
Sampling Freq. t=T 0.20 Representation 10%
Basis Functions:step; sgn; Frequency of Simplification
abs; times; plus; minus; divide and Parameter Optimization 50 gen.
Table 4.1: Properties of the applied force for both training and validation used in a SDOF system
with Coulomb friction and the corresponding evolutionary parameters.
Actual Identified
1% noise 5% noise
Values Values error Values error
m 2.000 - - - -
c 0.160 0.159 0.6% 0.154 3.8%
k 0.320 0.320 0.0% 0.318 0.6 %
Table 4.2: Properties of the oscillator with Coulomb friction, actual vs. identified. Mass assumed
known.
67
Type Restoring Force Equation Fitness Error
Actual r(x; _ x) = 0:160sgn( _ x) + 0:320x -
1% noise r(x; _ x) = 0:159sgn( _ x) + 0:320x 1.6%
r(x; _ x) = 0:318step( _ x) + 0:320x 0:158 1.6%
5% noise r(x; _ x) = 0:154sgn( _ x) + 0:318x 7.6%
r(x; _ x) = 0:306step( _ x) + 0:0319x 0:150 7.7%
Rational Polynomial r(x; _ x) =
51:212 _ x + 46:460 _ x
3
+ _ x
5
51:212 + 514:174 _ x
2
+ 66:293 _ x
4
+ 0:202 _ x
6
+ 0:320x 7.6%
Table 4.3: The actual and discovered equations of motion of the oscillator with Coulomb friction
and their associated RMS error for different noise levels.
States Training Data Validation Data
1% noise 5% noise 1% noise 5% noise
Displacementx 0.4% 5.0% 0.7% 9.0%
Velocity _ x 0.4% 5.0% 0.7% 9.0%
Acceleration x 0.5% 5.0% 0.7% 9.0%
Restoring Forcer(x; _ x) 0.5% 5.2% 0.7% 9.2%
Table 4.4: The RMS error associated with the estimated states in the training and validation phase
due to solving the GP-found ODEs corresponding to the oscillator with Coulomb friction.
68
−10 −5 0 5 10
−0.2
−0.15
−0.1
−0.05
0
0.05
0.1
0.15
0.2
Friction Force
Velocity

Rational Polynomial
Exact
Figure 4.3: Exact vs. rational polynomials, both obtained from GP.
0 18
−4
0
4
Excitation
t/T
(a)
0 2 4 6 8 10 12 14 16 18
−32
0
32
Displacement
t/T
(b)
0 2 4 6 8 10 12 14 16 18
−12
0
12
Velocity
t/T
(c)
−12 0 12
−11
0
11
Restoring Force
Velocity
(d)
−12 0 12
−0.2
0
0.2
Coulomb Friction
Velocity
(e)
−32 0 32
−11
0
11
Restoring Force
Displacement

Exact
1% noise
5% noise
(f)
Figure 4.4: Comparison of the actual and identified response of the oscillator with Coulomb
friction to a new exciting force for validation: (a) validation excitation; (b) time history of dis-
placement; (c) time history of velocity; (d) restoring force vs. velocity; (e) Coulomb friction vs.
velocity; (f) restoring force vs. displacement. Note that three superposed curves are shown in
each figure.
69
4.4 Identification of Tri-Linear Oscillator
Consisting of three lines joined at two sharp corners, the tri-linear oscillator is an ideal example
to unveil the capability of GP for the identification of complex nonlinear systems. The motion of
this example oscillator, i.e., presenting ideal “elasto-plastic” phenomena is governed by:
f(t) =m x +c _ x +f
s
(x) ; f
s
(x) =
8
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
:
kx
y
; if x
y
>x:
kx; if x
y
>x>x
y
:
kx
y
; if x>x
y
:
(4.3)
wheref(t) is the external excitation,m is the mass of the system,c is the viscous damping,f
s
(x)
is the stiffness force, and x, _ x andx are the acceleration, velocity, and displacement of the system,
respectively. Shown schematically in Fig. (4.5), as a result of the assumed time-invariant piece-
wise linear behavior, the stiffness forcef
s
(x) undergoes two trend change at the yielding points,
without encountering hysteretic loops. The parameters of the system and the applied force are
introduced in Table (4.5). Carrying out the same identification procedure as the oscillator with
Coulomb friction, the training data are noise-polluted displacement and velocity with 1% and 5%
noise as well as the restoring force. The unpolluted data are also employed to obtain insight into
the exact demonstration of the tri-linear behavior by a single mathematical expression.
The fittest individuals of the population for all three noise levels are shown in Table (4.6), and
the tree structure of the exact discovered ODE of the system that corresponds to the unpolluted
data, is shown in Fig. (4.6). Trained by unpolluted data, the exact solution is composed of several
step functions, each acting at a certain interval to create three different segments of the stiffness
force. Thus, careful investigation of the obtained expression leads to extracting the governing
70
multi-criterion differential equation, the range of applicability of each criteria and eventually the
key parameters of the system. Scrutinizing the revealed ODE and the plot of the stiffness force
in Fig. (4.7(i)), one may even rewrite the discovered restoring force equation presented in Table
(4.6) in the parametric form as follows:
g(x; _ x) =c _ x + 2f
y
step(x
y
+x)f
y
step(x
2
y
x
2
) +kxstep(x
2
y
x
2
)f
y
(4.4)
When polluted data are employed, the state equation is difficult to be interpreted physically,
however, it is accurate enough to be employed to approximate the response of the system under
the validation excitation as shown in Fig. (4.7). The associated error presented in Table (4.7)
indicates that the identified ODEs are acceptable, and state estimates are reasonably accurate.
It is also observed that, when the noise is 5%, the general trend of the oscillator is captured
effectively, but at some intervals the differential equation introduces error in the restoring force
estimate as shown in Fig. 4.7(e). However, the states are still estimated accurately, which is more
useful from the identification point of view, in most cases.
The key to the success of this approach is its ability to reveal the physical nature of the prob-
lem in an autonomous manner through intelligent search and direct identification of the functions
involved in the state equation, such as sgn and step. As a result, it removes the need to em-
ploy polynomial approximations and the consequent Gibbs phenomenon that generates erroneous
oscillations at unsmooth points. Moreover, the discovered solution captures the true nature of
the system being investigated, a crucial accomplishment that provides a physical insight into the
problem. This very important attribute distinguishes the identification by GP from similar black
box identification methods, such as Neural Networks, while providing a comparable automatic
scheme.
71
Excitation Parameters Values System Values Evolutionary Parameters Values
Mean 0 m 2.000 Population Size 100
Standard Deviation 1.00 c 0.080 Crossover Rate 0.80
Duration (Periods) 80T k 0.320 Mutation Rate 0.20
Sampling Freq. t=T 0.05 x
y
1.100 Representation 10%
Basis Functions: f
y
0.352 Frequency of Simplification
step; abs; times; plus; minus and Parameter Optimization 50 gen.
Table 4.5: Properties of the tri-linear oscillator, the applied excitations for both training and vali-
dation, and the evolutionary parameters.
Type Restoring Force Fitness Error
0% noise r(x; _ x) = 0:080 _ x + 0:704step(1:100 +x) 0:352step(1:210x
2
)
+0:320xstep(1:210x
2
) 0:352 0.00%
1% noise r(x; _ x) = 0:080 _ x + 0:352step(x) + 0:329xstep(1:000x
2
)
0:352step(xx
3
) 2.22%
5% noise r(x; _ x) = 0:080 _ x + 0:469step(x) 0:234step(0:534xx
3
)
0:117 7.98%
Table 4.6: The actual and discovered equations of the restoring force of the tri-linear oscillator
and their associated RMs error for different noise levels.
72
States Training Data Validation Data
0% noise 1% noise 5% noise 0% noise 1% noise 5% noise
Displacementx 0.2% 1.1% 11.6% 0.2% 1.6% 6.7%
Velocity _ x 0.2% 1.3% 16.2% 0.2% 1.6% 8.5%
Acceleration x 1.0% 1.5% 14.2% 0.2% 1.8% 9.2%
Restoring Forcer(x; _ x) 2.3% 3.8% 34.2% 0.6% 4.6% 23.1%
Table 4.7: The RMS error associated with the estimated states in training and validation phase
due to solving the GP-found ODEs corresponding to the tri-linear oscillator.
Figure 4.5: Behavior of the stiffness force in the tri-linear oscillator.
73
−0.352
0.080 X2
times
plus
0.704
1.100 X1
plus
step
times
plus
0.352
1.210
X1 X1
times
minus
step
times
minus
0.320
X1
1.210
X1 X1
times
minus
step
times
times
plus
Figure 4.6: Tree representation of GP-obtained actual differential operator of the tri-linear oscil-
lator.
74
0 80
−4
0
4
Excitation
t/T
(a)
60 70 80
−40
0
40
Displacement
t/T
(b)
60 70 80
−5.5
0
5.5
Velocity
t/T
(c)
60 70 80
−2
2
Acceleration
t/T
(d)
60 70 80
−1
0
1
Restoring Force
t/T
(e)
−40 0 40
−5.5
0
5.5
Velocity
Displacement
(f)
−40 0 40
−1
0
1
Restoring Force
Displacement
(g)
−5.5 0 5.5
−1
0
1
Restoring Force
Velocity
(h)
−40 0 40
−1
0
1
Stiffness Force
Displacement

Exact
0% noise
1% noise
5% noise
(i)
Figure 4.7: Comparison of the actual and identified response of the tri-linear oscillator to a new
exciting force for validation. (a) time history of validation excitation; (b) time history of displace-
ment; (c) time history of velocity; (d) time history of acceleration; (e) time history of restoring
force; (f) velocity vs. displacement; (g) restoring force vs. displacement;(h) restoring force vs.
velocity; (i) stiffness force vs. displacement. Note that four superposed curves are shown in each
figure.
4.5 Discussion
The computational benefits as well as advantages and disadvantages of the introduced identifica-
tion approach in this chapter are discussed in this section.
75
4.5.1 Computational Advantages of the Introduced Identification Technique
The proposed identification technique is composed of Genetic Programming (GP) as the main op-
timization engine, Genetic Algorithm (GA) for parameter optimization and Linear Algebra (LA)
rules for the symbolic simplification and complexity reduction. In this section the importance of
the inclusion of each technique in the identification process from a computational standpoint is
investigated. Therefore, four scenarios are considered for the identification scheme, and their per-
formances are compared: Using (1) GP only, (2) GP and GA, (3) GP and LA and (4) GP, GA and
LA. Due to the stochastic nature of the applied evolutionary optimization techniques in the iden-
tification procedure, inferring about the overall performance of such processes requires Monte
Carlo simulations. Since evolutionary approaches are extremely computationally expensive, only
10 realizations of the identification of the less complex Coulomb friction oscillator is considered
for the assessment. Each realization is initialized by a different seed number in Matlab. Thus,
the probabilistic evolutionary operations and random number generators are affected, and conse-
quently, result in a different outcome each time the evolution is restarted. However, the same seed
numbers are utilized when the Monte Carlo simulations are resumed for another scenario. The
optimization stops when the fitness of the best individual falls under 0.5 or 5000 generations are
reached. The processing unit has Intel quad core i7-3370 with hyper threading, whose CPU clock
is 3.40 Hz, and 16.0GB RAM. The program was run in parallel on all cores usingmatlabpool
command.
The realizations of all four scenarios are presented in the four columns of Fig. (4.8). The
fitness and the complexity of the best individual in the population during the evolution of 10
realizations and their average are shown for each scenario in the embedded panels of Fig. (4.8).
The fitness value indicates the RMS error between the reference restoring force response and the
76
model response, and the complexity value is the number of nodes in the tree. The first row in Fig.
(4.8) shows the evolution of the best individual’s fitness in the generation for 10 realizations and
the second row shows the complexity trend of the same individual. The first row shows that when
ordinary GP is used, almost half of the realizations do not meet the fitness-stop limit, while when
GA or LA are incorporated in the method, almost all of the realizations reach the fitness-stop
limit. The last panel on the right side of the first row shows how fast all the realizations converge
when both GA and LA are included in the procedure. Studying complexity changes in the second
row reveals that performing the identification only based on ordinary GP, results in very large
complex trees, and consequently, causes bloat. On the other hand, the second and the third panel
on the second row show the significance of GA and LA in decreasing the complexity. Limiting the
complexity and fast convergence of the realizations in the last panel confirm the effectiveness of
the incorporation of both GA and LA in the method under discussion. Comparing the realizations
of GP & GA and GP & LA reveals that fitness improvement is highly dependent on GA while
sharp changes in the complexity trend is mainly affected by LA. In order to show including GA
and LA doesn’t impose more computational cost, that will defeat the purpose, the time required
to complete all 10 realizations in each scenario is shown under the corresponding column. They
suggest that adding GA and LA to GP has computational benefits and the convergence time drops
significantly.
The main graph of Fig. (4.8) compares the evolution of the fitness, averaged over all 10
realizations, for all four scenarios. The fast convergence of the proposed methodology in this
chapter, that integrates GP, GA and LA, confirms that all the aforementioned techniques are re-
quired to provide the best overall performance among all combinations. When noise-free data are
processed, utilizing the optimization technique under discussion guarantees that the optimized
parameter values, unlike randomly generated numbers, provide the best fit to the data, and more
77
1 1000 2000 3000 4000 5000
0.5
1
5
10
20
Generations
Fitness

28974 gen. 8637 gen. 16844 gen. 1914 gen.
344 min. 263 min. 282 min. 132 min.
60% success 90% success 80% success 100% success
GP Only
GP and GA
GP and LA
GP, GA and LA
1 5000
0.5
1
10
20
Fitness
Realizations for GP
1 5000
10
50
100
200
300
Complexity
Generations
1 5000
0.5
1
10
20
Realizations for GP and GA
1 5000
10
50
100
200
300
Generations
1 5000
0.5
1
10
20
Realizations for GP and LA
1 5000
10
50
100
200
300
Generations
1 5000
0.5
1
10
20
Realizations for GP, GA and LA
1 5000
10
50
100
200
300
Generations
Figure 4.8: The comparison of different combinations of the utilized techniques to verify the
superiority of the proposed approach in discovering the restoring force equation of the oscillator
with Coulomb friction
importantly, reflect the physical parameters of the system. In addition, pruning trees by LA dur-
ing noise-free processes results in parsimonious physics-based models that describe the systems
under investigation more accurately. Therefore, incorporating both GA and LA, in addition to the
ordinary GP in the introduced identification approach, is essential and the key to the success of
the method.
78
4.6 Advantages and Limitations of GP Compared to other Identifi-
cation Methods
4.6.1 Advantages
The major benefits of the proposed method relative to other conventional non-parametric ap-
proaches for the identification of complex non-linear systems can be summarized as follows:
Polynomial based methods and NARMAX models
Intuitive assumptions, such as setting the maximum order of the polynomial or the maxi-
mum number or delays of NARMAX models to characterize the estimator, are not required
when GP is employed. In addition, when discontinuous functions are involved, polynomial
estimators cause unfavorable oscillations, and hence, introduce a significant error. Whereas
in the introduced hybrid approach, the employment of non-differentiable and discontinu-
ous basis functions for evolving structures not only does eliminate the need to erroneous
polynomial approximations, but also discloses the underlying physical traits of the system.
Neural Networks
A Neural Network uses sigmoid functions, and adjusts their parameters through filters, to
learn the behavior of the system and provide predictions on the basis of data only. However,
Neural Networks perform completely blindly and have the disadvantage of not providing
any information about the underlying physics of the system. Identification by GP, on the
other hand, is operated on a comparable automatic routine since the user only needs to enter
the data to the program. Yet, the outcome is the governing equation of the system that, in
many cases, unveils the physical nature and the key parameters associated with the system
under investigation.
79
4.6.2 Limitations
The convergence rate in the GP approach is influenced by various factors. The complexity of
the state equation (in the sense that the tree structure requires many elements), involvement of
discontinuous functions, high dimensionality of the search domain, and excessive functions in the
pool, are among the major factors that impose considerable delay in finding the solution. On the
other hand, if the pool of potential basis functions does not include the required building blocks,
the outcome will be only composed of the available functions, that may not be sufficiently accu-
rate, and do not reflect the physics of the problem, similar to the situation encountered with the
polynomial approximation of a discontinuous function in the case of Coulomb friction. More-
over, to capture the behavior of the system over the entire range of its motion, the training data
set needs to be rich and persistent, in terms of the characteristics and duration, respectively. For
instance, in systems with yielding behavior, low amplitude or short excitations may keep the yield
region hidden, in addition to causing over-fitting. Also, highly noise-polluted data may result in
significant distortion of the identified ODE, despite acceptable performance, and consequently,
inference about the physics of the problem may not be feasible.
4.7 Summary and Conclusions
By building on advancements in the field of Computational Intelligence, this chapter presents the
development and application of a hybrid approach that combines Genetic Programming, Com-
puter Algebra and Genetic Algorithms, to “discover” a parsimonious differential operator with
optimized parameters that represent an optimum match to the governing differential equation of
the target underlying complex nonlinear system (and even an exact solution when “clean” or
slightly polluted data are provided). Discontinuous basis functions, such as the Heaviside step
80
function, are included in the differential equations to capture the discontinuity or unsmooth cor-
ners of the response. No model is postulated up-front, and the identification scheme only requires
input and output data.
Different classes of nonlinear single-degree-of-freedom systems, spanning the range from
polynomial-type nonlinearity, to a system with discontinuity in its response, to piece-wise linear
behavior, are studied to investigate the utility and generalizability of the proposed approach. It is
shown that the method yields accurate estimates of the response, and is very effective in unveiling
the physical nature of the underlying process. It is also shown that the inclusion of discontinuous
basis functions removes the need for inaccurate polynomial approximations and the subsequent
Gibbs phenomenon. The advantages and limitations of the suggested approach compared to other
model-free identification schemes are evaluated and discussed.
81
Chapter 5
Identification of SDOF Systems with Hysteretic Behavior
5.1 Introduction
5.1.1 Background
S
IGNIFICANT effort has been devoted over the past decades to develop new data process-
ing algorithms that incorporate various techniques in machine learning, data mining, pattern
recognition, etc., to obtain useful information from large data sets acquired from complex non-
linear dynamical systems. However, despite the current advancements, a scarcity of useful tools
exists in the field of Applied Mechanics for the identification of hysteretic systems. Modeling
and identification of hysteretic dynamical systems on the basis of measured input and output have
been attracting significant attention over a along period of time. Hysteretic behavior is observed
in many fields such as mechanics, structures, electricity, magnetism, materials, among others.
Elasto-plasticity of the steel, the behavior of buildings under strong earthquakes, and the re-
sponse of aerospace systems incorporating joints are examples of systems exhibiting this history-
dependant phenomena. Estimators of such systems cannot be accurately founded merely based
on time-invariant expressions that involve instantaneous values of the state variables. Hence, the
82
consideration of the history of the evolving response should be taken into account in any valid
model of such systems.
5.1.2 Review of Previous Work
Noteworthy research on the identification of systems with memory-dependant characteristics in
the field of Applied Mechanics has been conducted by [40, 48, 78–93]. Bilinear hysteretsis is one
of the well-studied models that characterizes the behavior of many memory-dependent phenom-
ena. Much research has been conducted to formulate and identify the idealized bilinear hysteretic
model. These efforts include providing the exact solution of the system when subjected to a cer-
tain type of excitation [83], Equivalent Linearization [81, 94], Eigenvalues of Linearized Poincare
Map [95, 96], Extended Kalman Filter [97], Stochastic Averaging Technique [98], etc.
Within this context, first introduced by Bouc [99] in 1967, and then popularized by Wen [80]
in 1976, the Bouc-Wen model is a smooth model extensively utilized to characterize a history-
dependent phenomena based on a system of differential equations. This model incorporates a
separate differential equation for the restoring force to incorporate the nonlinear system’s mem-
ory. Over the past few decades, the parameter identification of hysteretic systems was extensively
focused on the identification of Bouc-Wen model, rather than the bilinear hysteretic oscillator —
the latter is governed by algebraic equations. Past studies on the Bouc-Wen model resulted in
the development of various techniques for tackling the problem, including algorithms based on
Gauss-Newton [88, 89], least squares estimation [68, 100, 101], Simplex [102, 103], Levenberg-
Marquardt [103, 104], Kalman Filters [103, 105–107], reduced gradient method [103], Genetic
Algorithms [108, 109], Differential Evolution [110], adaptive laws [101, 108], Neural Networks
[101, 111], etc. However, there still exists a scarcity of effective non-parametric approaches for
83
the identification of hysteretic systems.
5.1.3 Developments in Computational Intelligence Approaches
Over the past few decades, fast computing capabilities and computational intelligence approaches
have opened new horizons in data analysis, modeling and solving the questions that have been
challenging researchers for many years. One approach for tackling identification problems is
to intelligently search the entire domain of feasible solutions by employing optimization methods
such as gradient-based techniques [112, 113] and metaheuristic algorithms [114, 115]. Evolution-
ary algorithms, based on Darwin’s theory of natural selection [3], employ genetic operators, to
evolve a population of potential candidates to reach the optimal solution. Being a special variant
of Genetic Algorithms, Genetic Programming (GP) evolves a population of computer programs,
expressed in tree structures, with respect to a measure of quality. First introduced by John Koza in
1992 [13], GP has been utilized in various fields such as data analysis [18], control of mechanical
systems [22], music [16], stock selection [20], transportation [19], decision making [13], etc.
The aforementioned techniques for the identification of hysteretic systems are mainly based
on employing a model-class corresponding to the Bouc-Wen model, and the subsequent effort to
identify the parameters of the model, or monitor their variations. However, the non-parametric
evolutionary approach introduced for the identification of nonlinear systems in Chapter 4 can be
modified to yield accurate models and subsequent estimates for hysteretic systems under the dis-
cussion. This chapter shows that utilizing the non-parametric hybrid approach introduced in this
chapter, without postulating any model up-front, yields a parsimonious representation in the form
of systems of differential equations that uses input-output data only for the systems characterized
84
by the Bouc-Wen model and other models of hysteretic behavior such as the bilinear oscillator.
This approach adapts GP to discover a desirable structure and combines it with Genetic Algo-
rithms to optimize the numeric parameters of the tree structures during the course of evolution,
and consequently, to improve the fitness of the population. This hybrid technique also incorpo-
rates Linear Algebra to simplify evolving expressions symbolically, and subsequently to decrease
the complexity and control bloat. It will be shown in this chapter that the proposed technique
provides parsimonious differential equations that yield fairly accurate predictive capacity, even
when the validation excitation is substantially different from the one used for identification. This
method overcomes the limitations associated with conventional techniques utilized for the identi-
fication of hysteretic systems by: (1) eliminating the need to restrict the representative model to a
specific form, i.e., Bouc-Wen model, and (2) providing an automatic scheme that uses discontin-
uous basis functions to discover the differential operator of systems that exhibit discontinuity in
their response.
5.1.4 Scope
A Genetic Programming-based technique is presented and applied for the identification of hys-
teretic systems, normally encountered in the Applied Mechanics field, to unveil the stable reduced-
order and reduced-complexity differential equation that governs the underlying hysteretic system.
This chapter is organized as follows:
Section 5.2 presents the hybrid evolutionary-based technique for the identification of hys-
teretic systems; Sections 5.3 and 5.4 present the result of the identification of hysteretic systems
governed by the Bouc-Wen model and the bilinear hysteretic model respectively; and Section 5.5
provides the conclusion.
85
5.2 Identification Scheme Using Genetic programming
5.2.1 General Formulation
In order to overcome the challenge of incorporating the response history in the mathematical
model of hysteretic systems, a secondary differential equation for the restoring force has to be
coupled with the primary ODE. The generic form of the governing ODE of hysteretic systems
being studied is:
f(t) =m x +r(x; _ x)
_ r(x; _ x) =g(x; _ x;r)
(5.1)
where f(t) is the applied force, r(x; _ x) is the restoring force, _ r(x; _ x) is the derivative of the
restoring force,m is the mass, andx, _ x and x denote the displacement, velocity and acceleration
of the system, respectively. The final goal is to discover the differential equation of the “restoring
force” to be eventually coupled with the primary ODE of the system, presented in Eq. (5.1).
Assuming that the external excitationf(t), the massm and the acceleration of the system x are
available from measurements, the restoring force is determined byr(x; _ x) = f(t)m x. The
time history ofx and _ x can be obtained through careful integration of x. Finite difference method
can be exploited to obtain an estimate of the derivative of the restoring force _ r(x; _ x) from the
restoring force data:
_ r(x
i
; _ x
i
) =
r(x
i+1
; _ x
i+1
)r(x
i
; _ x
i
)
t
) r(x
i+1
; _ x
i+1
) = t: _ r(x
i
; _ x
i
) +r(x
i
; _ x
i
) (5.2)
wherer(x
i+1
; _ x
i+1
) andr(x
i
; _ x
i
) are two consecutive datum in the restoring force data set, and
t is the time step. Equation (5.2) demonstrates how the restoring force at a specific instance in
time depends on the previous instance, and the way the response history is incorporated in the
modeling.
86
5.2.2 Investigated Models
The Bouc-Wen model is extensively employed in the Applied Mechanics field for characteriz-
ing a wide class of hysteretic systems where the shape of hysteretic loops is determined by the
model parameters. The bilinear hysteretic model, governed by algebraic equations, is another
widely used model for describing systems with yielding behavior. While parametric modeling is
a popular approach for identifying the model parameters or monitoring the changes, the proposed
non-parametric scheme in this study provides an equivalent system of differential equations to
represent the hysteretic systems governed by the Bouc-Wen or bilinear formulation.
5.2.3 External Excitation
The applied external excitations are stationary Gaussian white noise with zero-mean and a cer-
tain standard deviation that drive the hysteretic system into its nonlinear range. In addition, the
excitation should be persistent enough to provide suitable data for the identification.
5.2.4 Noise
In the real world, data that are collected from installed sensors on systems, and then transmitted to
processing units for analysis. Thus, data pollution is inevitable. Hence, in order to investigate the
situation in a more realistic manner, stationary Gaussian white noise with zero mean and known
standard deviation is added to all the “measured” signals.
5.2.5 Identification Procedure
The author has previously introduced a technique in Chapter 3 for the identification of SDOF
systems that exhibit nonlinearity in their response. This technique was also employed for the
modeling of nonlinear systems with discontinuity in their response in Chapter 4. That technique
87
is improved to be capable of embracing SDOF hysteretic systems investigated in this chapter. Fig.
(5.1) summarizes the main steps of the application of this method to hysteretic systems, and even-
tually the validation of the obtained models. Two essential phases are involved in the course of
the identification: first, training, where synthetically generated data are utilized, and the evolution
is directed to find the governing differential equation of the restoring force of the system. Second,
validation, where the performance of the discovered equation is evaluated through applying new
excitations to access its range of validity.
Evolution Initialization
The first step in the training phase is to generate the training excitation and solve the differen-
tial equation of the system by means of standard numerical approaches to obtain the reference
response. The displacementx, velocity _ x and acceleration x are polluted with 1% and 5% noise.
Moreover, the restoring forcer and its approximate derivative _ r are computed. These elements
form the training data set that will be used next during the evolution.
In the evolution process, the first generation is built by randomly chosen basis functions,
variables, algebraic operations as well as randomly generated numbers, and applying ramped
half-and-half method [13]. Each individual represents the tree form of the differential equation
of the restoring force and is meant to map the states and the restoring force to the derivative of
the restoring force. Hence, the displacement x, velocity _ x and the restoring force r(x; _ x), to-
gether with the available basis functions, form a three-dimensional function space search domain,
whose variables are represented asX1, X2 andX3, respectively. The fitness of every individual
is assessed based on the input, the desired output, and the model output. Thus, the termination
criteria is checked, and if a desirable solution is not found, evolutionary operations manipulate
the genetic materials to produce a new population. The evolution stops when the the fitness of
88
Figure 5.1: Flowchart of the proposed hybrid evolutionary technique for the identification of
hysteretic systems.
89
the best individual doesn’t improve more than 0:1% after 10
5
generations. The stop criteria is
determined based on different test cases to guarantee reaching the solution in a reasonable time.
Important factors involved in the evolution process of the tree structures are listed in Table 5.1.
GP Parameters Values GA Parameters Values
Population Size 100 Method gaoptimset
Mutation 0.20 Population Size 30
Crossover 0.80 Initial Population From GP:c
Representation 10% PopInitRange [0:05c 20c]
Penalty on ComplexityP
n
0.30 TolFun 1e 2
Basis Functions:step; abs; times; plus; minus
Frequency of Simplification and Parameter Optimization: 50 generations
Table 5.1: Properties of GP and GA for discovering the optimum structure and optimum parame-
ters, respectively.
Calculation of Fitness
Adaptation of an individual to the environment reflects its quality, and is determined on the ba-
sis of a fitness criterion. A suitable criterion is crucial for the success of the evolution. Thus,
it should be chosen very carefully and thoughtfully, with the consideration of the purpose of the
optimization. In this study, the eventual aim of the evolution is to generate parsimonious repre-
sentative models for hysteretic systems. Thus, the fitness is defined as the root-mean-square of
the difference between the derivative of the reference restoring force response _ r and the model
response
b
_ r, denoted ase = _ r
b
_ r, normalized by the root-mean-square (RMS) of the reference
90
response _ r, and a penalty on the complexityp
n
based on the number of nodes of the treen. A
reasonable penalty weightp
n
is specified from studying different test cases.
=
k _ r
^
_ rk
k _ rk
+p
n
:n =
kek
k _ rk
+p
n
:n (5.3)
Bloat Control: Linear Algebra
When tree structures grow enormously, bloat occurs which delays the termination of the evolution
process and causes the exceedance of available computer memory [34–36]. In addition, parsimo-
nious models are less likely to suffer from overfitting, and are more reliable for generalizing to
future dynamic environments. Therefore, upon the creation of new trees, the bloat control module
comes into play to simplify trees and reduce the complexity of the individuals. This essential task
is carried out by transforming tree structures to mathematical expressions, and simplifying them
by means of computer algebra through symbolic manipulations in Matlab. Removing redundant
parts, grouping similar terms, and condensing biases shrink the tree structure, and consequently,
reduce the complexity. Moreover, experience shows that rejecting trees whose depth is higher
than 13 results in robust results.
Parameter Optimization: Genetic Algorithms
During the course of evolution, the most appropriate randomly generated numbers survive along
with the most suitable tree structures. In addition to relying on GP to find the optimum numbers
through the evolution, a parameter optimization module is coupled with GP to optimize the param-
eters, and subsequently, to improve the fitness. Since discontinuous functions are present in the
equations, many parameter optimization techniques such as gradient based methods, Nedler-Mead
91
Simplex approach, and direct search method are not suitable. Hence, Genetic Algorithms are uti-
lized in Matlab to optimize the parameters after the equations are simplified symbolically. The
parameters of GAs used in this study for parameter optimization are listed in Table (5.1). The em-
bedded numbers in the equations are used to specify the initial population range of [0:05c : 20c],
where c is the vector of the current unoptimized parameters. The population size is 30 and TolFun
is 0:01. The GA factors are chosen in a way that parameters of the equations improve significantly,
while at the same time, the computational cost is affordable with respect to the entire evolution
speed. The parameter optimization terminates when the change in the best fitness function value
over 50 generations is less than or equal to TolFun (MATLAB default). Then, expression trees
are reconstructed from the simplified mathematical equations through a complex process that re-
quires extensive amount of programming. The fitness values are updated and a new generation
is formed. When the search termination criterion is met, the evolution stops and the discovered
differential equation of the restoring force is coupled with the primary ODE to form the system of
differential equations that characterize the behavior of the hysteretic system under investigation.
Finally, the solution of these differential equations is compared to the reference response to assess
the performance of models.
Validation
For validation, a different excitation is used for testing the obtained ODEs that represent the Bouc-
Wen model to assess their admissibility. In addition, three dissimilar excitations with various
stimulation levels, that are stronger and weaker than that of the identification, are considered
for the bilinear hysteretic system to assess the generalizability of the discovered solution in a
substantially different dynamical environment. The results of the application of the proposed
technique to the investigated models are presented next.
92
5.3 Identification of SDOF Hysteretic Systems Governed by Bouc-
Wen Model
The differential equation of a system with the hysteretic Bouc-Wen behavior is characterized by:
f(t) =m x +r(x; _ x) (5.4)
_ r(x; _ x) =
1

[A _ x(j _ xjjr(x; _ x)j
n1
r(x; _ x)
_ xjr(x; _ x)j
n
)] (5.5)
wheref(t) is the applied force,m is the mass,r(x; _ x) is the restoring force andx, _ x and x are
displacement, velocity and acceleration, respectively. The rest of the parameters, including,A,
, ,
, and n, are responsible for producing different characteristics, such as softening, hard-
ening, as well as hysteresis loops [101]. The parameters of the system and the properties of the
white Gaussian excitation, for training and validation, and the key parameters of the evolution are
listed in Table (5.2). In order to find the restoring force and the estimate of its derivative using
finite difference, the massm has to be known a priori. This requirement is not a major handicap
in realistic situations.
The fittest individual of the population, obtained after the termination of the optimization proce-
dure, is presented in Table (5.3) for different noise levels. Although the fitness of the equations
is not relatively high, it will be shown that the entire system of differential equations is a very
good representative of the hysteretic system. These coupled differential equations are solved by
conventional numerical approaches [116], and the solution is compared to the reference response
of the system in the training phase in Fig. (5.2). Despite the discontinuity of the response and the
presence of sharp corners, the provided models are able to capture the phenomena, and result in
acceptable estimates.
93
A new excitation, that has the same properties as the training excitation, is generated with a
different seed number to validate the results. The responses in the validation phase are compared
in Fig. (5.3), and the associated error is presented in Table (5.4). Figure (5.2) and Figure (5.3)
reveal that the saturation of the restoring force is captured effectively by the discovered differential
equations, and the estimates of the displacement and acceleration are reasonably accurate, even
when the oscillator is subjected to an excitation not encountered in the training phase. The phase
plots of the restoring force versus velocity and displacement under the validation excitation are
enlarged in Fig. (5.4).
Excellent match of the actual and estimated response for the hysteretic model implies that
coupling a secondary differential equation that represents the restoring force, with the primary
one, is the key to incorporate the response history in the representative model of the system, and
leads to the accurate identification of the memory-dependent phenomena. Results reveal that the
discovered ODE that represents the governing equation of the hysteretic system, approximates
the Bouc-Wen model with reasonable accuracy, yet with fewer parameters, and is less complex.
Hence, the introduced approach is able to identify the hysteretic systems characterized by Bouc-
Wen model fairly accurately and eliminates the restrictions associated with parametric modeling,
and the employment of a restricted class of mathematical models for various hysteretic phenom-
ena.
5.4 Identification of SDOF Hysteretic Systems Exhibiting Bilinear
Behavior
The bilinear hysteretic oscillator is governed by:
94
Excitation Parameters Values Bouc-Wen Values
Mean 0 m 2.00
Standard Deviation
t
and
v
1.00 1.00
Duration (Periods) 20T A 1.00
Sampling Freq. t=T 0.05 1.00
Basis Functions: 2.00
step; abs; times; plus; minus
-0.50
n 1.00
Table 5.2: Properties of the Bouc-Wen model, the applied excitations for both training and vali-
dation.
Type Governing Differential Equations RMS Error
f(t) = 2:00 x +r(x; _ x)
Actual _ r(x; _ x) = 1:00 _ x 2:00r(x; _ x)j _ xj + 0:50 _ xjr(x; _ x)j –
f(t) = 2:00 x +r(x; _ x)
0% noise GP: _ r(x; _ x) = 0:93 _ x 2:32r(x; _ x)j _ xj 10.3%
f(t) = 2:00 x +r(x; _ x)
1% noise GP: _ r(x; _ x) = 0:92 _ x 2:31r(x; _ x)j _ xj 10.46%
f(t) = 2:00 x +r(x; _ x)
5% noise GP: _ r(x; _ x) = 0:90 _ x 2:26r(x; _ x)j _ xj 15.45%
Table 5.3: The actual and discovered equations of motion of the Bouc-Wen oscillator and their
associated RMS error when polluted and unpolluted data are used.
f(t) =m x +c _ x +f
s
(x); f
s
(x) =
8
>
>
>
<
>
>
>
:
kx Linear region

2
kx Plastic region
(5.6)
95
States Training Data Validation Data
0% noise 1% noise 5% noise 0% noise 1% noise 5% noise
Displacementx 2.5% 2.5% 2.3% 1.9% 2.1% 3.2%
Velocity _ x 1.2% 1.3% 1.4% 1.2% 1.2% 1.6%
Acceleration x 0.9% 0.9% 1.1% 0.8% 0.9% 1.1%
Table 5.4: The normalized error associated with the estimated states in the training and validation
phase due to solving the GP-found ODEs corresponding to the Bouc-Wen oscillator.
wheref(t) is the external force,m is the mass,c is viscous damping,x, _ x and x are the displace-
ment, velocity and acceleration, respectively, andf
s
(x) is the stiffness force that is determined
based on the following two conditions:
(i). f
s
(x) = kx : Elastic range of 2f
y
or when linear motion recovers from the plastic region
once the direction of the motion alters.
(ii). f
s
(x) =
2
kx : Plastic range beyond the allowable linear range of 2f
y
.
The behavior of this system is shown schematically in Fig. (5.5). More details about other char-
acteristics of this phenomena can be found in [83, 117, 118]. Since no single-expression explicit
mathematical equation exists to fully characterize the bilinear hysteretic behavior, it is the ideal
case to challenge the capability of the proposed approach in identifying systems that exhibit such
behavior. The properties of the system, the applied excitations for both training and validation,
and the evolutionary parameters are presented in Table (5.6). The differential operator of the
restoring force of the bilinear system is obtained from GP, and is presented in Table (5.7) for
different pollution levels. Note that, since a two-point finite difference approximation is applied
to estimate the derivative of the restoring force, discovered ODEs have a significant fitness error
in Table (5.7). However, when they are embedded in the system of differential equations, the
96
0 20
−3.5
0
3.5
Excitation
t/T
(a) Training excitation
0 10 20
−50
0
10
Displacement
t/T
(b) Time history of displacement
0 10 20
−2
2
Acceleration
t/T
(c) Time history of acceleration
0 10 20
−0.5
0
0.5
Restoring Force
t/T
(d) Time history of restoring force
−50 0 10
−0.5
0
0.5
Restoring Force
Displacement
(e) Restoring force vs. Displacement
−3 0 3
−0.5
0
0.5
Restoring Force
Velocity

Exact
0% noise
1% noise
5% noise
(f) Restoring force vs. Velocity
Figure 5.2: Comparison of the reference and estimated response of Bouc-Wen oscillator to the
exciting force used for training. Note 4 superposed curves are shown in each panel.
estimated response is reasonably accurate. The reliability of these surrogate models are shown
by comparing their solution to the reference responses in Fig. (5.6). The corresponding error,
reported in Table (5.5), also confirms their acceptance.
In order to confirm the admissibility of the discovered models for the validation and gen-
eralization, three validation excitations, the properties of which are listed in Table (5.6), were
generated using different seed numbers and with different standard deviations. The intensity of
97
0 20
−3.5
0
3.5
Excitation
t/T
(a) Validation excitation
0 10 20
−30
0
5
Displacement
t/T
(b) Time history of displacement
0 10 20
−2
2
Acceleration
t/T
(c) Time history of acceleration
0 10 20
−0.5
0
0.5
Restoring Force
t/T
(d) Time history of restoring force
−30 0 5
−0.5
0
0.5
Restoring Force
Displacement
(e) Restoring force vs. Displacement
−3 0 3
−0.5
0
0.5
Restoring Force
Velocity

Exact
0% noise
1% noise
5% noise
(f) Restoring force vs. Velocity
Figure 5.3: Comparison of the reference and estimated response of Bouc-Wen oscillator to the
exciting force used for validation. Note that 4 superposed curves are shown in each panel.
the validation excitations is chosen once to be much higher, and once much lower than that of
the identification. Figures (5.7(a)-5.7(c)) show the intensity of the validation excitations. The
intensity of case 2 is the same as the identification excitation. The time history of displacement,
acceleration and the restoring force of the system for the range of 20T 40T , that is close to
the end of the response, are shown in Fig. (5.7(d)-5.7(l)). The associated error is also reported
in Table(5.5). The range of the plots are changed accordingly to provide better comparison plots.
98
Though the similarity of the intensity of the identification and validation excitation in case 2 re-
sults in the most accurate response estimation in Fig. (5.7(e)), Fig. (5.7(h)) and Fig. (5.7(k)),
the rest of the models provide acceptable estimates even when the intensity of the validation ex-
citation is extremely different. The plots of the restoring force vs. displacement in Fig. (5.7(m),
Fig. (5.7(n)) and Fig. (5.7(o)) show that the post-yield stiffness of the system is also tracked
effectively by the discovered differential equations. The only scenario in which the discovered
models are not reasonably accurate is when highly-polluted data are used for the training and the
validation excitation is three times more intense than the identification excitation. This inaccurate
estimate is reflected in the response of case 3 in the right column of Fig. (5.7). In other scenarios,
models yield reasonable estimates. The phase plots of the restoring force versus displacement
under the validation excitations are enlarged in Figs. (5.8), (5.9) and (5.10). Therefore, the dis-
covered ODEs are, in most cases, reliable representatives of the hysteretic system even when the
excitations are substantially different from the identification excitations.
States Training Data Validation Data (V1,V2,V3)
Noise 0% 1% 5% 0% 1% 5%
x 8.3% 9.5% 15.5% (3.8%,7.4%,3.8%) (5.2%,8.3%,5.2%) (42.7%,15.6%,73.3%)
_ x 6.9% 7.91% 10.2% (3.8%,6.5%,3.5%) (4.9%,7.3%,4.9%) (35.4%,12.0%,74.7%)
x 5.0% 5.6% 7.1% (2.47%,4.7%,3.5%) (3.4%,5.3%,3.4%) (21.9%,7.8%,51.4%)
r 8.0% 9.0% 11.5% (3.6%,7.4%,3.6%) (5.0%,8.3%,5.0%) (33.5%,12.2%,87.9%)
Table 5.5: The normalized error associated with the estimated states in the training and validation
phase due to solving the GP-found ODEs corresponding to the bilinear hysteretic oscillator.
99
Excitation Parameters Values Oscillator Values
Gaussian stationary Mean 0 m 2.00
white noise for Sampling Freq. t=T 0.05 !
n
0.40
all excitations Duration (Periods) 50T c 0.08
Training Standard Dev.
t
1.00 0.40
Validation 1 (V1) Standard Dev.
v1
0.33 x
Y
1.00
Validation 2 (V2) Standard Dev.
v2
1.00 k 0.32
Validation 3 (V3) Standard Dev.
v3
3.00 f
y
0.32
Table 5.6: Properties of the bilinear hysteretic oscillator and the applied excitations for both
training and validation.
Type Governing differential equations RMS error
f(t) = 2:00 x +r(x; _ x)
0% noise GP: _ r(x; _ x) = 0:31 _ x + 0:13xabs( _ x)
+0:08 _ xabs( _ x)r(x; _ x)abs( _ x) 28.5%
f(t) = 2:00 x +r(x; _ x)
1% noise GP: _ r(x; _ x) = 0:30 _ x + 0:11xabs( _ x) + 0:06 _ xabs( _ x)
0:05r(x; _ x) 0:83r(x; _ x)abs( _ x) 28.8%
f(t) = 2:00 x +r(x; _ x)
5% noise GP: _ r(x; _ x) = 0:30 _ x + 0:06xabs( _ x) + 0:04 _ x
2
r(x; _ x)
0:06r(x; _ x) 0:52r(x; _ x)abs( _ x) 33.0%
Table 5.7: The actual and discovered equations of motion of the bilinear oscillator and their
associated RMS error, when polluted and unpolluted data are used.
100
5.5 Conclusion
A non-parametric identification approach is introduced to discover a system of differential equa-
tions that characterize certain hysteretic systems, without the need to postulate any parametric
model for the system up-front. This approach employs Genetic Programming as the main opti-
mization engine to evolve a population of differential equations stored in expression trees, and
incorporates it with parameter optimization methods using Genetic Algorithms, as well as linear
algebra techniques to improve the fitness and to simplify redundancies in trees, respectively. The
fitness of each tree is determined based on the RMS error between the derivative of the restoring
force, obtained using finite difference methods, and the model estimates. In addition, a penalty
on the complexity is considered to favor parsimonious solutions.
It is shown that the proposed hybrid technique results in stable and parsimonious models in
the form of differential equations with reasonably accurate response estimates for the Bouc-Wen
model where a differential equation governs the evaluation of the hysteretic force. This technique
also yields accurate models for the bilinear hysteretic oscillator, which is defined through a set
of algebraic expressions that are conditioned on the history of the dynamic system response.
The performance of the models are also tested under validation excitations that are chosen to
be, in one case more intense, and in another case less intense than the identification excitation.
Good matching between the reference and the estimated response confirms the validity of the
provided models. Thus, the presented non-parametric technique provides an equivalent system of
parsimonious differential equations to characterize the behavior of hysteretic systems effectively.
101
−30 0 5
−0.5
0
0.5
Restoring Force
Displacement
(a) Restoring force vs. Displacement
−3 0 3
−0.5
0
0.5
Restoring Force
Velocity

Exact
0% noise
1% noise
5% noise
(b) Restoring force vs. Velocity
Figure 5.4: Comparison of the phase plots of the reference and estimated restoring force of Bouc-
Wen oscillator to the exciting force used for validation. Note that 4 superposed curves are shown
in each panel.
102
Figure 5.5: Behavior of the stiffness force in the bilinear hysteretic oscillator.
0 50
−3.5
0
3.5
Excitation
t/T
(a) Training excitation,t= 1.00
20 40
−13
0
13
Displacement
t/T
(b) Displacement
20 40
−2
2
Acceleration
t/T
(c) Acceleration
20 40
−2
0
2
Restoring Force
t/T
(d) Restoring force
−13 0 13
−2
0
2
Restoring Force
Displacement
(e) Restoring force vs. Displacement
Figure 5.6: Comparison of the reference and estimated response of the bilinear hysteretic oscil-
lator to the exciting force used for training. Note that 4 superposed curves are shown in each
panel.
103
0 50
−10
0
10
Excitation
t/T
(a) V1 Excitation,v1= 0.33
0 50
−10
0
10
Excitation
t/T
(b) V2 Excitation,v2= 1.00
0 50
−10
0
10
Excitation
t/T
(c) V3 Excitation,v3= 3.00
20 40
−3.5
0
3.5
Displacement
t/T
(d) V1 Displacement
20 40
−13
0
13
Displacement
t/T
(e) V2 Displacement
20 40
−50
0
50
Displacement
t/T
(f) V3 Displacement
20 40
−0.75
0.75
Acceleration
t/T
(g) V1 Acceleration
20 40
−2
2
Acceleration
t/T
(h) V2 Acceleration
20 40
−7
7
Acceleration
t/T
(i) V3 Acceleration
20 40
−0.7
0
0.7
Restoring Force
t/T
(j) V1 Restoring force
20 40
−2
0
2
Restoring Force
t/T
(k) V2 Restoring force
20 40
−10
0
10
Restoring Force
t/T
(l) V3 Restoring force
−3.5 0 3.5
−0.7
0
0.7
Restoring Force
Displacement

Exact
0% noise
1% noise
5% noise
(m) V1: Restoring force vs. Dis-
placement
−13 0 13
−2
0
2
Restoring Force
Displacement

Exact
0% noise
1% noise
5% noise
(n) V2: Restoring force vs. Dis-
placement
−50 0 50
−10
0
10
Restoring Force
Displacement

Exact
0% noise
1% noise
5% noise
(o) V3: Restoring force vs. Dis-
placement
Figure 5.7: Comparison of the reference and estimated response of the bilinear hysteretic oscilla-
tor subjected to the exciting forces used for validation. Note that 4 superposed curves are shown
in each panel.
104
−3.5 0 3.5
−0.7
0
0.7
Restoring Force
Displacement

Exact
0% noise
1% noise
5% noise
Figure 5.8: Comparison of the phase plots of the reference and estimated restoring force of the
bilinear hysteretic oscillator subjected to the exciting force used for validation with
v1
= 0:33.
Note that 4 curves are superposed for different noises.
105
−13 0 13
−2
0
2
Restoring Force
Displacement

Exact
0% noise
1% noise
5% noise
Figure 5.9: Comparison of the phase plots of the reference and estimated restoring force of the
bilinear hysteretic oscillator subjected to the exciting force used for validation with
v2
= 1:00.
Note that 4 curves are superposed for different noises.
106
−30 0 30
−4
0
4
Restoring Force
Displacement

Exact
0% noise
1% noise
5% noise
Figure 5.10: Comparison of the phase plots of the reference and estimated restoring force of the
bilinear hysteretic oscillator subjected to the exciting force used for validation with
v3
= 3:00.
Note that 4 curves are superposed for different noises.
107
Chapter 6
Modeling of Human Spine
6.1 Introduction
T
HE evolution of living beings for millions of years resulted in the creation of humans,
which is believed to be the most intelligent living being on Earth. One of the challenges
of our time is to better understand the humans body to discover effective medicine to cure ill-
nesses, to invent new devices to improve the body’s health or replace malfunctioning organs, and
to design new exercises to enhance the body’s physical abilities. Our brain, as the most impor-
tant part of our body, plays the major role in controlling the behavior of different organs. The
human spine is a very important part of our skeleton on which our movement and stability highly
depends. In order to control the movement of the spine, brain orders are sent to the vertebrae
as electric signals through nerve channels of the nervous system. However, the vertebral column
may not be capable of obeying the brain orders due to birth defects or other issues. The sources
of spinal issues can be infections, injuries, tumors, bone changes due to aging and other special
conditions.
108
Disc disease is one of the most prevailing spinal problems that may lead to surgery, among
which disc degeneration is the one that can significantly influences one’s life. Various devices are
invented to treat the disc degeneration disease. Many posterior dynamic stabilization devices have
been utilized for the treatment of degenerative disc disease. X-Stop made by St Francis Medical
Technologies in Alameda, California, is extensively investigated for the treatment of degenerative
disc disease as an interspinous-process spacer [119–127]. Sangiorgio et al. [2] investigated the
range of motion of a spinal motion segment intact, injured, and treated using three different pos-
terior dynamic stabilization devices, that are the PrecuDyn, the X-stop and the Isobar posterior
stabilization. They showed that each implant has some advantages and improves the range of mo-
tion in a certain way. One approach to investigate the spinal motion under different conditions, is
to have a representative model for a healthy spine, and injured one and with implants to compare
their range of motion under a similar external stimulation. This model can also be used for the
prediction of the spine response subjected to a not-previously experienced stimulation. Hence,
the modeling of the spine can be used for the diagnosis of spine problems, the design of spine
related devices, and inventing new exercises for strengthening the spine muscles.
The evolutionary identification technique introduced in Chapter 5 is utilized in this chapter
for obtaining the differential equations that govern the behavior of a spinal motion segment under
three conditions: intact, injured and with X-Stop implant. The employed data set is obtained
from the experiments conducted by Sangiorgio et al. [2] on spinal motion segment under three
aforementioned conditions. More details about the experiment and the employed data set as well
as the method of analysis will be presented in the next section.
109
6.1.1 Review of Previous Work
Much research has been conducted to model and understand the biomechanical behavior of the
human spine. Notable investigations in the literature include considering a continuum model to
model the spine [128], studying the kinematics of the human spine [129, 130], providing a dis-
crete three-dimensional mathematical model for the spine [131, 132], modeling the muscle action
and the stability of the spine [133], introducing a geometric model [134] and using a viscoelastic
model to investigate the spine [135]. More studies on the mathematical modeling of the human
spine can be found in [136]. With the advances in processing capabilities of computers, more
sophisticated models can be used to represent the spine, and predict its biomechanical behavior.
Finite element methods are a class of approaches that are widely utilized for the modeling of hu-
man spine and its interaction with surrounding organs. Significant research has been conducted
in this area such as the work of [137–143]. However, modeling the spine using finite element
methods requires the consideration of many parameters associated with the spine and the adjacent
interacting body elements. Thus, due to the complexity of the inner body interactions, these mod-
els may not be able to predict the behavior of spine accurately, nor can be easily implemented in
biomedical devices for testing and diagnosis.
The goal of this study is to use the data obtained from real test cases to provided a fairly
simple mathematical models in the form of differential equations to characterize the behavior of
the human spine. Such a model can be used to examine the health of the spine through monitoring
its movement, or to detect spinal injuries and diseases. Simplicity of these models can also lead
to a better understanding of the problem, and provides insight to the phenomena.
110
6.1.2 Scope
The evolutionary based identification scheme that was introduced in chapter 5 is employed to
discover the differential equations that characterize the motion of spine due to certain excitations.
The materials in this chapter are organized in the following manner: Section 6.2 describes the
method of experimentation and the data set used for the analysis; Section 6.3 presents the general
identification procedure; Section 6.4 provides the identification results that include the discovered
models, and compare their estimates to the actual responses; Section 6.5 discusses the findings,
and Section 6.6 provides the summary and the conclusion of this chapter.
6.2 The Method of Experimentation
6.2.1 Experimental Materials
According to [2], nine lumbar spines were prepared to the experiments. The excessive fat tissues
and paravertebral muscles were removed, and the discs, facet joints, spinous processes, as well as
transverse processes were preserved. In addition, ligaments of transverse and spinous processes
were left intact. After conducting the test on intact specimens, the injured specimen was created
by making radial tears and rim lesions [2]. It was followed by the experimentation of X-stop
implants. The schematic posterior and lateral view of a single motion segment with and without
the X-stop implant is shown in Fig. (6.2). More detailed information on the testing procedure can
be found in the work of Sangiorgio et al. [2].
6.2.2 Experimental Setup
All experiments were conducted by an MTS 858 mini-bionix servohydrolic load frame that in-
corporated the Flextest system. The motion was recorded using Optotrak 3020 3-dimensional
111
motion tracking system [1]. A 6 degree-of-freedom load cell is also mounted on the apparatus to
record forces and moments continuously throughout the experiment. The test setup is shown in
Fig. (6.1) [1].
Figure 6.1: The test setup [1].
6.2.3 Data Set
The data set available for this investigation included the applied torsional torquet
r
to the spine, the
measured axial forcef
d
, the measured rotational torquer, and the measured axial displacement
112
Figure 6.2: The motion segment is shown schematically. The top picture is intact, and the bottom
picture is with X-stop implant [2].
d. The time history of all quantities that govern the kinematics and their normalized values are
plotted in Fig. (6.3) and Fig. (6.4), respectively. Due to the complexity of the problem and the
dependency of the movement in one direction to the other directions and the parameters of the
system, the phase plots of the actual and normalized data are plotted in Fig. (6.5) and Fig. (6.6),
respectively. They provide insight about the physics behind the problem, and help us with the
modeling. Data normalization is carried out by mapping all data to the range of [-1 1]. The
normalized vector s
n
can be calculated by shifting the average of data s to zero, and scaling them
properly, using the following equation:
s
n
=
s
smax+s
min
2
smaxs
min
2
(6.1)
The provided data are the response of a spinal motion segment under different conditions
subjected to identical torsional torques. These conditions are intact, injured, and with X-stop
implant. Most of the phase plots in Fig. (6.5) do not reveal the physics behind the problem.
113
0 50 100 150
0.5
1
1.5
2
Axial Displacement (mm)
Time (sec)
0 50 100 150
−5
0
5
Axial Force (N)
Time (sec)
0 50 100 150
−5
0
5
Torsional Rotation (Degree)
Time (sec)
0 50 100 150
−1
−0.5
0
0.5
1
x 10
4
Torsional Torque (N−mm)
Time (sec)

Intact
Injured
X−Stop
Figure 6.3: The time history of the data used for the modeling of spine.
However, the phase plot of the torsional torquet
r
vs torsional rotationr, that is repeated in Fig.
(6.7), is very similar to the hysteretic behavior of the mechanical systems that are investigated
in Chapter 5. Thus, the eventual aim is to find a differential equation that defines the derivative
of the torsional rotation, denoted as angular velocity _ r, in terms of the applied loads and the
measured states. The general formulation of the problem and the proposed identification scheme
is presented next.
114
0 50 100 150
−1
0
1
Normalized Axial Displacement
Time (sec)
0 50 100 150
−1
0
1
Normalized Axial Force
Time (sec)
0 50 100 150
−1
0
1
Normalized Torsional Rotation
Time (sec)
0 50 100 150
−1
0
1
Normalized Torsional Torque
Time (sec)

Intact
Injured
X−Stop
Figure 6.4: The time history of the normalized data used for the modeling of spine.
6.3 Non-Parametric Modeling of Human Spine Using Genetic pro-
gramming
6.3.1 Modeling Procedure
The identification technique discussed in Chapter 5 is employed to discover a model in the form
of differential equations that govern the behavior of the spine under three conditions of intact,
injured and with X-stop implant. The general form of the governing differential equations is
defined as follows:
_ r =f(t
r
;f
d
;r;d) (6.2)
115
1 1.5 2
−10
−5
0
5
10
Axial Force (N)
Axial Displacement (mm)
1 1.5 2
−10
−5
0
5
10
Torsional Rotation (Degree)
Axial Displacement (mm)
1 1.5 2
−1
−0.5
0
0.5
1
x 10
4
Torsional Torque (N−mm)
Axial Displacement (mm)
−5 0 5
−10
−5
0
5
10
Axial Force (N)
Torsional Rotation (Degree)
−5 0 5
−1
−0.5
0
0.5
1
x 10
4
Axial Force (N)
Torsional Torque (N−mm)
−5 0 5
−1
−0.5
0
0.5
1
x 10
4
Torsional Torque (N−mm)
Torsional Rotation (Degree)

Intact
Injured
X−Stop
Figure 6.5: The phase plots of the data used for the modeling of spine.
where the angular velocity _ r is defined in terms of the torsional torquet
r
, the axial forcef
d
, the
torsional rotationr and the axial displacementd. Since the data are not in a comparable range,
selecting variables and providing equations that fit the data reasonably accurately may not be fea-
sible. Thus, the data are normalized to [-1 1] using the equation introduced earlier in this chapter.
116
−1 0 1
−1
0
1
Normalized Axial Force
Normalized Axial Displacement
−1 0 1
−1
0
1
Normalized Torsional Rotation
Normalized Axial Displacement
−1 0 1
−1
0
1
Normalized Torsional Torque
Normalized Axial Displacement
−1 0 1
−1
0
1
Normalized Axial Force
Normalized Torsional Rotation
−1 0 1
−1
0
1
Normalized Axial Force
Normalized Torsional Torque
−1 0 1
−1
0
1
Normalized Torsional Torque
Normalized Torsional Rotation

Intact
Injured
X−Stop
Figure 6.6: The phase plots of the normalized data used for the modeling of spine.
The finite difference method was employed to obtain the estimate of the angular velocity _ r
from torsional rotationr measurements. Hence, the normalized recorded torsional rotationr and
the estimate of its derivative, the normalized torsional torquet
r
, the normalized axial forcef
d
and
the normalized axial displacementd construct the training data set provided to the identification
scheme described in Chapter 5. The properties of the excitation and the evolutionary parameters
117
−5 0 5
−1
−0.5
0
0.5
1
x 10
4
Torsional Torque (N−mm)
Torsional Rotation (Degree)

Intact
Injured
X−Stop
Figure 6.7: The phase plots of torsional torque vs. torsional rotation.
are presented in Table (6.1). Thus, the proposed hybrid technique is utilized to discover repre-
sentative models of the spinal motion segment under intact condition. Then, the coefficients of
this differential equation are optimized for the injured and X-stop scenarios.fmincon command
usedinterior-point algorithm for the parameter optimization. Note thatfmincon is not
the best option for solving highly nonlinear and complex optimization problems such as the the
problem under discussion, but the extreme computational cost of Genetic Algorithms and other
evolutionary based methods, drove us to use this method. The coefficients given by GP are con-
sidered to be the starting point. However, the range of variation during the optimization, should
be restricted to avoid the differential equations blow up. The parameter optimization factor are
also listed in Table (6.1). The models and the response estimate are presented in the next section.
118
Excitation Parameters Values GP Parameters Values Parameters Values
Peak of Sinusoidal 80000 Population Size 100 Method fmincon
Torsional Torquet
r
Nmm Crossover Rate 0.80 Parameters Variability 10%
Number of Cycles 5 Mutation Rate 0.20 TolFun 1e 10
Durationt 90 sec Representation 10% MaxFunEvals 10e10
Basis Functions:step; sgn; abs; times; plus; minus MaxIter 10e10
Frequency of Simplification and Parameter Optimization: 50 gen. TolX 1e 15
Table 6.1: Properties of evolutionary based identification approach for discovering the differential
equation that governs the motion of the spine.
6.4 Identification Results
The proposed GP-based identification approach yields differential equations that characterize the
behavior of the spinal motion segment described earlier. The discovered differential equations in
the normalized domain are:
(1) _ r =a
1
t
r
+rabs(a
2
r)a
3
ra
4
rabs(t
r
)
(2) _ r =a
1
+a
2
t
r
+a
3
t
2
r
r +a
4
t
3
r
a
5
ra
6
rabs(t
r
)a
7
t
r
abs(t
r
)a
8
t
2
r
rabs(t
r
)
(6.3)
where _ r is the angular velocity,r is the torsional rotation,t
r
is the torsional torque anda
1
a
8
are the numeric coefficients of the ODEs. It is worth mentioning that the data corresponding tof
d
andd were also included in the training data set for consideration but the applied technique didn’t
find the participation of these variables significant to include them in the governing differential
equations.
After the termination of the evolution, the coefficients of the obtained differential equations are
119
Type Condition Coefficients RMS error
a
1
a
2
a
3
a
4
a
5
a
6
a
7
a
8
k^ r rk
ODE 1
Intact 1.878 1.522 1.847 1.718 - - - - 12.37 %
Injured 2.600 2.034 2.138 2.647 - - - - 10.96%
X-stop 3.152 1.570 1.794 3.179 - - - - 11.04 %
ODE 2
Intact 0.001 3.891 4.372 4.433 1.503 0.872 7.104 1.71 12.76 %
Injured 0.0344 5.292 5.142 5.278 1.504 1.584 9.157 3.831 10.87 %
X-stop 0.001 5.334 3.147 4.602 1.267 1.751 7.897 2.618 10.36 %
Table 6.2: The optimized coefficients of the ODEs and the RMS error associated with the rotation
in the denormalized domain.
optimized using genetic algorithms to provide the best fit to the reference torsional rotationr for
intact, injured and X-stop scenarios. Hence, for the parameter optimization process,k^ r rk is the
cost function, the population size is 80 and the termination point is when the change in the best
fitness value, over 50 consecutive generations, doesn’t exceedTolFun = 0:01. These values are
obtained by experimental studies to ensure accurate results. Note that the parameter optimization
process, after the termination of the main course of search for ODEs, is computationally very
expensive. Because for each set of parameters in the population, the associated differential equa-
tions should be solved numerically to assess the corresponding fitness.
The coefficients of these equations and the corresponding error associated with them are presented
in Table (6.4). The bar graphs of the coefficients are also compared in Fig. (6.8).
The differential equations with optimized parameters are solved using standard numerical
methods, and the solution is transformed from the normalized domain to the denormalized do-
main. The results are plotted over the reference response in Fig. (6.9), Fig. (6.10) and Fig. (6.11)
120
a1 a2 a3 a4 a5 a6 a7 a8
0
5
10
ODE 1 Coefficients
Values
Coefficients
a1 a2 a3 a4 a5 a6 a7 a8
0
5
10
ODE 2 Coefficients
Coefficients
(a) The coefficients of ODE 1 in Equ. (6.3). (b) The coefficients of ODE 2 in Equ. (6.3).
Figure 6.8: (a) The bar graph of the coefficients of the differential equations for spinal motion
segment intact (blue), injured (green) and X-stop (red)).
to assess the accuracy of the representative models in intact, injured and X-stop scenarios, respec-
tively. Transforming the data from the normalized domain s
n
to the denormalized domain s is
carried out by reversing the normalization process as follows:
s = s
n
s
max
s
min
2
+
s
max
+ s
min
2
(6.4)
The results reveal that the discovered models can fairly accurately predict the response of
the spinal motion segment. The error is within a reasonable range, and the hysteretic loops are
captured effectively. Important findings are discussed more in the discussion section coming next.
6.5 Discussion
Due to the employment of an evolving population to obtain differential equations that character-
ize the systems based on the experimental data, the candidate equations are not unique. In the
previous section, the results of two candidate solutions, one of which was more complex, were
presented. The results show that the discovered ODEs can fairly accurately represent the system
121
0 50 100
−5
0
5
Actual Torsional Rotation
Time
−5 0 5
−1
−0.5
0
0.5
1
x 10
4
Actual Torsional Torque
Actual Torsional Rotation

Exact
ODE
0 50 100
−5
0
5
Actual Torsional Rotation
Time
−5 0 5
−1
−0.5
0
0.5
1
x 10
4
Actual Torsional Torque
Actual Torsional Rotation

Exact
ODE
(a) Short ODE 1 for intact sample (b) Long ODE 1 for intact sample .
Figure 6.9: The solution of the short and long ODE 1 against the reference response for the intact
spine segment motion.
under the investigation. The hysteresis loops were successfully identified. However, the slight
degrading behavior could not be captured effectively using a single differential equation.
Normalizing the variables and loads was the key to the success of the approach. Otherwise
the extreme difference in the range of motion and loads could prevent the identification technique
to give robust models. It was also observed that the more complex differential equation, denoted
as ODE 2 in Eq. (6.3) does not yield significantly improved estimates. Therefore, using a more
complex model is not necessary.
Another goal was to carry out statistical analysis to infer about the parameters of the differential
equations and their dependence on the spine condition through analyzing bar graphs. However,
122
0 50 100
−10
−5
0
5
10
Actual Torsional Rotation
Time
−10 −5 0 5 10
−1
−0.5
0
0.5
1
x 10
4
Actual Torsional Torque
Actual Torsional Rotation

Exact
ODE
0 50 100
−10
−5
0
5
10
Actual Torsional Rotation
Time
−10 −5 0 5 10
−1
−0.5
0
0.5
1
x 10
4
Actual Torsional Torque
Actual Torsional Rotation

Exact
ODE
(a) Short ODE 1 for injured sample (b) Long ODE 1 for injured sample
Figure 6.10: The solution of the short and long ODE 1 against the reference response for the
injured spine segment motion.
due to a lack of sufficient number of samples, we were not able to carry out an statistical analysis
on the parameters of the differential equations.
6.6 Summary and Conclusion
The GP-based hybrid identification method presented in Chapter 5 was used for the identifica-
tion of spinal motion segment intact, injured and with X-stop implant. The experimental data
obtained from different spinal conditions were provided by [2]. These data were utilized to dis-
cover a representative model in the form of differential equations. The eventual goal was to find
an expression that relates the angular velocity to the measured states and loads, involved in the
123
0 50 100
−10
−5
0
5
10
Actual Torsional Rotation
Time
−10 −5 0 5 10
−1
−0.5
0
0.5
1
x 10
4
Actual Torsional Torque
Actual Torsional Rotation

Exact
ODE
0 50 100
−10
−5
0
5
10
Actual Torsional Rotation
Time
−10 −5 0 5 10
−1
−0.5
0
0.5
1
x 10
4
Actual Torsional Torque
Actual Torsional Rotation

Exact
ODE
(a) Short ODE 1 for X-stop sample (b) Long ODE 1 for X-stop sample .
Figure 6.11: The solution of the short and long ODE 1 against the reference response for the spine
segment motion with X-stop.
kinematic of the spinal motion segment. It was shown that the proposed hybrid evolutionary iden-
tification technique discovers differential equations whose solution match the experimental data
fairly accurately. Two candidate solutions, with different number of terms, that provide the best
fit were presented, and their coefficients were optimized for various spinal conditions. The results
also reveal that more complex models do not necessarily increase the accuracy of the estimates.
Thus, parsimonious models have superiority over more complex models.
124
Chapter 7
Identification of MDOF Systems with Hysteretic Behavior
7.1 Introduction
7.1.1 Background
M
ODELING and analysis of complex nonlinear hysteretic systems is a broad research
area with applications in civil, aerospace and mechanical systems. The behavior of
structures in severe earthquakes, the response of aerospace devices involving joints, and the mo-
tion of bearings in machine tools are examples of such phenomena. Despite outstanding ad-
vancements in the processing capabilities of modern computers, there is a scarcity of suitable
approaches that benefit from high processing capacities for the identification of complex non-
linear multi-degree-of-freedom (MDOF) systems with hysteretic behavior, typically encountered
in Applied Mechanics field. For example, parametric identification is a class of techniques that
requires prior knowledge about underlying systems up-front, and excessive user dependency, to
provide a suitable model based on the physics of the problem. Data driven modeling through
black-box simulations is another class of techniques that carries out the identification blindly, and
lacks disclosing the physics of the studied phenomena [56]. The author has developed and im-
plemented a promising technique that utilizes computational intelligence approaches to eliminate
125
the aforementioned restrictions contained in conventional techniques to identify various single-
degree-of-freedom complex non-linear systems. That technique is improved and presented in
this chapter to be utilized for the identification of non-conservative dissipative multi-degree-of-
freedom systems.
7.1.2 Review of Related Studies
Over the past several decades, nonlinear systems incorporating non-conservative dissipative phe-
nomena have been investigated in many studies, mainly through parametric or nonparametric
modeling [56]. Parametric modeling is a class of techniques used to characterize hysteretic sys-
tems by attributing a suitable model to the system, and identifying the model parameters or detect-
ing their changes [63, 80, 88–90, 92, 99, 100, 110, 144–152]. However, the success of parametric
modeling highly relies on the user providing an accurate representative model based on a clear
understanding of the physical phenomena. Non-parametric modeling, on the other hand, is an-
other class of approaches which carries out the identification automatically, without setting any
major hypothesis about the system upfront. In this category, Neural Networks are powerful tools
of modeling complex systems involving memory-dependent characteristics [101, 111, 153–162].
Their data driven analysis is suitable when no additional information about the underlying sys-
tems is available, and they result in fairly accurate estimates. However, Neural Networks process
data in a completely blind manner, and result in black-box models, that are not suitable for reveal-
ing the physics behind complex systems.
Polynomial-based approaches are another class of nonparametric techniques that are extensively
employed for the identification of memory-dependent phenomena in dynamic environment [40,
48, 68, 163–167]. However, specifying the order of polynomials and the composition of basis
functions require insight into the dynamics of the problem.
126
Currently, there is a lack of automatic computational techniques that provide robust physics-based
computational models, which reflect hysteretic behavior of dynamical MDOF systems, exists in
the field of Applied Mechanics.
7.1.3 Advances in Computational Intelligence
Computational intelligence methodologies are inspired by natural phenomena to provide elegant
solutions to complex real-world problems in various fields. In this class of problem-solving tech-
niques, Genetic Programming (GP), is built on evolutionary algorithms, and offers a great poten-
tial for the identification of complex dynamical systems exhibiting non-conservative dissipative
behavior. GP maintains a population of computer programs, and evolves them by means of bio-
logical crossover and mutation operations. The performance of each individual in the domain is
evaluated using a fitness criteria. GP has proven to be successfully applicable to various classes of
problems [14–24, 168, 169]. To our knowledge, GP has never been adapted to be employed for the
identification of MDOF systems with hysteretic behavior. Therefore, this investigation explores
the potential of GP combined with nonlinear parameter optimizersto discover the mathematical
representation of nonlinear hysteretic MDOF systems.
7.1.4 Scope
A general procedure concerning the identification of complex single-degree-of-freedom systems
associated with challenging type of nonlinearities was presented in previous chapters. That proce-
dure is extended in this chapter and is implemented for the modeling and analyzing of non-linear
MDOF systems with memory-dependent dissipative characteristics. The eventual aim is to utilize
GP as the main problem-solving engine for discovering the suitable structure of the differential
127
equations that govern the behavior of the system of interest. GP is also coupled with nonlin-
ear parameter optimization techniques employing Genetic Algorithms to optimize the parameters
embedded in the equations, and consequently to improve the estimates.
It is shown that the suggested methodology provides reduced-complexity parsimonious systems
of differential equations that govern the dynamics of complex nonlinear MDOF systems with
non-conservative dissipative traits. Several test cases are provided to assess the applicability, reli-
ability, and validity of the methodology.
This chapter is organized as follows: an introduction to GP within the scope of the investiga-
tion in this chapter is presented in Section 7.2; assumptions, the formulation and the identification
procedure is presented in Section 7.3; the results of the identification of MDOF hysteretic systems
is presented in Section 7.4; statistical analysis of validation results and in-depth interpretation os
findings are discussed in Section 7.5, and the conclusion is provided in Section 7.6.
7.2 Introduction to Genetic Programming
7.2.1 Overview
Evolutionary Algorithms (EAs) are heuristic optimization methods premised on the theory of
natural selection. The fundamental idea behind EAs is to mimic the necessary processes for
evolution in nature —selection, mutation and crossover— to evolve a population of candidate
solutions. Defining a suitable fitness measure is also a critical aspect of the optimization process
to direct the evolution toward an acceptable result in a timely manner. Genetic Programming
(GP) is a branch of EAs, introduced by John Koza [13] which employs the same fundamental
128
principals. GP evolves any type of structure, presented in expression tress, whose performances
is quantifiable in the domain.
7.2.2 Population
In this study the population of evolving candidate solutions are the differential equations of the
restoring force at each degree-of-freedom (DOF), that will eventually take part in the system of
the differential equations that fully characterize the behavior of the MDOF system. In order to
obtain distinct and independent differential equations, a new population is formed for analyzing
each degree of freedom. No additional information is assumed to be known about the topol-
ogy of the system. The training data set is composed of all displacement and velocity states,
x
1
;x
2
;x
3
; _ x
1
; _ x
2
; _ x
3
as well as the compound restoring forcer
i
at a specific DOFi. Hence, they
form a 7 dimensional function space domain to estimate the derivative of the restoring force _ r
i
at that specific DOFi. In addition to the state variables, algebraic operations (+;;), as well
asabs,sgn and Heavisidestep functions are included in the library of essential building blocks
to construct the body of the evolving population. However, only appropriate elements will sur-
vive during the course of the evolution to form the most suitable structures, and subsequently,
the most accurate differential equations. The ramped half and half technique [13] is utilized to
construct the first population. Then evolutionary operators, described next, are used to advance
the evolution. The eventual goal is to discover independent and distinct differential equations for
the compound restoring forces of the system to form a complete system of differential operators
that fully characterize the system under discussion.
129
7.2.3 Evolutionary Operators
In order to imitate biological evolution, mutation and crossover operators are applied to the pop-
ulation. Mimicking mating in nature, crossover is an essential operation to direct the evolution.
Crossover operator interchanges two randomly selected limbs of parent trees to produce two off-
springs. Mutation operator replaces a randomly selected branch of the target tree by a new tree
grown at random. Similar to natural phenomena embracing evolution, mutation is a background
operation which is applied at a very low rate. It is primarily intended to generate small perturba-
tions to escape the optimization from local exterema. However, in GP this operation may account
for a significant transformation of trees. Mutation, in theory, provides a chance for GP to explore
the entire domain of the potential solutions. The adaptation of crossover and mutation operations
for evolving expression trees of differential equations is elaborated more extensively in Chapter
2.
Beside these two primary evolutionary operations, highly-fit individuals in the population are
preserved against not always beneficial evolutionary operators. Thus, a percentage of the fittest
individuals are transferred to the next generation intact. Selecting individuals is also required for
the implementation of evolutionary operators. The effectiveness of different selection schemes
is broadly studied in the literature. Many selection techniques, such as Truncation scheme [30],
roulette-wheel approach [31] and Stochastic Universal Sampling method [32], that do not give
a chance to not highly fit individuals to survive. However, Tournament selection [33] provides
an opportunity for less-fit individuals to be considered for future generations because they may
contain valuable genetic materials. Therefore, Tournament selection scheme is implemented in
this investigation.
130
7.2.4 Fitness Criterion
Establishing a suitable fitness criterion is critical for guiding the evolution toward an admissi-
ble solution in a timely manner. Inverse problem investigations, concerned in this investigation,
ideally result in discovering a representative model for the unknown system. Therefore, the sim-
ilarity of the response of the unknown system x and the response of the candidate model ^ x can
be a robust measure of the candidate model fitenss. Hence, the fitness error is calculated using
the root-mean-square (RMS) of the deviation array e = x ^ x, normalized by the RMS of the
reference response x:
=
kx ^ xk
kxk
(7.1)
7.3 Identification of MDOF Systems: General Procedure
7.3.1 Assumptions
Investigated Models
The identification method under discussion is implemented on a benchmark MDOF system ex-
hibiting non-conservative dissipative behavior. This structure has been employed for testing the
robustness of other identification schemes in the past [101, 170, 171]. This three degree-of-
freedom system in Fig. (7.1) is composed of three unequal lumped masses that are linked by six
arbitrary components to all other masses and to a fixed support.
131
Figure 7.1: The benchmark 3 degree-of-freedom system with interconnecting links and supports.
In Fig. (7.1), x
i
is the displacement, f
i
(t) is external excitation, and m
i
is the mass at
DOF i. The semi-physical Bouc-Wen model is well studied in the literature to represent non-
conservative dissipative systems. The interconnecting elements in the structure above, denoted as
g
i
, are governed by the Bouc-Wen model, and consequently, account for the hysteretic behavior
of the MDOF system. The Bouc-Wen model for a single degree-of-freedom system is formulated
as follows:
f(t) =m x +r(x; _ x) (7.2)
_ r(x; _ x) =
1

[A _ x(j _ xjjr(x; _ x)j
n1
r(x; _ x)
_ xjr(x; _ x)j
n
)] (7.3)
wheref(t) is the external excitation,m is the lumped mass, x is the acceleration,r(x; _ x) is the
restoring force,x is the displacement, and _ x is the velocity. The shape of the hysteretic loops is
132
ruled by the parameters of the model: ,A,,,
, andn [101]. Therefore, the MDOF system
under discussion is formulated as follows:
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
f
1
(t) =m
1
x
1
+g
1
(u
1
; _ u
1
) +g
2
(u
2
; _ u
2
) +g
3
(u
3
; _ u
3
) (7.4a)
f
2
(t) =m
2
x
2
+g
2
(u
2
; _ u
2
) +g
4
(u
4
; _ u
4
) +g
6
(u
6
; _ u
6
) (7.4b)
f
3
(t) =m
3
x
3
+g
3
(u
3
; _ u
3
) +g
4
(u
4
; _ u
4
) +g
5
( _ u
5
; _ u
5
) (7.4c)
_ g
1
(u
1
; _ u
1
) =
1

[A _ u
1
(j _ u
1
jjg
1
(u
1
; _ u
1
)j
n1
g
1
(u
1
; _ u
1
)
_ u
1
jg
1
(u
1
; _ u
1
)j
n
)] (7.4d)
_ g
2
(u
2
; _ u
2
) =
1

[A _ u
2
(j _ u
2
jjg
2
(u
2
; _ u
2
)j
n1
g
2
(u
2
; _ u
2
)
_ u
2
jg
2
(u
2
; _ u
2
)j
n
)] (7.4e)
_ g
3
(u
3
; _ u
3
) =
1

[A _ u
3
(j _ u
3
jjg
3
(u
3
; _ u
3
)j
n1
g
3
(u
3
; _ u
3
)
_ u
3
jg
3
(u
3
; _ u
3
)j
n
)] (7.4f)
_ g
4
(u
4
; _ u
4
) =
1

[A _ u
4
(j _ u
4
jjg
4
(u
4
; _ u
4
)j
n1
g
4
(u
4
; _ u
4
)
_ u
4
jg
4
(u
4
; _ u
4
)j
n
)] (7.4g)
_ g
5
(u
5
; _ u
5
) =
1

[A _ u
5
(j _ u
5
jjg
5
(u
5
; _ u
5
)j
n1
g
5
(u
5
; _ u
5
)
_ u
5
jg
5
(u
5
; _ u
5
)j
n
)] (7.4h)
_ g
6
(u
6
; _ u
6
) =
1

[A _ u
6
(j _ u
6
jjg
6
(u
6
; _ u
6
)j
n1
g
6
(u
6
; _ u
6
)
_ u
6
jg
6
(u
6
; _ u
6
)j
n
)] (7.4i)
where u
i
and _ u
i
are the relative displacement and the relevant velocity of the two ends of the
interconnecting linkg
i
determined according to the topology of the system as follows:
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
u
1
=x
1
(7.5a)
u
2
=x
2
x
1
(7.5b)
u
3
=x
3
x
1
(7.5c)
u
4
=x
3
x
2
(7.5d)
u
5
=x
5
(7.5e)
u
6
=x
6
(7.5f)
133
Excitation Parameters Values Bouc-Wen Values
Mean 0 n 1.00
Training S.D.*
t
1.00 1.00
Validation S.D.
v
0:50; 1:00; 1:50 A 1.00
Sampling Freq. t=T 0.05 1.00
Duration 80 Sec. 1.00
m
1
; m
2
; m
3
0.80, 2.00, 1.20
-0.50
*S.D.: Standard Deviation
Table 7.1: Properties of the Bouc-Wen model and the applied excitationsf
i
(t) for training and
validation.
The parameters of the Bouc-Wen models in all interconnecting components are identical, and
presented in Table (7.1)
Identification Formulation
The general system of differential equations that describes the three DOF system concerned in
this investigation, regardless of the type of the interconnecting componentsg
i
, is as follows:
134
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
f
1
(t) =m
1
x
1
+r
1
( _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
) (7.6a)
_ r
1
=q
1
(r
1
; _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
) (7.6b)
f
2
(t) =m
2
x
2
+r
2
( _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
) (7.6c)
_ r
2
=q
2
(r
2
; _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
) (7.6d)
f
3
(t) =m
3
x
1
+r
3
( _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
) (7.6e)
_ r
3
=q
3
(r
3
; _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
) (7.6f)
wherex
i
, _ x
i
and x
i
are the displacement, velocity, and acceleration, respectively,f
i
(t) is ex-
ternal excitation,m
i
is the mass, andr
i
and _ r
i
are the compound restoring force and its derivative
at DOFi. The quantityq
i
expresses the derivative of the restoring force _ r
i
in terms of the available
states. Two equations are dedicated for each DOFi in this system, and consequently, the set of all
equations together form a robust model for the system.
Excitation
Lumped masses of the MDOF system under investigation are stimulated by dissimilar zero mean
stationary white noise excitations. The applied training excitations at each DOF have identical
statistical properties which are presented in Table (7.1), and are produced using different seed
numbers. The external excitations stimulate the system sufficiently strong to disclose the yielding
region of the vibrating interconnecting components of the MDOF system.
135
Measurements
It is assumed that the external excitations, the lumped masses and their accelerations are known
from measurements. Hence, the displacements and velocities can be obtained by careful integra-
tion of accelerations. As the result, all statesx
i
, _ x
i
and x
i
,i = 1; 2; 3, are assumed to be available
from measurements.
7.3.2 Identification Procedure
The proposed methodology for the identification of nonlinear MDOF systems exhibiting non-
conservative dissipating behavior is summarized in the flowchart of Fig. (7.2). The identification
procedure consists of two major phases: training and validation.
136
Figure 7.2: The flowchart of the methodology presented for the identification of hysteretic MDOF
systems. The unidirectional arrows show the main stream of the training process.
137
Step 1 and 2: Stimulation
For training, first the synthetic training data set should be generated. Three different, but statisti-
cally similar, excitations stimulate the lumped masses at each degree of freedom.
Step 3: Response Calculation
The reference response of the structure under investigation is obtained using standard numerical
techniques, and the restoring force at each degree of freedom is calculated using the general
equation of motionr
i
=f
i
(t)m
i
x
i
. The estimate of the derivative of the restoring force is also
calculated using finite difference approximation as follows:
_ r
j
i
=
r
j+1
i
r
j
i
t
(7.7)
where _ r
j
i
is thejth datum in the array of the estimate of the derivative of the restoring forcer
i
, and
r
j+1
i
andr
j
i
are thej+1th andjth data in the array of the restoring forcer
i
, and t is the time step.
Step 4: Training Data Preparation
In order to determineq
i
in Eq. (7.6) that best describes the derivative of the restoring force _ r
i
,
the left hand side (LHS) of equations 7.6b, 7.6d and 7.6f, the states on the right hand side (RHS)
of equations 7.6b, 7.6d and 7.6f are input to GP. Therefore, the restoring force _ r
i
at DOFi, the
estimate of its derivative _ r
i
, along with displacement and velocity of all degrees of freedom form
the training data set. Hence, the evolution is specifically guided to obtain the differential operator
of the restoring force at each degree of freedom separately, and consequently should be repeated
138
3 times: each time one of the following data sets is utilized for training to attain a differential
operator:
(1) _ r
1
;r
1
; _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
(2) _ r
2
;r
2
; _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
(3) _ r
3
;r
3
; _ x
1
; _ x
2
; _ x
3
;x
1
;x
2
;x
3
(7.8)
The eventual aim is to discover three differential operators for the restoring forces _ r
i
, in terms
of the available states in the training data set to ultimately form the system of differential equa-
tions in Eq. (7.6f). Note that, because of the complexity of the dynamic of the system, and
the dependency of the response at each DOF on other elements, all the available displacement
and velocity states are included in each training data set. However, during the evolution, only
the most relevant state variables survive, and eventually construct the differential equations that
approximate the restoring force from the necessary building blocks.
Step 5: Genetic Programming Computation
GP uses available basis functions, algebraic operations and state variables to build a population,
evaluate them using the provided data set, and evolve them using evolutionary operators, to reach
an admissible solution.
Step 6: Un-optimized Model Formation
At the end of three rounds of evolution, three distinct differential equations are obtained for the
restoring forces at all degrees of freedom. Then, they are combined to form a system of differen-
tial equations. But it is still not optimized.
139
Step 7: Model’s Parameter Optimization
Due to the dependency of the system response on all coupled differential equations combined
together, the parameters are optimized using GAs, all together, to enhance the accuracy of the
response estimates. The computational burden at this stage is moderate, and acceptable. Thus,
the error between the model response and the reference response is the ideal choice of the cost
function. Hence, the cost function is defined as the summation of equally weighted errors between
the displacement, velocity and acceleration of the model response and the reference response at
every degree of freedom as follows:
=
3
X
i=1
(
kx
i
^ x
i
k
kx
i
k
+
k_ x
i

^
_ x
i
k
k_ x
i
k
+
k x
i

^
x
i
k
k x
i
k
) (7.9)
where x
i
; _ x
i
; x
i
are the displacement, velocity and acceleration, respectively, at DOFi and ^ x
i
;
^
_ x
i
;
^
x
i
are their estimates time-history arrays, respectively. The parameters of GA are listed in Table
(7.2).
GA Parameters Values
Method gaoptimset
Population Size 150
Initial Population From GP:c
PopInitRange [0:2c 5c]
Table 7.2: Properties of GA for parameter optimization.
140
Step 8: Final Optimized Model Completion
In the end, a system of differential equations is obtained that represent the MDOF system under
investigation.
Validation and Verification
The last phase of the identification scheme involves statistical analysis to predict the accuracy
of the discovered model when subjected to new excitations whose stimulations are substantially
different from the training excitations. This phase gages the generalizability of the discovered
model, and its applicability in new dynamical environments.
7.4 Identification of MDOF Hysteretic Systems Characterized by
Bouc-Wen Model
7.4.1 Training
The model shown in Fig. (7.1) is excited at all degrees of freedom by broad-band uncorrelated
forces to undergo horizontal motion. The excitations duration is 80 seconds, and their time histo-
ries are shown in the lower part of the time-history panels of Fig. (7.3), Fig. (7.4) and Fig. (7.5).
In order to show the complexity of the dynamical system under investigation, forces of the inter-
connecting linksg
i
are plotted against the relative displacement of their two ends in Fig. (7.7).
Note that, although 6 Bouc-Wen models are includded in this MDOF system, the estimate of the
governing system of differential equations in Eq. (7.6) involves only three equations for the com-
pound restoring forces at each DOF. Thus, the fairly simple proposed model is able to represent
141
the investigated complex system exhibiting nonlinear dissipative memory-dependent behavior. It
is also important to note that none of the information concerning the individual element restoring
forces (i.e.,g
i
) is used in the identification scheme. The plots in Fig. (7.7) are displayed to con-
firm that all of the 6 nonlinear elements are underlying significant hysteretic deformation. The
response of this stimulated 3-DOF system constructs 3 batches of training data sets according to
the list of Eq. (7.8). At the end of the training phase, the following system of differential equa-
tions is obtained, using the approach introduced earlier, to estimate the overall behavior of the
studied hysteretic MDOF system:
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
f
1
(t) = m
1
x
1
+r
1
_ r
1
= 2:80 _ x
1
_ x
2
r
1
abs( _ x
1
) 0:94 _ x
3
f
2
(t) = m
2
x
2
+r
2
_ r
2
= 2:75 _ x
2
0:19r
2
0:89 _ x
1
0:91 _ x
3
f
3
(t) = m
3
x
1
+r
3
_ r
3
= 2:80 _ x
3
_ x
2
r
3
abs( _ x
3
) 0:93 _ x
1
(7.10)
Solving these coupled differential equations by means of standard numerical techniques pro-
vides the response estimate. The time histories of the estimated response and the reference re-
sponse are shown in the time history panels of Fig. (7.3), Fig. (7.4) and Fig. (7.5), and the error
of the model response is reported in Fig. (7.6). The reasonable error associated with the estimates
confirm that the proposed identification technique yields reduced-order, reduced-complexity, op-
timized differential operators that effectively characterize the dynamics of the complex nonlinear
MDOF system with non-conservative dissipative traits.
142
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.3: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 1. The system is stimulated by 3 distinct external training excitations at all degrees-of-
freedom. The graphs of the reference and estimated response are plotted using solid and dotted
lines, respectively.
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.4: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 2. The system is stimulated by 3 distinct external training excitations at all degrees-of-
freedom. The graphs of the reference and estimated response are plotted using solid and dotted
lines, respectively.
143
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.5: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 3. The system is stimulated by 3 distinct external training excitations at all degrees-of-
freedom. The graphs of the reference and estimated response are plotted using solid and dotted
lines, respectively.
D3 D2 D1
0
0.1
0.2
0.3
0.4
0.5

Displacement
Velocity
Acceleration
Figure 7.6: The error of the estimated response of the 3-DOF hysteretic system at all 3 degrees-
of-freedom denoted as D1, D2, D3. The system is stimulated by 3 distinct external training
excitations at all degrees-of-freedom.
144
−1 1
−0.75
0.75
g
1
Force
x
1
(a) Force atg1 vs. displacementx1
−1 1
−0.75
0.75
g
2
Force
x
2
− x
1
(b) Force atg2 vs. relative displacementx2x1
−1 1
−0.75
0.75
g
3
Force
x
3
− x
1
(c) Force atg3 vs. relative displacementx3x1
−1 1
−0.75
0.75
g
4
Force
x
3
− x
2
(d) Force atg4 vs. relative displacementx3x2
−1 1
−0.75
0.75
g
5
Force
x
3
(e) Force atg5 vs. displacementx3
−1 1
−0.75
0.75
g
6
Force
x
2
(f) Force atg6 vs. displacementx2
Figure 7.7: Phase plots of the force at interconnecting elements of g
1
to g
6
vs. the relative
displacement of the corresponding ends when the 3-DoF hysteretic system is subjected to distinct
training excitations with
v
= 1:0 at all degrees of freedom.
145
7.5 Validation and Discussion: Statistical Analysis
Validation analysis is performed to verify the applicability of the proposed mathematical model
with unfamiliar dynamical environment, i.e., when the intensity of the stimulation and the time
history are different from which was used for training. Moreover, statistical analysis helps to ver-
ify the independence of the accuracy of the model estimates from the training excitation, and also
to evaluate the overall performance of the discovered model under different dynamical conditions.
To this end, 20 different excitation sets are generated using distinct seed numbers, for three dif-
ferent stimulation levels. The stimulation intensities are chosen to be less, equal, and more than
the intensity of the training excitation. The intensity level is controlled by the standard deviation
of the generated exciting forces
v
. the The statistical properties of these 60 excitation sets are
summarized in Table (7.1). The standard deviation of these validation excitations is denoted as

v
. Note that three different forces exist in each excitation set to be applied to all 3 degrees of
freedom of the MDOF structure in Fig. (7.1). The actual response of the structure excited by
these 360 excitations is compared to the mathematical model solution, and the histograms of
the associated error are presented.
Figure (7.8) shows the error distribution associated with the difference between the reference
response and the model response when the stimulation intensity of the exciting forces is low

v
= 0:50, as oppose to the training excitation
t
= 1:00. The error is calculated for displace-
mentx, velocity _ x and acceleration x responses of every DOFi separately. The associated normal
distribution of the data are also superposed on the histograms. The results reveal that the most
accurate estimate belongs to the acceleration and the least belongs to the displacement. How-
ever, considering the reduction of the validation intensity by half, and employing a reduced-order
reduced-complexity model of Eq. (7.10) for the hysteretic MDOF system under discussion, the
146
error is reasonable.
With the stimulation intensity of the training force at
t
= 1:00, similar histograms are presented
for the validation excitations with similar stimulation intensity at
v
= 1:00, and high stimula-
tion intensity at
v
= 1:50 in Fig. (7.9) and Fig. (7.10), respectively. The computed average
error in each scenario at every DOF is compared in Fig. (7.11). This error bar implies that, as
we would expect, the model yields the most accurate response estimates when the severity of the
training and validation excitations is comparable. However, the model still predicts the response
reasonably well even under lower or higher stimulation intensities.
147
0 100
0
8
DOF1
Distribution
Error
(a) Displacement of DOF 1
0 100
0
8
DOF2
Distribution
Error
(b) Displacement of DOF 2
0 100
0
8
DOF3
Distribution
Error
(c) Displacement of DOF 3
0 50
0
8
DOF1
Distribution
Error
(d) Velocity of DOF 1
0 50
0
8
DOF2
Distribution
Error
(e) Velocity of DOF 2
0 50
0
8
DOF3
Distribution
Error
(f) Velocity of DOF 3
0 20
0
8
DOF1
Distribution
Error
(g) Acceleration of DOF 1
0 20
0
8
DOF2
Distribution
Error
(h) Acceleration of DOF 2
0 20
0
8
DOF3
Distribution
Error
(i) Acceleration of DOF 3
Figure 7.8: Distribution of the percentage of the model response error at every degree of freedom
when subjected to 20 validation excitations with
v
= 0:5. Having
v
= 0:5, these excitations are
less intense than the identification excitation having
t
= 1:0. The vertical axes show the number
of samples and the horizontal axes show the corresponding error value of the same number of
excitations in the bin. To obtain an estimate about the average error, the normal distribution of the
data are also provided in each figure.
148
0 100
0
8
DOF1
Distribution
Error
(a) Displacement of DOF 1
0 100
0
8
DOF2
Distribution
Error
(b) Displacement of DOF 2
0 100
0
8
DOF3
Distribution
Error
(c) Displacement of DOF 3
0 50
0
8
DOF1
Distribution
Error
(d) Velocity of DOF 1
0 50
0
8
DOF2
Distribution
Error
(e) Velocity of DOF 2
0 50
0
8
DOF3
Distribution
Error
(f) Velocity of DOF 3
0 20
0
8
DOF1
Distribution
Error
(g) Acceleration of DOF 1
0 20
0
8
DOF2
Distribution
Error
(h) Acceleration of DOF 2
0 20
0
8
DOF3
Distribution
Error
(i) Acceleration of DOF 3
Figure 7.9: Distribution of the percentage of the model response error at every degree of freedom
when subjected to 20 validation excitations with
v
= 1:0. Having
v
= 1:0, these excitations
are as intense as the identification excitation having
t
= 1:0. The vertical axes show the number
of samples and the horizontal axes show the corresponding error value of the same number of
excitations in the bin. To obtain an estimate about the average error, the normal distribution of the
data are also provided in each figure.
149
0 100
0
8
DOF1
Distribution
Error
(a) Displacement of DOF 1
0 100
0
8
DOF2
Distribution
Error
(b) Displacement of DOF 2
0 100
0
8
DOF3
Distribution
Error
(c) Displacement of DOF 3
0 50
0
8
DOF1
Distribution
Error
(d) Velocity of DOF 1
0 50
0
8
DOF2
Distribution
Error
(e) Velocity of DOF 2
0 50
0
8
DOF3
Distribution
Error
(f) Velocity of DOF 3
0 20
0
8
DOF1
Distribution
Error
(g) Acceleration of DOF 1
0 20
0
8
DOF2
Distribution
Error
(h) Acceleration of DOF 2
0 20
0
8
DoF3
Distribution
Error
(i) Acceleration of DOF 3
Figure 7.10: Distribution of the percentage of the model response error at every degree of freedom
when subjected to 20 validation excitations with
v
= 1:5. Having
v
= 1:5, these excitations
are 1.5 times more intense than the identification excitation having
t
= 1:0. The vertical axes
show the number of samples and the horizontal axes show the corresponding error value of the
same number of excitations in the bin. To obtain an estimate about the average error, the normal
distribution of the data are also provided in each figure.
150
D1
D2
D3
D1
D2
D3
D1
D2
D3
0
0.1
0.2
0.3
0.4
0.5

Excitation3
Excitation2

Excitation1
Displacement
Velocity
Acceleration
Figure 7.11: The average of the validation error associated with the model response subject to
the validation excitations with respect to the reference response. The stimulation level of the
validation excitations are substantially different from the training excitation: Excitation 1 (
v
=
0:5) is less intense than the training excitation (
t
= 1:0), Excitation 2 (
v
= 1:0) is as intense
as the training excitation (
t
= 1:0) and Excitation 3 (
v
= 1:5) is more intense than the training
excitation (
t
= 1:0).
In each equal-in-intensity group, the data set whose error falls on the average of the error
across all data sets is selected to represent a typical validation excitation. Thus, three excitation
with
v
= 0:50,
v
= 1:00 and
v
= 1:50 are selected to investigate the accuracy of the predic-
tions in more detail. The reference response of the MDOF system subjected to these validation
excitations is compared to the model response.
Generated for validation exciting force with
v
= 0:50, Fig. (7.12) shows the actual forces of
the system in the linksg
1
g
6
. Compared to the training phase plots of Fig. (7.7), the hysteretic
151
loops of the less severe excitation
v
= 0:50 are smaller than that of the training. Figure (7.13),
Fig. (7.14) and Fig. (7.15) also compare the associated time histories of the states.
152
−0.5 0.5
−0.75
0.75
g
1
Force
x
1
(a) Force atg1 vs. displacementx1
−0.5 0.5
−0.75
0.75
g
2
Force
x
2
− x
1
(b) Force atg2 vs. relative displacementx2x1
−0.5 0.5
−0.75
0.75
g
3
Force
x
3
− x
1
(c) Force atg3 vs. relative displacementx3x1
−0.5 0.5
−0.75
0.75
g
4
Force
x
3
− x
2
(d) Force atg4 vs. relative displacementx3x2
−0.5 0.5
−0.75
0.75
g
5
Force
x
3
(e) Force atg5 vs. displacementx3
−0.5 0.5
−0.75
0.75
g
6
Force
x
2
(f) Force atg6 vs. displacementx2
Figure 7.12: Phase plots of the force at interconnecting elements of g
1
to g
6
vs. the relative
displacement of the corresponding ends when the 3-DOF hysteretic system is subjected to distinct
validation excitations with
v
= 0:5 at all degrees of freedom. These plots are obtained from
reference responses for providing insight to the hysteretic behavior of interconnecting elements
g
1
tog
6
.
153
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.13: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 1. The system is stimulated by 3 distinct external validation excitations with
v
= 0:5 at
all degrees-of-freedom which is less intense than the training excitation with
t
= 1:0. The graphs
of the reference and estimated response are plotted using solid and dotted lines, respectively.
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.14: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 2. The system is stimulated by 3 distinct external validation excitations with
v
= 0:5 at
all degrees-of-freedom which is less intense than the training excitation with
t
= 1:0. The graphs
of the reference and estimated response are plotted using solid and dotted lines, respectively.
154
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.15: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 3. The system is stimulated by 3 distinct external validation excitations with
v
= 0:5 at
all degrees-of-freedom which is less intense than the training excitation with
t
= 1:0. The graphs
of the reference and estimated response are plotted using solid and dotted lines, respectively.
Similarly, the time-history plots of the system and model response under the validation exci-
tations for
v
= 1:00 are presented in Fig. (7.16), Fig. (7.17) and Fig. (7.18), and for
v
= 1:50
are presented in Fig. (7.20), Fig. (7.21) and Fig. (7.22), respectively. The phase plots of the
actual forces of interconnecting links ofg
1
g
6
are also shown in Fig. (7.19) and Fig. (7.23),
respectively. The graphical illustrations imply that, when the intensity of the validation excitation
increases to
v
= 1:50, i.e., 50% higher than the training excitation
t
= 1:00, the hysteretic loops
enlarge significantly. However, the presented system of differential equations in Eq. (7.10) yield
reasonably accurate estimates even when the intensity of the stimulations increases significantly.
Hence, this model is reliable when the dynamic environment undergoes significant changes, and
future excitations will not necessarily stimulate the system within the range of the training excita-
tion, under which the model was obtained. It is worth mentioning that the training excitation was
intense enough to significantly stimulate the system to its yielding region to provide rich data for
155
training. Otherwise, slight perturbations can not reveal the true physical nature of the underling
structure.
In summary, the findings suggest that the proposed method yields reasonably accurate high
fidelity surrogate computational model in the form of systems of differential equations for the
hysteretic MDOF system with non-conservative dissipative behavior under investigation.
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.16: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 1. The system is stimulated by 3 distinct external validation excitations with
v
= 1:0 at
all degrees-of-freedom which is as intense as the training excitation with
t
= 1:0. The graphs of
the reference and estimated response are plotted using solid and dotted lines, respectively.
156
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.17: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 2. The system is stimulated by 3 distinct external validation excitations with
v
= 1:0 at
all degrees-of-freedom which is as intense as the training excitation with
t
= 1:0. The graphs of
the reference and estimated response are plotted using solid and dotted lines, respectively.
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.18: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 3. The system is stimulated by 3 distinct external validation excitations with
v
= 1:0 at
all degrees-of-freedom which is as intense as the training excitation with
t
= 1:0. The graphs of
the reference and estimated response are plotted using solid and dotted lines, respectively.
157
−1 1
−0.75
0.75
g
1
Force
x
1
(a) Force atg1 vs. displacementx1
−1 1
−0.75
0.75
g
2
Force
x
2
− x
1
(b) Force atg2 vs. relative displacementx2x1
−1 1
−0.75
0.75
g
3
Force
x
3
− x
1
(c) Force atg3 vs. relative displacementx3x1
−1 1
−0.75
0.75
g
4
Force
x
3
− x
2
(d) Force atg4 vs. relative displacementx3x2
−1 1
−0.75
0.75
g
5
Force
x
3
(e) Force atg5 vs. displacementx3
−1 1
−0.75
0.75
g
6
Force
x
2
(f) Force atg6 vs. displacementx2
Figure 7.19: Phase plots of the force at interconnecting elements of g
1
to g
6
vs. the relative
displacement of the corresponding ends when the 3-DoF hysteretic system is subjected to distinct
validation excitations with
v
= 1:0 at all degrees of freedom. These plots are obtained from
reference responses for providing insight to the hysteretic behavior of interconnecting elements
g
1
tog
6
.
158
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.20: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 1. The system is stimulated by 3 distinct external validation excitations with
v
= 1:0 at
all degrees-of-freedom which is as intense as the training excitation with
t
= 1:0. The graphs of
the reference and estimated response are plotted using solid and dotted lines, respectively.
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.21: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 2. The system is stimulated by 3 distinct external validation excitations with
v
= 1:0 at
all degrees-of-freedom which is as intense as the training excitation with
t
= 1:0. The graphs of
the reference and estimated response are plotted using solid and dotted lines, respectively.
159
0 80
−5
5
Excitation
Time (sec)
0 80
−6
6
Acceleration
Time (sec)
0 80
−1.5
1.5
Velocity
Time (sec)
0 80
−1.5
1.5
Displacement
Time (sec)
0 80
−2
2
Restoring Force
Time (sec)
Figure 7.22: Comparison of the reference and estimated response of the 3-DOF hysteretic system
at DOF 3. The system is stimulated by 3 distinct external validation excitations with
v
= 1:5
at all degrees-of-freedom which is 1.5 times more intense than the training excitation with
t
=
1:0. The graphs of the reference and estimated response are plotted using solid and dotted lines,
respectively.
160
−1.5 1.5
−0.75
0.75
g
1
Force
x
1
(a) Force atg1 vs. displacementx1
−1.5 1.5
−0.75
0.75
g
2
Force
x
2
− x
1
(b) Force atg2 vs. relative displacementx2x1
−1.5 1.5
−0.75
0.75
g
3
Force
x
3
− x
1
(c) Force atg3 vs. relative displacementx3x1
−1.5 1.5
−0.75
0.75
g
4
Force
x
3
− x
2
(d) Force atg4 vs. relative displacementx3x2
−1.5 1.5
−0.75
0.75
g
5
Force
x
3
(e) Force atg5 vs. displacementx3
−1.5 1.5
−0.75
0.75
g
6
Force
x
2
(f) Force atg6 vs. displacementx2
Figure 7.23: Phase plots of the force at interconnecting elements of g
1
to g
6
vs. the relative
displacement of the corresponding ends when the 3-DoF hysteretic system is subjected to distinct
validation excitations with
v
= 1:5 at all degrees of freedom. These plots are obtained from
reference responses for providing insight to the hysteretic behavior of interconnecting elements
g
1
tog
6
.
161
7.6 Summary and Conclusion
By building on the revolutionary advances in the field of evolutionary computing, this chapter
introduces a robust methodology that yields optimized computational models in the form of sys-
tems of differential equations that represent nonlinear MDOF systems with memory-dependent
energy-dissipative characteristics. The proposed technique utilizes Genetic Programming to dis-
cover suitable sets of differential operators through evolving their tree variant. This approach
also incorporates stochastic parameter optimization techniques involving Genetic Algorithms to
optimize the model parameters using the error of the model response as the cost function. The
performance of the discovered model in future dynamical situations is assessed through the sta-
tistical analysis of different validation forces having various time histories and stimulation levels.
Reasonably accurate estimates of the discovered model subjected to new external forces verifies
the robustness of the proposed approach. Although the validation forces were substantially dif-
ferent from the training excitations, and despite the complexity of the hysteretic phenomena, the
proposed methodology yields simple parsimonious computational models in the form of differ-
ential equations that are able to effectively characterize the behavior of nonlinear MDOF systems
with path-dependent dissipative traits.
162
Chapter 8
Summary and Conclusion
T
HIS study builds on advancements in the field of Computational Intelligence to develop a
generalized procedure for the identification of complex nonlinear phenomena. The pro-
posed methodology combines stochastic optimization methods utilizing Genetic Programming for
evolving the structure of differential equations, combined with nonlinear evolutionary optimiza-
tion methods using Genetic Algorithms for optimizing the parameters of differential equations,
and Computer Algebra techniques that involve symbolic manipulation of mathematical expres-
sions for reducing their complexity during the course of evolution, to “discover” parsimonious
differential operators that represent the system of interest. The proposed scheme requires input
and output data only, without postulating any model-class in advance. Employing the suggested
identification scheme, a variety of highly nonlinear and complex phenomena, ranging from me-
chanical to biological systems, are investigated, and the performance of the discovered models in
not-previously experienced dynamical environments is assessed.
163
Chapter 2 provides an overview of Genetic Programming and its adaptation for identification
purposes in the context of this research, an overview of Genetic Algorithms for parameter opti-
mization, and an overview of Linear Algebra techniques for bloat control. This chapter founds
the basis for presenting the proposed identification scheme in the following chapters.
Chapter 3 provides the proposed approach for the identification of a class of dynamical sys-
tems that exhibit polynomial nonlinearities in their response. The details of the identification
scheme and the associated flowchart are presented. The Duffing, van der Pol and noisy-Duffing
van der Pol oscillators are identified using the proposed technique. It is shown that this tech-
nique yields computational models in the form of differential operators that provide an optimum
match to the governing differential equations of the target complex nonlinear systems, and sub-
sequently disclose the correct nature of the investigated systems. Findings verify the reliability
of the proposed identification scheme, and the robustness of the resulting models. Therefore, the
method of this chapter provides a robust methodology for developing high-fidelity reduced-order,
reduced-complexity, computational models that reflect the correct “physics” of the underlying
phenomena.
Chapter 4 provides the proposed scheme for the identification of dynamical phenomena that
exhibit discontinuity in their response. The details of the identification procedure and the asso-
ciated flowchart are presented. Incorporating discontinuous and non-differentiable functions in
the mathematical representation, the suggested algorithm can discover differential operators that
match the exact governing differential equations to characterize the behavior of nonlinear sys-
tems with sharp discontinuities in their response. Moreover, this technique reveals an accurate
single-expression, with direct physical interpretation, that represents the governing multi-region
(response domain) equations of systems that incorporate certain classes of nonlinear phenomena
164
(such as yielding). Yet, unlike many conventional non-parametric techniques whose approxi-
mations result in undesirable oscillations around unsmooth points, automatic incorporation of
discontinuous basis functions in this approach eliminates the need of such approximations and
their concomitant errors. The capabilities and generalization extent of the proposed technique is
assessed even when polluted data are employed. Due to the stochastic nature of the proposed
identification technique, Monte carlo simulation is performed to confirm the necessity of the in-
tegration of all three essential components of the proposed approach i.e., Genetic Programming,
Genetic Algorithms, and Linear Algebra, in an efficient evolutionary process to reach robust par-
simonious models.
Chapter 5 provides the robust evolutionary based hybrid method for discovering the differen-
tial equations of nonlinear systems that exhibit hysteretic behavior. The details of the proposed
approach and its flowchart are presented. Utilizing discontinuous basis functions, the proposed
method provides parsimonious models to represent governing differential equations of the under-
lying hysteretic oscillators. This technique is utilized to identify the Bouc-Wen model and the
bilinear hysteretic oscillator — both exhibit abrupt change in their memory-dependent response.
The representative models are subjected to validation excitations that are substantially different
from the probing signals to confirm the generalizability of the models for different dynamical phe-
nomena. The results verify the effectiveness of the approach and the accuracy of the subsequent
differential operators that characterize the behavior of the studied hysteretic systems.
Chapter 6 presents the nonparametric modeling of human spine using the technique proposed
herein by means of experimental data. The eventual goal is to discover a differential equation
that relates the angular velocity to the measured loads and states of the spine. The method of the
experimentation, and the properties of the available data sets are described. It is shown that the
165
proposed hybrid technique yields differential operators with optimized parameters that accurately
estimate the spine’s movement.
Chapter 7 presents the proposed scheme for the identification of multi-degree-of-freedom sys-
tems exhibiting non-conservative dissipative behavior. The effectiveness of the approach is tested
using a benchmark hysteretic MDOF system governed by the Bouc-Wen models. Genetic Pro-
gramming yields distinct differential equations for the compound restoring force at each degree-
of-freedom. Then, Genetic Algorithms are employed to optimize the parameters of the entire
system of differential equations that represent the MDOF structure under discussion. The cost
function at this stage is defined as the error of the response estimates at all degrees of freedom. A
statistical analysis is performed to investigate the performance of the discovered optimized model
in future dynamical environments. Variety of validation excitations with different intensity lev-
els are applied to assess the generalizability of the discovered model, when not-previously-seen
excitations stimulate the system, more than or less than the excitations employed for training.
Findings suggest that the proposed methodology results in simple and robust models in the form
of a system of differential equations that can effectively represent nonlinear MDOF systems with
memory-dependent dissipative traits in future dynamical environments, and provide reasonably
accurate estimates.
166
Bibliography
[1] Sophia N. Sangiorgio, Sean L. Borkowski, Richard E. Bowen, Anthony A. Scaduto,
Nathan L. Frost, and Edward Ebramzadeh. Quantification of increase in three-dimensional
spine flexibility following sequential ponte osteotomies in a cadaveric model. Spine De-
formity, 1(3):171 – 178, 2013.
[2] Sophia N Sangiorgio, Hormoz Sheikh, Sean L Borkowski, Larry Khoo, Christopher R
Warren, and Edward Ebramzadeh. Comparison of three posterior dynamic stabilization
devices. Spine, 36(19):E1251–E1258, 2011.
[3] Charles Darwin. The origin of species by means of natural selection or the preservation of
favored races in the struggle for life. Murray, London, 1859.
[4] Riccardo Leardi. Genetic algorithms in chemometrics and chemistry: a review. Journal of
chemometrics, 15(7):559–569, 2001.
[5] B. Engel. Problem solving with genetic algorithms. Mathematics and Computer
Education, 34(3):277 – 88, 2000. problem solving;genetic algorithms;complex prob-
lems;evolution;potential solution;chromosome;optimal solution;Digital Dimmer Switch
problem;computational procedures;mathematical function;volume;cardboard box;BASIC
code;selection;mating;mutation;.
[6] G´ abor Renner and Anik´ o Ek´ art. Genetic algorithms in computer aided design. Computer-
Aided Design, 35(8):709–726, 2003.
[7] Carlos Manuel Mira da Fonseca. Multiobjective genetic algorithms with application to
control engineering problems. Citeseer, 1995.
[8] Ricardo Salem Zebulum, Marco Aurelio Pacheco, and Marley Maria Be Vellasco. Evolu-
tionary electronics: automatic design of electronic circuits and systems by genetic algo-
rithms, volume 22. CRC press, 2001.
[9] Shu-Heng Chen. Genetic algorithms and genetic programming in computational finance.
Springer, 2002.
[10] Andrew Horner and D Goldberg. Genetic algorithms and computer-assisted music compo-
sition. Urbana, 51(61801):14, 1991.
167
[11] Thomas B¨ ack and Hans-Paul Schwefel. An overview of evolutionary algorithms for pa-
rameter optimization. Evolutionary computation, 1(1):1–23, 1993.
[12] Mitsuo Gen and Runwei Cheng. Genetic algorithms and engineering optimization, vol-
ume 7. John Wiley & Sons, 2000.
[13] Jhon. Koza. Genetic Programming: On the Programming of Computers by Means of Nat-
ural Selection. MIT Press, Cambridge, MA, 1992.
[14] W.A. Tackett. Genetic programming for feature discovery and image discrimination. pages
303 – 9, San Mateo, CA, USA, 1993.
[15] F. Gruau. Automatic definition of modular neural networks. Adaptive Behavior, 3(2):151
– 83, 1994.
[16] L. Spector and A. Alpern. Induction and recapitulation of deep musical structure. Montreal,
Quebec, Canada, 1995.
[17] Simon Handley. Predicting whether or not a nucleic acid sequence is an e. coli promoter
region using genetic programming. pages 122 – 127, Herndon, V A, USA, 1995.
[18] Man Leung Wong, Kwong Sak Leung, and J.C.Y . Cheng. Discovering knowledge from
noisy databases using genetic programming. Journal of the American Society for Informa-
tion Science, 51(9):870 – 81, 2000.
[19] D. Howard and S.C. Roberts. The prediction of journey times on motorways using genetic
programming. pages 210 – 21, Berlin, Germany, 2002.
[20] Ying Becker, Peng Fei, and Anna Lester. Stock selection: An innovative application of
genetic programming methodology. In Rick Riolo, Terence Soule, and Bill Worzel, editors,
Genetic Programming Theory and Practice IV, Genetic and Evolutionary Computation,
pages 315–334. Springer US, 2007.
[21] Gregory. S. Hornby, Jason D. Lohn, and Derek S. Linden. Computer-automated evolution
of an x-band antenna for nasa’s space technology 5 mission. Evol. Comput., 19(1):1–23,
2011.
[22] Dimitris Dracopoulos and Dimitrios Effraimidis. Genetic programming for generalised
helicopter hovering control. In Genetic Programming, volume 7244 of Lecture Notes in
Computer Science, pages 25–36. Springer Berlin / Heidelberg, 2012.
[23] Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental
data. science, 324(5923):81–85, 2009.
[24] Sara Silva, Orlando Anunciac ¸˜ ao, and Marco Lotz. A comparison of machine learning meth-
ods for the prediction of breast cancer. In Evolutionary Computation, Machine Learning
and Data Mining in Bioinformatics, pages 159–170. Springer, 2011.
168
[25] R. A. FISHER. The Genetical Theory of Natural Selection. Diver, N.Y ., 1958.
[26] L. J. Fogel, A. J. Owens, and Walsh M. J. Artificial intelligence through simulated evolu-
tion. Wiley, New York, 1966.
[27] H. P. Schwefel. Kybernetische evolution als strategie der experimentellen forschung in der
stromungstechnik, 1973.
[28] I. Rechenberg. Evolutionsstrategie: Optimierung technischer systeme nach prinzipien der
biologischen evolution. Frommann-Holzboog, 1973.
[29] J. H. Holland. Adaptation in natural and artificial systems. The University of Michigan
Press, Ann Arbor, 1975.
[30] D. Thierens and D. Goldberg. Convergence models of genetic algorithm selection schemes.
pages 119 – 29, Berlin, Germany, 1994.
[31] P.J. Fleming and R.C. Purshouse. Evolutionary algorithms in control systems engineering:
a survey. Control Engineering Practice, 10(11):1223 – 41, 2002.
[32] J.E. Baker. Reducing bias and inefficiency in the selection algorithm. pages 14 – 21,
Hillsdale, NJ, USA, 1987.
[33] D.E. Goldberg and K. Deb. A comparative analysis of selection schemes used in genetic
algorithms, pages 69–93. Foundations of genetic algorithms. 1991.
[34] W. B. Langdon and R. Poli. Genetic programming bloat with dynamic fitness. In
W. Banzhaf, R. Poli, M. Schoenauer, and T. C. Fogarty, editors, The First European Work-
shop on Genetic Programming, volume 1391 of LNCS, page 96112. Springer Paris, 1998.
[35] W. Banzhaf and W.B. Langdon. Some considerations on the reason for bloat. Genetic
Programming and Evolvable Machines, 3(1):81 – 91, 2002.
[36] S. Silva and E. Costa. Dynamic limits for bloat control: variations on size and depth. pages
666 – 77, Berlin, Germany, 2004.
[37] S. Luke and L. Panait. Lexicographic parsimony pressure. In W. B. Langdon et al., editor,
GECCO 2002: Proceedings of the Genetic and Evolutionary Computation Conference,
pages 829–836, New York, USA, 2002. Morgan Kaufmann Publishers.
[38] C. Gagne, M. Schoenauer, M. Parizeau, and M. Tomassini. Genetic programming, valida-
tion sets, and parsimony pressure. pages 109 – 20, Berlin, Germany, 2006.
[39] P Ibanez. Identification of dynamic parameters of linear and non-linear structural models
from experimental data. Nuclear Engineering and Design, 25(1):30–41, 1973.
[40] S.F. Masri and T.K. Caughey. A nonparametric identification technique for nonlinear dy-
namic problems. Journal of Applied Mechanics, 46(2):433 – 447, 1979.
169
[41] S.F. Masri, H. Sassi, and T.K. Caughey. Nonparametric identification of nearly arbitrary
nonlinear systems. V 49(N 3):619 – 628, 1982.
[42] S.F. Marsi, G.A. Bekey, H. Sassi, and T.K. Caughey. Non-parametric identification of a
class of non-linear multidegree dynamic systems. Earthquake Engineering and Structural
Dynamics, 10(1):1 – 30, 1982.
[43] H. B. Yun, F. Tasbighoo, S.F. Masri, J.P. Caffrey, R.W. Wolfe, N. Makris, and C. Black.
Comparison of modeling approaches for full-scale nonlinear viscous dampers. Journal of
Vibration and Control, 14(1-2):51 – 76, 2008.
[44] H. B. Yun, Sami F. Masri, Raymond W. Wolfe, and Gianmario Benzoni. Data-driven
methodologies for change detection in large-scale nonlinear dampers with noisy measure-
ments. Journal of Sound and Vibration, 322(1-2):336 – 357, 2009.
[45] C.H. Kim, A.R. Mijar, and J.S. Arora. Development of simplified models for design and op-
timization of automotive structures for crashworthiness. Structural and Multidisciplinary
Optimization, 22(4):307 – 21, 2001.
[46] M.S. Allen, H. Sumali, and D.S. Epp. Piecewise-linear restoring force surfaces for semi-
nonparametric identification of nonlinear systems. Nonlinear Dynamics, 54(1-2):123 – 35,
2008.
[47] Edward F. Crawley and Allan C. Aubert. Identification of nonlinear structural elements by
force-state mapping. AIAA journal, 24(1):155 – 162, 1986.
[48] Francesco Benedettini, Danilo Capecchi, and Fabrizio Vestroni. Identification of hysteretic
oscillators under earthquake loading by nonparametric models. Journal of Engineering
Mechanics, 121(5):606 – 612, 1995.
[49] K. Shin and J.K. Hammond. Pseudo force-state mapping method: incorporation of the
embedding and force-state mapping methods. Journal of Sound and Vibration, 211(5):918
– 22, 1998.
[50] M. Haroon, E.A. Douglas, Y .W. Luk, and A.A. Ferri. A time and frequency domain ap-
proach for identifying nonlinear mechanical system models in the absence of an input mea-
surement. Journal of Sound and Vibration, 283(3-5):1137 – 55, 2005.
[51] M. Ajjan Al-Hadid and J.R. Wright. Developments in the force-state mapping technique
for non-linear systems and the extension to the location of non-linear elements in a lumped-
parameter system. Mechanical Systems and Signal Processing, 3(3):269 – 290, 1989.
[52] G. Box. Time series analysis: Forecasting and control. Holden-Day, San Francisco, 1970.
[53] I.J. Leontaritis and S.A. Billings. Input-output parametric models for non-linear systems.
part i: Deterministic non-linear systems. International Journal of Control, 41(2):303 –
328, 1985.
170
[54] I.J. Leontaritis and S.A. Billings. Input-output parametric models for non-linear systems.
part ii: Stochastic non-linear systems. International Journal of Control, 41(2):329 – 344,
1985.
[55] Li Haitao and Akira Mita. Damage indicator defined as the distance between arma models
for structural health monitoring. Structural Control and Health Monitoring, 15(7):992 –
1005, 2008.
[56] G. Kerschen, K. Worden, A.F. Vakakis, and J.-C. Golinval. Past, present and future of
nonlinear system identification in structural dynamics. Mechanical Systems and Signal
Processing, 20(3):505 – 92, 2006.
[57] Tomoo Saito and James L. Beck. Bayesian model selection for arx models and its appli-
cation to structural health monitoring. Earthquake Engineering and Structural Dynamics,
39(15):1737 – 1759, 2010.
[58] S.F. Masri, A.G. Chassiakos, and T.K. Caughey. Structure-unknown non-linear dynamic
systems: identification through neural networks. Smart Materials and Structures, 1(1):45
– 56, 1992.
[59] S.F. Masri, A.G. Chassiakos, and T.K. Caughey. Identification of nonlinear dynamic sys-
tems using neural networks. Transactions of the ASME. Journal of Applied Mechanics,
60(1):123 – 33, 1993.
[60] K. Worden and G.R. Tomlinson. Modeling and classification of non-linear systems using
neural networks-i. simulation. Mechanical Systems and Signal Processing, 8(3):319 – 56,
1994.
[61] K. Worden, G.R. Tomlinson, W. Lim, and G. Sauer. Modeling and classification of non-
linear systems using neural networks .ii. a preliminary experiment. Mechanical Systems
and Signal Processing, 8(4):395 – 419, 1994.
[62] A.G. Chassiakos and S.F. Masri. Modelling unknown structural systems through the use of
neural networks. Earthquake Engineering and Structural Dynamics, 25(2):[d]117 – 128,
1996.
[63] A.W. Smyth, S.F. Masri, A.G. Chassiakos, and T.K. Caughey. On-line parametric iden-
tification of mdof nonlinear hysteretic systems. Journal of Engineering Mechanics,
125(2):133 – 142, 1999.
[64] S.F. Masri, A.W. Smyth, A.G. Chassiakos, T.K. Caughey, and N.F. Hunter. Application
of neural networks for detection of changes in nonlinear systems. Journal of Engineering
Mechanics, 126(7):666 – 676, 2000.
[65] S.F. Masri, A.W. Smyth, A.G. Chassiakos, T.K. Caughey, and N.F. Hunter. Application
of neural networks for detection of changes in nonlinear systems. Journal of Engineering
Mechanics, 126(7):666 – 676, 2000.
171
[66] J. S. Pei and Andrew W. Smyth. New approach to designing multilayer feedforward neural
network architecture for modeling nonlinear restoring forces. i: Formulation. Journal of
Engineering Mechanics, 132(12):1290 – 1300, 2006.
[67] J. S. Pei and Andrew W. Smyth. New approach to designing multilayer feedforward neural
network architecture for modeling nonlinear restoring forces. ii: Applications. Journal of
Engineering Mechanics, 132(12):1301 – 1312, 2006.
[68] S.F. Masri, J.P. Caffrey, T.K. Caughey, A.W. Smyth, and A.G. Chassiakos. Identification
of the state equation in complex non-linear systems. International Journal of Non-Linear
Mechanics, 39(7):1111 – 27, 2004.
[69] G.a. Baker and P. Graves-Morris. Pad Approximations. Cambridge University Press, 2nd
edition, 1996.
[70] L. Emmel, S.M. Kaber, and Y . Maday. Pade-jacobi filtering for spectral approximations of
discontinuous solutions. Numerical Algorithms, 33(1-4):251 – 64, 2003.
[71] S.M. Kaber and Y . Maday. Analysis of some pade-chebyshev approximants. SIAM Journal
on Numerical Analysis, 43(1):437 – 454, 2005.
[72] J.S. Hesthaven, S.M. Kaber, and L. Lurati. Pade-legendre interpolants for gibbs recon-
struction. Journal of Scientific Computing, 28(2-3):337 – 59, 2006.
[73] J. Foster and F. B. Richards. The gibbs phenomenon for piecewise-linear approximation.
The American Mathematical Monthly, 98(1):47–49, 1991.
[74] T.A. Driscoll and B. Fornberg. A pade-based algorithm for overcoming the gibbs phe-
nomenon. Numerical Algorithms, 26(1):77 – 92, 2001.
[75] B.D. Shizgal and Jae hun Jung. Towards the resolution of the gibbs phenomena. Journal
of Computational and Applied Mathematics, 161(1):41 – 65, 2003.
[76] T. Chantrasmi, A. Doostan, and G. Iaccarino. Pade-legendre approximants for uncer-
tainty analysis with discontinuous response surfaces. Journal of Computational Physics,
228(19):7159 – 80, 2009.
[77] TL Lew, AB Spencer, F Scarpa, K Worden, A Rutherford, and F Hemez. Identification
of response surface models using genetic programming. Mechanical Systems and Signal
Processing, 20(8):1819–1831, 2006.
[78] T.K. Caughey. Random excitation of system with bilinear hysteresis. American Society of
Mechanical Engineers – Transactions – Journal of Applied Mechanics Series E, 27(4):649
– 652, 1960.
[79] W.D. Iwan. A distributed-element model for hysteresis and its steady-state dynamic re-
sponse. Transactions of the ASME. Series E, Journal of Applied Mechanics, 33(4):893 –
900, 1966.
172
[80] Yi-Kwei Wen. Method for random vibration of hysteretic systems. 102(2):249 – 263, 1976.
[81] W.D. Iwan and L.D. Lutes. Response of the bilinear hysteretic system to stationary random
excitation. Journal of the Acoustical Society of America, 43(3):545 – 552, 1968.
[82] W. D. Iwan. A generalization of the concept of equivalent linearization. International
Journal of Non-Linear Mechanics, 8(3):279 – 287, 1973.
[83] S.F. Masri. Forced vibration of the damped bilinear hysteretic oscillator. Journal of the
Acoustical Society of America, 57(1):106 – 12, 1975.
[84] Wilfred D. Iwan and Arturo O. Cifuentes. Model for system identification of degrading
structures. Earthquake Engineering and Structural Dynamics, 14(6):877 – 890, 1986.
[85] B.F. Spencer and L.A. Bergman. On the reliability of a simple hysteretic system. Journal
of Engineering Mechanics, 111(12):1502 – 1514, 1985.
[86] O. Vinogradov and I. Pivovarov. Vibrations of a system with non-linear hysteresis. Journal
of Sound and Vibration, 111(1):145 – 52, 1986.
[87] C.-Y . Peng and W.D. Iwan. Identification methodology for a class of hysteretic structures.
Earthquake Engineering and Structural Dynamics, 21(8):695 – 712, 1992.
[88] Mohammad Yar and Joseph K. Hammond. Modeling and response of bilinear hysteretic
systems. Journal of Engineering Mechanics, 113(7):1000 – 1013, 1987.
[89] M. Yar and J.K. Hammond. Parameter estimation for hysteretic systems. Journal of Sound
and Vibration, 117(1):161 – 72, 1987.
[90] J.B. Roberts and P.D. Spanos. Random Vibrations and Statistical Linearization. Wiley,
New York, 1990.
[91] S.F. Masri, R.K. Miller, M.-I. Traina, and T.K. Caughey. Development of bearing friction
models from experimental measurements. Journal of Sound and Vibration, 148(3):455 –
75, 1991.
[92] Chin-Hsiung Loh and Sheng-Tsai Chung. Three-stage identification approach for hys-
teretic systems. Earthquake Engineering and Structural Dynamics, 22(2):129 – 150, 1993.
[93] A.G. Chassiakos, S.F. Masri, A. Smyth, and J.C. Anderson. Adaptive methods for identifi-
cation of hysteretic structures. volume 3, pages 2349 – 53, Evanston, IL, USA, 1995.
[94] Koichiro Asano and Wilfred D Iwan. An alternative approach to the random response of
bilinear hysteretic systems. Earthquake engineering & structural dynamics, 12(2):229–
236, 1984.
[95] R Pratap, S Mukherjee, and FC Moon. Dynamic behavior of a bilinear hysteretic elasto-
plastic oscillator, part i: free oscillations. Journal of Sound and Vibration, 172(3):321–337,
1994.
173
[96] R Pratap, S Mukherjee, and FC Moon. Dynamic behavior of a bilinear hysteretic elasto-
plastic oscillator, part ii: oscillations under periodic impulse forcing. Journal of Sound and
Vibration, 172(3):339–358, 1994.
[97] Masaru Hoshiya and Etsuro Saito. Structural identification by extended kalman filter. Jour-
nal of Engineering Mechanics, 110(12):1757–1770, 1984.
[98] Mohammad Yar and Joseph K Hammond. Modeling and response of bilinear hysteretic
systems. Journal of engineering mechanics, 113(7):1000–1013, 1987.
[99] R. Bouc. Forced vibration of mechanical systems with hysteresis. Prague, Czechoslovakia,
1967.
[100] Sashi K. Kunnath, John B. Mander, and Lee Fang. Parameter identification for degrading
and pinched hysteretic structural concrete systems. Engineering Structures, 19(3):224 –
232, 1997.
[101] A.W. Smyth, S.F. Masri, E.B. Kosmatopoulos, A.G. Chassiakos, and T.K. Caughey. Devel-
opment of adaptive modeling techniques for non-linear hysteretic systems. International
Journal of Non-Linear Mechanics, 37(8):1435 – 51, 2002.
[102] R.H. Sues, S.T. Mau, and Y .-K. Wen. Systems identification of degrading hysteretic restor-
ing forces. Journal of Engineering Mechanics, 114(5):833 – 846, 1988.
[103] Haochuan Zhang, Greg C. Foliente, Yongmin Yang, and M. Fai. Seismic response of
self-centring hysteretic sdof systems. Earthquake Engineering and Structural Dynamics,
31(5):1131 – 1150, 2002.
[104] Y .Q. Ni, J.M. Ko, and C.W. Wong. Identification of non-linear hysteretic isolators from
periodic vibration tests. Journal of Sound and Vibration, 217(4):737 – 56, 1998.
[105] Jeen-Shang Lin and Yigong Zhang. Nonlinear structural identification using extended
kalman filter. Computers and Structures, 52(4):757 – 64, 1994.
[106] Meiliang Wu and A.W. Smyth. Application of the unscented kalman filter for real-time
nonlinear structural system identification. Structural Control and Health Monitoring,
14(7):971 – 90, 2007.
[107] E.N. Chatzi and A.W. Smyth. The unscented kalman filter and particle filter methods
for nonlinear structural system identification with non-collocated heterogeneous sensing.
Structural Control and Health Monitoring, 16(1):99 – 123, 2009.
[108] Jih-Lian Ha, Ying-Shieh Kung, Rong-Fong Fung, and Shao-Chien Hsien. A comparison
of fitness functions for the identification of a piezoelectric hysteretic actuator based on the
real-coded genetic algorithm. Sensors and Actuators, A: Physical, 132(2):643 – 650, 2006.
174
[109] N.M. Kwok, Q.P. Ha, M.T. Nguyen, J. Li, and B. Samali. Bouc-wen model parameter
identification for a mr fluid damper using computationally efficient ga. ISA Transactions,
46(2):167 – 179, 2007.
[110] A. Kyprianou, K. Worden, and M. Panet. Identification of hysteretic systems using the
differential evolution algorithm. Journal of Sound and Vibration, 248(2):289 – 314, 2001.
[111] E. B. Kosmatopoulos, A. W. Smyth, S. F. Masri, and A. G. Chassiakos. Robust adaptive
neural estimation of restoring forces in nonlinear structures. Journal of Applied Mechanics,
68(6):880–893, 2001.
[112] R. Fletcher. Practical methods of optimization, Volume 1. Wiley, 1987.
[113] J. A. Snyman. Practical mathematical optimization: an introduction to basic optimization
theory and classical and new gradient-based algorithms. Springer-Verlag, 2005.
[114] Christian Blum and Andrea Roli. Metaheuristics in combinatorial optimization: Overview
and conceptual comparison. ACM Computing Surveys, 35(3):268 – 308, 2003. Combina-
torial optimization;Metaheuristics;.
[115] J. Antonio Parejo, Antonio Ruiz-Cortes, Sebastian Lozano, and Pablo Fernandez. Meta-
heuristic optimization frameworks: A survey and benchmarking. Soft Computing,
16(3):527 – 561, 2012.
[116] D. R. Kincaid and E. W. Cheney. Numerical Analysis: Mathematics of Scientific Comput-
ing. Brooks/Cole, Pacific Grove, CA, 1991.
[117] R. Pratap, S. Mukherjee, and F.C. Moon. Dynamic behavior of a bilinear hysteretic elasto-
plastic oscillator, part i: free oscillations. Journal of Sound and Vibration, 172(3):321 –
337, 1994.
[118] W.D. Iwan. The dynamic response of the one degree of freedom bilinear hysteretic system.
3WCEE, New Zealand, 1965.
[119] Rolf Sobottke, Klaus Schl¨ uter-Brust, Thomas Kaulhausen, Marc R¨ ollinghoff, Britta
Joswig, Hartmut St¨ utzer, Peer Eysel, Patrick Simons, and Johannes Kuchta. Interspinous
implants (x stop R
, wallis R
, diam R
) for the treatment of lss: is there a correlation between
radiological parameters and clinical outcome? European Spine Journal, 18(10):1494–
1503, 2009.
[120] James F Zucherman, Ken Y Hsu, Charles A Hartjen, Thomas F Mehalic, Dante A Im-
plicito, Michael J Martin, Donald R Johnson, Grant A Skidmore, Paul P Vessa, James W
Dwyer, et al. A multicenter, prospective, randomized trial evaluating the x stop interspinous
process decompression system for the treatment of neurogenic intermittent claudication:
two-year follow-up results. Spine, 30(12):1351–1358, 2005.
175
[121] JF Zucherman, KY Hsu, CA Hartjen, TF Mehalic, DA Implicito, MJ Martin, DR John-
son II, GA Skidmore, PP Vessa, JW Dwyer, et al. A prospective randomized multi-center
study for the treatment of lumbar spinal stenosis with the x stop interspinous implant: 1-
year results. European Spine Journal, 13(1):22–31, 2004.
[122] Ken Y Hsu, James F Zucherman, Charles A Hartjen, Thomas F Mehalic, Dante A Im-
plicito, Michael J Martin, Donald R Johnson, Grant A Skidmore, Paul P Vessa, James W
Dwyer, et al. Quality of life of lumbar stenosis-treated patients in whom the x stop inter-
spinous device was implanted. Journal of Neurosurgery: Spine, 5(6):500–507, 2006.
[123] Johannes Kuchta, Rolf Sobottke, Peer Eysel, and Patrick Simons. Two-year results of
interspinous spacer (x-stop) implantation in 175 patients with neurologic intermittent clau-
dication due to lumbar spinal stenosis. European Spine Journal, 18(6):823–829, 2009.
[124] Olaf J Verhoof, Johannes L Bron, Frits H Wapstra, and Barend J van Royen. High failure
rate of the interspinous distraction device (x-stop) for the treatment of lumbar spinal steno-
sis caused by degenerative spondylolisthesis. European Spine Journal, 17(2):188–192,
2008.
[125] Giuseppe MV Barbagallo, Giuseppe Olindo, Leonardo Corbino, and Vincenzo Albanese.
Analysis of complications in patients treated with the x-stop interspinous process decom-
pression system: proposal for a novel anatomic scoring system for patient selection and
review of the literature. Neurosurgery, 65(1):111–120, 2009.
[126] Manal Siddiqui, Efthimios Karadimas, Malcolm Nicol, Francis W Smith, and Douglas
Wardlaw. Effects of x-stop device on sagittal lumbar spine kinematics in spinal stenosis.
Journal of spinal disorders & techniques, 19(5):328–333, 2006.
[127] Vikram Talwar, Derek P Lindsey, Amy Fredrick, Ken Y Hsu, James F Zucherman, and
Scott A Yerby. Insertion loads of the x stop interspinous process distraction system de-
signed to treat neurogenic intermittent claudication. European Spine Journal, 15(6):908–
912, 2006.
[128] H.J. Cramer, Y .King Liu, and D.U. von Rosenberg. A distributed parameter model of the
inertially loaded human spine. Journal of Biomechanics, 9(3):115 – 130, 1976.
[129] AUGUSTUS A WHITE III and MANOHAR M PANJABI. The basic kinematics of the
human spine: a review of past and current knowledge. Spine, 3(1):12–20, 1978.
[130] Gary Monheit and Norman I Badler. A kinematic model of the human spine and torso.
Technical Reports (CIS), page 746, 1990.
[131] Manohar M Panjabi. Three-dimensional mathematical model of the human spine structure.
Journal of Biomechanics, 6(6):671–680, 1973.
176
[132] T.B. Belytschko, T.P. Andriacchi, A.B. Schultz, and J.O. Galante. Analog studies of forces
in the human spine: Computational techniques. Journal of Biomechanics, 6(4):361 – 371,
1973.
[133] Marek Dietrich, Krzysztof Kedzior, and Tomasz Zagrajek. Modeling of muscle action and
stability of the human spine. In Multiple Muscle Systems, pages 451–460. Springer, 1990.
[134] Roy L. Adler, Jean-Pierre Dedieu, Joseph Y . Margulies, Marco Martens, and Mike Shub.
Newton’s method on riemannian manifolds and a geometric model for the human spine.
IMA Journal of Numerical Analysis, 22(3):359–390, 2002.
[135] C.T. Terry and Verne L. Roberts. A viscoelastic model of the human spine subjected to +gz
accelerations. Journal of Biomechanics, 1(2):161 – 168, 1968.
[136] Marko de Jager. Mathematical modelling of the human cervical spine: a survey of the
literature. WFW-93.027, 1993.
[137] CE Aubin, JL Descrimes, J Dansereau, W Skalli, F Lavaste, H Labelle, et al. [geometrical
modeling of the spine and the thorax for the biomechanical analysis of scoliotic deformities
using the finite element method]. In Annales de chirurgie, volume 49, pages 749–761,
1994.
[138] Narayan Yoganandan, Srirangam Kumaresan, Liming V oo, and Frank A Pintar. Finite
element applications in human cervical spine modeling. Spine, 21(15):1824–1834, 1996.
[139] Wayne Z Kong and Vijay K Goel. Ability of the finite element models to predict response
of the human spine to sinusoidal vertical vibration. Spine, 28(17):1961–1967, 2003.
[140] Michael AK Liebschner, David L Kopperdahl, William S Rosenberg, and Tony M Keaveny.
Finite element modeling of the human thoracolumbar spine. Spine, 28(6):559–565, 2003.
[141] Anne Polikeit, Lutz Peter Nolte, and Stephen J Ferguson. The effect of cement augmenta-
tion on the load transfer in an osteoporotic functional spinal unit: finite-element analysis.
Spine, 28(10):991–996, 2003.
[142] E Dall’Ara, R Schmidt, D Pahr, P Varga, Y Chevalier, J Patsch, F Kainberger, and P Zysset.
A nonlinear finite element model validation study based on a novel experimental technique
for inducing anterior wedge-shape fractures in human vertebral bodies¡ i¿ in vitro¡/i¿. Jour-
nal of biomechanics, 43(12):2374–2380, 2010.
[143] Li-Xin Guo, Yi-Min Zhang, and Ming Zhang. Finite element modeling and modal analysis
of the human spine vibration configuration. Biomedical Engineering, IEEE Transactions
on, 58(10):2987–2990, 2011.
[144] Jeng-Wen Lin, Raimondo Betti, Andrew W. Smyth, and Richard W. Longman. On-line
identification of non-linear hysteretic structural systems using a variable trace approach.
Earthquake Engineering and Structural Dynamics, 30(9):1279 – 1303, 2001.
177
[145] Zhishen Wu, Bin Xu, and Koichi Yokoyama. Decentralized parametric damage detec-
tion based on neural networks. Computer-Aided Civil and Infrastructure Engineering,
17(3):175–184, 2002.
[146] Jann N Yang and Silian Lin. On-line identification of non-linear hysteretic structures
using an adaptive tracking technique. International Journal of Non-Linear Mechanics,
39(9):1481–1491, 2004.
[147] SJ Li, H Yu, and Y Suzuki. Identification of non-linear hysteretic systems with slip. Com-
puters & structures, 82(2):157–165, 2004.
[148] Fayc ¸al Ikhouane and Jos´ e Rodellar. On the hysteretic bouc–wen model. Nonlinear Dy-
namics, 42(1):79–95, 2005.
[149] Seyed Ali Ashrafi, Andrew W. Smyth, and Raimondo Betti. A parametric identification
scheme for non-deteriorating and deteriorating non-linear hysteretic behaviour. Structural
Control and Health Monitoring, 13(1):108–131, 2006.
[150] Jann N Yang, Hongwei Huang, and Silian Lin. Sequential non-linear least-square estima-
tion for damage identification of structures. International Journal of Non-linear mechanics,
41(1):124–140, 2006.
[151] Fayc ¸al Ikhouane and Oriol Gomis-Bellmunt. A limit cycle approach for the parametric
identification of hysteretic systems. Systems & Control Letters, 57(8):663–669, 2008.
[152] Eleni N Chatzi, Andrew W Smyth, and Sami F Masri. Experimental application of on-line
parametric identification for nonlinear hysteretic systems with model uncertainty. Struc-
tural Safety, 32(5):326–337, 2010.
[153] Khaldoon Bani-Hani and Jamshid Ghaboussi. Nonlinear structural control using neural
networks. Journal of engineering mechanics, 124(3):319–327, 1998.
[154] Jyh-Da Wei and Chuen-Tsai Sun. Constructing hysteretic memory in neural networks.
Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on, 30(4):601–
609, 2000.
[155] Chih-Chieh Huang and Chin-Hsiung Loh. Nonlinear identification of dynamic systems
using neural networks. Computer-Aided Civil and Infrastructure Engineering, 16(1):28–
41, 2001.
[156] CY Kao and Shih-Lin Hung. Detection of structural damage via free vibration re-
sponses generated by approximating artificial neural networks. Computers & structures,
81(28):2631–2644, 2003.
[157] J-S Pei, AW Smyth, and EB Kosmatopoulos. Analysis and modification of volterra/wiener
neural networks for the adaptive identification of non-linear hysteretic dynamic systems.
Journal of sound and vibration, 275(3):693–718, 2004.
178
[158] Jin-Song Pei and Andrew W Smyth. New approach to designing multilayer feedforward
neural network architecture for modeling nonlinear restoring forces. i: Formulation. Jour-
nal of engineering mechanics, 132(12):1290–1300, 2006.
[159] Jin-Song Pei and Andrew W Smyth. New approach to designing multilayer feedforward
neural network architecture for modeling nonlinear restoring forces. ii: Applications. Jour-
nal of engineering mechanics, 132(12):1301–1312, 2006.
[160] Zhao Tong, Yonghong Tan, and Xianwen Zeng. Modeling hysteresis using hybrid method
of continuous transformation and neural networks. Sensors and Actuators A: Physical,
119(1):254–262, 2005.
[161] Ruili Dong, Yonghong Tan, Hui Chen, and Yangqiu Xie. A neural networks based model
for rate-dependent hysteresis for piezoceramic actuators. Sensors and Actuators A: Physi-
cal, 143(2):370–376, 2008.
[162] Xinlong Zhao and Yonghong Tan. Modeling hysteresis and its inverse model using neu-
ral networks based on expanded input space method. Control Systems Technology, IEEE
Transactions on, 16(3):484–490, 2008.
[163] F Badrakhan. Dynamic analysis of yielding and hysteretic systems by polynomial approx-
imation. Journal of sound and vibration, 125(1):23–42, 1988.
[164] YQ Ni, JM Ko, and CW Wong. Nonparametric identification of nonlinear hysteretic sys-
tems. Journal of Engineering mechanics, 125(2):206–215, 1999.
[165] S.F. Masri, F. Tasbihgoo, J.P. Caffrey, A.W. Smyth, and A.G. Chassiakos. Data-based
model-free representation of complex hysteretic mdof systems. Structural Control and
Health Monitoring, 13(1):365 – 387, 2006.
[166] Sami F Masri, Roger Ghanem, Felipe Arrate, and John Caffrey. Stochastic nonparametric
models of uncertain hysteretic oscillators. AIAA journal, 44(10):2319–2330, 2006.
[167] Farzad Tasbihgoo, John P. Caffrey, and Sami F. Masri. Development of data-based model-
free representation of non-conservative dissipative systems. International Journal of Non-
Linear Mechanics, 42(1):99 – 117, 2007.
[168] Amir Hossein Alavi, Amir Hossein Gandomi, Mohammad Ghasem Sahab, and Mostafa
Gandomi. Multi expression programming: a new approach to formulation of soil classifi-
cation. Engineering with Computers, 26(2):111–118, 2010.
[169] Amir Hossein Gandomi and Amir Hossein Alavi. Multi-stage genetic programming: a new
strategy to nonlinear system modeling. Information Sciences, 181(23):5227–5239, 2011.
[170] S.F. Masri, R.K. Miller, A.F. Saud, and T.K. Caughey. Identification of nonlinear vibrat-
ing structures. ii. applications. Transactions of the ASME. Journal of Applied Mechanics,
54(4):923 – 50, 1987.
179
[171] SF Masri, JP Caffrey, TK Caughey, AW Smyth, and AG Chassiakos. A general data-based
approach for developing reduced-order models of nonlinear mdof systems. Nonlinear Dy-
namics, 39(1-2):95–112, 2005.
180

Abstract (if available)

Abstract This study builds on major advances in the field of Computational Intelligence to develop a state‐of‐the‐art data‐driven methodology that provides parsimonious optimized computational models in the form of systems of differential equations that characterize the behavior of complex nonlinear phenomena observed in mechanical and biological systems. The proposed hybrid identification scheme integrates various stochastic optimization methods and computer algebra techniques, such as Genetic Programming and Genetic Algorithms, to evolve structures of differential equations, to optimize their parameters, and to reduce their complexity for controlling bloat. The investigated scenarios include systems that exhibit polynomial‐type nonlinearities in their response, systems that show discontinuity in their nonlinear behavior, systems with memory‐dependent and dissipative characteristics, as well as the human spine. The investigations are conducted by processing input and output data obtained from synthetic simulations as well as experiments. It is shown that the proposed technique yields reduced‐order, reduced‐complexity, optimized differential equations, that accurately characterize the behavior of the investigated systems, and provide accurate estimates. The generalization extent of the discovered models is scrutinized by assessing their performance in new dynamical environments through applying validation excitations that are substantially different from the excitations employed for training. Findings reveal that the resulting models provide reasonably accurate estimates, even when models are subjected to new stimulations with various intensities. Thus, the proposed approach of this study presents a robust data‐driven methodology based on evolutionary computation techniques that provides elegant computational models to represent variety of complex nonlinear systems.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Studies into data-driven approaches for nonlinear system identification, condition assessment, and health monitoring

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Analytical and experimental studies in the development of reduced-order computational models for nonlinear systems

PDF

Characterization of environmental variability in identified dynamic properties of a soil-foundation-structure system

PDF

Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures

PDF

Analytical and experimental studies in modeling and monitoring of uncertain nonlinear systems using data-driven reduced‐order models

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Structural nonlinear control strategies to provide life safety and serviceability

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Analytical and experimental studies of modeling and monitoring uncertain nonlinear systems

PDF

Optimal clipped linear strategies for controllable damping

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Theoretical and computational foundations for cyber‐physical systems design

PDF

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

PDF

On the synthesis of dynamics and control for complex multibody systems

PDF

Vision-based studies for structural health monitoring and condition assesment

PDF

New approaches in modeling and control of dynamical systems

PDF

On the synthesis of controls for general nonlinear constrained mechanical systems

Asset Metadata

Creator Bolourchi Yazdi, Seyed Ali (author)

Core Title Studies into computational intelligence approaches for the identification of complex nonlinear systems

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering (Structural Engineering)

Publication Date 04/26/2016

Defense Date 11/13/2013

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag computational intelligence,computer algebra,differential equations,dimensionality reduction,evolutionary computation,experimental data,expression trees,genetic algorithms,genetic programming,human spine,hysteretic systems,intelligent selection,large volume of data,nonlinear systems,OAI-PMH Harvest,parsimonious models,system identification

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Masri, Sami F. (committee chair), Ghanem, Roger G. (committee member), Mendel, Jerry M. (committee member), Wellford, L. Carter (committee member)

Creator Email bolourch@gmail.com,bolourch@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-386427

Unique identifier UC11296124

Identifier etd-BolourchiY-2426.pdf (filename),usctheses-c3-386427 (legacy record id)

Legacy Identifier etd-BolourchiY-2426.pdf

Dmrecord 386427

Document Type Dissertation

Format application/pdf (imt)

Rights Bolourchi Yazdi, Seyed Ali

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA