.fltin nth
, . .5 v). I
v ‘ ”19.1.?
A4 1.. . § . O :I
if iv. .2. ‘ .

c...» n ,n.

1 iv.
.. ..
$5.3... .
‘ :32}

.1...)

E s».
, Aug;
.m a .msfifi
. a . ‘ )
4w}.
1 9
$3.“ 1....
L: .. Jam.
. .J.
345w“.
. I} I . u
,, ‘ Wis:
12....-. .5:
.533 $1., 555.
J}. . £512.

:3?

.233;
.3 1
51 471!

, i550”. "
9:35:
.1 L. Luna... 3. . x
I, Z. .
. .. I 3!... r. ..

x2355. 1%,
vmvuh, tngiiifi
k.r.....n...
. .5

 

b. 5 a...
.4.“ fl...

1 .91 {Q‘s

is...

.Zil‘:
.t..1|: :\\

353345
. a. i ..l. .. .. a.
. . L many-fluid. I:

3‘!”

. . w... .2:
. , ..
‘ 4 5.3.1.;
. in. ._ . . am. . a.
. x ,7: V , V . ’u‘l
y. p I V «M?! .i. .
.9313 Pet. , a 1..., .6: .73.». .
. i .. :3...
. . :1, 1,... v1. .vtal...ll..
{Huuuipfié nun; .a . . ‘ 53-111.qu .... .13
‘ , z... .3751: 5‘3. (.1 £15..
IVEIQX 01:? .1 .2 ‘

.

.. an:
3 i... t»
.iZJI. .
. 7.3.11 . .
A. . . ‘ .

.5 . ,u.
33.“ ‘ ..

 

"zl 2!. '4
it) .1 .-
$3..

a; 4......
... .¥.Q.1.I..ufbuiivk
...a .If..- ‘ k ‘ ..

29114

This is to certify that the
dissertation entitled

PARAMETRIC IDENTIFICATION OF
CHAOTIC/NONLINEAR SYSTEMS AND
REDUCED ORDER MODELS BASED ON

PROPER ORTHOGONAL
DECOMPOSITION

presented by

Yang Liang

has been accepted towards fulfillment
of the requirements for the '

 

 

 

 

Ph.D. degree in Mechanical Emeering
Major [Profess73é Signature
l 2 / l / 2 o o f
Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

[ LIBRARY
Michigan State
I University

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

'JAN 1 2 2009i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p'lClRC/DateDue.indd-p.1

 

 

PARAMETRIC IDENTIFICATION OF
CHAOTIC / NONLINEAR SYSTEMS AND
REDUCED ORDER MODELS BASED ON

PROPER ORTHOGONAL DECOMPOSITION

By

Yang Liang

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2005

ABSTRACT

PARAMETRIC IDENTIFICATION OF
CHAOTIC/ NONLINEAR SYSTEMS AND REDUCED
ORDER MODELS BASED ON PROPER ORTHOGONAL
DECOMPOSITION

By
Yang Liang

In the dissertation, the parametric identification and proper orthogonal modes
(POMS) were investigated on chaotic/nonlinear systems such that a reliable and
general-purpose process can be developed to reconstruct the mathematical models
of nonlinear and/or chaotic systems.

First, a parametric identification method was examined for different chaotic sys-
tems, e. g. whirling, multi-degree-of-freedom, and strong nonlinearity, and by sim-
ulation and experiments. The original parametric identification method of chaotic
systems is a hybrid time and frequency domain method based upon the harmonic
balance method applied to unstable periodic orbits (UPOs), and solved by least mean
squares. A chaotic base-excited single pendulum system was simulated. The identi-
fication method was modified for the whirling system with measured data of angular
displacement. The nonlinearity was also approximated by two types of function se—
ries: linear interpolation functions and harmonic functions. Poincare sections showed
that the identified system and the original simulation system had similar chaotic

behavior. Then, the identification method was applied to an experimental chaotic

double pendulum under vertical base excitation. Several noise reduction techniques
were applied to reduce the identification errors due to the noise contamination in
the experimental data and the strong nonlinearity. Meanwhile, an error optimiza-
tion process, based on the linear regression and statistics, was proposed to improve
the identified parameters by selecting sub-harmonics of the unstable periodic orbits.
An energy balance method, as a second-step identification after the harmonic bal-
ance method, was applied to give a more accurate estimation of the small damping
coefficients.

It was also found that any chaotic orbit can be an approximated representation
of some UPO, and the longer the orbit the better the approximation. Thus, there
is an option between collecting UPOs for the extraction, and using longer arbitrary
data segments. The identification process can be simplified to a frequency domain
method. Examples were examined to show the success of the simplification.

The study then goes to parametric identification and the POM for building re-
duced order models of unknown systems. The special interest here is systems with
strong nonlinearity, where the reduced order models from experimental POM was
limited in simulating the unknown systems other than the neighborhood set plane
where the POM data is in the phase space. An added-constraint method was pro—
posed and examined by two simple systems: a two-mass system with nonlinear spring
and a mass-pendulum system. It showed that the added constraint can improve ap-
plicability of the POM reduced order model as well as increasing the accuracy of
the simulation result. Nevertheless, the method is to be tested by high dimensional
nonlinear systems, to which the proposed added constraint method can really make

a big difference in applications.

Copyright © by
Yang Liang
2005

To Xiaobei.

ACKNOWLEDGMENTS

I wish to acknowledge the National Science Foundation (OMS-0099603) and NASA
(N AG-1-01048 with Dr. Walt Silva) for their finaicial support for this research.

I would like to express my gratitude to my advisor Professor Brian Feeny for the
guidance, help and encouragement provided during the course of my research.

I am please to acknowledge my colleagues at the Dynamics Lab who created a
pleasant working atmosphere. I would like to thank particularly Brian Olson for
helpful discussions.

I would also like to acknowledge the help of Professor Steven Shaw on discussions
of nonlinear normal modes, and Professor Michael Frazier on talks of the wavelet
analysis.

Finally, I would like to express my love and gratitude to my wife Xiaobei, my
mom, dad and my brother for their enormous support and understanding during my

study at Michigan State University.

vi

TABLE OF CONTENTS

LIST OF TABLES ix

LIST OF FIGURES x

1 Introduction 1

1.1 Motivation ................................. 1

1.2 Literature review ............................. 4

1.3 Outline of the dissertation ........................ 6

2 Simulation Study of a Base-Excited Single Pendulum 10

2.1 Introduction ................................ 10

2.2 Horizontally excited single pendulum & identification algorithm . . . 11
2.3 Introduction to phase plane reconstruction and extraction of unstable

periodic orbits ............................... 17

2.3.1 Choosing the time delay and embedding dimension ...... 18

2.3.2 Extracting unstable periodic orbits (UPOs) .......... 18

2.4 Numerical simulation of the horizontally excited pendulum ...... 19

2.4.1 Phase space reconstruction .................... 19

2.4.2 Parametric identification ..................... 20

2.5 Discussion and validation ......................... 27

2.6 Concluding remarks ............................ 28

3 Experimental Study of a Base-Excited Double Pendulum 29

3.1 Description of the double pendulum system .............. 29

3.2 Method .................................. 32

3.3 Experiment description .......................... 37

3.4 Results and validation .......................... 39

3.4.1 Phase space reconstruction & UPO extraction ......... 39

3.4.2 Identified parameters ....................... 45

3.4.3 Friction issue ........................... 46

3.4.4 High frequency noise in unstable periodic orbits ........ 48

3.4.5 Digital differentiation and error reduction ........... 54

3.4.6 Validation methods ........................ 62

3.5 The energy balance method for identification of the damping coefficients 65

3.6 Concluding remarks ...... , ...................... 68

vii

4 Simplification to Hequency Domain Method 70

4.1 System identification by frequency method ............... 70
4.2 More chaotic properties and the improved frequency domain method 72
4.3 Calculation examples ........................... 76
4.3.1 Base-excited single pendulum .................. 77
4.3.2 Experimental double pendulum ................. 79
4.4 Concluding remarks ............................ 80

5 The Reduced Order Models of Nonlinear Systems Based on Proper
Orthogonal Modes 85
5.1 Need for reduced order models ...................... 85
5.2 Proper orthogonal decomposition introduction ............. 86
5.3 Improving POD for strong nonlinear systems .............. 90
5.3.1 Limitation of the POD ...................... 90
5.3.2 The example of a two-mass system with nonlinear spring . . . 92
5.3.3 The example of a mass-pendulum system ............ 99
5.4 Concluding remarks ............................ 112
6 Conclusions and Future Work 114
6.1 Main contributions ............................ 114
6.2 Directions for future work ........................ 117

BIBLIOGRAPHY 120

viii

LIST OF TABLES

2.1 The nk, Pk and qk values when harmonic basis is applied. ......

3.1 Physical properties of the double pendulum. ..............
3.2 Experimental settings. ..........................
3.3 Identified parameters by applying low-pass filter, FFT convolution,
Fourier series expansion of UFO, harmonic balance method with sub-
harmonics whose frequencies are S excitation frequency fe=5 Hz, and
sub—harmonics optimization. .......................
3.4 Identified non—dimensional parameters by optimized identification pro-

45

cess and noise reduction provided that the friction coefficients are known. 48

3.5 Comparison of .7:(y) with and without filter added. ..........
3.6 Optimized identification with no digital filter applied ..........
3.7 Comparison of calculated values by Fourier Series (FS) method and
Digital Differentiation (DD) method when error tolerance is set as 8%
with filtering and optimization process applied. ............
3.8 Identified parameters with the Fourier series method and low-pass fil-
tering, but without optimization to sub-harmonics set. ........
3.9 Comparison of the identified values and the true parametric settings of

the simulated system with low-pass filtering and optimization applied.

3.10 Natural frequencies .............................
3.11 Identified damping coefficients by the energy balance method.

4.1 Comparison of the true values and the identified parameters by new
frequency method. ............................

5.1 Identified parameters of the POM-constraint model (5.16) under the

first three conditions ...........................
5.2 Physical parameters used in the simulation ...............
5.3 Estimated polynomial coefficients of the constraint equation .....

ix

52
54
55
60
64

65
67

101

2.1

2.2

2.3
2.4

2.5

2.6

2.7

2.8

3.1
3.2
3.3

3.4

3.5
3.6
3.7
3.8
3.9
3.10

3.11

LIST OF FIGURES

A period-4 orbit (a) sampled signal; (b) actual continuous signal me;
(e) -- — constant rotation part wckt, — oscillatory part $var- ......
Reconstructed phase plane (a) and Poincaré section plot (b); 6(t) is
represented by theta(t) in the plots ....................
Mutual information of the data; the first minimum is at dt=5. . . . .
(a)A period-1 orbit; (b) A period-2 orbit; (c) A period-2 orbit; (d) A
period-2 orbit ................................
The "k curve for linear interpolation basis; (*) identified values, (—)
true function. ...............................
The pk curve for linear interpolation basis; (*) identified values, (—)
true function. ...............................
The qk curve for linear interpolation basis; (0) identified values, (——)
true function. ...............................
Poincaré section plots of identified systems: (a) by linear interpolation
basis, (b) by Harmonic basis. ......................

A sketch of the double pendulum .....................
Mutual information 1(dt). ........................
Phase portrait of experimental data with dt=24, (a) 01(t)—01(t + dt),
(b) 61(t)-02(t), (c) 02(t)—02(t + dt); 0 is represented by theta in the
figure. ...................................
A period 4 UPO; thetal and theta2 in the plot represent 91 and 62;
dt=24; same for the following figures. ..................
A period 9 UPO. .............................
A period 12 UPO ..............................
A period 12 UPO ..............................
A period 15 UPO ..............................
A period 15 UPO ..............................
Phase portrait of the simulated system with C11 = 1.02 x 10-3,
011:0.0366 and Cg = 0.0485; (#:24 ...................
F F T amplitude of the 0.1 and 01 with and without low pass filter ap-
plied; the continuous line is FFT of signals with filter of 1 / 5 f3 cut-off
frequency; the dotted line is FFT of signals without filter; It is the
order of sub-harmonics ...........................

14

20
21

22

24

25

26

27

30
39

40

42
42
43
43
44
44

49

52

3.12

3.13

3.14

3.15

4.1
4.2

4.3

4.4

4.5

5.1

5.2

5.3

5.4
5.5

5.6

5.7

Signals without and with low pass filter applied; (a) and (b) are 0.1
without and with filter; (c) and (d) are y without and with filter applied
to each convolution components; k is the kth sampled point of the orbit. 53
Obtained 91 and 03 orbits by Fourier series method (a) and (b) by dig-
ital differentiation method (c) and (d); the recurrence impulses occur
at k=800 in (a) and (b). ......................... 57
Identification residual 62 with each dot representing a subharmonic; the
horizontal position of each dot is the predicted value of 6'2, the vertical
position of each point is the identification residual; identification was
done with low-pass filtering, the Fourier series expansion and the sub-

harmonic set whose frequency 3 the excitation frequency. ...... 59
Identification residual 52 after optimization with each dot representing
a subharmonic. .............................. 62
The left shift one dimensional map ................... 73
Identified functions by interpolation functions on noise-free data, (a)
nk, (b) pk, (c) qk .............................. 79
Comparison of Poincare sections, (a) the original system, (b)the iden-
tified system ................................ 80
Identified functions by interpolation functions on 10% noise contami-
nated data, (a) nk, (b) pk, (c) qk ..................... 81
The residual plot of the first matrix equation of the double pendulum
system ................................... 82

Phase portraits of the original system (‘+’), the POD only model (the
dashed line), and the POD-constraint model (the dotted line) under

the 5th initial condition. ......................... 92
First NNM displacement shapes in y1, y2 coordinates, number 1—5
correspond to the five initial conditions. ................ 93

Comparison of the invariant manifold: (a) under small amplitude (the
first three conditions), (b) under all five conditions; Surf 1 is the em-
pirical surface from POD, Surf 2 is Shaw’s surface, the continuous lines
are the orbits under five conditions; :1: axis is y'2, y axis is 312, z axis is yl. 97
A sketch of the mass-pendulum system .................. 100
The shapes of the lst N N M of the mass-pendulum system under 5 i. c.
: 1. (0125,0385), 2. (0250,0733), 3. (0.375,1.10), 4. (05001.54),
5. (0575,1892) .............................. 102
Comparison of Model I (the dotted lines) to the original mass-
pendulum system (the continuous lines) simulated under lst initial
condition, (a) phase plane, (b) mode shape. .............. 105
Comparison of Model I (the dotted lines) to the original mass-
pendulum system (the continuous lines) simulated under 3rd initial
condition, (a) phase plane, (b) mode shape. .............. 106

xi

5.8

5.9

5.10

5.11

5.12

5.13

Comparison of Model I (the dotted lines) to the original mass-
pendulum system (the continuous lines) simulated under 5th initial
condition, (a) phase plane, (b) mode shape. ..............
Comparison of Model II (the dotted lines) to the original mass-
pendulum system (the continuous lines) simulated under lst initial
condition, (a) phase plane, (b) mode shape. ..............
Comparison of Model II (the dotted lines) to the original mass-
pendulum system (the continuous lines) simulated under 3rd initial
condition, (a) phase plane, (b) mode shape. ..............
Comparison of Model 11 (the dotted lines) to the original mass-
pendulum system (the continuous lines) simulated under 5th initial
condition, (a) phase plane, (b) mode shape. ..............
A time domain comparison between the original system (continuous
line), the Model-I (dashed line), the Model—II (dotted line) under the
3rd initial condition. ...........................
Natural frequency comparison under different initial conditions (rep-
resented by the maximum potential energy) between the original sys-
tem (continuous line), the Model-I (dashed line), the Model—II (dotted
line); the ‘*’ and ‘0’ mark the five initial conditions for Model-I and
Model-II ...................................

xii

111

CHAPTER 1

Introduction

1.1 Motivation

System identification, as an inverse question of dynamical systems analysis and quite
often a prerequisite of dynamics analysis, has been widely used and investigated in
industry, such as for the control and modeling of complex structures. There are

basically two classifications of system identification:

1. Non-parametric identification requires no actual governing dynamical equations
of a system, but rather an approximated model by a series of approximation
functions. Thus, no physical paremeters are identified, but coefficients of the
approximation functions. Some of the non-parametric identification methods
are based upon differential equations, which share many common features with
the parametric identification that we are going to discuss. On the other hand,
many more are based upon time-delay models, which are to be used in control
applications to approximate short-term behaviors of systems. The time-delay
models are also quite easy for instant, real time approximation of the dynamical

behavior.

2. Parametric identification is based upon governing differential equations, or re-

duced order models of the systems. Therefore, the identified parametric systems
can be better representations of the original unknown systems, and are often

used for dynamics analysis of the systems.

Our focus here is to use parametric identification to reconstruct a model which can
simulate an original unknown system. The advantages of a nonlinear mathematical
model over a linear one are quite obvious: more accurate description of the real phys-
ical systems and thus, more precise prediction of dynamical behavior, which satisfies
the requirements of modern technologies. For the nonlinear systems, chaotic systems
are among the most mysterious. Due to their random-like responses, they were once
considered stochastic systems, until “chaos theory” and chaotic data observation tech-
niques were developed. Chaos is generally considered as aperiodic time-asymptotic
behavior in a deterministic system which exhibits sensitive dependence on initial con-
ditions [1]. The phase-space trajectories of chaotic systems do not settle down to fixed
points or periodic orbits. Many developed identification methods are not suitable for
systems with chaotic responses. Therefore, this study here is aimed at developing
and perfecting an identification process which can be applied to general chaotic (also
nonlinear) systems. The systems in our scope include strongly nonlinear systems,
whirling systems, and high dimensional systems, because of their special features and
difficulty in accurate identification.

A general expression of a class of nonlinear oscillation systems is

Mite) + are) + Km) + fNL(:E, 55) + pr(:E, f)fp(t) = f(t), (1.1)

where f is a vector variable, matices M, K, C are mass, spring and damping matri-
ces, vectors f N L and f N p are nonlinear functions of displacement and velocity, and
fp(t) is the parametric excitation term. For simple systems, we can derive the math-

ematical model, and identify the actual physical parameters. However, for many high

dimensional systems, such as fluids and continuous structures, we tend to avoid the
exact differential equation of the original system, but use a simplified reduced order

model based upon modal analysis. The benefits of reduced order models (ROMS) are:
1. ROMS are simple and require much fewer parameters to identify.

2. Generally, only a few of the low order normal modes of that system show up in

the dynamics of an actual system.

3. Identification of partial differential equations (PDES) of continuous systems
is much more complex than ordinary differential equations (ODES) and also

restricted by experimental data.
4. The ODES are handier in simulation than the PDE, and thus, easier for analysis.

We can simplify nonlinear continuous or high dimensional systems by using either
assumed modes, linear normal modes (LN MS), nonlinear normal modes (N N Ms) [2]
or proper orthogonal modes (POMS) [3] to project the PDES into the form of a series
of ordinary differential equations in modal coordinates. Then, Since usually only few
of the normal modes (or orthogonal modes if POMS are used) are typically active
with some level of Significance in the dynamical behavior of a system, a truncated
series of equations can be used to approximate the original system.

For linear systems, both NNMS and the weighted POMS [4] become the linear
normal modes (LNMs). But, it Should be noticed that the NNM method is more
fit to theoretical analysis which assumes known systems’ structures, and governing
equations. On the other hand, POM method is more of an experimental method,
which deals with orthogonal modes obtained from experimental or testing data, by
which the governing differential equations can be reconstructed.

Hence, the scheme of the dissertation involves two parts: parametric identification

and reconstruction of reduced order models by POD.

1 .2 Literature review

Till now, most of the methods of parametric identification of nonlinear systems fo-
cused mainly on free vibration, random excitation or periodically forced steady state
vibration behavior. Mohammad [5] introduced a direct identification method by
using time-domain displacement, velocity and acceleration Signals. The method is
Simple, but also requires all of the displacement, velocity and acceleration informa-
tion, and is not noise-resistant. In [6], nonlinear resonances by random excitations
were utilized and part of the parameters were identified. Chen and Tomlinson [7]
proposed a time series acceleration, velocity and displacement model (AVD model)
with a narrowed investigation scope: dry Coulomb friction, viscous damping and
nonlinear stiffness. Also, many other methods are effective within limited dynami-
cal systems [8, 9, 10, 11, 12, 13, 14, 15]: weak non-linearity, complex algorithm and
time-consuming, non-chaotic behavior, Single or two degree-of—freedom systems.
Some of the identification methods, though not originally proposed or examined
under chaotic systems, have potential application to be applied to chaotic nonlinear
systems. Plakhtienko [16, 17] introduced a special weight method for parametric iden-
tification, which is derived from orthogonal basis functions. If Fourier series functions
are applied as the orthonormal basis, the method then becomes the harmonic balance
method. The inverse harmonic balance method has been applied to various nonlin-
ear systems identification [18, 19, 20, 21] which are under forced excitation and have
steady-state periodic responses. Thothadri [22] also developed a similar harmonic
balance identification on multi-degree—of-freedom nonlinear systems. A wavelet-based
approach was discussed by Ghanem and Romeo [23]. Both the discrete wavelet trans-
form and Fourier transform are based upon Signal reconstruction by an orthonormal
basis, except that the wavelet transform is a localized time-frequency reconstruction.

Meanwhile, the harmonic balance method is, in essence, a frequency domain method,

and much easier to apply than the wavelet method. Due to this property, Feeny and
Yuan [24, 25] proposed a method for chaotic systems, which exploits the harmonic
balance of extracted unstable periodic orbits (UPOS), because one of the fundamen-
tal properties of deterministic chaos is that the chaotic set of a dynamical system
contains infinite number of unstable periodic orbits. It is the primary method used

in this the dissertation. There are several reasons for choosing this method:

a The identification process can be theoretically applied to any dynamical sys-
tem. For linear systems, it becomes the frequency response function method if
multiple excitation frequencies are tested. For non-chaotic nonlinear systems,

there have been successful applications of the harmonic balance method.

0 The method exploits the existence of unstable periodic orbits, and has generated

satisfying identification results.

0 Some possible developments exist for further improvement, including the friction

issue, experiments, and error reduction.

Parametric identification of large order systems is greatly facilitated in reduced or-
der models. A popular tool for reduced order modeling is proper orthogonal decompo-
sition. Proper orthogonal decomposition (POD), or Karhunen-Loeve decomposition,
has been widely used for empirical modal analysis. As a statistical method utilizing
the correlation matrix derived from a set of measurement histories, POD produces
the empirical modes, which are known as proper orthogonal modes (POMS), from
the eigenvectors of the correlation matrix, and the corresponding energy levels of the
modes, known as proper orthogonal values (POVS), from the matrix eigenvalues.

For linear systems with evenly distributed mass, POMS are actually the linear
normal modes [26]. Therefore, it can be used as modal analysis tool [27, 28, 29, 4, 30,

31]. Meanwhile, since its introduction, POD has been widely used in fluid, turbulence

areas as modal reduction tool [3, 32, 33], and recently in structural dynamics [34, 35,
36, 37, 38, 39, 40, 41, 42, 43]. Some used it to estimate the active dimension of a
dynamical system [33, 37, 38].

For reduced order models, combining POD with system identification is a very
promising area in the analysis of systems dynamics. For chaotic and nonlinear sys-
tems, the POM based reduced order models can approximate the systems’ behavior
by fewer modes than the models based upon linear normal modes [41], which implies
that the POM models are more efficient than the LN M models for reducing nonlinear
systems. POD’s applications in system identification and reduced order modeling
for nonlinear/chaotic systems has also gained some improvements in applications of
systems with Simple nonlinearity [39, 44, 40, 42, 43]. In [39] the POM itself is also

part of the system identification process.

1.3 Outline of the dissertation

The purpose of this dissertation is to generate accurate and robust models of non—
linear and/or chaotic vibration systems. We can either do it by direct parametric
identification provided that the governing equations of the systems are handy, or
for large-order systems, by obtaining approximated reduced order models based upon
POMS or N N Ms of the systems, and identifying parameters of the reduced order mod-
els. Consequently, the dissertation consists of two major parts. The first part focuses
on parametric identification of single and two degree-of-freedom (d. o. f.) systems.
The systems studied are whirling pendulum systems. The second part is on the para-
metric identification of ROMS. POD is used for the reduced-order modeling, and we
focused attention on broadening the applicability of identified models to conditions
beyond those for which the POD was performed. Details are outlined below.

First, we will concentrate the focus on parametric identification method, e. g.

a harmonic balance based identification process for chaotic and nonlinear systems
[24]. Since the method has been tested on Simple single (1. o. f. systems, the goal of
this study is to test it on systems with more complicated nonlinearity and multiple
d. o. f. systems. The pendulums, due to their whirling nature and trigonometric
nonlinearity, are distinctive from ordinary oscillatory systems, and thus are chosen
for the investigation. The pendulums are also easily examined by experiments due
to their simple structures. A Simulated single pendulum was first examined. Then, a
base-excited double pendulum was selected for experimental application to whirling
and multi—degree-of-freedom systems. Once successfully tested, the follow-up goal is
to Simplify and improve the current method. Hence, the identification investigation
will be divided into three separate chapters, and each with a different focus.
Chapter 2 deals with a simulation study of a chaotic single pendulum system. The
main contribution of this chapter is to apply the identification method to whirling
chaotic systems. A brief introduction will be given to explain special features, dynam-
ical behaviors and tools for analysis of chaotic systems. The purpose is to examine
the method on a whirling, chaotic and parametrically excited single degree-of-freedom
system. It is assumed that the system is unknown to us except that the pendulum
is parametrically excited with some known frequency. Meanwhile, to make the simu-
lation resemble an experiment setup, only the angular displacement signal was used.
An embedding technique was then applied to obtain the reconstructed phase space,
which is essential to extract the UPOS from the chaotic data. The extracted orbits
are actually not the true unstable periodic orbits, but their approximations, whose
accuracy can be measured by the tolerance error during the extraction. It is also
hoped that many UPOS can be obtained so as to increase the reliability of the iden-
tified results. Based upon the UPO, the harmonic balance method can be applied
to build the identification matrix equation, and the parameters can be identified si-

multaneously. Another aspect of the identification study is the use of function series

to approximate the unknown nonlinear functions. Two sets of function series were
used: the Fourier series and the linear interpolation functions. Both approximations
Showed good matches with the true values.

Chapter 3 upgrades the study to a experimental double pendulum system. This in-
vestigation was to examine the identification method by a experimental multi-degree—
of—freedom system, to develop noise reduction, optimization techniques to improve
the identification results, and to estimate damping terms more accurately by energy
balance method. Again, the system is chaotic, whirling and parametrically excited.
Since the collected angular displacements were noise contaminated, noise reduction
techniques were applied. The strong nonlinearity of the double pendulum was also
found to contribute to the errors in the identified parameters, especially those non-
linear terms with velocity and/ or acceleration Signals. To quantify and reduce error
in the identification results, an identification error was defined to indicate quantita-
tively how ‘good’ the identification is, and based upon it, an optimization process
was developed. The core of the Optimization process is to obtain a smart choice of
sub-harmonics of the UPOS for harmonic balancing. Based upon the results of the
harmonic balance method, the energy method was developed as a follow-up to give
more accurate estimation of the damping parameters.

Chapter 4 describes a Simplified identification process. The underlying reason
for the improvement is because UPOS become less available for extraction, as the
degrees of freedom of a system goes up. As a result, the identification result can
deteriorate with fewer UPOS. The idea comes from the unstable periodic orbits inside
the chaotic data. Since the UPOS are the ‘Skeleton’ of the chaotic set, every chaotic
orbit is considered to be an approximated periodic orbit of large period. Thus, the
extraction step can be skipped, making the whole identification process look quite like
a frequency domain method for linear systems. The importance of the Simplification

can not be underestimated. It actually makes it possible to identify all kinds of

systems in the frequency domain.

Ultimately, as the d. o. f. of systems increases, the need for reduced-order mod-
eling becomes mandatory. Thus, the second part of the dissertation is to combine
POD and parametric identification to build reduced order models for nonlinear sys-
tems. For high dimensional systems, reduced order models are effective replacements
of the original governing equations. POD provides a simple and theoretically appli-
cable vehicle for reduced order modeling through experimental data. Our focus is
on systems with strong nonlinearity, where the POM based models are limited to a
small neighborhood of the initial conditions where the POM comes from due to the
nonlinearity. A new added-constraint method was developed and tested to improve
the reconstructed reduced order models. Due to the complexity of this topic, the
study was only a beginning of a series of future studies.

Chapter 5 poses the problem of reduced—order modeling with strongly nonlinear
systems. By borrowing ideas from N NM analysis, a geometric constraint method was
proposed to accompany the POM models. The idea was tested on two simple systems:
a two-mass system with nonlinear Spring and a mass-pendulum system. The analysis
was based upon simulation. It turned out that the added constraint method can not
only increase the accuracy of the POM models, but extend the applicability of those
models. Although it is not clear if this method will be as effective on high (1. o. f.
systems as it is to the analyzed low (1. o. f. examples, it gives us a possible direction
to improve the accuracy of the reduced order models.

Finally, Chapter 6 will conclude the dissertation with possible future research on
this topic. The improvements could be real-time identification of chaotic systems,
further identification of small parameters (i. e. frictions terms), POD for whirling
systems, and experiments or simulations of reduced order models of nonlinear high

dimensional systems.

CHAPTER 2

Simulation Study of a

Base-Excited Single Pendulum

2. 1 Introduction

A Simulated, horizontally base-excited pendulum is investigated as a first step for
whirling, chaotic system identification. A harmonic balance parametric identification
method is examined here, because of its Simplicity and capacity for handling chaotic
data.

Pendulum systems are among the most thoroughly investigated dynamic systems
in chaotic, nonlinear system research [45, 46, 47, 48, 49, 50, 51] for their simplicity in
both theoretical expression and experimental validation. Water [45] studied the un-
stable periodic orbits of a vertically excited system. Jeong and Kim [46] investigated
the bifurcation phenomena and routes to chaos of a horizontally excited system. Fur-
thermore, Bishop [47] and Dooren [48] studied the regions of chaos of a parametrically
excited pendulum in parametric space. These works provide the general relationship
between chaotic behavior and the influence of parameters on the system, and facilitate
our investigation of the pendulum with horizontal base-excitation.

Meanwhile, based upon the harmonic balance method [18, 19, 20, 21], Feeny and

10

Yuan [24, 25] proposed a method for chaotic systems which exploits the extracted
unstable periodic orbits. This method was successfully applied to single degree of
freedom systems, and is theoretically applicable to chaotic or non-chaotic, strongly
nonlinear multi-degree-of-freedom systems.

Our purpose here is to apply the harmonic balance algorithm to the horizontally
excited pendulum system such that we can examine the algorithm’s applicability on
systems with chaotic whirling behavior and parametric excitation. The identification
method has been examined on smooth excitation single (1. o. f. systems [24]. However,
the pendulum is a more complicated case. Since the single pendulum is both whirling
and under parametric excitation, the algorithm must be modified and can now be
improved partially. In the following sections, the Single pendulum system and the
modified method will be introduced first. Pre-requisites for applying this method
will be discussed, primarily regarding the phase plane reconstruction. Then, the

Simulation results will be presented and discussed.

2.2 Horizontally excited single pendulum 8:: iden-
tification algorithm

Single pendulum systems have simple structures, but strong non-linearity due to
their whirling property. The governing differential equations can be simulated eas-
ily, and the experimental verification is also feasible. Based upon these advantages,
horizontally base excited single pendulums are chosen for investigation here. The

non-dimensional form of the differential equation is

é+2g/ré+1/r2sin0—fsintcos0=0, (2.1)

11

where t = wr, w is the angular excitation velocity, and 'r is the actual time; 5 is
the viscous damping coefficient; r = w/wN, w N = mge/ J is the natural frequency
of the linearized system, m is the pendulum mass, 6 is the the centroid offset from
the pendulum hinge, J is the moment of inertia about the joint point. Meanwhile
coefficient f = mea/ J is the non-dimensional excitation amplitude, and a is the
excitation amplitude. For simplicity, we denote c1- : 2E/r as the new non-dimensional
damping coefficient. The function f sintcosfl is the nonlinear parametric excitation
term, and 1/r25in0 is the autonomous nonlinear part. The angular displacement
9 is in 5'1 space (one dimensional Sphere Space), whereas the angular velocity and
acceleration are in RI space.

In the Simulation study, it is assumed that we know little of the system except
for the parametric identification. To apply the identification algorithm [24] to the
pendulum system, the following more general expression of a single d. o. f. system

can be assumed:
Grit + kLL‘ + fanl($,i‘) + {d1 + fpnl(17, i‘)}p(t) = —:i, (2.2)

where k is the linear stiffness parameter, fanl(.r,:i:) is the autonomous nonlinear
part, fpnl(x, 1'2) is the time-independent part of the parametric excitation term, and
coefficient (11 is the amplitude of external excitation force. If all is non-zero and fpnl
equal to zero, the system becomes a nonlinear system with external excitation force.
On the contrary, if (11 is zero and fpnl is non-trivial, it is then a parametrically excited
system. For the examined single pendulum, parameter k = 0 and d1 = 0. Meanwhile,
Since the non-linearity of the pendulum is caused by geometry, it is convenient to
assume fan] and fpnl in equation (2.2) as functions of only displacement, such that

it becomes

C7413 + fanttv) + fpn,l(x)p(t) = —'T (23)

12

Based on equation (2.3), the harmonic excitation term p(t) = (11 cos wt + [31 sin wt

with excitation frequency (12. Unknown functions fpnl and fan) can be approximated

by
N
fpnl % Z Pia-(x), (2.4)
i=1
N
fanl z 2 Qi¢i($)a (2-5)
i=1

where {oz-(:2: )} is a set of orthogonal basis functions, and P1, q,- are unknown parame-
ters. Since the choice of basis functions can affect accuracy of the identified functions,
two sets of basis functions were tested in the Simulation: linear interpolation functions
and Fourier series functions.

Substituting (2.4), (2.5) and the harmonic excitation into (2.3), we have

N
cri: + Z q,¢>,(:z: )+ 2 mega: x)cost + Z 171%“? :r)sint = —:'z}, (2.6)
i=1 i=1

where n; = Piezl and p,- = P1131. With equation (2.6), harmonic balance identifica-
tion [24] can be applied after extracting unstable or stable periodic orbits from the
collected displacement data. For a period k orbit, if displacement :1: is in R1 (one
dimensional real space), the displacement, velocity and acceleration signals can be

approximated by the following truncated Fourier series expansions:

l ' ' t
Ik(t) ~-a0k+ 3:101 “003th +bj ,kS SinL:) (2-7)
M .
lo“; 37:3(a — aJ-ksinJ—w—kt+bj kc cosj%t), (2.8)
j—l
t
m ~: —j2w2/k2(ajkcos]—-m +b,,,sin’“’), (2.9)
j— 1 k k

where, for non-dimensional equation (2.6), w = 1. However, for :c in 5'1 (one di-

13

 

 

 

 

 

 

 

 

 

2 2
1 0 ’ 1
Q 2 2 ' 4
g) 0. U)
C
m m 4 - ]
1
6 .
2 t 8 .
3 x . 10 1 . 1
0 (15 1 1.5 2 0 (15 1 1.5 2
peflod pefiod
(a) (b)
10

 

 

 

 

 

32
C)
E

5 .

10 '

15 * ‘ ‘

0 0.5 1 1.5 2
peflod

(C)

Figure 2.1. A period-4 orbit (a) sampled signal; (b) actual continuous signal 136; (c)
- - - constant rotation part wckt, - oscillatory part rum».

mensional sphere space, e. g. the angular displacement 0 in the Single pendulum),
Equations (2.7) — (2.9) are not applicable for whirling. The harmonic balance method
can still be applied here, but with proper modification. The issue is that the velocity
and acceleration are R1 signals, and cannot be derived directly from the sampled S1
displacement signal.

Consider an experiment in which whirling angle 6 Signal is collected with an en-
coder such that the output :I: has values in the interval [—7r, 7r). Thus, a: is discontin-

uous. According to the illustration in Figure 2.1 of a period four orbit of the single

14

pendulum system, the sampled data {as} should be converted to a continuous Signal

{sec} in R1, and then decomposed as

xck = wckt + (Evankv

where wckt represents the constant rotating part of a whirling periodic orbit, and
3‘00ch is the oscillatory part of the orbit signal, which can be approximated by the
truncated Fourier series expansion

1 M

jwt jwt
xvar,k(t)~ §U0k+zl( ujkcos— 1: +11, ksin k).
J=l

Based upon xck, the corresponding velocity and acceleration can be expressed as

' wt ° t
t) zwck + Z— J——w wksinj—w +vj kcos— Jw ), (2.10)
j=1 k( uj k k

cat t
+vj ks1n‘7w (2.11)

M
site) a Z —12w2/k2(u,cok >

SJ__W
i=1 k
Meanwhile, for a period-k orbit, if ¢;(r) : S1 ——> C[——1, 1], as in the case of the
pendulum, the continuous functions (bl-(3%), ¢i($k) cos wt and (25,021.) sin wt can also

be approximated by

M . .
wt
o,- (xk)~ ~ CtO k/2 + 2(C CU- kcos J:——t + d-- jksinJT), (2.12)
J'2 —1
M J'wt J'wt
o (:rk)cost~ ~ eiOk/z + :(e ez-J-Jccos— +fij ksin— k ), (2.13)
j:1
M . .
t t
d>-(:1:k)sint~g,-Ok/2+ Z( (gijkcosfll: +11, j,ksin 2:— ). (2.14)
j=1

15

Substituting (2.10)—(2.14) into (2.6), the differential equation can be transformed

into the following matrix equation for a period-k orbit

 

 

 

 

 

 

{ wet 6101/2 9101/2 0101/2 €N0,1~/2 9N0,k/2 emit/2]
fume €11,1c 911,1c 011,1: 8N1,k 9N1,k CN1,k
Jfflu: f11,k h11,k d11,k fN1,k thJc leJc
$11141 e1M,k 911m ClM,k eNMJc 9NM,k CNMJc

\‘A‘iguMJc f1M,k h1M,k d1M,k fNM,k hNMJc dNM,k/

( Cr \

721 f 0 \
P1 w2u1,k/k2

(11 (flunk/[C2

"N M2w2'uMJc/k2
pN \M2w2vMJc/k2 )
\‘IN/

(2.15)

For Simplicity, A]: is denoted as the left-hand side matrix,
)‘=(C7‘ "1 P1 91 "N PN 4N )T’
and
>1".

Bk=(0 w2u1,k/k2 w2vl,k/k2 M2wQUM,k/k2 M2w2vMJc/k2

Equation (2.15) can then be expressed as AkA = Bk.

For multiple periodic orbits, by combining the Single periodic orbit matrices to-

16

gether, we can obtain the following matrix equation

211:5, (2.16)
where
= T T T T
A (A1 Ak AK),
T
,8=(B¥.‘ ... ’62; ... [81:19) ,

and A is a K(2M +1) x (3N +1) identification matrix . When K(2M + 1) > 3N +1,
the matrix can be solved by least mean square method. The least squares solution of
(2.16) is

A = (ATA)‘1ATB. (2.17)

The error of the identification algorithm can be affected by several factors, includ-
ing accuracy of unstable periodic orbits, noise contamination, system nonlinearity

and the choice of basis functions of approximation.

2.3 Introduction to phase plane reconstruction
and extraction of unstable periodic orbits

Usually in an experiment, only limited Signals can be acquired accurately, e. g. the
angular displacement signal in this Simulation. However, for the extraction of unstable
periodic orbits, phase plane information is necessary. It is then that the phase plane
reconstruction technique [52] is applied. Suppose that 3(k) is the sampled smooth
Signal from dynamical systems. Smooth dynamics in an n dimensional phase Space

could be approximately represented by an embedding dimension space

S(k) = {s(k),s(k + Td), . . .,S(k +(d—1)Td)},

17

where d is the dimension of the reconstructed phase Space and Td is the time delay
of the embedding dimension. Although the reconstructed phase Space is a distorted
appearance of the real phase plane, it provides us information of phase orbits, at-
tractors, periodic orbits and other chaotic characteristics. To the pendulum system,
s(k) = 0(k), 6 E [-7r,7r) is actually a non—smooth observable. But the method can
also be applied here Since the functions of angular displacement, Speed and accelera-

tion are all smooth.

2.3.1 Choosing the time delay and embedding dimension

The time delay value Td can be obtained from average mutual information function

[52], which is expressed as

P a

123m = game. WOMW, (2.18)
where t is the time lag, a. is the sampled data 3(1), A is the set of {3(1)}, 1: = z' + T,
b is thus s(z' + t), and B is considered as the set of {3(2' + t)}. Meanwhile, PACE) is
the probability function of observing .7: out of a set A, PB(y) is the probability of
observing y out of a set B, and PAB(:c,y) is defined as the probability of observing
function of (:r, y) out of the set of A n B. After calculating I (t), the ‘best’ time delay

Td is chosen when I (t) reaches its first minimum.
To determine the value of d E for minimum required delay dimensions, the false
nearest neighbors method was applied [52]. The false nearest neighbors method iden-

tifies the embedding dimension for which false trajectory crossings do not occur.

2.3.2 Extracting unstable periodic orbits (UPOS)

There are numerous unstable periodic orbits in a chaotic Signal. The UPO set is a

dense set within the chaotic orbits. The reconstructed phase space can be used to ex-

18

tract the hidden UPOS. Since the horizontally excited Single pendulum is excited by
a harmonic Signal with period T, the UPOS will have the periodicities of integral mul-
tiples of the excitation period T. Though theoretically, no exact UPO can be found
through the collected data, the theory proved that some very close approximation of
the UPO exists provided that the data is long enough. With the error tolerance 6
set for UPO extraction, we say that an approximated UPO of period It: is extracted if
|S(n — kT) — S (n)| < e, where S (n — kT) is the starting point of the extracted orbit,
and S(n) is the end point. Usually e is set as 1—5% of the Span of the data [24], and
the smaller 6 value is always desirable if possible, since the corresponding UPO will
be closer to the real periodic orbits.

Meanwhile, a recurrence error between the starting and end points will occur in
the extracted UPO by the choice of e, which is an exhibition of the accumulated error
compared to the ‘true’ periodic orbit. Generally, to acquire more accurate UPOS,

smaller 8 is used.

2.4 Numerical Simulation of the horizontally ex-

cited pendulum

2.4.1 Phase space reconstruction

The Simulation was based upon the non-dimensional differential equation (2.1). A
data set of 30000 points was gathered with a sampling rate f3 = 25/T. Displayed in
Figure 2.2 (a) is the 2-D reconstructed phase Space. The coefficients were chosen as
f = 1.52, c = 0.03, r = 0.8333. The actual minimum embedding dimension needed
for a complete phase space reconstruction was found to be four. However, Since our
purpose of reconstructing phase Space is to extract UPOS, which compare points with

[CT time interval (implying the addition of the SI dimension of time), the possibility

19

 

 

   

 

 

 

 

 

 

4 4
2 2
s E
s 0 E 0
3 o
5
-2 -2
-4 -4
-4 - o 4 - — o 2 4
theta(t) theta(t)
(a) (b)

Figure 2.2. Reconstructed phase plane (a) and Poincaré section plot (b); 0(t) is
represented by theta(t) in the plots.

of encountering a false nearest neighbor is greatly reduced. Hence, two embedding
coordinates are adequate for this case, and are also simpler for display purposes.
According to Figure 2.3 of the average mutual information, the calculated best
time delay was Td = 5. The phase portrait and the Poincaré section are displayed
in Figure 2.2 (a) and (b). The Poincaré section plot, obtained by plotting the cross-
section of the phase flow at a fixed periodic time position T = t mod 27r since the
system is under periodic excitation, provides us a handy tool for visually comparing
chaotic properties of different systems. Similar Poincare section plots give evidence

of similar dynamical behaviors of the systems.

2.4.2 Parametric identification

For Simplicity, the data used in the identification process was noise-free, thus excluding
the noise-generated error in the identified parameters. The UPOS were extracted

before applying the identification algorithm. An error tolerance of e = 3% was

20

Mutual Information
2.5 l I l I

 

1.5r

I(dt)

0.5 *-

 

 

 

 

Figure 2.3. Mutual information of the data; the first minimum is at dt=5.

applied in the extraction. The error tolerance can be made smaller if longer sampled
data is available. Consequently, around 25 UPOS were extracted from period 1 to
period 16. Figure 2.4 Shows a example of periodic orbits of period one and two. Some
are whirling orbits, whereas others are oscillating orbits. With the UPOS extracted,
in the identification process, two sets of basis functions were tested to approximate

the unknown nonlinear functions in the governing equation.

21

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 4
3
1
2
E E
31 s 0
1
2
2
3 3
3 2 1 0 1 2 4 2 0 2 4
theta(t) theta(t)
(a) (b)
4 4
2 2
€ 55‘
.t. .t
E 0 E °
0 o
5 5
2 2
4 4
4 2 0 2 4 4 2 0 2 4
theta(t) theta(t)
(d)

(C)

Figure 2.4. (a)A period-1 orbit; (b) A period-2 orbit; (c) A period-2 orbit; (d) A
period—2 orbit.

Harmonic basis functions
If basis functions are set as ¢2k_1(r) = cos kz, ¢2k(:r) = sin kx, k=1, 2, 3, ..., the

differential equation (4.12) is then converted to

crrb + ZkK=1(n2k_1 cos k3: + "2k sin km) cost + Zf=1(p2k_1 cos k2: + P2]: sin kz) sint

+ 25:1(‘12k—1 cos k2: + ‘12]: sin kz) = —i.
(2.19)

22

Table 2.1. The "k1 Pk and qk values when harmonic basis is applied.

 

 

- 7% n2 Ink-fit 19;, pk lpk-fikl 1);, (1;, lqk-ékl

 

coskO 0.009 0.0 0.009 -1.526 -1.52 0.006 -0.004 0.0 0.004

 

sin/c0 0.003 0.0 0.003 -0.0 0.0 0.0 1.434 1.44 0.006

 

 

 

 

 

 

 

 

 

 

 

 

For the investigated pendulum, it is convenient to set K = 1, i. e. include only the
first harmonics, Since the differential equation is simple and consists of harmonic type
non-linearity.

For systems like the Single pendulum, where the non-linearity comes from angular
rotation and the displacement variables belonging to the sphere space, a harmonic
basis can be a very good choice due to its capability of representing rotation and
periodic behaviors. Displayed in Table 2.1 are the identified coefficients of (2.19)
with fik, 13k, 91: representing the estimated values, and nk, pk, qk representing the
true values. The estimated or value is 0.0339. Compared to the real values, all of the

parameters are quite accurate with minor errors less than 0.01.

Linear interpolation basis functions

Linear interpolation functions can be expressed as

Sig—1E, (k—1)d<.r<kd
21(22): @331‘1—“9, kdgx<(k+1)d,

0, otherwise

where d is the distance interval. Had we not known the apprOpriate basis func-
tions, e. g. harmonics in this case, interpolation functions are good alternatives [21].
For better expression of the nonlinear functions, generally more interpolation points

(functions) are needed [25]. However, that implies more unknown parameters and

23

I

0.4

0.2 ~

nk
O
1
pt
I»
*-

 

__1 1 1 1 1 1
-4 -3 -2 -1 o 1 2 3 4
theta position

.-
h.
1—

Figure 2.5. The "k curve for linear interpolation basis; (*) identified values, (—-—) true
function.

requires more periodic orbits for the least mean square method. Eight points were
utilized in this identification. The interval [—7r, 71) was divided evenly into eight equal
sub-intervals. The unknown curves can be represented by the piece-wise linear curve
connected between sub-interval nodes.

Figures 2.5, 2.6, and 2.7 Show the piece-wise curves of {11k}, {pk}, {qk}. Mean-
while, the identified friction coefficient cr = 0.0348. It can be seen that {nk} curve
is nearly zero, {pk} curve is a cosine curve, and the {qk} curve is a Sine curve, which
are consistent with the differential equation (2.1). The eight-point approximation is
still a rough representation of the real curve, and consequently introduces more error
than the harmonic basis case. In the next section, we will discuss the error’s influence

on the dynamic characteristics of the system.

24

1.5~

*2

 

_2 1 1 m 1 1 1 1

-4 -3 -2 -1 1 2 3

#1

0
theta position

Figure 2.6. The Pk curve for linear interpolation basis; (*) identified values, (—) true
function.

25

I

0.5

-1.5L 0

 

_2 l l l 1 L l I j
-4 -3 -2 -1 o 1 2 3 4
theta position

Figure 2.7. The qk curve for linear interpolation basis; (0) identified values, (—) true
function.

26

2.5 Discussion and validation

Direct comparison showed that the identification algorithm was satisfying for the
Single degree of freedom system. Furthermore, it was the harmonic basis that gener-
ated the more precise identification result, due to the fact that the unknown nonlinear
functions consist of harmonics. However, when investigating unknown systems, direct
comparison is unavailable. In this case, it is then natural to examine the dynamical

behavior, e. g. chaotic characteristics, of the identified system for validation purposes.

 

 

   

 

 

 

 

 

 

4 4
2 2
f. 3
E O E 0
0 0
5 5
-2 '2
'4 '4
-4 —2 o 2 4 -4 —2 o 4
theta(t) theta(t)

(a) (b)

Figure 2.8. Poincaré section plots of identified systems: (a) by linear interpolation
basis, (b) by Harmonic basis.

Displayed in Figure 2.8 (a) and (b) are Poincaré sections of the identified systems
with harmonic basis and linear interpolation basis. The Poincaré sections look similar
to the original pendulum, which implies the dynamical similarity of these three chaotic
systems. Thus, for the purpose of chaotic characteristics, the linear interpolation
approximation, though a coarse representation, is still acceptable.

For more validation methods of the unknown systems, we can also refer to the

linearized model method, comparison of fractal dimensions, bifurcation diagrams, and

27

Lyapunov exponents. But for this simulation case, where the parameters are known,

direction comparison is the best tool.

2.6 Concluding remarks

In this chapter, the simulated parametric identification procedure was investigated
for a chaotic horizontally base-excited pendulum. The identification was assumed to
apply to a Single-degree-of-freedom system with a known parametric excitation fre-
quency and unknown nonlinearity. The identification algorithm [24, 25] was modified
for application to the whirling pendulum system. Two types of basis functions were
applied to approximate the unknown functions. Then, the simulations of the identi-
fied models were presented to verify the effectiveness of the algorithm. Both the direct
comparison and the validation methods showed the accuracy of the identified param-
eters and the Similarity between the original and the identified systems. Although
both of the methods gave Similar results, the harmonic basis functions matched the
form of the true nonlinearity, and therefore in this case, gave a more precise result.
On the other hand, the linear interpolation was also valuable since it is a good choice
for a general unknown type of non-linearity in a Single variable. According to the
simulation result, the modified parametric identification process was successful in the
single degree of freedom whirling system. However, the Simulation was limited to
single degree of freedom systems with no noise contamination. For more general ap-
plications, multi—degree of freedom systems and experimental verifications are needed

for full investigation of the proposed method.

28

CHAPTER 3

Experimental Study of a

Base-Excited Double Pendulum

The approximated periodic orbits can be extracted by phase space reconstruction
techniques [53, 54, 52] with a certain tolerance of error for extraction. The approxi-
mated UPOS made little influence on the identification of the single pendulum system.
However, it is not clear whether it is so for multi—degree-of-freedom systems, or sys-
tems with strong nonlinearity. Also, it is necessary to examine the harmonic balance
based identification method in experimental systems in which the collected signals
are noise contaminated. These are the reasons for the continued research on an ex-
perimental multi-degree-of-freedom system: the double-pendulum system. Another
reason that the double-pendulum system was chosen for experiment is that it is easier
to adjust parameters of the double pendulum system to exhibit chaotic behavior than

the Single pendulum system.

3.1 Description of the double pendulum system

A schematic diagram of the double pendulum is shown in Figure 3.1. The first arm

has mass ml, centroid offset 61, arm length 11, and angular inertia J61 based upon

29

Vertical base
excitation

    

Central arm

Figure 3.1. A Sketch of the double pendulum.

the arm centroid point. The second arm has mass m2, centroid offset 82, arm length
12, and angular inertia ch. 91 is the absolute angular displacement of the first arm
and 62 is relative angular displacement of the second arm.

The two arms of the pendulum are supported and connected by low friction bear-
ings. The bearings are assumed to have two types of friction: dry Coulomb friction
and viscous damping. To the specific double pendulum that was used in the exper-
iment, tests indicated that the first arm bearings had dominantly Coulomb friction
because of no oil lubrication, and the second arm bearings had dominantly viscous

friction due to full oil lubrication. With these known properties, we can then obtain

30

the non-dimensional governing differential equations of this system:

f

$1+ b11$2 COS ($2 — $1) - b11$§ sin ($2 — $1) + 012 8111 $1
+1313 Sin 0511? + 61182971051) - 0120152 - €51)= 0 (31)
(152+ b21051 COS ((132 - (151) + 52145? sin (<22 - 41) + 1922 sin 42

+b23 Sin (2523? + 02(452 — (151) = 0

 

where 451 = 01 and $2 = 01 +02 are absolute angular deflections, 311 = m2e2l1/J01,
1112 = 9(m181 + m211)/J011 513 = (W161 +m211)/Jo1 1 321 = m282’1/J021 1’22 =
9m262/Jo21b23 = m262/412 1012 = Cr2/J01162 = C'r'2/‘1021 J01 = Jc1+m1€f+m2lf

, and J02 = J62 + 172283; y is the excitation displacement of the support with known

7

1, :r > 0
frequency fe; and f (.13) = sign(:z:) = l 0, :1: = 0 is a Sign function representing
—1, a: < 0

 

Coulomb friction. Function f (:12) is valid if there are no sticks. For convenience of
the analysis, a non-dimensional form of the governing differential equation is desired.
By letting 7' = 27rth and Q = 27rfe, then, 52 = if; - 931% = 27rfejd; = {2%. Under a
Sinusoidal excitation y = acosr, where a is the excitation amplitude, equation (3.1)

can be expressed as

f

451" + 311%” COS ($2 — Q51) - 1311452'2 sin (¢2 - 451) + 312 Bin $1
< +313 sin my” + C118ign(¢>1') - 012(42' - 41') = 0 (3 2)
$2" + 321051" COS (¢2 — <21) + le¢>1’2 sin (42 - 051) + 322 sin $2

( +323 Si114521!” + C2(¢2' - 451') = 0

 

WllCI‘C $11, = ad;$i,Bi1 = 02:1, 822 = 012/92 for i=1, 2, 823 = —bi3a, C11 = 611/92,
012 = 012/!) and C2 = c2/w. The differential equation is non-autonomous because
of the time dependence of y(r).

For identification purposes, the excitation Signal y can be approximated by the

31

Fourier series expansion of
J
y(r) = 2: (EJ- cosjr + Fj sinjr), (3.3)
J =1
where the coefficients Ej and F J- are unknown. For the present system with harmonic
excitation, y = E1 cos 7' + F1 sin 1' is adequate for identification. Since all the param-

eters are unknown with the only exception of the excitation frequency, equation (3.2)

can then be transformed into

31140” COS (4)2 - ¢1)- 13114228in ($2 - <01) + 1312 Bin $1
+313 Sin (1)1 eosr + 814 sin 421 sinr + Cllsign(d>1') - 012(cb2' — 451’) = —¢1”
321¢1IICOS (452 — $1) + 1321151,2 Sin ((152 — 451) + 322 Sin 432

+853 Sin $2 cos 1' + 8&4 sin (131 Sin 7 + C2(¢2’ _ 4,1!) ____ _¢2//

1

 

(3.4)
where Bg3 = 313E1 and Bg4 = Bi3F1 for i=1, 2. For Simplicity, we will use 323 and
82-4 to replace 8:3 and 8:4 in the later parts of this report. Equation (3.4) is then

the desired form for identification.

3.2 Method

The identification process is similar to the one used in the simulated parametrically ex-
cited Single pendulum (Chapter 2), which contains data acquisition/post-processing,
phase plane reconstruction / extracting unstable periodic orbits, formation of the iden-
tification matrix and the solution by the least mean square method. Due to apparatus
limitations, some signal noise occurred during digitization. Therefore certain digital
filtering techniques were applied to the acquired Signals. Details of the post-processing
will be introduced in next section. Other modifications are also applied to the double

pendulum system for improvement of the identification.

32

Due to the complexity of the non-linearity of the double pendulum equation (3.4),
our unknown parameters contain only these coefficients of terms of the differential
equations, which will greatly Simplify the identification matrix. There are totally
11 unknown parameters in the two differential equations, namely, Bu, 812, ...,
824, Cu, 012 and C2. Similar to the case in the Single pendulum simulation, the
angular displacements 61 and 02 are both variables in S1 (one dimensional sphere
Space). However the angular Speeds, accelerations and Sin (b,- (2' =1,2) belong to R1
(one dimensional continuous real Space). Hence, for any period-k orbit (there may
be multiple orbits of same periodicity k, the Fourier series expression of the periodic

orbits is

jw , jwt
$1.1,tt()~91,,,,k1t+00k1/2+XX 03,,-—sz08 k t'+b 31151117) (35)
j: —1

and

jw , jwt
$2,kl(t)~ 92Mk1t+60k1/2+ Z(c cjkjcos k t+djk,ls1n—E-), (3.6)
j=1

where 91,“, is the average rotation speed per cycle for the lth orbit of period k. By

doing so, the following equations can be obtained through the Fourier expansion of

periodic orbits:

jw

 

- m jw , jwt
4511,10 )~91,,kz+ Z 73(— “3,2,1005— k t'+b 32181117). (37)
i=1
- jw . jwt
¢2,k,ll’"z() “9211+ 2:3]: -C3',k,tCOS—— k +61 k3,,zsm-k—L (38)
j: —-1
m -2 2
-' J w J'w J'w wt
1512,10)“ '21— ,2 (03210087 t+b3,,kzsin T) (39)
J:
m 2 2 -
u w
$2,k,l(t) 2 Z -] 26(jleOSJ—k— +djk1Sln k ), (310)
J=1

  

33

jwt

 

 

 

 

.. €1,0,k,l J'wt .
¢-i,k,l(t) C05 (¢2,k,z ‘ $1.191) ’e“ 2 + Z (ei,j,k,l COS T+fi,j,k,l 3111 k ), (3-11)
i=1
-2 . 91,0,k,z m jwt , jwt
4523M“) 31“ (<92,k,z — $1M) * 2 + Z (9i,j,k,l COS T + hi,j,l,k sm T)’ (3-12)
i=1
. 3 191,0,“ m jwt , jwt
szgn(q§i,k,l(t)) z 2 + Z (pi,j,k,l cos —k— + qz‘,j,k,l 8111 T)’ (3.13)
i=1
. 7‘1,0,k,l m jwt , jwt
Sln ¢i,kl z 2 + Z (ri,j,k,l COS T + Si,j,k,l s1n T), (3.14)
i=1
. U1 k m 'wt . 'wt
sm (25,-,ch cost a —%—’—l + Z (“1'03“ cos 1;- + ”23ij sm J—k—), (3.15)
1:1
. . w m 'wt . 'wt
srn¢i,k,331nt a: Jig-lg + Z (wi,j,k,l cos 37 + Zi,j,lc,l sm '77), (3.16)
1:1

where k=1, 2, ..., K, and i=1, 2, is the corresponding period of the orbit; K is the
maximum periodicity. Meanwhile, for the non-dimensional differential equations, time
t here is actually T in equation (3.4), and fundamental frequency is 1. If incremental
encoders are used to sense the angular displacements, then velocity and acceleration
are not directly measured. Their Fourier transforms can be handily generated by the
displacement’s Fourier transform for periodic response. However, for noise contam—
inated displacement signals, noisy errors may be amplified in the obtained Fourier
spectrum of velocity and acceleration.

Also, Fourier series expansions of some terms in the differential equations, such as
(1322,]: sin (432$, — 431,16) , are not obtained by direct Fourier series expansion of their time
domain signal, but by convolution of known Fourier expansion components with a low
pass filter applied. The purpose of applying a low pass filter to each signal component
before convolution is to avoid noise amplification. Reasons will be explained in the

following sections. Substituting (3.5—3.16) into (3.4), and equating the coefficients of

34

terms with identical harmonic order, we obtain two matrix equations:

01‘

l (6 — 9)2,0,1,1
(6 - 9)2,1,1,1
(f - h)2,1,1,1

(8 — 9)2,0,k,l
(8 - 9)2,1,k,z
(f - h)2,1,k,z

 

7'1,0,1,l ul,0,l,1
Tl,1,1,1 ul,l,l,1

31,1,1,1 U1,1,1,l

T1,O,k,l “1,0,k,l
Tl,1,k,l “1,1,k,l

Sl,l,k,l v1,l,k,l

(Blll

312
313
B14
C11
\ 012 )

 

 

w1,0,1,1

w1,1,1,1

zl,l,1,1

w1,0,k,l

wl,l,k,l

 

21,1,k,1

01,1,1

51,1,1

9 ,

81
H

II
:91

35

pl,0,1,1
pl,l,l,1

q1,1,1,1

p1,0,k,l
P1,1,k,l

ql,l,k,l

\.

 

91,1,1- Q2,1,1
d1,1,1 - b1,1,1

-Cl,1,1+ a1,1,1

91M - 92M
d1,k,l - b1,k,l

-C1,k,z + “I.“

\

 

(3.17)

(3.18)

and

 

 

 

 

 

 

{ (e-g)1,o,1,1 7‘2,0,1,1 'U2,0,1,1 w2,0,1,1 -91,1,1+92,1,1 )
(6-9)1,1,1,1 1"2,1,1,1 “2,1,1,1 w2,1,1,1 d1,1,1-191,1,1 (B \
(f-h)1,1,1,1 S2,1,1,1 v2,1,1,1 22,1,1,1 -(01,1,1-ai,1,1) 11
. . . . . 312
313
(6-9)1,o,k,z 72,031,: ”2,0331 w2,0,k,z 91M —92,k,z B14
(e-9)1,1,k,z ’"2,1,k,l “2,13,: w2,1,k,l d1,k,l -b1,k,z \ C
U -h)1,1,k,z 32,1,k,l “02,1131 22,1,k,z -(61,k,z —01,k,z) 2 j
( s 2 2 s s )
( 0 )
C1,1,1
d1,1,1
_ 0 ,
CLka—z
deJk—2
\ l (3.19)
01'
A252 = (3'2, (3.20)

where 51 and 52 are vectors of unknown parameters. Given a set of periodic orbits, it
is adequate to equate the coefficients of the first several orders of sub-harmonics since
they are usually the major components of the periodic orbits and less contaminated by
noise. We then truncate the Fourier series expansion, and take the first M harmonics,
such that (2M + 1) - K > NC, where (2M + 1) - K is the total number of rows in

the matrix A.,-, and NC is the number of unknown coefficients in 551 or :32. With

36

these conditions satisfied, the two equations can be solved by the least mean square
method:

5‘51 =(A1TA1)‘1A1Tc71 (3.21)

and

52 = (AgA2)—1A§If2. (3.22)

After the non-dimensional parameters are identified, the physical properties can
then be restored according to non-dimensional parameters’ definition if part of the
physical parameters are known prior to identification. In the experiment, mg, 82, [1
were treated as known since these physical properties were easily evaluated. Also, by

defining 325 = 8,23 + 82.24, i=1, 2, we can obtain

9
Pb = lBiZ/BiSI = (ZerPa- (3.23)

Hence Pb is a constant for a constant excitation amplitude a and can be used as an

indicator of the accuracy of the identified results.

3.3 Experiment description

Figure 3.1 shows a sketch of the double-pendulum that was used in the experiment.
Two optical encoders (US Digital) were separately attached to the central arm (US
Digital H58—1024-157H) and the second arm (US Digital HSS-1024) to measure the
relative angular displacements 01 and 62. Both of the encoders had a resolution
of 1024, and were therefore capable of detecting a minimum angular difference of
0.35160. The two encoders sent out TTL square waves, which are noise-resistant, to
two EDAC (Encoder Digital to Analogue Converter) converters, which transformed

the TTL waves into analogue signals. After that, a data acquisition terminal trans-

37

lated all the signals into computer-acceptable digital signals.

Table 3.1. Physical properties of the double pendulum.

 

 

 

 

 

 

 

 

 

 

ml (kg) 0.1362 mg (kg) 0.040
e1 (m) 0.0127 e2 (m) 0.0267
11 (m) 0.0635 12 (m) 0.0534
J1(kgxm2)* 6.99x10-4 J2 (111;,me 4.033x10—5
Cu 1.01 x 10-3 Cl * 0.0485
02 * 0.00366 — —

 

 

For validation purposes, Table 3.1 lists all the physical properties of the double
pendulum. The asterisks in the table denote that some properties are not directly
measured, but estimated from other dynamic methods, which implies that those pa—
rameters could have small errors. To estimate the mass moment of inertia, a small-
amplitude free vibration was tested on the double pendulum. By evaluating the two
natural frequencies of the system through the FFT, the mass moment of inertia values

were calculated. Some parameters related to sampling and experimental settings are

listed in Table 3.2.

Table 3.2. Experimental settings.

 

Sampling rate (f3)

500 Hz

Excitation freq. ( f e)

5H2

 

 

Cut-off freq. (fc)

 

 

80 KHz

Excitation amplitude (a)

 

1.15 cm

 

With all these settings, data were obtained during a three-hour-long chaotic vi-

bration. The data record lasted 22 minutes.

38

 

3.4 Results and validation

3.4.1 Phase space reconstruction & UPO extraction

The embedding dimension was chosen to be eight by the false nearest neighbors
method. Mutual information [52] of the signal was used for choosing adequate time

delay Td-

Mutual Information
0. 1 8 I I I I I I I I I

 

0.16

0.14

0.12

0.1

I(dt)

0.08

0.06

0.04

0.02

 

 

 

 

0 5 10 15 20 25 30 35 40 45 50

Figure 3.2. Mutual information I (dt).

It can be seen from Figure 3.2 that there are weak minima of I (dt) at dt =5,
15, 18 and 24. But, Td=24 in a driving period of 100 samples is somewhat close

to a quarter period, the ideal delay for a sinusoidal signal. The reconstructed phase

39

 

 

 

 

 

 

 

 

 

 

 

 

    

  
 

   
  

 
  
   

 

              

 

 

 

4 4
Ad
A 2 A 2
3 ‘13
E 0 § 0 5(2)
3 E 34
2 2 '
7 |/' /
4 4
4 2 o 4 4 o
theta1 (t) theta2(t)
(a)
4
1 \M 1“?)
2 ifdllflwfik‘m'lillll
A r /~ ywun' v '
a 41!!!!M/Alhzllzfillél.
... I ‘ ' .
. (1|)Mmlliiml
‘ .31.!le 1|
n, ’§. 7 V "1
4
4 2 o 4
theta1(t)
(C)

Figure 3.3. Phase portrait of experimental data with dt=24, (a) 01(t)—91(t+ dt), (b)
01(t)—02(t), (c) 92(t)—62(t + dt); 9 is represented by theta in the figure.

portrait is plotted in Figure 3.3 with Td = 24. The portrait shows that the central
arm, represented by 01, oscillated in small angles most of the time with occasional
large-angle whirling, which implies relatively larger noise in the 01 signal due to the
limitation of the optical encoders. The second arm displacement, represented by 92,
consisted mainly of whirling vibration.

After choosing the embedding dimension and time delay, we used the reconstructed
phase plane to extract unstable periodic orbits (UPOS). The error tolerance of extrac-

tion e was set at 5%. Supposing that the sampled signal is s,(t) = [91“) 92'“ + lel

40

for i=1, 2, an approximated UPO with error tolerance 19 is extracted if
||s,-(t) — 3,-(t + kT)” < 6132' = 1,2, (3.24)

where T is the excitation period and k is periodicity of the recurrence. For the
given experimental settings, f3=500Hz, fe=5Hz, T=fs/fe=100. Then, for example,
for a given k = 4, if inequality (3.24) is satisfied, a period 4 orbit is then said to
be extracted. One data set of 670,000 points was used in the analysis. For a 5%
tolerance, 6 distinct orbits were extracted.

Figures 3.4—3.9 show all of the distinct extracted UPOS. In all cases, the small
arm whirled. In Figures 3.5 and 3.7, the central arm whirled, whereas in other figures,
the central arm oscillated without whirling. Compared to [24, 25], there were much
fewer extracted UPOS in this double pendulum system than those in the single d. o.
f. systems, which even had tolerance error much smaller than 5%. This is because
recurrences are less frequent in higher dimensional spaces for small periodicity 16.
Suppose the number of boxes of size 1" needed to cover an attractor of dimension
01 is N1 ~ Alra, where A1 some constant determined by the chaotic attractor.
Roughly assuming even distance of data over an chaotic attractor, the probability of
recurrently hitting a given box 191 is p1 ~ 1 / N1 ~ 1 /A1r—0’. Comparing this with

another chaotic attractor of dimension 6, we have

In £2 fi—a
N T ,
P2 A1

which implies that the higher the dimension of a chaotic attractor, the smaller the

probability of finding a UPO with tolerance of error 1".

41

 

Phase Portrait of UPO. period4 PhasePortrait of UPO. period4
v 4

 

 

 

 

 

 

 

 

1.5
1
A A 2
6 05 1‘5
2: .t
a . 3
3'5 1»
5 0. 5 5
2
1 {/3
1. 5 4
2 0 1 4 2 0 2
theta1 (t) that a2 (t)

(a)

Figure 3.4. A period 4 UPO; theta1 and theta2 in the plot represent 91 and 02;
dt=24; same for the following figures.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Portrait of UPO. period9 Phase Portrait of UPO. period!)
4 4
/ Z

2 2 /
5 6 y
.t I.
:: 0 RT 0
.‘S S
a) o
5 5

2 l 2 /\ /% /7

I
4 4
4 0 2 4 4 2 0 2
theta1 (t) theta2(t)
(b)

(a)

Figure 3.5. A period 9 UPO.

42

theta1(t+dt)

theta1 (t+dt)

Portrait of UPO. period12

 

 

 

 

0
theta1 (t)
(a)

theta2(t+dt)

O
\\

Phase Portrait of UPO, period12

 

 

 

/

 

/

/

\\)‘
\\

 

 

 

 

 

 

 

V

 

Figure 3.6. A period 12 UPO.

Portrait of UPO. period12

 

 

 

 

 

2 0 2
theta1 (t)

(a)

theta2(t+dt)

O
theta2(t)

Phase Portrait of UPC, period12

 

[lg/l

 

 

2

Figure 3.7. A period 12 UPO.

43

o
theta2(t)
(b)

2

 

 

 

Portrait of UPO. period15

 

o .3
in _I in

theta1(t+dt)
.0
U! D

 

 

 

theta2(t+dt)

N

0

Phase Portrait of UPC. period 1 5

 

 

 

 

 

 

 

2 2

0
theta2(t)
(b)

Figure 3.8. A period 15 UPO.

 

 

 

 

1. 5
2 0 2
theta1(t)
(a)
Portrait of UPO, period15
15
1
7.5" 05
it.
I: 0
S
d)
5 0. 5
1
1. 5
2 1 0 1 2
theta1(t)

(a)

theta2(t+dt)

Phase Portrait of UPC. period 1 5

 

 

 

 

o
theta2(t)
(b)

Figure 3.9. A period 15 UPO.

44

 

 

3.4.2 Identified parameters

Table 3.3 lists all of the identified parameters by Fourier series expansion of UPOS
with these enhancements: a low-pass filtering, FFT convolution of nonlinear terms
involving angular velocity and acceleration, the harmonic balance method and the
optimization. The reasons for the modifications will be explained in the following

sections.

Table 3.3. Identified parameters by applying low-pass filter, F F T convolution, Fourier
series expansion of UPO, harmonic balance method with sub-harmonics whose fre-
quencies are S excitation frequency fe=5 Hz, and sub-harmonics optimization.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

— Identified values ’h‘ue values Errorx 100%
Bu 0.1167 0.1131 3.2%
812 0.0781 0.0711 9.8%
313 0.0809 — -
314 0.0264 — -
315 _ 33%;, + Bl4 0.0849 0.0820 3.6%
1321 1.6149 1.6816 4.0%
822 0.2562 0.2630 3.6%
823 0.2777 - '-
324 0.0715 - -
325 = 3333 + 324 0.2868 0.2913 1.6%
(111 * 0.0001 1.02x10-3 -
012 * 0.0007 0.00366 -
02 * 0.0106 0.0485 -
J01 (kg-m2)* 5.8113 x 10-4 5.99 x 10-4 2.0%
J02 (kg-m2)* 4.1995 x 10—5 4.033 x 10"5 3.9%
a (cm) 1.18 1.15 2.6%
mlel (kg-m) 1.73x10—3 1.60x10—3 8.5%

 

 

The sub-harmonic terms whose frequencies were smaller than (or equal to) the
excitation frequency were selectively included by an optimization method discussed

later in the identification procedure, because the noise is actually reduced in this

45

frequency interval. The optimization process utilizes linear regression techniques
and will be introduced later. The friction coefficients have relatively much larger
error comparing to the errors of other identified values. Part of the reason for this
inaccuracy is due to the fact that both of the Coulomb friction and viscous damping
factors are much smaller than other factors in the same matrix equations. So an error
deemed small for a stiffness/ mass parameter may cause large errors of the damping
parameters. As such, the damping coefficients will be estimated later by other means.
However, despite the friction factors inaccuracy, most of other parameters match
with the actual values within an error range of 10%, which is generally satisfying for
experimental data.

Meanwhile, as we have mentioned before, some of the ‘true’ values (marked with
‘*’) listed in Table 3.3 were obtained by an indirect dynamical method, e. g. small
angle free vibration, and thus those ‘true’ physical parameters and related non-
dimensional parameters may also have some error. Thus, more verification methods

were examined (see Section 3.4.6).

3.4.3 Friction issue

Our first result has shown that friction coefficients may not be precisely identified due
to their small values. Previous researches [24, 25] also indicate that friction parame-
ters were more difficult to identify accurately than other parameters in experiments
for a single degree of freedom system. The inaccuracy can result from mainly several

reasons:
1. experimental noise in the sampled data;

2. inadequate sensor sensitivity (the resulting error can also be treated as a noise

component);

46

3. inaccurate model of friction, e. g. viscous friction, dry friction, or their combi-

nafion;
4. recurrence errors of the extracted UPOS;
5. strong nonlinearity of the investigated system.

Meanwhile, since the damping factors are much smaller than others, and the method
of least mean squares is used for solving, an error deemed small for other parameters
can cause relatively large errors for the damping terms. Inadequate sensor sensitivity
can actually be considered as one of the noise sources during data acquisition process.
The noise decreases the accuracy of the extracted UPOS. As a result of the above
reasons, the identified friction coefficients are erroneous, whereas other estimated
parameters show only small discrepancies from the real values. On the other hand,
if the friction coefficients can be determined prior to identification, the identification
could likely be improved. Here, one would like to know the roles of friction parameters
in the identification process, and whether the friction errors have influence on other
identified parameters in the two-degree-of-freedom double pendulum system, which
is lightly damped.

We can assume that friction parameters are already known, and we want to iden-
tify the other parameters in our harmonic balance method. In this work, we de-
termined the friction in each pendulum bearing by a small amplitude free vibration
method. The free vibration test indicated that in the central arm, solely Coulomb
friction was involved, and in the second arm, solely viscous friction was involved. The
non-dimensional form of the friction parameters are listed in Table 3.1. As such, the
friction parameters were estimated by using free vibration decrements. Hence, we
applied the identification algorithm with known friction coefficients.

The result in Table 3.4 shows that the friction parameters, due to their much

smaller values compared to other parameters (less than 1/5 of other parameters),

47

Table 3.4. Identified non-dimensional parameters by optimized identification process

and noise reduction provided that the friction coefficients are known.

 

 

 

 

 

 

 

 

 

 

 

 

— Identified values 'h‘ue values Errorx 100%
Bu 0.1178 0.1131 4.2%
312 0.0774 0.0711 8.8%
313 0.0829 — -
814 0.0185 — -—
315 _ 33123 + B14 0.0852 0.0820 3.9%
321 1.6040 1.6816 4.6%
B22 0.2646 0.2630 0.6%
B23 0.2840 — '-
B24 0.0483 - -
325 = 3833 + 3%,, 0.2880 0.2913 1.2%

 

 

 

 

 

have little influence on the overall result. This was also shown in [25]. It implies that
the small errors in other parameters can cause a large percent error in the friction
terms. Hence, the first result in Table 3.3 is believed to be reliable for the coefficients
of conservative and parametric excitation terms.

For further verification, a simulation of the double pendulum system based on
equation (3.2) and the identified parameters was examined. However, the simulated
double pendulum system is extremely sensitive to parameter settings, e. g. friction
parameter values. The simulation was done in Matlab by digital integration. The
simulation result was obtained (shown in Figure 3.10) under Cu = 1.02 x 103,
Cu = 0.00366 and C2 = 0.0485 and other parameters set as the identified values in

Table 3.3.

3.4.4 High frequency noise in unstable periodic orbits

In the previous application of the chaotic system identification process ([24, 25], the

single pendulum simulation) in which single (1. o. f. systems were examined, the

48

 

 

N

theta1 (t+dt)
O

theta2(t+dt)
O

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 2
4 4
4 2 o 2 4 4 o 4
theta1 (t) theta2(t)
(a)
4
2 m /l
a, g
S 0
Q)
5
2
1 IV I ‘i
4 A 4; A
4 2 0 2 4
theta1(t)
(a)

Figure 3.10. Phase portrait of the simulated system with Cu 2 1.02 X 10’3,
011:0.0366 and C2 = 0.0485; dt=24.

49

identification algorithm was noise resistant. However, a similar hypothesis could not
be applied to the present two (1. o. f. system. The reason is simply due to the strong
non-linearity of the double pendulum, and it can be explained by the difference in
the governing differential equations of the investigated single and double pendulums.

The governing equation of the horizontally excited single pendulum in Chapter 2
is

8+c/r6+1/rzsin0— fsintcosO=0. (3.25)

Although the parametric excitation term is nonlinear, the 6 and 8 terms are, on
the other hand, linear. The sinG term in (3.25) can also be regarded as linear in
terms of harmonic functions. Suppose the contaminated signal is composed of the
real signal and noise n(t): 0(t) = 0t(t) + n(t). The angular velocity and acceleration
of periodic orbits can be obtained by equations similar to (3.7)—(3.10). Hence, by
equations (2.8) and (2.9), the noise in the obtained velocity and acceleration jth sub-
harmonic is actually amplified by jw/k and j 2w2 / 192, though n(t) is rather small in the
displacement signal. The overall signals of speed and acceleration are then considered
to be contaminated mainly by high frequency noise. Apparently, to the system (3.25),
where speed and acceleration terms are all linear, since in the identification matrix,
only the first K terms of harmonic order are used, the high frequency noise whose
frequency larger than K + 1th harmonic term will be automatically filtered out, which
explains why the identification process is noise resistant for systems like (3.25).
However, in the double pendulum system (3.4), the identification process ap-
pears to be less noise resistant than the previous examples since velocities and ac-
celerations do not appear linearly in differential equations. Specifically, it is the
(03;)2 sin ((152 - (351) and ¢,~” cos (432 — (251) terms that contribute most to the noise in-

accuracy of the result. Other high order terms also have similar problems of noise

50

amplification. It can be explained, for example, by the following F F T equation
3519’)? sin ((22 - (151)] = 350502 ® f[Sin (CD2 - 451)], (326)

where 27(2) represents the Fourier transform of :r, and operator <8) represents convo-

lution. The convolution operation can be expressed by
+00
12(3) 45 0(6) = / F(0)G(a — (5)65, (3.27)
—00

where F (a) and (1(0) are two integrable functions. By convolution, the high frequency
noise in each component is, therefore, mixed into the lower-frequency components
of the nonlinear terms. Furthermore, it actually amplifies the noisy influence of the
angular displacement, e. g. high frequency noise, since in the algorithm, the frequency
components of é,- and 6,- are obtained by equations (3.7—3.10). Hence, the truncation
of the first K harmonic terms could not reduce the disturbance of the noise.

One way to avoid noise contamination is to filter out the high frequency noise of the
displacement signals after UPO extraction and before convolution. By doing so, the
noise component can be effectively controlled. The algorithm was incorporated in the
identification process, and it turned out to be effective. For g = ((2592 sin ((132 — (251),
Table 3.5 compares the difference of the first four orders of Fourier series coefficients
of y with low-pass filtering before convolution and without filtering. The discrepancy
between the two methods is quite obvious and intolerable. It shows the necessity of
applying an adequate low—pass filter. On the other hand, we may also refer to some
optimization method such that only sub-harmonies with less noise contamination can
be selected into the identification process. The cut-off frequency was set to be 1/5f3
Hz. Figure 3.11 shows the velocity and acceleration frequency spectra of 91. It can

be seen that the acceleration’s high frequency noise is quite intolerable, and even

51

with the low pass filter of 1 / 5 f3 cut-off frequency, there is still considerable noise

remaining.

Table 3.5. Comparison of .7: (y) with and without filter added.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i Coefl'icients of cos (ix) Coefficients of sin (31:)
without filter with filter without filter with filter
1 -3.2486 -0.3087 0.0328 0.0000
2 -34.0713 -27.1330 17.7370 17.7403
3 -15.6946 -8.6732 -23.6908 -23.2652
4 -31.2079 -24.1696 -10.2571 -10.6034
01 0.5 '
8 0.08 '55 0.4»
8 8
U) (U
B 0.06 E 0.3-
m 2
:5: 0.04 3: 02i
O o :12: ......
E 0.021 ti- o.1i 3f
0 1 1 i 7" l. '1" L. . 0 . L . . .
0 100 200 300 400 500 0 100 200 300 400 500
k k

(a) (b)

Figure 3.11. FFT amplitude of the 6.1 and 01 with and without low pass filter applied;
the continuous line is F FT of signals with filter of 1/ 5 f3 cut-off frequency; the dotted
line is FFT of signals without filter; It is the order of sub-harmonics.

Displayed in Figure 3.12 are also signals with and without the low pass filter for
the period-9 UPO. The low pass filtered signals are more smooth and assumed to be

closer to real signals. Meanwhile, an identification process without any filtering was

52

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 1
g 05 g 0.5
‘5 ‘5
33 33
g a:
~ 0 '5 OA
0. 5 A A . . o. 5 A A A A
o 200 400 600 800 1000 o 200 400 600 800 1000
k k
(a) (b)
1 fl 1 . . .
2 3
>5 >1
1 m A A A A . A A
o 200 400 600 800 1000 o 200 400 600 800 1000
k k

(C) (d)

Figure 3.12. Signals without and with low pass filter applied; (a) and (b) are 6.1
without and with filter; (c) and (d) are y without and with filter applied to each
convolution components; It is the kth sampled point of the orbit.

53

examined by utilizing the same periodic orbits displayed in Figures 3.4—3.9. Listed
in Table 3.6 are the identified parameters. Coefficient B21 displays a larger error
because of the strong non-linearity of the related term and the relatively larger noise
contamination level (due to 01’s small oscillation amplitude and the sensitivity of

Optical encoders).

Table 3.6. Optimized identification with no digital filter applied.

 

 

 

 

 

 

 

 

 

 

 

Parameters Identified values True values Errorx 100%
Bu 0.1125 0.1131 0.5%
812 0.0790 0.0711 11.1%
815 0.0825 0.0820 0.6%
B21 1.1611 1.6816 31.0%
B22 0.2553 0.2630 3.0%
Bg5 0.2564 0.2913 12.0%
011 0.0003 1.02 x 10‘3 —

Clg -0.0002 0.00366 -
Cg 0.0164 0.0485 —

 

 

 

 

 

3.4.5 Digital differentiation and error reduction
Recurrence tolerance

Due to the limited length of the experimental data, not many periodic orbits were
extracted by setting the extraction error tolerance small, e. g. less than 5% in the
present experiment. Thus, one would naturally tend to increase the error tolerance
such that more periodic orbits can be extracted. However, with the increased error
tolerance and hence more, but less accurate, periodic orbits, the identification results
turned out to get worse for the double pendulum system. The observation is consistent

with [55]. Table 3.7 lists the identified parameters (Fourier series method) with 8%

54

error tolerance of extraction (in this case, 58 different orbits extracted), and showed
large errors. For the present system, the coefficients Bu and B21 are parameters
of strong nonlinear terms, and therefore, identification of these two parameters are
usually less stable and more prone to error due to the fact that the strong nonlinear
terms can magnify high frequency noise and the recurrence error is a major source of

high frequency noise.

Table 3.7. Comparison of calculated values by Fourier Series (FS) method and Digital
Differentiation (DD) method when error tolerance is set as 8% with filtering and
optimization process applied.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Parameters FS ID FS Errorx100% DD ID DD Errorx100%

B11 0.0126 88.5% 0.1224 8.2%
312 0.0730 2.7% 0.0760 6.9%
315 0.0908 10.7% 0.0908 10.7%
B21 1.4703 12.6% 1.6103 4.3%
322 0.2607 0.9% 0.2559 3.7%
325 0.2851 2.1% 0.2772 4.8%
011 0.0014 — -00003 —

012 0.0005 — 0.0016 —

C2 0.0053 — 0.0093 —

 

 

The large recurrence error resulted from the error tolerance of the UPO extraction.
Figure 3.13 (a) and (b) display the calculated angular velocity and acceleration curves
of a period-9 (Figure 3.5) orbit by means of the Fourier series expansion method. An
impulse caused by the recurrence error occurs in the periodic velocity and acceleration
curves, which is not true for the real periodic orbits. If we look at the velocity curve,

it can be expressed as

A a

6(t) = 9(1) + 66(1—16) + 13(1), (3.28)

where 6(t) is the calculated velocity, 0(t) is the real velocity, ad (t — tc) is the impulse

55

with amplitude a proportional to the recurrence error, to is time delay, and 17(t) is
noise other than recurrence, which is considered to be small due to low-pass filtering.

Thus, after digital filtering, and applying the Fourier transform, we obtain

2(9) 2* 77(6) + aexp(— 3101.)),

where a is the expression of the impulse function in frequency domain and is a white
noise. Furthermore, this noise contaminates all sub—harmonics of the velocity curve,
which could not be eliminated by the low pass digital filter and could hence generate

large error in the identification result.

Digital differentiation

Through the previous analysis, for a large recurrence error, one would naturally con-
sider that a similar scenario would happen to other dynamic systems with strong
non—linearity. In this case, digital differentiation could be applied given adequate
sampling points per cycle, 0. g. a high sampling rate fs. In this experiment, there
were 100 points per excitation period. Hence, a five-point differentiation algorithm
was applied to obtain the derivatives and double derivatives of angular displacements:

= 8if(:v + h) — f0: — hi] - f(22 + 210+ for - 2h)
12h

 

35(5) + 0014) (3.29)

and

: 16[f(:r+h)+f(:r—h) —2f(:r:)] ——f(:r+2h) ——f(:1:—2h) +2f(:1:)

12712 + 0(h4).

(3.30)

fix)

 

The errors of these algorithms can be reduced with a smaller time interval h. The
five-point algorithm can also reduce the influence of high frequency noise since it

involves the neighboring four points. The obtained curves are in Figure 3.13 (c) and

56

 

   

 

 

 

 

 

 

 

 

 

 

 

 

1 0.5
3; 0.5 g
a ‘6
2 5.
a)
5 o g
o. 5 A A 0. 5
o 200 400 600 800 1000 0 200 400 600 800 1000
k k
12 05
1
0.8 A
1‘ as
>5 0.6 ‘6
8 g 0
‘6 0.4 *6
5 C
02 “
0
0.2 o. 5
0 200 400 600 800 1000 0 200 400 600 800 1000
k k
(C) (d)

Figure 3.13. Obtained 0.1 and 01 orbits by Fourier series method (a) and (b) by digital
differentiation method (c) and (d); the recurrence impulses occur at k=800 in (a) and

(b).

(d) for the extracted period—9 orbit. It is apparent that the recurrence impulses are
eliminated for both of the cases, and much less high frequency noise is displayed in the
velocity and acceleration curves. The corresponding identified parameters are listed
in Table 3.7, which is more precise compared to the result of Fourier series method.
However, the algorithm for digital differentiation itself introduces calculation er-
rors in equation (3.29) and (3.30). It could not predict precisely the small value
parameters, e. g. friction parameters, according to an identification test based on the

simulated noise-free double pendulum system. On the other hand, the Fourier series

57

method can identify all of the parameters with satisfying accuracy if error tolerance
of extraction is small enough (2% in this case). It shows that the digital differenti-
ation method is more stable, but not more accurate compared to the Fourier series

expansion algorithm.

Choice of sub-harmonics or harmonics

The identification result varied when a different choice of sub-harmonics was applied.
The problem was not so troublesome in previous applications of the harmonic balance
method [24, 25] where nonlinearity was simple and not as strong as the double pen-
dulum case. However, in the present experiment, different sub-harmonic sets led to
quite different estimated parameters. It is then necessary for us to seek some general
rules for judging the result.

In equation (3.18) and (3.20), the first M sub-harmonics of each UPO (sub-
harmonics are functions of sin (ix/k) or cos (ix/k) where 2' and k represent the ith
term in Fourier series of a period I: orbit. Here, we denote M k as the truncated sub-
harmonics for a period-k orbit. The choice of Mk remains an issue. For a periodic
orbit whose period is a multiple k of the excitation period, if Mk 3 k, i. e. including
only sub-harmonic frequencies less than or equal to the driving frequency, the result
was tested to be the best for the examined pendulum system. One reason is that
these sub-harmonics consisted of a large part of the energy of the displacement sig-
nal, and hence, contained a relatively small portion of the noise contamination (see
Figure 3.11). Also, the noise components in the velocity and acceleration signals were
reduced according to equations (3.7—3.10).

The remaining question is whether it is possible to measure the identification
error and use it as an indicator of how ‘true’ the identification is. To quantify the

identification error, we refer to linear regression techniques, borrow some concepts in

58

 

residue

 

 

l l l

5 .1
-10 -5 0 5 10 15 20
predicted values of theta2 acoel.

 

Figure 3.14. Identification residual é'g with each dot representing a subharmonic; the
horizontal position of each dot is the predicted value of (72, the vertical position of each
point is the identification residual; identification was done with low-pass filtering, the
Fourier series expansion and the sub-harmonic set whose frequency 5 the excitation
frequency.

statistics, and transform (3.18)-(3.22) into
6, = A,§:’, —-(7,-,z' = 1,2, (3.31)

where é,- is the residual vector. Then, we can define the identification error 5,- as

_ lleitlloo

i— 7 i
ll‘h‘lloo

where (7,- : 14,-3f, is the predicted vector of (7,. With the identification error defined,

the rule of thumb for judging a good identification is 5,- < Q}, where 6c is a positive

59

critical value. cc = 10% was used in the experiment. Figure 3.14 displays the residuals
52 and the corresponding identification error 52 is 18.4% when all 140 sub-harmonics
were included in the identification, whose frequencies were less than or equal to the
driving frequency. Apparently, the results for sz for j=1, 2, 5, are not satisfying,
and the comparison in Table 3.8 and Table 3.3 also corroborates the rule of thumb
since B21 has a 11.2% error. On the other hand, the Blj parameters have 51 = 8.8%,

and are quite consistent with the result after optimization in Table 3.3.

Table 3.8. Identified parameters with the Fourier series method and low-pass filtering,
but without optimization to sub-harmonics set.

 

 

 

 

 

 

 

 

 

 

 

Parameters Identified 'D‘ue values Errorx100%
Bu 0.1178 0.1131 4.2%
812 0.0782 0.0711 10.0%
815 0.0861 0.0820 5.0%
321 1.4935 1.6816 11.2%
322 0.2562 0.2630 0.5%
325 0.2828 0.2913 3.4%
011 0.0001 1.02 x 10-3 —

012 0.0007 0.00366 —
C2 0.0147 0.0485 —

 

 

 

 

 

Our problem is now how to optimize the identification process so as to improve the
accuracy, i. e. minimize the identification error. From the statistics point of view, we
have a good linear regression curve if the resulting residuals are distributed evenly and
randomly around the predicted values. For this purpose, an optimization algorithm
was developed to exclude the sub-harmonics terms which result in large residuals, and
retain the terms which consist of most of the sub-harmonics and should have small
residuals. The sub-harmonics having large residuals are very likely contaminated

heavily by noise. The optimization process involves the following steps:

60

1. Given the set of sub-harmonics, do the identification process and find the max-

imum absolute residual value emax.

2. For a level of significance 0, which is a small value, remove from the sub-

harmonic set those sub-harmonics whose corresponding residual e > (1 —

fi)ema:rc

3. Repeat the first step by using the remaining sub-harmonics, and compute the

identification error 5,.

4. If a, < 10%, then stop the optimization process and assume that the desired
result has been obtained; if not, go back to step 1 with the remaining sub-

harmonics set and repeat the optimization process.

For the investigated system, after three optimization processes with 0 = 5% (by
excluding 10 erroneous sub-harmonics), the identification errors were reduced to 81 =
8.1% and 52 = 9.4%. The corresponding results are listed in Table 3.3. Displayed
in Figure 3.15 is the residual distribution of e'é after optimization. Compared to the
identified values without optimization in Table 3.8 and Figure 3.14, the Optimized ones
have smaller errors, smaller residuals, and are more accurate. However, the proposed
optimization can not work well for all cases. After a few cycles of optimization,
if the identification errors are still undesirable, we should either use more precise
UPOS, or re-select the set of sub—harmonics before optimization such that the noise
contamination can be minimized. In order to obtain more precise UPOS, besides the
UPO extraction with smaller tolerance of error, we can also refine the extracted UPOS
[55]. But the refinement process may generate erroneous results for quasi-periodic
orbits and may not be adequate for limited data sets. For better selection of sub-
harmonics, one may choose only those sub-harmonics that have the largest amplitude

in the UPO acceleration FFT spectrum (also avoid those noise contaminated high

61

frequency harmonics, since it could only be noise). Furthermore, the selection could
be simplified by choosing the harmonics instead of sub-harmonics if most of the orbits

are composed mainly of the harmonics of the driving frequency fe.

 

resndue
O
l
\
\
‘
‘D
'-
a"
Ofl
\
L

 

 

1 l l

5 l
-10 -5 0 5 10 15 20
predicted values of theta2 accel.

 

Figure 3.15. Identification residual 62 after optimization with each dot representing
a subharmonic.

3.4.6 Validation methods

Generally, most of the parameters of a nonlinear system are unknown to us. The
direct comparison discussed in previous sections is not available for most applications.
Besides, lots of the ‘true’ parameters in Table 3.1 were also estimated. Hence, we do

not know the exact errors for the identified values. Two other methods were applied

62

here to verify the identified parameters and the effectiveness of the identification
algorithm. The first one is to verify our method by identifying the simulated double
pendulum system such that comparisons can be made based on identification results
and the phase portraits of the experimental and simulated systems. The second

method involves the linearized system properties, e. g. natural frequencies.

Identification of the simulated double pendulum system

To verify the effectiveness of the identification process, a simulated system was also
examined with parameters set as the identified values listed in Table 3.3 except for
the friction terms (see Section 3.4.3). The error tolerance of extracted unstable peri-
odic orbits is 5%. The phase portrait has been shown in Figure 3.10. Some similarity
observed in the simulated system compared to the experimental phase portrait (dis-
played in Figure 3.3), e. g. chaotic behavior. Many detailed chaotic characteristics
were not available due to the inadequate experiment data. However, it turned out
that more and different UPOS could be extracted from the equally large, sampled
data of the simulated system, than from the experimental data set. The difference
can be explained by the sensitivity of the chaotic systems to parameter settings and
initial conditions. Meanwhile, for the simulated system, the comparison in Table 3.9
between the identified values and the parameter settings shows that all the param-
eters including the friction coefficients are identified correctly, which confirms the
effectiveness of this algorithm. The friction coefficients, probably due to their much
smaller values and the weakness of the least mean square method, are still hard to
calculate very accurately, and thus, have larger error percentages. With little noise in
the simulated data, the error can only come from the recurrence error of the extracted
periodic orbits. Also, it confirms the difficulty of identifying small friction parameters

in the experiment, since more noise contamination occurred in the experimental data.

63

Table 3.9. Comparison of the identified values and the true parametric settings of the
simulated system with low-pass filtering and optimization applied.

 

 

 

 

 

 

 

 

 

 

 

 

 

— Identified Parameters setting Errorx100%
J01 (kg.m2 5.767 x 10‘4 5.757 x 10—4 0.2%
J02 (kg.m2 4.310 x 10“5 4.228 x 10‘5 1.0%
B11 0.1166 0.1167 0.1%
812 0.0787 0.0781 0.8%
315 0.0848 0.0849 0.1%
011 0.0005 0.0 -
012 0.0038 0.00366 3.8%
B21 1.5943 1.6149 1.3%
B22 0.2568 0.2562 0.2%
825 0.2866 0.2868 0.1%
02 0.0429 0.0485 11.5%

 

 

 

 

 

Linear properties

The linearized pendulum can also be applied to validate the identification results, e. g.
by comparing the natural frequencies of the linearized system. Suppose the pendulum
has only small-angle oscillations without an excitation force. By discarding the higher

order terms and neglecting dry friction, equation (3.4) can be simplified as

<21” + 311452" + 312951 - C12(¢2' - (61’) = 0
(22" + 321¢1”+ 322951 — 02(¢2' - 451') = 0

(3.32)

Since our goal is to examine the natural frequencies, by neglecting the damping terms,

equation (3.32) can be further simplified to the form of

<21" = —Bl2/(1 — 311321l¢1+ B22311/(1— 811321)” (3 33)

$2" = B12321/(1— B11321)¢>1+ —B22/(1 - 311320922

64

and the characteristic matrix of equation (3.33) is

( 0 1 0 0)
—B12/(1—B B) 0 B B /(1—B B )0
A = 11 21 22 11 11 21 . (3.34)
0 0 0 1

 

 

(312321/(1-311321) 0 -B22/(1-B11321) 0)

The eigenvalues of matrix A are the natural frequencies of the linearized system
in non-dimensional form. The natural frequencies of the identified system were found
to be 1.33GHz and 2.904Hz. Through FFT analysis, the natural frequencies obtained
by the experiment of free vibration are 1.25Hz and 3.00Hz. Comparison in Table
3.10 shows that the natural frequencies match the FFT result from experiment and
the ’true’ values calculated by the settings of the double pendulum. It suggests that
the identified parameters, excluding the friction terms, are reliable for the purpose of

system linearization.

Table 3.10. Natural frequencies.

 

— Identified Experimental Error x 100% ‘true’
fn1(hz) 1.336 1.25 6.9% 1.292
fn2(hz) 2.904 3.00 3.2% 2.940

 

 

 

 

 

 

 

 

 

3.5 The energy balance method for identification
of the gamping coefficients

In chapters 2 and 3, the identification results were satisfying except for the small

friction coefficients. The inaccuracy of friction coefficients became worse for noisy

65

experiment data. The reason is due to solution method, i. e. least mean square
method, which tends to give more inaccurate solution for smaller parameters, e. g.
the small friction coefficients in this study. Also, it showed that the errors in friction
terms, since they are much smaller than other parameters, have little influence on
other estimated values. Consequently, we apply an energy balance method to estimate
the damping coefficients.

Suppose we have a single (1. o. f. system with a governing differential equation

= -1... — can — re) — fp($)P(t) + 90). (3.35)

where k, c, f (.23) and fp(:r) are unknown, and c is a small damping coefficient. First, we
can apply harmonic balance method to identify the parameters. Multiplying equation
(3.35) by cos 75-1”!- and integrating it over one cycle of a periodic orbit, it can then be

expressed as

ftt+T(:iE + kw) cos %dt =

(3.36)
—- ftt+T{c:i: + f(:z:)} cos #8: — ftt+T{fp(a:)p(t) + g(t)} cos grim.

Similarly, we can also obtain equations for sin 1%: functions. The unknown functions
f (11:) and fp(:1:) can be approximated by function series with unknown coefficients (see
Chapter 1), and the integrals become the corresponding sub-harmonics in the FFT
spectrum. Then, the equation (3.36) can be converted to a matrix equation similar
to equation (2.15). One of the advantages of the harmonic balance method is that it
can identify all the parameters simultaneously. However, smaller parameters are also
vulnerable to noise influence from other terms.

To eliminate the noise from other terms, with the parameters partially known
from the harmonic balance identification, we can refer to energy balance method [56]

for a second-step identification. For some periodic orbit D of (3.35), we want to find

66

the energy generated or dissipated by each term in the equation over one cycle, and

we obtain

f0 idx = —kjg:rd:c — ch idz — fD f(a:)dx — .710 fp(:r)p(t)d:l: + fD g(t)d:i:. (3.37)

The 56', km and f (11:) terms represents the conservative terms. Hence, theoretically, the
work done by the conservative forces within one cycle is zero. Equation (3.37) then

becomes

ch :irdx = — ij fp(:r)p(t)d:r + f1) g(t)d:r, (3.38)

where the left hand side of the equation is the energy dissipated by the damping
term, the right hand side is the energy input into the system by force and parametric
excitation. The integrals in (3.38) can be calculated. Thus, if we can obtain multiple
periodic orbits of the system, we can then solve the damping coefficients by the least
mean square method. Since we have obtained the quite accurate estimation of those
unknown functions in the harmonic balance identification, the damping coefficient

can then be identified more accurately by the energy balance method.

Table 3.11. Identified damping coefficients by the energy balance method.

 

 

 

 

 

— Identified Experimental Error x 100%
— C 11 012 C2
The identified -0.0004 0.0045 0.0396
True values 0.000101 0.00366 0.0485
error — 23% 18%

 

 

 

 

 

 

For the single pendulum system in Chapter 1, the identified damping coefficient
is cr = 0.0353, and the corresponding error is 2.1% which is much smaller than

result error of the harmonic balance identification. For multi—d. o. f. systems, we can

67

apply a similar method to each differential equation, and obtain the more accurate
estimation of the damping terms. Table 3.11 listed the identified damping terms by
energy balance method. Apparently, the second-step identification obtained much
better results. But the error is still much larger than 10%. The reasons could be
noise contamination in the signals, or inaccurate estimation of the ‘true’ values which

were also estimated by the small amplitude free vibration.

3.6 Concluding remarks

A double pendulum experiment was examined for chaotic system identification. The
investigated system was a multi-degree-of-freedom system with strong non-linearity,
e. g. mixed (0,, ¢,' and ¢,-’ non-linearity with whirling in 51 space. However, only
displacement signals were directly measured. To adapt to these new challenges, some

modifications were added to the harmonic balance identification algorithm:

1. The identification appeared to be less noise-resistant in this case, mainly due
to the strong nonlinear term of ((1),!)2 sin ((02 — (01) and 05,-” cos (0)2 - (1)1). The
high frequency noise contaminated the strong nonlinear terms at lower frequen-

cies without adequate low pass filtering of each component before convolution.

2. A digital differentiation algorithm was developed and applied to the experimen-
tal data in order to make the identification results more robust even with large
recurrence errors in the extracted orbits. It could be of use for limited data
sets. However, the digital differentiation algorithm also introduced differentia—

tion error, and thus, did not give accurate values of friction terms.

3. Choices of sub-harmonic terms also had influence on the identified parameters.
Inappropriate selection of sub-harmonics can generate poor results. To avoid

poor results, the key factor was to avoid noise contaminated sub-harmonics. For

68

the present system, the sub-harmonics of frequencies less than the excitation

frequency were selected to avoid noise.

4. Linear regression techniques were applied to quantify the identification error 6,,
which reflected the error of the identified parameters by examining the residues
and the predicted values. Based upon the identification error, an optimization
algorithm was proposed to improve the result. An identification error less than
10% indicated an rather satisfying result. However, optimization is limited by
its statistical properties, and can not work for all data to satisfy the rule of

thumb.

Friction was a problem in the identification process. For slightly damped systems,
since the friction factors were much smaller than other parameters, the identification
could not produce accurate values of the friction terms. Noise and recurrence error
were the two factors that contributed to this error. However, the harmonic balance
algorithm is robust, and the validation process showed that the friction error had
little effect on other identified parameters. Based upon the identified non-damping
parameters, the damping coefficients were identified again in the second step of the
identification—the energy balance method which exploits the energy balance and dis-
sipation of the vibration system. The results were tested to be improved. Meanwhile,
the convergence of the identification algorithm was not fully discussed in the exper-
imental study, mainly due to the direct comparison between the identified and the
measured parametersln lien of convergence, the model was verified using residuals,
comparisons to known parameters and linearized properties. Through the experi-
ment, it can be concluded that the examined identification method can be applied
to systems with chaotic response, strong non-linearity and multiple degrees of free-
dom. With adequate modification, the identification result could be improved, and

the quality of the result can be quantified by the identification error.

69

CHAPTER 4

Simplification to Frequency
Domain Method

4.1 System identification by frequency method

Pure frequency based parametric identification method has long been limited to non-
chaotic systems. The reason is that, generally, frequency domain methods are based
upon frequency response functions, or the harmonic balance method. For linear sys-
tems, their responses to any kind of dynamical load (which can be expressed as the
sum of simple loads, e. g. harmonic loads), by superposition, is equal to the sum of the
responses of each simple load acting separately. Thus, for linear systems, the harmonic
balance identification is equivalent to using frequency response functions (FRFs) if
multiple frequencies of excitation were used. The FRF is perfect for linear system
identification, since the FRF can be easily generated by an impact response. For
weak nonlinearity, Richard [12] and Kerschen [11] developed the conditional reverse
path (CRP) method based upon the FRFs, which seperates the response into linear
and nonlinear response parts. In [18, 19, 20, 21, 22], the harmonic balance method
was proposed to identify non-chaotic nonlinear systems, because it requires steady

state responses, i. e. periodic orbits. Similarly, Plakhtienko [16, 17] and Ghanem [23]

70

developed respectively the special weight method and the wavelet analysis method.
All of the three methods have a common feature: a basis of orthonormal functions.
Furthermore, they all require a periodic response, though it is not an obvious require-
ment in the latter two methods. Discrete signal reconstruction by the orthonormal
basis requires periodic orbits. It inspired Yuan and Feeny [24, 25] to extract unstable
periodic orbits (UPOS) from data and then apply the harmonic balance method for
identification. This method was examined by whirling and experimental systems in
previous chapters. It can be regarded as a hybrid time-frequency (HTF) method.
The idea comes from the knowledge that hyperbolic chaotic systems, though having
no steady state periodic orbits, have chaotic responses with dense sets of unstable
periodic orbits inside.

In previous research [24, 29] and previous chapters, the HTF method has been
tested successfully on various chaotic systems, either by simulation or by experiment.
Although the method is capable of being applied to any chaotic system, the process
could be time-consuming since it requires extracting UPOS from the chaotic data. For
better results, it requires more accurately extracted UPOS and more orbits. Therefore,
the method requires a large amount of the collected data. Furthermore, for large d.
o. f. systems, recurrences are less frequent, meaning UPOS of low periodicity are
visited less often. These drawbacks make it hard or impossible for fast parametric
identification, or real-time identification. Hence, one would like to investigate if it is
possible to skip the process of extracting unstable periodic orbits by applying more
chaotic properties, such that the harmonic balance method can be more conveniently

applied to nonlinear/ chaotic systems.

71

4.2 More chaotic properties and the improved fre-
quency domain method

The UPOS are extracted by satisfying the recurrence error tolerance,
IISUC + nT) - 509)” < 6, (4.1)

where S(k) = [ mag) $094.77,” $0; +de)] is the displacement signal in
an m dimensional embedding space [52], T is the excitation period, Td is the time
delay for embedding dimension, and 6 is the tolerance of recurrence error. Hence,
we say a close approximation of the UPO is found when equation (4.1) is satisfied.
The extraction method is usually targeted for short term periodic orbits. For long
periodic orbits (n large), the extraction can be much simpler.

As we know that the set of unstable periodic orbits is the dense skeleton of the
hyperbolic chaotic attractor, it actually implies that for any chaotic orbit (of some
hyperbolic chaotic attractor with UPOS inside), a UPO can be found in any neigh-
borhood of the chaotic orbit, and the extracted UPOS can be utilized for chaos in-
vestigation [57, 53, 58, 59, 60]. Here, we will use this property and also add a claim
that:

Given a bounded hyperbolic chaotic set, assuming unstable periodic orbits inside
the attractor, any chaotic orbit within some time interval nT, where n is a large
integer, can be ‘approzimated’ to some extent by an unstable periodic orbit with nT
periodicity with a similar frequency spectrum. The similarity increases as it gets
larger.

The word ‘approximate’ here generally can not be examined by the rule in (4.1) for
extracting UPOS, rather by the frequency spectra: the end-points may have significant

error, but the most of the orbit is a close approximation to a period-n orbit.

72

 

 

 

 

 

1 I I l T
0.9 - 4
0.8 — .
0.7 — .
0.6 - .
E 0 5 ~ ,

:5

0.4 - ]
0.3 — .
0.2 - -
0.1 ~ -

0 l 1 1 4L
0 0.2 0.4 0.6 0.8 1

X0‘)

Figure 4.1. The left shift one dimensional map

The claim comes from the fact that every chaotic set can be converted to equivalent
symbolic dynamics [52, 61]. With symbolic dynamics, the hidden properties of chaos
emerge clearer before our eyes. For instance, we can take a look at a one dimensional

map (shown in Figure 4.1) expressed by

22:, 0 S a: < 0.5
:r(k+1) = . (4.2)
2:1: - 1, 0.5 S :r < 1

This map is conjugate to the left shift operation of binary numbers in the interval

[0, 1) in symbol dynamics,

$(k +1) = :i:(k) x [10]), mod 1, (4.3)

73

where the subscript b represents binary numbers. The uncountable set of irrational
numbers in [0, 1) forms the chaotic attractor. On the other hand, the countable set
of rational numbers in [0, 1) form the set of the unstable periodic orbits, and it is
dense in [0, 1). Since the set of irrational numbers is a set of Lebesgue measure 1, the
probability of observing a unstable periodic orbit is zero. However, for any chaotic
orbit :r(t; 1‘0) of length n starting from an irrational number 1:0, there is always some
periodic orbit y(t; yo) with periodicity 71. starting from rational number yo, such that
”yo - 17:0” < 2‘", which implies the first n binary digits of $0 and yo are identical.
Comparing :r(t; $0) and y(t; yo), we observe that, for the first it points of these two
orbits,

|1'(k; x0) — y(k; yoll < 2—<"—k>.vxc s n. (4.4)

The error between the chaotic orbit and the periodic orbit increases as t increases.
For t = n, since y(t;y0) is a period n orbit, y(ngyo) = yo z :00, whereas, :r(n;:r:0)
becomes some number other than 320, and thus, forms a ‘recurrence’ error. For the
extraction of short periodic orbits, we can use the rule of (4.1) by choosing chaotic
orbits with small recurrence errors. However, the extraction of UPOS is harder for
high (1. o. f. systems, and can only utilize a small part of the data for identification.

Indeed, it has been observed by Yuan [62] that long chaotic orbits were good
for parametric identification by the harmonic balance method, though they did not
satisfy the UPO extraction rule. It then raises the question whether the long chaotic
orbits are usable for UPO approximation. For any long orbit :6, i. e. large n number,
given an error level 6, according to inequality (4.4), there exists a K, such that for
the first K points of that chaotic orbit, the error between the chaotic orbit and some

period-n orbit satisfies

||:L‘(k;1:0)— y(k; yo)” < 2-("-—k> < wk 3 K < n.

74

It implies that the first K points of a sampled chaotic orbit :1:(t; 2:0) of length n are
very close to the corresponding points of some unstable periodic orbit y(t; yo), and

the rest of a: may deviate from y. Hence, we can express a:(t) by

1705:1130) = y(t; 90) + 610) + 6205), (4-5)

I(tn‘ )-y(t;y ), ts K
where el and e2 are error functions, and e1(t) = O 0 ,
0, t > K

0, t S K
62(t) = . Function cl represents the first K points

1:03:20) - y(t;yo), t > K
of the chaotic orbit where the approximation error is small and [e1(t)| < 6. Function

e2 represents the tail part of the chaotic orbit where the error to the periodic orbit
tends to increase beyond the tolerance of error 6, but is still limited [82(t)] < 1. In

the frequency domain, the discrepancy between :1: and y is

NW; $0) - y(t; 310)} = 70910)} + jr{f«’2(t)}- (4-6)

Obviously, if we can manage to reduce the error caused by the tail part, i. e.
(n — K) << n, the majority of the chaotic orbit can be a good representation of some
periodic orbit. As for the tail part, in the frequency domain, the recurrence error
generated by e2(t) can also be neglected if (n — K) < n. Hence, in the frequency
domain, the chaotic orbit :r(t;;ro) becomes a close approximation to y(t;y0). The
solution to this problem is simply a large n, i. e. a long chaotic orbit. For instance,
suppose the tolerance of error 6:0.01, and n=100. There are only six tail points
out of the length-100 orbit with errors greater than 6, and in frequency domain, the
amplitude of the recurrence error e2 is less than 0.02, which is acceptable. Therefore,
every long chaotic orbit in the left-shift map is a good approximation of some unstable

periodic orbit, and the approximation is better if n is larger.

75

This simple example based upon the shift map on the unit interval can be gen-
eralized since all the chaotic systems are equivalent to some corresponding symbol
dynamics. Thus, the claim is supported.

The claim simplifies the investigated systems’ identification method, and turns it
into a pure frequency method. Assuming a single chaotic data {S(k)} was sampled,
the chaotic data itself is a very long approximated unstable periodic orbit in the
frequency domain, and can be used directly for identification purposes. To acquire
more UPOS within the limited data and improve the identification results, we can
divide the orbit into pieces by some interval 7' = nT. For each piece, we can apply
the harmonic balance method and get the same equations as those in (2.15). However,
without inequality (4.1) satisfied, the recurrence error can cause jump noise in the
frequency spectrum, and thus generate inaccuracy for identification process. This
shortcoming can be overcome by choosing nT large enough such that the jump noise
influence to the whole power spectrum is reduced to an acceptable level.

It is now clear that a pure frequency method is justified. The pure frequency
method for chaos is intimately tied to our understanding of unstable periodic orbits

and the hyperbolic chaotic set.

4.3 Calculation examples

We will use two examples here to examine the improved method. The first one is
a simulated single pendulum with horizontal base excitation. The second example
is a double pendulum experiment. In both of the systems, a linear regression based

optimization method was applied to improve the identification results.

76

4.3.1 Base-excited single pendulum

The simulation was based upon a horizontally excited single pendulum. Its non-

dimensional differential equation is
6+ 2€/rfl + 1/r2 sin0 — fsintcosd = 0, (4.7)

where t = wr, w is the angular excitation velocity, and r is the actual time; 6 is
the viscous damping coefficient; r = w/wN, (0 N = mge/ J is the natural frequency
of the linearized system, m is the pendulum mass, e is the the centroid offset from
the pendulum hinge, J is the moment of inertia about the hinge point. Meanwhile
coefficient f = mea/ J, and a is the excitation amplitude. For simplicity, we denote
cr = 25 / r as the new non-dimensional friction coefficient. The function f sintcos 0 is
the nonlinear parametric excitation term, and 1/r2 sin0 is the autonomous nonlinear
part. In the simulated system, or = 0.036, f = 1.52 were used. Assuming viscous

damping and displacement nonlinearity, the differential equation is

crab + fanl(:r) + fnnl(:c) cost + fpnl(a:) sint = —:ii, (4.8)

where fan] is the unknown autonomous part, fnnl and fpnl are the unknown para-

metric excitation parts. The unknown functions can be approximated by

N

fanl a: Z Qi¢i($)a (4.9)
i=1
N

fnnl z 2 ni¢i($)a (4.10)
i=1
N

fpnl z E Pi¢i($)a (4.11)
i=1

77

where ¢,(:L') is a linear interpolation function:

W, (k—1)d<a:<kd
iiilg—d—‘E, kdgx<(k+1)d-

0, otherwise

(MSU)

 

\

Hence, the differential equation is to be identified in the form of

N N N
crzi; + k2 + Z q,¢,-(:r) + Z n,(,b.,-(2:) cost + Z p,d),-(:1:) sint = —:'13. (4.12)

Meanwhile, since angular displacement is a S1 signal, we used 2 = exp (id) to solve
for 0 since it is much easier than the original method (see 2.2). Two sets of data were
examined: the noise-free data and the data with 10% random noise which was a noise
source with uniform distribution and the amplitude equal to 10% of the amplitude
of the signal. Compared to noise with normal distribution, the uniformly distributed
noise has more high frequency components and easier for computer generation. The
long data was cut into equal-length pieces with r = 60T. The 1' value can be selected
randomly as long as the periodicity is large enough for UPO approximation (see
previous section for error estimation) such that the recurrence impulse is small enough
in the frequency domain. Meanwhile, we should also avoid large 7' such that adequate
number of UPOS can be obtained for identification purposes. For the noise free data,
the identified cr = 0.0352, and the unknown functions were displayed in Figure 4.2.
Compared to the HTF method in Chapter 2, the error is even smaller. Figure 4.3 also
compares the Poincare’ sections of the original system and the identified system. The
chaotic property of the two systems is quite similar. For the data with 10% random
noise, the identification algorithm shows strong noise-resistance. The estimated or is

0.0343, and the identified stiffness functions are displayed in Figure 4.4.

78

0.5' 2.

 

 

 

 

 

 

"c‘ 0* *+*—.—*—*—*—*— 3f 0
-1 .
_05 A . . . _2 . . i .
—4 -2 O 2 4 -4 -2 0 2 4
theta position theta position
(a) (b)
2 .
1 .
5 0
-1 .
_2 l . 4 1
-4 -2 0 2 4
theta position

Figure 4.2. Identified functions by interpolation functions on noise-free data, (a) nk,

4.3.2 Experimental double pendulum

The second example is a double pendulum system with vertical base excitation. It
is both a whirling and a multi-degree-of-freedom system. A 100T orbit length was
selected. During identification process, the high frequency noise was filtered out by
a digital filtering technique, and linear-regression-based optimization was applied for
a better result. The identified parameters are listed in Table 4.1 in comparison with

the ‘true’ values (the values were measured by other methods, and thus, may not be

79

 

 

theta(t’rdt)
O

theta(H-dt)
O

l
M
l
N

   

 

 

 

 

 

 

_4 -4
4 0 0
theta(t) theta(t)

(a) (b)

Figure 4.3. Comparison of Poincaré sections, (a) the original system, (b)the identified
system

the exact parameter values, see Chapter 3).

Shown in Figure 4.5 is the residual plot of the first matrix equation. It can be seen
that the residuals are much less than 10% of the maximum predicted value, implying
a good identification result by the 10% rule of thumb. Actually, the small residual
values were obtained even without using the optimization process. It supports the

claim that any long chaotic orbit is close to some unstable periodic orbit.

4.4 Concluding remarks

In this chapter, we improved and simplified the parametric identification method
for chaotic systems [24]. Assuming hyperbolic chaos for the investigated system,
due to properties of a chaotic attractor, the set of unstable periodic orbits forms
a dense skeleton of the attractor. Also, a long chaotic segment is arbitrarily close
to a UPO, except at the tail or head part, and is approximately the same in its

frequency spectrum. Thus, balancing harmonics of an arbitrary long segment of

80

0.5' 2.

 

 

 

 

 

 

-1
-0.5 ‘ + * ‘ -2 * ‘ 1 ‘
-4 -2 O 2 4 -4 —2 0 2 4
theta position theta position
(a) (b)
2 .
1 .
5 0*
-1 >
..2 . A - 1
-4 -2 O 2 4
theta position

(C)

Figure 4.4. Identified functions by interpolation functions on 10% noise contaminated
data, (a) 71),, (b) 19),, (c) qk.

chaotic data simplifies the identification method from a complicated one involving
harmonic balance method and extraction of UPOS to a pure frequency method applied
to a long chaotic segment. Simulation and experiment were examined to verify this
improvement. The method was seen to be quite noise resistant.

With this simplification, previous research on whirling systems, and the optimiza-
tion of the identification by refining the set of sub-harmonics, the harmonic balance
method has been tested to be a very powerful and simple tool for the identification

of strongly nonlinear systems. The new discovery also unites the identification pro-

81

 

 

 

 

5 T I I l I
4 *— -1
3 r a
2 r- .4
9 1~ .
E 1
8
o: 0 r ' ‘
-1 _ ' ' a
-2 - -
-3 _ "i
_4 1 1 1 1 1
“0.03 '0.02 ”0.01 0 0.01 0.02 0.03

Predicted value of theta1 acceleration

Figure 4.5. The residual plot of the first matrix equation of the double pendulum
system

82

Table 4.1. Comparison of the true values and the identified parameters by new

frequency method.

 

 

 

 

 

 

 

 

 

 

 

Parameters Identified values True values Errorx100%
311 0.1125 0.1131 7.8%
812 0.0770 0.0711 8.3%
815 0.0821 0.0820 0.1%
321 1.6091 1.6816 4.3%
822 0.2539 0.2630 4.5%
825 0.2839 0.2913 3.5%
(111 0.0003 1.02 x 10"3 —

C12 -0.0010 0.00366 —
Cg 0.0077 0.0485 -

 

 

 

 

 

cess of periodically excited linear and nonlinear systems under the commonly used
frequency domain method, i. e. the harmonic balance method.

Meanwhile, the research prompts a great amount of future work. Compared to
some other previous frequency methods of nonlinear systems [12], the present method
has no limitations to strongly nonlinear or chaotic systems. However, at present, the
harmonic balance method can be applied to nonlinear systems only if under periodic
excitation. For linear systems, the present method can be simplified to frequency
response functions (FRF) method, and can be applied to free vibration (under im-
pact) and random vibration (under random excitation, e. g. ergodic random signals)
cases. For nonlinear systems, identification under impact or random excitation are,
at present, still limited to weak nonlinearity, and are not available to the harmonic
balance method.

Nevertheless, the requirement of periodic orbits could be changed in the future if
we look at the periodic response in the frequency domain. We can say linear systems
have stable frequency response functions (F RF) under different excitations, due to

the fact that their stable responses are always stable periodic orbits with no change of

83

shape. On the contrary, the nonlinear systems do not have structurally stable FRFs
under different excitations. Their FRFs are liable to change for different excitation
amplitudes, and are not defined for random excitation. The extreme representation
of the instability is chaotic systems, in which the stable responses are no longer stable
periodic orbits, but chaos with embedded UPOS. Actually, we assume that a single
free vibration after a single impact is not adequate for nonlinar system identification
unless we know exactly what type of nonlinearity that the system has. Instead, we
would refer to a series of impact vibration signals, e. g. systems under periodic impact
excitation. For randomly excited systems, the case is more complicated. We believe
that some analogue can be made between chaotic systems under periodic excitation

and nonlinear systems under random (ergodic) excitation.

84

CHAPTER 5

The Reduced Order Models of
Nonlinear Systems Based on

Proper Orthogonal Modes

5.1 Need for reduced order models

As the degree of freedom increases, the numbers of parameters increases, and the
complexity of the analysis and identification is magnified. A good way to overcome
these issues is with reduced order models. As such, with increasingly large system
order, at some point, reduced order modeling becomes mandatory.

In this chapter, we use proper orthogonal decomposition (POD) as the vehicle
for reduced order modeling. The parametric identification is then performed on the
reduced order model itself, thus marrying the POD with parametric identification.

In what follows, we review the POD: its strengths and shortcomings. We then
propose an enhancement of POD reduced order modeling that allows the identified

system to be applied to conditions beyond those for which the POD was performed.

85

5.2 Proper orthogonal decomposition introduc-
tion

Proper orthogonal decomposition, or Karhunen-Loeve decomposition, has been widely
used tool for empirical modal analysis. As a statistical method utilizing the correlation
matrix derived from a set of measurement histories, POD leads to proper orthogonal
modes (POMS), which are taken as empirical modes, from the eigenvectors of the
correlation matrix, and the corresponding energy level of each mode, known as proper
orthogonal values (POVS), from the matrix eigenvalues. For linear systems with
evenly distributed mass, POMS actually represent the linear normal modes (LNMS)
[26]. Therefore, POD can be used as modal analysis tool [27, 28, 29, 4, 30, 31]. In
[31], a weighted POD method was proposed to show the common nature between
POMS and LNMS. Meanwhile, since its introduction, POD has been widely used in
turbulence as modal reduction tool [3, 32, 33], and recently in structural dynamics
[34, 35, 36, 37, 38, 39, 40, 41, 42, 43].

The process of finding POMS is very simple. For a random field 0(2, t) defined on

a spatial domain 9, the field can be decomposed into two parts by
9(1', t) = #(IL‘) + 19(1“. t), (5-1)

where 71(1) is the mean value, and 19(r,t) is the oscillation part. The snapshot of
19(x,t) in time tk, k=1, 2, ..., N, is denoted as 19,,(2‘). The objective is to find a
function p(r) which is most similar to the set of 19,,(22) on average. It can be obtained

by solving the following integral eigenvalue problem:

[,(21.(:r>'z9k<z>>¢<z>dz = we). (5.2)

86

jn ¢2(x)d:r = 1, (5.3)

where (19k(:r)19k(z)) is the averaged cross correlation function, and the equation (5.3)
is an added normalization condition to make the solution unique.

Solving (5.2) and (5.3), we can obtain eigenfunctions (b, (:13), and corresponding
eigenvalues A,. The ¢,(;r) is the ith proper orthogonal mode, and the corresponding
A, is the proper orthogonal value of this mode under the collected data set. The
mode indices i are ordered by the values of the POVS, from largest to smallest. The
POV A, also shows how much energy lies in that mode. The total energy is defined
by 23:1 Aj, and the share of energy in each mode is A, / 23:1 Aj.

The Signal 79k(r) can now be represented in terms of ¢,(r) by

l
191,013) z X aj¢j(:v), (5-4)
j=1

where the coefficient a, = (19(x,t),¢j(1:)). The number of significant modes l is
generally determined by 23:1 )1 j greater than, for example, 99.9% of the total energy.
In practice, the data is discretised in both time and space. Thus, for m observa-

tions of an n dimensional vector :5, the ensemble of the signal becomes a matrix

$11 xlm

lex'i ml= E 2 . (5.5)

xnl - . . xnm

If :E has zero mean (:3 is the varying part of some field variable), then the covariance
matrix S is defined as

s = ~1—xxT. (5.6)
m

The POMS and POVS are now the eigenvectors and eigenvalues of the matrix S. For

high dimensional systems like turbulence, the matrix can be huge, and a direct solu-

87

tion for eigenvectors and eigenvalues is quite difficult. In this case, another method
was introduced by Sirovich [63], but the details are not to discuss in this dissertation.

With the POMS as empirical modes, we can then apply them to build reduced
order models (ROMS). The dissertation research on POD and parametric identifica-
tion is quite different from the POD method in fluid turbulence research, the latter
of which is mainly a tool for model order reduction from PDES to ODEs by POMS
and the Garlekin method [33]. The investigation in the dissertation is a further step
in building reduced order structural models by POD without necessarily requiring
knowledge of the continuous system. Assuming no, or partial, knowledge of the
governing equations of the systems, we try to reconstruct the systems by building

differential equations, based upon the empirical modes, expressed as
Ms} + or + K55 + 13.107. 5') + 5.207. one) = rut). (5.7)

where all of the coefficients are unknown, functions fnl and fnz are the nonlinear
functions, and y}(t) is the external excitation force with known excitation frequency,
y‘fit) is the parametric excitation. We will then need a parametric identification
method to find the best suitable values for the unknown coefficients. For simplicity,
the systems investigated here are free vibration systems with no damping. The reason
is that this is the initial investigation of a vast issue. The free vibration systems
without damping can also display nonlinear normal modes under bounded initial
conditions [64, 65, 2], which is a good tool for comparison with POMS. Therefore,

equation (5.7) is reduced to
M5 + K55 + fn1(a‘:’) = 0. (5.8)

The linear part of the model involves relatively few problems. However, there are

88

questions regarding nonlinearity. First, what kind of nonlinearity Should be assumed?
Even though we know the physical nonlinearity of the original system, it is not clear
what form it will have when transmitted into the new modal coordinates. Thus,
function series, such as polynomials, are needed to approximate the nonlinearity.
Polynomial functions can be functions of displacement or velocity. Generally, it is
necessary to know what kind of nonlinearity is involved in the unknown systems.
This is why we treat the reconstruction process as a ‘graybox’. Second, for the
approximation function series, any unnecessary term should be discarded based upon
a priori information. The major reason is not for simplicity, but for their influences
on the dynamical responses, such as chaos, bifurcation and other nonlinear behaviors
sensitive to parameters. Sometimes, the nonlinearity approximation needs guesses
and trials, at the present research progress to approximate the dynamical responses
of the original system. Third, for the reconstructed reduced order models, one would
like to know to what extent it can represent the original system. It is a problem
especially for systems with strong nonlinearity. The reduced order models may not
work well under different initial conditions. We will discuss these problems later.
Even with many uncertainties and unknowns, POD is popular in modeling and
modal analysis because it is simple and empirical. For any system, linear or nonlinear,
it gives empirical orthogonal modes of the dynamical behavior, which can be used to
set up empirical governing equations, or for detecting the active modes in the vibra-
tion data. It is especially useful for analysis of complex structures or fluids. Structural
dynamics and fluid turbulence models are continuous and high dimensional. For the
turbulence case, the dynamics is also chaotic [66]. Usually, the dynamical governing
differential equations of those systems are hard to obtain, or hard to analyze. Thus,
one has to refer to approximation methods for analysis, either by discretisation meth-
ods like finite element analysis, or by reduced order models, e. g. the POM method.

The finite element method is more direct and usually accurate. However, it requires

89

complex processing, and sometimes, can not give accurate simulation if the systems
are complicated, or chaotic. On the other hand, POM reduced order modeling is
empirical, as it obtains proper orthogonal modes through experimental data and re-
sults in a small number of the governing equations based upon the modal coordinates
of the empirical POM. The POM method treats a system like a gray box (known
basic dynamical mechanism, unknown global governing equations and parameters),
and it gives us functional models that can approximate the dynamical behavior of the
unknown systems. We can see that the large scale discretization and ROM methods
both have advantages and have different applications. However, in many cases, the
reduced order models are not inferior to the other discretisation methods, and those

models are much easier to simulate Since they are composed of ODEs, instead of

PDEs

5.3 Improving POD for strong nonlinear systems

5.3.1 Limitation of the POD

Proper orthogonal decomposition is basically a linear statistical method. When ap-
plied to nonlinear systems, POD has its own limitations. Unlike linear systems,
displacement-response configurations of nonlinear normal modes vary with the initial
conditions, implying that POMS of nonlinear systems vary with the initial conditions.
Therefore, the reduced order model from POM under one data set may only be suit-
able for initial conditions within a small neighborhood in the phase space of those
that generated the sampled data. Beyond this neighborhood, the approximation de-
teriorates due to nonlinearity. This can be a shortcoming of the POM-based models,
especially for systems with strong nonlinearity.

In contrast, N N Ms are varying curved surfaces instead of planes in the phase space.

90

One may consider N N M based reduced order models for more precise prediction of
system behavior. The N N M is a useful tool for theoretical modal analysis [64, 65, 2] of
nonlinear systems and model order reduction in N NM modal coordinates. However,
compared to POM model order reduction, it is complex and time consuming, and
requires clear knowledge of the governing equations of systems. For unknown systems
with strong nonlinearity, it is then necessary to find a way to improve POM based
models for better prediction and applicability under varying initial conditions, which
provides robustness of the reconstructed reduced order model. To investigate this, we
focused on two nonlinear systems: a two-mass system and a mass-pendulum system,
both of which were investigated under the initial conditions exciting only the first
nonlinear normal mode, such that Single (1. o. f. reduced-order models can be built,
instead of reconstructing the original two (1. o. f. systems. For the former system,
the study focused on the limitation of the reduced-order model based only on one
POM, and then developing a new modeling method with an added constraint. The
purpose of this example is to Show whether the added-constraint method can improve
the accuracy of the reduced-order model as well as expanding the applicability of the
model under more initial conditions. For the latter example, the study focused on
how to treat unknown nonlinearity in the new models: nonlinearity approximation
with the constraint coordinate included, or nonlinearity with the constraint excluded,
which simplifies the nonlinearity approximation since there is only one degree-of-
freedom involved for this example. By the two examples, the POM-only models and
the proposed POM added-constraint models were investigated by low dimensional
systems. Investigation on these simple systems will provide a theoretical basis for

further study on high dimensional systems.

91

5.3.2 The example of a two-mass system with nonlinear

spring

We consider a two-mass system with the differential equations

5131 = —2:l31 +1122 —:lL":13

3

(152 = 1‘1 —2:z:2

(5.9)

where m = 1, and k = 1. In [2], the NNMS and modal differential equations of

equation (5.9)were solved for small amplitude vibration.

 

30
20 -
10 ~

0 ..

dy2/dt

1
V!

-10

 

 

 

I
0.)
O

l.

-20 -10 0 10

20

Figure 5.1. Phase portraits of the original system (‘+’), the POD only model (the
dashed line), and the POD-constraint model (the dotted line) under the 5th initial

condition.

To build a reduced-order model for this system, we will only consider the first

92

 

‘L5

 

 

 

 

-10 0 10 20

Figure 5.2. First N N M displacement shapes in y1, y2 coordinates, number 1—5 corre-
spond to the five initial conditions.

mode with different initial conditions. Assuming that the system is unknown to us

except that it has cubic nonlinearity, the goal is to build a Single d.o.f. model for the

first mode, and make it applicable for general initial conditions on the first NN M. We

used five different initial conditions for the simulation. With zero initial velocities, the

five initial displacements were (1, 1.3833), (1.5,2.9250), (1.75,4.0365), (2.0,5.4667),

(3.0, 17.0). The reduced order model based on POD was done under the fifth initial

condition only.

For the two d. o. f. system, there are two POMS: p'i and p3. The transformation

matrix A can be defined as

A=Wzfi 92)

93

Hence, the modal coordinates can be defined as

37 = Af, (5.10)
y :1:
1 = AT 1 . (5.11)
92 332
For the Simulated example, we obtained 1)] = ( —0.990 0.141 ) and P2 =

( 0,141 0,990 ). The POM p3 is the dominant POM, and consists of 99.7% of the

—0.990 0.141
total energy. Hence, the matrix A = . The modal coordinates yl,

0.141 0.990
yg obtained were y1 = —.990:1:1 + 0.141132 and y2 = 0.14121 + 990232. We can then

use new coordinates to transform the original differential equations into

.9" =9 9.9
1 1(1 2) . (5.12)

92 = 9201.92)

As y2 accounts for a greater portion of the total energy, according to the original
POD, the yl would become significant, and the two (1. o. f. system could be reduced
3 single (1. o. f. system in terms of y2 in equation (5.12). The dynamical behavior of

the original two-mass system could be approximated by yz only, using
= A . (5.13)

However, if the yg mode is not sufliciently dominant, the reduced order model based
only on the linear 3);; and equation (5.13) can not predict the system behavior satis-
factorily.

Suppose the coefficients of the ODES, such as equation (5.9), are unknown. Then,

94

the transformation of the coefficients in equation (5.12) will produce complicated
expressions of the original parameters. To bypass this complication, the POM modal
projection and the parameter estimation can be coordinated as follows. The form of

the dominant equation can be expressed by

3 .
92 = (€312 + Z njyé. (5.14)
j=1

where the parameters are unknown before identification. Assuming displacement
Signals measured, a parametric identification method can then be applied. For this
example, the direct identification method [5] was used, since the simulation signals
were noise free and in free vibration. Digital differentiation was applied to obtain
velocity and acceleration Signals.

Displayed in figure 5.1 are phase portraits of the original two-mass system and the
one d.o.f. identified POM model under the fifth initial condition. A discrepancy be-
tween the two systems is visible in the behavior of yg. On the other hand, under these
initial settings, the behavior of the two-mass system is apparently one dimensional.
Meanwhile, Figure 5.2 shows the changing modal displacement shapes under these
conditions. The displacement mode Shape of the fifth condition is the most differ-
ent among them. Obviously, the modal surface is quite curved, and the POM model
(5.14) under the fifth condition can not be applied to other initial settings. This raises
the question of how to improve our model without introducing extra dimensions to
make the POM model more applicable to other initial conditions.

Traditionally, for systems like this one, POMS have to be re-calculated to recon-
struct the model under other conditions. To increase the applicability to a range of
initial conditions, the proposed solution is to add a constraint to the modal coordi-
nates. Its inspiration comes from solving for individual N N Ms [2]. Basically, the idea

is to introduce a constraint coordinate yl obtained from the POD. However, unlike

95

an extra d.o.f., it is a coordinate dependent on modal coordinate 3);; and velocity ()2,

and it can be expressed by

M M 'k
91: f(y2 92) ~23 Z (Unfit/2- (5.15)
j=0k=0

The constraint coordinate yl increases the accuracy of the POM model in addition
to enlarging its applicable domain in the phase space.

In the Simulation, a cubic polynomial was used to approximate the constraint
surface, and the coefficients a“, were estimated by the least mean-square method.
The model differential equation was also modified by introducing y1:

3 3 Wk
9 =2— 2 .. a- 939
1 3-0 k—O M 2 2 , (5.16)

372 = klyl + km + n(—0.990y1 + 0.141y2)3

where k1 and k2 are unknown coefficients of the linear terms, and n is the coefficient

of the known 11:13 nonlinearity. substituting (5.15) into equation (5.16), we can obtain
3 3 ' k 3 ' k 3
y2= l( (92 (5)1520 Z 0),),9239'2 + (7292 + n( 2 Z aj,ky2jy'2 + 0-14192) ,
M: k=0 j=0k=0
(5.17)

which is a simple single d. o. f. nonlinear oscillator. The new model can be regarded

as the added—constraint model.

Table 5.1. Identified parameters of the POM-constraint model (5.16) under the first
three conditions

 

— k1 k2 n
Identified -0.9932 -1.7069 -0. 1607
'Ii‘ue values -0.9532 -1.7212 —0. 1412
Errors 4.2% 0.9% 13.8%

 

 

 

 

 

 

 

 

 

96

 
  

_,,"- -
o‘--
o'~_~--

    

    
  

, 'l'pII'-‘l‘4' 4. o
a" - I’ll; ’0 .‘O ’ s:
. I ’:,l”lllll’.l,”',l."’o1:.$.‘
"I O I'I'IIIIII’II'I'.’
u” l',,;.llllllllllo
II ,1 -...qulluuunlo,
0 _5 fl 11,!" III! I '

  
 

  
  
     
   
  

 

 

     
  
    

   

“
s‘ O
\\\‘ a ‘0‘ .‘
\\ o \
\att‘xxttsmt‘n‘vtv' «:o.
\\\‘ - “\\ ‘s .s
\ \\\\\“\\\'. \\ \“““. \.

    

y-dy2

Figure 5.3. Comparison of the invariant manifold: (a) under small amplitude (the
first three conditions), (b) under all five conditions; Surf 1 is the empirical surface

from POD, Surf 2 is Shaw’s surface, the continuous lines are the orbits under five
conditions; 2: axis is y'g, y axis is yg, z axis is yl.

97

First, the coefficients were estimated under the first three conditions. Displayed
in Table 5.1 are identified parameters of the y2 differential equation, and it shows
a close match with the true values calculated from the linear transformation of the
two-mass equation (5.9) except the nonlinear term. The discrepancy in the coefficient
of the nonlinear term may come from the small amplitude vibration where nonlinear
effect is not so apparent, and the error from digital differentiation. Upon obtaining
the identified equation, Figure 5.3(a) displays the constraint surface (surf 1) fit for
the first three initial conditions which are of smaller oscillation amplitude. The con-
tinuous lines are the corresponding three orbits (the orbits are partially covered by
the surfaces and are hard to see). Surf 2 is the surface obtained by Shaw’s polynomial
approximation of NNM, which is a good approximation of the NNM surface under
small amplitude vibration. For small amplitude vibration, the two surfaces match
well. However, as the vibration amplitude increases, Shaw’s surface deviates from the
true orbits. Instead, the empirical POM model, which was fit from data that was not
‘Small’ in asymptotic terms, shows a better match for more orbits in the phase space.

Figure 5.3(b) Shows the surface fit (Surf 1) for the identified system under all of
the five initial conditions. In Figure 5.1, the dynamical orbit under the fifth initial
condition is also Shown in the phase space (the dotted line), and shows a closer
match with the original system compared to the POD-only model. N oticeably, the
fifth orbit is much larger than the other four orbits. Shaw’s surface shows large
discrepancy from the real orbits for large amplitude vibration. On the other hand,
noticing that the identified model changed after using a larger data set, the empirical
surface changes to fit the fifth orbit, but also shows discrepancy for small amplitude
orbits. The discrepancy results from an inadequate polynomial approximation as well
as the nature of least mean square method. We can adjust our overall approximation
by using a different data set consisting of more records of small amplitude orbits,

or a weighted data set which gives more weight to the small orbits, such that the

98

discrepancy over the whole surface can be more evenly distributed.

Through the analysis above, it is clear that if yl is neglected, i. e. neglecting
equation (5.15) when the energy in yl is negligible comparing to the dominant POM,
the ROM becomes a POM-only model (5.14). The constraint model, however, does
not neglect yl, but incorporates the curvature of the true N NM, to achieve the form
of equation (5.15). Effectively, this puts the model into the form of equation (5.17),
which is a differential equation with nonlinearity in displacement and velocity. By
its construction, it is made to accommodate dynamics of a wide range of initial

conditions.

5.3.3 The example of a mass-pendulum system
The differential equation and the nonlinearity

In the previous example, the added-constraint method was successful in enhancing the
adaptability of the POM reduced order model for a nonlinear system. Nevertheless,
under assumptions of known nonlinearity, the examined system is quite simple. In
practice, the nonlinearity is quite often not clear to us. It is necessary to use function
series to approximate the nonlinearity.

We tested this problem by examining a more complex system: the mass-pendulum
system. Figure 5.4 shows a Sketch of this system, where the mass block has a horizon-
tal displacement 1:, and the arm has an angular displacement 6. The corresponding
differential equation is

:f: + 01117 + b12{flcos6 —- 92 sin 6} = 0

.. ,, , (5.18)
0 + b21 Sind + bggci'cosd = 0

where bu = k/m1+m2, b12 = Time/{(7% + m2)l}, b21 = "l’zgel/JOZ) 922 =

171.26 / J02, ml and 771.2 are the masses of the block and the pendulum (the mass of the

99

 

 

 

Figure 5.4. A sketch of the mass-pendulum system.

bar included), J02 is the mass moment of inertia of the pendulum with respect to the
pendulum joint, e is the distance of the pendulum centroid to the joint, and :i: = z/l
is the non-dimensional displacement. For convenience, we drop the hat of the non-
dimensional displacement (it. The nonlinear terms are strong, including interaction
terms of the mass and the pendulum, mixed nonlinear terms involving displacement,
velocity and acceleration, and trigonometric functions. As a result, it is hard to de-
termine what kind of nonlinear format will show in the POM reconstructed reduced
order model. In the simulation, the initial conditions were chosen to excite only the
first N N M of the system. The POM extraction and the model identification process
were similar to the previous example. However, the two NNMS of the mass-pendulum
system were found only in the absence of whirling.

Now, we denote 21 as the major POM which consists of most of the energy of
the data, and 22 as the minor POM which is orthogonal to 21. The added-constraint

model will be composed of 21 and 22. Two types of the added constraint models were

100

tested. The first one is
Z2 = fc(21, Z1) , (5.19)
Z1 + kziZ1+ P(Z1) = 0
where the differential equation of zl consists of only 21. On the contrary, the second
model, considering 22 as a dependent variable of 21, can be expressed as

Z2 = fc(Z1, Z1) (5 20)

Z1+ kz1Z1 + kz2Z2 + p(Z1, Z2) = 0

In both of the models, the constraint equations are functions of 2'1 and 21. But, in
the first model, the POM differential equation has only one independent variable 21,
and the nonlinearity is approximated by a polynomial of 21. In contrast, the POM
differential equation of the second model has both 21 and 22, and the nonlinearity
is approximated by a polynomial of 21, 22. The second model appears to be a bet-
ter approximation model to the system dynamics, because, much like the previous
example, the 22, as a function of 21 and 2'1, is also introduced in the differential
equation and is generally closer to the actual nonlinearity. It also brings more com-
plexity in the nonlinearity identification, and the problem may be worse for high (1.
o. f. systems. Besides, extra terms in the nonlinearity approximation may unstable
dynamical behavior. For this simulation study, the unknown parameters of equation
(5.19) and (5.20) were identified by the harmonic balance method discussed in the

previous three chapters of the dissertation.

Table 5.2. Physical parameters used in the simulation

 

ml (kg) 0.1 mg (kg) 0.1
k(N/m) 100 I(m) .04

 

 

 

 

 

 

 

101

Simulation setting and the POD

The parameter values used in the Simulation are listed in Table 5.2, and the mass
of the pendulum bar was neglected. Five initial conditions were tested, in which
the initial velocities were zero, and the initial displacements for :1: and 6, from small
to large amplitude, were: (0125,0385), (0250,0733), (0.375,1.10), (0.500,1.54),
(0575,1892). The resulting dynamical behavior corresponded the first N NM of this

system, under five different initial conditions on that NNM.

 

 

theta
O

 

 

 

 

 

 

 

Figure 5.5. The shapes of the lst NNM of the mass-pendulum system under 5 i. c. :
1. (0125,0385), 2. (0250,0733), 3. (0.375,1.10), 4. (0.500,1.54), 5. (0575,1892)

The dominant POM obtained under the first initial condition was 2] =
[ 0,2838 0.9589 ]T, which consists of 99.93% of the total energy. The minor POM

was 23 = [ —0.958902838 ]T. Figure 5.5 shows the five displacement shapes of the

102

first NN M under these non-whirling conditions. It is clear that for large amplitude
vibration (curves 4, 5) the normal mode is quite twisted. If a POM based reduced
order model is built on 21, it is accurate for the first initial condition. However, the
reduced order model under 2] is not applicable to the other initial conditions listed
here. Therefore, one would like to test if the added-constraint model can adapt to

those initial conditions.

Table 5.3. Estimated polynomial coefficients of the constraint equation

 

a1 0.00132 a2 5.96x10“3 (1;; -00105
a4 -000327 a5 1.21x10—4 0.6 2.49

 

 

 

 

 

 

 

 

 

The added constraint and the weighted least mean square solution

In the study, the added constraint was estimated by a cubic polynomial:
_ - 3 - 3 2 - - 2
2:2 — alzl + 0221 + a3zl + a421 + 052121 + 0.62121 , (5.21)

where a, are the unknown coefficients. For better approximation of the small am—
plitude oscillation, the coefficients were solved by the weighted least mean square
method based on the data of all five initial conditions.

We can denote the following equation as the unweighted least mean square matrix
equation:

.45 : q", (5.22)

where (if is the vector of the unknown coefficients, (7 is some known terms without

103

unknown coefficients,

and the small matrices A, are data matrices from free responses to different initial
conditions. Vector f can be solved by this unweighted least mean square method.

Then, we can transform the matrix A to a weighted matrix

w1A1
Aw = 1

w545

where w, is the weight. Correspondingly, we also denote

w1‘71

1% <75

Replacing A, (7 by Aw, (7w, equation (5.22) becomes

which is the weighted least means square method. The advantage of the weighted
solution method is that we can adjust the approximated solution :5 by giving different
weights to subsets. We can give large weights to some subsets of the data if we want
the approximated solution to have small discrepancy with the corresponding subsets.
In the case of the N N M approximation, larger weights were given to the data sets
from smaller oscillations such that the approximated constraint surface can have an

evenly distributed error over all vibration amplitudes. Thus, the result is better

104

 

 

   

 

 

 

 

 

 

10 . 0.015
0.01 r
5 ,
0.005 .
E
1; 0 ‘N 0
‘o
-0.005 r
_5 b
-0.01 >
-10 4 -0.015 A
-0.5 0 0.5 -0.5 0 0.5
21 21

(a) (b)

Figure 5.6. Comparison of Model I (the dotted lines) to the original mass-pendulum
system (the continuous lines) simulated under lst initial condition, (a) phase plane,
(b) mode shape.

than the estimation from the unweighted least mean square method which tends to
accommodate more with the large amplitude orbits.
To the investigated system, the estimated values of the coefficients are listed in

Table 5.3. The constraint were Shared by both of the reduced order models.

Model of type I

Substituting a cubic polynomial into model-1 (5.19), we obtain

3 3 ‘° k
z = E:-_ E: _ a- z 32

21+ kz121+ “11.1213: 0
where 162] = 0.0571, an, 2 —-000468 in the investigation. The nonlinear term is
represented by only a cubic polynomial, and is simple for estimation.
Displayed in figures 5.6, 5.7, 5.8, are the dynamical behaviors of the original
system and the approximated model of the first type. The approximated model was

built on the data from five initial conditions. It shows that the approximated model

105

 

 

   

 

 

 

 

 

 

20, 1 0.15
0.1
10 » i 0.05 »
§ 0’
.— 0
g -0.05+
-10. -0.1
—0.15» <
— L 4
20 A A A A A -0.2 A A A
-1.5 -1 -0.5 0 0.5 1 1.5 -2 -1 0 1 2
21 21

Figure 5.7. Comparison of Model I (the dotted lines) to the original mass-pendulum
system (the continuous lines) simulated under 3rd initial condition, (a) phase plane,
(b) mode shape.

 

 

 

 

 

 

 

 

 

40 r 0 4 -
20 r
'2‘
1;, 0
t:
-20 P
-40 * * -0.4
-2 -1 0 1 2 -2 -1 0 1 2
21 21

Figure 5.8. Comparison of Model I (the dotted lines) to the original mass-pendulum
system (the continuous lines) simulated under 5th initial condition, (a) phase plane,
(b) mode shape.

106

 

 

0.015

0.01 .

1 0.005 '

dz1/dt

-0.005 ’

 

-0.01 .

 

 

 

 

 

 

-0.015 ‘
0.5 -0.5 0 0.5

 

(a) (b)

Figure 5.9. Comparison of Model II (the dotted lines) to the original mass-pendulum
system (the continuous lines) simulated under lst initial condition, (a) phase plane,
(b) mode shape.

is a good representation of the original system. The approximation is better for the
small amplitude vibration conditions, e. g. the first initial condition. As the vibration
amplitude increases, the phase plane curves Show increasing difference between the
full and reduced order models. Part of the reason could be that the cubic nonlinearity
of 21 is inadequate to represent the nonlinearity of the original system. Meanwhile,
by mode shape comparison in figures 5.6, 5.7, 5.8, the vibration mode shapes of the
identified model are close to the real N NM Shapes, and the minor displacement 22
is properly approximated. Noticeably, the energy in 22 coordinate increases as the
oscillation amplitude increases. This model has the same differential equation of 21
as the POM-only model, implying same dynamic behavior of 21. The only difference
is that the POM-only model consider 22 = 0, whereas, in the present mode, 22 is a
curve which can accommodate the change of initial conditions. Obviously, comparing
to the POM-only model, the added-constraint made the response of the ROM closer
to the actual first N N M of the original system, and thus improved the approximation

of the dynamic behavior of the original system.

107

 

 

 

 

 

 

 

 

 

0.15 - -
2° ’ 0.1
10 , 4 0.05 ’
5 0 .
1- 0 t
8 —0.05 -
-10 ’ ‘ -0.1 >
_ l
_20 , 0.15
—0.2 A
-2 -1 0 1 2 -2 -1 0 1 2
21 21
(a) (b)

Figure 5.10. Comparison of Model II (the dotted lines) to the original mass-pendulum
system (the continuous lines) simulated under 3rd initial condition, (a) phase plane,
(b) mode Shape.

 

 

 

 

 

 

 

 

 

40 - - a . - 0.4
20 L 1
fi
1'. 0 ‘
'O
—20 . 1
—40 A A -0 4
-3 -2 -1 0 1 2 3 -2 -1 0 1 2

Figure 5.11. Comparison of Model II (the dotted lines) to the original mass-pendulum
system (the continuous lines) simulated under 5th initial condition, (a) phase plane,
(b) mode shape.

108

Model of type II

A more complex polynomial was used to estimate the nonlinearity in model-2 (5.20),
and it was
Z2 = 23320 233:0 aj.kZ1jZ'1‘° , (525)
2’1 + kzlzl + 152222 + dlzi’ + dngzg + (132le = 0
where the identified values are kzl = 0.0403, 1622 = 0.383, (11 = 0.00218, d2 = —0.146,
and d3 = —0.453 for the mass-pendulum system. The nonlinear terms are up to the
cubic order of 21 and 22. Theoretically, the nonlinearity approximation is better than
that of the model-I.

N otieeably, the 2% term was not included in the 21 differential equation of (5.25). It
was found that there was unstable behavior of the identified system if 2% nonlinearity
was included, and the corresponding response tended to be quasi-periodic. It looks
like a trial and error process, but it is not. Actually, the original systems also showed

3 term

the bifurcation to quasi-periodic response as the initial offset increases. The 22
made the response of the reduced-order model quasi-periodic at a smaller initial offset.
The identification process is performed on state values during responses, and not on
stability characteristics. Therefore, in sensitive systems, while a model is identified
for response values, it may be close to a bifurcation set, and may not have the right
stability. Thus, the model of equation (5.25) was selected for better approximation
of the original system. It shows that the nonlinear terms in the identification model
can bring unexpected results in the simulation. To the added constraint model, since
another dependent variable 22 is included, the problem becomes more complicated.
For this example, the ‘full’ set of nonlinear terms actually deteriorated the system
dynamics of the second model. However, the dilemma is that one does not know in
advance what terms are unnecessary, or undesired. The problem has to be solved by

reducing the higher order terms of the approximation functions, which increases the

difficulty of system identification.

109

1.5:

 

 

, f -. I: -
1 .’ \ /.-' i:
‘ l 1 .- 1:
i I '.
‘_ i I \ \2
N 0 " , l \. l: 1..
l 1 °_ . .
_ l \, I. \.
-0.5- ' ‘ I; '- " i:
l - \1 .
, ’: . ’ \1
-1 \‘ )5 ).\ l \
_15 1 1 1 . .
0 20 40 60 80 100
time (s)

Figure 5.12. A time domain comparison between the original system (continuous
line), the Model—I (dashed line), the Model—II (dotted line) under the 3rd initial
condition.

Figures 5.9, 5.10, 5.11 show the comparison between model type II (5.25) and
original system by phase planes and vibration mode Shapes. Due to the effect of more
nonlinear terms, by comparing the phase plane orbits, this model better approximates
to the original system in large amplitude vibration (Fig. 5.10, 5.11) than the first
model. On the other hand, the approximation of the small amplitude vibration is
not so satisfying comparing to the model of the first type. Possible reason can be the
influence of the nonlinear terms of 22.

Displayed in Figure 5.12 is a time domain signal comparison between the original
system and the two models. The time domain responses of the three systems are
very similar. But, the vibration frequencies are lightly different. In Figure 5.13, the
variation of the natural frequencies is compared among the three systems. It shows

that the reduced order model can also display inaccuracy in frequency domain. The

110

0.038

 

0.037 r

_o
o
00
a)

0.035 '

0.034 -

natural freq. (hz)

0.033 '

 

 

 

 

0.032 * ‘
0 0.02 0.04 0.06

potential energy (J)

Figure 5.13. Natural frequency comparison under different initial conditions (repre-
sented by the maximum potential energy) between the original system (continuous
line), the Model—I (dashed line), the Model—II (dotted line); the ‘*’ and ‘0’ mark the
five initial conditions for Model-I and Model-II.

111

natural frequency change of the original system is a hardening curve as the vibration
amplitude (maximum potential energy) increased. Whereas, the curve of model-I is
a softening curve which resulted from the negative identified value of an, in (5.24),
and the corresponding frequency discrepancy is large. The model of the second type
is somewhat a better approximation of the original system in the natural frequency

comparison.

5.4 Concluding remarks

In this chapter, the POM reduced order models of nonlinear systems was investigated.
Unlike the model order reduction techniques in turbulence [3, 32, 33], the complete
process of the proposed idea is to reconstruct and identify the governing differential
equations based upon proper orthogonal modes. The advantages of the method are
the abilities to simplify and identify an unknown complex system simultaneously. But,
it also brings the difficulty of identifying correctly the nonlinearity of the original
systems. Moreover, for systems with strong nonlinearity, the POM based reduced
order model can not give accurate prediction of the dynamics in conditions other than
the given initial conditions upon which the proper orthogonal modes were extracted.
Hence, the dissertation study generalized the POM based method by introducing
an added-constraint such that the reduced order model can be more accurate in
prediction of system dynamics as well as applicable to different initial conditions.
Two examples were examined: a two-mass system with a hardening spring and
a mass-pendulum system. Both of the two examples show success of the added-
constraint model, in which, due to the strong nonlinearity, the basic POM models
are not adequate in Simulating system behaviors under different initial conditions.
Meanwhile, in the Simulation study of the mass-pendulum system, two types of added-

constraint models were investigated: one with only the POM coordinates and the

112

other one with POM coordinates and added-constraint coordinates involved. The
former one was found to be more suitable for small amplitude vibration. In this
example, however, the latter one was shown to be more accurate in predicting large
amplitude vibration.

The study was successful in the two examples of low-dimensional systems, which
are under free vibration with no damping. However, more experimental and sim-
ulation studies are needed for forced vibration systems, high-dimensional complex
systems, and chaotic systems. The questions involve what POMS are needed for
reduced order modeling, how to identify nonlinear terms in multi-degree-of-freedom
systems, how to apply the added-constraint method to improve the accuracy and
applicability of the POM model, and determining which type of models are best for
the investigated systems. If all of the obstacles can be solved, the added-constraint

POM model can be very powerful in system identification and simulation.

113

CHAPTER 6

Conclusions and Future Work

6.1 Main contributions

The purpose of this dissertation was to provide a handy and useful tool for parametric
identification of nonlinear/chaotic systems, and also to reconstruct reduced order
models of complex systems based upon proper orthogonal modes and parametric
identification methods.

The parametric identification method is based upon harmonic balance method
[24], and focused on chaotic systems. Due to the nature of chaos—no stable steady
state behavior—the identification process first refers to the extraction of the unstable
periodic orbits, which form a dense set existing within the chaotic set. Then the
harmonic balance method can be applied upon the UPOS to identify the unknown
parameters. The dissertation study did a thorough investigation on the application
of this method: simulation and experimental study, multi-degree-of-freedom systems,

and whirling systems. The major achievements include

1. The method was examined successfully on whirling systems, e. g. single-
pendulum and double-pendulum systems. In both of the cases, only angular
displacement signals were supposed to be known, consistent with the limita-

tions of practical experiments. The special feature of whirling systems was

114

that the angular displacements were not continuous signals, but signals in 31
space [—1r,7r). Thus, special modifications were made in the process of UPO

extraction and the approximation of velocity and acceleration signals.

. Noise reduction techniques were applied in the experimental investigation of
multi-degree-of-freedom systems, e. g. the double pendulum system with base
excitation. There were three reasons for applying noise reduction: the noise
contamination in the signals, the recurrence error of the extracted UPO, and
the strong nonlinearity which magnified the noise influence and propagated
it from high to low frequencies. The noise reduction techniques include low-
paSS-filtering, digital differentiation, and choice of sub-harmonics by which the
less noise-contaminated sub-harmonics of the UPOS were selected for harmonic

balance identification.

. Based upon the understanding of the least mean square method, an error quan-
tification method was proposed to examine the accuracy of the estimated param-
eters. The quantification ideas were based upon linear regression and statistical
methods. Also, by the error quantification process, an optimization process was
established to improve the identification results by an Optimized selection of
sub-harmonics. The indicators of a good identification were small identification

error value, and an even and random distribution of the estimation errors.

. With more understanding of the chaotic set and the unstable periodic orbits, un-
der the assumption of hyperbolic chaotic attractor and unstable periodic orbits
contained within, the identification process which originally involved the extrac-
tion of UPOS and the harmonic balance method was Simplified to a nearly pure
frequency method - simple harmonic balance identification. The simplification
was due to the fact that any chaotic orbit, if long enough, is an approximation

of some UPO. The longer the chaotic orbit, the better the approximation. This

115

property was justified with simple symbol dynamics. The calculation examples
showed that the simplification could identify the parameters with satisfying re-
sults. The meaning of this simplification should not be under-estimated. The
most direct result is to make the identification process fast and requires rela-
tively small data set, i. e. it is more efficient. Hence, it is possible to apply this
method to cases of real-time identification. Meanwhile, the simplification also
makes the identification of all nonlinear systems possible by only one method -

the harmonic balance identification.

With the progress in identification techniques, the focus of the dissertation then
switched to the reduced order models by means of proper orthogonal modes, e. g. the
reduced order models for systems with strong nonlinearity. The POMS are experimen-
tal optimized orthogonal modes. They have been useful for model order reduction.
The dissertation study goes further in the direction of model order reduction, i. e.
reconstruction of reduced order models for unknown systems. Two simple examples
were given to test the method of added-constraint models. The special feature of the
added-constraint method iS that extra constraint equations are added to minor mode
such that the constraint and the differential equation of the dominant POM together

form the new reduced order models. Through two example systems, it revealed:

1. In both the studies, the nonlinear added-constraint method showed better ap-
proximation of the systems’ dynamics than the POM-only reduced order models.
Also, the method made the reduced order models applicable to more general ini-
tial conditions. Compared to the nonlinear normal modes, the added-constraint
method is empirical and able to adjust the approximation model by different

data, i. e. different initial conditions.

2. The identification of the nonlinear terms was more difficult if little knowledge

was known for the nonlinearity of the original systems. Two types of nonlin-

116

ear models were examined in the mass-pendulum model: the model with no
constraint modal coordinates in the governing differential equations and the
one with constraint modal coordinates in the differential equations. The first
model was simpler, and more accurate for small amplitude oscillation. On the
contrary, the second model was more feasible for large amplitude oscillation,
and showed the difficulty in finding appropriate nonlinear terms in the model

to represent the nonlinearity of the original systems.

6.2 Directions for future work

There were achievements through this study on both parametric identification method
and the added-constraint POM reduced order models. However, many questions still

remain:

o In the experimental study, the identification of friction coefficients were erro-
neous. The reason was due to the small values of those parameters. It indicates
that the small parameters are more vulnerable to noise contamination, and thus,
hard to identify. Therefore, the energy balance method was applied to identify
the friction parameters. For non-damping small terms, e. g. small geometric
terms, similar procedure can also be developed. Suppose we have a single d. o.

f. vibration system expressed as
5 + k2: + c5 + 6,,53 -_— f(t), (6.1)

where f (t) is a periodic excitation, and an is a small nonlinear coefficient com-
pared to linear coefficient k and c. Like the identification to small damping
terms in previous examples, the estimation of an will be hard in harmonic iden-

tification process, due to the influence of other large terms. To reduce the noise

117

from other terms, multiplying (6.1) and integrating over one cycle of a periodic

orbit, we obtain

t+T __ t+T 2 t+T . t+T 4 t+T
ft rsrdt+ ft k5 dt+ ft crxdt-i—an [t a: dt= ft f(t):rdt.
(6.2)
N oticeably, the integration of the velocity term over one cycle of a periodic orbit

is zero. Hence, equation (6.2) can be transformed to
t—l—T 4 t+T ,2 t+T 2 t+T
an/t :1: dt :/t :1: dt _/t kit dt+/t f(t):1:dt. (6.3)

The an can be solved by least mean square method if there are more than one
periodic orbit. The advantage of equation (6.3) is that there are fewer, and
only displacement and velocity terms. If the displacement Signal is directly
measured, it is believed to be quite accurate. On the other hand, the velocity
and acceleration Signals are to be obtained by the Fourier series method or
digital differentiation, which have been proved to have more high-frequency
noise. The high frequency noise is the worst for the acceleration signal. Hence,
we are quite confident that the estimation of an by (6.3) is more accurate
than the harmonic balance identification. For multi-degree-of—freedom systems,
similar methods are expected to be developed and improve the accuracy of

nonlinearity identification.

The simplification of the harmonic balance identification method makes it pos-
sible for real-time identification of chaotic systems. As we know that every
chaotic orbit can be an approximated orbit of some UPO, real-time identifica-
tion can be possible if the length of the UPO can be properly chosen. The goal
is to identify relatively slowly varying parameters of the governing equation,

or to assess damage of systems. Meanwhile, due to the similarity between the

118

harmonic balance method and the wavelet analysis, the real-time identification
of chaotic systems can also be done by wavelet analysis. The questions are how

fast can the process be and the issue of noiseresistance.

It is also possible that the identification method can be applied to nonlin-
ear/chaotic systems under random or impulse excitation. The idea lies in the
analogue between the random excitation and the impulse excitation. The basic
root of this idea still comes from re-thinking of the unstable periodic orbits and
chaotic set. But it Should be considered more in the frequency domain, instead

of the time domain.

The POD process for whirling systems requires study. The whirling Signals, i. e.
the angular displacements, are not applicable to the traditional definition of the
proper orthogonal decomposition, since these signals are in a spherical Space.
Hence, the reduced order models are hard to apply to whirling systems. The
possible solution can be to use the complex variables, instead of real variables,

to represent the whirling signals.

Nonlinear reduced order modeling is difficult. The approximation method can
be troublesome since unnecessary nonlinear terms can bring undesired dynam-
ical behavior for strongly nonlinear systems. The problem will become more
difficult for complex structures where the nonlinearity is not clear. An idea

might be to incorporate stability measures into identification process [67].

The added-constraint model, though successfully explained in the two examples,
is still in development, and needs to be examined by complex, high-dimensional
structures. The problems involve the determination of POM coordinates and
the constraint constraints, and how far can the method improve the performance

of the reduced order models.

119

BIBLIOGRAPHY

[1] K. T. Alligood and et a1. CHAOS An Introduction to Dynamical Systems. New
York: Springer—Verlag, 1997.

[2] S. W. Shaw and C. Pierre, “Normal modes for non-linear vibratory systems,”
Journal of Sound and Vibration, vol. 164, no. 1, pp. 85—124, 1993.

[3] J. L. Lumley, Stochastic Tools in Turbulence. New York: Academic Press, 1970.

[4] S. Han and B. F. Feeny, “Enhanced proper orthogonal decomposition for the
modal analysis of homogeneous structures,” Journal of Vibration and Control,
vol. 8, pp. 19—40, 2002.

[5] K. S. Mohammad and et al., “Direct parameter estimation for linear and non-
linear structures,” Journal of Sound and Vibration, vol. 153, no. 3, pp. 471—499,
1992.

[6] M. R. Hajj, J. Fung, and et al., “Parametric identification of nonlinear dynamic
systems,” Computers 8 Structures, vol. 20, no. 3, pp. 487—493, 1985.

[7] Q. Chen and G. R. Tomlinson, “Parametric identification of systems with dry
friction and nonlinear stiffness a time series model,” Journal of Vibration and
Acoustics, vol. 118, no. 2, pp. 252—263, 1996.

[8] R. K. Kapania and S. Park, “Parametric identification of nonlinear structural
dynamic systems using time finite element method,” AIAA Journal, vol. 35,
no. 4, pp. 719—726, 1997.

[9] G. Verbeek, A. Dekraker, and D. H. Vancampen, “Nonlinear parametric identifi-
cation using periodic equilibrium states - application to an aircraft landing gear
damper,” Nonlinear Dynamics, vol. 7, no. 4, pp. 499—515, 1995.

[10] M. R. Hajj, J. Fung, and et al., “Damping identification using perturbation tech-
niques and higher order spectra,” Nonlinear Dynamics, vol. 23, no. 2, pp. 189—
203, 1986.

[11] G. Kerschen, V. Lenaerts, and et al., “Identification of a continuous structure
with a geometrical non-linearity. part i: Conditioned reverse path method,” Jour-
nal of Sound and Vibration, vol. 262, no. 4, pp. 889—906, 2003.

120

[12] C. M. Richard and R. Singh, “Identification of multi-degree—of-freedom non-linear
systems under random excitations by the reverse-path spectral method,” Journal
of Sound and Vibration, vol. 213, no. 4, pp. 673—708, 1998.

[13] O. Gottlieb and M. Feldman, “Application of a hilbert-transform based algorithm
for parameter estimation of a nonlinear ocean system roll model,” Journal of
Offshore Mechanics and Arctic Engineering, vol. 119, pp. 239—243, 1997.

[14] M. Feldman, “Non-linear free vibration identification via the hilbert transform,”
Journal of Sound and Vibration, vol. 208, no. 3, pp. 475—489, 1997.

[15] M. Feldman, “Non-linear system vibration analysis using hilbert transform—i. free
vibration analysis method ‘freevib’,” Mechanical Systems and Signal Processing,
vol. 8, no. 2, pp. 119—127, 1994.

[16] N. P. Plakhtienko, “Methods of identification of nonlinear mechanical vibrating
systems,” International Applied Mechanics, vol. 36, no. 12, pp. 1565—1594, 2000.

[17] N. P. Plakhtienko, “A method of special weight-functions in the problem of
parametric identification of mechanical systems,” DOPOV AKAD NA UK A,
vol. 8, pp. 31—35, 1983.

[18] K. Yasuda, S. Kawamura, and et al., “Identification of nonlinear multi-degree-
of-freedom systems (presentation of a technique),” JSME International Journal,
Series 111, vol. 31, pp. 8-14, 1988.

[19] K. Yasuda and S. Kawamura, “A nonparametric identification technique for non-
linear vibratory systems (proposition of a technique),” JSME International Jour-
nal, Series 111, vol. 31, pp. 365—372, 1989.

[20] K. Yasuda and K. Kamiya, “Identification of a nonlinear beam (proposition of
an identification technique),” JSME International Journal, Series III, vol. 33,
pp. 535—540, 1990.

[21] K. Yasuda and K. Kamiya, “Experimental identification technique of nonlinear
beams in time domain,” Nonlinear Dynamics, vol. 18, pp. 185—202, 1997.

[22] M. Thothadri and et al., “Nonlinear system identification of multi-degree—of—
freedom systems,” Nonlinear Dynamics, vol. 32, pp. 307—322, 2003.

[23] R. Ghanem and F. Romeo, “A wavelet-based approach for model and param-
eter identification of nonlinear systems,” International Journal of Non-linear
Mechanics, vol. 36, pp. 835—959, 2001.

[24] C. M. Yuan and B. F. Feeny, “Parametric identification of chaotic systems,”
Journal of Vibration and Control, vol. 4, no. 4, pp. 405—426, 1998.

[25] B. F. Feeny, C. M. Yuan, and et al., “Parametric identification of an experimental
magneto-elastic oscillator,” Journal of Sound and Vibration, vol. 247, no. 5,
pp. 785—806, 2001.

121

[26] B. F. Feeny and R. Kappagantu, “On the physical interpretation of proper or-
thogonal modes in vibrations,” Journal of Sound and Vibration, vol. 214, no. 4,
pp. 607—616, 1998.

[27] B. F. Feeny, “On the proper orthogonal modes and normal modes of continuous
vibration systems,” Journal of Vibration and Acoustics, vol. 124, no. 1, pp. 157-
160, 2002.

[28] B. F. Feeny, “On proper orthogonal coordinates as indicators of modal activity,”
Journal of Sound and Vibration, vol. 255, no. 5, pp. 805—817, 2002.

[29] B. F. Feeny and Y. Liang, “Interpreting proper orthogonal modes of randomly
excited vibration systems,” Journal of Sound and Vibration, vol. 265, no. 5,
pp. 953-966, 2003.

[30] G. Kerschen and J. C. Golinval, “Physical interpretation of the proper orthogonal
modes using the Singular value decomposition,” Journal of Sound and Vibration,
vol. 249, no. 5, pp. 849—865, 2002.

[31] S. Han and B. Feeny, “Application of proper orthogonal decomposition to struc-
tural vibration analysis,” Mechanical Systems and Signal Processing, vol. 17,
no. 5, pp. 989—1001, 2002.

[32] J. L. Lumley, “The structure of inhomogeneous turbulent flow,” Atmospheric
Turbulence and Radio Wave Program, pp. 166—178, 1967. (A. M. Yaglom and
V. I. Tatarski, editors) Moscow: Nauka.

[33] G. Berkooz and et al., “The proper orthogonal decomposition in the analysis of
turbulent flows,” Annual Review of Fluid Mechanics, vol. 25, pp. 539—575, 1993.

[34] W. Z. Lin and et al., “Proper orthogonal decomposition and component mode
synthesis in macromodel generation for the dynamic Simulation of a com-
plex mems device,” Journal of Micromechanics and Microengineering, vol. 13,
pp. 646—654, 2003.

[35] M. S. Riaz and B. F. Feeny, “Proper orthogonal modes of a beam sensed with
strain gages,” Journal of Vibration and Acoustics, vol. 125, pp. 129—131, 2003.

[36] A. Sarkar and M. P. Paidoussis, “A cantilever conveying fluid: coherent modes
versus beam modes,” International Journal of Non-linear Mechanics, vol. 39,
pp. 467—481, 2004.

[37] J. P. Cusumano and B. Y. Bai, “Period-infinity periodic motions, chaos, and
spatial coherence in a 10-degree-of-freedom impact oscillator,” Chaos, Solutions,
and Fractals, vol. 3, pp. 515—535, 1993.

[38] J. P. Cusumano and et al., “Spatial coherence measurements of a chaotic flexible-
beam impact oscillator,” Aerospace Structures: Nonlinear Dynamics and System
Response, ASME-AD, vol. 33, pp. 13—22, 1993.

122

[39] V. Lenaerts and et al., “Proper orthogonal decomposition for model updating
of non-linear mechanical systems,” Mechanical Systems and Signal Processing,
vol. 15, no. 1, pp. 31-43, 2000.

[40] I. T. Georgiou and I. B. Schwartz, “Dynamics of large scale coupled struc-
tural / mechanical systems: a singular perturbation/ proper orthogonal decompo-
sition approach,” Journal of Applied Mathematics, vol. 59, no. 4, pp. 1178—1207,
1999.

[41] R. Kappagantu and B. F. Feeny, “An ‘optimal’ modal reduction of a system with
frictional excitation,” Journal of Sound and Vibration, vol. 224, no. 5, pp. 863—
877, 1999.

[42] R. V. Kappagantu and B. F. Feeny, “Part1: dynamical characterization of a
frictionally excited beam,” Nonlinear Dynamics, vol. 22, pp. 317—333, 2000.

[43] R. V. Kappagantu and B. F. Feeny, “Part2: proper orthogonal modal modeling
of a frictionally excited beam,” Nonlinear Dynamics, vol. 23, pp. 1—11, 2000.

[44] G. Kerschen, B. F. Feeny, and et al., “On the exploitation of chaos to build
reduced-order models,” Computer Methods in Applied Mechanics and Engineer-
ing, vol. 192, no. 13-14, pp. 1785—1795, 2003.

[45] W. V. dc Water and M. Hoppenbrouwers, “Unstable periodic orbits in the para-
metrically excited pendulum,” Physical Review A, vol. 44, no. 10, pp. 6388—6397,
1991.

[46] J. J eong and S.-Y. Kim, “Bifurcations in a horizontally driven pendulum,” Jour-
nal of Korean Physical Society, vol. 35, no. 5, pp. 393—398, 1999.

[47] S. R. Biship and M. J. Clifford, “Zones of chaotic behavior in the parameterically
excited pendulum,” Journal of Sound and Vibration, vol. 189, pp. 142—147, 1996.

[48] R. V. Dooren, “Chaos in a pendulum with forced horizontal support motion: a
tutorial,” Chaos, Solutions and Fractals, vol. 7, no. 1, pp. 77—90, 1996.

[49] V. B. Raybov, “Using lyapunov exponents to predict the onset of chaos in non-
linear oscillators,” Physical Review E', vol. 66(1), no. 016214, 2002.

[50] W. Szemplinska-Stupnicka and E. Tyrkiel, “The oscillation-rotation attractors
in the forced pendulum and their peculiar properties,” International Journal of
Bifurcation and Chaos, vol. 12, no. 1, pp. 159—168, 2002.

[51] M. Clerc, P. Coullet, and E. Tirapegui, “The stationary instability in quasi-
reversible systems and the lorenz pendulum,” International Journal of Bifurca-
tion and Chaos, vol. 11, no. 3, pp. 591—603, 2002.

[52] H. D. Abarbanel, R. Brown, and et al., “The analysis of observed chaotic data
in physical systems,” Reviews of Modern Physics, vol. 65, pp. 1331-1392, 1993.

123

[53] D. Auerbach and et al., “Exploring chaotic motion through periodic orbits,”
Physical Review Letters, vol. 58, pp. 2387—2389, 1987.

[54] D. P. Lathrop and E. J. Kostelich, “characterization of an experimental strange

attractor by periodic orbits,” Physical Review A, vol. 40, pp. 4028—4031, 1989.

[55] Z. Al-Zamel and B. F. Feeny, “Improved estimations of unstable periodic orbits
extracted from chaotic sets,” 2001 ASME Design Engineering Technical Confer-
ences, Proceedings of Det’01, 2001.

[56] J. W. Liang and B. F. Feeny, “Identifying damping in forced oscillators using
a harmonic energy balance,” IMECE-2002, 11 2002. New Orleans, IMECE2002
CD-ROM.

[57] J. P. Eckmann and D. Ruelle, “Ergodic theory of chaos and strange attractors,”
Reviews of Modern Physics, vol. 57, no. 3, pp. 617—656, 1985.

[58] O. Biham and W. Wenzel, “Characterization of unstable periodic orbits in
chaotic attractors and repellers,” Physical Review Letters, vol. 63, no. 8, pp. 819—
822, 1987.

[59] P. Cvitanovié, “Invariant measurement of strange sets in terms of cycles,” Phys-
ical Review Letters, vol. 61, no. 24, pp. 2729-2732, 1987.

[60] G. H. Gunaratne and I. Procaccia, “Organization of chaos,” Physical Review
Letters, vol. 59, no. 13, pp. 1377—1380, 1987.

[61] K. Mischaikow and et al., “Construction of symbolic dynamics from experimental
time series,” Physical Review Letters, vol. 82, no. 6, pp. 1144—1147, 1999.

[62] C. M. Yuan, Parametric identification of chaotic systems, PHD dissertation.
East Lansing, MI 48824: Michigan State University, 1995.

[63] L. Sirovich, “Turbulence and the dynamics of coherent structures, part i: Coher-
ent structures,” Quarterly of Applied Mathematics, vol. 45, pp. 561—571, 1987.

[64] I. T. Georgiou and I. B. Schwartz, “Invariant manifolds, nonclassical normal
modes, and proper orthogonal modes in the dynamics of the flexible spherical
pendulum,” Nonlinear Dynamics, vol. 25, pp. 3-31, 2001.

[65] E. Pesheck and et al., “Nonlinear modal analysis of structural systems using
multi-mode invariant manifolds,” Nonlinear Dynamics, vol. 25, pp. 183-205,
2001.

[66] D. Ruelle and F. Takens, “On the nature of turbulence,” Communications in
Mathematical Physics, vol. 20, pp. 167-192, 1971.

[67] M. Thothadri and F. C. Moon, “Nonlinear system identification of systems with

periodic limit—cycle response,” Nonlinear Dynamics, vol. 39, no. 1—2, pp. 63-77,
2005.

124