# Time-Dependent Statistical Analysis of Wide-Area Time-Synchronized Data

Authors: A. R. Messina; P. Esquivel; F. Lezama
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

Characterization of spatial and temporal changes in the dynamic patterns of a nonstationary process is a problem of great theoretical and practical importance.  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 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, 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 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 time-frequency dynamics of oscillatory processes.
The efficiency and accuracy of the developed procedures for extracting localized information of power system behavior from time-synchronized phasor measurements of a real event in Mexico is assessed.


---

## Body

## 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 low-scale 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 wide-area 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, data-driven 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 wide-area 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 time-synchronized 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 real-time 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, data-driven framework that integrates the use of a wavelet-based 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. Extensions to current approaches to estimating propagating and standing features in near-real-time 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 time-synchronized 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, time-varying modes shapes, and time instants of intermittent or irregular transient behavior associated with abrupt changes in system topology or operating conditions in a near-real-time setting.


## 2. Statistical Characterization of Measured Data

### 2.1. Empirical Orthogonal Function 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–4]. When data is available at multiple locations, a spatio-temporal model is required that represents system behavior. Figure 1 shows a conceptual representation of PMU data collected at different spatial locations.Figure 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 atn spatial locations (measurement locations)  defined by xj,j=1,…,n, at  N instants in time, tk,k=1,…,N. Let  u(xj,tk) be a space-time scalar field representing a time trace, where xj,  is a set of spatial variables (i.e., measurement locations) on a space  Ωk, and tk,  is the time at which the observations are made.Using this notation, the set of data can be represented by anN×n-dimension ensemble (observation) matrix, F(x,t),  of the form(2.1)F(x,t)=[f1(x1,t)f2(x2,t)⋯fn(xn,t)]=[u̅(x1,t1)u̅(x2,t1)⋯u̅(xn,t1)u̅(x1,t2)u̅(x2,t2)⋯u̅(xn,t2)⋮⋮⋱⋮u̅(x1,tN)u̅(x2,tN)⋯u̅(xn,tN)].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](2.2)F̂(x,t)=∑i=1maiφi*(x),
where ai are the temporal amplitudes or eigenvectors, φi(x) are the spatial component maps or eigenfunctions, and m are the statistical modes or POMs, where  ai=F(x,t)φi. Each of these maps represents a standing oscillation, and the temporal coefficients, ai, 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(2.3)C=1NFT(x,t)F(x,t).In standard real-EOF analysis, matrixC is real and symmetric and possesses a set of orthogonal eigenvectors with positive (real) eigenvalues. The vectors φi(x) 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. Wavelet-Based EOF Analysis

In real-world 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 wavelet-based 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 wavelet-based 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.AppendixA 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, near-real-time formulation with the ability to resolve localized information is then proposed.


## 2.1. Empirical Orthogonal Function 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–4]. When data is available at multiple locations, a spatio-temporal model is required that represents system behavior. Figure 1 shows a conceptual representation of PMU data collected at different spatial locations.Figure 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 atn spatial locations (measurement locations)  defined by xj,j=1,…,n, at  N instants in time, tk,k=1,…,N. Let  u(xj,tk) be a space-time scalar field representing a time trace, where xj,  is a set of spatial variables (i.e., measurement locations) on a space  Ωk, and tk,  is the time at which the observations are made.Using this notation, the set of data can be represented by anN×n-dimension ensemble (observation) matrix, F(x,t),  of the form(2.1)F(x,t)=[f1(x1,t)f2(x2,t)⋯fn(xn,t)]=[u̅(x1,t1)u̅(x2,t1)⋯u̅(xn,t1)u̅(x1,t2)u̅(x2,t2)⋯u̅(xn,t2)⋮⋮⋱⋮u̅(x1,tN)u̅(x2,tN)⋯u̅(xn,tN)].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](2.2)F̂(x,t)=∑i=1maiφi*(x),
where ai are the temporal amplitudes or eigenvectors, φi(x) are the spatial component maps or eigenfunctions, and m are the statistical modes or POMs, where  ai=F(x,t)φi. Each of these maps represents a standing oscillation, and the temporal coefficients, ai, 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(2.3)C=1NFT(x,t)F(x,t).In standard real-EOF analysis, matrixC is real and symmetric and possesses a set of orthogonal eigenvectors with positive (real) eigenvalues. The vectors φi(x) 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. Wavelet-Based EOF Analysis

