Multi-variate statistical data analysis techniques offer a powerful tool for analyzing power system response from measured data. In this thesis, a statistically-based, data-driven framework that integrates the use of complex wavelet empirical orthogonal function (EOF) analysis and the method of snapshots is proposed to identify and extract dynamically independent spatio-temporal patterns from time synchronized data. The procedure allows identification of the dominant spatial and temporal patterns in a complex data set and is particularly well suited for the study of standing and propagating features that can be associated with electromechanical oscillations in power systems. It is shown that, in addition to providing spatial and temporal information, the method improves the ability of conventional correlation analysis to capture temporal events and gives a quantitative result for both the amplitude and phase of motions, which are essential in the interpretation and characterization of transient processes in power systems. The efficiency and accuracy of the developed procedures for capturing the temporal evolution of the modal content of data from time synchronized phasor measurements of a real event in México is assessed. Results show that the proposed method can provide accurate estimation of nonstationary effects, modal frequency, time-varying mode shapes, and time instants of intermittent or irregular transient behavior associated with abrupt changes in system topology or operating conditions.
Centro de Investigación y de Estudios Avanzados del
Instituto Politécnico Nacional Unidad Guadalajara
Extracción de Patrones Dinámicos de
Mediciones de Área Amplia mediante
Funciones Ortogonales Empíricas
Tesis que presenta:
Pedro Esquivel Prado
para obtener el grado de:
Doctor en Ciencias
en la especialidad de:
Ingeniería Eléctrica
Director de Tesis:
Dr. Arturo Román Messina
CINVESTAV del IPN Unidad Guadalajara, Guadalajara, Jalisco, Agosto de 2010.
Centro de Investigación y de Estudios Avanzados del
Instituto Politécnico Nacional Unidad Guadalajara
Extraction of Dynamic Patterns from
Wide-Area Measurements using
Empirical Orthogonal Functions
A thesis presented by:
Pedro Esquivel Prado
to obtain the degree of:
Doctor in Science
in the subject of:
Electrical Engineering
Thesis Advisors:
Dr. Arturo Román Messina
CINVESTAV del IPN Unidad Guadalajara, Guadalajara, Jalisco, August 2010.
Extracción de Patrones Dinámicos de
Mediciones de Área Amplia mediante
Funciones Ortogonales Empíricas
Tesis de doctorado en Ciencias
Ingeniería Eléctrica
Por:
Pedro Esquivel Prado
Maestro en Ciencias en Ingeniería Eléctrica
CINVESTAV del IPN Unidad Guadalajara 2005-2007
Becario del CONACYT, expediente no. 172551
Director de Tesis:
Dr. Arturo Román Messina
CINVESTAV del IPN Unidad Guadalajara, Agosto de 2010.
Extraction of Dynamic Patterns from
Wide-Area Measurements using
Empirical Orthogonal Functions
Doctor of Science Thesis
In Electrical Engineering
By:
Pedro Esquivel Prado
Master in Science in Electrical Engineering
CINVESTAV del IPN Unidad Guadalajara 2005-2007
Scholarship granted by CONACYT, No. 172551
Thesis Advisors:
Dr. Arturo Román Messina
CINVESTAV del IPN Unidad Guadalajara, August 2010.
v
Acknowledgments
My parents: Ignacio Esquivel Alcaraz and Lilia Prado Robles for giving me all
their love, comprehension and for motivating me to study.
My brothers: María Inés, José Alejandro, Ezequiel, Olivia, Mauricio, María
Reyna Leticia and to Patricia for giving me all their support.
I would like to express my sincere thanks to my thesis supervisor, Ph. D. Arturo
Román Messina for providing me with the necessary support, advice, facilities
and the enthusiasm required to successfully complete this thesis; the search for
relevant literature is very much appreciated.
I wish to thank my professors and friends at CINVESTAV Guadalajara, as well
as to my friends at my birth place (El Papalote, Nayarit, México).
Finally, acknowledgements are also due to CONACYT for its support.
vi
Resumen
Las técnicas de análisis de datos estadísticos multi-variables ofrecen una
poderosa herramienta para analizar la respuesta en sistemas de potencia a
partir de datos medidos. En esta tesis, un marco estadístico basado en los datos
que integra el uso del análisis de las funciones ortogonales empíricas complejas
y el método de snapshots es propuesto para identificar y extraer patrones
espacio-temporal independientes dinámicamente a partir de mediciones
fasoriales sincronizadas. El procedimiento permite la identificación de los
patrones espaciales y temporales dominantes en un conjunto de datos
complejos y es particularmente conveniente para el estudio de características
estacionarias y de propagación que pueden ser asociadas con oscilaciones
electromecánicas en sistemas de potencia. Se muestra que, además de
proporcionar información espacial y temporal, el método mejora la habilidad
del análisis de la correlación convencional para capturar eventos temporales y
proporciona un resultado cuantitativo tanto para la amplitud como para la fase
del movimiento; las cuales son esenciales en la interpretación y caracterización
de procesos transitorios en sistemas de potencia. La eficiencia y exactitud del
procedimiento desarrollado para capturar la evolución temporal del contenido
modal de los datos a partir de medidas fasoriales sincronizadas en tiempo de
un evento real en México es evaluado. Los resultados muestran que el método
propuesto puede proporcionar estimación exacta de efectos no estacionarios,
frecuencia modal, forma del modo variante en el tiempo, y los instantes de
tiempo del comportamiento transitorio irregular o intermitente asociado con
cambios abruptos en la topología del sistema o condiciones de operación.
vii
Abstract
Multi-variate statistical data analysis techniques offer a powerful tool for
analyzing power system response from measured data. In this thesis, a
statistically-based, data-driven framework that integrates the use of complex
wavelet empirical orthogonal function (EOF) analysis and the method of
snapshots is proposed to identify and extract dynamically independent spatio-
temporal patterns from time synchronized data. The procedure allows
identification of the dominant spatial and temporal patterns in a complex data
set and is particularly well suited for the study of standing and propagating
features that can be associated with electromechanical oscillations in power
systems. It is shown that, in addition to providing spatial and temporal
information, the method improves the ability of conventional correlation
analysis to capture temporal events and gives a quantitative result for both the
amplitude and phase of motions, which are essential in the interpretation and
characterization of transient processes in power systems. The efficiency and
accuracy of the developed procedures for capturing the temporal evolution of
the modal content of data from time synchronized phasor measurements of a
real event in México is assessed. Results show that the proposed method can
provide accurate estimation of nonstationary effects, modal frequency, time-
varying mode shapes, and time instants of intermittent or irregular transient
behavior associated with abrupt changes in system topology or operating
conditions.
viii
Index
Chapter 1 ... 1
Introduction
1.1 Background and Motivation ... 2
1.2 Statement of the Problem ... 3
1.3 A Brief Review of Previous Work ... 5
1.3.1 Model-Based Approaches ... 5
1.3.2 Measurement-Based Techniques ... 6
1.3.3 Spatio-Temporal Statistical Models ... 7
1.4 Thesis Objectives ... 8
1.5 Contributions ... 9
1.6 Summary and Outline of the Dissertation ...10
1.7 Publications...11
1.7.1 Refereed conferences ...11
1.7.2 Refereed journals ...11
1.7.3 Book chapters ...12
1.8 References ...12
Chapter 2 ... 15
Empirical Orthogonal Function Analysis
2.1 Notation ...16
2.2 Theoretical Development ...16
2.3 Discrete Domain Representation ...20
2.4 The Method of Snapshots...22
ix
2.5 Energy Relationships ...24
2.6 Interpretation of EOFs Using Singular Value Decomposition ...25
2.6.1 Singular Value Decomposition ...26
2.6.2 Relation with the Eigenvalue Decomposition ...27
2.6.3 Geometric Interpretation Using SVD ...29
2.7 Relation Between POMs and Normal Modes ...29
2.7.1 The Notion of Normal Modes ...29
2.7.2 Physical Interpretation of POMs in Terms of the Eigenvalues ...31
2.7.3 Physical Interpretation of POMs Using Singular Value Decomposition ...33
2.8 References ...36
Chapter 3 ... 38
Complex Empirical Orthogonal Function
3.1 Complex EOF Analysis ...39
3.2 Theoretical fundamentals of Empirical Orthogonal Functions...41
3.3 Complex-Vector EOF Analysis ...44
3.4 Analysis of Propagating Features ...46
3.4.1 Spatial Amplitude Function,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
S
S
...47
3.4.2 Spatial Phase Functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
II
I
...47
3.4.3 Temporal Amplitude Functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
R
R
...47
3.4.4 Temporal Phase Functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
...48
3.4.5 Spatio-Temporal Analysis ...49
3.4.6 Motivating Examples: Modeling of Propagating Waves Using the Covariance
Matrix ...52
3.4.7 Coherency Identification from Measurements Data ...59
3.4.8 Mode Shape Identification ...60
3.5 Numerical Computation of POMs: The Algorithm ...60
3.5.1 Numerical Example ...61
3.6 Concluding Remarks ...68
3.7 References ...68
x
Chapter 4 ... 70
Near-Real-Time of Measured Data Using Complex EOF
Analysis
4.1 Motivation ...71
4.2 Analytical Aspects of Time-Dependent EOFs ...71
4.3 Near-Real-Time EOF Analysis ...72
4.3.1 Time-Dependent Covariance Matrix ...74
4.3.2 Spectral Decomposition of Time-Dependent Covariance Matrix ...77
4.3.3 Time-Dependent Basis Functions ...78
4.3.4 Feature Extraction ...79
4.3.5 Wave Components ...79
4.4 Characterization of Spatio-Temporal Behavior ...82
4.4.1 Temporal Characterization ...82
4.4.2 Spatial Characterization...82
4.5 Proper Orthogonal Decomposition of Nonstationary Random Processes ...83
4.5.1 Statistics of Random Processes ...84
4.5.2 Stationarity of Random Process ...85
4.5.3 Moments of Random Variables ...85
4.5.4 Expectation Value and Standard Deviation ...86
4.5.5 Recursive Algorithm ...88
4.6 Multiscale EOF Analysis ...90
4.6.1 Wavelet-Based EOF Analysis ...90
4.6.2 Empirical Mode Decomposition-Based EOF Analysis ...91
4.7 Conclusions ...93
4.8 References ...94
Chapter 5 ... 97
Application
5.1 Complex Empirical Orthogonal Function Analysis of Power-System Oscillatory
Dynamics ...98
xi
5.1.1 Application to Time Synchronized Measured Data ...98
5.1.2 Construction of POD Modes via the Method of Snapshots ... 101
5.2 Spatiotemporal Analysis of Measured Data ... 102
5.2.1 Temporal Properties ... 105
5.2.2 Frequency Determination from Instantaneous Phases ... 107
5.2.3 Mode Shape Estimation ... 108
5.2.4 Energy Distribution ... 109
5.3 Extraction in near-real-time of propagating Features of Power-System Oscillatory
Dynamics ... 110
5.4 Discussion ... 116
5.5 References ... 117
Chapter 6 ... 119
Conclusions and Recommendations for Future Work
6.1 General Conclusions ... 119
6.2 Recommendations for Future Work ... 121
Appendix A ... 123
EOFs in the Presence of Measurement Error
A.1 References...125
Appendix B ... 126
Optimal Linear Predictor
B.1 References...129
Appendix C ... 130
Correlation Function Using Linear Instantaneous Estimates
xii
C.1 References...131
xiii
Index of Figures and Tables
Figures
Chapter 3
Fig. 3. 1. Conceptual diagram illustrating the phenomenon of temporal and
spatial energy in a field.
... 43
Fig. 3. 2. Conceptual view of the proposed algorithm.
... 61
Fig. 3. 3. Space-time evolution of propagating wave on proposed field.
... 62
Fig. 3. 4. Standing wave estimated from original field.
... 63
Fig. 3. 5. Frequency spectrum.
... 64
Fig. 3. 6. Amplitude of the standing wave
... 64
Fig. 3. 7. Travelling wave estimated from field original.
... 65
Fig. 3. 8. Frequency spectrum.
... 66
Fig. 3. 9. Wave constant.
... 67
Chapter 4
Fig. 4. 1. Visualization of PMU data in terms of spatial and temporal formation
...73
Fig. 4. 2. Conceptual diagram of the sliding window-based approach...75
Fig. 4. 3. Proposed algorithm to the complex EOF identification...89
xiv
Chapter 5
Fig. 5. 1. Schematic of the MIS system showing the location of the observed
oscillations. Measurement locations are indicated by shaded circles.
... 99
Fig. 5. 2. Time traces of recorded bus frequency swings recorded January 1, 2004
and detail of the oscillation build up.
... 100
Fig. 5. 3. Time-Varying means and fluctuating components of the recorded bus
frequency signals. The time series are smoothed using wavelet shrinkage.
... 104
Fig. 5. 4. Comparison of the Fourier transform spectrum of the original signal
and the spectra constructed with the instantaneous mean removed.
... 104
Fig. 5. 5. Energy captured as a function of the number of modes. The percentage
of energy located in the jth mode is measured by
¦
n
i
i
i
E
1
100
.
... 104
Fig. 5. 6. Reconstruction of the original data using the three leading POMs. Solid
line represent the original time series and dotted lines represent the composite
oscillation obtained by adding the temporal modes.
... 105
Fig. 5. 7. Temporal patterns of variability associated with the dominant mode.
... 106
Fig. 5. 8. Instantaneous frequency of POM 1. Time interval 06:27:42-06:29:39.
107
Fig. 5. 9. Mode shape of POM 1 for various time intervals of interest.
... 109
Fig. 5. 10. Energy distribution.
... 110
Fig. 5. 11. Time trace of bus frequency swings recorded on January 1, 2004 and
detail of the oscillation build up.
... 111
Fig. 5. 12. Reconstruction of the original data with the traveling wave mode
identified using the conventional (off-line) and recursive (on-line,
2
W
) EOF
analysis.
... 113
xv
Fig. 5. 13. (Upper) average energy of travelling wave mode using the
conventional (off-line) EOF analysis, and (lower) instantaneous energy of
travelling wave mode using the recursive (on-line,
2
W
) EOF analysis.
... 114
Fig. 5. 14. Instantaneous frequency of dominant mode.
... 114
Fig. 5. 15. (Upper) mode shape of travelling wave mode using the conventional
(off-line) EOF analysis, and (lower) time-varying mode shape of travelling wave
mode using the recursive (on-line,
2
W
) EOF analysis.
... 115
Fig. 5. 16. Spatial patterns of the leading mode showing temporal variability.
116
Tables
Chapter 5
Table 5. 1. RMS error in EOF approximation...113
1
Chapter 1
Introduction
Many physical processes in power systems involve variability over both space
and time. This variability can be quite complicated, including variability caused
by protective control actions, automatic control, and load characteristics.
This introductory chapter presents an outline of the research work in this
thesis, and defines concepts linked to power system dynamic behavior.
Furthermore, the motivation for the study of spatial and temporal variability in
power system regimes is presented. The discussion will apply equally well to
numerical or measured data.
The objectives of the work are stated and the main contributions are
summarized. An outline of the structure of the thesis is also given.
2
1.1 Background and Motivation
Many dynamic processes arising in power system behavior involve phenomena
that vary in time and space in complicated ways [1-4]. Examples of such
processes include inter-area oscillations between generators or interconnected
electrical systems and frequency stability [5-8].
Inter-area oscillations of electromechanical nature have been frequently
observed in large interconnected power systems involving areas or groups of
generators swinging against each other. Wide-area monitoring systems offer the
possibility to analyze and characterize complicated oscillations involving the
exchange of kinetic energy between machines or electrical areas in large,
interconnected power systems [3-13]. As the number of measured signals
increase, accurate characterization of relevant modal behavior becomes
difficult.
Complex data sets obtained from synchronized measurements at various
system locations tend to produce a large amount of information. Moreover,
different faults may result in different dynamic patterns thus requiring the use
of correlation statistically-based techniques to examine dynamic trends and
phase relationships.
Characterization of spatial and temporal changes in the dynamic pattern of
system oscillations is a problem of great theoretical and practical importance
that has been barely studied.
In the study of complicated inter-area oscillations it is of importance to
analyze not only the temporal behavior of the underlying phenomena, but also
the physical (spatial) distribution of the phenomena across the system. This
phenomenon has some similarities with the analysis of multivariate data such
as that obtained from measurements of wave phenomena, climate variations or
turbulence in complex flow. Further, propagating features should be identified
in order to assess control requirements and the extent and degree of
nonlinearity.
3
Multivariate data analysis techniques provide a useful means to cope with
increasing complexity and information from wide-area monitoring schemes, to
allow statistical data to predict future performance, and to deal with existing
and new uncertainties [1,4,8-13].
Modelling and controlling such phenomena is a complex spatio-temporal
dynamical problem with important practical applications for which few
analytical models have been proposed. The computation time required to solve
large interconnected systems is large and many become prohibitive for practical
systems. To reduce the complexity of the problem, several simplifications are
commonly used that may result in a poor characterization of global system
behavior.
The study of analytical techniques to extract the most significant information
from simultaneous measurements of system oscillations constitutes the main
focus of this work.
1.2 Statement of the Problem
As noted earlier, phenomena observed in power system oscillatory dynamics
are diverse and complex. Power system ambient data measurements are known
to exhibit noisy, nonstationary fluctuations resulting from random changes in
load and other random phenomena. Forced oscillations resulting from major
events, on the other hand, may contain motions with a broad range of temporal
scales and can be highly nonlinear and time-varying.
Low-frequency system oscillations have been experienced in many power
systems posing a major obstacle to normal operation and severely limiting
transmission capabilities.
System studies on the nature of these oscillations are based on modal
analysis of analytical models of system behavior or in sequential analysis of
time series obtained from simulation of the nonlinear model or measured data
[1]. The problem of selecting the most significant modes is of considerable
interest and has been studied intensively by several researchers.
4
Approaches such as Prony analysis, eigenrealization algorithms, and auto-
regressive processes, among others, have been successfully used to characterize
modal behavior based on a single data set. Linear analysis of power system
models can also be used to extract modal information as well as to design
system controllers [6,14-20]. However, such an approach is computationally
intensive since a large number of studies is needed to estimate various
contingency scenarios and cannot be directly applied to multivariate data.
Due to the size and dimension of realistic power system models, the volume
of information that is required to characterize the inter-area mode phenomenon
from measured data is very high. This includes the determination of dominant
oscillatory modes, the analysis of phase relationships between key system
signals and the extraction of dynamic trends and nonlinearities.
A second challenge of phasor measurements unit (PMU) applications is how
to utilize the large volumes of data. Ways are needed to conveniently
summarize and visualize the large amount of data generated by dynamic
system recorders, and extract the information relevant to the phenomenon of
interest.
Further, data gaps, and degraded spatial resolutions could result in biased
or poor estimates of power system electromechanical modes.
When simultaneously measured responses throughout a system are
available, dominant linear modes should be extracted using correlation
techniques rather that individual analysis of the system response. This provides
a global picture of system behavior and enables statistical characterization of
the observed phenomena [21-22].
Large stressed interconnected power systems exhibit nonlinear dynamic
performance when are subjected to large disturbances. The extraction of
significant dynamics from multiple measurements is a challenging problem.
Factors such as noise, trends and nonlinearities may affect modal extraction or
preclude the identification of phase relationships among output signals.
5
Due to its large number of degrees-of-freedom, and the sparse geographical
distribution, power systems exhibit highly complex phenomena including
modal interactions, temporarily chaotic vibrations, and intermodulation [6].
An understanding of these phenomena is necessary if the dynamics and
ultimately the propagating nature of these oscillations are to be predicted.
To determine the mechanisms governing these physical variations, it is
essential to characterize the large-scale interactions between the system modes.
Today wide-area analyses of inter-system oscillations require a systematic
approach that extracts significant information. Multivariate data analysis
techniques offer a powerful tool to analyze sets of data obtained from
simulations or measurements that supplements information on conventional
analysis techniques.
An appealing feature of these approaches is their ability to extract more
energy fluctuations than conventional linear decomposition techniques. Also
the modal representations of the observed data are physically interpretable and
provide meaningful insights into the observed dynamics.
We provide here a brief summary of literature describing the application of
statistical techniques to nonlinear oscillatory problems in power systems.
1.3 A Brief Review of Previous Work
1.3.1 Model-Based Approaches
Model-based approaches have been successfully used to extract modal
information as well as to synthesize transfer functions and to design system
controllers.
These approaches are based on analytical approximations obtained from the
series expansion of the nonlinear power system model around a given
operating condition. This theory has been employed to determine both, the
interacting modes of oscillation and the degree of nonlinearity.
6
Recently, nonlinear interaction between the fundamental modes of
oscillation has been studied using normal form theory and bilinear analysis
techniques [8,10,23].
Because these techniques are based on analytical models, reliable and
accurate methods for modal analysis of large system models are critically
needed for determining dynamic information. In addition, models needed to
capture nonlinear phenomena require large, complex analytical formulations
that makes unsuitable for the study of practical systems or characterization of
random behavior.
With the increasing utilization of measurement devices throughout the
system, measurement-based techniques are required that supplement
information to analytical models.
1.3.2 Measurement-Based Techniques
Measured data provide a useful complement to model-based methods in
detecting temporal and frequency variability directly from time-synchronized
measurements.
Signal processing methods for estimating power system electromechanical
modes can, in turn, be broadly classified into two categories: ringdown
algorithms and mode-meter methods. The second approach addresses the
problem of estimating system behavior from power system ambient noise
measurements.
A large collection of methods have been investigated for extracting modal
information from measured system response. These methods include such
diverse strategies as Ritz vector reduction, rational approximation, and proper
orthogonal decomposition. Despite the considerable collection of research for
extracting modal information, there remain significant barriers to the extraction
of relevant information, specifically for the analysis of real-time measurements.
There are three significant drawbacks to such an approach. The first is that
existing algorithms are restricted to essentially linear signals. The second is the
large amount of data to be processed for accurate mode extraction.
7
In addition, event-based oscillations cannot be evaluated in isolation from
other events. It is clear from the above review that statistical analysis tools are
needed with the ability to treat large data sets, and enable extraction of
dominant patterns.
1.3.3 Spatio-Temporal Statistical Models
Statistical techniques have been widely used in many engineering and science
applications. These include Geophysical research, climate variations and
vibratory problems [24-31]. Among these techniques, the proper orthogonal
decomposition (POD) method (also called empirical orthogonal function
analysis and Karhunen-Loève technique) has been shown to be capable of
representing complicated phenomena with a handful of degrees of freedom
[12].
Recently, several methods have been proposed for spatio-temporal analysis
of measured data. Among the several methods available of analysis, empirical
orthogonal function (EOF) analysis is a particularly useful tool in studying large
quantities of multi-variate data [32]. The decomposition generates a physically
motivated basis for the data. More recently, these techniques have gained wide
popularity in applications related to data analysis and reduced-order modelling
of various physical processes or models.
Applications, for example, to unsteady fluid flow, turbulence, aerospace,
optimal control, structural dynamics, microstructural design, solution of
stochastic partial differential equations, heat transfer and non-destructive
testing and system identification have been reported [1,14-15,22,24-32].
In recent years there have been several studies in which empirical
orthogonal function analyses were used to highlight potential physical
mechanisms associated with system variability. One of the unresolved issues is
the estimation of propagating and standing features that can be associated with
observed or measured data and their application to real-time data.
A limited number of studies on statistical analysis of mode propagation in
nonlinear power system models exist. In this context a number of recent studies
have sought to identify such interactions and their associated time scales, by
8
analyzing both observational records and the results of power system
simulations [14-15].
In this thesis, we extend these formulations, and use POD-based methods to
extract modal information from measured data.
The conceptual framework developed provides rigorous bases for the
analysis, detection, and simplification of complex oscillations, and enables the
simultaneous study of synchronized measurements.
The advantage of this technique is that is easy to implement, efficient for
short-time records, and is based on the covariance structure of the data. The
approach is general and can be used to study both, measured data and
simulated data.
Primary applications of statistical methods for studying oscillatory dynamic
systems include modal extraction and management of disturbance recording.
1.4 Thesis Objectives
The primary objective of this research is to develop an analytical, multiscale
data analysis framework for identifying the dominant spatial, temporal and
frequency variability regimes among critical modes from real-time measured
data. Following the above problem, the specific objectives of this research are:
1. To develop a systematic, data-driven low-dimensional representation of
complex inter-system oscillations based on proper orthogonal
decomposition for the analysis, characterization and simplification of
oscillatory processes. An approach to the joint estimation of temporal
frequency and phase characteristics is to this end developed.
2. To extend existing approaches to characterize the spatio-temporal dynamics
of inter-area oscillations using near-real-time information.
3. To analyze spatio-temporal complexity and the nature of system
oscillations. Unlike other spectral analysis techniques which attempt to fit a
spectrum to a presupposed model function, the proposed technique makes
9
no assumptions concerning spectral characteristics and can treat both spatial
and temporal information.
Other identified objectives are:
x To evaluate the performance of the method under noisy, nonlinear
conditions using real-time information.
x To assess the ability of EOF-based models for efficient dynamic modelling of
large data sets.
x To address numerical issues associated with the analysis of large, complex
informating resulting from transient oscillations.
1.5 Contributions
The main original contributions of this research are:
1. The development of a mathematical and computational framework for data-
driven quantification of linear and nonlinear phenomena in physical
processes.
2. The derivation of an analytical framework to extract modal information.
Specifically the development of efficient algorithms to extract mode shapes,
and phase relationships between selected variables or signals. An approach
to space-time prediction that achieves dimension reduction and that uses a
statistical model is also presented.
3. The generalization of existing approaches to account for complex
formulations of the proper orthogonal decomposition and its numerical
implementation.
4. The evaluation of alternative formulations based on singular value
decomposition analysis to characterize standing and propagating features in
complex phenomena. An extension of this method handles the temporal
variation commonly present in propagating waves.
10
5. The derivation of analytical criteria to characterize physical behavior in
terms of physically meaningful modes in mechanical systems.
1.6 Summary and Outline of the Dissertation
This dissertation is divided into six chapters and three appendices. A brief
summary of each chapter is described below.
Chapter 1 contains the background and the overview of this thesis. The
organization of the remainder of this dissertation is as follows.
In Chapter 2, an overview of the basic concepts and approximations leading
to the empirical orthogonal function and singular value decomposition
methods is given, together with an introduction to numerical aspects for proper
orthogonal decomposition.
Chapter 3 proposes a complex POD formulation that combines proper
orthogonal decomposition with Hilbert analysis. The methodology for POD
analysis is developed in detail, and its relation to other techniques is described.
In Chapter 4 a statistically-based data-driven framework that integrates the
use of wavelet-based EOF analysis and a sliding-window based method is
proposed to identify and extract, in near-real-time, dynamically independent
spatio-temporal patterns from time-synchronized measurements.
Different
orthogonal systems are investigated including orthogonal wavelets, Hilbert
analysis, and empirical orthogonal functions. The developed methods combine
the ability of empirical orthogonal function analysis to decorrelate complex
vector fields with that of wavelet analysis to extract deterministic features.
Wavelet-EOF analysis is particularly appropriate for the modelling of transient
data containing contributions from events whose behavior changes over time
and frequency.
Chapter 5 discusses the practical application of the developed
methodologies to the study of a practical power system. The salient features of
the proposed methodology are demonstrated by applying it on measured data
sets.
11
Finally, some concluding remarks, suggestions for future work and
significant contributions, are presented in Chapter 6.
Certain details of the implementation of the methodology are relegated to an
appendix. The appendices cover the analysis of EOFs in the presence of
measurement error, optimal linear prediction, the analysis of correlation
functions using linear estimates.
1.7 Publications
The following publications have resulted from this work:
1.7.1 Refereed conferences
1. P. Esquivel, and A. R. Messina, "Complex Empirical Orthogonal Function
Analysis of Wide-Area System Dynamics," IEEE Power and Energy Society,
General Meeting (PES)
, Pittsburgh, Pennsylvania USA, 2008.
2. P. Esquivel P., "Wide-Area Wave Motion Analysis Using Complex
Empirical Orthogonal Functions," International Conference on Electrical
Engineering, Computing Science and Automatic Control (ICEEE)
, Toluca, Estado
de México, México, 2009.
3. P. Esquivel, and A. R. Messina, "Time-Dependent EOF Analysis of Power
System Measurements," International Council on Large Electric Systems
(CIGRÉ), Ontario, Toronto Canada, 2009.
4. F. Lezama Z., A. L. Rios, P. Esquivel, and A. R. Messina, "A Hilbert-Huang
Based Approach for On-Line Extraction of Modal Behavior from PMU
Data," North American Power Symposium (NAPS), Starkville, Mississippi
USA, 2009.
1.7.2 Refereed journals
1. A. R. Messina, P. Esquivel and F. Lezama, "Time-Dependent Statistical
Analysis of Wide-Area Time-Synchronized Data," Journal of Mathematical
Problems in Engineering
(MPE), Special Issue on Short Range Phenomena:
Modeling, Computational Aspects, and Applications, 2010.
12
2. P. Esquivel, and A. R. Messina, "Spatio-Temporal Decorrelation based in
Statistical Analysis from Wide-Area Measurements," to be submitted to
International Journal of Electric Power System Research
, 2010.
1.7.3 Book chapters
1. P. Esquivel, E. Barocio, M. A. Andrade, and F. Lezama, "Complex Empirical
Orthogonal Function Analysis of Power-System Oscillatory Dynamics", in:
"Analysis of Nonlinear and Non-stationary Oscillations: A Time-Frequency
Perspective," Book Chapter, Springer Verlag, 2009.
1.8 References
[1] A. R. Messina, "Inter-area Oscillations In Power Systems: A Nonlinear and
Nonstationary Perspective," Springer Verlag, 2009.
[2] T. Hiyama, D. Kojima, K. Ohtsu, and K. Furukawa, "Eigenvalue-based Wide area
stability monitoring of power systems," Control Engineering Practice, vol. 13, pp.
1515-1523, 2005.
[3] D. Hu, and V. Venkatasubramanian, "New wide-area algorithms for detecting
angle instability using synchrophasors," Proceedings of the Western Protective Relay
Conference
, pp. 1-37, 2006.
[4] R. Avila-Rosales, and J. Giri, "Wide-area monitoring and control for power system
grid security," 15
th
[5] D. Y. Wong, G. J. Roger, B. Porretta, and P. Kundur, "Eigenvalue analysis of very
large power systems," IEEE Transactions on Power Systems, vol. 3, no. 2, pp. 472-
480, May 1988.
PSCC, Liege
, pp. 1-7, 2005.
[6] P. Kundur, "Power system stability and control," New York: McGraw-Hill, Inc.,
1994.
[7] V. Vittal, N. Bhatia, and A. A. Fouad, "Analysis of the interarea mode phenomena
in power systems following large disturbances," IEEE Transactions on Power
Systems
, vol. 6, no. 2, pp. 15151521, Nov., 1991.
[8] D. Ruiz-Vega, A. R. Messina, and G. Enriquez-Harper, "Analysis of inter-area
oscillations via non-linear time series analysis techniques," 15
th
[9] J. Arroyo, E. Barocio, R. Betancourt, and A. R. Messina, "A bilinear analysis
technique for detection and quantification of nonlinear modal interaction in power
systems", in Proc. 2006 IEEE Power Engineering Society General Meeting, pp. 18.
PSCC
, pp. 1-7,
Aug., 2005.
13
[10] Pacific Northwest National Laboratory, U. S. Department of Energy, "The role of
synchronized wide area measurements for electrical power grid operations,"
Networked Embedded Control for Cyber Physical Systems
(HCSS-NEC4CPS), 2006.
[11] H. Ukai, K. Nakamura, and N. Matsui, "DSP- and GPS-based synchronized
measurement system of harmonics in wide-area distribution system," IEEE
Transactions on Industrial Electronics
, vol. 50, no. 6, pp. 1159-1164, Dic., 2003.
[12] A. R. Messina, and V. Vittal, "Extraction of dynamic patterns from wide-area
measurements using empirical orthogonal functions," IEEE Transactions on Power
Systems
, vol. 22, no. 2, pp. 682-692, May 2007.
[13] A. R. Messina, V. Vittal, D. Ruiz-Vega, and G. Enriquez-Harper, "Interpretation
and visualization of Wide-area PMU measurements using Hilbert analysis," IEEE
Transactions on Power Systems
, vol. 21, pp. 1763-1771, no. 4, Nov., 2006.
[14] A. Hannachi, I. T. Jolliffe, and D. B. Stephenson, "Review. Empirical orthogonal
functions and related techniques in atmospheric science: a review," International
Journal of Climatology
, vol. 27, pp. 1119-1152, May 2007.
[15] E. Campbell, "A review of methods for statistical climate forecasting," Technical
Report
06/134, Aug., 2006.
[16] D. J. Trudnowski, J. M. Johnson, and J. F. Hauer, "Making Prony analysis more
accurate using multiple signals," IEEE Transactions on Power Systems, vol. 14, no. 1,
pp. 226-231, Feb., 1999.
[17] R. F. Jeffers, "Techniques for wide-area state estimation in power systems," A
thesis submitted to the faculty of the Virginia polytechnic institute and state
university, 2006.
[18] P.-F. Bach, "Wide area control and long distance transfer of electric power," IEA
Workshop
, pp. 1-7, Nov., 2004.
[19] J. F. Hauer, and J. G. DeSteese, "A tutorial on detection and characterization of
special behavior in large electric power systems," Pacific Northwest National
Laboratory
, July 2004.
[20] J. Arroyo, R. Betancourt, A.R. Messina, E. Barocio, "Development of Bilinear
Power System Representations for Small Signal Stability Analysis", Electric Power
Systems Research
, vol. 77, no. 10, pp. 12391248, Aug., 2007.
[21] P. Esquivel and A. R. Messina, "Complex Empirical Orthogonal Function Analysis
of Wide-Area System Dynamics," Power and Energy Society, General Meeting IEEE,
pp. 1-7, July 2008.
[22] P. Holmes, J. L. Lumley, and G. Berkooz, "Turbulence, coherent structures,
dynamical systems and symmetry," New York: Cambridge University Press, 1996.
14
[23] J. Arroyo, E. Barocio, R. J. Betancourt, A. R. Messina, "Quantifying nonlinearity in
power systems using normal forms theory and higher-order statistics", In Proc.
2005
IEEE PES General Meeting
, pp. 18.
[24] J. F. Hauer, N. B. Bhatt, K. Shah, and S. Kolluri, " Performance of WAMS East in
providing dynamical information for the North East blackout of august 14, 2003,"
The Bonneville Power Administration
, pp. 1-6, 2005.
[25] T. P. Barnett, "Interaction of the Monsoon and pacific trade wind system at
interannual time scales. Part 1: the equatorial zone," Monthly weather review, vol.
111, pp. 756-773, 1983.
[26] M. A. Merrifield, and R. T. Guza, "Notes and correspondence. Detecting
propagating signals with complex empirical orthogonal functions: a cautionary
note," Journal of physical oceanography, pp. 1628-1633, Oct., 1990.
[27] J. Terradas, R. Oliver, and J. L. Ballester, "Application of statistical techniques to
the analysis of solar coronal oscillations," The Astrophysical Journal, vol. 614, pp.
435-447, Oct., 2004.
[28] D. H. Chambers, R. J. Adrian, P. Moin, D. S. Stewart, and H. J. Sung, "Karhunen-
Loève expansion of Burgers´ model of turbulence," Physics Fluids, vol. 31, no. 9,
pp. 2573-2582, Sep., 1988.
[29] S. Han, and B. Feeny, "Application of proper orthogonal decomposition to
structural vibration analysis," Mechanical systems and signal processing, vol. 17, no.
5, pp. 989-1001, 2003.
[30] S. Toh, "Statistical model with localized structures describing the spatial-temporal
chaos of Kuramoto-Sivashinsky equation," Journal of the Physical Society of Japan,
vol. 56, no. 3, pp. 949-962, Mar., 1987.
[31] D. J. Trudnowski, "Estimating Electromechanical Mode Shape from
Synchrophasor measurements," IEEE Transactions on Power Systems, vol. 23, no. 3,
pp. 1188-1195, August 2008.
[32] R. P. Pacheco, and V. Steffen Jr., "On the identification of non-linear mechanical
systems using orthogonal functions," International Journal of Non-linear Mechanics,
vol. 39, pp. 1147-1159, 2004.
15
Chapter 2
Empirical Orthogonal Function
Analysis
Empirical orthogonal function (EOF) analysis is a statistical method of finding
optimal distributions of energy from an ensemble of multidimensional
measurements. The essential idea is to generate an optimal basis for the
representation of an ensemble of data collected from measurements or
numerical simulations of a dynamic system.
Given an ensemble of measured data, the technique yields an orthogonal
basis for the representation of the ensemble, as well as a measure of the relative
contribution of each basis function to the total energy with no a priori
assumption on either spatial or temporal behavior. The method can serve two
purposes, namely order reduction by projecting high-dimensional data into a
lower dimensional space and feature extraction by revealing relevant, but
unexpected, structure hidden in the data.
The following sections provide a review of some aspects of the qualitative
theory of empirical orthogonal functions that are needed in the analysis of low-
dimensional models derived from the technique.
We start by introducing the method in the context of statistical correlation
theory.
16
2.1 Notation
Most of the notation is standard. The Euclidean norm of a vector
x
is denoted
by
x
. Other norms are indicated by double bars:
.
. Vectorial quantities are
denoted by boldface letters and scalar quantities by italic letters.
The inner (scalar, o dot) product of two elements, u, v in
n
is given by
¦
n
i
i
i
v
u
1
)
,
(
v
u
2
L
denotes the space of square integrable real or complex-valued functions.
Thus,
]
1
,
0
[
2
L
denotes the space of functions defined over the unit interval
1
0
d
d x
. The inner product space is defined by
³
:
dx
x
g
x
f
g
f
)
(
)
(
)
,
(
*
where * denotes the complex conjugate.
Other symbols are defined in the text.
2.2 Theoretical Development
The proper orthogonal decomposition (POD) method is an optimal technique of
finding a basis that spans an ensemble of data, collected from an experiment or
numerical simulations [1-4]. More precisely, assume that
)
,
(
k
j
t
x
u
,
n
j
,
,
1
"
,
N
k
,
,
1
"
, denotes a sequence of observations on some domain
:
x
where
x
is a vector of spatial variables, and
]
,
0
[ T
t
k
is the time at which the
observations are made.
Without loss of generality, the time average of the time sequence
¦
N
k
k
k
m
t
u
N
t
u
u
1
)
,
(
1
)
,
(
)
(
X
X
X
is assumed to be zero [3]. Generalizations to this approach are discussed below.
17
The POD procedure determines empirical orthogonal functions (EOFs),
f
,...
1
,
)
(
i
x
i
, such that the projection onto the first p EOFs (a low order
representation)
N
k
n
j
x
t
a
t
x
u
p
i
i
i
k
j
,
,
1
,
,...,
1
,
)
(
)
(
)
,
(
^
1
"
¦
M
(2.1)
is optimal in the sense that the average least-squares truncation error,
j
H
N
p
x
t
a
t
x
u
p
i
i
i
k
j
j
d
¦
,
2
1
)
(
)
(
)
,
(
M
H
(2.2)
is minimized for any
N
p
d
, where
!
.
denotes ensemble average,
2
/
1
,
!
f
f
f
, and . denotes the
2
L
norm over
: . The
i
a
's are time-
dependent coefficients of the decomposition to be determined so that (2.1)
results in a maximum for (2.2). These special orthogonal functions are called the
proper orthogonal modes (POMs) of the reduced basis for
)
,
(
k
j
t
x
u
. The
computation of the EOFs in the presence of measurement error is given in
Appendix A.
Following [1], assume that the field
)
,
(
k
j
t
x
u
is decomposed into a mean
value
)
,
(
t
x
u
j
m
, and a fluctuating part
)
,
(
j
t
x
u
)
,
(
)
,
(
)
,
(
k
j
k
j
m
k
j
t
x
u
t
x
u
t
x
u
(2.3)
More formally, let
2
L
denote the space of square integrable functions. It
follows that, a normalized basis function
M is optimal if the average projection
of
u
onto
M is maximized, i.e., [2]
1
to
subject
|
)
),
,
(
(
|
max
2
2
])
1
,
0
([
2
t
x
u
k
j
L
(2.4)
where the inner product is defined as
V
U
V
U
u
!
,
, and
¦
!
m
j
j
T
1
2
2
,
M
M
M
M
M
M
18
The optimization problem can be recast in the form of a functional for the
constrained variational [3]
1
)
1
(
|
)
),
,
(
(
|
]
[
2
2
t
x
u
J
k
j
(2.5)
where
O
is a Lagrange multiplier. A necessary condition for the extrema is that
the Gateaux derivative vanishes for all variations
])
1
,
0
([
2
L
G\
M
,
G
. This
can be expressed as
>
@
)
(
,
0
2
0
:
L
d
dJ
(2.6)
Consider now the Hilbert space of all pairs
g
f ,
where
f
and g are
functions of
]
1
,
0
[
2
L
, i.e., square integrable functions of the space variable
x
on
the interval
]
1
,
0
[
, where .,. denotes the standard inner product on
2
L
defined
by
dx
x
g
x
f
g
f
³
1
0
*
)
(
)
(
,
(2.7)
and
dx
³
:
!
2
2
,
(2.8)
where
: is the domain of interest over which
)
(x
u
and are defined and the
asterisk
(*)
denotes conjugate transpose.
It follows immediately from (2.6) that
@
@
u
u
u
u
d
dJ
d
dJ
,
,
,
Re
2
,
,
,
0
]
[
0
0
(2.9a)
where use has been made of the inner product properties.
1
Given a function to maximize,
)
(P
f
, subject to the constraint
0
)
(
P
g
, the Lagrange function
can be defined as
)
(
)
(
)
,
(
P
g
P
f
j
P
F
O
O
.
19
Noting from the commutability of the averaging operator and spatial
integral, that
0
)
(
)
(
'
)
'
(
)
'
(
)
(
)
(
)
(
'
)
'
(
)
'
(
)
(
)
(
,
,
,
*
1
0
1
0
*
1
0
*
1
0
*
1
0
*
»¼
º
«¬
ª
³ ³
³
³
³
dx
x
x
dx
x
x
u
x
u
dx
x
x
dx
x
u
x
dx
x
x
u
u
u
\
OM
M
\
M
O
M
\
\
M
O
M
\
(2,9b)
the condition for the extrema reduces to
)
(
'
)
'
(
)
(
)
(
*
x
dx
x
x
u
x
u
³
:
(2.10)
Equation (2.10) has a finite number of orthogonal solutions
)
(x
i
M
(the
proper orthogonal modes) with corresponding real and positive eigenvalues
i
O
.
They are consequently called empirical orthogonal functions.
Defining
'
)
'
(
)
,
'
(
)
,
(
*
dx
x
t
x
u
t
x
u
³
:
R
(2.11)
where
¦
N
i
k
k
t
x
u
t
x
u
N
x
x
R
1
)
,
'
(
)
,
(
1
)
'
,
(
(2.12)
The problem of minimizing (2.9) becomes that of finding the largest
eigenvalue of the eigenvalue problem
R
, subject to
1
2
.
In practice, the observations that form the data are only available at discrete
spatial grid points herein called snapshots. In this case, the kernel
)
'
,
(
x
x
R
can
be written in matrix notation as [5-6]
»
»
»
¼
º
«
«
«
¬
ª
)
,
(
)
,
(
)
,
(
)
,
(
1
1
1
1
n
n
n
n
x
x
R
x
x
R
x
x
R
x
x
R
#
%
#
)
x'
R(x,
where n indicates the number of measurement points, and
20
n
j
i
t
x
u
t
x
u
N
x
x
R
N
k
k
j
k
i
j
i
,...,
1
,
,
)
,
(
)
,
(
1
)
,
(
1
¦
(2.13)
In other words, the optimal basis is given by the eigenfunctions
i
M
of (2.13)
whose kernel is the autocorrelation function
)
,
'
(
)
,
(
)
'
,
(
k
j
k
j
t
x
u
t
x
u
x
x
R
.
2.3 Discrete Domain Representation
As highlighted in the previous section, time series are usually recorded in
discrete form even though the underlying process itself is continuous. For
discretely sampled measured data, the integral time-average can be
approximated by a sum over the set of sampled data points [1]. In this case, the
vectors
n
j
t
x
u
t
x
u
t
x
u
t
x
T
N
j
j
j
k
j
j
j
,...,
1
,
)]
,
(
...
)
,
(
)
,
(
[
)
,
(
2
1
u
x
(2.14)
represent a set of snapshots obtained from the observed data at
n
locations
}.
,...,
,
{
2
1
n
x
x
x
The set of data can then be written as the N×n-dimension ensemble matrix,
X [6]
»
»
»
¼
º
«
«
«
¬
ª
)
,
(
)
,
(
)
,
(
)
,
(
]
[
1
1
1
1
1
N
n
N
n
n
t
x
u
t
x
u
t
x
u
t
x
u
#
%
#
x
x
X
(2.15)
where each column corresponds to the measured response at a specific time.
Typically,
N
n
z
, so X is generally rectangular. Under these assumptions,
the actual integral (2.10) can be written as
C
, where
¦
N
k
k
j
k
i
ij
t
x
u
t
x
u
N
1
)
,
(
)
,
(
1
C
.
Assuming the PODs to be of the form
¦
N
l
l
i
l
i
w
1
x
M
, where
i
l
w
is a
coefficient to be determined, the problem of minimizing (2.2) can be recast as
the problem of finding the largest eigenvalue of the linear equation
M
O
M
C
(2.16)
21
where
C
is the autocorrelation (covariance) matrix defined as
)
(
)
(
1
1
1
1
2
1
2
2
2
1
2
1
2
1
1
1
m
i
T
n
i
m
i
n
T
n
T
n
T
n
n
T
T
T
n
T
T
T
T
u
x
u
x
N
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
N
N
»
»
»
»
»
¼
º
«
«
«
«
«
¬
ª
¦
#
%
#
#
X
X
C
(2.17)
The resulting covariance matrix
C
, is a real, symmetric (
ji
ij
C
C
), positive
and semi-definite matrix. Consequently, it possesses a set of orthogonal
eigenvectors
n
i
i
,...,
1
,
M
, i.e.,
¯
®
z
j
i
j
i
ij
j
T
i
,
0
,
G
M
M
Using standard linear algebra techniques, the covariance matrix can be
expressed in the form
T
V
U
C
(2.18)
where
U
and
V
are the matrices of right and left eigenvectors and
]
[
diag
2
1
n
O
O
O
are the ordered eigenvalues of the correlation matrix,
n
t
t
2
1
.
The eigenvalues computed from (2.18) are real, nonnegative and can be
ordered such that
0
2
1
t
t
t
t
n
O
O
O
. The eigenvectors of
C
are called the
proper orthogonal modes (POMs), and the associated eigenvalues are called the
proper orthogonal values (POVs).
For practical applications, the number of snapshots
N
can be rather large,
leading to a very large eigenvalue problem. There are at least two methods to
solve the eigenvalue problem (2.16) [6]: the direct method and the method of
snapshots.
The direct method attempts to solve the eigenvalue problem involving the
n
N
u
matrix directly using standard numerical techniques. This can be
computationally intensive if the number of observations is larger than the
number of observing locations or grid points.
22
The method of snapshots, on the other hand, is based on the fact that the
data vectors
i
u
and the POD modes
l
span the same linear space. The latter is
explored here.
The next section describes in more detail the nature of the approximation
employed here to construct the statistical representation.
2.4 The Method of Snapshots
The method of snapshots is based on the fact that the data vectors
i
u
and the
POD modes span the same linear space [6]. In this approach, we choose the
eigenfunctions
M to be a linear combination of the snapshots:
¦
N
l
i
i
l
i
w
1
x (2.19)
where the coefficients
i
l
w
are to be determined such that
u
maximizes
¦
!
n
j
j
x
N
1
,
,
1
max
(2.20)
These l functions are assembled into an
N
n
u matrix, , known as the
modal matrix. Equation (2.19) in matrix form becomes
XW
(2.21)
where
»
»
»
¼
º
«
«
«
¬
ª
p
p
p
p
n
n
n
n
n
M
M
M
2
1
;
»
»
»
¼
º
«
«
«
¬
ª
p
p
p
p
n
n
n
n
n
x
x
x
2
1
X
;
»
»
»
¼
º
«
«
«
¬
ª
p
p
p
p
n
n
n
n
n
w
w
w
2
1
W
in which
»
»
»
»
»
¼
º
«
«
«
«
«
¬
ª
1
1
2
1
1
1
N
w
w
w
#
w
,
»
»
»
»
»
¼
º
«
«
«
«
«
¬
ª
2
2
2
2
1
2
N
w
w
w
#
w
, ...,
»
»
»
»
»
¼
º
«
«
«
«
«
¬
ª
n
N
n
n
n
w
w
w
#
2
1
w
23
Substituting the expression (2.19) into the eigenvalue value problem (2.16)
gives
¦
¦
N
l
i
i
l
N
l
i
i
l
w
w
1
1
x
x
C
(2.22)
where
)
,
(
)
/
1
(
j
i
ij
u
u
N
C
.
This can be written as the eigenvalue problem of dimension n
:
CW
(2.23)
where
@
n
w
w
w
2
1
w
and
is a diagonal matrix storing the eigenvalues
i
O
of the covariance matrix
C. In words, the first-order necessary optimality condition for
M to provide a
maximum in (2.20) is given by (2.16). This completes the construction of the
orthogonal set
^
`
n
M
M
M
2
1
.
Once the modes are found using these equations, the flow field can be
reconstructed using a linear combination of the modes
¦
f
1
)
(
)
(
)
,
(
^
i
i
i
k
j
k
x
t
a
t
x
u
(2.24)
for some
)
(t
a
i
, where the
)
(t
a
i
are the time-varying amplitudes of the POD
modes
)
(x
i
.
The truncated POD of
u
is
x
t
a
t
x
u
p
i
i
i
k
j
k
¦
1
)
(
)
(
)
,
(
^
(2.25)
where p is the number of dominant modes and is an error function. Having
computed the relevant eigenmodes, the temporal behavior of each mode can be
evaluated as the inner product of the eigenmode (the POD mode,
i
) and the
24
original data
.
To ensure uniqueness of the solution, the normalization condition
of
i
i
M
M ,
=1 is imposed.
The temporal coefficients are then
expressed as
i
i
i
i
a
M
M
M
,
,
x
(2.26)
Note that the temporal modes are uncorrelated in time, i.e.,
j
ij
j
i
t
a
t
a
))
(
),
(
(
, where
1
ij
G
for
j
i
, 0
else, and that the system (2.25) is
optimal in the sense that minimizes the error function
¦
¦
p
l
p
i
i
i
k
j
x
t
a
t
x
u
1
1
)
(
)
(
)
,
(
)
(
M
M
H
(2.27)
It should also be stressed that no conditions are imposed on the data set: the
data can be a sample of a stationary process or a sample of a nonstationary
process.
Equation (2.24) is called the Karhunen-Loève decomposition and the set
j
M
are called the empirical basis [5].
2.5 Energy Relationships
The use of the POD method leads naturally to a discussion of truncation
criteria. Several techniques to derive truncated expansions have been proposed
in the literature. Here, we choose to reduce the residual terms,
R
, such that
the mean square value
¸
¸
¸
¸
¹
·
¨
¨
¨
¨
©
§
¦
N
i
i
p
O
R
1
1
(2.28)
be as small as possible.
25
Defining the total energy E , by
)]
,
(
Trace[
m
i
E
x
x
R
, one obtains
¦
p
i
i
E
1
(2.29)
The associated percentage of total energy contributed by each mode can
then be expressed as
¦
n
i
i
i
i
E
1
(2.30)
Typically, the order p of the reduced basis
M such that a predetermined
level of the total energy E , of the ensemble of snapshots is captured, i.e., 99%.
The p -dominant eigenfunctions are then obtained as
}
)
(
:
)
(
min{
arg
to
subjet
%
99
100
)
(
%
0
1
1
E
p
E
p
E
p
E
n
j
j
p
i
i
t
u
¦
¦
(2.31)
where
0
E is an appropriate energy level, and
100
)
(
%
0
d
d
p
E
is the percentage
of energy that is captured by the POD basis. By neglecting modes
corresponding to the small eigenvalues a reduced-order model can be
constructed. In Appendix B, an optimal linear predictor to determine the lower
integer p is proposed.
The key advantage of this technique is that allows extracting modal
information from short and often noisy time series without prior knowledge of
the dynamics affecting the time series.
2.6 Interpretation of EOFs Using Singular Value
Decomposition
A useful alternative method for estimating modal characteristics can be
developed based on the analysis of the response matrix X in (2.15).
26
Before outlining the procedure for singular value decomposition (SVD)
analysis, we introduce some background information on singular value
analysis.
2.6.1 Singular Value Decomposition
Let A be a real
n
m
u
matrix. The singular value decomposition theorem states
that A can be decomposed into the following form [7]
T
V
U
A
66
(2.32)
where
]
[
col
2
1
m
u
u
u
U
is an
m
m
u
orthonormal matrix (
1
U
V
T
), is
an
n
m
u
pseudo-diagonal and semi-positive definite matrix with diagonal
entries containing the singular values, and
]
[
col
2
1
n
v
v
v
V
is an
n
n
u
orthonormal matrix (
1
U
U
T
). The columns of
U
and
V
are called left and
right singular vectors for
A
.
Matrix has the form
n
m
m
»
»
»
»
¼
º
«
«
«
«
¬
ª
for
0
0
0
0
0
0
0
0
0
0
0
0
2
1
#
%
#
%
#
#
#
66
or
n
m
m
!
»
»
»
»
»
»
»
»
»
¼
º
«
«
«
«
«
«
«
«
«
¬
ª
for
0
0
0
0
0
0
0
0
0
0
0
0
2
1
#
%
#
#
#
%
#
#
66
Throughout this research, we will, consider only the case when
n
m
!
. The
diagonal entries of , i.e., the
i
ii
V
, can be arranged to be nonnegative and in
order of decreasing magnitude
0
2
1
t
t
t
t
m
V
V
V
.
27
Equivalently, we can express the model as
@
»
»
»
»
»
»
»
»
¼
º
«
«
«
«
«
«
«
«
¬
ª
»
»
»
»
¼
º
«
«
«
«
¬
ª
T
n
T
k
T
k
T
k
m
k
k
v
v
v
v
u
u
u
u
#
#
#
#
%
#
1
1
1
1
1
0
0
0
0
0
0
A
(2.33)
or
@
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
A
A
V
V
D
U
U
A
where D is the diagonal matrix of nonzero singular values, and
U
and
V
are
the matrices of left and right singular vectors, respectively, corresponding to the
nonzero singular values.
It is clear that only the first
r
of the
u
's and
v
'
s make any contribution to A
, and can be expressed as an outer product expansion
¦
r
i
T
i
i
i
1
)
( v
u
A
(2.34)
where the vectors
i
u
and
i
v
are the columns of the orthogonal matrices
U
and
V
, respectively. Techniques to compute the POMs based on singular value
decomposition are next discussed.
2.6.2 Relation with the Eigenvalue Decomposition
An interesting interpretation of the POD modes can be obtained from the
singular value analysis of the response matrix X .
Using the notation in section 2.2 let the response matrix be given by
»
»
»
¼
º
«
«
«
¬
ª
)
,
(
)
,
(
)
,
(
)
,
(
1
1
1
1
N
n
N
n
t
x
u
t
x
u
t
x
u
t
x
u
#
%
#
X
(2.35)
28
where the columns correspond to a response at time. The SVD of the response
matrix X may be written in compact form as
T
V
U
X
(2.36)
where
U
is an orthonormal
N
N
u
matrix whose columns are the left singular
vectors of X , is
n
N
u
matrix containing the singular values of X along the
main diagonal and zeros elsewhere, and
V
is an
n
n
u
orthogonal matrix whose
columns correspond to the right singular vectors of X . The response matrix, X ,
is complex and symmetric and possesses a set of orthogonal singular vectors
with positive singular values.
In terms of the notation above for SVD, it can be seen directly from (2.18)
that the correlation matrix defined previously is given by
T
T
T
U
U
V
U
V
U
XX
2
66
6
6
(2.37)
and
T
T
T
V
V
V
U
V
U
X
X
2
(2.38)
Hence, (2.36) becomes
»
»
»
»
»
¼
º
«
«
«
«
«
¬
ª
»
»
»
»
¼
º
«
«
«
«
¬
ª
2
2
2
2
1
2
1
n
n
T
%
%
U
U
U
XX
(2.39)
It follows immediately from (2.37)-(2.39), that the singular values of X are
the square roots of the eigenvalues of
T
XX or
X
X
T
[2,8]. In addition, the left
and right eigenvectors of X are the eigenvectors of
T
XX and
X
X
T
, respectively.
Also of interest, the trace of (2.39) is given by
¦
r
i
i
T
1
2
XX
The POMs, defined as the eigenvectors of the sample correlation matrix
C
are thus equal to the left singular vectors of X . The POVs, defined as the
29
eigenvalues of matrix
C
are the squares of the singular values divided by the
number of samples
N
.
2.6.3 Geometric Interpretation Using SVD
The SVD of matrix X , seen as a collection of column vectors provides
important insight into the oriented energy distribution of this set of vectors.
The energy of the vector sequence
k
x
building a matrix
m
N
u
is defined by
the square root of sums of squares of all the elements (Frobenius norm)
¦
¦
¦
¦
N
j
r
k
k
r
k
k
m
n
j
n
F
m
N
r
N
E
1
1
1
2
1
2
)
(
2
.
)
,
min(
,
)
1
(
)
(
x
X
C
(2.40)
Thus, the first columns of the matrix
V
give an optimal orthonormal basis
for approximating the data [5]. So that the energy of a vector sequence is equal
to the energy in its singular spectrum
2
2
1
,
,
r
V
V
or alternative in its
eigenspectrum
r
,
,
1
. A key property of the SVD is that the extreme in this
oriented energy distribution occur at each left singular direction.
2.7 Relation Between POMs and Normal Modes
Significant physical insight into the physical interpretation of POMs can be
obtained from the study of mechanical systems.
In this section we review the theory behind the notion of normal modes in
mechanical systems and show that, under certain conditions, the POMs from
correlation analysis converge to the normal modes.
2.7.1 The Notion of Normal Modes
Consider an
m
-degree-of-freedom nonlinear mechanical system. The equation
of motion describing the forced vibration of the system may be written as
follows:
,
f
y
y
(2.41)
30
where
is the mass matrix,
is the stiffness matrix, and
y is the
displacement vector. Using the transformation
x
y
2
1
, Eq. (2.41) takes the
form
,
^
f
x
x
(2.42)
where
2
1
2
1
. We also assume that
1
exists.
The general solution that represents the free natural response of the
dynamical systems (2.41) is [9]
,
)
(
^
)
(
)
(
1
2
0
t
t
t
i
0
f
x
e
f
x
(2.43)
where
)
(t
f
x
is a solution of the free motion
)
0
^
(
f
, and can be expressed in
terms of the normal modes
n
v
of the free motion
¦
¦
m
n
m
n
n
n
n
n
n
n
n
n
n
t
t
t
t
1
1
,
)
sin(
))
sin(
^
)
cos(
^
(
)
(
v
v
D
C
A
x
f
(2.44)
or,
,
)
(
)
(
1
¦
p
n
n
n
t
t
v
a
x
f
(2.45)
where
,
)
sin(
)
(
n
n
n
n
t
t
A
a
(2.46)
are the temporal amplitudes, and the
n
v
is the right eigenvector associated with
the
n
-th system mode; the
n
A
are constants depending on the initial conditions.
The modal shapes are determined by the eigenvector
n
v
.
The function
)
sin(
n
n
t
operates simply as a mapping scale. The resulting
vibration modes of the oscillations are stationary and one thus refers to such
vibrations as stationary oscillations. Such a natural vibration mode repeats itself
with a period
n
n
2
T
, where
n
are the natural frequencies (in
s
rad
).
31
Since the determinant of the
m
m
u
matrix
)
(
2
0
is a polynomial of
degree
m
2
in
Z
it has
m
2
roots and therefore there are
m
2
normal modes.
Equation (2.45) can be equivalently expressed as
,
)
(
)
(
)
(
)
(
)
(
2
1
2
1
2
22
21
1
12
11
t
t
a
t
a
t
a
v
v
v
v
v
v
v
v
v
t
m
mm
m
m
m
m
a
x
f
v
»
»
»
»
¼
º
«
«
«
«
¬
ª
»
»
»
»
¼
º
«
«
«
«
¬
ª
#
#
%
#
#
(2.47)
where
@
m
v
v
v
v
2
1
is the modal matrix of
m
normal modes
n
v
and
@
T
m
t
a
t
a
t
a
t
)
(
)
(
)
(
)
(
2
1
a
is the vector of modal coordinates with elements
)
sin(
)
(
n
n
n
n
t
t
A
a
representing the time modulation of the natural modes
[10].
2.7.2 Physical Interpretation of POMs in Terms of the Eigenvalues
The POMs have a physical interpretation of interest in terms of the eigenvalues
of the correlation matrix. Let the response matrix,
, of the mechanical system
be expressed in the form
.
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
1
1
2
1
1
1
1
1
2
1
1
1
1
1
0
1
2
0
1
1
0
»
»
»
»
»
»
»
»
¼
º
«
«
«
«
«
«
«
«
¬
ª
¦
¦
¦
¦
¦
¦
¦
¦
¦
m
n
nm
N
n
m
n
n
N
n
m
n
n
N
n
m
n
nm
n
m
n
n
n
m
n
n
n
m
n
nm
n
m
n
n
n
m
n
n
n
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
#
%
#
#
(2.48)
Now defining the vectors
a
as
1
u
m
arrays of the functions
)
(t
a
n
evaluated
N
times, one can form an ensemble matrix
,
@
,
)
sin(
)
sin(
)
sin(
)
sin(
)
sin(
)
sin(
)
sin(
)
sin(
)
sin(
2
2
2
1
1
1
1
2
1
2
2
1
1
1
1
0
2
0
2
2
1
0
1
1
2
1
»
»
»
»
¼
º
«
«
«
«
¬
ª
m
N
m
m
N
N
m
m
m
m
m
m
m
t
t
t
t
t
t
t
t
t
A
A
A
A
A
A
A
A
A
a
a
a
#
%
#
#
(2.49)
from which it follows that
32
T
m
m
T
T
T
m
T
T
T
v
v
v
v
v
v
A
AV
X
a
a
a
2
2
1
1
2
1
]
[
, (2.50)
where the
n
v
are the column vectors formed by the natural modes of the
system. These are orthogonal to each other in the metrics of the mass and
stiffness matrices.
It thus follows that the correlation matrix
C
, may be written as
.
)
(
)
(
1
1
2
2
1
1
2
2
1
1
T
m
m
T
T
T
T
m
m
T
T
T
N
N
v
v
v
v
v
v
C
a
a
a
a
a
a
(2.51)
To verify that a modal vector is a POM, we post-multiply the matrix
C
by
n
v
to obtain
.
v
)
(
)
(
1
)
(
1
^
2
2
1
1
2
2
1
1
n
T
m
m
T
T
T
T
m
m
T
T
n
T
n
N
N
v
v
v
v
v
v
v
Cv
C
a
a
a
a
a
a
(2.52)
Using the orthogonality condition
1
n
T
n
v
v
for
k
n
, the above expression
simplifies to
@
.
1
^
1
1
m
k
T
m
k
k
T
k
k
T
N
v
v
v
C
a
a
a
a
a
a
(2.53)
Further, noting that
,
)
sin(
)
sin(
)
(
sin
,
)
sin(
)
sin(
1
1
2
2
2
1
1
¦¦
¦¦
m
l
m
k
k
n
l
n
k
n
k
k
l
k
k
k
k
T
k
m
l
m
k
k
n
l
n
k
n
n
k
T
n
t
t
t
N
t
t
A
A
A
A
A
A
A
A
a
a
A
A
a
a
v
v
(2.54)
Equation (2.53) can now be rewritten in the form
.
^
k
k
k
T
k
v
v
C
a
a
(2.55)
In such case, as the frequencies of the modes are distinct, each term
N
n
k
T
n
v
a
a
will tend to zero as
N
tends towards infinity, except for the term
N
k
k
T
k
v
a
a
which converges approximately to
k
v
[11]. In other words, the
eigenvectors (POMs) of
C
converge to the modal vector as expected.
33
In what follows we provide a similar interpretation of POMs in the
framework of singular value decomposition analysis.
2.7.3 Physical Interpretation of POMs Using Singular Value
Decomposition
Consider again the system (2.41). When the displacements are sampled
N
times, we can form a displacement-history array. The observed field is then
represented by the data matrix:
@
,
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
^
^
^
^
2
1
2
2
2
2
1
1
1
2
1
1
2
1
»
»
»
»
¼
º
«
«
«
«
¬
ª
N
m
N
N
m
m
m
t
x
t
x
t
x
t
x
t
x
t
x
t
x
t
x
t
x
#
%
#
#
x
x
x
(2.56)
where
¦
m
n
n
n
k
k
t
a
v
x
1
)
(
)
(
^
^
^
.
From (2.56), it follows that the measurement matrix can be expressed in the
form
,
^
^
)
(
^
^
)
(
^
^
)
(
^
^
)
(
^
^
)
(
^
^
)
(
^
^
)
(
^
^
)
(
^
^
)
(
^
^
v^
1
1
2
1
1
1
2
1
2
2
1
1
2
1
1
1
2
1
1
1
1
T
m
n
mn
N
n
m
n
n
N
n
m
n
n
N
n
m
n
mn
n
m
n
n
n
m
n
n
n
m
n
mn
n
m
n
n
n
m
n
n
n
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
v
t
a
»
»
»
»
»
»
»
»
¼
º
«
«
«
«
«
«
«
«
¬
ª
¦
¦
¦
¦
¦
¦
¦
¦
¦
#
%
#
#
(2.57)
where
,
^
^
^
^
^
^
^
^
^
2
1
2
22
21
1
12
11
v^
»
»
»
»
¼
º
«
«
«
«
¬
ª
mm
m
m
m
m
v
v
v
v
v
v
v
v
v
#
%
#
#
(2.58)
is the modal matrix.
Defining
34
,
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
)
(
^
^
2
1
2
2
2
1
2
1
2
1
1
1
»
»
»
»
¼
º
«
«
«
«
¬
ª
N
m
m
m
N
N
t
a
t
a
t
a
t
a
t
a
t
a
t
a
t
a
t
a
#
%
#
#
(2.59)
the system (2.57) can be rewritten as
,
^
]
)
(
^
)
(
^
)
(
^
][
^
^
^
[
^
v^
2
1
2
1
T
m
m
t
t
t
v
v
v
a
a
a
(2.60)
Further insight into the nature of this model can be obtained by expressing
(2.60) in the form
,
]
^
][
[
^
v^
T
T
9
U
Q
(2.61)
with the orthogonality conditions
U
U
UU
*
*
, where
v^
is the modal matrix,
,
is an
)
(
m
m
u
identity matrix, is an
))
(
(
m
N
m
u
matrix full of zeros, and Q
is an
))
(
(
m
N
N
u
matrix.
We remark that matrix Q does not influence
^
since it is multiplied by a
matrix of zeros [8]. Note also that matrix Q affects only matrix
^
and that the
scaling not affect the coefficients of
^
.
In order for (2.60) to admit a decomposition of the form (2.61), matrices
v^
and
]
^
[
Q should be orthogonal.
According to (2.61), the response matrix may then be written as
,
]
)
^
(
^
][
)
^
(
[
]
^
^
^
[
^
1
2
1
T
n
n
m
diag
diag
Q
v
v
v
a
a
(2.62)
or, in extended form,
.
^
)
(
^
^
)
(
^
0
0
^
0
0
0
0
0
0
^
]
^
^
[
^
1
1
1
1
1
»
¼
º
«
¬
ª
»
»
»
¼
º
«
«
«
¬
ª
Q
v
v
a
a
a
a
a
a
t
t
m
m
m
#
%
(2.63)
The eigenvectors of
^
can now be selected such that the columns of the
third matrix on the right hand side of (2.63) are orthogonal.
35
From linear algebra, the columns of
1
^
^
n
diag
a
, will be orthogonal if and
only if
.
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
^
)
(
^
)
(
^
^
)
(
^
)
(
^
^
1
1
1
1
1
1
1
1
»
»
»
»
»
»
¼
º
«
«
«
«
«
«
¬
ª
»
»
¼
º
«
«
¬
ª
¸¸¹
·
¨¨©
§
»
»
¼
º
«
«
¬
ª
¸¸¹
·
¨¨©
§
m
m
m
T
m
m
T
m
m
m
T
T
n
n
T
n
n
t
t
diag
t
t
diag
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
#
%
#
(2.64)
Inspection of (2.64) reveals that each off-diagonal term can be expressed in
the form
°
°
°
°
°
¯
°
°
°
°
°
®
z
¸¸¹
·
¨¨©
§
¦
¦
¦
¦
¦
.
)
^
sin(
^
)
^
sin(
^
)
^
sin(
^
)
^
sin(
^
,
)
^
sin(
^
)
^
sin(
^
^
^
^
^
1
1
1
2
1
1
2
k
n
for
t
t
t
t
k
n
for
t
t
N
k
j
k
k
N
j
n
j
n
n
k
j
k
k
N
n
j
n
n
N
j
n
j
n
n
N
j
n
j
n
n
k
n
k
T
n
j
j
A
A
A
A
A
A
a
a
a
a
(2.65)
Thus, if the natural frequencies
n
Z
are distinct, it then follows that if the
columns of
1
^
^
n
diag
a
are orthogonal if we consider an infinite set of sampled
values. In particular, as
N
increases the off-diagonal terms vanish, except for
the diagonal terms.
Formally, we may write
.
,
,
0
^
^
^
^
k
n
to
N
f
k
k
n
n
z
f
o
o
i
a
a
a
a
(2.66)
Therefore, it can be concluded that POMs converge to the modal vectors of
an undamped and unforced linear system for
N
sufficiently large.
Once the columns of
1
^
^
n
diag
a
are made orthogonal, the columns of Q
can be computed in order that they are orthogonal to those of
1
^
^
n
diag
a
. The
analysis of more complicated situations follows along the same lines.
36
The following observations can be made about the singular value
decomposition of the response matrix
^
.
x The columns of
U
represent the eigenmodes of the system.
x The first
n
columns of
V
are normalized time modulation of the modes.
x For a finite length time series, the columns of
U
are modal vectors.
On this basis, a comprehensive numerical approach to extract dynamic
trends from measured data is proposed in Chapter 3.
2.8 References
[1] P. Holmes, J. L. Lumley, G. Berkooz, "Turbulence, coherent structures, dynamical
systems and symmetry," New York: Cambridge University Press, 1996.
[2] G. Kerschen, J. C. Galinval, Alexander F. Vakakis and Lawrence A. Bergman, "The
method of proper orthogonal decomposition for dynamical characterization and
order reduction of mechanical systems: An overview," Nonlinear Dynamics, vol. 41,
pp. 147-169, 2005.
[3] S. S. Ravindran, "Reduced-order controllers for control of flow past an airfoil,"
International Journal for Numerical Methods in Fluid
, vol. 50, no. 5, pp. 531-554, 2006.
[4] A. R. Messina, V. Vittal, "Extraction of dynamics pattern from wide area
measurements using empirical orthogonal functions," IEEE Transactions on Power
Systems,
vol. 22, no. 2, pp. 682-692, 2007.
[5] G. Kerschen, J.-C. Golinval, A. F. Vakakis, and L. A. Bergman, "The method of
proper orthogonal decomposition for dynamical characterization and order
reduction of mechanical systems: An overview," Nonlinear Dynamics, vol. 41, pp.
147-169, 2005.
[6] L. Sirovich, "Turbulence and the dynamics of coherent structures part 1: coherent
structure," Quarterly of Applied Mathematics, vol. 95, no. 3, pp. 561-571, Oct., 1986.
[7] D. Kalman, "A singularly valuable decomposition: the SVD of a matrix," The
College Mathematics Journal
, vol. 27, no. 1, pp. 2-23, Jan., 1996.
[8] G. Kerschen, and J. C. Golinval, "Physical interpretation of the proper orthogonal
modes using the singular value decomposition," Journal of Sound and Vibration, vol.
249, no. 5, pp. 849-865, 2002.
[9] J. Argyris, and H. P. Mlejner, "Dynamical of structures," Text on Computational
Dynamics
, North Holland, Amsterdam, 1991.
37
[10] S. Han, "Application of proper orthogonal decomposition to structural vibration
analysis," Mechanical systems and signal processing, vol. 17, no. 5, pp. 989-1001, 2003.
[11] B. F. Freeny, and R. Kappagantu, "On the physical interpretation of proper
orthogonal modes in vibrations," Journal of Sound and Vibration, vol. 211, no. 4, pp.
607-616, 1998.
38
Chapter 3
Complex Empirical Orthogonal
Function
Many power system processes involve variability over both space and time.
These processes may contain moving patterns, travelling waves of different
spatial scales and temporal frequencies. The identification of dynamic
structures can be time-consuming and difficult because of the large amount of
information that needs to be processed and compared.
The conventional empirical orthogonal eigenfunction analysis is primarily a
method of compressing time and space variability of a data set into the lowest
possible number of spatial patterns. Each of these patterns is composed of
standing modes of variability and modulated by a time function. As such, real-
EOF analysis detects only standing waves, not travelling waves. The key point
to observe is that real-EOF analysis cannot deal with propagating features and
it only uses spatial correlation of the field, it is necessary to use both spatial and
time information in order to identify such features.
In this chapter, we extend the conventional empirical orthogonal function
analysis to the study and detection of propagating features in nonlinear
processes. A new modelling technique that accurately describes the
fundamental dynamic of large data sets through the use of reduced-order
modelling is proposed.
This approach is ideally suited to the study of propagating and standing
features that can be associated with observed oscillations and is effective at
identifying and isolating nonstationary behavior.
39
3.1 Complex EOF Analysis
Empirical orthogonal function analysis of data fields is commonly carried out
under the assumption that each field can be represented as a spatially fixed
pattern of behavior. This method, however, cannot be used for detection of
propagating features because of the lack of phase information. To fully utilize
the data, a technique is needed that acknowledges the nonstationarity of the
time-series data [1-11].
Let
)
,
(
k
j
t
x
u
be a space-time scalar field representing a time series, where
j
x
,
n
j
,
,
1
, is a set of spatial variables on a space
k
:
, and
k
t
,
N
k
,
,
1
is the
time at which the observations are made.
Provided
u
is simple and square integrable, it has a Fourier representation
of the form [1]
¦
f
1
)
(
)
(
)]
sin(
)
(
)
cos(
)
(
[
)
,
(
m
m
j
m
j
j
t
m
b
t
m
a
t
x
u
(3.1)
where
)
(
)
(
Z
m
j
a
and
)
(
)
(
Z
m
j
b
are the Fourier coefficients defined as
³
³
)
(
)
(
sin
)
,
(
1
cos
)
,
(
1
d
t
m
t
x
u
b
d
t
m
t
x
u
a
j
m
j
j
m
j
This allows description of travelling waves propagating throughout the
system. Equation (3.1) can be rewritten in the form
¦
f
1
)
(
)
(
)
,
(
m
t
im
m
j
j
c
e
c
t
x
u
(3.2)
where
)
(
)
(
)
(
)
(
)
(
)
(
Z
Z
Z
m
j
m
j
m
j
ib
a
c
[1] and
1
i
is the unit complex number.
Expanding (3.2) and collecting terms gives
40
^
`
^
`
)
,
(
)
,
(
)]
sin(
)
(
)
cos(
)
(
[
)]
sin(
)
(
)
cos(
)
(
[
)
,
(
1
)
(
)
(
1
)
(
)
(
t
x
u
i
t
x
u
t
m
a
t
m
b
i
t
m
b
t
m
a
t
x
u
j
H
j
m
m
j
m
j
m
m
j
m
j
j
c
¦
¦
f
f
(3.3)
where the real part of
c
u
is given by (3.1) and the imaginary part is the Hilbert
transform of
)
,
(
t
x
u
j
.
In formal terms, the Hilbert transform of a continuous time series
)
,
(
t
x
u
j
is
defined by the convolution [12]
³
f
f
dy
y
t
y
u
t
x
u
j
H
)
(
1
)
,
(
, (3.4)
where the integral is taken to mean the Cauchy principal value.
The most well-known classical methods for computing the Hilbert transform
are derived from the Fourier transform. However, this transform has a global
character and hence, it is not suitable for characterization of local signal
parameters. Alternatives for local implementation of the Hilbert transformation
which are based on local properties are developed and tested in this analysis.
For discretely sampled data, the Hilbert transform can be derived, in time
domain, by applying a rectangular rule to (3.4). It can be shown that
¦
f
f
|
k
j
H
k
k
t
u
t
x
u
1
2
)
)
1
2
(
(
2
)
,
(
W
S
, (3.5)
where
W
is the step size. When (3.5) is applied to a discrete time series
)
,
(
k
j
t
x
u
,
,
1
,
0
r
k
, one gets [8]
¦
¦
t
f
f
0
)]
1
2
(
)
1
2
(
[
1
2
1
2
1
2
)
1
2
(
2
)
,
(
k
k
k
j
H
k
t
u
k
t
u
k
k
k
t
u
t
x
u
. (3.6)
41
In previous formulations, the Hilbert transform was estimated by truncating
the series (3.6) [11]. This truncation was approximated using a convolution filter
as
¦
f
L
L
k
j
k
j
H
L
h
t
x
u
t
x
u
A
A
A
,
)
(
)
,
(
)
,
(
(3.7)
where
h
is a convolution filter with unit amplitude response and 90° phase
shift. In this research, it has been found that a simple filter that has the desired
properties of approximate unit amplitude response and /2 phase shift is given
by [13]
°¯
°
®
z
0
,
0
0
,
)
2
(
sin
2
)
(
2
l
l
l
l
l
h
(3.8)
where
L
l
L
d
d
. We omit the calculations.
As
f
o
L
, equation (3.7) yields an exact Hilbert transform. This represents a
filtering operation upon
)
,
(
t
x
u
j
in which the amplitude of each Fourier spectral
component remains unchanged while its phase is advanced by /2. Hannachi et
al. [8] found that
25
7
d
d L
provides adequate values for the filter response.
In what follows, we discuss the extension of the conventional EOF analysis
using the above approach to compute standing and propagating features of
general nonstationary processes where the eigenvectors of the covariance
matrix are complex and can be expressed alternatively as a magnitude and
phase pair.
3.2 Theoretical fundamentals of Empirical Orthogonal
Functions
The relationship between spatial and temporal behavior can be obtained by
noting that spatio-temporal information can be mapped into a time and space
grid. We suppose that we have a gridded data set composed of a space-time
field,
)
,
( t
x
u
representing the value of the field u , at time
t
, and spatial position
42
x . Figure 3.1 shows a conceptual representation of spatial and temporal
variability captured in the energy field. Rows in this plot capture spatial
information, whilst columns capture temporal information.
The method of complex EOF analysis is an optimal technique of
biorthogonal decomposition to find a spatial and temporal basis that spans an
ensemble of data collected from experiments or numerical simulations [1-3].
The method essentially decomposes a fluctuating field into a weighted linear
sum of spatial orthogonal modes and temporal orthogonal modes such that the
projection onto the first few modes is optimal. Before introducing the method,
the theoretical background for the biorthogonal decomposition of spatial and
temporal variability is introduced.
The fundamental idea of EOF analysis is to find a basis
M for linear, infinite-
dimensional Hilbert space
])
1
,
0
([
2
L
, that maximizes the averaged projection of a
response matrix, X , suitably normalized; the corresponding functional for this
constrained variational problem is solved and reduced in section 2.2 to
0
)
(
)
(
'
)
'
(
)
'
(
)
(
*
1
0
1
0
*
»
¼
º
«
¬
ª
³ ³
dx
x
x
dx
x
x
u
x
u
\
OM
M
(3.9)
where, if
0
)
(
*
x
\
, the optimal basis is given by the eigenfunctions
j
M
of the
integral equation
)
(
'
)
'
(
)
'
(
)
(
1
0
*
x
dx
x
x
u
x
u
OM
M
³
whose kernel is the averaged autocorrelation function
)
'
,
(
)
(
)
(
*
x
x
x
u
x
u
def
C
[4-
8]. They are consequently called empirical orthogonal functions (EOFs). Now,
we assume that if
0
)
(
*
z
x
\
, then (3.9) can be rewritten as
dx
x
x
dx
dx
x
x
u
x
u
x
)
(
)
(
'
)
(
)
'
(
)
(
)
'
(
*
1
0
*
1
0
1
0
*
\
O
M
\
M
³
³³
(3.10)
such that the inner product
0
))
(
)
'
(
(
*
z
x
x
\
M
C
with orthogonal eigenvectors,
\
M, , i.e.,
43
¯
®
z
¯
®
z
j
i
j
i
j
i
j
i
ij
j
T
i
ij
j
T
i
,
0
),
(
and
,
0
),
(
\
G
\
\
M
G
M
M
From (3.10) it can be seen that if there exists an arbitrary variation (spatial),
0
)
(
*
z
x
\
, then the original field can be reconstructed using two optimal
orthogonal basis,
\
M, , given from (3.10).
a) Three-dimensional view of energy distribution.
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
u [x
1,
x
2 ,
x
3,
...
x
n
]
t
1
t
N
.
.
.
Spatial
Energy
Temporal
Energy
...
...
...
...
...
...
...
...
...
...
...
...
...
b) Two-dimensional representation.
Fig. 3. 1. Conceptual diagram illustrating the phenomenon of temporal and spatial energy in a
field.
44
Based on this notion, an efficient technique to find the optimal orthogonal
basis using complex EOF analysis is developed.
3.3 Complex-Vector EOF Analysis
Drawing on the above approach, an efficient formulation to compute complex
time-dependent POMs has been derived. Following Feeny et al. [3], assume that
I
R
c
iX
X
X
is a n (spatial points) by N (temporal points) ensemble complex
matrix. From the preceding results, it follows that
H
I
I
I
H
R
R
R
c
c
i
t
isin
t
cos
c
c
V
U
V
U
X
X
¦
¦
)]
(
)
(
[
(3.11)
where
c
and
c
X
T are the magnitude and phase of
c
,
H
R
V
and
H
I
V
are the
conjugate transpose of
R
V
and
I
V
respectively, the superscript H denotes a
Hermitian matrix and the subscripts
I
R
c
,
,
, indicate the complex, real and
imaginary vector.
Implicit in this model is the assumption that
I
V
V
I
V
V
I
U
U
I
U
U
I
H
I
R
H
R
I
H
I
R
H
R
,
,
From equation (2.17), the autocorrelation matrix
C
for the case complex data
can be written in the form
I
R
iC
C
C
. It can be easily verified that
)]}
(
)
(
)
(
)
(
[
)]
(
)
(
)
(
)
(
{[
)]
(
)
(
)][
(
)
(
[
)]}
(
)
(
)
(
)
(
[
)]
(
)
(
)
(
)
(
{[
)]
(
)
(
)][
(
)
(
[
t
cos
t
sin
t
sin
t
cos
i
t
sin
t
sin
t
cos
t
cos
t
isin
t
cos
t
isin
t
cos
t
sin
t
cos
t
cos
t
sin
i
t
sin
t
sin
t
cos
t
cos
t
isin
t
cos
t
isin
t
cos
c
T
c
c
T
c
c
T
c
c
T
c
c
c
T
c
T
c
T
c
c
T
c
c
T
c
c
T
c
c
T
c
T
c
c
c
c
T
c
c
T
c
c
H
c
T
c
c
T
c
c
H
c
c
X
X
X
X
X
X
X
X
X
X
X
X
(3.12)
which translates into
45
]
[
]
[
]
[
]
[
]
[
]
[
]
[
]
[
R
R
T
I
T
I
I
I
T
R
T
R
I
I
T
I
T
I
R
R
T
R
T
R
R
T
I
I
T
R
I
T
I
R
T
R
c
H
c
T
I
T
I
R
R
T
R
T
R
I
I
T
I
T
I
I
I
T
R
T
R
R
R
T
I
R
T
R
I
T
I
I
T
R
R
H
c
c
i
i
i
i
V
V
V
V
V
V
V
V
U
U
U
U
U
U
U
U
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
¦
(3.13)
where,
T
R
and
T
I
denotes the transpose of
R
and
I
, respectively.
As is apparent from (3.11), the columns of
R
U
and
I
U
are the eigenvectors
of the real and imaginary parts of
H
c
c
, and the columns of
R
V
and
I
V
are the
eigenvectors of real and imaginary parts of
c
H
c
, respectively. The n singular
values on the diagonal of
R
and
I
are the square roots of the nonzero
eigenvalues of the real and imaginary parts of both
H
c
c
and
c
H
c
divided
by the number of samples
N
, where n is the rank of the complex matrix
c
X
.
It can be seen from (3.12)-(3.13) that the autocovariance matrix
I
R
iC
C
C
,
where the real part is a symmetrical matrix, that is to say,
7
R
R
C
C
and
I
C
is a
skew-symmetric matrix i.e. say,
I
I
C
C
7
. If the size of
I
C
is odd then the
determinant of
I
C
will always be zero. Because the symmetrical matrix is a
particular case of the Hermitian matrix, then all its eigenvectors are real.
Further, the eigenvalues of the skew-symmetric matrix are all imaginary pure
and is a normal matrix; its eigenvectors are complex conjugates.
From the decomposition given in (3.12) can be seen that the imaginary part
is zero when time is in phase with the extremum of the cosine or sine, that is,
the sum of the two components is zero; at this time instant both are symmetrical
matrices. The imaginary part measures the degree of asymmetry when the sum
of both matrices is different from zero; this is used to define the existence of
arbitrary variations into the space,
0
)
(
*
z
x
\
; this feature is used to define the
existence of travelling wave components.
We also remark from the eigenvalue decomposition in (3.13) that the basis of
the empirical orthogonal functions are defined as
)
(x
R
to the real part of the
field complex, and
)
(x
I
as the imaginary part; the approximation to the
original field is written as the sum of pure standing and travelling wave
components denoted by the subscripts
swc
and
twc
, as
46
twc
swc
c
X
X
X
(3.14)
In an attempt to address these issues in the context of EOF analysis, an
analytical framework that treats both standing and travelling features is next
introduced.
3.4 Analysis of Propagating Features
The time-dependent complex coefficients associated with each eigenfunction
can be conveniently split into their amplitudes and phases.
Following the notation defined earlier, the ensemble of data can be
expressed as the complex expansion
)
(
)
(
)
(
)
(
^
1
)
(
)
(
1
)
(
)
(
x
t
i
x
t
q
j
j
I
j
I
p
j
j
R
j
R
c
¦
¦
+
+
A
A
X
(3.15)
where
)
(
)
(
t
j
R
A
and
)
(
)
(
t
j
I
A
are obtained as the projection of the basis
)
(
)
(
x
j
R
and
)
(
)
(
x
j
I
respectively into complex field
c
X
of the form
)
(
)
(
)
(
)
(
)
(
)
(
x
x
j
I
c
j
I
j
R
c
j
R
X
A
X
A
(3.16)
The time-dependent complex coefficients associated with each eigenfunction
can be conveniently split into their amplitudes and phases. From complex EOF
analysis in (3.15), the ensemble of data can be expressed as [5,12-15]
)
2
(
)
(
)
2
(
)
(
)
(
)
(
)
(
)
(
)
(
1
)
(
)
(
)
(
)
(
1
)
(
x
t
x
t
j
I
j
I
j
I
q
j
j
I
j
R
j
R
j
R
p
j
j
R
c
¦
¦
I
I
S
R
S
R
X
,(3.17)
where
)
(
)
(
t
j
R
R
,
)
(
)
(
t
j
I
R
are the temporal amplitude function, and
)
(
)
(
x
j
R
S
, and
)
(
)
(
x
j
I
S
are the spatial amplitude functions associated with the wave
decomposition. This effectively decomposes the data into spatial and temporal
modes.
47
Four measures that define moving features in
)
,
( t
x
u
can then be defined:
1. Spatial distribution of variability of each eigenmode
2. Relative phase fluctuation
3. Temporal variability in magnitude
4. Variability of the phase of a particular oscillation
The succeeding sections describe the properties of these representations.
3.4.1 Spatial Amplitude Function,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
S
S
The spatial amplitude functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
S
S
, are defined as
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
t
t
t
t
t
t
j
I
j
I
j
I
j
R
j
R
j
R
S
S
+
+
(3.18)
This function shows the spatial distribution of variability associated with
each eigenmode.
3.4.2 Spatial Phase Functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
I
I
This function shows the relative phase fluctuation among the various spatial
locations where u is defined and is given by
°¿
°
¾
½
°¯
°
®
°¿
°
¾
½
°¯
°
®
)
(
[
)]
(
[
tan
)
(
)
(
[
)]
(
[
tan
)
(
)
(
)
(
1
)
(
)
(
)
(
1
)
(
t
re
t
im
t
t
re
t
im
t
j
I
j
I
j
I
j
R
j
R
j
R
I
I
(3.19)
This measure, for which an arbitrary reference value must be selected, varies
continuously between 0 and
S
2 . Together these equations give a measure of
the space distribution of energy and can be used to identify the dominant
modes and their phase relationships. Further, for each dominant mode of
interest, a mode shape can be computed by using the spatial part of (3.15).
3.4.3 Temporal Amplitude Functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
R
R
Similar to the description of the spatial amplitude function, the temporal
amplitude functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
R
R
, can be defined as
48
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
t
t
t
t
t
t
j
I
j
I
j
I
j
R
j
R
j
R
A
A
R
A
A
R
+
+
(3.20)
These functions give a measure of the temporal variability in the magnitude
of the modal structure of the field,
u
. As described later, the general form of
these equations is very amenable to computational analysis.
3.4.4 Temporal Phase Functions,
)
(
),
(
)
(
)
(
t
t
j
I
j
R
This temporal variation of phase associated with
)
,
(
t
x
u
is given by
°¿
°
¾
½
°¯
°
®
°¿
°
¾
½
°¯
°
®
)
(
[
)]
(
[
tan
)
(
)
(
[
)]
(
[
tan
)
(
)
(
)
(
1
)
(
)
(
)
(
1
)
(
t
re
t
im
t
t
re
t
im
t
j
I
j
I
j
I
j
R
j
R
j
R
A
A
A
A
(3.21)
For the simple case of a sinusoidal wave with fixed frequency and wave
number,
t
t
i
Z
)
(
. In the more general (and interesting case), the derivative of
the phase and frequency of the modal components can be calculated from
i
i
i
i
i
i
i
k
dt
d
dx
d
k
c
I
(3.22)
where
i
c
is the phase velocity of the function.
Equations (3.15)-(3.22) provide a complete characterization of any
propagating effects and periodicity in the original data field which might be
obscured by standard cross-spectral analysis. Finally, it might be remarked that,
in the special case of real analysis, these expressions simplify to the normal
definitions. We now turn our attention to the analysis of spatio-temporal
behavior.
49
3.4.5 Spatio-Temporal Analysis
In order to investigate travelling features, let Equation (3.17) be expressed as
)
(
)
(
1
)
(
)
(
)
(
1
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
i
j
I
q
j
j
I
i
j
R
p
j
j
R
c
j
I
j
I
j
R
j
R
e
x
t
e
x
t
¦
¦
II
I
S
R
S
R
X
(3.23)
where the real physical field is reconstructed by taking the real part of the
before complex expansion,
)
cos(
)
(
)
(
)
cos(
)
(
)
(
)
(
)
(
)
(
1
)
(
)
(
)
(
)
(
1
)
(
x
t
x
t
j
I
j
I
j
I
q
j
j
I
j
R
j
R
j
R
p
j
j
R
¦
¦
II
I
S
R
S
R
X
(3.24)
This equations defines how the dynamical characterization can provide
information on the variability of the j -th mode into original field as [15-17]
)
cos(
)
(
)
(
)
cos(
)
(
)
(
)
(
)
(
)
(
1
)
(
)
(
)
(
)
(
1
)
(
x
t
x
t
x
t
x
t
j
I
j
I
j
I
q
j
j
I
j
R
j
R
j
R
p
j
j
R
¦
¦
S
R
S
R
X
(3.25)
In terms of the complex representation, the wave components are given by
[5]:
1. Wave number
)
( j
R
,
)
( j
I
. The wave numbers are given by:
dx
d
dx
d
j
I
j
I
j
R
j
R
)
(
0
)
(
)
(
)
(
)
(
)
(
I
I
(3.26)
where units are in
1
.
m
rad
, where
0
)
(
j
R
. Since
)
( j
R
I
takes values from 0 or
;
)
( j
R
I
is the angle of the base function given from the eigenvectors of the
matrix
R
C
obtained in (3.13). Given that
R
C
is a symmetrical matrix, all of its
eigenvalues are real; it can then be seen from (3.25) that in the first term, all
parts represent stationary wave modes, whereas in the second term, they
represent travelling wave modes. The wave number is only defined for
travelling waves; therefore (3.25) can be rewritten as
50
)
cos(
)
(
)
(
)
cos(
)
(
)
(
)
(
)
(
)
(
1
)
(
)
(
)
(
1
)
(
x
t
x
t
t
x
t
j
I
j
I
j
I
q
j
j
I
j
R
j
R
p
j
j
R
¦
¦
S
R
S
R
X
(3.27)
2. Angular frequencies
)
( j
R
,
)
( j
I
. Represent the angular frequency in
s
rad /
and are defined as
dt
d
dt
d
j
I
j
I
j
R
j
R
)
(
)
(
)
(
)
(
)
(
)
(
(3.28)
3. Average phase speeds
)
( j
R
c
,
)
( j
I
c
. They are computed from
)
(
)
(
)
(
)
(
)
(
)
(
0
j
I
j
I
j
I
j
R
j
R
j
R
c
c
(3.29)
physical units are
s
m / .
These equations provide a complete characterization of any propagating
features and periodicity in the original data field.
From (3.17) and (3.27) it can be seen that the term associated with the j -th
travelling wave can be expressed in the form
)]}
2
/
(
)
2
/
{cos[(
)
cos(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
x
t
x
t
j
I
j
I
j
I
j
I
j
I
j
I
j
I
j
I
S
R
S
R
(3.30)
Using the trigonometric identities
)
sin(
)
sin(
)
cos(
)
cos(
)
cos(
B
A
B
A
B
A
and
)
sin(
)
cos(
)
cos(
)
sin(
)
sin(
B
A
B
A
B
A
, where
)
2
/
(
)
(
S
t
A
j
I
and
)]
2
/
(
)
(
x
B
j
I
, and noting that
)
cos(
)
2
/
sin(
)
cos(
)
2
/
sin(
)
sin(
)
2
/
cos(
)
sin(
)
2
/
cos(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
x
x
t
t
x
x
t
t
j
I
j
I
j
I
j
I
j
I
j
I
j
I
j
I
S
S
then (3.30) can be re-expressed in the form
51
)}
cos(
)
cos(
)
sin(
)
{sin(
)
cos(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
x
t
x
t
x
t
j
I
j
I
j
I
j
I
j
I
j
I
j
I
j
I
j
I
j
I
S
R
S
R
.(3.31)
From (3.31), it can be seen that the travelling wave components are also
identified as the sum of two standing wave components. To obtain the
decomposition of the original field in its pure standing wave components, it is
necessary to compute the difference with the pure travelling wave components
of the form
)
cos(
)
(
)
(
)
(
)
(
1
)
(
t
x
t
j
R
j
R
p
j
j
R
twc
swc
S
R
X
-
X
X
¦
(3.32)
where
swc
X
represents the decomposition of the original field given by the pure
standing wave components, and
)
cos(
)
(
)
(
)
(
)
(
)
(
1
)
(
x
t
x
t
j
I
j
I
j
I
q
j
j
I
twc
¦
S
R
X
, (3.33)
represents the decomposition of the original field given by the pure travelling
wave components, where the damping of each mode is given by its amplitude.
A practical criterion for choosing the number of relevant modes is given by
the energy percentage contained in the p and q modes of the form
%
99
100
2
1
)
(
1
)
(
u
¦
¦
F
q
j
twc
j
p
j
swc
j
X
(3.34)
where
2
.
F
denotes the Frobenius norm.
Equations (3.27) through (3.34) provide a complete characterization of any
propagating effects and periodicity in the original data field which might be
obscured by normal cross-spectral analysis. Finally, it might be remarked that,
in the special case of real analysis, these expressions are simplified to the
normal definitions.
Later, we discuss refinements to the above model, so that travelling features
can be accurate characterized using real-time information.
52
3.4.6 Motivating Examples: Modeling of Propagating Waves Using
the Covariance Matrix
In an effort to better understand the mechanism of wave propagation using
complex EOF analysis we consider, as first example, a nondispersive plane
wave.
For simplicity, consider a nondispersive plane wave propagating at phase
speed
c
and wavenumber defined as
c
k
/
, past an array of sensors at
positions j given by [5]
¦
j
j
j
t
kx
t
kx
t
u
)]
sin(
)
(
)
cos(
)
(
[
)
(
(3.35)
Expanding (3.35) we have
¦
Z
Z
Z
E
Z
D
Z
Z
E
Z
D
)}
sin(
)]
cos(
)
(
)
sin(
)
(
[
)
cos(
)]
sin(
)
(
)
cos(
)
(
{[
)
(
t
kx
kx
t
kx
kx
t
u
j
j
j
j
j
(3.36)
Now, using the identities
)
cos(
)
(
)
sin(
)
(
)
(
)
sin(
)
(
)
cos(
)
(
)
(
j
j
j
j
j
j
kx
kx
b
kx
kx
a
(3.37)
(3.35) can be rewritten as
^
`
¦
j
j
j
t
b
t
a
t
u
)
sin(
)
(
)
cos(
)
(
)
(
(3.38)
Where, for this example,
u
is a white, band-limited signal given by
¯
®
!
d
d
$
1
2
2
1
2
2
2
,
,
0
,
)
(
)
(
if
if
(3.39)
To obtain phase information between stations, a complex representation of
(3.38) is invoked:
^
`
¦
t
i
j
j
j
e
ib
a
t
u
)]
(
)
(
[
)
(
(3.40)
The covariance matrix of (3.40) is given as
53
t
k
j
jk
t
u
t
u
)
(
)
(
C
*
(3.41)
where
t
denotes time averaging and the asterisk complex conjugation.
Substations (3.40) onto (3.41) and expanding gives
^
` ^
`
^
`
^
`
¦
¦
¦
¦
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
Z
)]
(
)
(
)
(
)
(
[
)]
(
)
(
)
(
)
(
[
,
)]
(
)
(
[
*
)]
(
)
(
[
)]
(
)
(
[
*
)]
(
)
(
[
C
k
j
k
j
k
j
k
j
k
k
j
j
t
i
k
k
t
i
j
j
jk
a
b
b
a
i
b
b
a
a
ib
a
ib
a
e
ib
a
e
ib
a
(3.42)
Using (3.37) and simplifying, (3.42) reduces to
¦
Z
Z
E
Z
D
Z
E
Z
D
)]}
sin(
)
(
)
sin(
)
(
[
)]
cos(
)
(
)
cos(
)
(
{[
C
2
2
2
2
j
k
j
k
k
j
k
j
jk
kx
kx
kx
kx
i
kx
kx
kx
kx
(3.43)
For calculation purposes it is better to use the exponential form. Using the
condition give in (3.39),
C
can be rewritten as
^
`
¦
¦
'
$
$
2
1
)
(
2
2
)
sin(
)
cos(
C
x
ik
k
j
k
j
jk
jk
e
kx
kx
i
kx
kx
(3.44)
where
k
j
jk
x
x
x
'
.
Replacing the summation with an integral in (3.44) and multiplying by
'
'
, yields
7
'
'
$
³
'
d
e
x
ik
jk
jk
2
to
,
C
2
1
)
(
2
(3.45)
This equation is amenable to integration by parts. Defining
)
(
)
(
)]
(
[
)
(
,
)
(
d
d
k
d
x
i
d
x
k
i
jk
jk
¿
¾
½
¯
®
'
'
(3.46)
where
)
(
d
is given by
54
)
(
)]
(
[
)
(
)
(
)]
(
[
)
(
)
(
d
k
d
x
i
d
d
k
d
x
i
d
d
jk
jk
¸
¸
¹
·
¨
¨
©
§
'
'
(3.47)
and substituting (3.46) into (3.45), one has
)
(
)]
(
[
)
(
2
C
2
1
2
d
e
k
d
x
i
d
jk
jk
³
¸
¸
¹
·
¨
¨
©
§
'
7
$
(3.48)
Integrating (3.48) by parts, and using the above results we obtain
@
2
1
)]
(
[
)
(
2
C
2
jk
jk
e
k
d
x
i
d
¸
¸
¹
·
¨
¨
©
§
'
7
$
(3.49)
which gives
@
@
jk
jk
jk
x
ik
x
ik
jk
x
ik
jk
jk
e
e
k
d
x
i
d
e
k
d
x
i
d
'
'
'
¸
¸
¹
·
¨
¨
©
§
'
7
$
¸
¸
¹
·
¨
¨
©
§
'
7
$
)
(
)
(
2
)
(
2
2
1
2
1
)]
(
[
)
(
2
,
)]
(
[
)
(
2
C
(3.50)
After some algebra, we can show that
»
»
»
¼
º
«
«
«
¬
ª
¸
¸
¹
·
¨
¨
©
§
'
7
$
'
'
'
'
'
'
'
'
2
)
(
2
)
(
2
)
(
2
)
(
2
)
(
2
)
(
2
)
(
2
)
(
2
1
1
2
2
2
2
1
1
)]
(
[
)
(
2
C
jk
jk
jk
jk
jk
jk
jk
jk
x
k
i
x
k
i
x
k
i
x
k
i
x
k
i
x
k
i
x
k
i
x
k
i
jk
jk
e
e
e
e
e
e
e
k
d
x
i
d
(3.51)
¸
¸
¹
·
¨
¨
©
§
¸
¸
¹
·
¨
¨
©
§
¸
¸
¹
·
¨
¨
©
§
'
7
$
'
'
'
'
'
'
2
)
(
2
)
(
2
)
(
2
)
(
2
)
(
2
)
(
2
2
1
2
1
2
1
)]
(
[
)
(
2
C
jk
jk
jk
jk
jk
jk
x
k
i
x
k
i
x
k
i
x
k
i
x
k
i
x
k
i
jk
jk
e
e
e
e
e
e
k
d
x
i
d
(3.52)
or
2
)]
(
)
(
[
2
)]
(
)
(
[
2
)]
(
)
(
[
2
1
2
1
2
1
2
)]
(
[
)
(
2
C
jk
jk
jk
x
k
k
i
x
k
k
i
x
k
k
i
jk
jk
e
e
e
k
d
x
i
d
'
'
'
¸
¸
¹
·
¨
¨
©
§
¸
¸
¹
·
¨
¨
©
§
'
7
$
(3.53)
55
2
)]
(
)
(
[
2
)]
(
)
(
[
2
)]
(
)
(
[
2
1
2
1
2
1
2
2
)]
(
[
)
(
C
jk
jk
jk
x
k
k
i
x
k
k
i
x
k
k
i
jk
jk
e
i
e
e
k
d
x
d
'
'
'
¸
¸
¸
¹
·
¨
¨
¨
©
§
¸
¸
¹
·
¨
¨
©
§
'
7
$
(3.54)
Finally, using the complex identity
i
e
e
i
i
2
)
(
)
sin(
in (3.54), we
obtain
2
)]
(
)
(
[
1
2
1
2
1
2
2
2
)]
(
)
(
[
1
2
2
1
2
1
2
2
)]
(
)
(
[
2
)]
(
)
(
[
sin
2
]
[
,
2
)]
(
)
(
[
sin
)]
(
[
)
(
C
jk
jk
x
k
k
i
jk
jk
x
k
k
i
jk
jk
jk
e
x
k
k
x
k
k
e
x
k
k
k
d
x
d
'
'
'
¸¸¹
·
¨¨©
§
'
7
$
¸¸¹
·
¨¨©
§
'
¸
¸
¹
·
¨
¨
©
§
'
7
$
(3.55)
or, more simply,
jk
x
k
i
jk
jk
e
x
k
c
'
¸¸¹
·
¨¨©
§ '
'
0
$
2
sin
C
2
(3.56)
where k
' is the wavenumber bandwidth, x
' is the array length, and
1
)
0
(
sin
,
0
)
sin(
)
(
sin
,
2
)
(
,
2
)
(
)
(
,
)
(
)
(
1
2
1
2
1
2
7
0
'
c
y
if
y
y
y
c
k
k
k
k
k
k
with
)
(
sin
)
(
sin
y
c
y
c
.
General algebraic expression for the eigenvalues and eigenfunctions of
C
for an arbitrary number of sensors
)
(n
are very difficult to determine and are
not pursued here.
For the case of two sensors, the above model can be reduced to
56
2
,
1
,
,
2
)
(
sin
2
)
(
sin
C
2
)
(
1
2
2
)
(
2
1
2
2
1
2
2
1
»
»
»
»
¼
º
«
«
«
«
¬
ª
0
$
¸
¹
·
¨
©
§
'
0
$
¸
¹
·
¨
©
§
'
0
$
0
$
k
j
e
x
x
k
c
e
x
x
k
c
x
x
k
i
x
x
k
i
jk
(3.57)
Equation (3.57) illustrates some important properties of
C
. From the relation
above it follows that the eigenvalues of
C
are given by
]
C
det[
,k
j
,
, i.e.,
»
»
»
»
¼
º
«
«
«
«
¬
ª
0
$
¸
¹
·
¨
©
§
'
0
$
¸
¹
·
¨
©
§
'
0
$
0
$
2
)
(
1
2
2
)
(
2
1
2
2
1
2
2
1
2
)
(
sin
2
)
(
sin
det
e
x
x
k
c
e
x
x
k
c
x
x
k
i
x
x
k
i
(3.58)
2
2
2
2
2
2
)
(
2
1
2
)
(
2
1
2
2
2
2
2
2
2
sin
)
(
2
2
)
(
sin
2
)
(
sin
)
(
det
2
1
2
1
¸¸¹
·
¨¨©
§
¸
¹
·
¨
©
§ '
'
0
$
0
$
0
$
¸¸¹
·
¨¨©
§
¸
¹
·
¨
©
§
'
0
$
¸¸¹
·
¨¨©
§
¸
¹
·
¨
©
§
'
0
$
0
$
0
$
0
$
x
k
c
e
x
x
k
c
e
x
x
k
c
x
x
k
i
x
x
k
i
(3.59)
It is easy to see further that
2
2
sin
)
(
4
)
2
(
2
2
2
2
2
2
2
2
2
,
1
»
»
¼
º
«
«
¬
ª
¸¸¹
·
¨¨©
§
¸
¹
·
¨
©
§ '
'
0
$
0
$
0
$
r
0
$
x
k
c
(3.60)
and also,
»
¼
º
«
¬
ª
¸
¹
·
¨
©
§ '
'
r
0
$
¸
¹
·
¨
©
§ '
'
0
$
r
0
$
¸¸¹
·
¨¨©
§
¸
¹
·
¨
©
§ '
'
0
$
0
$
0
$
r
0
$
2
sin
1
2
sin
2
2
sin
4
)
(
4
)
(
4
2
2
2
2
2
2
2
2
2
2
2
2
,
1
x
k
c
x
k
c
x
k
c
O
(3.61)
It then follows that the eigenvectors of
0
]
C
[
,
%
,
k
j
, are given by
57
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
»
»
»
¼
º
«
«
«
«
¬
ª
0
$
¸
¹
·
¨
©
§
'
0
$
¸
¹
·
¨
©
§
'
0
$
0
$
0
0
2
)
(
sin
2
)
(
sin
1
,
2
1
,
1
2
1
)
(
1
2
2
)
(
2
1
2
2
1
1
2
2
1
e
x
x
k
c
e
x
x
k
c
x
x
k
i
x
x
k
i
(3.62)
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
»
»
»
¼
º
«
«
«
«
¬
ª
0
$
»
¼
º
«
¬
ª
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
0
$
»
¼
º
«
¬
ª
¸
¹
·
¨
©
§ '
'
0
$
'
'
0
0
2
sin
1
2
sin
2
sin
2
sin
1
1
,
2
1
,
1
2
2
2
2
2
2
x
k
c
e
x
k
c
e
x
k
c
x
k
c
x
k
i
x
k
i
(3.63)
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
»
»
»
¼
º
«
«
«
«
¬
ª
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
'
'
0
0
2
sin
2
sin
2
sin
2
sin
1
,
2
1
,
1
2
2
2
2
x
k
c
e
x
k
c
e
x
k
c
x
k
c
x
k
i
x
k
i
(3.64)
which can be simplified to
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
¼
º
«
¬
ª
'
0
0
0
0
1
1
,
2
1
,
1
x
k
i
e
(3.65)
where
x
k
i
x
k
i
e
if
e
'
'
7
%
7
%
%
%
1
,
1
1
,
2
1
,
2
1
,
1
,
,
(3.66)
2
2
1
with
,
1
,
2
1
,
1
x
k
i
x
k
i
x
k
i
e
e
e
7
7
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
(3.67)
Finally, solving for
2
,
1
% and
2
,
2
% , yields
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
»
»
»
¼
º
«
«
«
«
¬
ª
0
$
¸
¹
·
¨
©
§
'
0
$
¸
¹
·
¨
©
§
'
0
$
0
$
0
0
2
)
(
sin
2
)
(
sin
2
,
2
2
,
1
2
2
)
(
1
2
2
)
(
2
1
2
2
2
1
2
2
1
e
x
x
k
c
e
x
x
k
c
x
x
k
i
x
x
k
i
(3.68)
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
»
»
»
¼
º
«
«
«
«
¬
ª
0
$
»
¼
º
«
¬
ª
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
0
$
»
¼
º
«
¬
ª
¸
¹
·
¨
©
§ '
'
0
$
'
'
0
0
2
sin
1
2
sin
2
sin
2
sin
1
2
,
2
2
,
1
2
2
2
2
2
2
x
k
c
e
x
k
c
e
x
k
c
x
k
c
x
k
i
x
k
i
(3.69)
58
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
»
»
»
¼
º
«
«
«
«
¬
ª
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
¸
¹
·
¨
©
§ '
'
0
$
'
'
0
0
2
sin
2
sin
2
sin
2
sin
2
,
2
2
,
1
2
2
2
2
x
k
c
e
x
k
c
e
x
k
c
x
k
c
x
k
i
x
k
i
(3.70)
and
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
»
¼
º
«
¬
ª
'
0
0
0
0
1
2
,
2
2
,
1
x
k
i
e
(3.71)
where, in the above,
x
k
i
x
k
i
e
e
'
'
7
%
7
%
%
%
2
,
1
2
,
2
2
,
2
2
,
1
,
if
,
(3.72)
2
2
1
with
,
2
,
2
2
,
1
x
k
i
x
k
i
x
k
i
e
e
e
7
7
»
¼
º
«
¬
ª
»
¼
º
«
¬
ª
%
%
(3.73)
The following observations can be made from this analysis:
x As can be seen from (3.67) and (3.73), the method yields complex conjugate
eigenvectors and an average wave number, k .
x The value
j
x
k
gives the mode shape; this can be used for detection of wave
propagating into field original. We remark that the performance of the
complex EOF analysis, as measured by the percentage of variance given in
(3.61), depends on the spread in wave number relative to the array size; as
the parameter
x
k
'
'
decreases, more of the variance is contained in the
lowest complex EOF mode.
The development given in this section indicates that modal spatial patterns
from a time domain complex EOF analysis may be computed in a
straightforward manner.
A point of particular interest is that, as a standard technique for describing
coherent variability in spatial data, a relatively wide wave number bandwidth
)]
2
(
0
[
S
!
'
' x
k
results from (3.56). The importance of the parameter
x
k
'
'
in
explaining the variability is analogous to the importance of phase wobble when
computing the coherence in a spectral band.
59
The more the phase varies over a band, the lower the coherence for that
band. In the expression for
C
from (3.56),
jk
x
k
' is the phase analog and
jk
x
k
'
'
is the phase spread between any two sensors. Thus, as
jk
x
k
'
'
increases,
jk
C
decreases. Note further that
X
k
'
'
is a measure of the largest range in phase
computed over the array (between the two furthest sensors) while
)
1
/(
'
'
N
X
k
represents the minimun phase range between adjacent sensors.
In the next subsections we discuss the practical computation of coherency
and mode shape identification in relation to the context of singular values and
eigenvalues decomposition from measurements data.
3.4.7 Coherency Identification from Measurements Data
The procedure adopted to identify dynamic coherency from a set of observed
data using singular value decomposition is summarized as follows:
1. Given a set of measurements, assume that the data matrix
is an n (state
variables) by
N
(instants of time) matrix.
2. Compute the principal component analysis of
based on the singular
value decomposition.
3. Solve the singular value decomposition
T
V
U
X
. Compute the vectors
U
and
T
V
W
.
4. Reconstruct the
matrix through the full principal component analysis
decomposition
W
.
5. Variation of
is truncated to the first terms. Being as this is captured by the
first principal components.
Here,
has the coordinates
k
n
n
n
t
t
t
,
2
,
1
,
,...,
,
, in a n -dimentional space which
are called scores plot. Similar components have similar
t
-coordinates.
Therefore, such groups form clusters. This information of the clusters, in turn,
can be used to identify the coherence groups in vast interconnected systems.
These can be represented graphically as discussed in our numerical application
of the method (Chapter 5).
60
3.4.8 Mode Shape Identification
Several practical procedures for computing modal information are possible.
Based on the above ideas we proposed a simple algorithm to compute the mode
shape of measured data in the context of eigenvalues decomposition.
1. Perform a eigenvalues decomposition of the correlation matrix
C
T
N
1
.
2. The mode shapes are the columns of the matrix
U
. Sort the vectors
according to their relative sign.
Following a similar reasoning, it is now clear that the POMs, defined as the
eigenvectors of the correlation matrix or the right singular vectors of SVD,
determine the mode shape of a particular oscillation mode.
The above approach provides a clear picture of swing oscillation energy in
the network and allows us to identify important information of interest in the
application of power system stability analysis.
3.5 Numerical Computation of POMs: The Algorithm
In this section, a step-by-step description of the algorithm used to extract modal
information is presented. The procedure adopted to compute the POMs can be
summarized as follows [15]:
1. Given an ensemble of measurements of a nonstationary process, compute
the response matrix X . Form the complex time series matrix
H
jX
X
X
^
,
where
H
X
is the Hilbert transform of X
2. Compute the singular value decomposition
T
V
U
X
. Compute the
singular vectors
U
,
V
and the corresponding singular values
V
3. Determine the time evolution of the temporal modes,
i
a
. Extract standing
and propagating features using the complex SVD formulation.
Figure 3.2 illustrates the proposed algorithm.
61
Fig. 3. 2. Conceptual view of the proposed algorithm.
For different events recorded at the same location, statistical averaging can
be employed to take advantage of the statistics of the data. In this case, the
snapshots can be thought to be realizations of random fields generated by some
kind of stochastic process. The processing steps are detailed in the next chapter.
A numerical example is given now to visualize the use and application of
the above formulation.
3.5.1 Numerical Example
As an illustrative example, consider a field
)
,
( t
x
f
composed of standing and
travelling components denoted by
)
(
1 swc
f
,
)
(
2 twc
f
of the form
)
,
(
)
,
(
)
,
(
2
1
t
x
f
t
x
f
t
x
f
(3.74)
where
)
cos(
)
,
(
)
cos(
)
,
(
2
2
2
2
1
1
1
1
t
x
k
A
t
x
f
t
x
k
e
A
t
x
f
t
Z
Z
D
with constant wave number, angular frequency and amplitude.
The model is evaluated at 64 equally spaced points from
1
x
to
m
64 and
sampled at
s
t
1
'
intervals. Data are available at
s
t
400
,
,
2
,
1
; the mapped
grid is of order
64
400
u
.
62
For the purposes of the present example, we choose the following
parameters:
u
p
A
.
2
1
,
u
p
A
.
1
2
, damping factor
u
p.
02
.
0
D
, wave numbers,
0
1
k
,
32
2
2
S
k
in
1
.
m
rad
and angular frequency,
16
2
1
S
Z
,
64
2
2
S
Z
in
s
rad
, with evenly spaced contours having constant slopes equal to
5
.
0
2
2
k
Z
.
As can be seen from (3.74) for the given data, the function
1
f
is a damped
standing wave whereas
2
f
is a travelling wave. Figure 3.3 shows a plot of the
space-time evolution of
)
,
( t
x
f
. In this numerical example, complex EOF
analysis was applied from which a standing mode and a travelling mode are
estimated. Each time series is augmented with an imaginary component, where
the Hilbert transform is used to obtain the complex representation in order to
provide phase information.
Much insight into the nature of standing and propagating features can be
found by examining the space-time evolution of selected variables.
In the results that follow, the standing and travelling wave component were
obtained to validate the proposed methodology. The standing wave mode
contains the 19.5 % of the energy while the travelling wave mode contains the
79.5 %, and the remainder 1% is distributed in approximated error.
Fig. 3. 3. Space-time evolution of propagating wave on proposed field.
63
A) Standing wave: This example serves to illustrate the analysis of standing
wave functions from simulated (measured) data. The importance of
determining the standing wave component can be understood by noting that it
is the stationary wave present at different locations of the original field, that is
to say, it represents a local oscillation.
The estimated standing wave is shown in Fig. 3.4. Figure 3.5 compares the
Fourier transform spectrum of standing wave computed using the fast Fourier
transform (FFT) and the instantaneous frequency given from proposed
algorithm.
Both methods are found to be in good agreement. It can be seen from these
plots that the frequency,
1
Z , and amplitude given in Fig. 3.6 coincide with
expected behavior, thus showing the correctness of the proposed approach.
Fig. 3. 4. Standing wave estimated from original field.
64
a) Proposed technique.
b) FFT algorithm.
Fig. 3. 5. Frequency spectrum.
Fig. 3. 6. Amplitude of the standing wave
65
B) Travelling wave: The analysis of spatial changes in dynamic patterns of
oscillatory systems is a problem of great theoretical and practical importance.
The extraction of dynamic patterns from wide-area measurements using EOF
analysis has been the main objective in the proposed method. Knowledge of the
travelling wave component is of great utility for the study of propagating
phenomena in power systems. The detection and characterization of travelling
wave features is used to diagnostic and detect interactions into different areas.
Knowledge of travelling features may help to some types of devices such as
flexible alternating current transmission system (FACTS) controllers to be tuned
over a range of frequencies. The travelling wave component computed from
numerical example is shown in Fig. 3.7. As can be seen in these plots, the results
obtained from proposed method are consistent with the observed oscillations
thus giving validity to the results.
Fig. 3. 7. Travelling wave estimated from field original.
A point of particular interest is the ability of the proposed technique for the
detection of instantaneous features such as the instantaneous frequency. In Fig.
3.8 is plotted the frequency spectrum of the travelling wave obtained using the
FFT algorithm and the proposed technique. Note that the computed frequency
coincides with
2
Z in (3.74), the approximated instantaneous frequency is
computed as in subsection 3.4.5.
66
Also of relevance is the wave number
)
(k
through space; this is equivalents
to the angular frequency in the time domain. Figure 3.9 compares the frequency
spectrum using the FFT algorithm and the instantaneous wave number
computed using the proposed method.
a) Proposed technique.
b) FFT algorithm.
Fig. 3. 8. Frequency spectrum.
67
a) Proposed technique.
b) FFT algorithm.
Fig. 3. 9. Wave constant.
These findings provide basic insights into the nonstationary behavior of
general time series. This may be useful for understanding and predicting the
dynamic performance of measured data.
68
3.6 Concluding Remarks
In this chapter, a new method of temporal representation of nonstationary
processes in power systems has been presented. Complex empirical orthogonal
function analysis provides an efficient and accurate way of looking at the
temporal variability of transient processes while at the same time providing
spatial information about each one of the dominant modes with no a priori
assumption on either spatial or temporal behavior. The main advantage of this
approach is its ability to compress the variability of large data sets into the
lower possible number of temporal modes.
Complex empirical orthogonal function analysis is shown to be a useful
method for identifying standing and travelling patterns in wide-area system
measurements. Numerical results show that the proposed method provides
accurate estimation of nonstationary effects, modal frequency, time-varying
mode shapes, and time instants of intermittent transient responses. This
information may be important in determining strategies for wide-area control
and special protection systems. The identified system modes from the
decomposition may also serve to reveal relevant, but unexpected structure
hidden in the data such as that resulting from short-lived transient episodes.
Other issues such as the effect of numerical approximations on modal estimates
will be investigated in subsequent chapters.
3.7 References
[1] H. Dankowicz, P. Holmes, G. Berkooz, and J. Elezgaray, "Local models of spatio-
temporally complex fields," Physica D, vol. 90, pp. 387-407, 1996.
[2] N. Aubry, R. Guyonnet, and R. Lima, "Spatiotemporal analysis of complex signals:
theory and applications," Journal of Statistical Physics, vol. 64, nos. 3/4, 1991.
[3] B. F. Feeny, "A complex orthogonal decomposition for wave motion analysis,"
Journal of Sound and Vibration
, vol. 310, pp. 77-90, 2008.
[4] A. R. Messina, V. Vittal, D. Ruiz-Vega, and G. Enriquez-Harper, "Interpretation
and visualization of Wide-area PMU measurements using Hilbert analysis," IEEE
Transactions on Power Systems
, vol. 21, pp. 1763-1771, no. 4, Nov., 2006.
69
[5] M. A. Merrifield, and R. T. Guza, "Notes and correspondence. Detecting
propagating signals with complex empirical orthogonal functions: a cautionary
note," Journal of physical oceanography, pp. 1628-1633, Oct., 1990.
[6] S. Toh, "Statistical model with localized structures describing the spatial-temporal
chaos of Kuramoto-Sivashinsky equation," Journal of the Physical Society of Japan,
vol. 56, no. 3, pp. 949-962, Mar., 1987.
[7] J. F. Hauer, and J. G. DeSteese, "A tutorial on detection and characterization of
special behavior in large electric power systems," Pacific Northwest National
Laboratory
, July 2004.
[8] A. Hannachi, I. T. Jolliffe, and D. B. Stephenson, "Review. Empirical orthogonal
functions and related techniques in atmospheric science: a review," International
Journal of Climatology
, vol. 27, pp. 1119-1152, May 2007.
[9] G. Kerschen, J.-C. Golinval, A. F. Vakakis, and L. A. Bergman, "The method of
proper orthogonal decomposition for dynamical characterization and order
reduction of mechanical systems: An overview," Nonlinear Dynamics, vol. 41, pp.
147-169, 2005.
[10] P. A. Parrilo, S. Lall, and F. Paganini, "Model reduction for analysis of cascading
failures in power systems," Proceedings of the American Control Conference, pp. 4208-
4212, 1999.
[11] A. R. Messina, "Inter-area Oscillations In Power Systems: A Nonlinear and
Nonstationary Perspective," Springer Verlag, 2009.
[12] S. L. Hahn, "Hilbert transforms in signal processing," The Artech House, Signal
Processing Library, 1996.
[13] T. P. Barnett, "Interaction of the Monsoon and pacific trade wind system at
interannual time scales. Part 1: the equatorial zone," Monthly weather review, vol.
111, pp. 756-773, 1983.
[14] R. D. Susanto, Q. Zheng, and X.-H. Yan, "Complex singular value decomposition
analysis of equatorial waves in the pacific observed by TOPEX/Poseidon
altimeter," Journal of Atmospheric and oceanic technology, vol. 15, pp. 764-774, 1998.
[15] C. Wolter, M. A. Trindade, and R. Sampaio, "Obtaining mode shapes through the
Karhunen-Loève expansion for distributed-parameter linear systems," Dynamics
and Vibration
, pp. 1-21, 2000.
[16] J. D. Horel, "Complex principal component analysis: Theory and examples,"
Journal of climate and Applied Meteorology
, vol. 23, pp. 1660-1673, 1984.
[17] D. J. Trudnowski, "Estimating Electromechanical Mode Shape from
Synchrophasor measurements," IEEE Transactions on Power Systems, vol. 23, no. 3,
pp. 1188-1195, August 2008.
70
Chapter 4
Near-Real-Time of Measured Data
Using Complex EOF Analysis
Extracting the essential dynamics in the underlying process of complicated
phenomena is a complex problem; this is particularly true in the study of the
wide-area oscillatory phenomena where many data are necessary.
In this chapter, a statistically-based, data-driven framework that integrates
the use of EOF analysis and a sliding window-based method is proposed to
identify and extract in near-real-time, dynamically independent spatio-
temporal patterns from time synchronized data. The approach is nonparametric
and multiscale and deals with the information in space and time
simultaneously, which result in a more compressed representation of the data.
This can lead to a greater understanding of the mechanisms which generate the
time series and result in a better diagnosis of specific events leading to the
observed oscillatory behavior.
The method is expected to have various applications in dynamic analysis
including on-line analysis of mode propagation, ambient noise measurements
and dynamic coherency. More importantly, the technique has the potential to
be applied for near-real-time in wide-area monitoring and analysis.
71
4.1 Motivation
With the growing utilization of measurement devices, the need for
measurement-based techniques that allow the extracting of spatial and
temporal features of measured data is becoming more demanding [1-6].
Power system response is primarily driven by random load variations,
generation variations, unforeseen system outages and protective relaying which
make its analysis and control difficult. This random nature of the system
behavior has been largely disregarded or inadequately assessed in applications
of modal analysis techniques for oscillatory problems in power systems.
Because of these inherent problems, deterministic approaches are not adequate
to accommodate continuous changes in system planning and operation, and can
result in unrealistic estimates of system behavior.
As discussed in the introductory section, event-based oscillations cannot be
evaluated in isolation from other events. Moreover, deterministic procedure for
power system analysis and design used today may provide a conservative
estimate of system operation, and lead to inadequate design of controllers [7-
11]. To establish operating limits and allows efficient design of system
controllers, the random nature of these variations needs to be realistically
modeled and simulated.
This is especially important in the analysis and characterization of wide-area
system oscillations in which many dynamic phenomena may take place, often
involving different time scales [12-13].
4.2 Analytical Aspects of Time-Dependent EOFs
On-line monitoring of large-scale power systems by means of time-
synchronized phasor measurement units (PMUs) provides the opportunity to
analyze and characterize inter-system oscillations. Wide-area measurements
sets, however, are often relatively large, and may contain phenomena with
differing temporal scales. Extracting the relevant dynamics from these
measurements is a difficult problem. As the number of observations of real
72
events continues to increase, statistical techniques are needed to help in identify
relevant temporal dynamics from noise or random effects in measured data.
In the last decade, extensive research has been carried out for developing
modal extraction algorithms that can accurately monitor transient response.
Among these approaches, the EOF analysis has been used to extract dynamic
patterns from measured data.
The EOF analysis is not only useful as method of constructing optimal
dynamical models, but also provides a simple diagnostic tool for compressing
the extensive information required for a complete characterization of data of the
process into a manageable set of modes.
Difficulties associated with interpreting the swing dynamics of wide-area
have led researchers to developed tools to overcome these difficulties such as
extensions of existent techniques.
Conventional EOF analysis assumes stationariness of the underlying
dynamic process. Such an approach can only provide information about
standing waves and cannot isolate portions of the phenomena in which
dynamical changes take place. In order to overcome the above limitations, a
sliding window-based method is used to systematically analyze the
observational data in near-real-time. This allows to isolate the portions of data
where oscillations are present and improves numerical efficiency and accuracy.
4.3 Near-Real-Time EOF Analysis
Time-synchronized phasor measurements collected by PMUs or other dynamic
recorders can be interpreted conveniently in terms of statistical models
involving both temporal and spatial variability [1]. When data is available at
multiple locations, a spatio-temporal model is required that represents system
behavior. Figure 4.1 shows a conceptual representation of PMU data collected
at different spatial locations.
73
Fig. 4. 1. Visualization of PMU data in terms of spatial and temporal information.
Assume, in order to introduce the more general ideas that follow, that
measured data is available at n spatial (measurement locations) defined by
j
x
,
n
j
,...,
1
at
N
instants in time,
k
t
,
N
k
,...,
1
. Let
)
,
(
k
j
t
x
u
be a space-time scalar
field representing a time trace, where
j
x
, is a set of spatial variables (i.e.,
measurement locations) on a space
:
, and
k
t
, is the time at which the
observations are made. Using this notation, the set of data can be represented
by an
n
N
u
-dimension ensemble (observation) matrix,
)
,
( t
x
X
, of the form
@
»
»
»
»
¼
º
«
«
«
«
¬
ª
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
2
1
2
2
2
2
1
1
1
2
1
1
2
1
N
n
N
N
n
n
n
t
x
u
t
x
u
t
x
u
t
x
u
t
x
u
t
x
u
t
x
u
t
x
u
t
x
u
t
x
t
x
t
x
t
x
#
%
#
#
u
u
u
X
In this formulation, each column corresponds to the system response at a
specific time, and can represent the system response to a single event or
represent an ensemble of time responses to multiple events measured at a
single location.
As it was highlighted in the previous section, conventional EOF analysis
assumes stationarity of the underlying dynamic process and is therefore not
74
suitable for an on-line implementation. In addition, such an approach, can only
provide information about standing waves and cannot isolate portions of the
phenomena in which dynamic changes take place [1,7].
In order to overcome the above limitations, a sliding window-based method
is used to systematically analyze the observational data in near-real-time. This
allows to isolate and extract the portions of data where oscillations are present
and improves numerical efficiency and accuracy.
4.3.1 Time-Dependent Covariance Matrix
Time series are recorded in discrete form even though the underling process
itself is continuous [14,15]. For discretely measured data, a sliding window-
based method is proposed here that allows recursive computation of the
covariance matrix.
Assume that
)]
,
(
,
),
,
(
),
,
(
[
)
,
(
2
1
k
n
k
k
k
t
x
u
t
x
u
t
x
u
t
x
u
denotes a sequence of
observations, collected at time instant
k
t
, to
f
,...,
1
k
, evenly spread throughout
the time period, where
n
x
is set of spatial variables (measurement locations).
Assume, in order to introduce these ideas, that
)
,
(
k
t
x
u
is a
n
u
1
vector-file
with complex or real values that represents the snapshots at time instant
k
t
, to
f
,
,
2
,
1
k
, on some domain
:
x
. Hereafter,
)
,
(
k
t
x
u
is referred as
)
(
k
x
t
u
to
enhance the time-dependent analysis. Then,
)]
(
,
),
(
),
(
[
)
(
2
1
k
x
k
x
k
x
k
x
t
u
t
u
t
u
t
n
u
where
n
x
is set of spatial variable defined over the range
)
,
1
[
f
.
A sliding window-based approach has been combined with EOF analysis to
resolve localized information. In this approach, a sliding window frame of fixed
size, say
W
, is shifted regularly throughout the data span from the beginning of
the record to the end of the data as is shown in Fig. 4.2. The analysis window is
then slid the same distance repeatedly and the covariance matrix is computed
recursively. It should be stressed that
W
can range from a single sample to the
entire length of the record.
75
Fig. 4. 2. Conceptual diagram of the sliding window-based approach.
Referring to Fig. 4.2, let matrix
)
(
k
x
t
X
be an ensemble of measured data of
size
n
u
W
, where
W
is the window size, n is the number of sensors or (PMUs),
and k indicates the time instant at which the response matrix is computed.
Using the theory in section 3.3, the time-dependent response matrix of size
n
u
W
,
can then be defined as
k
to
t
u
t
u
t
u
t
u
t
u
t
u
t
u
t
u
t
u
t
n
n
n
x
x
x
x
x
x
x
x
x
k
x
,...,
2
,
1
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
2
1
2
1
2
1
2
2
2
1
1
1
»
»
»
»
»
¼
º
«
«
«
«
«
¬
ª
W
W
W
W
#
%
#
#
X
(4.1)
More formally, given a set of available sample time series, the temporal
autocorrelation matrix,
)
(
k
x
t
C
, of
)
(
k
x
t
X
extending over a window length is
defined as
¦
'
k
j
j
H
x
j
x
k
H
x
k
x
k
x
t
t
k
k
t
t
t
1
)
(
)
(
1
)
(
)
(
)
(
u
u
X
X
C
, (4.2)
or, equivalently,
.
)
(
)
(
1
)
(
)
(
)
(
1
¦
'
k
j
j
x
j
H
x
k
x
k
H
x
k
x
t
t
k
k
t
t
t
u
u
X
X
C
(4.3)
where, in the above, the covariance matrix can be real or complex.
76
In our formulation, matrix
)
(
k
x
t
C
, is positive semidefinite matrix of size
n
n
u , and possesses a set of complex time-dependent orthogonal eigenvectors,
)
(
)
(
x
k
t
i
M
,
n
i
,
,
2
,
1
, that capture propagating features. Extensions of this
approach to consider the case where the autocorrelation matrix is nonintegrable
are discussed in [16,17]. This is, however, not considered here.
An interesting and useful interpretation of the covariance matrix can be
obtained by rewriting (4.2) as
)
(
)
(
)
(
)
(
)
(
1
1
k
H
x
k
x
k
j
j
H
x
j
x
k
x
t
t
t
t
t
k
u
u
u
u
C
¦
(4.4)
Now if we let
¦
1
1
1
)
(
)
(
1
1
)
(
k
j
j
H
x
j
x
k
x
t
t
k
t
u
u
C
, we can
)
(
)
(
1
)
(
1
)
(
1
k
H
x
k
x
k
x
k
x
t
t
k
t
k
k
t
u
u
C
C
, (4.5)
where the increment
)
(
' of
C
at time instant
k
is given as
)
(
)
(
1
)
(
1
1
k
x
k
x
k
x
t
t
k
t
k
t
u
u
C
C
'
'
. It then follows that, in the limit,
0
lim
'
'
of
t
k
C
.
This analysis suggests that the covariance matrix at time step
k
can be
computed recursively using information from both, the previous time step
)
1
(
k
and the current measurement using only elementary operations. The
stationary test for the recursive covariance matrix can be applied to
demonstrate the stationarity of the dataset [18].
The above formulation can be extended in straightforward manner to the
real-time estimation of temporal variations in the data.
Similar to the previous developed, it can be shown that for a fixed window
size,
W
, (4.2) and (4.3) can be expressed as
)
(
)
(
1
)
(
1
)
(
1
k
H
x
k
x
k
x
k
x
t
t
t
t
u
u
C
C
W
W
W
, (4.6)
77
where
¦
1
1
1
)
(
)
(
1
1
)
(
W
W
j
j
x
j
x
k
x
t
t
t
u
u
C
, besides that, its increment is given as
)
(
)
(
1
)
(
1
1
k
x
k
x
k
x
t
t
t
t
u
u
C
C
W
W
'
'
, leading to
0
lim
'
'
of
t
C
W
.
Now, it is clear that the information at time instant
k
depends on the
information at the previous. This suggests that the temporal behavior of a
harmonic wave can be accurately characterized by sliding the window at each
sampling period.
Equation (4.6) provides an efficient method to approximate the time-
dependent covariance matrix: since the method is based on local information
the technique is well-suited for real-time applications. The selection of an
optimal window size is an important but difficult problem and will be
addressed in future research. In the analysis of highly nonstationary signals,
large values of may obscure the analysis of temporal changes occurring at
segments of the observed record. As discussed in our numerical simulations,
short-width windows provide improved localized information but result in
enhanced computational effort. We have found that, the best linear estimate for
the autocorrelation matrix
C
can be obtained using the averaged two-point
correlation function [19]. The computation of the covariance matrix using linear
estimate is given in Appendix C.
We next examine the properties of the temporal covariance matrix and
introduce a numerical algorithm to extract modal information.
4.3.2 Spectral Decomposition of Time-Dependent Covariance Matrix
Recent papers have examined the structural properties of the covariance matrix.
In this section, we extend this theory to the case of time-varying system
formulations. Techniques to extract the dominant modes of behavior in vector
data sets are then discussed.
In an effort to extend standard real-EOF analysis to deal with propagating
phenomena a complex formulation is adopted. Here, the real part is augmented
with an imaginary obtained from the Hilbert transform of each time series.
78
Assume the that
)
(
)
(
)
(
)
(
)
(
k
x
I
k
x
R
k
x
t
i
t
t
X
X
X
is a complex matrix with
1
i
, where the subscripts
I
R,
, indicate the real and imaginary vectors. Under these
assumptions, the autocorrelation matrix
)
(
k
x
t
C
in (4.2) and (4.3) becomes
)
(
)
(
)
(
)
(
)
(
k
x
I
k
x
R
k
x
t
i
t
t
C
C
C
(4.7)
where
k
R
I
I
R
k
I
I
R
R
k
x
k
H
x
k
I
R
R
I
k
I
I
R
R
k
H
x
k
x
i
t
t
i
t
t
]
[
]
[
)
(
)
(
]
[
]
[
)
(
)
(
T
T
T
T
T
T
T
T
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
. (4.8)
Given a model of the form (4.7), it is straightforward to show that
)
(
)
(
T
)
(
)
(
k
x
R
k
x
R
t
t
C
C
is a symmetrical matrix, and that
)
(
)
(
k
x
I
t
C
is a skew-
symmetric matrix, this is,
)
(
)
(
)
(
T
)
(
k
x
I
k
x
I
t
t
C
C
. Since the symmetrical matrix is a
particular case of the Hermitian matrix, all of its eigenvectors are real; the
elements of the skew-symmetric matrix are all purely imaginary and its
eigenvectors are complex conjugate. Note also that the bases for the complex
autocorrelation matrix
)
(
k
x
t
C
, are now defined as
)
(
)
(
x
k
t
R
for the real part and
)
(
)
(
x
k
t
I
for the imaginary part.
4.3.3 Time-Dependent Basis Functions
As discussed earlier in this chapter, the fundamental idea of the time-
dependent proper orthogonal decomposition is find a instantaneous basis
)
(x
k
t
for linear, infinite-dimensional Hilbert space
])
1
,
0
([
2
L
that maximize the
averaged projection of a function
)
(
k
x
t
u
on some domain
:
x
, suitably
normalized
1
)
(
subject to
)
(
|
))
(
),
(
(
|
max
2
2
2
])
1
,
0
([
2
x
x
x
t
k
k
k
t
t
t
k
x
L
u
(4.9)
Similar to the previous development, it follows that, the optimal basis are
given by the eigenfunctions
)
(x
k
t
of the equation (4.7). They are hereafter
called the time-dependent empirical orthogonal functions [19].
79
4.3.4 Feature Extraction
Many power system phenomena derive from nonlinear interactions between
travelling waves of different spatial scales and temporal frequencies. This
section extends the developed model to detect propagating phenomena in
nonstationary processes.
Once the spatial eigenvectors associated with the real and imaginary part of
C
are calculated, the original field can be approximated by a spatio-temporal
model [1]. Consider a field
)
(
k
x
t
u
composed of standing and travelling
components, denoted by
)
(
),
(
)
(
)
(
k
x
twc
k
x
swc
t
t
u
u
of the form
)
(
)
(
)
(
)
(
)
(
k
x
twc
k
x
swc
k
x
t
t
t
u
u
u
(4.10)
From EOF analysis [1], the complex field can be expressed as the complex
expansion
¦
¦
q
j
j
I
k
j
I
j
I
k
j
I
k
x
twc
p
j
k
j
R
j
R
k
j
R
k
x
swc
x
t
x
t
t
t
x
t
t
1
)
(
)
(
)
(
)
(
)
(
1
)
(
)
(
)
(
)
(
)
cos(
)
(
)
(
)
(
)
cos(
)
(
)
(
)
(
S
K
S
R
u
S
R
u
(4.11)
for
n
q
p
d
,
, with q even. In Eq. (4.11), the unknown instantaneous coefficients
, are the temporal and spatial amplitude functions, respectively, and ,
K
are the frequency functions and wave component to be determined [9,20].
Unlike the stationary case, the time-dependent expressions in (4.11) need to
incorporate temporal information.
4.3.5 Wave Components
Let now the complex field (4.10) be expanded in terms of a truncated EOF basis
of
p
and
q
modes as [1]
¦
¦
q
j
H
j
I
k
j
I
p
j
H
j
R
k
j
R
k
x
x
t
i
x
t
t
1
)
(
)
(
1
)
(
)
(
)
(
)
(
)
(
)
(
)
(
MM
M
A
A
u
(4.12)
80
where the complex, time-dependent coefficients,
)
(
)
(
k
j
R
t
A
and
)
(
)
(
k
j
I
t
A
are
given by projecting the bases functions onto the original vector field [20,21].
That is,
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
x
t
t
x
t
t
k
k
t
I
k
x
k
x
I
t
R
k
x
k
x
R
u
A
u
A
(4.13)
Following a similar approach to that described in section 3.4, Equation (4.12)
can now be recast in the form
¦
¦
q
j
i
j
I
k
j
I
p
j
i
j
R
k
j
R
k
x
j
I
j
I
j
R
j
R
e
x
t
e
x
t
t
1
(
)
(
)
(
1
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
II
T
I
T
S
R
S
R
u
(4.14)
where
)
(t
R
and
)
(t
S
are the temporal and spatial amplitude functions
associated to the wave component of the decomposition, respectively, and
)
(t
TT
and
)
(t
II
are the temporal and spatial phase functions, respectively. The original
time series can then be reconstructed from the real part of (4.14). The details are
given in section 3.4.
Now, the temporal parameters of (4.14) are computed as
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
x
x
x
x
x
x
t
t
t
t
t
t
k
k
k
k
t
I
t
I
I
t
R
t
R
R
k
x
I
k
x
I
k
I
k
x
R
k
x
R
k
R
S
S
A
A
R
A
A
R
+
+
+
+
(4.15)
Similarly, the angular frequency
R
,
I
is defined as
t
t
t
t
t
t
t
t
k
I
k
I
k
I
k
R
k
R
k
R
)
(
)
(
)
(
,
)
(
)
(
)
(
1
1
(4.16)
and wave number
)
(
)
(
x
k
t
I
x
x
x
x
k
k
k
t
I
t
I
t
I
)
(
)
(
)
(
)
(
)
(
)
(
1
I
I
(4.17)
81
The temporal phase function
)
(
k
R
t
,
)
(
k
I
t
and the spatial phase function
)
(
)
(
x
k
t
I
II
given in (4.16) and (4.17) are computed from
°¿
°
¾
½
°¯
°
®
°¿
°
¾
½
°¯
°
®
°¿
°
¾
½
°¯
°
®
)]
(
[
)]
(
[
tan
)
(
)]
(
[
)]
(
[
tan
)
(
)]
(
[
)]
(
[
tan
)
(
)
(
)
(
1
)
(
)
(
)
(
1
)
(
)
(
1
x
re
x
im
x
t
re
t
im
t
t
re
t
im
t
k
k
k
t
I
t
I
t
I
k
x
I
k
x
I
k
I
k
x
R
k
x
R
k
R
A
A
A
A
I
(4.18)
In the developed models, a criterion for choosing the number of relevant
modes is given by the energy percentage contained in the
p
and
q
modes.
From the above it is easy to see that
%
99
100
)
(
)
(
)
(
)
(
2
1
)
(
1
)
(
u
¦
¦
F
k
x
q
j
k
twc
j
p
j
k
swc
j
k
x
t
t
t
t
E
u
O
O
(4.19)
where
2
.
F
denotes the Frobenius norm. From (4.19), the increment
)
(
k
swc
t
O
and
)
(
k
twc
t
O
are given in difference as
)
(
)
(
)
(
and
)
(
)
(
)
(
1
1
k
twc
k
twc
k
twc
k
swc
k
swc
k
swc
t
t
t
t
t
t
O
O
O
O
O
O
(4.20)
If a sliding window is used to systematically analyze the observational data
in near-real-time, then
)
(
)
(
)
(
)
(
1
1
'
W
O
O
O
O
k
swc
k
swc
k
swc
k
swc
t
t
t
t
(4.21)
and
)
(
)
(
)
(
)
(
1
1
'
W
O
O
O
O
k
twc
k
twc
k
twc
k
twc
t
t
t
t
(4.22)
These equations provide a complete characterization of any propagating
effects and periodicity into original field which might be obscured by normal
cross-spectral analysis.
82
4.4 Characterization of Spatio-Temporal Behavior
In the previous sections, it was shown that dynamic behavior from observed
data can be characterized by four measures given in the proposed method. The
space-time instantaneous description provides a compact way of looking the
spatial and temporal variability of the observed data [22-24].
In what follows we discuss the spatial and temporal properties used to the
modal characterization.
4.4.1 Temporal Characterization
The variation in the instant of time of the modal structure is given by means of
the description of its temporal amplitude and phase function [20].
As discussed in Chapter 3, several measures that define possible moving
features in the vector field can now be defined.
x Temporal amplitude functions,
)
(
k
R
t
R
,
)
(
k
I
t
R
. In analogy to the real,
stationary case, a useful measure of temporal variability in the magnitude of
the modal structure of the field is obtained from (4.15).
x Temporal phase function,
)
(
k
R
t
,
)
(
k
I
t
. This function describes the
temporal variation of phase associated with
)
(
k
x
t
u
.
x Temporal energy: Their evolution dynamic permits identify transient
phenomena present in the field.
When combined with temporal information, these parameters completely
define modal structure. The temporal characterization defines the temporal
evolution of modal structure; it can be used to estimate deviations from steady-
state performance at a particular location. In other words, it provides a measure
of the standing wave components present in the original field.
4.4.2 Spatial Characterization
The position of any type of wave in the space field can be measured by its
spatial amplitude and phase function [20]. More formally,
83
x Spatial amplitude functions,
)
(
)
(
x
k
t
R
S
,
)
(
)
(
x
k
t
I
S
: Shows the spatial
distribution of variability in the field associated in the spatial structure.
x Spatial phase function,
)
(
)
(
x
k
t
I
II
: Shows the relative phase of fluctuations
among the various spatial locations.
x Spatial energy: Their dynamical evolution permits identify too transient
phenomena like in the case of the temporal energy, this uses the spatial
relation.
The spatial characterization is used in this document to define the travelling
wave component; it can be used to detect the extent of modal interaction at any
physical locations, their spatial amplitude functions give the distribution in the
field.
4.5 Proper Orthogonal Decomposition of Nonstationary
Random Processes
Since the EOF analysis depends on the decomposition of a time-dependent
data-driven covariance matrix, a technique is needed to estimate this matrix in
near-real-time. Conventionally, this is done based on the moment estimation
procedure. In this section, a moving average-based computing method for real
time computation of the instantaneous mean is proposed.
We recall from subsection 4.3.1 that the covariance matrix for the discrete
case can be written as
)
(
)
(
1
)
(
1
)
(
1
k
H
x
k
x
k
x
k
x
t
t
k
t
k
k
t
u
u
C
C
(4.23)
where
)
(
k
x
t
u
is assumed to have zero mean.
To complete the statistical characterization of system behavior in near-rea-
time it is now necessary to determine an appropriate measure of the temporal
mean behavior [25,26]. We first examine some statistical properties.
84
4.5.1 Statistics of Random Processes
In order to introduce the definition and property statistics of random processes,
a random experiment with a finite or infinite number of outcomes from a
sample space
:
x
, each occurring with a probability
)
)
(
(
x
t
P
d
u
, is considered
[2-3].
A stochastic process is a noncountable infinity of random variables, one for
each
t
. For a specific
t
,
)
(t
u
is a random variable with distribution
}
)
(
{
)
,
(
x
t
P
t
x
F
d
u
(4.24)
This function depends on
t
, and is equal to the probability of the event
}
)
(
{
x
t
d
u
consisting of all outcomes
[ such that, at the specific time
t
, the
samples
)
,
(
[
t
u
of the given process do not exceed the number x . Of central
importance in this methodology will be random variables that are Gaussian, or
equivalently, normal distributed. The function
)
,
( t
x
F
in (4.24) will be called the
first-order distribution of the process
)
(t
u
. Its derivative with respect to x
is:
x
t
x
F
t
x
f
w
w
)
,
(
)
,
(
(4.25)
gives the first-order density of
)
(t
u
. A random variable is Gaussian if its density
function is given by
»
»
¼
º
«
«
¬
ª
¸¸¹
·
¨¨©
§
2
)
(
2
/
1
2
)
(
)
(
2
1
exp
)
)
(
2
1
(
)
,
(
t
t
u
x
t
t
x
f
x
x
m
x
V
SV
(4.26)
where
)
(
)
(
)
(
x
E
t
u
x
m
represents the expectation value or mean, and
]
))
(
[(
)
(
2
x
E
x
E
t
x
V
is the standard deviation, which represents a level of added
white noise; note that these two parameters characterize the Gaussian density.
The normal density function has the familiar bell shape. The normal
distribution function in (4.24) is then given by
[
V
[
SV
d
t
t
u
t
t
x
F
x
x
x
m
x
³
f
»
»
¼
º
«
«
¬
ª
¸¸¹
·
¨¨©
§
2
)
(
2
/
1
2
)
(
)
(
2
1
exp
)
)
(
2
1
(
)
,
(
(4.27)
85
A random variable u is a process of assigning a number
)
(
[
u
to every
outcome
[ . The resulting function mush satisfies the following two conditions
but is otherwise arbitrary:
I. The set
}
{
x
d
u
is an event for every x .
II. The probabilities of the events
0
}
{
f
u
P
and
0
}
{
f
u
P
.
We shall say that the statistics of a random variable u are known if we can
determine the probability
}
{
R
P
u
that u is in a set R of the x axis consisting
of countable unions or intersections of intervals. The following sections provide
the nature of the approximation to construct the time-dependent statistical
representation using (4.26) and (4.27).
4.5.2 Stationarity of Random Process
Consider a random process,
)
(
i
t
u
, of length k (
k
i
,
,
2
,
1
) and other group of
samples displaced in the time a quantity from first group, the second group is
)
(
W
i
t
u
(see Fig. 4.2). These second group of
k
random variables is characterized
by their joint probability density functions (PDFs)
)
,
,
;
,
,
(
W
W
W
W
k
i
k
i
t
t
x
x
f
. The
joint PDF of the groups of random variables can be or not identical. If the
functions are identical, then
)
,
,
;
,
,
(
)
,
,
;
,
,
(
W
W
W
W
k
i
k
i
k
i
k
i
t
t
x
x
f
t
t
x
x
f
(4.28)
for all and
k
. It is said that the process is stationary in strict sense. In other
words, the statistical properties of a stationary random process are invariants
above a translation of time [14,15].
4.5.3 Moments of Random Variables
Given a stationary random process,
)
(
k
t
u
, the first moment of the process can
be expressed in terms of the expectation value as
³
f
f
))
(
(
))
;
(
)
(
))
(
(
k
k
k
k
k
t
u
d
t
x
f
t
u
t
E u
(4.29)
86
where
))
(
(
k
t
E u
represents the expectation value, and the value of the moment
depend of the instant
k
t
. The first moment is the mean of the signal.
If the statistical correlation into
)
(
k
t
u
and
)
(
1
k
t
u
is measure through of joint
moment
³ ³
f
f
f
f
))
(
(
))
(
(
)
,
;
,
(
)
(
)
(
))
(
)
(
(
1
1
1
1
1
k
k
k
k
k
k
k
k
k
k
t
u
d
t
u
d
t
t
x
x
f
t
u
t
u
t
t
E
u
u
(4.30)
then, it is noted that this joint moment depend of the time instants
k
and
1
k
,
(4.30) is called the autocorrelation function of the random process. This
implicates that the autocorrelation function is dependent upon the
autocorrelation values between times
k
and
1
k
[14,15]. If the process is
stationary, then,
))
(
)
(
(
))
(
)
(
(
1
1
1
W
k
k
k
k
t
t
E
t
t
E
u
u
u
u
(4.31)
It can also be stated that, if the process is nonstationary with the property
that the mean value is a constant, then the autocorrelation function can be
related to the autocovariance function, and is given by
)]}
(
)
(
)][
(
)
(
{[
1
1
k
m
k
k
m
k
t
u
t
t
u
t
E
u
u
C
(4.32)
where
))
(
(
)
(
k
k
m
t
E
t
u
u
and
))
(
(
)
(
1
1
k
k
m
t
E
t
u
u
are the mean value of
)
(
k
t
u
and
)
(
1
k
t
u
respectively.
4.5.4 Expectation Value and Standar Deviation
The proposed real-time formulation described above requires the efficient
computation of the temporal expectation value and standard deviation [2-3].
This can be done recursively as discussed below.
The discrete-time standard derivation on the interval
k
i
d
d
1
for the j-th
field can be expressed as
¦
k
i
k
x
m
i
x
k
x
k
x
k
x
t
u
t
u
k
t
E
t
E
t
j
j
j
j
j
1
2
)
(
2
))
(
)
(
(
1
1
]
)))
(
(
)
(
[(
)
(
u
u
V
(4.33)
Now, it is easily verified that (4.33) can be rewriten as
87
2
)
(
1
1
2
)
(
))
(
)
(
(
))
(
)
(
(
)
(
)
1
(
k
x
m
k
x
k
i
k
x
m
k
x
k
x
t
u
t
u
t
u
t
u
t
k
j
j
j
j
j
¦
V
(4.34)
where the recursive form to (4.33) is
2
)
(
1
))
(
)
(
(
1
1
)
(
1
2
)
(
k
x
m
k
x
k
x
k
x
t
u
t
u
k
t
k
k
t
j
j
j
j
V
V
(4.35)
or into arbitrary windows wide ,
2
)
(
1
))
(
)
(
(
1
1
)
(
1
2
)
(
k
x
m
k
x
k
x
k
x
t
u
t
u
t
t
j
j
j
j
W
V
W
W
V
(4.36)
We remark that, in equations (4.35) and (4.36), it is necessary to compute the
temporal expectation or mean value at time instant
k
t
. Assume to this end, that
the cumulative mean value of a signal
)
(
i
x
t
j
u
in the time interval
k
i
d
d
1
is
approximated by
¦
k
i
i
x
k
x
k
x
m
t
u
k
t
E
t
u
j
j
j
1
)
(
)
(
1
1
))
(
(
)
(
u
for
n
j
N
k
,
,
2
,
1
,
,
2
,
1
(4.37)
where for convenience of notation, we will rewrite this expression as
)
(
))
(
(
)
(
)
(
))
(
)(
1
(
1
)
(
1
1
)
(
k
x
k
x
m
k
x
k
i
i
x
k
x
m
t
u
t
u
k
t
u
t
u
t
u
k
j
j
j
j
j
¦
(4.38)
Equation (4.37) can then be replaced by
)
(
1
1
))
(
(
1
)
(
1
)
(
)
(
k
x
k
x
m
k
x
m
t
u
k
t
u
k
k
t
u
j
j
j
(4.39)
Similarly to (4.39), the cumulative mean value for a window size , can be
calculated as the average from each previous value through the recursive
relation
)
(
1
1
))
(
(
1
)
(
1
)
(
)
(
k
x
k
x
m
k
x
m
t
u
t
u
t
u
j
j
j
W
W
W
. (4.40)
The relation above has been used as basis to a method of estimation to the
cumulative mean value and standard deviation at each time instant k .
88
The development given in this section indicates that the computation for
(4.26) and (4.27) are estimated by means of the proposed time-varying statistical
formulation, consequently, these equations are called time-dependent
distribution and density function respectively. The next section discuss in more
detail the nature of the approximations used here to compute propagating
features using real-time information.
4.5.5 Recursive Algorithm
The algorithm used to obtain a time-varying representation of modal behavior
can be summarized in the following way:
1. Let
)
(
k
x
t
u
represent the time evolution of a nonlinear time series at time
instant k . Obtain the complex form of the POD formulation. This is
implemented in two steps:
a) Form the covariance matrix,
)
(
k
x
t
C
b) Add the complex part to each column of
)
(
k
x
t
C
2. The second step is the real-time computation procedure involves the
descomposition of the data into a time-varying mean and a fluctuating
component. This is done using the theory in subsection 4.5.4.
3. Compute the fluctuation of the signal by calculating
)
(
)
(
)
(
)
(
k
x
m
k
x
k
x
t
u
t
u
t
u
j
j
j
, where
)
(
j
x
m
u
is the time-varying mean of the
signal.
4. Compute the complex field representation using (3.6).
5. Determine time-varying covariance matrices using (4.6).
6. Decompose the complex-vector field into standing and travelling waves,
and
7. Extract dominant modes of behavior.
A flowchart of the proposed approach is shown in Fig. 4.3. In this plot,
dashed boxes indicate the use of information from the previous time step
)
1
(
k
, and
)
(
E
indicate the expectation value or time-varying mean.
89
This approach makes feasible the analysis and characterization of transient
processes using real-time information. Variations to these approaches that
extend their practical use to the real of near-real-time stability assessment and
control are being investigated in Chapter 5.
The size of the window frame is a crucial parameter of the algorithm given
in Fig. 4.3. This parameter gives the number of data used in the average value;
large values of
W
may obscure the analysis of temporal changes occurring at
segments of the observed record. This is an issue that warrants further
investigation and will be investigated further in future stages of this work.
u
x
(t
k
)
Complex field
Complex EOF analysis
PMU
E
k-1
(u
x
(t
k-1
))
C
k-1
)
(
)
(
1
1
1
k
H
x
k
x
k
k
t
t
k
k
k
u
u
C
C
))
(
(
)
(
)
(
k
x
k
x
k
x
t
E
t
t
u
u
u
)
(
)
(
k
x
swc
t
u
)
(
)
(
k
x
twc
t
u
))
(
(
)
(
)
(
)
(
)
(
)
(
k
x
k
x
twc
k
x
swc
k
x
t
E
t
t
t
u
u
u
u
)
(
)
(
1
1
1
k
H
x
k
x
k
k
t
t u
u
C
C
W
W
W
Fig. 4. 3. Proposed algorithm to the complex EOF identification
The use of the EOF analysis to determine the instantaneous parameters of
one or any number of nonstationary time series is a problem treated in this
section. In the following are given two methodologies that in combination with
the EOF analysis are used to extract the relevant scales into time series.
90
4.6 Multiscale EOF Analysis
In real-world applications data are multiscale due to events occurring at
different locations with different localization in time and frequency scales [27-
29]. To address the problem of multiscale modeling in data obtained from time-
synchronized phasor measurements, a multiscale EOF approach has been
developed. Wavelet transformation (WT) and conventional empirical
orthogonal function (EOF) analysis are combined tentatively; the method
combines the ability of EOF to decorrelate the variables for extracting a linear
relationship with that of wavelet analysis to extract deterministic features and
approximately decorrelate autocorrelated measurements. In this approach, the
method computes the dominant EOFs of wavelet coefficients at each scale and
then combines the results at relevant scales.
The method has potential for analyzing complicated signals associated with
data containing contributions from events whose behavior changes over time
and frequency, and inherits the advantages of both WT and EOF.
4.6.1 Wavelet-Based EOF Analysis
As highlighted in previous sections, conventional EOF is best for analyzing data
collected from a steady-state process, containing linear relationships between
the variables. In this section, a multiscale EOF analysis methodology based on
wavelet analysis is introduced. The method consists of decomposing each
variable of interest on a selected family of wavelets. The EOF model is then
determined independently for the coefficients at each scale; the models at
important scales are then combined in an efficient scale-recursive manner to
yield the model for all scales together.
The wavelet analysis in the context of the proposed formulation is next
summarized.
Consider a function
)
(t
f
belonging to a finite energy function, that is
f
³
f
f
dt
t
f
2
)
(
(4.41)
91
The continuous wavelet transform (WT) of
)
(t
f
with respect to a mother
wavelet is a convolution integral given as [30-32]
dt
a
b
t
t
f
a
b
a
W
³
f
f
¸
¹
·
¨
©
§
\
*
)
(
1
)
,
(
(4.42)
where
)
(t
\ is the wavelet basis function or mother wavelet,
a
is a scale
parameter and b is a time parameter. The superscript (*) denotes complex
conjugation. It is assumed that
)
(t
\ is also a finite energy function satisfying the
condition
³
f
f
\
f
Z
Z
Z
\
d
C
)
(
(4.43)
in which
)
(
Z
\
is the Fourier transform of
)
(t
\ defined as
³
f
f
Z
\
S
Z
\
dt
e
t
t
i
)
(
2
1
)
(
(4.44)
The function
)
(t
f
can be reconstructed from
)
,
( b
a
W
by using the double
integral representation
db
a
da
a
b
t
b
a
W
a
C
t
f
³ ³
f
f
f
f
\
¸
¹
·
¨
©
§
\
S
)
,
(
1
1
)
(
2
(4.45)
where the scale parameter
a
is restricted to positive values only.
Four different wavelet bases functions can be used in (4.42); practical details
in applying wavelet analysis are taken from [30-31]. The use of this special
bases functions into EOF analysis prompt further investigation of the origin of
mechanism generating nonlinear behavior [32].
4.6.2 Empirical Mode Decomposition-Based EOF Analysis
The empirical mode decomposition (EMD) technique is a systematic method for
numerically decomposing a time dependent signal
)
(t
g
into its own intrinsic
mode functions (IMFs). This decomposition is expressed in the form
92
)
(
)
(
)
(
1
1
t
r
t
c
t
g
n
n
j
j
¦
(4.46)
where the components
}
,...,
1
),
(
{
n
j
t
c
j
represent the oscillatory functions and
)
(t
r
n
is the residue, extracted through a process called sifting. The dominants
IMFs represent the oscillatory behavior of interest directly related to physical
characteristics, while the rest of the terms are associated with featureless, often
non-oscillatory characteristics. A nonlinear mechanism can give rise to the
presence of different time-scales in a signal, which would then be separated as
different IMFs. The next key idea is that each IMF has its very own time-scale of
variations, with oscillations that are symmetric about the local zero mean. One
of the biggest benefits of the EMD method is that it can be used as a novel filter.
[30-31] provide details of the empirical mode decomposition.
Traditionally, filtering is carried out in frequency space only, but there is a
great difficulty in applying the frequency filtering when the data are either
nonlinear, nonstationary, or both, and in general, harmonics are generated at all
frequencies. By performing the EMD decomposition, the different coexisting
modes of oscillation in the signal are separated. The advantage of this time-
frequency filtering is that the results preserve the full nonlinearity and
nonstationarity.
In the multiscale analysis proposed, the empirical orthogonal functions of
the coefficients computed using wavelet or EMD analysis is applied at each
time-scale and then combined the results at relevant scales. The approach is
referred to as multiscale WEOF analysis and involves three main steps:
Identifying scales where significant events are detected.
1. Computing the wavelet or EMD decomposition for each column (time
series) in the data matrix,
2. Computing the covariance matrix of the wavelet or EMD coefficients for
each scale of interest, and
3. Extracting dynamic features from the selected scales.
93
This methodology has the advantage of being able of analyzing the time
series in individual manner; its near-real-time application is limited for the
multiscale decomposition of the time series. The multiscale decomposition in
real-time of time series already has been treated in wavelet and EMD analysis
[35-36].
4.7 Conclusions
Wide-area real-time monitoring prove to be invaluable in power system
dynamic studies by giving a quick assessment of the damping and frequency
content of dominant system modes after a critical contingency.
In this chapter, an alternative technique based on time-dependent complex
EOF analysis of measured data is proposed to resolve the localized nature of
transient processes and extract dominant temporal and spatial information. It is
shown that, in addition to providing spatial and temporal information, the
method improves the ability of conventional correlation analysis to capture
temporal events. Beside, the method gives a quantitative result for both the
amplitude and phase of motions, which are essential in the interpretation and
characterization of transient processes in power systems.
By combining a recursive method for identification of nonstationary effects
in transient processes based on a sliding window method and complex EOF
analysis, a novel framework, for on-line characterization of temporal behavior,
is proposed and attempts to consider the influence of both, spatial and temporal
variability. The basic idea of the recursive EOF is that the instantaneous
fluctuation of a harmonic wave form is measured by moving the window at
each sampling period. This recursive algorithm benefits the real-time
computation at the high-speed sampling.
94
4.8 References
[1] K. Karhunen, "Spectral theory of Stochastic processes," Mathematica-Physica, pp. 1-
7, 1946.
[2] D. G. Manolakis, V. K. Ingle, and S. M. Kogon, "Statistical and adaptive signal
processing: spectral estimation, signal modelling, adaptive signal filtering and
array processing," Norwood, Massachusset: Artech House, Inc., 2005.
[3] A. Papoulis, "Probability, random variables, and stochastic processes," New York:
McGraw-Hill
, 1984.
[4] F. Wang, and R. Scott, "On the prediction of linear stochastic systems with a low-
order model," Tellus, vol. 1, no. 57A, pp. 12-20, 2005.
[5] P. A. Parrilo, S. Lall, and F. Paganini, "Model reduction for analysis of cascading
failures in power systems," Proceedings of the American Control Conference, pp. 4208-
4212, 1999.
[6] S. S. Ravindran, "Reduced-order controllers for control of flow past an airfoil,"
International Journal for Numerical Methods in Fluid
, vol. 50, pp. 531-554, 2006.
[7] M. B. Galin, "Study of the dynamics of ultralong waves using time-dependent
empirical orthogonal functions," Izvestiya, Atmospheric and Oceanic Physics, vol. 38,
no. 1, pp. 27-39, 2002.
[8] B. F. Farrell, and P. J. Ioannou, "A theory for the statistical equilibrium energy
spectrum and heat flux produced by transient Baroclinic wave," Journal of the
Atmospheric Sciences
, vol. 51, no. 19, pp. 2685-2698, 1994.
[9] J. Terradas, R. Oliver, and J. L. Ballester, "Application of statistical techniques to
the analysis of solar coronal oscillations," The Astrophysical Journal, vol. 614, pp.
435-447, Oct., 2004.
[10] J. Q. Restrepo, "A real-time wide-area control for mitigating small-signal
instability in large electrical power systems," Electrical Engineering and Computer
Sciences
, May 2005.
[11] S. Toh, "Statistical model with localized structures describing the spatial-temporal
chaos of Kuramoto-Sivashinsky equation," Journal of the Physical Society of Japan,
vol. 56, no. 3, pp. 949-962, Mar., 1987.
[12] R. Avila-Rosales, and J. Giri, "Wide-area monitoring and control for power system
grid security," 15
th
[13] M. A. Merrifield, and R. T. Guza, "Notes and correspondence. Detecting
propagating signals with complex empirical orthogonal functions: a cautionary
note," Journal of physical oceanography, pp. 1628-1633, Oct., 1990.
PSCC, Liege, pp. 1-7, 2005.
[14] J. G. Proakis, and D. G. Manolakis, "Digital treatment of signals," Pearson Ed. 2007.
95
[15] V. Oppenheim, and R. W. Shafer, "Discrete-Time Signal Processing," 2
nd
[16] M. Li, and S. C. Lim, "Modeling network traffic using generalized Cauchy
process," Physica A, vol. 387, no. 11, pp. 2584-2594, 2008.
Edition,
Prentice Hall
, 1998.
[17] M. Li, "Generation of teletraffic of generalized Cauchy type," Physica Scripta, vol.
81, no. 2, pp. 1-10, 2010.
[18] M. Li, W.-S. Chen, and L. Han, "Correlation matching method of the weak
stationarity test of LRD traffic," Telecommunication Systems, vol. 43, no. 3-4, pp. 181-
195, 2010.
[19] P. Holmes, J. L. Lumley, and G. Berkooz, "Turbulence, coherent structures,
dynamical systems and symmetry," New York: Cambridge University Press, 1996.
[20] T. P. Barnett, "Interaction of the Monsoon and pacific trade wind system at
interannual time scales. Part 1: the equatorial zone," Monthly weather review, vol.
111, pp. 756-773, 1983.
[21] A. Hannachi, I. T. Jolliffe, and D. B. Stephenson, "Review. Empirical orthogonal
functions and related techniques in atmospheric science: a review," International
Journal of Climatology
, vol. 27, pp. 1119-1152, May 2007.
[22] H. Dankowicz, P. Holmes, G. Berkooz, and J. Elezgaray, "Local models of spatio-
temporally complex fields," Physical D, vol. 90, pp. 387-407, 1996.
[23] N. Aubry, R. Guyonnet, and R. Lima, "Spatiotemporal analysis of complex signals:
theory and applications," Journal of Statistical Physics, vol. 64, no. 3-4, pp. 683-739,
1991.
[24] Z. Gamm, I. I. Golub, L. S. Irkutsk, "Determination of locations for FACTS and
energy storage by the singular analysis," IEEE, pp. 411-414, 1998.
[25] E. Campbell, "A review of methods for statistical climate forecasting," Technical
Report
06/134, Aug., 2006.
[26] K.-Y. Kim, and G. R. North, "EOF-Based linear prediction algorithm: examples,"
Journal of Climate
, vol. 12, pp. 2076-2092, July 1999.
[27] S. Uhlig, "On the complexity of Internet traffic dynamics on its topology",
Telecommunication Systems
, vol. 43, no. 3/4, pp. 167-180, 2010.
[28] M. Li, W.-S. Chen, and L. Han, "Correlation matching method of the weak
stationarity test of LRD traffic", Telecommunication Systems, vol. 43, no. 3/4, pp. 181-
195, 2010.
[29] M. Li, "Generation of teletraffic of generalized Cauchy type", Physica Scripta, vol.
81, no. 2, pp. 1-10, 2010.
[30] M. Farge, "Wavelet transforms and their applications to turbulence," Annu. Rev.
Fluid Mech
., vol. 24, pp. 395-457, 1992.
96
[31] C. Torrence, and G. P. Compo, "A practical guide to Wavelet analysis," Bulletin of
the American Meteorological Society
, pp. 61-78, 1998.
[32] D. Ruiz-Vega, A. R. Messina, and G. Enriquez-Harper, "Analysis of inter-area
oscillations via non-linear time series analysis techniques," 15
th
[33] A. R. Messina, and V. Vittal, "Extraction of dynamic patterns from wide-area
measurements using empirical orthogonal functions," IEEE Transactions on Power
Systems
, vol. 22, no. 2, pp. 682-692, May 2007.
PSCC
, pp. 1-7,
Aug., 2005.
[34] P. Flandrin, G. Rilling, and P. Goncalves, "The empirical mode decomposition as a
filter bank," IEEE Signal Process. Lett., vol. 11, no. 2, pp. 112-114, Feb., 2004.
[35] A. R. Messina, V. Vittal, G. T. Heydt, and T. J. Browne, "Nonstationary approaches
to trend identification and denoising of measured power system oscillations,"
IEEE Transactions on Power Systems
, pp. 1-10, 2009.
[36] S.-W. Chen, H.-C. Chen, and H.-L. Chang, "A real-time QRS detection method
based on moving-average incorporating with wavelet denoising," Comput. Meth.
Prog. Biomed
., vol. 82, pp. 187-195, 2006.
97
Chapter 5
Application
This chapter examines the application of the proposed techniques to assess
swing oscillations patterns in power systems. Attention is focused on three
main aspects, namely the identification of critical modes and the associated
areas involved in the oscillations, the use of empirical orthogonal function
analysis to determine coherency, and the extraction of modal characteristics.
The evaluation of oscillation patterns in power systems is mainly done using
both, conventional EOF analysis and the complex EOF method. To facilitate
comparison between both simulation approaches and allow full comparison
between static and dynamic approaches, statistical techniques are developed to
extract modal characteristics directly from time domain simulations. The
approach can be used to analyze simulated and measured data.
Application studies on a large-scale power system model show that results
determined from the reduced-order representation are in qualitative and
quantitative agreement with those of conventional methods even in conditions
far from linearity. The results obtained from EOF analysis are compared with
conventional time domain analysis techniques for verification purposes.
It is shown that the method can predict the spatial location of damping with
good accuracy, and also gives some indication of the correct mechanism of
damping.
98
5.1 Complex Empirical Orthogonal Function Analysis of
Power-System Oscillatory Dynamics
5.1.1 Application to Time Synchronized Measured Data
To test the ability of the method to analyze complex oscillations, we analyze
data from time-synchronized measurements. The data used for this study were
recorded by phasor measurement units (PMUs) over a 400 ms window during a
real event in northern México. A brief description of the data follows. More
detailed information on system measurements can be found in [1,2].
At local time 06:27:42 in the early hours of 1 January 2004, undamped
oscillations involving frequency, voltage and power were observed throughout
the northern systems of the Mexican interconnected system (MIS). The main
event that originated the oscillations was a failed temporary interconnection of
the Northwest regional systems (Hermosillo (H) and Mazatlan Dos (MZD)
substations) to the MIS through a 230 kV line between Mazatlan Dos and Tres
Estrellas (TTE) substations. It is noted that, prior to this oscillation incident, the
northwestern system operated as an electrical island.
Oscillations in the northern systems with periods about 0.61 Hz, 0.50 Hz and
0.27 Hz persisted for approximately 120 ms before the northwestern system was
disconnected from the Mexican interconnected system. During the time interval
06:27:42-06:28:54 the system experienced severe fluctuations in frequency,
power, and voltage resulting in the operation of protective equipment with the
subsequent disconnection of load, independent generation and major
transmission resources.
Figure 5.1 shows a geographical diagram of the MIS showing the PMU
locations and the location of the initiating event. For demonstration purposes,
three buses spread across the system are selected. To allow comparison with
previous work, bus frequency signals are used in the analysis.
Figure 5.2 is an extract from PMU measurements of this event showing the
observed oscillations of selected bus frequencies.
99
For the purpose of further comparison to EOF analysis, the relevant time
interval of concern is zoomed in. System measurements in this plot demonstrate
significant variability suggesting a nonstationary process in both space and
time.
The most prominent variations occur in the interval during which the
oscillation start at 06:27:42 and the interval in which the operating frequency is
restored to the nominal condition (60 Hz) by control actions (06:28:21).
As discussed later, during this period the system experiences changes in
frequency (amplitude) content and mode shapes.
Fig. 5. 1. Schematic of the MIS system showing the location of the observed oscillations.
Measurement locations are indicated by shaded circles.
100
Fig. 5. 2. Time traces of recorded bus frequency swings recorded January 1, 2004 and detail of
the oscillation build up.
As a first step towards the development of a POD basis, the observed
records are placed in a complex data matrix,
)
(t
x
X
, as
»
»
»
¼
º
«
«
«
¬
ª
)
(
)
(
)
(
)
(
)]
(
),
(
),
(
[
)
(
1
1
1
1
N
x
N
x
x
x
TTE
MZD
H
x
t
u
t
u
t
u
t
u
t
t
t
t
j
j
#
%
#
u
u
u
where
3
,
2
,
1
j
is the spatial position index (grid location),
t
is time, and N is
the number of data points in the time series. For our simulations, 2021
snapshots are available, representing equally spaced measurements at three
different geographical locations. Each time series is then augmented with an
imaginary component to provide phase information and the EOF method is
employed to approximate the original data by a general, nonstationary, and
correlated frequency model.
101
The following subsections describe the application of complex empirical
orthogonal function analysis to examine the temporal and spatial variability of
measured data.
5.1.2 Construction of POD Modes via the Method of Snapshots
The method of snapshots was applied to derive a spatiotemporal model of the
oscillations. In order to improve the ability of the method to capture temporal
behavior, the individual time series are separated into their time-varying mean
and fluctuating components. By separating the data into their mean and
fluctuating components, EOF analysis is able to selectively determine the
temporal behavior of interest. The method may also help reduce the detrimental
effects of crossing events, spatial aliasing and random contamination.
In the application of the proposed method to measured data, we assume
that each signal can be construed as a superposition of fast oscillations on top of
a slow oscillation (the time-varying instantaneous mean) [3]. The slow
oscillation essentially captures the nonlinear (and possibly time-varying) trend
while the slow oscillations are the fluctuating parts.
A two-stage analysis technique based on wavelet shrinkage is proposed to
determine the temporal properties of time-synchronized information. In the
first step, the original system time histories are decomposed into their time-
varying mean speeds and fluctuating speeds through wavelet shrinkage. More
formally, the recorded time series are decomposed into their time-varying mean
frequencies,
)
(
)
(
k
x
m
t
u
j
, and nonstationary fluctuating components,
)
(
k
x
t
u
j
as
follows:
)
(
)
(
)
(
)
(
k
x
k
x
m
k
x
t
u
t
u
t
u
j
j
j
(5.1)
In the second stage, complex-EOF analysis is applied to the fluctuating field
to decompose the spatio-temporal data into orthogonal temporal modes.
Following Donoho [4], the two-stage approach can be used to reconstruct an
unknown function f from noisy data
i
i
i
z
t
f
d
V
)
(
, i=0,...,n-1, where
i
d
is the
observed data point (the noise contaminated measurement point),
i
z
is a
102
standard Gaussian white noise,
i
t
are equispaced data points, and
V
is a noise
level. The first stage is to find the estimates
(.)
)
,
(
(.)
^
dy
y
T
f
, where
(.)
)
,
(
dy
y
T
is a suitable reconstruction formula with spatial smoothing parameters
G
, and
dy
is data-adaptive choice of the spatial smoothing parameter . The separating
procedure is carried out using the wavelet shrinkage method because of its
ability to capture fast changes in the data.
The adopted separation/identification procedure can be summarized as
follows:
1. Expand the data
)
(
k
x
t
u
j
into wavelet series. Using the wavelet
decomposition structure, estimate level dependent thresholds for signal
compression using suitable thresholding approach
2. Obtain a denoised compressed function,
)
(
j
x
m
u
, representing the
instantaneous mean, via the wavelet shrinkage, and
3. Compute the fluctuation of the signal by calculating
)
(
)
(
)
(
)
(
k
x
m
k
x
k
x
t
u
t
u
t
u
j
j
j
, where
)
(
j
x
m
u
is the time-varying mean of the signal. Among the existing
techniques, we use the Birge-Massart strategy in [5].
The algorithm is simple to implement and the computational requirements
are small.
5.2 Spatiotemporal Analysis of Measured Data
Based on the analytical procedure outlined in section 3.5, the complex-EOF
method was used to determined dynamic trends and to analyze phase
relationships. In this procedure, the measured data are augmented with an
imaginary component defined by its Hilbert transform, and the temporal
patterns are extracted by using the procedure in section 3.4.
Complex EOF analysis was performed on the original time series, and the
nonlinearly detrended (wavelet shrinkage model) time series. The analysis
presented here uses the Daubechies wavelet with a fixed decomposition level.
103
A level-dependent threshold is then obtained using a wavelet coefficients
selection rule based on the Birge-Massart strategy [5].
As seen in Fig. 5.3, this model very effectively describes the long-term
behavior of the data while also capturing transient fluctuations. This, in turn,
results in improved characterization of system behavior.
Spectral analysis results for the leading POM in Fig. 5.4 show that the main
power is concentrated in oscillations with frequencies about 0.61, 0.50 and 0.27
Hz, which can be associated with major inter-area modes in the system [2].
Small peaks in Fig. 5.4 may indicate nonlinear interactions between frequency
components. From Fig. 5.4 it is also seen that the time-varying instantaneous
mean (bold gray line) acts as a high pass-filter. This approach can be used to
separate the slow from the fast components as well as to improve the numerical
accuracy of the method. This is discussed further in the next two subsections.
Figure 5.5 presents the corresponding eigenspectrum of these modes
computed to capture 99.9% of the signal's energy. In these plots, the horizontal
axis shows the number of modes required to attain 99% of the averaged total
energy; the vertical axis shows the energy captured in (3.34) by each POM. The
data set being analyzed corresponds to bus frequency measurements observed
with the PMUs at various geographical locations.
Complex EOF analysis of the synchronized measurements of frequency in
Fig. 5.5 shows that wide-area system dynamics is well represented by three
modes; the two leading modes together account for 96.5 % of the total energy.
Individually, these modes account for 72, 24.5 and 3.5 % of the energy (see Fig.
5.5 caption for details).
104
Fig. 5. 3. Time-Varying means and fluctuating components of the recorded bus frequency
signals. The time series are smoothed using wavelet shrinkage.
Fig. 5. 4. Comparison of the Fourier transform spectrum of the original signal and the spectra
constructed with the instantaneous mean removed.
Fig. 5. 5. Energy captured as a function of the number of modes. The percentage of energy
located in the jth mode is measured by
¦
n
i
i
i
E
1
100
.
105
Figure 5.6 compares the reduced solution using the POD basis functions to
the full solution obtained from the measured data for the Hermosillo signal in
Fig. 5.3. As observed in this plot, using only three modes we are able to
accurately approximate the measured data over the entire observation period.
The agreement between the reduced order model and the observed behavior
illustrates the high degree of accuracy that is possible with a simplified model.
Fig. 5. 6. Reconstruction of the original data using the three leading POMs. Solid line represent
the original time series and dotted lines represent the composite oscillation obtained by adding
the temporal modes.
On the basis of these results, we conducted detailed analysis aiming at
disclosing hidden information in the data. For clarity of exposition, the analysis
of temporal and spatial patterns will be presented separately.
5.2.1 Temporal Properties
In order to reveal the hidden wave signatures in the time series, we examine the
spatial, temporal and frequency patterns of system behavior in the light of
complex orthogonal function analysis. Figure 5.7 illustrates the temporal
evolution (amplitude and phase) of the dominant mode.
These functions display a number of interesting features. As discussed
below, the analysis identifies two periods of interest; a transient period
associated with the interconnection of the systems (06:27:42-06:28:21), and a
106
nearly stationary period in which the frequency of the interconnected system is
restored to its normal value (06:28:54-06:29:39).
The first interval manifests particularly strong temporal activity as can be
seen in Fig. 5.7a. In interpreting these results, we remark that break or changes
in the temporal functions may signal different physical regimes or control
actions.
An examination of the temporal phase in Fig. 5.7b, on the other hand,
reveals a nonstationary behavior in which the phase (frequency) content
changes with time. Here, the slope of the spatial phase function represents the
instantaneous frequency. The slowly increasing trends indicate periods of
essentially constant frequencies.
a) Temporal amplitude
b)
Temporal phase
Fig. 5. 7. Temporal patterns of variability associated with the dominant mode.
107
5.2.2 Frequency Determination from Instantaneous Phases
Additional insight into the frequency variability of the observed oscillations can
be obtained from the analysis of instantaneous frequencies. Recognizing that
the instantaneous frequency is the time derivative of the temporal phase
function,
T
, the instantaneous frequencies can be estimated from (3.28) for each
mode of concern.
The study focuses on POM 1 which is the mode that captures most of the
variability in the signal. Figure 5.8 gives the instantaneous frequency of POM1
for the interval of interest in this study. Also plotted, is the instantaneous mean
frequency (nonlinear trend) determined using wavelet shrinkage analysis
above.
Nonstationary features are evident in this plot. Analysis of Fig. 5.8 shows
two modal components: a 0.27 Hz component associated with the steady-state
behavior of the system, and a 0.64 Hz component associated with the transient
system fluctuation following the system interconnection. The 0.27 Hz
component captures the slow ambient swings previous to the onset of system
oscillations and the steady behavior of the system. The results are consistent
with those based on nonstationary analysis of the observed oscillations giving
validity to the results.
Fig. 5. 8. Instantaneous frequency of POM 1. Time interval 06:27:42-06:29:39.
108
5.2.3 Mode Shape Estimation
One of the most attractive features of proposed technique is its ability to detect
changes in the shape properties of critical modes arising from topology changes
and control actions. Changes in mode shape may indicate changes in topology
or changes in load/generation and may be useful for control decisions and the
design of special protection systems. This is a problem that has been recently
addressed using spectral correlation analysis [6].
Using the spatial phase and amplitude, the phase relationship between key
system locations (the mode shapes) can be determined. In this analysis, we
display the complex value as a vector with the length of its arrow proportional
to eigenvector magnitude and direction equal to the eigenvector phase.
Figure 5.9 shows the mode shape for the three intervals of interest above
(06:27:42-06:28:06, 06:28:06-06:28:21 and 06:28:21-06:28:54) computed using the
spatial function (3.15). It is interesting to note that the dominant POM mode
shape changes with time. The effect is more pronounced for the time interval
06:28:06-06:28:21 in which several control actions take place in the system. This
information may be useful to identify the dominant generators involved in the
oscillations, and ultimately devise control mechanisms to damp the observed
oscillations.
These results are in general agreement with previously published results
based on real EOF analysis and Prony results [2]. The new results, however,
provide clarification on the exact phase relationships between key system
measurements as a function of time.
109
Fig. 5. 9. Mode shape of POM 1 for various time intervals of interest.
5.2.4 Energy Distribution
In the previous section it was shown that a linear combination of an individual
eigenmode can accurately reconstruct the temporal behavior of simultaneous
measurements at different geographical (spatial) locations. A key related
question of interest is that of finding a small number of measurements that will
provide a good estimation of the entire field of interest.
Based on the decomposed EOFs, complex EOF analysis was used to
determine the locations with the most energy. Figure 5.10 shows the
participation of each location to the total energy of the record. The x-axis shows
spatial sensor location and the y-axis shows the energy value. From Fig. 5.10, it
is evident that modes 1 and 2 are quite prominent at the Mazatlan Dos and
Hermosillo substations while mode 3 is more strongly evident at the Tres
Estrellas substation. This is consistent with conventional analysis (not shown).
However, the proposed approach provides an automated way to estimate mode
shapes without any prior information of the time intervals of interest.
110
Fig. 5. 10. Energy distribution.
The analysis of propagating features using the data-driven framework that
integrates the use of EOF analysis and the sliding window-based method
proposed in Chapter 4 is given in the next section to analyze the recorded data
from the real event in northern México recollected by phasor measurement
units.
5.3 Extraction in near-real-time of propagating Features of
Power-System Oscillatory Dynamics
To test the ability of the proposed method given in Chapter 4 to isolate the
portion of data where oscillation are present and improves numerical efficiency
and accuracy in dynamical characterization of oscillations in power systems, we
analyze time-synchronized data recorded by phasor measurement units during
a real event in northern México. The records of selected signals are given in Fig.
5.2, this with the objective of compare the results obtained in section 5.2 with
the proposed recursive formulation in Chapter 4 to the estimation of
instantaneous parameters.
Using the algorithm proposed in subsection 4.5.5 to the estimation of
instantaneous parameters, in Fig. 5.11 is shown the fluctuating components of
111
the individual time series separated into their time-varying mean. The time-
varying mean to the selected signal was computed using (4.39) for a fixed
window size of 2 samples. As is observed, in this plot demonstrate significant
variability suggesting a nonstationary process; the most prominent variations
occur in the interval during which the oscillations starts at 06:27:42 and the
interval in which the operating frequency is restored to the nominal condition
(60 Hz) by control actions at about 06:28:54.
Fig. 5. 11. Time trace of bus frequency swings recorded on January 1, 2004 and detail of the
oscillation build up.
Based on the frequency data collected from the fluctuating components
given in Fig. 5.11, each time series is augmented with an imaginary component
using Hilbert analysis. Therefore, a near-real-time complex EOF analysis was
applied to reveal the ability of the proposed method to isolate the portion of
112
data where oscillations are present and improves numerical efficiency and
accuracy.
Using the proposed approach, two statistical modes describing standing and
traveling features were identified. Mode 1, which accounts for 98.5 % of the
total energy describes the main (standing) feature in the records. The second
mode contains approximately 1 % of the energy and describes propagating
features.
Figure 5.12 compares the signal reconstructed using the conventional (off-
line) EOF and the proposed recursive (on-line) formulation for a fixed window
size of 2 samples, while Table 1 gives root mean square
)
(rms
relative error
using different window size. For simplicity, only the traveling wave mode is
used in the computation.
Simulation results clearly show that near-real-time approximation enhance
the ability of the method to capture local information.
From Table 5.1 it can be seen that large values of
W
obscure the analysis of
temporal changes occurring at segments of the observed record. The
rms
relative error is computed of the form
³
³
dt
f
dt
f
f
rms
2
2
^
(5.2)
where f is the real scalar field and
f^
represents the estimated value.
113
Fig. 5. 12. Reconstruction of the original data with the traveling wave mode identified using the
conventional (off-line) and recursive (on-line,
2
W
) EOF analysis.
Table 5. 1 RMS error in EOF approximation.
Signal
Conventional (rms)
Recursive (rms)
k=1,2,...,2002
Hermosillo
0.3591
0.3458
0.3406
0.2614
Mazatlan Dos
0.4928
0.4923
0.4871
0.4008
Tres Estrellas
0.7536
0.5620
0.5575
0.5477
Also of interest, Fig. 5.13 compares the amplitude of modal components
provided by conventional EOF analysis and that obtained from the proposed
formulation. Note that conventional analysis can only provide average
information. In sharp contrast with this, the near-real-time implementation
enables to single out the times at which sudden changes in system behavior
take place. This can be confirmed through the instantaneous frequency in Fig.
5.14, where the nonlinear trend was computed from the recursive EOF
formulation using a 5-point moving average.
114
Fig. 5. 13. (Upper) average energy of travelling wave mode using the conventional (off-line)
EOF analysis, and (lower) instantaneous energy of travelling wave mode using the recursive
(on-line,
2
W
) EOF analysis.
Fig. 5. 14. Instantaneous frequency of dominant mode.
One of the most attractive features of proposed technique is its ability to
detect changes in the shape properties of critical modes arising from topology
115
changes [1], and control actions. Figure 5.15 shows the mode shape of dominant
electromechanical modes obtained using both approaches. For simplicity
attention is constrained to the time interval 133 to 135 ms. Changes in the mode
shape of electromechanical modes may indicate changes in topology or changes
in operating patterns and may be useful for control decisions and the design of
special protection systems [1,6]. In this analysis, complex values are displayed
as vectors with the length of the vector proportional to the eigenvector
magnitude and the direction equal to the eigenvector phase. A phase of 0° is
represented as a vector directed upward.
The proposed approach provides an automated way to estimated mode
shape without any prior information of the time intervals of interest.
Fig. 5. 15. (Upper) mode shape of travelling wave mode using the conventional (off-line) EOF
analysis, and (lower) time-varying mode shape of travelling wave mode using the recursive
(on-line,
2
W
) EOF analysis.
116
Finally, Fig. 5.16 shows the time evolution of the dominant travelling wave
at different system locations obtained using the wavelet empirical orthogonal
function analysis described in section 4.6. For all simulations, the Morlet
wavelet mother function has been used. The analysis depicts the fill evolution
of temporal dynamics and provides details about the local behavior of
multiscale data. Here, the shaded areas in the plot delineate the regions where
the energy is significant. This provides detailed information of time instants
associated with transient system behavior whose behavior changes over time
and frequency. The results stand in contrast with current modeling done at a
single scale.
Fig. 5. 16. Spatial patterns of the leading mode showing temporal variability.
5.4 Discussion
Wide-area, real-time monitoring may prove invaluable in power system
dynamic studies by giving a quick assessment of the damping and frequency
content of dominant system modes after a critical contingency. In this chapter,
an alternative technique based on time-dependent complex EOF analysis of
117
measured data is used to resolve the localized nature of transient processes and
to extract dominant temporal and spatial information.
A method of spatially decomposing oscillation patterns in near-real-time
into their standing and travelling parts is presented. By combining a sliding
window approach with complex EOF, a novel framework for on-line
characterization of temporal behavior is proposed that attempts to consider the
influence of both, spatial and temporal variability.
Numerical results show that the proposed method provides accurate
estimation of nonstationary effects, modal frequency, time-varying mode
shapes, and time instants of intermittent transient responses. This information
may be important in determining strategies for wide-area control and special
protection systems.
It is shown that, in addition to providing spatial and temporal information,
the method improves the ability of conventional correlation analysis to capture
temporal events and gives a quantitative result for both the amplitude and
phase of motions, which are essential in the interpretation and characterization
of transient processes in power systems.
Extensions to this approach to determine the optimal locations to place
PMUs are underway. Applications to real-time system monitoring, protection
and control will be addressed in future work.
5.5 References
[1] A. R. Messina, M. A. Andrade, and E. Barocio, "Wide-area monitoring and
analysis of inter-area oscillations using the Hlbert-Huang transform," pp. 285-316,
Book Chapter in: "Leading-Edge Electric Power System Research," Cian M.
O'Sullivan (editor)
, Nova Science Publishers, Inc., 2008. New York.
[2] A. R. Messina, and V. Vittal, "Extraction of Dynamic Patterns from Wide-Area
Measurements using Empirical Orthogonal functions," IEEE Transactions on Power
Systems
, vol. 22, no. 2, pp. 682-692, May 2007.
[3] L. Chen, and C. W. Letchford, "Proper orthogonal decomposition of two vertical
profiles of full-scale nonlinearity downburst wind speeds [lzcl]," Journal of Wind
Engineering and Industrial Aerodynamics
, vol. 93, pp. 187-216, 2005.
118
[4] D. L. Donoho, and I. M. Johnstone, "Ideal spatial adaptation by wavelet
shrinkage," Biometrical Trust, vol. 81, no. 3, pp. 425-455, August 1994.
[5] M. Misiti, Y. Misiti, G. Oppenheim, and J. Poggi, "Wavelet Toolbox for use with
Matlab," The MathWorks Inc., Natick, Massachusetts, 1996.
[6] D. J. Trudnowski, "Estimating Electromechanical Mode Shape from
Synchrophasor measurements," IEEE Transactions on Power Systems, vol. 23, no. 3,
pp. 1188-1195, August 2008.
119
Chapter 6
Conclusions and
Recommendations for Future Work
6.1 General Conclusions
Wide-area monitoring schemes offer the possibility to analyze and characterize
complicated oscillations involving the exchange of kinetic energy between
machines or electrical areas in large interconnected power systems. In this
research, time-synchronized PMU measurements collected at various systems
locations are used here to propose a multi-variate statistical data analysis
technique with the ability to capture accurately and efficiently independent
spatio-temporal patterns from time synchronized data.
The work developed in this thesis, presents a new approach to dynamic
characterization of spatial and temporal changes in wide-area measurements of
system oscillatory based on complex empirical orthogonal function (EOF)
analysis. The method improves the ability of conventional EOF analysis in
interpretation and characterization of nonstationary effects, modal frequency,
time-varying mode shapes and time instant of intermittent or irregular transient
behavior associated with abrupt changes in system topology or operating
conditions.
This approach enables the systematic computation of dynamic trends and
characteristics of simultaneously observed oscillations and allows the practical
application of the formulation to a wide variety of phenomena. Further, because
of its simplicity and generality, the proposed method has several advantages
over other methods that make it suitable for analyzing complex phenomena.
120
The modal characteristics of the response matrix are estimated using an
alternative method based on singular value decomposition that improves
efficiency and results in a more efficient representation of the data. The ability
to compress the variability of large data sets into the lower possible number of
statistical modes makes the technique adequate for the study of multivariate
data.
The use of conventional EOF analysis in combination with Hilbert analysis
provides an efficient mean to compute standing and propagating features of
general nonstationary processes and enables to identify important information
of interest in the analysis of power system dynamic phenomena.
A sliding window-based approach has been combined with EOF analysis to
systematically analyze the observational data in near-real-time and resolve
localized information; this allows isolate the portions of data where oscillations
are present and improves numerical efficiency and accuracy to the proposed
method. Furthermore, an interesting and useful interpretation of the covariance
matrix is introduces that allows extracting modal information using recursive
relations. A rigorous statistical formulation that leads to a recursive estimation
of the covariance matrices is also given.
Other main conclusions obtained from this work can be summarized as
follows:
¾ The conceptual framework developed provides a rigorous basis for the
analysis, detection and simplification of complex oscillations and enables the
simultaneous study of synchronized measurements.
¾ The use of data-driven proper orthogonal decomposition makes it is
possible to identify the amplitude and phase relationships of general
arbitrary processes. Unlike existing approaches, the bases of the linear
approximations are given by the data itself.
¾ The advantage of this technique is that is easy to implement and is based on
the covariance structure of the data.
121
¾ The method needs only measurements from the test system and can be
applied to systems with different types of oscillations.
¾ Experience with the analysis of a complex system shows that the method
can be efficiently (and accurately) applied to the analysis of wide-area
phenomena involving a large number of machines and time scales.
¾ Complex empirical orthogonal function analysis shows promise for the
study of spatio/temporal behavior rapidly varying oscillations. This is an
aspect of the thesis that deserves further investigation.
¾ A truncation criterion is proposed for selecting the most significant
oscillatory modes in base to its kinetic energy (variance). The key advantage
of this criterion is that allows extracting information from noisy time series
without prior knowledge of the dynamics affecting the time series.
6.2 Recommendations for Future Work
The recommendations for future work can be summarized as follows
1. The extension of the proposed method to the frequency domain, where a
complex representation is also necessary. The imaginary part is computed
from the time series using the fast Fourier transform (FFT) and to estimate
the real part it is necessary develop some method. The advantage that is
hoped with this future work is the study of transient phenomena sampled in
different time step and frequency scales. Furthermore, extracting
information from noisy time series without prior knowledge of the
dynamics affecting the time series.
2. It is well known that the covariance matrix C , can often result in an ill-
conditioned matrix due the large-scale wide-area analysis. Column scaling
which balance the Euclidean norm of the columns of C to unity can be
applied to improve its condition number. The column scaling can be
formally expressed by introducing a diagonal matrix D in the eigenvalues
122
decomposition, where the analytical development is considered as a
suggestion for future work.
3. The near-real-time spatio-temporal characterization from time series
synchronized data estimated using the recursive relation of proposed
method may be important in determining strategies for control and special
protection systems. The extensions of the proposed technique to design
system controllers like future work prove to be invaluable in wide-area real-
time monitoring.
4. A theory on how complex oscillations are propagated in geographic
coordinate system is necessary. It is thought that complex EOF analysis may
provide insight into the nature of this behavior. Again, this is an aspect of
the work that requires further work.
5. The use of many time synchronized phasor measurement units to wide-area
monitoring over space and time can be time-consuming and difficult
because of the large amount of information that needs to be processed and
compared, the use of block monitoring can be applied as an alternative
method to characterize the variability in wide-area. These blocks can be
analyzed in synchronized manner; this might reflect the real-time
computation at the high-speed sampling.
6. The method is particularly promising, and may be expanded into other
areas of analysis including the study of voltage-related phenomena, heat
transfer, fluid flow, propagating wave in tsunamis, structural dynamics, etc.
The suggestions for future work is now extend into other areas.
123
Appendix A
EOFs in the Presence of
Measurement Error
In this appendix, the EOF decomposition is extended to the analysis of temporal
behavior in the presence of measurement error. Assume that
)
,
( t
x
u
denote a
sequence of observations on some domain
:
x
, where
x
is a vector of spatial
variables and
t
is the time at which the observations are made. The approach
for
)
,
( t
x
u
can be computed as the sum of a smooth signal
)
,
( t
x
y
and additive
noise
)
,
( t
x
h
as [1-3]
)
,
(
)
,
(
)
,
(
t
x
h
t
x
y
t
x
u
(A.1)
Taking the variance [.]
E
, of
)
,
( t
x
u
gives
)
)]
,
(
)
,
(
([
2
t
x
u
t
x
u
E
m
h
y
u
/
/
/
(A.2)
where
)
,
( t
x
u
m
represent the mean of
)
,
( t
x
u
. Then,
)
)]
,
(
)
,
(
([
2
t
x
y
t
x
y
E
m
y
/
)
)]
,
(
)
,
(
([
2
t
x
h
t
x
h
E
m
h
/
and
)
(
)
(
x
x
H
u
h
y
u
M
M
/
C
C
C
(A.3)
where
u
C
,
y
C ,
h
C
are the covariance matrices defined to
)
,
( t
x
u
,
)
,
( t
x
y
,
)
,
( t
x
h
respectively, and
)
(x
M
represent an optimal basis for the modal decomposition
124
of the ensemble of data; the
)
(x
j
M
are the eigenfunctions of the linear equation
OM
M
C
.
If we let that the added measurement error for
)
,
( t
x
u
be white noise with
variance
2
h
V , then
)
(
)
)(
(
2
x
I
x
H
h
u
h
u
y
M
V
M
/
C
C
C
(A.4)
where
n
n
I
u
is the identities matrix and the letter ( H ) denotes the conjugate
transpose or Hermitian matrix. The j-th original time series is expressed as the
sum of the mode functions [4]
¦
n
j
j
j
x
t
a
t
x
y
1
)
(
)
(
)
,
(
M
(A.5)
where the time series are given as
)
(
)
,
(
)
(
x
t
x
y
t
a
H
j
M
(A.6)
Hence,
I
t
a
t
a
E
h
H
2
)]
(
)
(
[
V
/
(A.7)
where the diagonal matrix,
/ , contain the eigenvalues
}
,
,
,
{
2
1
n
O
O
O
.
Furthermore
T
n
t
a
t
a
t
a
)]
(
,
),
(
[
)
(
1
Thus from (A.4) and (A.7), the noisy
)
,
( t
x
u
process and the smooth
)
,
( t
x
y
have the same eigenfunctions,
)}
(
{
x
j
M
, but different eigenvalues (i.e.,
}
{
j
O and
}
{
2
h
j
V
O
, respectively). In addition, the time series amplitude for both
processes has the same mean (zero), but the variance of the
)
,
( t
x
y
process
amplitude time series is less (as expected). Thus, if white measurement error is
present (and it always is with observational data), then the spatial patterns (i.e.,
the
)}
(
{
x
j
M
) are not affected, but the variance explained by each pattern must
adjusted accordingly. This is seldom (if ever) done in practice.
125
A.1. References
[1] C. Kim Wikle, "Spatio-temporal statistical models with application to atmospheric
processes," Statistics, Geological and Atmospheric Sciences, 1996.
[2] H. Jazwinski, "Stochastic processes and filtering theory," Dover Edition, 2007.
[3] J.-P. Chiles, Geostatistics: Modeling Spatial Uncertainty, New York: Wiley, 1999.
[4] P. Holmes, J. L. Lumley, and G. Berkooz, "Turbulence, coherent structures,
dynamical systems and symmetry," New York: Cambridge University Press, 1996.
126
Appendix B
Optimal Linear Predictor
Assume we have data at n locations
}
,
,
,
{
2
1
n
x
x
x
and that we are interested in
predicting
)
,
( t
x
y
from
)
,
( t
x
u
. The truncated EOF expansion in the discrete
framework is given by
n
m
to
x
t
a
t
x
u
m
j
j
j
¦
1
)
(
)
(
)
,
(
^
M
(B.1)
In practice, it is common to assume that this truncation can be used as a
prediction of the smooth process
)
,
( t
x
y
[1-3]. that is
)
,
(
^
)
,
(
t
x
u
t
x
y
{
(B.2)
It is natural to ask if (B.2) is an optimal predictor. To investigate the
optimality of (B.2), we consider the space-only process analogous to (A.1), this
is
)
,
(
)
,
(
)
,
(
t
x
h
t
x
y
t
x
u
(B.3)
where
)
,
( t
x
y
is a smooth signal and
)
,
( t
x
h
is additive noise.
If we assumed that
)
,
( t
x
u
contain Gaussian additive noise and that we are
interested in predicting
)
,
( t
x
y
from
)
,
( t
x
u
, then the optimal predictor is just a
linear predictor. This optimal predictor is given by [1]
n
j
to
t
x
u
t
x
u
t
x
y
E
t
x
u
H
t
x
y
j
j
,...,
1
)
,
(
)]
,
(
|
)
,
(
[
1
)
,
(
)
,
(
C
C
(B.4)
127
where we have assumed that the means of
)
,
( t
x
u
and
)
,
(
t
x
y
j
are known to be
zero, and we define
)
)]
,
(
)
,
(
([
))
,
(
var(
)
))]
,
(
(
)
,
(
))][
,
(
(
)
,
(
([
))
,
(
),
,
(
cov(
2
)
,
(
)
,
(
t
x
u
t
x
u
E
t
x
u
t
x
y
E
t
x
y
t
x
y
E
t
x
y
E
t
x
y
t
x
y
m
t
x
u
T
j
j
j
t
x
y
j
C
C
(B.5)
Now, if we let
2
)]
,
(
)
,
(
[
h
k
j
t
x
h
t
x
h
E
V
for
k
j
, and equal to zero otherwise,
then we can write
2
)
,
(
)
,
(
)
,
(
)
,
(
h
t
x
y
t
x
h
t
x
y
t
x
u
V
C
C
C
C
(B.6)
now, since
)
,
( t
x
y
C
is real and symmetric, we can write the following
decomposition:
H
y
t
x
y
M
M
/
)
,
(
C
(B.7)
where the eigenvectors and eigenvalues are:
)
,...,
(
)
,...,
(
)
(
)
1
(
1
n
y
y
y
n
diag
O
O
M
M
M
/
(B.8)
then, using the orthogonality
I
H
H
M
M
M
M
, we obtain
H
u
H
h
H
y
t
x
u
M
M
M
M
V
M
M
/
/
2
)
,
(
C
(B.9)
where
I
h
y
u
2
V
/
/
(B.10)
Now, (B.4) can be written
n
j
to
t
x
u
t
x
u
t
x
y
E
H
u
H
t
x
y
j
j
,...,
1
)
,
(
)]
,
(
|
)
,
(
[
1
)
,
(
/
M
M
C
(B.11)
furthermore, since
¦
n
k
j
k
i
k
k
y
j
i
x
x
t
x
y
t
x
y
E
1
)
(
)
(
)
(
)]
,
(
|
)
,
(
[
M
M
O
(B.12)
then
128
)
(
)
(
)
,
(
j
y
t
x
y
x
x
j
M
M
/
C
(B.13)
where
n
j
to
x
x
x
T
j
n
j
j
,...,
1
)]
(
),...,
(
[
)
(
1
M
M
M
(B.14)
In addition, as for the EOF case [4]
¦
n
j
j
j
t
a
x
t
x
u
1
)
(
)
(
)
,
(
M
(B.15)
and
)
,
(
)
(
)
(
t
x
u
x
t
a
H
j
j
M
(B.16)
We can then write the optimal predictor as
n
j
to
x
t
a
t
a
x
t
x
u
t
x
y
E
n
k
j
k
k
h
k
y
k
y
H
H
y
j
j
,...,
1
)
(
)
(
)
(
)
(
)]
,
(
|
)
,
(
[
1
2
)
(
)
(
1
/
/
¦
M
V
O
O
M
M
M
M
M
(B.17)
Now, the truncated predictor of interest is given by
¦
m
j
j
j
x
t
a
t
x
u
1
)
(
)
(
)
,
(
^
M
(B.18)
Thus, the optimal predictor (B.17) is equivalent to the truncated predictor
(B.18) in the trivial case where
0
2
h
V
and
0
)
(
)
(
1
¦
n
m
l
l
l
x
t
a
M
. In general,
setting (B.17) equal to (B.18) gives
¦
¦
¦
m
k
j
k
k
n
m
l
j
k
l
h
l
y
l
y
m
k
j
k
k
h
k
y
k
y
x
t
a
x
t
a
x
t
a
1
1
2
)
(
)
(
1
2
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
M
M
V
O
O
M
V
O
O
(B.19)
which result in the equality constraint
n
j
to
x
t
a
x
t
a
n
m
l
j
k
l
h
l
y
l
y
m
k
j
k
k
h
k
y
k
y
,...,
1
)
(
)
(
)
(
)
(
1
1
2
)
(
)
(
1
2
)
(
)
(
¸
¸
¹
·
¨
¨
©
§
¦
¦
M
V
O
O
M
V
O
O
(B.20)
129
Clearly, (B.20) shows that if
0
2
h
V
, then optimality is achieved if
0
)
(
t
a
l
or
0
)
(
l
y
O
for
n
m
l
,
,
1
. Thus, in the absence of measurement error, if the
truncation parameter
m
is large so that the
)
(l
y
O are very close to zero, then the
truncated EOF smoothing technique is nearly optimal. In the presence of
measurement error, the discrete orthogonality of the eingevectors in (B.20)
implies optimality when
0
)
(
t
a
k
for
m
k
,
,
1
and
0
)
(
t
a
l
or
0
)
(
l
y
O
for
n
m
l
,
,
1
. Consequently, the smoothed quantity
)
,
( t
x
y
should be viewed as
a method for removing noise.
B. 1. References
[1] C. Kim Wikle, "Spatio-temporal statistical models with application to atmospheric
processes," Statistics, Geological and Atmospheric Sciences, 1996.
[2] H. Jazwinski, "Stochastic processes and filtering theory," Dover Edition, 2007.
[3] J.-P. Chiles, Geostatistics: Modeling Spatial Uncertainty, New York: Wiley, 1999.
[4] P. Holmes, J. L. Lumley, and G. Berkooz, "Turbulence, coherent structures,
dynamical systems and symmetry," New York: Cambridge University Press, 1996.
130
Appendix C
Correlation Function Using Linear
Instantaneous Estimates
Given
)
(
1
k
x
t
u
, we seek an estimate for
)
(
k
x
t
u
of the form
)
(
)
,
(
1
1
k
x
k
k
x
t
u
t
t
C
, by
requiring the function
)
,
(
1
k
k
x
t
t
C
to minimise the expression [1]
2
1
1
)
(
)
,
(
)
(
k
x
k
k
x
k
x
t
u
t
t
t
u
C
(C.1)
Following a similar procedure to that in section 2.2, a necessary condition is
that, for any
)
,
(
1
k
k
x
t
t
V
0
)
(
)]
,
(
)
,
(
[
)
(
0
2
1
1
1
G
G
G
k
x
k
k
x
k
k
x
k
x
t
u
t
t
t
t
t
u
d
d
V
C
(C.2)
Noting that the expression inside the averaging brackets is equal to:
)}
(
)]
,
(
)
,
(
[
)
(
)}{
(
)]
,
(
)
,
(
[
)
(
{
1
*
1
1
*
1
1
1
k
x
k
k
x
k
k
x
k
x
k
x
k
k
x
k
k
x
k
x
t
u
t
t
t
t
t
u
t
u
t
t
t
t
t
u
*
*
V
C
V
C
G
G
differentiating with respect to
G
, evaluating at
0
G
and equating to zero we
get:
]
)
(
)
(
)
,
(
)
,
(
Re[
2
]
)
(
)
(
)
,
(
Re[
2
1
*
1
1
1
1
*
1
k
x
k
x
k
k
*
x
k
k
*
x
k
x
k
x
k
k
*
x
t
u
t
u
t
t
t
t
t
u
t
u
t
t
C
V
V
(C.3)
Since
)
,
(
1
k
k
*
x
t
t
V
is an arbitrary variation this implies that
)
(
)
(
)
,
(
)
(
)
(
1
*
1
1
1
*
k
x
k
x
k
k
x
k
x
k
x
t
u
t
u
t
t
t
u
t
u
C
(C.4)
and therefore
131
)
(
)
(
)
(
)
(
)
,
(
1
*
1
1
*
1
k
x
k
x
k
x
k
x
k
k
x
t
u
t
u
t
u
t
u
t
t
C
(C.5)
The fluctuations typically present in a turbulent system make the
assumption that
)
(
)
(
1
*
1
k
x
k
x
t
u
t
u
is invertible at every point
1
k
t
a reasonable one,
so that (C.5) will be well defined [2-3].
C. 1. References
[1] P. Holmes, J. L. Lumley, and G. Berkooz, "Turbulence, coherent structures,
dynamical systems and symmetry," New York: Cambridge University Press, 1996.
[2] H. Jazwinski, "Stochastic processes and filtering theory," Dover Edition, 2007.
[3] J.-P. Chiles, Geostatistics: Modeling Spatial Uncertainty, New York: Wiley, 1999.
Excerpt out of 146 pages
- scroll top
- Quote paper
- Pedro Esquivel Prado (Author), 2010, Extraction of Dynamic Patterns from Wide-Area Measurements using Empirical Orthogonal Functions, Munich, GRIN Verlag, https://www.grin.com/document/280564
Look inside the ebook
-
-
-
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X. -
-
Upload your own papers! Earn money and win an iPhone X.