Research Article  Open Access
TimeDependent Statistical Analysis of WideArea TimeSynchronized Data
Abstract
Characterization of spatial and temporal changes in the dynamic patterns of a nonstationary process is a problem of great theoretical and practical importance. Online monitoring of largescale power systems by means of timesynchronized Phasor Measurement Units (PMUs) provides the opportunity to analyze and characterize intersystem oscillations. Widearea measurement sets, however, are often relatively large, and may contain phenomena with differing temporal scales. Extracting from these measurements the relevant dynamics is a difficult problem. As the number of observations of real events continues to increase, statistical techniques are needed to help identify relevant temporal dynamics from noise or random effects in measured data. In this paper, a statistically based, datadriven framework that integrates the use of waveletbased EOF analysis and a sliding windowbased method is proposed to identify and extract, in nearrealtime, dynamically independent spatiotemporal patterns from time synchronized data. The method deals with the information in space and time simultaneously, and allows direct tracking and characterization of the nonstationary timefrequency dynamics of oscillatory processes. The efficiency and accuracy of the developed procedures for extracting localized information of power system behavior from timesynchronized phasor measurements of a real event in Mexico is assessed.
1. Introduction
Phenomena observed in power system oscillatory dynamics are diverse and complex. Remotely sensed measured data are known to exhibit noisy, nonstationary fluctuations resulting primarily from small magnitude, random load changes in load, driven by lowscale motions or nonlinear trends originating from control actions or other changes in the system. Extracting from these sets indices that capture significant spatial and temporal dynamics is very challenging [1–4].
Recent improvements in widearea monitoring schemes have led to renewed investigation of nonlinear and nonstationary behavior of system oscillations [2, 3]. In the last decade, extensive research has been carried out for developing modal extraction algorithms that can accurately monitor transient response. Multivariate statistical data analysis techniques offer a powerful tool for analyzing power system response from measured data [1]. By seeing the snapshots of system data as a realization of random field generated by some kind of stochastic process, datadriven statistical can be applied to investigate propagating phenomena of different spatial scales and temporal frequencies [5–7]. Existing analysis methods, however, do not take into account the multiscale nature of measured system dynamics arising from events occurring at different locations with different localization in time and frequency [6, 7]. This, in turn, may obscure visualization of transient widearea phenomena or limit the ability of the method to deal with data measured at different sampling rates or containing missing information.
Analysis and characterization of timesynchronized system measurements requires mathematical tools that are adaptable to the varying system conditions, accurate and fast, while reducing the complexity of the data to make them comprehensible and useful for realtime decisions. This is particularly true in the study of large datasets of dynamic processes whose energy changes with time or frequency.
In this paper, a statistically based, datadriven framework that integrates the use of a waveletbased empirical orthogonal function (EOF) analysis and the method of snapshots is proposed to identify and extract dynamically independent spatiotemporal patterns from timesynchronized data. Extensions to current approaches to estimating propagating and standing features in nearrealtime that can be associated with observed or measured data are discussed and numerical issues are addressed.
The procedure allows identification of the dominant spatial and temporal patterns in a complex dataset and is particularly well suited for the study the temporal evolution of critical modal parameters. 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 motion, and modal content, 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 timesynchronized phasor measurements of a real event in Mexico is assessed. Results show that the proposed method can provide accurate estimation of nonstationary effects, modal frequency, timevarying modes shapes, and time instants of intermittent or irregular transient behavior associated with abrupt changes in system topology or operating conditions in a nearrealtime setting.
2. Statistical Characterization of Measured Data
2.1. Empirical Orthogonal Function Analysis
Timesynchronized 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–4]. When data is available at multiple locations, a spatiotemporal model is required that represents system behavior. Figure 1 shows a conceptual representation of PMU data collected at different spatial locations.
Assume, in order to introduce the more general ideas that follow, that measured data is available at n spatial locations (measurement locations) defined by , at instants in time, . Let be a spacetime scalar field representing a time trace, where is a set of spatial variables (i.e., measurement locations) on a space , and is the time at which the observations are made.
Using this notation, the set of data can be represented by an N×ndimension ensemble (observation) matrix, , of the form
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.
Empirical orthogonal function (EOF) analysis provides a basis for the modal decomposition of an ensemble of data in terms of the smallest possible number of basic modes or proper orthogonal modes (POMs, most energetic global functions) [5–8]. This decomposition is of the form [5] where are the temporal amplitudes or eigenvectors, are the spatial component maps or eigenfunctions, and are the statistical modes or POMs, where . Each of these maps represents a standing oscillation, and the temporal coefficients, , represent how this pattern oscillates through time [9–14].
The temporal and spatial components are calculated from the eigenvectors and eigenfunctions of the covariance matrix
In standard realEOF analysis, matrix is real and symmetric and possesses a set of orthogonal eigenvectors with positive (real) eigenvalues. The vectors are called empirical orthogonal functions (EOFs); the map associated with each eigenvector represents a pattern which is statistically independent of the others and spatially orthogonal to them.
Two important properties of EOFs are that the spatial distributions are orthogonal and their time series are uncorrelated. The method accounts for spatial and temporal changes and can be used to extract dynamic patterns from measured data. Measured data, however, may exhibit quite different dynamics at each system location or exhibit abrupt changes in modal quantities that cannot be captured using existing stationary models [11].
2.2. WaveletBased EOF Analysis
In realworld applications data are multiscale due to events occurring at different locations with different localization in time and frequency. To address the problem of multiscale modeling in data, a waveletbased EOF approach has been developed. In this approach, the empirical orthogonal functions of the wavelet coefficients are computed at each scale and then combined the results at relevant scales. The approach is referred to as waveletbased EOF analysis [15, 16] and involves three main steps:(1)computing the wavelet decomposition for each column in the data matrix, (2)computing the covariance matrix of wavelet coefficients for each scale of interest, and combining the results at the dominant scales of interest, (3)extracting dynamic features from the selected scales.
Appendix A briefly summarizes wavelet analysis in the context of the proposed formulation. In what follows we review the theory behind EOF decomposition and present some results and extensions to conventional analysis to treat changes in space and temporal variability.
A complex, nearrealtime formulation with the ability to resolve localized information is then proposed.
3. NearReal Time EOF Analysis
As it was highlighted in the previous section, conventional EOF analysis assumes stationarity of the underlying dynamic process and is therefore not suitable for an online 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].
In order to overcome the above limitations, a sliding windowbased method is used to systematically analyze the observational data in nearrealtime. This allows to isolate and extract the portions of data where oscillations are present and improves numerical efficiency and accuracy.
3.1. TimeDependent Covariance Matrix
Assume that denotes a sequence of observations, collected at time instants , to , evenly spread throughout the time period, where is set of spatial variables (measurement locations).
A sliding windowbased approach has been combined with EOF analysis to resolve localized information. In this approach, a sliding window frame of fixed size, say , is shifted regularly throughout the data span from the beginning of the record to the end of the data as shown in Figure 2. The analysis window is then slid the same distance repeatedly and the covariance matrix is computed recursively. It should be stressed that can range from a single sample to the entire length of the record.
Referring to Figure 2, let matrix be the covariance matrix of size , where is the window size, is the number of sensors or PMUs, and k indicates the time instant at which the response matrix is computed.
More formally, given a set of available sample time series, the temporal autocorrelation matrix, , of extending over a window length is defined as or, equivalently,where, in the above, the covariance matrix can be real or complex.
In our formulation, matrix is a positive semidefinite matrix of size , and possesses a set of complex timedependent orthogonal eigenvectors, , that capture propagating features. Extensions of this approach to consider the case where the autocorrelation matrix is nonintegrable are discussed in [17, 18]. This is, however, not considered here.
An interesting and useful interpretation of the covariance matrix can be obtained by rewriting (3.1a) as
Now if we let , we can write
This analysis suggests that the covariance matrix at time step can be computed recursively using information from both, the previous time step and the current measurement using only elementary operations [19]. The stationarity test for the recursive covariance matrix can be made to demonstrate the stationarity of the dataset [20].
Similar to the previous development, it can be shown that for a fixed window size, , (3.1a) and (3.1b) can be rewritten as
Equation (3.4) provides an efficient method to approximate the timedependent covariance matrix: since the method is based on local information the technique is wellsuited for realtime 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 specific segments of the observed record. As discussed in our numerical simulations, shortwidth windows provide improved localized information but result in enhanced computational effort. We have found that, the best estimate for the autocorrelation matrix C can be obtained using the averaged twopoint correlation function [5].
3.2. Complex Formulation
In an effort to extend standard realEOF 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. Refer to Esquivel and Messina for more details [11, 21, 22].
Assume then that is a complex matrix with , where the subscripts , indicate the real and imaginary vectors. Under these assumptions, the autocorrelation matrix in (3.1a) and (3.1b) [11, 23] becomes
Given a model of the form (3.5), it is straightforward to show that is a symmetrical matrix, and that is an asymmetric matrix or hemisymmetric matrix, that is, . Since the symmetrical matrix is a particular case of the Hermitian matrix, all of its eigenvectors are real; the elements of the asymmetrical matrix are all purely imaginary and its eigenvectors are complex conjugate. Note that the bases for the complex autocorrelation matrix , are now defined as for the real part and for the imaginary part.
Appendix B contains further discussion of this subject.
3.3. Propagating Features
Many power system phenomena derive from nonlinear interactions between traveling 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 are calculated, the original field can be approximated by a spatiotemporal model [4, 8, 10]. Consider a field composed of standing and traveling components, denoted by of the form
From EOF analysis [11], the complex field can be expressed as the complex expansion where are the temporal and spatial amplitude functions, respectively, and , are the frequency functions and wave component to be determined.
Let now the complex field (3.6) be expanded in terms of a truncated EOF basis of and modes as [12] where the complex, timedependent coefficients, and are given by
Equation (3.8) can now be recast in the form where and are the temporal and spatial amplitude functions associated to the wave component of the decomposition, respectively, and and are the temporal and spatial phase functions, respectively, [6, 7]. The original time series can then be reconstructed from the real part of (3.10). The details are omitted.
In the developed models, a criterion for choosing the number of relevant modes is given by the energy percentage contained in the and modes: where denotes the Frobenius norm. From (3.11), and are given in difference as
Equations (3.7)–(3.12) provide a complete characterization of any propagating effects and periodicity in the original data field which might be obscured by normal crossspectral analysis.
It might be remarked that, in the special case of real analysis, these expressions are simplified to the standard definitions.
4. Recursive Algorithm
A flowchart of the proposed approach is shown in Figure 3. In this plot, dashed boxes indicate the use of instantaneous information from the previous time step , and indicates the expectation value. This approach makes feasible the analysis and characterization of transient processes using realtime information. Variations to these approaches that extend their practical use to the realm of nearrealtime stability assessment and control are being investigated. These features will be discussed in future research.
The following sections describe the application of multivariate statistical techniques to measured data including the estimation of instantaneous parameters.
5. Simulations
5.1. Records of Selected Signals
The data used for this study were recorded by multiple phasor measurement units during a real event in northern Mexico. See [21, 22] for more details about this event.
The main event that originated the oscillations was a failed temporary interconnection of the Northwestern regional system to the Mexican interconnected system through a 230 kV line between Mazatlan Dos and Tres Estrellas substations. As a result of topological changes and load shedding, the observed oscillations exhibit highly complex phenomena including transient motions characterized by changing frequency content and variations in the mode shapes of critical electromechanical modes.
Among the existing PMU locations, frequency measurements at three major substations round the north, northwestern and northeastern systems were selected for study: Hermosillo (H), Mazatlan Dos (MZD) and Tres Estrellas (TTE). Figure 4 is an extract from PMU measurements of this event showing the observed oscillations of selected bus frequencies.
Measurements were recorded over 400 ms collected at a rate of 0.20 samples per second.
5.2. Analysis of Propagating Features
System measurements in this plot demonstrate significant variability suggesting a nonstationary process in both space and time. As observed in these plots, the most prominent variations occur in the interval during which the oscillation 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 : 21.
Based on the frequency data collected by PMUs, EOF analysis was applied to reveal spatial and temporal dynamics. As a first step towards the development of an empirical basis, the covariance matrix was formed by the ensembles of the frequency observations at different system locations, that is,
To visualize the complex temporal and spatial dynamics that takes placed in the system following the failed interconnection, each time series is augmented with an imaginary component using Hilbert analysis and the complex EOF method is employed to approximate the original data. Further, in order to improve the ability of the method to capture temporal behavior, the individual time series are separated into their timevarying mean and fluctuating components; complexEOF analysis is applied to the fluctuating field to decompose the spatiotemporal data into orthogonal modes associated to the standing and traveling wave components.
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 compares the signal reconstructed using conventional EOF and the proposed online formulation, whilst Table 1 gives the rms error. A twosample data analysis window was used in the online analysis of system characteristics. For simplicity, only the traveling wave mode is used in the computation. Simulation results clearly show that nearrealtime approximations enhance the ability of the method to capture local information.

