.-

z: .u

.I .H. 1:.7.

.»...

 

 

 

1 .fifiv...‘ V...” .

i _. r. . ., . . . .3: .1,3‘. . . ., .., . ‘ ‘ ‘
,mmowaa‘ , _ ‘ . , 3%w3w3w33... A... a» , A . . , . . . ‘ . . ‘ 1:. Km.

. 3,. 4::

 

THESIS

I
30m

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the

thesis entitled

Phase Space Reconstruction by Alternative Methods

presented by
Guojian Lin

has been accepted towards fulfillment
of the requirements for

MS . Mcdnnicnl Engineering
degree 1n

' I
Major prof74r

 

Date 8/’ (9/0/

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

 

 

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6101 cJCIRC/DatoDuost-pis

 

PHASE SPACE RECONSTRUCTION BY ALTERNATIVE METHODS
By

Guojian Lin

AN ABSTRACT OF A THESIS
Submitted to

Michigan State University

In partial fulfillment of the requirements
for the Degree of

MASTER OF SCIENCE

Department of Mechanical Engineering
2001

Professor Brian Feeny

ABSTRACT
PHASE SPACE RECONSTRUCTION BY ALTERNATIVE METHODS
By

Guojian Lin

Two alternative methods in phase space reconstruction are set forward and
studied: a fractional derivative method and a filter method. The fractional derivative
method is tested on sinusoidal signals first, and then applied to three systems: the Duffing
oscillator, the Lorenz system and a two-well experimental system. The derivative leakage
problem caused by fast Fourier transforms (FFT) is studied and a criterion is developed to
measure the distortion of a one-derivative signal calculated by the PFI‘ method. Based on
this criterion, 3 simple but effective method—truncation method is put forward.
Fractional derivatives are then used as pseudo phase vectors. The embedding dimension,
unstable periodic orbits and correlation dimension are studied for the three systems. The
filter method is tested on sinusoidal signals first, and applied to three systems: the
Duffing oscillator and two Coulomb oscillators. Optimal values for the filter-parameter
are determined. The filtered signal and delayed signal compose the pseudo space vectors.
The embedding dimension is computed by the false nearest neighbor (FNN) method. The
study of this thesis shows the two alternative methods can be effective phase space

reconstruction tools.

ACKNOWLEDGEMENTS

I would like to express my deepest appreciation to the following people and
institutes for having made the completion of this thesis possible.

To the National Science Foundation for funding the research in this thesis (NSF,
grant No. DMI-9800323)

To NASA, grant No. NAG-l-01048. Thanks to Dr. Walt Silva of the NASA
Langley Research Center, the contract monitor,

To my advisor Dr. Brian Feeny, for his guidance and invaluable assistance,
without which this thesis would not have been possible.

To my other committee members, Dr. Ranjan Mukherjee and Dr. Hassan Khalil
for their kind help and useful advice.

To Joe Cusumano of Penn State University, for the experimental data of the two-
well system.

To J. W. Liang of the Ming-Chi Institute of Technology, for the Fortran 77
program to calculate FNN.

To C. M. Yuan of Chung Shan Institute of Technology, for the Fortran 77
program to calculate correlation dimension and MATLAB program to calculate unstable
periodic orbits.

To all my other instructors and colleagues at MSU, for their skillful instructions in
the classroom and lab, but more importantly, for their friendship and constant

encouragement.

iii

TABLE OF CONTENTS

LIST OF FIGURES ................................................................................. vi

CHAPTER 1

INTRODUCTION .................................................................................... 1
1.1 Thesis Statement ........................................................................ 1
1.2 Nonlinear System and Dynamics ...................................................... 1
1.3 Phase Space Reconstruction ........................................................... 3
1.4 FFI‘ (Fast Fourier Transform) ......................................................... 8
1.5 Overview ................................................................................ 12

CHAPTER 2

SYSTEMS AND TOOLS .......................................................................... 14
2.1 Introduction ............................................................................. 14
2.2 System Simulation ..................................................................... 14
2.3 Duffing System ........................................................................ 14
2.4 Coulomb Friction Oscillator System (Oscillator l) ............................... 15
2.5 Another Friction Oscillator System (Oscillator 2) ................................. 17
2.6 Lorenz System ......................................................................... 18
2.7 Two-well System ...................................................................... 20
2.8 Average Mutual Information ......................................................... 21
2.9 False Nearest Neighbors .............................................................. 22
2.10 Unstable Periodic Orbits ............................................................ 23
2.11 Correlation Dimension .............................................................. 23

CHAPTER 3

PHASE SPACE RECONSTRUCTION BY MEANS OF FRACTIONAL

DERIVATIVES ..................................................................................... 24
3.1 Motivation .............................................................................. 24
3.2 Test on Sinusoidal Time Series ...................................................... 27
3.3 Leakage Problem ....................................................................... 28
3.4 Look for a “Good” Value for a ...................................................... 35
3.5 FNN Analysis .......................................................................... 47
3.6 Extraction of Unstable Periodic Orbits .............................................. 52
3.7 Correlation Dimension ................................................................ 56

CHPATER 4

RECONSTRUCTION WITH A FILTER OPERATOR ....................................... 58
4.1 Motivation .............................................................................. 58
4.2 Basic Idea ............................................................................... 58
4.3 Theoretical Analysis ................................................................... 59
4.4 Test on Sinusoidal Time Series ....................................................... 60
4.5 Look for “good” Value for a ......................................................... 62

iv

4.6 FNN Analysis .......................................................................... 72

Chapter 5

CONCLUSIONS .................................................................................... 76
5.1 Conclusions ............................................................................. 76
5.2 Future Work ............................................................................ 77

BIBLIOGRAPHY ................................................................................... 78

1.1

1.2

1.3

2.1

2.2

2.3

2.4

2.5

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8

3.9

3.10

3.11

3.12

3.13

LIST OF FIGURES

Schematic representation of embedding process ........................................ 6
Derivative by FFT (no leakage) .......................................................... 11
Derivative by FFT (leakage) .............................................................. 12
Duffing oscillator system .................................................................. 15
Coulomb friction oscillator 1 ............................................................. 16
Coulomb friction oscillator 2 ............................................................. l7
Lorenz system example (waterwheel) ................................................... l9
Two-well system example ................................................................ 20
Bode diagram of s“ (a=1) ................................................................. 27
Fractional derivative Lx(n) and original signal x(n) ................................... 28
Rectangle window .......................................................................... 29
Near-rectangle window .................................................................... 30
Kaiser-Bessel window ..................................................................... 30
Leakage error of the Duffing system .................................................... 33
Leakage error of the Lorenz system ...................................................... 34
Leakage error of the two-well system ................................................... 34
Ivs. h for the Duffing system........................... ................................. 35
I vs. a of the Duffing system .............................................................. 36
2-D pseudo phase space of the Duffing system ........................................ 37
True 2-D phase space of the Duffing system ........................................... 38
2-D pseudo phase space of the Duffin g system ....................................... 39

vi

3.14

3.15

3.16

3.17

3.18

3.19

3.20

3.21

3.22

3.23

3.24

3.25

3.26

3.27

3.28

3.29

3.30

3.31

3.32

3.33

3.34

3.35

4.1

I vs. h of the Lorenz system ............................................................... 40

I vs. a of the Lorenz system ............................................................... 40
2-D pseudo phase space of the Lorenz system ......................................... 41
True 2-D phase space of the Lorenz system ............................................ 42
2-D pseudo phase space of the Lorenz system ......................................... 43
I vs. h of the Two-well system ............................................................ 44
I vs. a of the Two-well system ............................................................ 44
2-D pseudo phase space of the Two-well system ...................................... 45
True 2-D phase space of the Two-well system ........................................ 46
2-D pseudo phase space of the Two-well system ...................................... 47
FNN vs. d (Duffing) ........................................................................ 48
FNN vs. (I (Duffing, MOD) ............................................................... 49
FNN vs. d (Lorenz) ........................................................................ 49
FNN vs. d (Lorenz, MOD) ................................................................ 50
FNN vs. d (Two-well) ..................................................................... 51
FNN vs. d (Two-well, MOD) ............................................................. 51
Recurrent points (Duffing) ................................................................ 53
Recurrent points (Duffing, true phase space) .......................................... 53
Recurrent points (Lorenz) ................................................................. 54
Recurrent points (Lorenz, true phase space) ............................................ 55
Recurrent points (Two-well) .............................................................. 55
Recurrent points (T wo-well, true phase space) ........................................ 56
Bode diagram of a/(a+s) .................................................................. 60

vii

4.2

4.3

4.4

4.5

4.6

4.7

4.8

4.9

4.10

4.11

4.12

4.13

4.14

4.15

4.16

4.17

4.18

4.19

4.20

Filtering sinusoidal signal ................................................................. 61

Filtering square signal ..................................................................... 62
I vs. h (Duffing) ............................................................................ 63
I vs. a (Duffing) ............................................................................ 64
Pseudo phase space of the Duffing system ............................................. 65
I vs. h (Oscillator 1) ........................................................................ 66
I vs. a (Oscillator 1) ........................................................................ 66
Separation of the stick-slip signal ....................................................... 67
Pseudo phase space of [x,Lx] for the first oscillator ................................... 68
True phase space of the first oscillator .................................................. 68
I vs. h (Oscillator 2) ........................................................................ 69
I vs. a (Oscillator 2) ........................................................................ 7O
Separation of the stick-slip signal for the second oscillator (velocity) .............. 70
Separation of the stick-slip signal for the second oscillator (displacement) ........ 71
Pseudo phase space of [x,Lx] for the second oscillator ............................... 72
True phase space of the second oscillator ............................................... 72
FNN vs. d and h (Duffing) ................................................................ 73
FNN vs. d and h (Oscillator 1) ............................................................ 74
FNN vs. d and h (Oscillator 2) ............................................................ 75

viii

CHAPTER 1

INTRODUCTION

1.1 Thesis Statement

In this thesis, two alternative methods in phase space reconstruction are set
forward and studied: a fractional derivative method and a filter method. These methods
are preposed as supplements to current toolbox of phase-space reconstruction methods,
dominated by the method of delays. The method of delays is powerful, yet known to
introduce distortions of the phase-space, and have problems in discontinuous systems.
The study of this thesis shows the two alternative methods can be effective phase space
reconstruction tools. Whether the fractional derivative method is advantageous is not yet
clear, but the door is open to future investigation. The filter method is shown to improve
reconstructions in stick-slip systems.

In the rest of this chapter, we provide background on the phase-space
reconstruction of nonlinear systems, and also on the fast Fourier transform (FFT) as it is

relevant to our work.

1.2 Nonlinear system and dynamics
We may categorize all existing systems as linear and nonlinear systems. The
essential difference is that linear systems can be broken down into parts, solved
separately and recombined to get the answer. Thus standard tools such as normal modes,

Laplace transforms and Fourier analysis can be applied to simply the linear systems.

But many things do not have the property of superposition. For example,
listening to two favorite songs at the same time does not double the pleasure you get! [23]
The familiar tool of Fourier analysis fails to simplify nonlinear systems. The Fourier
transform of a linear system changes differential equations into an algebraic problem,
which makes the analyzing process much easier. Fourier analysis of nonlinear systems
turns differential equations in the time domain into convolutions among the Fourier
transforms of the independent variables in the frequency domain, which is not an
improvement. Therefore, the methods to analyze linear and nonlinear systems are
substantially different.

Also we can categorize all existing systems into known and unknown systems. A
known system has known equations of motion (ODEs or PDEs) are known. For these
systems, we know how many active states they have and what the states are. When we
study the dynamics of known systems, we calculate the state vectors by computer
simulation or by analytical methods, given initial conditions. This way we can get phase
space vectors and apply additional nonlinear dynamics methods.

Unfortunately, models are not exact. Moreover, there are many unknown systems
people are exploring. For these systems, we do not know the equations of motion.
Therefore, we cannot get state vectors from computer simulations or analytical methods.
The only thing we have is experimental data. In most-cases, we need to use lots of
sensors to get all the state vectors of one system. But in practice, typically, we cannot
measure each state in an experimental system. In fact, it may not be known what the
states are and how many active states there are. Therefore, for nonlinear systems,

experimentalists often make use of the phase-space reconstruction.

1.3 Phase space reconstruction
Here we use the term state space or phase space for the space in which orbits of
dynamical systems evolve. State space is more commonly used in the engineering
literature, while phase space is more common in the physics literature [34]. In this thesis,
we use both terms.

We assume the system can be modeled by continuous differential or discrete-time

equations. For the continuous case, we use u(t)=[ul (t),u2 (t),...,u, (t)] to denote the state

vector of the system. The equation of motion is then written as

du(t)

dt = 004(1)).

 

where G(u) is continuous and differentiable in its variables. For the discrete case, which
is the realistic situation, the equation of motion is written as u(n)=u( to +721), a map
from R’ to R’ , where z,is the sampling time interval [34]. The two cases are
considered very similar in the sense that either the discrete time series can approximate

the time derivative, as in

dug) u(to + (n+l)z,, —u(t0 + nr,)
dr 7 °

8

 

 

or the map can be on the Poincare section of the flow.

It is typically not possible to measure all of the (states in a system. The phase
space reconstruction is a representation of the full active phase space from a limited
number of measured observables, which in most cases is one single observable. The
reconstructed phase space is then called the pseudo phase space, compared to the original
true one. In other words, the goal of a phase-space reconstruction is to take, typically, a

single sampled measurement history and view it in a higher-dimensional space, which

can somehow accommodate the dynamics that generated the signal. The fully
reconstructed state space is useful in system characterization, nonlinear prediction, and in
estimating bounds on the size of the system [10].

Mathematically, the phase space reconstruction is an embedding process. An
embedding is a smooth map, (it, from the manifold M to a space U such that its image
¢(M)CU is a smooth sub manifold of U and that d) is a diffeomorphism between M and
¢(M). In other words, the embedding of M in U is a “realization” of M as a sub manifold
within U [4]. In particular, the fact that the embedding gives a diffeomorphism hem/ecu
two manifolds means that we have an important prerequisite with which to set up a
differentiable equivalence relation. A general existence theorem for embeddings in
Euclidean spaces was given by Whitney [1] who proved that a smooth (C 2) m-
dimensional manifold (which is compact and Hausdorff) may be embedded in R2“.
This theorem is the basis of reconstruction techniques for phase portraits from time series
measurements proposed by Packard et a1. [2] and by Takens [3,4]

Several researchers first suggested that the full dynamics might be recovered from
a time series. Packard et a1 [2] and Takens [3] suggested that a phase portrait, equivalent
in some sense to that of the underlying dynamical system, could be reconstructed from
time derivatives formed from the data.

Suppose we have scalar measurements s(n)=s( 10 + n1,) . If values can be acquired

for the time derivatives of the measurement s(n) such as ii , etc., then we can plot the

phase space and find the relationships among the derivative vectors and the measurement.

The first-order derivative of the measurement is

ds(t) I
dt t=to+nr, ’

so» =
We need to approximate the derivative, since we only have discrete measurements

sampled every t, . Naturally the approximation would be

~ s(t0 -I-(n+1)‘z,)-s(to +1125)
1'

S

S(n)

 

This is an approximate value for the first-order derivative, which means there is some
error. When we apply the same method to second-order derivative, we will have an even

bigger error. Actually, the derivatives are just approximated by different combination of

s(n) and delayed measurement s(n+kz,) [34].

In approximately 1980 Packard et al. [2], Takens [3] and David Ruelle
simultaneously and independently introduced the idea of using time-delay coordinates to
reconstruct the phase space of an observed dynamical system. The Santa Cruz group [2]
implemented the method in an attempt to reconstruct the time derivatives of s(t). The
main idea is that we really do not need the derivatives to form a coordinate system in
which to capture the structure of orbits in phase space, but that we could directly use the

lagged variables s(n+T)=s( to +(n+T)z,), where T is some integer to be determined.

Then using a collection of time lags to create a vector in d dimensions,
y(n) = [s(n), s(n + T), s(n + 2T),..., s(n + (d -1)T)] ,
we would have provided the required coordinates [34].
Takens [3] provided theorems to justify the use of the method of delays and the
method of derivatives in smooth systems. Takens’ theorem for MOD is summarized as

follows.

Let M be a compact manifold of dimension m. For pairs (F,v), F a smooth (i.e.
C 2) vector field and v a smooth function on M, it is a generic property that

¢F,v (Y) : M -—> R2” , defined by

(0;... (y) = (V(y).V(¢1(y))....,V((02,,.(y)))T
is an embedding, where (a, is the flow of F. Here we can think v(y) as the m-dimensional

measurement of the system.

The concept discussed here is illustrated in Fig. 1.1 (Broomhead and King [4]).

(1)
True phase space

 

‘ V r (2)observable(1-dimension)

N f‘.

 

 

  
 

(3)
pseudo phase space

Figure 1.1 Schematic representation of embedding process
In Fig. 1.1, (1) is the true phase space of the dimension m, (2) is the measured
observable, which in most cases is in Ldimension, and (3) is the embedded pseudo phase

space in the dimension n.>. 2m+1. The measurement (2) is a function of the true state

variables in (1). In most cases, what we have hand is the measurement (2). The
embedding (P is the process from the true phase space (1) to the pseudo phase space (3)
via the observable v(y).

MOD was further refined relatively recently (Sauer et al., 1991). So far, time-
delay embedding is, perhaps, the only systematic method for going from scalar data to
multidimensional phase space, and it has been the most explored in the literature.

In MOD, the problem is how to choose the lag T and the dimension d. If the time
lag T is too large, the delay coordinates may be statistically uncorrelated to each other
and the pseudo phase space will be random. If T is too small, the delay coordinates may
be very similar to each other so that the pseudo vectors tend to lie on a diagonal in the
reconstructed space.

Average mutual information and correlation methods can be used to choose a
good delay. The average mutual information is really a kind of generalization to the
nonlinear world from the correlation function in the linear world.

Many methods are put forward for determining the embedding dimension d: the
singular-value analysis, the saturation of system invariants, the false nearest neighbors
and the true vector fields

Also, Broomhead and King [4] suggested the use of a form of singular value
decomposition (referred to hereafter as singular systems analysis or SSA) as a means of
overcoming many aspects of MOD. Use of a basis obtained by truncating the full SSA
basis to the first N most significant vectors (i.e. with the N largest eigenvalues) is claimed
to offer a number of enhancements to MOD, including an element of data-adaptive noise

reduction and filtering, provided N is large enough to capture the deterministic dynamics.

Gibson et al. [5] demonstrated that the relationship between delays, derivatives
and principal components consists of a rotation and a rescaling. Consequently, from
Gibson’s point of view, statements about the nature of the equivalence between the
original and the reconstructed phase portraits would not depend on the coordinate system.

To determine embedding dimension, Kennel et a1 [6] put forward the method of
false nearest neighbors (FNN). Cao et al. [7] proposed the method of averaged false
nearest neighbors.

In the paper of Sauer et al. [8], mathematical formulations of the embedding
methods commonly used for the reconstruction of attractors from data series are
discussed. Embedding theorems, based on previous work by H. Whitney and F. Takens,
are established for compact subsets A of Euclidean space R" .

In 1998, Gilmore [35] presented an integral-differential phase space embedding.

He suggested using the vector (x, 2,3: ,. . .) or the vector (I x, x, in...) for reconstruction. He

pointed out the serious weakness: noise will be amplified a lot when high order
derivatives are used.

It is widely recognized that a time series of measurements from more than one
probe is more desirable if practical considerations permit. In the paper of P. L. Read [9],
multivariate SSA (M-SSA) is shown to offer significantly improved reconstruction in
terms of signal-to-noise ratio and structural uniformity 'of the attractor over earlier

methods.

1.4 FFI‘ (Fast Fourier Transform)

The expression for the DFT is, setting the sampling time t, = 1for simplicity,

X (k) = ”gimme-1'2“" ’”

n=0

The inverse (IDFT) is written:

N-l
x(n)=2 x(k)e+j2flnk/N
k=0

Since straightforward evaluation of DFI‘s is excessively time-consuming, the fast
Fourier transform (FFI‘), which is a family of algorithms for computing DPT, provides a
means of greatly speeding up DFI‘ computations.

As we mentioned before, Fourier transform is an effective tool in linear system
analysis. The Fourier transform of a linear system changes differential equations into an
algebraic problem. In nonlinear system analysis, Fourier transform may also be useful,
but not enough to analyze the complicated nonlinear systems.

Measured scalar data is always contaminated by noise. For linear and some
nonlinear systems, we can use Fourier methods to identify and remove noise. However,
in some nonlinear systems, the measured quantities are irregular or chaotic in time
domain. These time series have been considered bothersome “noise” in the sense they are
both irregular and broadband in frequency domain. Fourier analysis will not be of much
assistance in making the separation.

The spectrum of a sine wave of infinite duration peaks at a single frequency. But
if the FFT is applied to the sine wave over a finite interval, the single spectral peak is
spread into a range of frequencies. This is called spectral leakage. Leakage is due to the
finite data record, or window, which effectively multiplies the signal by a rectangle and

thus convolves each spectral peak with sin(f)lf. The peaks are broadened, making spectral

separation of nearby frequencies difficult, and the sidelobes can obscure small-amplitude
spectral peaks at other frequencies. For example, suppose the measured scalar data is
x(k), where k=1,2,...,N. Then we apply FFI‘ to x(k). Because x(k) is in finite length N,
the PET result corresponds to x(n)*o)(n), where n=1,2,...,oo and (s(n) is a rectangle
window function, u)(n)=1 for 15 n S. N , and 0 otherwise. Let x1(n)=x(n)* (s(n), then
X1(ju))=X(iO)) ® W00», i.e., X1(ju)) is the convolution of X00» and W(iu)).

A signal observed for a finite interval of time (window) may have distorted
spectral information in the Fourier transform due to the ringing of the sin(f)lf spectral
peaks of the rectangular window. To avoid or minimize this distortion, a longer signal
record can be used or a signal is multiplied by a window-weighting function before the
DFT is performed. Window choice is crucial for separation of spectral components that
are near one another in frequency or where one component is much smaller than another.

There are some popular window functions: the rectangle window, triangular
window, Gaussian window, Cosine bell window, Hamming window, Blackman window,
and Kaiser-Bessel window.

Here we take an example of a sinusoidal signal x(t)=sin(6t+l), with t0=0.

Fig. 1.2 shows the no-leakage case, where the maximum time tw = 212 *100/ 6 and
sampling interval “:3:th IN (N=10,000). In this figure, the solid line is the original

sinusoidal signal x(n), the dotted line is the derivative lx(n) calculated by FFI‘ and IFFI‘,
and the pointed line is the derivative dx(n) calculated analytically. In this case, there are

exactly 100 periods in x(n). Thus there is no leakage, and lx(n) and dx(n) overlapped.

10

x, Ix and dx(solid—x. dotted-4x, pointed--dx)

 

 

 

 

Al'-

4 :L. ...... J .......... a . . _ .

. . . . ooooooo J nnnnnnnn ho .

r ........ . ....... k ...... g----'-'- --.'--.' .“H””'.
II. . o . L- IIIII ” IIIII 0"“..I.

. . ............ s IIIIIIII a--- . .

. IIIIIIII Wat . . . . .

..r iiiiiii .— . . . . .

. . IIIIIIII 1 nnnnnnnnnnn . . .

. . . o tIoJ tttttttt Lilia .
I‘ ....... .1 UUUUUUU J--- -J ........ J ........ ‘8---““I'hl

. . . . nnnnnn r uuuuuuuu a sssssss .

. nth iiiiiiiiiiiiiiii «it . . .

P llllll c u o u c u

. tttttttt r ooooo . . . . .

. . 00-1 tttttttttttttttt .Ioo . .
If ........ r ....... L ''''' ah ........ ‘-.8I."I.’“" """"" "'I

. . . . . . tutu;

. - '-. ........ ’ ........ ‘-‘-- .

. sssss L nnnnnnnnnnnnn a . . .

«Ho . . . . . .

. iiiiiii f tttttttt . llllll . . . .
l.‘ ....... n. ....... JD----.-I ”nuunuuuu..l....-h'-...""".1
. . . . . . tttttttt H
. . ...... o. nnnnnnnn a ........ u .

. ccccccc T ooooooo J \ . . . .
“Ill! . . . a c .

. our-J iiiiiiii L rrrrrrrr . . .
Ir ........ r ....... L ....... ""--u.‘ ........ ‘- ..... .1
. . . . . . IIIIII d
c . . a o IIIIIII P llllllll a
. . tttttttttttt a IIIIIIII no . .

. lllllll JI ll . g c . .
i ttttttt . . . . . .

c IIIIIIIIIIIII
U‘ ........ u ........ “DUDDHIIN . ....... dc......."”"" ...... «I

. . . . . . turn“.

. . . . IIIIIIII e lllllll .

. o IIIIIIIIIIIIIII J llllllll .u . o

a llllllll u . a o c a

an llllllll .llll . o a c c
I' '''''''' TI--.-..‘ ......... “..””"J" ....... O ........ .l

. . . . iiiiii o rrrrrr .

. . . . . . tun.

. . M iiiiiiii L IIIIIIII 4 IIIII .

. llll. IIIIIIIIIIIIIII o u u .

NHHTI . . . . . .
I".-D'HHHHA".'.."'J."..’”U.II cl..." ..... * UUUUUUU * ....... *1

. . . . a. iiiiiiii 1 IIIIIIII .

. . . . . r* tttttttt a

o . llllll L llllllll L lllllll .

o IIIIIIII cl IIIIIIII o . c .
Put . . . o . .
1'-|"~“""I"".l‘I..I...I..'I' ....... ‘ ....... . ........ .l
. . . nnnnnnnnnn r tttttttt f .

. . . . . . ttttttt .

. . . . . i talk tttttttt J
. . nnnnnnnnnnn f tttttttt .000 i . .

n rrrrrrrr 1 rrrrr . . . . .
I'KI'I'II . ....... o ....... c ........................ l
H u -I--.' ............. *1 ...... w ..... w
. . . . . . ”nu!

. c . - I. llllllll LIIOII .

o h IIIIIIIIIIIIIII Q IIIIIII . u .
IIIIIIIII . g g . .
I? """" O. ......... 1'... U‘ ........ ‘ ........ . ....... .1
I q ...... ‘ IIIIIIIIIIIIIIII u c u
. . . . 5 cccccccc T aaaaaa .

_ _ _ _ _ lurpr ....... J
5 4 2 U 7.. 1.4 5

.
.6 cc... x. x

97 98 99100

92939495%

91

t

Figure 1.2 Derivative by FFT (no leakage)

Here we choose another value for the frequency so that the total time isn’t integral

multiplication of periods. Then there would be some leakage. Fig. 1.3 shows the leakage

100 and sampling interval ts=tml N

case, where the maximum time t1m

(N=10,000). The solid line is the original sinusoidal signal x(n), the dotted line is the

derivative lx(n) calculated by FFI‘ and IFFT, and the pointed line is the derivative dx(n)

calculated analytically. We can see leakage causes distortions in the beginning and the

end of the differentiated signal.

11

x. Ix and dx(solid—x. dotted-4x, pointedudx)

 

 

 

 

 

 

- I , ........ I ........... r ..................... I .......... j. _

5 P i,»“””:' M... : : i z :

,c'. : N1; 2 2 i I

4--s ----------- a ----------- a "N; ------ ----------- a ----------- r-

: i i "c . l i i

I I I fl. I I I I

2 —-r ----------- : ----------- : --------- '- ;.'- —-: ----------- : ----------- r~

: ' . a. : w:

' : : : 3, : : :

o _-; ---------- 4 ----------- r ---------- . ----------- :—

i 2 i i ‘q 1 I i

x : : : : *-,: : :
3 '2 "I """""" : """""" r """"" r """"" a; """"" : """""" r“
C I I I I I t, I I
a : : : : : '3... : . :
—’-‘ .4 —-: ---------- 4 ----------- 'r ---------- i- ---------- :------' .4": -------- 4»;-
" I I I I I s I I' I
X I I I I I '3" I a f. 1" I
: : : : : in." .~,:; :

.5 --s ---------- ~: ----------- ;. ---------- ;. ---------- ~: ----------- aura-ere

s a s a a a i": s

'8 '"i """"" 1 """""" I”""""'.' """"" ‘. """""" '2' """" 1'?“
-1o—-s ---------- ~: ----------- ;. ---------- ;. ---------- ~: ---------- ~:— -------- -‘:-:~
-12--i ---------- 1' ----------- 'r ---------- I ---------- i ---------- + ---------- 'r
994 995 996 99.7 998 999 100

t
Figure 1.3 Derivative by FFI‘ (leakage)
1.5 Overview

A proposed fractional derivative method is tested on sinusoidal signals first, and
then applied to three systems: the Duffing oscillator, the Lorenz system and a two-well
experimental system. The derivative leakage problem caused by fast Fourier transforms
(FFT) is studied and a criterion is developed to measure the distortion of a one-derivative
signal calculated by the FFT method. Based on this criterion, a simple but effective
method—a truncation method—is put forward. Fractional derivatives are then used as
pseudo phase vectors. The embedding dimension, unstable periodic orbits and correlation
dimension are studied for the three systems.

A filter method is proposed for problems with stick-slip. The method is tested on

sinusoidal signals first, and applied to three systems: the Duffing oscillator and two

12

Coulomb oscillators. Optimal values for the filter-parameter are determined. The filtered
signal and delayed signal compose the pseudo space veCtors. The embedding dimension
is computed by the false nearest neighbor (FNN) method.

This thesis is organized as follows: Chapter 2 introduces some systems used in
this thesis. Chapter 3 focuses on phase space reconstruction by means of fractional
derivatives. Chapter 4 addresses a filter method for phase space reconstruction.

Conclusions are drawn in Chapter 5.

13

CHAPTER 2

SYSTEMS AND TOOLS

2.1 Introduction
In this thesis, we develop reconstruction techniques and test them on chaotic
systems. We choose well-known models to test the fractional derivative method, and
stick-slip systems to test the filter idea. These reconstruction methods are tested using
standard nonlinear time series analysis tools. The systems and tools used are introduced

in this chapter.

2.2 System simulation
For a given system, if we know the equations of motion, then the system can be
simulated by computers and thus we get the state vectors by numerical integration.
In this thesis, data of the Duffing system, Coulomb friction oscillator systems, and
the Lorenz system are from computer simulation. Data of a two-well oscillator is from

experiments conducted by Joe Cusumano at Penn State University.

2.3 Duffing system

The Duffing system in this thesis is given by:

Jr+arx+x+x3 =flcosax ,

with parameters Ot=0.2, 8:27, and (0:133. [10]

14

The sampling time interval is At=0.053 and there are 94.4 points per driving

period.

 

 

 

 

Fm

Figure 2.1: Duffing oscillator system

Physically, there are many Duffing oscillators corresponding to the equation
above. One simple example (shown in Fig. 2.1) is a unit mass attached to a nonlinear
spring with restoring force R(x)=— x — x3 and external force F(t)=Bcosu)t. The spring gets
stiffer as the displacement x incmases——this is called a hardening spring.

The Duffing oscillator is a nonautonomous system for which a multitude of

nonlinear dynamics analysis and experiments have been applied [26,27].

2.4 Coulomb friction oscillator system (oscillator 1)
This is an example of a stick-slip chaotic oscillator. When stuck, the sliding mass has
zero velocity, and hence the displacement is momentarily constant. This system is prone

to reconstruction failure when the method of delays is used [10].

15

 

i . att)

 

 

 

 

° ' Noon
Figure 2.2 Coulomb friction oscillator 1

This model is shown in Fig. 2.2 [10]. Neglecting the transverse component of the

normal load, the nondimensional equation of motion is of the form

5c+ 2:5” x + n(x)f(x) = Acoth ,

where x is the displacement, I; is the damping ratio, and Acoth represents a harmonic

excitation. rm represents the coefficient of friction, and, using the Coulomb model, is
given by

f(x) = mign(1r),x is o, —1 s for) s 1,5: = o.
The normal load is given by n(x)=1+k,, x for x>-1/k,, , and n(x)=0 for x<—1/kn when the

contact is lost. If the model had a bilateral contact, there would always be a nonnegative
normal load, and the friction model would be such as in the oscillator studied by

Anderson and Ferri [28].

16

Here tanh(50x) is used in place of f(.r) so it’s a smooth system. The parameter

values are A=1.9, 9.=1.25rad/s, §=0, kn =1.5. The sampling time interval is At =0.053 and

there are 100.5 points per driving period.

2.5 Another friction oscillator system (oscillator 2)

This stick-slip chaotic system also leads to reconstruction failures when the
method Of delays is used. When stuck, the mass rides on the belt with constant velocity in

contrast to the zero velocity of the previous oscillator.

 

 

Figure 2.3 Coulomb friction oscillator 2

This system is shown in Fig. 2.3. It is a forced, belt-driven oscillator studied by

Stelter and Popp [37]. The equation of motion is:

17

Sc+x+x, +f,(v,) = Ycos(Q.t)

Here, v represents the belt velocity, v,=x-v is the relative velocity between the
mass and the belt, and Y and .Q is the amplitude and frequency of base motion. The

friction model used was of the form f,=Nu( v, )sign( v, ), where N is the ratio between

)+ a, + av,2 . Finally,

 

 

the normal load and spring stiffness, and |.I(v,)=(;10 — [(1)/(1+ ,1 v,
x, = ,u(v)N is the static equilibrium position. The parameters used in this study are

,uo=0.4, p,=o.1, A=1.42 s/m, 0:0.0132/m2, N=10m, v=1.0m/s, Y=0.05m, and

9:2.1625 rad/s.

The sampling time interval is At=0.03s and there are 97.4 points per driving

period.

2.6 Lorenz system

This is a paradigm for chaos in an autonomous system. The equations of motion

i=0(y-X)

y rx - y - xz
a = xy - bz
Here 0’, r, b>0 are parameters. Ed Lorenz (1963) derived this three-dimensional

system from a drastically simplified model of convection rolls in the atmosphere. A

mechanical model of the Lorenz equations was invented by Willem Malkus and Lou

18

Howard at MIT in the 19705. The simplest version is a toy waterwheel with leaky paper

cups suspended from its rim, as shown in Figure 2.4 [23].

 

Figure 2.4 Lorenz system example (waterwheel)

Water is poured to the paper cups from the top and leaked from bottom holes.
When the water input flow rate is too low, the top cups cannot hold enough water to
overcome friction, and the wheel remain static. When the flow gets faster, the top cups
gain enough water to overcome fi'iction and the wheel starts to rotate (figure a). At last
the wheel rotates steadily in one direction or the other (figure b). The probabilities of both
directions are equal and it depends on the initial condition.

By increasing the flow rate further, stability would be lost. Some cups may be so
heavy that they force the wheel to slow down or even reverse the rotating direction
(figure c). Then it may spin in the other direction for a while. The motion becomes
chaotic [23].

We choose the parameters o=lO, r=28, b=8/3. Sampling time interval is At=0.01 s.

19

.I
I.

 

rF‘

2.7 Two-well system

 

Figure 2.5 Two-well system example

A two—well system is shown in Fig. 2.5 [29]. The magneto-elastic beam
experiment was developed at Penn State University. The experiment consisted of a
stiffened beam buckled by two magnets. Two permanent magnets were placed on the
base of the frame holding the beam to create the two-well potential. The uncovered
portion of the beam near the clamped end acted as an elastic hinge from which the
position of the beam was measured by strain gauges. The frame was then fixed through a
rigid mount to an electromagnetic shaker. A periodic driving signal was fed through a

power amplifier to the shaker to provide the external forcing function.

20

The driving frequency was set at 7.5 Hz and data was collected at the sampling
frequency of 187.5 Hz for 7000 periods of excitation. At this sampling frequency and

driving frequency, there are 25 samples per driving period.

2.8 Average mutual information
One of the nonlinear dynamics tools we will use in Chapter 3 and Chapter 4 is the
average mutual information (AMI). The AMI can be used to measure nonlinear
independence of two reconstruction coordinates. This is computed by plotting the two
coordinates, and dividing the resulting two-dimensional pseudo space into square bins.
Denoting A as the set corresponding to one axis, and B as the set corresponding to the
other, we can then set event a to represent a bin in A, and event b to represent a bin in B.

Then PA (a)is the probability of a in A, which is the number of data in bin 0 divided by
the total number of data, and P, (b) is defined likewise. PM (a,b) is the joint probability,

which equals the number of data that are in a and also in b, divided by the total number of

data. In other words, the bins form a grid on the two-dimensional plot. Bin 0 is a vertical

strip, and bin b is a horizontal strip. PA (a) and P, (b) are obtained by counting data in
the vertical and horizontal strips, and PM (a,b) is obtained by counting data in the

intersection of the strips. Given these probabilities, then

PM (a,b)

[AB (a,b) = logz m
A B

The average mutual information is then obtained by averaging this quantity over

all the bins [10].

21

In the MOD, the delay h is chosen as the lowest value which produces a local

minimum in the AMI between x, and xM. We note that the full reconstruction of
dimension d will also include elements xi,” , k=2,3,. ..,d-1. It is possible that the average

mutual information between x, and x“, is not minimal.

2.9 False nearest neighbors
Here is another tool for nonlinear dynamical analysis. The false nearest neighbors
computation [Kennel and Abarbanel] is a method of determining if there are false
crossings of trajectories. In the full reconstruction of the phase space, trajectories should
not cross. When the reconstruction dimension is too small, the phase flow is projected
onto a smaller space, and trajectories cross in the projection. In the FNN computation, the
nearest neighbor of each point in the reconstruction is determined. If the distance between
the point and its nearest neighbor increases by more than a factor of 10 as the
reconstruction dimension is increased, then a false crossing of trajectories is deemed to
occur, and the nearest neighbor is said to be false. For a reconstruction dimension for
which there are no false nearest neighbors, it is said that there are no false crossings of
trajectories. Hence, the reconstruction is one-to-one with the true phase space, and the

reconstruction is an embedding.
Since a reconstruction is built from pseudo state variables (i.e., delays, fractional
derivatives, or other operations on the observable), the embedding dimension usually

differs slightly from the dimension of the true state space. For example, in the Duffing

system, the 3-D phase space geometry is 9i’ x S1 for (x,5r and t mod T) coordinates, and

the delay reconstruction tends to be 4-5 dimensional.

22

2.10 Unstable periodic orbits
A chaotic invariant set is the closure of infinitely many unstable (saddle-type)
periodic orbits. A chaotic trajectory continually “visits” these UPOs. The UPOs provide a
“skeleton” which characterizes the system topologically. They can be used for the
estimation of Lyapunov exponents [39, 40], computation of fractal dimensions [41, 42],
controlling chaos [43], template analysis and knot theory [44, 45,46], evaluation of signal
deterrninacy [34], and the identification of nonlinear systems [47] and system parameters

[48, 49, 50, 30].

2.11 Correlation dimension
The correlation integral can be calculated by counting number of pairs of points

that have a distance smaller than r, that is,

C(r) = N1? * (number of pairs (i,j) with x(i,j)<r), or can be put:
1 N N
C(r) = F2 290- | x(i) — x( j) l) , where N is the number of points, and 0(.) is
is] j=i+l
the Heaviside step function, i.e.,

0(r-|x(i)-x(j)|) is 1 when r>|x(i)-x(j)| and 0 when r<|x(i)-x(i)|.

Then the correlation dimension is dam = limw. (see, e.g., [34,43])

H0 log r

23

CHAPTER 3

PHASE SPACE RECONSTRUCTION BY MEANS OF

F RACTIONAL DERIVATIVES

3.1 MOTIVATION

For most applications, the MOD is a dependable way of producing the pseudo
phase space. While points on a true phase trajectory represent an instant in time, the
MOD produces d-dimensional pseudo phase points that involve a time span of (d-1)T.
This “smearing” of time can cause distortions in the reconstructed phase space [31]. Here
we try reconstructing the phase space without incorporating delays.

In this chapter, we will try fractional derivatives as the reconstruction vectors.
Packard et al. [2] first suggested using derivatives as the reconstruction vectors.
However, the traditional derivation method has a big problem--noise amplification. For
example, the observable displacement vector x(t) has some noise n(t) when measuring,
where n(t) changes abruptly but has small amplitudes. Because of its small amplitude, the
noise n(t) can be omitted in most cases. But when we use derivatives (velocity,
accelerate, etc.) of the displacement from the contaminated data x(t), the noise
component might be so big that the resulting derivatives we get is far away from the real
derivative values and the derivative data become meaningless. To reduce the noise
amplification but at the same time generate a sufficient set of independent coordinates for

a phase space reconstruction, we use fractional derivatives.

24

As is known, in the time domain, the first-order derivative of some variable x(t) is

. .. 2
= E, and the second-order derivative is x = i; =_d_(_d_x) , and so on. We use the
dt dt dt dt

Operator D"x for the derivatives x, x,x, where a is the order of the derivative.

Then D'x=x, D2x=x, In calculus, we have only integers for a. What if we use a

fractional number, like V2, for a instead of an integer? We define fractional derivatives as

 

the operators D" x= where a is a fractional number.

For ordinary derivatives (a is an integer), we surely know how to do it (many
engineering people do it every day). However, calculating fractional derivatives is
different. It seems very difficult to do it in the time domain. Thus we do it in the
frequency domain.

The procedure is:

1) FFI‘ of x(n) to Xfiw)
2) Multiply X(j(n) by (jw)"a to get F0m)=FFI’( D“ x)
3) Inverse FFI‘ of F000) to D“ x

The quantities [x,D"x,D2"x,...,D“’”"’x] will be used to build the reconstruction
vectors.

We still have noise amplification. But when we choose a “small” value for a,
noise amplification will not be a “big” problem. For example, if we want to use five
derivative vectors to reconstruct the phase space, we would use the vectors
D'x,sz,D3x,D‘x,D’x if integer derivatives were used. D5 x is used and it might have

a very severe noise amplification problem and become meaningless. Now we could use

25

fractional derivatives, say, the vectors D°"’x,D°'8x,D"2x,D“°x,sz. For vectors of the

like dimension, the applications of fractional derivatives will make the noise
amplification much smaller than that in the traditional case.

The next step in the reconstruction is to normalize each axis of the data, such that

x, D‘x,

R.’ R.

 

Y=[ ,...] where R j is the span of the unnormalized coordinate R j =max D’“x,.-

min Djaxi. The reason for doing this is for the time series analysis metric methods that

involve finding data with small “balls”.

 

Let us examine the Bode diagram of s“ (d in the time domain):

dr"
=2010g|G(iw)|=2010g|(ju))"a|=20a*logo),

Thus the Bode diagram of a fractional derivative term is a straight line that intersects the

(0 axis at (0:1 with a slope of 20a per decade. The Bode diagram of the a=l case ( s' = s)

is shown in Fig. 3.1.

When we choose a small value for a, the slope 20*a would be small too.
Compared to traditional derivatives (s, s2, s3 ), it is better in respect of noise amplification

though it is still a high pass filter.

26

 

40 Y I 7 UVYVVU V I IrlYVIU
D ..II.

 

 

 

I
I
I
I
I I I I I I I II I I I
I I I I I I I II I I I I I I I I I I I I I I II
I I I I I I I II I I I I I I I I I I I I I I I I
30 -------------------- I ............................ ' ................... 'uLl—l
I I I I I I I II I I I I I I I I I I I I I I I I
I I I I I I I II I I I I I I I II I I I I I I I I
I I I I II I I I I I I I I I I I I I I I I I
I
I
I I
I I I I I I I II I I I I I I I I I I I I I I I II
I I I I I I I II I I I I I I I I I I I I I I I I I
20 ..................................................................... -'—LA_.
I I I I I I I II I I I I I I I I I I I I I I I I
I I I I I I I II I I I I I I I I I I I I I I I I
I I I I I I I II I I I I I I I I I I I I I I I I
A I 0
g ' I
I
U . g
I I I I I I I II I I I I I II I I I I I I I I I
10 ................................... -J--0--|-AJ-'-fl.------k---L--L-A-L-'-LA_
I I I I I I I II I I I I I I I I I I I I I I I I
I I
I I
I I
I I
1 ' I
I I
O ' I
....................... ...; .1 .I .I .. .I I. 5 -I. h.
I I
I I
I I
I I
I
I
I
-1o -
...................................................................... A—
I
I
I
I
I I I I I I I II I
I I I I I I I II I I I I I I II I I I I I I I I I
-20 I I I I I I I II I I I I I I III I I I I I I I I
A L L A A A AA A A A A AJAA A A A L LLA

F roquoncy (rodboc)

Figure 3.1 Bode diagram of s" (a=1)
3.2 Test on sinusoidal time series

Before applying the fractional derivative method into actual systems, we tested it

on sinusoidal time series. Here we take x(t)=sin(6t)+l with the sampling time interval

1, =0.0ls.

Below, we analyze the procedure in page 25 for the sinusoidal signal for the case

of a=1. Applying the Fourier transform to x(t), we have

X(jw)=-;- j[b'(w+ 6) - 6(a) - 6)] + 6(a))

Then, the Fourier transform of the derivative of x(t) is

F(jw) = X(jw)* H(ja)) = -%w[6(w+ 6)—5(w-6)]+ ja)6(a)) = 36(w-6)+36(w+ 6)

Inverting the transform, we have

Dx(t)=6cos(6t)

27

The result Dx(t) is the same as the traditional first order derivative

3:1:- = %(sin(6t) + 1) = 6cos(6t) . Fig. 3.2 shows the fractional derivative signals Dx(t) for

different a values.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

' :i? E: '1

g ', .1 r: w .4 i

'05 "‘ 4 I A ' ':: '1 :I' :I
40 Q 44 Is 48 so '40 42 u Is 48 so

 

BLILII'o-‘Io-Inwu
abb‘t'JDmem

‘2 M 5 48 5]

Figure 3.2 Fractional derivative Lx(n) and original signal x(n) (solid—x(n),

dotted—Lx(n))

3.3 Leakage problem

The derivative is based on FFI‘s that are prone to leakage, caused when the signal
of total time (length of the signal) is not periodic with period T (i.e., the endpoints do not
match up). In the time domain, leakage causes distortions in the beginning and the end of

the differentiated signal.

28

As we said in Chapter 1, a window method can reduce leakage. The most popular
window functions are rectangle window, triangle window, Gaussian window, Cosine bell
window, Hamming window, Blackman window, Kaiser-Bessel window, etc., which can
be found in reference [11]. Two typical window functions are studied here: near-
rectangle window and Kaiser-Bessel window.

The rectangle window and its frequency distribution are shown in Fig. 3.3. Fig.
(a) is in the time domain, and Fig. (b) is the Fourier transform. The window function is
shown as N-point samples with N=50, which sets the scale for the frequency axis
6 = 27zk/N. The range, —71 <6 <72 , corresponds to the Nyquist range,

— f, /2 < f < fI /2 , or the dimensionless frequency range — N/2 < k < N/2 [ll].

 

 

 

 

 

 

 

 

 

 

It 0:13
1 125 H .20
‘ 100 ,. l l
H W ‘. "’ ll ”"1“.
|| . i
|| ..
ll .
-25 -2o .15 -1o .5 o 2 10 15 20 25 -'7r I I I ' 3 I I I fig
3. Rectanglewindow. b. Log-magnitude(rectangle). I

Figure 3.3 Rectangle window (from reference [11])

In order to make the rectangle window smoother, sinusoidal functions are applied
at both ends. The resulting window, dubbed here the near-rectangle window, is shown in

Fig. 3.4.

29

 

 

0.9 -
0.0 -
0.7 -
0.5 ~
0.5 - ,
0.4 -
0.3 -
0.2 ~
0.1 ~

 

 

1

 

0100020003ti30400050in600070000000900010000
Figure 3.4 Near-rectangle window

It
I°[”a)/l_(N/2)2]
The Kaiser-Bessel window is defined as w(n)=

I 0 (7m)

 

, where

.. k 2
10(X) = 2 (XI/('2) . This window and its frequency distribution are shown in Fig. 3.5
k-O -

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

[11].
o 12.5
m 1.00
{I
fiqilllflllmu 1 1““!I'llrhjn ' V '
-252015-10-50510152025 -x o. a
m. Kaiser-Bessel window. 0. Log-magnitude(l(aisIr-Beasel.a-2.5).

Figure 3.5 Kaiser-Bessel window [11]
The amount of leakage is vague. In linear signal analysis, windowing is very

effective since it reduces distortions found in the frequency domain. However, nonlinear

30

signal analysis tends to be geometric, with tools often applied in the time domain.
Applying windowing functions across the time record can distort the geometry of the
signal. We expect it to have detrimental effects on nonlinear prediction, unstable periodic
orbit extraction, and fractal dimension computation. To measure signal distortion, we
compute

= max I Lx(i) - Dx(i) I ._

ERR . ,1—l,2,...,N
max I Dx(r)|

 

where Lx is the derivative calculated by FFI‘ with leakage, Dx is the true derivative
without leakage and N is the length of the sampled series Lx and Dx.

To reduce the leakage error, here we put forward a simple but effective method—
truncation method. Truncation method means to chop off some points at both ends. For
example, Lx(n) has 10,000 points, we can truncate N, =500 points at both ends with 9000
points remaining in Lx(n).

The leakage error ER is calculated for the near-rectangle window, the Kaiser-
Bessel window and the truncation method for the sinusoidal signals. The results are

shown in Table 1. Here “Original” means the case before truncation.

 

 

 

 

 

 

 

 

 

 

 

 

ERR*max[Dx(i)| a=0.2 a=0.4 a=0.6 a=0.8 a=l.0 a=l.2
Original 0.9289 2.5452 4.9770 7.8754 9.0585 20.5177
Near rectangle 2.1746 2.7021 3.4154 4.5277 5.9930 9.2721
Kaiser-Bessel 1.6612 2.1549 2.9677 4.1960 6.0088 8.5697
Truncation 0.0102 0.0107 0.0084 0.0070 0.0128 0.0357

 

Table 1: Leakage errors for different methods

31

 

For the results in Table 1, the true value of D‘x was estimated by adjusting the
sampling interval such that the sinusoidal fit perfectly in the time window. The sampled
D‘x signal was then converted to the sampling rate used in L°x by interpolation.

The truncation result in the table is achieved when the truncation number N, =500
and total number of points left is N-2* N, =10000-2*500=9000.

From Table 1, we can see that the truncation method can decrease the leakage
error ERR a lot more than the window functions. This is due to the fact that the time-
domain leakage distortions mainly exist in the beginning and ending points of the
differentiated signal. In the middle section, the leakage error is pretty small compared to

both ends.

In the following, a criterion is proposed for determining how many points to be
truncated, that is, how to choose an appropriate value for N, .
In Fig. 3.6, we plotted the leakage error y, = le, —Dx,| for the Duffing

oscillator, where Lx is the first-order derivative by the FFI‘ method (with leakage) and
Dx is the true valued first-order derivative from the system simulation (without leakage).
The total number of the data set is N=57344. (a) is for the beginning points l~ 250 and

(b) is for the ending points 57095~57344. Based on Fig. 3.6, we can select an appropriate

value for the truncation number N, , given an upper limit a for the leakage error. In this
case, for example, N, must be greater than 50 if 8:0.1.
(1) The Duffing system—we used N,=150 and the maximal value of [Dxl is

14.0720.

32

le-Dxl

(3deth

ILx-Dxl

 

 

 

 

100 200 300

(a)

 

 

 

 

J

 

'73?

5.71

5.72 5.73 5.74
" 3:10‘

0))

Figure 3.6 Leakage error of the Duffrng system

(2) The Lorenz system—we used N, =150 and the maximal value of [Dxl is

149.9532.

ILx-Dxl

500

400 r

300
200
1 00

 

 

 

1 00 200 300
00

33

 

500 - - -
400 »
300 »
200’ 4

100 i
- J
g7 5.71 5.72 5.73 5.74
n x 104

le-Dxl

 

 

 

 

(b)
Figure 3.7 Leakage error of the Lorenz system
(3) The two-well system—we used N, =300 and the maximal value of [Dxl is

3009.4.

x10

 

W

ILx-Dxl
O M h 0') CD

 

 

o 1 06 200 300

X113

 

M OJ #-

le-Dxl

 

 

 

[5.7 5.71 5.72 5.73 5.74
n x104

(b)

Figure 3.8 Leakage error of the two-well system

34

3.4 Look for a “g ” value of a
In the method of delays (MOD), one effective way to determine the optimal delay
h is the average mutual information (AMI) method [32]. We apply this method to look
for a good value for the fractional derivative order a such that the coordinate axes in the
pseudo space have optimal independence.
The AMI can be used to make sure D“ coordinates are “reasonably nonlinearly

independent”.

1) Duffing system
First, we compute the average mutual information I between the vectors x(n) and
x(n-r-h) for different values of delay h. Here x(n) is the original signal and x(n+h) is the
delayed signal. The result of I vs. h is shown in Fig. 3.9. In MOD, we choose the value of

h where I reaches its first minimum. In this case h=28 and I is about 0.22.

 

2.2

AA
man

.I
h

I (average mutual information) i
D .I
to _. M

.0
or

.0
I.

 

 

 

.0
M

 

h (delay number)

Figure 3.9 I vs. h for the Duffing system

35

Then we calculate the average mutual information I between the vectors x(n) and
L“ x(n) for different derivative orders a, where x(n) is the original signal and L" x(n) is
the differentiated signal. Different a will result in different L“ x(n), thus different I. Fig.

3.10 shows the result of I.

 

...)
(D

..s .3
b 01

.5
N

l (average mutual information)

 

 

 

 

0 0.2 0.4 0.8 0.8 1 1.2 1.4 1.6 1.8 2
a (fractional)

Figure 3.10 I vs. a of the Duffing system
The prescribed usage of the AMI is to choose the reconstruction parameter, in this
case a, at the minimum. The minimum occurs where a=l. Thus a full derivative leads to
coordinates that are most nonlinearly independent. It turns out that at 0:1, the value of
0.22 for the AM] is the same as the minimal value for delay coordinates. While we want
to choose a to minimize the AMI, we do not want a to be too large since it will amplify

noise.

36

 

 

Lindabrbonumrn

 

 

 

 

 

 

 

 

 

Figure 3.11 2-D pseudo phase space of the Duffing system (horizontal axis—x(n),
vertical axis— L“ x(n))

In Fig. 3.11, we plotted the two-dimensional pseudo phase space [x, L“ x] for

different a values. For comparison, the true two-dimensional phase space [x,rc] is plotted
in Fig. 3.12. From the figures, we can see that our pseudo phase space resembles the true
phase space, which gives us confidence to calculate other nonlinear characteristics such
as false nearest numbers (FNN), unstable periodic orbits (UPO) and correlation

dimensions.

37

 

 

 

 

 

Figure 3.12 True 2-D phase space of the Duffing system (horizontal axis—
displacement x(n), vertical axis—velocity y(n))
For comparison of the derivative method and method of delays, two-dimensional
pseudo phase space [x(n),x(n+h)] is plotted in Fig. 3.13 for different delay h values. The
delay h is selected such that the average mutual information between x(n) and x(n+h)

approximately equals to that between x(n) and L‘x(n) for different a values in Fig. 3.11.

38

 

 

11M?!) 2 x(mh)

 

 

 

 

 

 

 

 

 

2
“M, mm
0

 

 

 

 

 

 

   

Figure 3.13 2-D pseudo phase space of the Duffing system (horizontal axis—x(n),

vertical axis—x(n+h)

2) Lorenz system
For data of Lorenz system, we apply the same method for choosing a good value
for the fractional derivative order a. For the delay reconstruction, the result of I vs. h is
shown in Fig. 3.11. Again, we choose the value of h where I reaches the first minimum.

In this case, h=19 and I=0.85.

39

 

 

 

 

€235.25 .355 $233 _

N.

40

30

h (delay)

10

_ _ d _ _
. u n u u
- - - . n
u . . . a
g c n . o
o u . - n
. a - - -
. - . . -
g g o - c
- o c . c
- g o . c
v IIIIIIII ‘ IIIIIIIII ' ........ L IIIIIIIII '
- g o . n
o . o . .
u — c . -
- . . o
- - u n c
— . u c -
- c u c u
— - - . .
- o c . u
- . - . u
- c a . —
' ........ ‘ ......... ' ........ ‘ ........ 7
- n o . -
— c o . -
- u - c u
n o n c u
. o a . c
. g c . u
. o o . c
. u - . .
c c a . n
- - - . .
. - - . c
v IIIIIIII ‘ IIIIIIIII ' ........ ‘ IIIIIIII '
c a - g c
. a o . p
. . o . u
- . - . u
- . o c c
- . o u u
. u . - -
- c . - -
u o - g o
- u . c -
o - . c .
f IIIIIIII b IIIIIIIII f 00000000 L IIIIIIII f
c . o c -
. - a u u
. — . . u
- u c . u
a o . . p
u - c . n
. c . o p
- - . .
c o o a .
c g . _ .
. c . . u
I ........ r ........ ‘ ......... r ...... ‘ IIIIIIIII r
. u g . u
. c c . o
u o . . .
- c - .
. p . u .
u o o o -
o a u n —
. . c . u
n n . o
. . _ . .
g - u -
n p — h n
5 2 5. 1 5 o
2 1 0

Figure 3.14 I vs. h of the Lorenz system

Also, we calculated I vs. a, which is shown in Fig. 3.12. The optimal AMI is

0.8 at a=0.4. In this case, a fractional derivative leads to coordinates that are most

I

 

 

   

   

--'------‘------‘-------P------,------

   

-----pc-

     

'------‘------

 

 

nonlinearly independent.

. o g
- u .
o o g
a . n
o o .
IIIII+ IIIII r ..... . ..... Or ..............
. u -
n u u
- . n
. o o
u u -
. - -
IIIOI} ..... P IIIII fl IIIII .- IIIIIIIIIIIIIIII
. u .
a o .
o o c
o u c
c u u
u a o
u o a .
If...) IIIIII q ..... fl ...... c ................
o c o
o - .
a u u
c o u
. g g
. g -
IIII.* ..... ‘ ..... C ..... J. ...............
n o -
u c c
_ c c
c . c
u o —
o . o
'I...* ..... v ..... O ..... * ...............
u u o p
c u u o
g - u c
g u o .
u o o c
c . n -
IICOCLI ..... P IIIII F ...... n. IIIIIIIIIIIIIIIIIIIII
n . - .
- g g u
. . - o
c o u u
n u c a
c o u -
o .
I ..... C ...... * ..... * ...... n ..................
. - . .
. n o -
c o . .
c o - o
a g - -
- c - -
'I.II‘ ..... ‘ ..... ‘ ..... * .....
. . g
u c c -
. c n a
c c o u
p . - .
c o .
1|..I‘I IIIII V IIOI IIIIIIIIII ‘ .....
o u u u
n c u n
g g . o g
. o g u o
o o c n u c
n n g - . c
p h h n p n
2 a. e 4. 2. 1 a.
1 1 4| 4| 0

eofiées .355 omeoé _

0.4

1.4 1.6 1.8

1.2

0.8

0.6

0.4

0.2

0.2

a (fractional)

Figure 3.15 I vs. a of the Lorenz system

40

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.16 2-D pseudo phase space of the Lorenz system (horizontal axis—
x(n), vertical axis— L“ x(n))

In Fig. 3.16, we plotted the two-dimensional pseudo phase space [x, L" x] for

different a values. For comparison, the true two-dimensional phase space [x,ic] is plotted
in Fig. 3.17. From the figures, we can see that our pseudo phase space resembles the true

phase space.

41

 

u-L
U1
l

 

 

 

 

Figure 3.17 True 2-D phase space of the Lorenz system (horizontal axis—
displacement x(n), vertical axis—velocity y(n))
For comparison of the derivative method and method of delays, the two-
dimensional pseudo phase space [x(n),x(n+h)] is plotted in Fig. 3.18 for different delay h
values. The delay h is selected such that the average mutual information between x(n)

and x(n+h) is approximately equal to that between x(n) and L“x(n) for different a values

in Fig. 3.16.

42

 

 

10
. x(n+h) o

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.18 2-D pseudo phase space of the Lorenz system (horizontal axis—x(n), vertical

axis—x(n+h)

3) Two-well system
For data from the two well experiment, we apply the same ideas for evaluating the
fractional derivative order a. For comparison, we look at I vs. h for the method of delays.

The result of I vs. h is shown in Fig. 3.13. We choose h=8 and 1:06.

43

 

25 I I I I 17

...
01

d

l (average mutual information)

.0
or

 

 

 

 

h (delay)
Figure 3.19 I vs. h of the Two-well system

Also, we calculated I vs. a, which is shown in Fig. 3.14. We choose I=0.58 at

 

l (average muual information)

 

 

 

 

 

 

 

Figure 3.20 I vs. a of the Two-well system

 

 

 

 

 

 

Figure 3.21 2-D pseudo phase space of the Two-well system (horizontal axis—x(n),
vertical axis— L" x(n))

In Fig. 3.21, we plotted the two-dimensional pseudo phase space [x, L“ x] for

different a values. For comparison, the true two-dimensional phase space [x,x] is plotted
in Fig. 3.22. From the figures, we can see that our pseudo phase space resembles the true

phase space.

45

 

 

 

 

30910500

 

Figure 3.22 True 2-D phase space of the Two-well system (horizontal axis—
displacement x(n), vertical axis—velocity y(n))
For comparison between the derivative method and method of delays, two-
dimensional pseudo phase space [x(n),x(n+h)] is plotted in Fig. 3.23 for different delay h
values. The delay h is selected such that the average mutual information between x(n)

and x(n+h) is approximately equal to that between x(n) and L‘x(n) for different a values

in Fig. 3.21.

46

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.23 2-D pseudo phase space of the Two-well system (horizontal axis—x(n),

vertical axis—x(n+h)

3.5 FNN analysis—embedding dimension d5
Now that we have at hand all the vectors x(n), D“ x(n), ..., we can reconstruct
pseudo phase space [x, D“x,D2‘x,D3“x, ..., D““"'x]. To determine the embedding
dimension d E, false nearest neighbor (FNN) numbers are computed for different

reconstruction dimensions d.
In order to reduce noise amplification, the highest derivative order (d-l)a should
be “small”. Instead of choosing the values for a that minimize the AMI, we choose a=2/7

in each case so that the highest derivative order (d-l)a=(8-l)*2/7=2 when d=8. From the

47

pseudo phase portraits (Fig. 3.11, 3.16 and 3.21) and from the AMI values (Fig. 3.10,
3.15 and 3.20), we assert that this value of a produces independent reconstruction
coordinates.
l) Duffing system
For the Duffing system reconstructed with fractional derivatives, Fig. 3.24 shows
the resulting FNN vs. d where d ranges from one through eight. Note the total number of

vectors x(n) is N=57044.

 

Fllll

 

 

 

|+3eriesil56433|24819|1152I 15 | o l o I 0 | n ]

 

 

Figure 3.24 FNN vs. d (Duffing) —fractional derivatives

From Figure 3.24, we conclude the embedding dimension d£=4 or 5, which

makes sense for a 3-dimensional (t,x,:t) system expressed in terms of pseudo phase
coordinates.

From Fig. 3.10, we see that the AMI is slightly above 1, comparable to the value
for a delay of 4 in Fig. 3.9. The FNN result using MOD (h=4) is shown in Fig. 3.25. In

MOD, the embedding dimension d 5 =3 or 4.

48

 

Fllll

 

 

I:o—Series1l56905l21193|1272| o | o I o I o I o I

 

 

 

 

Figure 3.25 FNN vs. d (Duffing, MOD)
2) Lorenz system
For the Lorenz system, Fig. 3.26 shows the resulting FNN vs. d where d ranges
from one through eight. From Fig. 3.15, we see that the AMI is slightly above 1 when
a=2/7, comparable to the value for a delay of 13 in Fig. 3.14. Note the total number of

vectors x(n) is N=57044.

 

FNN

 

 

I+Series1|5$95WBIMI o I o I o I o I o I

 

 

 

 

Figure 3.26 FNN vs. d (Lorenz)—fractional derivatives

49

From Figure 3.25, we conclude the embedding dimension d 5 :4.
From Fig. 3.15, we see-that the AMI is slightly above 1, comparable to the value
for a delay of 13 in Fig. 3.14. The FNN result using MOD is shown in Fig. 3.27. In

MOD, the embedding dimension d5 =3.

 

FNN

 

 

6
+Seriesil56853l4184l o Io I o I o I o I o I

 

 

 

 

Figure 3.27 FNN vs. d (Lorenz, MOD)
3) Two-well system
For the Two-well system, Figure 3.28 shows the resulting FNN vs. d, where d

ranges from one through eight. Note the total number of vectors x(n) is N=57044.

50

 

Fllll

 

6
l—o—Series1I56461I32032I 396] e I o o I o I o I

 

 

 

Figure 3.28 FNN vs. d (T wo—well)—-fractional derivatives
From Figure 3.28 we conclude the embedding dimension dE=4 or S, which is
reasonable for a 3-dimensional autonomous system in 9I3expressed in pseudo phase
coordinates.
From Fig. 3.20, we see that the AMI is slightly above 1, comparable to the value
for a delay of 4 in Fig. 3.19. The FNN result using MOD is shown in Fig. 3.29. In MOD,

the embedding dimension d E = 4 or 5

 

Fllll

 

 

 

l+Series1I56952I31999I4103I e I 2 I 2 I o I o I

Figure 3.29 FNN vs. d (Two-well, MOD)

 

 

 

51

3.6 Extraction of unstable periodic orbits (UPOs)
Since the periodicity of UFOs directly involves time, it is possible for distortions
to occur in the MOD. This was observed in large DOF systems by Kappagantu [33].
Thus, it is of interest to examine the UFO extraction in our reconstruction.
Here we calculate recurrent points for different period number k’s. We say x(i) is
a recurrent point and i is the recurrent index, if |x(i)-x(i+kn)|<e, where x is the
reconstruction phase space vector and 8 is a pre-determined upper limit value for the
distance between points. For autonomous systems, n is 1. For non-autonomous systems, n
is the number of points (samples) per driving period.
1) Duffing system
Fig. 3.30 shows the number of recurrent points for the Duffing system with the
reconstruction vector [x, D‘x,D2"x,D3"x, ..., D“'""'x]. Our embedding dimension is
m=5 and the fractional derivative order is a=2/7. The x-axis is possible period numbers
and the y-axis is the number of recurrent points corresponding to different period
numbers. Note that the total number of points is N=57044 and there are approximately 94

samples per driving period.

52

 

number of recurrent points
8 of 8

—I
(II
I
r

 

 

 

0 5 10 15 20 25 33
period number

Figure 3.30 Recurrent points (Duffing)-—fractional derivatives

For comparison, we plotted the number of recurrent points in Fig. 3.31 in the true

phase space [x,x]. The fractional derivative reconstruction produces data with similar

densities of UFOs as in the original state space.

 

 

 

 

40 r T
35" -1
mm!- '1
.E
8.
E ..
€25?
8
2
~520- 4
E
E
315» -
10- -
5 1 1 l 1 J
0 5 10 15 20 25 33

period number

Figure 3.31 Recurrent points (Duffing, true phase space)

53

2) Lorenz system

Fig. 3.32 shows the number of recurrent points for the Lorenz system with the
reconstruction vector [x, D“x,D2"x,D3‘x, ..., D“""'x]. Our embedding dimension is
m=4 and the fractional derivative order is a=2/7. Note that the total number of points is

N=57044 and the system is autonomous (n=1).

 

 

 

 

25 r r
20» -
(0
E
8.
E15. -1
e:
3
2
‘8
310- ~
.9
E
3
5)- .4
0 r 7 l r r in, 1

 

 

0 50 100 150 203 250 303 350 4C!) 450 500
period number

Figure 3.32 Recurrent points (Lorenz)-—fractional derivatives

For comparison, we plotted the number of recurrent points in Fig. 3.33 in the true

phase space [x, x , it ]. The fractional derivative reconstruction produces data with similar

densities of UFOs as in the original state space.

54

 

6‘. B
I I
A l

d
O
U

number of recurrent points

 

 

 

 

 

 

 

5 M AL
0 1 l l l L
303 350 400 450 503

0 50 1C!) 150 200 250
period number

Figure 3.33 Recurrent points (Lorenz, true phase space)
3) Two-well system
Fig. 3.34 shows the number of recurrent points for the Two-well system with the
reconstruction vector [x, D“x,D2‘x,D3‘x, ..., D“""x]. Our embedding dimension is
m=5 and the fractional derivative order a=2l7. Note that the total number of points is

N=57044 and there are exactly 25 samples per driving period.

 

 

 

 

25 r
20h -
an
.§
8.
E15)- "
3
‘é‘
“5
.210" -I
E
C
5" -l
o l L L l l
0 20 40 60 80 1C!) 120

 

period number

Figure 3.34 Recurrent points (Two-well)—-fractional derivatives

55

For comparison, we plotted the number of recurrent points in Fig. 3.35 in the true

phase space [x,x]. The fractional derivative reconstruction produces data with similar

densities of UFOs as in the original state space.

 

 

 

 

 

 

 

 

22 j fii Y
20- -
1e- -
.‘E‘s' ‘
O
.214- 4
E f -
2
'8
h10" :-
‘3
2 5' "
61- -1
4- .1
2 l l l l L
0 20 40 50 m 100 120

period number

Figure 3.35 Recurrent points (Two-well, true phase space)

3.7 Correlation dimension
Different dimensions are calculated for the Duffing system and the two-well

system. For the correlation dimension, we applied two different algorithms.

 

 

 

 

 

The Duffing system The two-well system

True FD Delay True FD Delay

[x,y] [x,Lx} [x,.xm. ] [Ly] [XLXI [x,.x,... 1
Do 1.2359 1.1809 1.2907 1.6566 1.6387 1.6346
DI 1.3143 1.2694 1.4662 1.5025 1.3712 1.4527

 

 

 

 

 

 

 

 

56

 

D2 * 1.5495 1.5578 1.6475 1.6539 1.6887 1.6018

 

 

D2 ** 1.6436 1.6778 . 1.9155 1.6446 1.6101 1.7679

 

 

 

 

 

 

 

 

Table 2. Results of different dimensions for the Duffing system and the two-well
system.

In Table 2, D0 is the box dimension, D1 is the information dimension and D2 is
the correlation dimension. “True” means it is calculated in the true phase space [x,y];
“”FD means it is calculated in the pseudo phase space [x,Lx], where Lx is the fractional
derivative (FD); “delay” means it is calculated in the pseudo delay phase space [x, ,xm ].
D2 * was computed by box-covering for generalized dimensions, according to the
algorithm of reference [51]. D2 ** was determined by computing the correlation integral
C(r) for various ball sizes, r, then finding the slope in C(r) versus r [34, 43]. The C(r)

were computed using the algorithm of Yuan [50].

57

CHAPTER 4

RECONSTRUCTION WITH A FILTER OPERATOR

4.1 Motivation

The two Coulomb friction models in Chapter 2 are systems with stick-slip
dynamics. A “stick” is defined as an event in which a relative velocity is zero for a finite
interval of time [10]. In 1997, Feeny and Liang [10] examined chaotic single-degree-of—
freedom stick-slip friction models and used the method of delays to reconstruct the phase
space. They illustrated that the MOD reconstruction can cause pseudo trajectories to
collapse and thus the reconstructed attractor is not topologically similar to the real
attractor. The problem was cured by adding an additional observable. However, an
additional observable may not always be available. Hence we are seeking alternative

reconstruction methods to separate the pseudo signal collapse.

4.2 Basic idea
In chapter 3, we studied the use of fractional derivatives D“ for reconstructing
the phase space. In this chapter, we augment the method of delay with an operator term.
That is, if we have is discrete signal x(n), the pseudo state variables are y(n)=[x(n), Lx(n),
x(n+h), x(n+2h),...]. At first, we tried to represent Lx with an integral(Gilmore, [35]).
That is, for a measured variable x(t), we reconstruct the phase space using the vector

y(t)=[x(t), Lx(t), x(t+h), x(t+2h), ...], where Lx(t)=_Ix(t)dt and h is the delay number.

However, method of integrals is not feasible for any measured signal x(n) having

a constant item (zero-frequency component). Integration in time domain corresponds to

58

Us in Laplace transform or 1/ju) in Fourier transform. For the constant component, 5 or

ju) is 0, which makes Us or l/ja) meaningless. Therefore we propose using a first-order

filter al(a+s) in the frequency domain.

Similar as chapter 3, the procedure to calculate Lx(n) is listed below:
1) FFT of x(n) to XGw)
2) Multiply Xfiw) by al(a+jm) to get LXGw)

3) Inverse FFT of LXGw) to Lx(n)

4.3 Theoretical analysis

First, we analyze the process theoretically for the sinusoidal signal x(t)=cos( we t).

Then,

X(jw) = é-JM) -a)o) +%5(w+wo) , the Fourier transform of the filter is

 

H(ja))= a. ,hence
a+1w

LXGw)=X(jw)H(ito). Transforming to the time domain, we have

Lx(t) = —-l-—sin(ax + tan " (3))
Ja2 + (02 to
Here we look at the Bode diagram of the filter al(a+s).

a
a+jw

 

dB = 2010g I G(jw) |= 201og| I= 20*(log | a | -log I a + ja) I)

With (0 increasing, the magnitude decreases. The Bode diagram of the low pass

filter Lx(s)=a/(a+s) when a=l is shown in Fig. 4.1.

59

 

 

 

 

 

 

-5 2
n-1o ‘
E E
EH15 E
5 ..
-20 '1

'25 . ......:'r L . ......1

1o-1 1o° 1°

Frequencyaedlsec)

Figure 4.1 Bode diagram of al(a+s)

4.4 Test on sinusoidal and square time series
Before applying the derivative method into actual systems, we tested it on

sinusoidal time series. Here we take x(t)=sin(4t)-0.5sin(2t) with the initial value x(t=0)=0
at the sampling time interval 1, =21d256=0.02454s (128 samples per period n).

Fig. 4.2 shows the filtered signals Lx(n) for different a values. The solid line is
the original signal x(n) and the dotted line is the filtered signal Lx(n) with the filter

a/(s+a).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.2 Filtering sinusoidal signals

Considering that we are going to apply this filter method to the stick-slip signals,

we tested it on square waves. Here we take the initial value x(t=0)=l and the sampling

time interval 1, =21rl7683=0.00818s (384 samples per period 1:).

Fig. 4.3 shows the filtered signals Lx(n) for different a values. The solid (thin)
line is the original signal x(n) and the dotted (thick) line is the filtered signal Lx(n) with

the filter a/(s+a).

61

 

 

 

 

 

 

 

 

 

a=0. 1 a=0.7
1 a - 1 ~ a - f
o s » « 0.5 J
o \/\/ + 0
as 1 1 0.5 >
.1 L .1
o l 2 3 4 s s 7 o 1 2 3 4 5 s 1

 

 

 

 

 

 

 

v7

 

 

 

 

 

 

 

 

 

 

 

 

 

1115523621 'oiiii’éé:
Figure 4.3 Filtering square signal

From Fig. 4.3, we can see, unlike the delayed signal x(n+kh), which is still a

stick-slip signal, the filtered signal Lx(n) is not a stick-slip signal any more. Thus it may

overcome the pseudo signal collapse problem in MOD phase space reconstruction.

4.5 Look for “g ” value of a

In the filter a

 

, a is an important parameter. If a is very small, it is like an
8+0

integrator. Low frequencies are emphasized. The value of a is the break-frequency and
the bandwidth of the ln-order filter. If a is large, it is a low-pass filter, with the signal
unchanged across its bandwidth. We want the signal to be distinctive without too much

low-frequency drift.

62

In the method of delays (MOD), one effective way to determine the optimal delay
h is average mutual information (AMI) method. We apply this method to look for a good
value for the filter parameter a.

(l) Duffing system

First, we compute the average mutual information I between the vectors x(n) and
x(n+h) for different delay h’s, where x(n) is the original signal and x(n+h) is the delayed
signal. The result of I vs. 11 is shown in Fig. 4.4. In MOD, we choose h value where I
reaches a minimum value for the first time, i.e., in this case h=28 and I is about 0.22.
There are about 94 samples per driving period in this system. The optimal value of h is

close to ll3 of sample numbers per driving period.

 

 

 

 

2.2 . I I I 1
2 ----------r-----------§--v---------é .................................. —
18 P- --------- r ----------- 'r ----------- i ----------- + ----------- + ---------- -l
E‘ E E
.2 1.5 - ---------------------------------- —
a : :
E : :
e 1.4 -- ---------------------------------- -
.5 . .
g 12 b- --------- r -------- :r """"""" i ----------- v ---------- r ---------- -i
E 3 i
g, 1 .. ......... I. ........... g. ........... Q ........... s ........... A ---------- —
53 i 2
a: ' I
508 r——--- -----~-----------u---------—-: -------------------------------- --
06 - --------- r ----------- r ----------- 3 ----------- , ---------- . --------- —
04 _'. .......... u. ....... a. ........... 1 ..................... L ......... .a
0.2 I 1 ' 1 1
0 1D 20 33 40 50 50

h (delay number)

Figure 4.4 I vs. h (Duffing)
Then we calculate the mutual information I between the vectors x(n) and Lx(n)

for different filter parameters a, where x(n) is the original signal and Lx(n) is the filtered

63

signal. Different a will result in different Lx(n), thus different I. I vs. a is shown in Fig.
4.5. There is no local minimum. I gets smaller as a gets smaller. In regards to our
concerns about excessively small a, we choose I=O.22 and a=0.2, which is close to the
optimal mutual information for MOD at h=28. The ratio of the break frequency to the

excitation frequency is 0.2/l.333=0.15.

 

09 I I r T 1

 

as L -------------- :- .- 3 s -: -------------------- -

 

07

 

o
m

 

 

 

 

 

 

 

I (av otaqo mutual information)

 

 

 

a, 1 1 1 1 i

 

Figure 4.5 I vs. a (Duffing)
Fig. 4.6 below shows the phase space of [x,Lx] when we choose different values
for a. It can be compared with the true phase space of [x,y] in Chapter 3, where x is the

displacement and y is the velocity.

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.6 Pseudo phase space of the Duffing system (horizontal—original signal x(n),
vertical—filtered signal Lx(n))
(2) Friction oscillator 1
For data sets from simulation of the first friction oscillator system, we apply the
same idea to choose a good value for the parameter a. Fig. 4.7 shows the result I vs. h.
We choose h=33 and I=0.75. The optimal value of h is close to 1/3 of the number of

samples per driving period 100.

65

 

.MNN
NMbO’l

.5
0

..s
b

I (average mutual information)
- 'u or

F3
m

 

 

 

 

 

 

9
m

0 1 O 20 30 40 50 60
h (delay)

Figure 4.7 I vs. h (Oscillator 1)
Also, we calculated I vs. a, which is shown in Fig. 4.8. There is no local
minimum. I gets smaller as a gets smaller. In regards to our concerns about excessively
small a, we choose I=0.8 and a=0.5, which is close to the optimal mutual information for

MOD at h=33. The filter bandwidth frequency a is 2/5 the excitation frequency.

 

1.6

o.

      

 

.
.
--°-1
"1

..

otooooood

------J-oonco-Oc ...................................

.5
b
l

'u
r

ooouooofioooooo..uoooo

 

.5

.0
on

...-...,.......?.......'.....I.' .0...-

 

P
0)
1
I
I
l
l
‘
I
I
I
I
l

I (average mutual information)

 

 

 

 

 

 

Figure 4.8 I vs. a (Oscillator 1)

66

 

 

 

 

l
x0.5 ,» ..
..l
U
S
x 0 ..

-1 . l l l l L l l l l
50 5 10 15 2O 25 30 35 40 45 50

time

 

 

 

 

Figure 4.9 Separation of the stick-slip signal (solid thin line—the original
displacement signal x(t), dotted thick line—the signal Lx after being filtered with the
filter parameter a=0.5)
Fig. 4.9 shows the filter we applied can “separate” the stick-slip signal for this
oscillator.
Fig. 4.10 shows the phase space of [x,Lx] when we choose a=0.5. For comparison,

true phase space of the system is shown in Fig. 4.11.

67

1.4 I I I I I I

 

0.8 -

05 -

0.4 -

Lx

 

 

 

-0.2

I

1

414

 

 

415
-1.5

 

2.5

 

1.5~

 

 

 

 

-1
-1.5 -1 415 O 0.5 1 1.5 2

Figure 4.11 True phase space of the first oscillator (x—displacement, y-velocity)

68

(3) Friction oscillator 2

For data sets from simulation of the second friction oscillator system, we apply
the same idea to choose a good value for the parameter a. Fig. 4.12 shows I vs. h. Again,
for minimal mutual information, we choose h=24 at I=l.6. The optimal value of h is

approximately % of the sample number per driving period of 94 samples.

 

2.2

M

l (average mutual information)

...h
0

..s
0')

 

 

 

 

.3
h

h (delay)

Figure 4.12 I vs. h (Oscillator 2)
Also, we calculated I vs. a, which is shown in Fig. 4.13. In regards to our
concerns about excessively small a, we choose I=l.55 and a=1, which is close to the
Optimal mutual information for MOD at h=24. The ratio between a and the excitation

frequency is 0.4624.

69

 

l (average mutual information)

.5
0|

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.4
1.3 I l l l l
o 0.5 1 1.5 2 2.5 3
a
Figure 4.13 I vs. a (Oscillator 2)
1 __ I I I I r ' -
0.5 - -
D h v 1
if 415 ~
U
5
X
-1 .
f f
-1.5 >-
-2 .A -
U 5 10 15 '20 25 30
time

Figure 4.14 Separation of the stick-slip signal for the second oscillator (solid thin line—

the original velocity signal x(n) , dotted thick line—the signal Lx(n) after being filtered

with the filter parameter a=1.0)

7O

 

 

 

 

.6
01

 

 

 

 

 

 

 

l I l l

0 5 10 15 20 25 30
time

on.
-

 

Figure 4.15 Separation of the stick-slip signal for the second oscillator (solid thin
line—the original displacement signal x(t), dotted thick line—the signal Lx after being
filtered with the filter parameter a=l.0)

Fig. 4.14 shows the velocity signal and that after filtering. Fig. 4.15 shows the
displacement signal and that after filtering. Fig. 4.14 shows the filter we applied can
“separate” the stick-slip signal for this oscillator.

Fig. 4.16 shows the phase space of [x,Lx] when we choose a=l.0. For comparison,

true phase space of the system is shown in Fig. 4.17.

71

 

2 .
1.5- ~
1_ _
0.5- a

:1 0- -
-o.5e -
-1 - ..
-15» -

 

 

 

 

 

 

 

 

Figure 4.17 True phase space of the second oscillator (x—displacement, y-velocity)

4.6 FNN analysis—embedding dimension E

72

Now that we have at hand all the vectors x(n), Lx(n), x(n+h), ..., we can
reconstruct pseudo phase space [x(n), Lx(n), x(n+h), x(n+2h), ..., x(n+(d-2)h)]. To
determine the embedding dimension d E, false nearest neighbor (FNN) numbers are
computed for different delay h’s and different dimension d’s.

(l) Duffing system

For the Duffing system, as mentioned above, we choose a=0.2. Fig. 4.18 depicts
results for FNN vs. d and h, where d ranges from one through nine and h ranges from one

through sixty-one. Note the total number of vectors x(n) is N=10000.

 

h (delay) 0 ’10 d (dimension)

Figure 4.18 FNN vs. d and h (Duffing). The phase-space reconstruction includes an Lx
axis
From the figure above, the embedding dimension d 5:4 or 5, which still makes

sense for the 3-dimensional (t,x, or) system for most values of delay h.

Also, the figure shows that FNN is rather robust in indicating d 5 as h varies.

73

(2) Oscillatorl:

Stick-slip in this oscillator can cause problems in the MOD phase space
reconstruction, in which the pseudo phase trajectories collapse. This collapse produces
anomalies in the FNN calculations, such that the indicated embedding dimension varies
greatly with the delay [10].

The addition of the Lx term to the delay vector is meant to “separate” this
collapse. Successful separation should produce a more robust FNN computation.

For Oscillator 1 system, as we mentioned above, we choose a=0.5. Fig. 4.19
depicts results for FNN vs. d and h where d ranges from one through nine and h ranges
from one through sixty-one. Note the total number of vectors x(n) is N=10000.

From Fig. 4.19, the embedding dimension d 5 :3 or 4, which makes sense for the

3-dimensional (t,x, 3:) system for most values of delay h.

Also, the figure shows that FNN is rather robust in indicating d E as h varies.

4O

-5

 

h (delay) 0 '10 d (dimension)
Figure 4.19 FNN vs. d and h (Oscillator l) The phase-space reconstruction includes an Lx

axis.

74

(3) Oscillator 2:

Stick-slip causes anomalies in the FNN computation based on delay
reconstructions, such that the indicate embedding dimension varies greatly with h [10].

For Oscillator 2 system, as mentioned above, we choose a=1. Fig. 4.20 depicts
results for FNN vs. d and h where d ranges from one through nine and h ranges from one
through sixty—one. Note the total number of vectors x(n) is N=10000.

In Fig. 4.20, the indicated embedding dimension d E varies considerably with the

delay h. For small h, the embedding dimension is d 5 :3 or 4, which makes sense for the

3-dimensional (t,x,Jr) system. The lack of robustness in the FNN computation for larger
h may be an indicator that a better value of a may exist. In this oscillator, the waveform
has a nearly harmonic appearance at the sticking phase. Its higher harmonics are weaker
than those of the more square waveform of oscillator 1. Hence, the filtered waveform
may not be as distinctive as in Oscillator 1. Another Lx might be designed to improve the

FNN results (for example, another value of a, or a higher-order filter).

-6

 

h (delay) 0 ’10 - d (dimension)
Figure 4.20 FNN vs. d and h (Oscillator 2). The phase-space reconstruction includes an

Lx axis.

75

CHAPTER 5

CONCLUSIONS AND FUTURE WORK

5.] Conclusions

We have developed methods for reconstructing the phase space of nonlinear
system. We proposed a fractional derivative method, in which the pseudo phase space is
built from an observable and its successive fractional derivatives, which are easily
computed in the frequency domain. The method was tested on the Duffing and Lorenz
systems and also on data from a two-well beam experiment. Standard time-series
analyses were successfully performed.

We also proposed the use of filters as pseudo phase space coordinates. The idea of
the filter was to generate a distinctive signal in cases of stick-slip, for which intervals of
the motion (sticks) lead to indistinguishable points in delay reconstructions.

A summary of conclusions is as follows.

1. The fractional derivative method offers a derivative-based phase space
reconstruction without excessive noise amplification.

2. In FNN analysis, UPO analysis and correlation dimension calculations, the
fractional derivative method gives very good results, which shows it is probably a viable
alternative method for phase space reconstruction.

a

 

3. Using the FFT as a tool, original signal and that filtered by H(s)= along

9
0+3

with sufficient delay terms, produces a one-to-one map of the data from the phase space

76

to the embedding space, and thus phase space reconstruction using the signals will not
have drift problems.

a

 

4. The filter H(s)= can “separate” stick signals, which makes it possible to
a

+8

separate the pseudo signal collapse.

5.2 Future Work
The following future work is suggested

1. Standardize a procedure for determining a in both fractional derivative and filter
methods.

2. Look deeply at the effect of a and drift. Will drift occur as a-90 in FFT based
methods?

3. Qualify distortion and evaluate reconstruction methods.

77

BIBLIOGRAPHY
[1] Whitney, H., 1936, Differentiable Manifolds, Ann. Math. 37 pp645

[2] Packard, N.H., Crutchfield, J.P., Farmer J.D., and Shaw, RS, 1980, Geometry from a
time series, Phys. Rev. Lett. 45, pp712-716

[3] Takens, F., 1980, Detecting strange attractors in turbulence, in: Dynamical Systems
and Turbulence, eds. D. Rand and L. S. Young, Lecture Notes in Mathematics, Vol. 898
Springer, Heidelberg, pp366-38l.

[4] Broomhead, D.S. and King, G.P., 1986, Extracting qualitative dynamics from
experimental data, Physica 20D, pp 217-36

[5] Gibson, J.P., Farmer, J.D., Casdagli, M. and Eubank, S., 1992, An analytic approach
to practical state space reconstruction, Physica 57D, ppl-30

[6] Kennel, M.B., Brown, R. and Abarbanel, H.D.I., 1991, Determining embedding
dimension for phase-space reconstruction using a geometrical construction, Physical
Review A, Vol 45, No. 6, pp3403-34ll

[7] Cao, L, 1997, Practical method for deterrning the minimum embedding dimension of
a scalar time series, Physica 110D, pp43-52

[8] Sauer, T., Yorke J.A., and Casdagli, M., 1991, Embedology, Journal of Statistical
Physics, Vol 65, pp579-6l6

[9] Read, P.L., 1993, Phase portrait reconstruction using multivariate singular systems
analysis, Physica D. 69, pp353-365

[10] Feeny, BF. and Liang, J.W., 1997, Phase-Space Reconstruction and Stick-Slip,
Nonlinear Dynamics l3, pp39-57

[l 1] Higgins, R.J., 1990, Digital Signal Processing in VLSI, Prentice Hall, ppl39-l43
[12] Gouesbet, G. and Maquet, J., 1992, Construction of phenomenological models from
numerical scalar time series, Physica 58D, pp202-15

[13] Hignett, P., 1985, Characteristics of amplitude vacillation in a differentially heated
rotating fluid annulus, Geophys. Astrophys. Fluid Dyn. 31, pp247-281.

78

[14] Jacobs, O.L.R., 1974, Introduction to Control Theory
[15] Koffman, E.B., Friedman, P.L., 1997, Koffman-Friedman FORTRAN

[16] Wiercigroch, M., Kraker, B.D., 2000, Applied Nonlinear Dynamics and Chaos of
Mechanical Systems with Discontinuities

[l7] Ibrahim, R.A., Soom, A., 1992, Friction-induced Vibration, Chatter, Squeal, and
Chaos

[l8] Higgins, Richard J., 1990, Digital Signal Processing in VLSI

[19] Elliott, DR, 1987, Handbook of Digital Signal Processing Engineering Application
[20] Cadzow, J .A., Landingham, H.F.V., 1985, Signals, Systems and Transforms

[21] Carlson, GE, 1992, Signal and Linear System Analysis

[22] Peitgen J urgens Saupe, 1992, Chaos and Fractals

[23] Strogatz, Steven H., 1994, Nonlinear Dynamics and Chaos

[24] Potapov, A., 1997, Distortions of reconstruction for chaotic attractors, Physica D
101, pp207-226

[25] Moon, EC, 1987, Chaotic Vibrations, Wiley-Interscience, New York

.[26] Guckenheimer, J. and Holmes, P.J., 1983, Nonlinear Oscillations, Dynamic Systems
and Bifurcations, Springer-Verlag, New York

[27] Moon, RC, 1992, Chaotic and Fractal Dynamics—An Introduction for Applied
Scientists and Engineers, John Wiley and Sons, New York

[28] Anderson, J .R. and Ferri, A.A, 1990, Behavior of a single-degree-of-freedom system
with a generalized friction law, Journal of Sound and Vibration, 140(2), pp287-304

[29] Feeny, B. F., Yuan, C. M. and Cusumano, J. P., 1999, Parametric Identification of a
Two-well Oscillator, DETC 99, Las Vegas, on CD-ROM

[30] Yuan, C. M. and Feeny, B. F., 1998, Parametric Identification of Chaotic Systems,
Journal of Vibration and Control 4 (4), pp405-426

[31] Malinetskii, G.G., Potapov, A.B., 1993, Limitations of Delay Reconstruction for

Chaotic Dynamical Systems, Physical Review E: Statistical Physics, Plasmas, Fluids and
Related Interdisciplinary Topics, v48, n2, pp904

79

[32] Fraser, A.M., 1989, Reconstructing Attractors from Scalar Time Series: a
Comparison of Singular System Analysis and Redundancy Criteria, Physica D 34, pp39l-
404

[33] Kappagantu, R., 1997, An Optimal Modal Reduction for Frictionally Excited
Systems, PhD thesis, Michigan State University

[34] Abarbanel, H.D.I., 1993, The Analysis of Observed Chaotic Data in Physical
Systems, Review of Modern Physics, Vol.65, No. 4. pp1331-1392

[35] Gilmore, Robert, 1998, Topological analysis of chaotic dynamical systems, Review
of Modern Physics, Vol. 70, No. 4, ppl455-1526

[36] Zamel, Z.Al and Feeny, BR, 2001, Improved Estimations of Unstable Periodic
Orbits Extracted From Chaotic Sets, Proceedings of DETC’Ol, 2001 ASME Design
Engineering Technical Conferences, September 9-12, 2001, Pittsburgh, Pennsylvania,
USA

[37] Popp, K and Stelter, P., 1990, Stick-slip vibrations and chaos, Philosophical
Transactions of the Royal Society of London, A, 332, pp89-105

[38] Tseng, C.C., Pei, SC, and Hsin, SC, 2000, Computation of fractional derivatives
using Fourier-transform and digital FIR differentiator, Signal Processing, 80(1), pp151-
159

[39] Auerbach, D., Cvitanovic, P., Eckmann, J .P., Gunaratne, G., and Procaccia, I., 1987,
Exploring Chaotic Motion through Periodic Orbits, Physical Review Letters 58, pp2387-
2389

[40] Lathrop, DP. and Kostelich, E.J., 1989, Characterization of an experimental strange
attractor by periodic orbits, Physical Review A 40, pp4028-4031

[41] Grebogi, C., Ott, E., and Yorke, J., 1988, Unstable periodic orbits and the
dimensions of multifractal chaotic attractors, Physical Review A 37(5), pp17 l 1-1724

[42] Van de Water, W., Hoppenbrouwers, M., and Christiansen, F., 1991, Unstable
Periodic Orbits in the Parametrically Excited Pendulum, Physical Review A 44(10),
pp6388-6398

[43] Ott, E., 1993, Chaos in Dynamical Systems, Cambridge University Press, NY.

[44] Tufillaro, N.B., Abbott, T., and Reilly, J., 1992, An Experimental Approach to
Nonlinear Dynamics and Chaos, Addison-Wesley

[45] Tufillaro, N ., Solari, H., and Gilmore, R., 1989, Physical Review Letters 64, pp2350

80

[46] Tufillaro, N., Solari, H., and Gilmore, R., 1990, Relative Rotation Rates—
Fingerprints for Strange Attractors, Physical Reviews A 41, pp5717

[47] Hammel, S., and Heagy, J., 1992, Chaotic system identification using linked
periodic orbits, presented at the SIAM Conference on Applications of Dynamics
Systems, Snowbird, Utah, October 15-19

[48] Kesaraju, RV. and Noah, S.T., 1994, Characterization and detection of parameter
variations of nonlinear mechanical systems, Nonlinear Dynamics 6, pp433-457

[49] Van de Wouw, N ., Verbek, G., and Van Campen, DH, 1995, Nonlinear parametric
identification using chaotic data, Journal of Vibration and Control 1, pp291-305

[50] Yuan, C.M., 1995, A Method of Parametric Identification for Chaotic Systems, PhD
Thesis, Michigan State University, East Lansing

[51] Feeny, BR, 2000, Fast Multifractal Analysis by Recursive Box Covering,
International Journal of Bifurcation and Chaos, 10(9) pp2277-2287

81