In real-world 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 wavelet-based 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 wavelet-based 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.AppendixA 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, near-real-time formulation with the ability to resolve localized information is then proposed.


## 3. Near-Real 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 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].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.


### 3.1. Time-Dependent Covariance Matrix

Assume thatfk(x)=[fk(x1),fk(x2),…,fk(xn)] denotes a sequence of observations, collected at time instants tk, to k=1,2,…,∞, evenly spread throughout the time period, where xn is set of spatial variables (measurement locations).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τ, 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.Figure 2
Visualization of PMU data in terms of spatial and temporal information.Referring to Figure2, let matrix Ck(x) be the covariance matrix of size  τ×n, where τ 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.More formally, given a set of available sample time series, the temporal autocorrelation matrix,Ck, of Fk(x) extending over a window length is defined as(3.1a)CkΔ̲̲FkH(x)Fk(x)k=1k∑i=1kfiH(x)fi(x),
or, equivalently,(3.1b)CkΔ̲̲Fk(x)FkH(x)k=1k∑i=1kfi(x)fi(x)H,where, in the above, the covariance matrix can be real or complex.In our formulation, matrixCk, is a positive semidefinite matrix of size n×n, and possesses a set of complex time-dependent orthogonal eigenvectors, φi(t),i=1,…,n, 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(3.2)kCk=∑i=1k-1fi(x)fiH(x)+fk(x)fkH(x).Now if we letCk-1=(1/(k-1))∑i=1k-1fi(x)fi(x)H, we can write(3.3)Ck=k-1kCk-1+1kfk(x)fkH(x).This analysis suggests that the covariance matrix at time stepk can be computed recursively using information from both, the previous time step (k-1) 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(3.4)Ck=τ-1τCk-1+1τfk(x)fkH(x).Equation (3.4) provides an efficient method to approximate the time-dependent co-variance 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 specific 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 estimate for the autocorrelation matrix C can be obtained using the averaged two-point correlation function [5].


### 3.2. Complex Formulation

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. Refer to Esquivel and Messina for more details [11, 21, 22].Assume then thatFk=[FR+jFI]k is a complex matrix with j=-1, where the subscripts R,I,  indicate the real and imaginary vectors. Under these assumptions, the autocorrelation matrix Ck in (3.1a) and (3.1b) [11, 23] becomes(3.5)FkFkH=[FRFRT+FIFIT]k+j[FIFRT-FRFIT]k,FkHFk=[FRTFR+FITFI]k+j[FRTFI-FITFR]k.Given a model of the form (3.5), it is straightforward to show that CR=CRT is a symmetrical matrix, and that CI is an asymmetric matrix or hemisymmetric matrix, that is, CIT=-CI. 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 Ck, are now defined as φR(x) for the real part and φI(x) for the imaginary part.AppendixB 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 ofC are calculated, the original field can be approximated by a spatio-temporal model [4, 8, 10]. Consider a field fk composed of standing and traveling components, denoted by fswc,ftwc of the form
(3.6)fk=[fswc+ftwc]k.From EOF analysis [11], the complex field can be expressed as the complex expansion(3.7)[fswc]k=[∑i=1pRR(i)(t)SR(i)(x)cos(ωR(i)t)]k,[ftwc]k=[∑i=1qRI(i)(t)SI(i)(x)cos(ωI(i)t+KI(i)x+π)]k,
where R,S are the temporal and spatial amplitude functions, respectively, and ω, K 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 p and q modes as [12](3.8)f̂k=[∑i=1pAR(i)(t)φR(i)H(x)+j∑i=1qAI(i)(t)φI(i)H(x)]k,
where the complex, time-dependent coefficients, AR(i)(t) and AI(i)(t) are given by(3.9)AR(i)=[fφR(i)(x)]k,AI(i)=[fφI(i)(x)]k.Equation (3.8) can now be recast in the form(3.10)f̂k=[∑i=1pRR(i)(t)SR(i)(x)ej(θR(i)+φR(i))]k+[∑i=1qRI(i)(t)SI(i)(x)ej(θI(i)+φI(i)+π)]k,
where R(t) and S(x) are the temporal and spatial amplitude functions associated to the wave component of the decomposition, respectively, and θ(t) and φ(x) 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 thep and q modes:
(3.11)∑i=1p[Δλ(i)swc]k+∑i=1q[Δλ(i)twc]k∥fk∥F2=99%,
where ∥·∥F2 denotes the Frobenius norm. From (3.11), (Δλswc)k and (Δλtwc)k are given in difference as
(3.12)(Δλswc)k=(λswc)k+1-(λswc)k,(Δλtwc)k=(λtwc)k+1-(λtwc)k.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 cross-spectral analysis.It might be remarked that, in the special case of real analysis, these expressions are simplified to the standard definitions.