(a)
(b)
(c)
Also of interest, Figure 6 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 nearrealtime implementation enables to single out the times at which sudden changes in system behavior take place. Further, Figure 7 shows the mode shape of dominant electromechanical modes obtained using both approaches.
(a)
(b)
(a)
(b)
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 [9, 23], and control actions. 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, 22]. 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.
The proposed approach provides an automated way to estimate mode shapes without any prior information of the time intervals of interest. The computational effort for analysis of large blocks of data is directly related to the number of spatial variables and the size of the observation window considered. In our numerical simulations the data analysis window corresponds to ongoing that research is being conducted to estimate the appropriate window size for realtime applications.
Finally, Figure 8 shows the time evolution of the dominant travelling wave at different system locations obtained using the wavelet empirical orthogonal function analysis described above. 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.
6. Conclusion
Widearea, realtime 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 paper, an alternative technique based on timedependent complex EOF analysis of measured data is proposed to resolve the localized nature of transient processes and to extract dominant temporal and spatial information. A method of spatially decomposing oscillation patterns in nearrealtime into their standing and travelling parts is presented.
By combining a sliding window approach with complex EOF, a novel framework for online 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, timevarying mode shapes, and time instants of intermittent transient responses. This information may be important in determining strategies for widearea 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 realtime system monitoring, protection and control will be addressed in future work.
Appendices
A. Wavelet Transform
Consider a function belonging to a family of finite energy functions, that is
The continuous wavelet transform (WT) of is defined as [24] where is the basis wavelet function or mother wavelet, a is a scale parameter and b is a time parameter. The superscript denotes complex conjugation. It is assumed that is also a finite energy function satisfying the condition in which is the Fourier transform of defined as
The function can be reconstructed from by using the double integral representation where the scale parameter a is restricted to positive values only.
B. Singular Value Decomposition
This section reviews the singular value decomposition and its features that relevant in the context of the proposed formulation.
Using the SVD decomposition, a complex matrix can be represented as where and are the magnitude and phase of , respectively, are left singular vectors and are singular values.
Using (B.1), we can rewrite (3.5) as or As it is seen from (B.1), the columns of and are the eigenvectors of real and imaginary parts of , and the columns of and are the eigenvectors of real and imaginary parts of , respectively, [10, 11]. The n singular values on the diagonal of and are the square roots of the nonzero eigenvalues of the real and imaginary parts of both and divided by the number of samples .
From the decomposition given in (B.2), it 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, for this instant both are symmetrical matrixes.
The imaginary part measures the grade of asymmetries when the sum of both matrixes is different from zero; this is used to define the existence of traveling wave components.
References
 A. R. Messina, Interarea Oscillations In Power Systems: A Nonlinear and Nonstationary Perspective, Springer, Berlin, Germany, 2009.
 J. F. Hauer, W. A. Mittelstadt, R. Adapa, M. K. Donelly, and W. H. Litzenberger, “Power system dynamics and stability—section 8: direct analysis of widearea dynamics,” in CRC Electric Power Engineering Handbook, chapter 11, CRC Press, Cleveland, Ohio, USA, 2000. View at: Google Scholar
 J. F. Hauer and J. G. DeSteese, “A tutorial on detection and characterization of special behavior in large electric power systems,” Tech. Rep. PNL14655, Pacific Northwest National Laboratory, Cambridge, Mass, USA, July 2004. View at: Google Scholar
 J.P. Chilès and P. Delfiner, Geostatistics: Modeling Spatial Uncertainty, Wiley Series in Probability and Statistics: Applied Probability and Statistics, John Wiley & Sons, New York, NY, USA, 1999. View at: MathSciNet
 P. Holmes, J. L. Lumley, and G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge Monographs on Mechanics, Cambridge University Press, Cambridge, Mass, USA, 1996. View at: MathSciNet
 T. P. Barnett, “Interaction of the monsoon and Pacific trade wind system at interannual time scales. Part I: the equatorial zone,” Monthly Weather Review, vol. 111, no. 4, pp. 756–773, 1983. View at: Google Scholar
 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,” American Meteorological Society, vol. 15, pp. 764–774, 1997. View at: Google Scholar
 N. Aubry, R. Guyonnet, and R. Lima, “Spatiotemporal analysis of complex signals: theory and applications,” Journal of Statistical Physics, vol. 64, no. 34, pp. 683–739, 1991. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 C. Wolter, M. A. Trindade, and R. Sampaio, “Obtaining mode shape through the KarhunenLoeve expansion for distributedparameter linear systems,” Vibrations and Acoustic, vol. 10, pp. 444–452, 2001. View at: Google Scholar
 H. Dankowicz, P. Holmes, G. Berkooz, and J. Elezgaray, “Local models of spatiotemporally complex fields,” Physica D, vol. 90, no. 4, pp. 387–407, 1996. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 P. Esquivel and A. R. Messina, “Complex empirical orthogonal function analysis of widearea system dynamics,” in Proceedings of IEEE Power and Energy Society General Meeting (PES '08), July 2008. View at: Publisher Site  Google Scholar
 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, no. 1–3, pp. 147–169, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 A. Hannachi, I. T. Jolliffe, and D. B. Stephenson, “Empirical orthogonal functions and related techniques in atmospheric science: a review,” International Journal of Climatology, vol. 27, no. 9, pp. 1119–1152, 2007. View at: Publisher Site  Google Scholar
 J. Terradas, R. Oliver, and J. L. Ballester, “Application of statistical techniques to the analysis of solar coronal oscillations,” Astrophysical Journal, vol. 614, no. 1, pp. 435–447, 2004. View at: Publisher Site  Google Scholar
 D. Mwale, T. Y. Gan, K. Devito, C. Mendoza, and U. Silins, “Precipitation variability and its relationship to hydrologic variability in Alberta,” Hydrological Processes, vol. 23, pp. 3040–3056, 2009. View at: Google Scholar
 S. Uhlig, “On the complexity of Internet traffic dynamics on its topology,” Telecommunication Systems, vol. 43, no. 34, pp. 167–180, 2010. View at: Google Scholar
 M. Li and S. C. Lim, “Modeling network traffic using generalized Cauchy process,” Physica A, vol. 387, no. 11, pp. 2584–2594, 2008. View at: Publisher Site  Google Scholar
 M. Li, “Generation of teletraffic of generalized Cauchy type,” Physica Scripta, vol. 81, no. 2, pp. 1–10, 2010. View at: Google Scholar
 J. G. Proakis and D. G. Manolakis, Digital Treatment of Signals, Pearson Education, Delhi, India, 2007.
 M. Li, W.S. Chen, and L. Han, “Correlation matching method of the weak stationarity test of LRD traffic,” Telecommunication Systems, vol. 43, no. 34, pp. 181–195, 2010. View at: Google Scholar
 A. R. Messina and V. Vittal, “Extraction of dynamic patterns from widearea measurements using empirical orthogonal functions,” IEEE Transactions on Power Systems, vol. 22, no. 2, pp. 682–692, 2007. View at: Publisher Site  Google Scholar
 A. R. Messina, V. Vittal, D. RuizVega, and G. EnríquezHarper, “Interpretation and visualization of widearea PMU measurements using Hilbert analysis,” IEEE Transactions on Power Systems, vol. 21, no. 4, pp. 1763–1771, 2006. View at: Publisher Site  Google Scholar
 B. F. Feeny, “A complex orthogonal decomposition for wave motion analysis,” Journal of Sound and Vibration, vol. 310, no. 12, pp. 77–90, 2008. View at: Publisher Site  Google Scholar
 C. Torrence and G. P. Compo, “A practical guide to wavelet analysis,” Bulletin of the American Meteorological Society, vol. 79, no. 1, pp. 61–78, 1998. View at: Google Scholar
Copyright
Copyright © 2010 A. R. Messina et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.