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.

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 [14].

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 [57]. 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 [14]. 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.

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 space-time 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×n-dimension 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) [58]. 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 [914].

The temporal and spatial components are calculated from the eigenvectors and eigenfunctions of the covariance matrix

In standard real-EOF 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. 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.

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, 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 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 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.

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 time-dependent 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 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 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 spatio-temporal 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, time-dependent 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 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 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 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.

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 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.

Figure 5 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.

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 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.

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, 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

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.

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.