## 3.1. Time-Dependent Covariance Matrix

Assume thatfk(x)=[fk(x1),fk(x2),…,fk(xn)] denotes a sequence of observations, collected at time instants tk, to k=1,2,…,∞, evenly spread throughout the time period, where xn is set of spatial variables (measurement locations).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τ, 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.Figure 2
Visualization of PMU data in terms of spatial and temporal information.Referring to Figure2, let matrix Ck(x) be the covariance matrix of size  τ×n, where τ 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.More formally, given a set of available sample time series, the temporal autocorrelation matrix,Ck, of Fk(x) extending over a window length is defined as(3.1a)CkΔ̲̲FkH(x)Fk(x)k=1k∑i=1kfiH(x)fi(x),
or, equivalently,(3.1b)CkΔ̲̲Fk(x)FkH(x)k=1k∑i=1kfi(x)fi(x)H,where, in the above, the covariance matrix can be real or complex.In our formulation, matrixCk, is a positive semidefinite matrix of size n×n, and possesses a set of complex time-dependent orthogonal eigenvectors, φi(t),i=1,…,n, 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(3.2)kCk=∑i=1k-1fi(x)fiH(x)+fk(x)fkH(x).Now if we letCk-1=(1/(k-1))∑i=1k-1fi(x)fi(x)H, we can write(3.3)Ck=k-1kCk-1+1kfk(x)fkH(x).This analysis suggests that the covariance matrix at time stepk can be computed recursively using information from both, the previous time step (k-1) 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(3.4)Ck=τ-1τCk-1+1τfk(x)fkH(x).Equation (3.4) provides an efficient method to approximate the time-dependent co-variance 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 specific 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 estimate for the autocorrelation matrix C can be obtained using the averaged two-point correlation function [5].


## 3.2. Complex Formulation

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. Refer to Esquivel and Messina for more details [11, 21, 22].Assume then thatFk=[FR+jFI]k is a complex matrix with j=-1, where the subscripts R,I,  indicate the real and imaginary vectors. Under these assumptions, the autocorrelation matrix Ck in (3.1a) and (3.1b) [11, 23] becomes(3.5)FkFkH=[FRFRT+FIFIT]k+j[FIFRT-FRFIT]k,FkHFk=[FRTFR+FITFI]k+j[FRTFI-FITFR]k.Given a model of the form (3.5), it is straightforward to show that CR=CRT is a symmetrical matrix, and that CI is an asymmetric matrix or hemisymmetric matrix, that is, CIT=-CI. 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 Ck, are now defined as φR(x) for the real part and φI(x) for the imaginary part.AppendixB 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 ofC are calculated, the original field can be approximated by a spatio-temporal model [4, 8, 10]. Consider a field fk composed of standing and traveling components, denoted by fswc,ftwc of the form
(3.6)fk=[fswc+ftwc]k.From EOF analysis [11], the complex field can be expressed as the complex expansion(3.7)[fswc]k=[∑i=1pRR(i)(t)SR(i)(x)cos(ωR(i)t)]k,[ftwc]k=[∑i=1qRI(i)(t)SI(i)(x)cos(ωI(i)t+KI(i)x+π)]k,
where R,S are the temporal and spatial amplitude functions, respectively, and ω, K 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 p and q modes as [12](3.8)f̂k=[∑i=1pAR(i)(t)φR(i)H(x)+j∑i=1qAI(i)(t)φI(i)H(x)]k,
where the complex, time-dependent coefficients, AR(i)(t) and AI(i)(t) are given by(3.9)AR(i)=[fφR(i)(x)]k,AI(i)=[fφI(i)(x)]k.Equation (3.8) can now be recast in the form(3.10)f̂k=[∑i=1pRR(i)(t)SR(i)(x)ej(θR(i)+φR(i))]k+[∑i=1qRI(i)(t)SI(i)(x)ej(θI(i)+φI(i)+π)]k,
where R(t) and S(x) are the temporal and spatial amplitude functions associated to the wave component of the decomposition, respectively, and θ(t) and φ(x) 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 thep and q modes:
(3.11)∑i=1p[Δλ(i)swc]k+∑i=1q[Δλ(i)twc]k∥fk∥F2=99%,
where ∥·∥F2 denotes the Frobenius norm. From (3.11), (Δλswc)k and (Δλtwc)k are given in difference as
(3.12)(Δλswc)k=(λswc)k+1-(λswc)k,(Δλtwc)k=(λtwc)k+1-(λtwc)k.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 cross-spectral 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 Figure3. In this plot, dashed boxes indicate the use of  instantaneous information from the previous time step (k-1), and E(·) indicates the expectation value. 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 realm of near-real-time stability assessment and control are being investigated. These features will be discussed in future research.Figure 3
Algorithm for near-real-time computation of standing and propagating features from measured data.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). Figure4 is an extract from PMU measurements of this event showing the observed oscillations of selected bus frequencies.Figure 4
Time traces of recorded bus frequency swings recorded January 1, 2004 and detail of the oscillation buildup.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,(5.1)Fk=[fH(tk),fMZD(tk),fTTE(tk)].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 time-varying mean and fluctuating components; complex-EOF analysis is applied to the fluctuating field to decompose the spatio-temporal 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.Figure5 compares the signal reconstructed using conventional EOF and the proposed on-line formulation, whilst Table 1 gives the rms error. A two-sample data analysis window was used in the on-line analysis of system characteristics.  For simplicity, only the traveling wave mode is used in the computation. Simulation results clearly show that near-real-time approximations enhance the ability of the method to capture local information.Table 1
RMS error in the approximation.
SignalsOff-line (rms)On-line (rms)Hermosillo0.35910.2614Mazatlan Dos0.49280.4008Tres Estrellas0.75360.5477Reconstruction of the original data with the traveling wave mode identified using the off-line and on-line analysis.
(a)(b)(c)Also of interest, Figure6 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. Further, Figure 7 shows the mode shape of dominant electromechanical modes obtained using both approaches.(a) average energy of traveling wave mode using the off-line analysis, and (b) instantaneous energy of traveling wave mode using the on-line analysis.
(a)(b)(a) mode shape of traveling wave mode using off-line analysis, and (b) time-varying mode shape of traveling wave mode using the proposed on-line analysis method.
(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 on-going that research is being conducted to estimate the appropriate window size for real-time applications.Finally, Figure8 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.Figure 8
Spatial pattern of the leading mode showing temporal variability.


## 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). Figure4 is an extract from PMU measurements of this event showing the observed oscillations of selected bus frequencies.Figure 4
Time traces of recorded bus frequency swings recorded January 1, 2004 and detail of the oscillation buildup.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,(5.1)Fk=[fH(tk),fMZD(tk),fTTE(tk)].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 time-varying mean and fluctuating components; complex-EOF analysis is applied to the fluctuating field to decompose the spatio-temporal 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.Figure5 compares the signal reconstructed using conventional EOF and the proposed on-line formulation, whilst Table 1 gives the rms error. A two-sample data analysis window was used in the on-line analysis of system characteristics.  For simplicity, only the traveling wave mode is used in the computation. Simulation results clearly show that near-real-time approximations enhance the ability of the method to capture local information.Table 1
RMS error in the approximation.
SignalsOff-line (rms)On-line (rms)Hermosillo0.35910.2614Mazatlan Dos0.49280.4008Tres Estrellas0.75360.5477Reconstruction of the original data with the traveling wave mode identified using the off-line and on-line analysis.
(a)(b)(c)Also of interest, Figure6 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. Further, Figure 7 shows the mode shape of dominant electromechanical modes obtained using both approaches.(a) average energy of traveling wave mode using the off-line analysis, and (b) instantaneous energy of traveling wave mode using the on-line analysis.
(a)(b)(a) mode shape of traveling wave mode using off-line analysis, and (b) time-varying mode shape of traveling wave mode using the proposed on-line analysis method.
(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 on-going that research is being conducted to estimate the appropriate window size for real-time applications.Finally, Figure8 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.Figure 8
Spatial pattern of the leading mode showing temporal variability.


## 6. Conclusion

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 paper, an alternative technique based on time-dependent 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 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.



---

*Source: 751659-2010-06-20.xml*