Mode Shape Estimation using Complex Principal Component Analysis and k-Means Clustering

TL;DR

This paper introduces an empirical approach combining complex PCA, Hilbert Transform, and k-means clustering to identify low damped modes and mode shapes from frequency data, validated on simulated and real power system data.

Contribution

The paper presents a novel two-step method for mode shape estimation using complex PCA and clustering, applicable to real power system measurements.

Findings

Method accurately estimates dominant modes in simulated data.

Initial real data tests show promising results.

Approach effectively groups mode shape estimates into meaningful clusters.

Abstract

We propose an empirical method for identifying low damped modes and corresponding mode shapes using frequency measurements from a Wide Area Monitoring System. The method consists of two main steps: Firstly, Complex Principal Component Analysis is used in combination with the Hilbert Transform and Empirical Mode Decomposition to provide estimates of modes and mode shapes. The estimates are stored as multidimensional points. Secondly, the points are grouped using a clustering algorithm, and new averaged estimates of modes and mode shapes are computed as the centroids of the clusters. Applying the method on data resulting from a non-linear power system simulator yields estimates of dominant modes and corresponding mode shapes that are similar to those resulting from modal analysis of the linearized system model. Encouraged by the results, the method is further tested with real PMU data at transmission grid level. Initial results indicate that the performance of the proposed method is promising.

Equations

$$\Delta \dot{\mathbf{x}}=\mathbf{A}\Delta \mathbf{x}+\mathbf{B}\Delta \mathbf{u}$$

$$\begin{array}{rl}\Delta \mathbf{x}(t)& =\sum _{j=1}^n{\mathbf{\Phi }}_j{c}_j{e}^{{\lambda }_{j}t}\end{array}$$

$$\mathbf{X}=\left[\begin{array}{c}{\mathbf{x}}_1\\ {\mathbf{x}}_2\\ ⋮\\ {\mathbf{x}}_M\end{array}\right]=\left[\begin{array}{ccccc}{x}_1(t_1)& x_1(t_2)& \cdots & x_1(t_N)& \\ x_2(t_1)& x_2(t_2)& \cdots & x_2(t_N)& \\ ⋮& & \ddots & & \\ x_M(t_1)& x_M(t_2)& & x_M(t_N)& \end{array}\right]$$

$$\mathbf{S}=\left[\begin{array}{c}{\mathbf{s}}_1\\ {\mathbf{s}}_2\\ ⋮\\ {\mathbf{s}}_M\end{array}\right]={\mathbf{U}}^T\mathbf{X}$$

$$\begin{array}{r}\mathbf{C}=\frac{1}{1-N}\mathbf{X}{\mathbf{X}}^T\Rightarrow \left(1-N\right)\mathbf{C}=\mathbf{X}{\mathbf{X}}^T\end{array}$$

$$\begin{array}{rl}\left(1-N\right)\mathbf{\Lambda }=& \mathbf{S}{\mathbf{S}}^T\\ =& {\mathbf{U}}^T\mathbf{X}{\left({\mathbf{U}}^T\mathbf{X}\right)}^T\\ =& {\mathbf{U}}^T\mathbf{X}{\mathbf{X}}^T\mathbf{U}\\ =& {\mathbf{U}}^T\mathbf{CU}\end{array}$$

$$\mathbf{C}{\mathbf{u}}_j={\lambda }_j{\mathbf{u}}_j,j=1\dots M$$

$$\begin{array}{rl}& \mathbf{X}=({\mathbf{U}}^T{)}^{-1}\mathbf{S}=\mathbf{US}\\ & \Rightarrow {\mathbf{x}}_i=\sum _{j=1}^M{u}_{ij}{\mathbf{s}}_j\approx \sum _{j=1}^{M_{PC}}{u}_{ij}{\mathbf{s}}_j\end{array}$$

$${\mathbf{y}}_i={\mathbf{x}}_i+jH\left({\mathbf{x}}_i\right)$$

$$\mathbf{Y}=\left[\begin{array}{c}{\mathbf{y}}_1\\ {\mathbf{y}}_2\\ ⋮\\ {\mathbf{y}}_M\end{array}\right]$$

$$\left(1-N\right)\tilde{\mathbf{C}}=\mathbf{Y}{\mathbf{Y}}^H$$

$$\begin{array}{rl}& \mathbf{Y}=({\mathbf{V}}^H{)}^{-1}\mathbf{Z}=\mathbf{VZ}\\ & \Rightarrow {\mathbf{y}}_i=\sum _{j=1}^M{v}_{ij}{\mathbf{z}}_j\approx \sum _{j=1}^{M_{CPC}}{v}_{ij}{\mathbf{z}}_j\end{array}$$

$$\mathbf{S}={\mathbf{U}}^T\mathbf{X}$$

$${\mathbf{S}}^{\prime }=\mathbf{S}-\mathbf{R}$$

$$\mathbf{Y}={\mathbf{S}}^{\prime }+jH\left({\mathbf{S}}^{\prime }\right)$$

$$\mathbf{Z}={\mathbf{V}}^H\mathbf{Y}={\mathbf{V}}^H\left({\mathbf{S}}^{\prime }+jH({\mathbf{S}}^{\prime })\right)$$

$$\begin{array}{rl}\mathbf{X}& =\mathbf{US}\\ & =\mathbf{U}\left({\mathbf{S}}^{\prime }+\mathbf{R}\right)\\ & =\mathbf{U}\left(\mathrm{Re}\left(\mathbf{Y}\right)+\mathbf{R}\right)\\ & =\mathbf{U}\left(\mathrm{Re}\left(\mathbf{VZ}\right)+\mathbf{R}\right)\end{array}$$

$$\mathbf{X}=\mathrm{Re}\left(\mathbf{UVZ}\right)=\mathrm{Re}\left(\mathbf{WZ}\right)$$

$${\mathbf{x}}_i=\sum _{j=1}^M{w}_{ij}{\mathbf{z}}_j\approx \sum _{j=1}^{M_{CPC}}{w}_{ij}{\mathbf{z}}_j$$

$$\begin{array}{r}\mathbf{p}=[f,\sigma ,\mathrm{Re}(G_1),\mathrm{Im}(G_1),\mathrm{Re}(G_2),\mathrm{Im}(G_2)...\\ ...\mathrm{Re}(G_M),\mathrm{Im}(G_M)]\end{array}$$

$$\begin{array}{r}\mathbf{P}=\left[\begin{array}{c}{\mathbf{p}}_1\\ {\mathbf{p}}_2\\ ⋮\\ {\mathbf{p}}_Q\end{array}\right]=\left[\begin{array}{ccccc}f& \sigma & \mathrm{Re}(G_1)& \dots & \mathrm{Im}(G_M)\\ f& \sigma & \mathrm{Re}(G_1)& \dots & \mathrm{Im}(G_M)\\ ⋮& & & & \\ f& \sigma & \mathrm{Re}(G_1)& \dots & \mathrm{Im}(G_M)\end{array}\right]\end{array}$$

Full text

©2019 IEEE. DOI: 10.1109/SGSMA.2019.8784556

I Introduction

Power system operators are facing the increasingly evident challenge of operating the grid securely during complex, uncertain and rapidly changing generation and demand. This challenge can and must be met by providing the operators with better tools for monitoring and control of the grid. In this paper we focus on the monitoring part and how the operators can improve their situational awareness and gain better knowledge about stability properties and dynamic phenomena occurring in the grid.

Worldwide, Phasor Meaurement Units (PMUs) are being installed within transmission grids at increasingly faster rates. It is expected that the information from PMUs will replace the current remote terminal units (RTUs) as source of voltage and current measurements to SCADA/EMS systems. This will change the way operators can monitor the system state. First of all, Wide Area Monitoring Systems opens new possibilities for getting fast and precise information about system dynamics. While WAMS have been in use many places for offline studies and fault analysis, the development of applications for online use has not yet matured. One reason for this may be that it takes time for operators to get used to and take advantage of the new technology. Research is therefore needed in order to develop applications that are readily accepted and understood by the operators. This paper aims to contribute to that development.

The method proposed in this paper takes advantage of the high number of frequency (or potentially voltage angle) measurements available in the transmission grid in order to extract precise information about low damped oscillatory modes and their observability.

Multiple empirical methods for estimating information about modes have been proposed in the litterature. In[1], Prony analysis is used to extract frequencies and damping ratios from measured oscillations. This method, however, does not provide information about mode shapes. In [2], a method for estimating modes and mode shapes using various mode decomposition techniques and Power Spectral Density analysis is proposed. In [3], a Bayesian approach is proposed for monitoring of electromechanical modes.

The proposed method presented in this paper is based on using Principal Component Analysis (PCA) and k-Means clustering for estimating modes and corresponding mode shapes from PMU frequency measurements. The method consists of two main parts: In Part I, the oscillatory behavior measured in a number of measurements is decomposed into a few main oscillatory components using a variant of Complex Principal Component Analysis (CPCA). The decomposition allows the oscillatory behavior to be reduced to a set of parameters, which can be stored in the form of a multidimensional point, referred to as an observation. Finally, in Part II of the method, the observations from Part I are clustered using the k-Means algorithm, to determine which of the observations to include when computing averaged modes and mode shapes.

CPCA is well described in [4], where it is aimed at extracting dynamic patterns in geophysical data sets. Variants of CPCA has been applied to power system analysis in[5, 6, 7]. The implementation used in this paper is slightly different, but functions in much the same way. An important contribution from this paper is the addition of the additional layer, referred to as Part II of the method, where the information provided by CPCA is processed further into more easily understandable information, i.e. by filtering and clustering observations and computing averaged quantities.

We present and describe the necessary background theory in section II, and describe the proposed method in section III. Application examples on simulated as well as real PMU data are presented in section IV. Finally, discussion and conclusions are given in sections V and VI.

II Background theory

This section provides a short introduction to modal analysis, and presents theory on the most important tools used by the proposed method.

II-A Modal Analysis

Modal analysis of power systems provides useful knowledge of the system dynamics, and can be used for multiple purposes, including assessment of small signal stability, and design and tuning of power oscillation dampers [8].

Fundamental to modal analysis of power systems is the linearization of the dynamic equations of the system. The state space form, assumed familiar to the reader, is given as follows[9]:

[TABLE]

Through eigendecomposition of the system matrix $\mathbf{A}$, the dynamic modes of the system are determined from the eigenvalues, and the observability of the modes (mode shapes) are contained within the elements of the right eigenvectors. It can be shown that the response of the system can be written

[TABLE]

where

[TABLE]

The element $\Phi_{ij}$ of $\mathbf{\Phi}_{j}$ describes the amplitude and phase with which mode $j$ appears in state $\Delta x_{i}$.

When performing modal analysis analytically, difficulties arise when modelling large systems due to insufficient knowledge of the system topology, component characteristics and controller tuning. Therefore, methods are developed for online identification of electromechanical modes that are based on measurements only. However, whether the methods are model based utilizing a state estimator or purely based on measurements, the concept of modes is inherently related to linear theory. Power systems, however, are nonlinear by nature and what appears as critical modes may change rapidly during a disturbance or any changes in operation.

II-B Principal Component Analysis

Principal Component Analysis (PCA, sometimes referred to as Empirical Orthogonal Functions) is a powerful dimensionality reduction technique, which allows the variance in a set of time series to be decomposed into a few orthogonal Principal Components (PCs) that explain the main share of the variance. A brief introduction of the method is given below, based on the thorough description of the method given in [10].

Assuming we have $M$ measurement series, each containing $N$ samples, assembled in a $M\times N$ matrix as follows:

[TABLE]

Here, the notation $x_{i}(t_{k})$ denotes the value of series $i$ at time $t_{k}$. We assume that each of the series has zero mean (if this is not the case, the mean is subtracted prior to the analysis).

We want to transform the correlated series $\mathbf{x}_{1},\mathbf{x}_{2}\dots\mathbf{x}_{M}$ into the uncorrelated series $\mathbf{s}_{1},\mathbf{s}_{2}\dots\mathbf{s}_{M}$ using the the linear transformation

[TABLE]

To find the matrix $\mathbf{U}$, we start by establishing the correlation matrix of $\mathbf{X}$:

[TABLE]

Similarly, the covariance matrix of $\mathbf{S}$ is given by

[TABLE]

Since we want the series $\mathbf{S}$ to be uncorrelated, we want the matrix $\mathbf{\Lambda}$ to have non-zero elements only along the diagonal. Thus, we want to choose the matrix $\mathbf{U}$ such that $\mathbf{C}$ is diagonalized. The diagonalization is carried out by performing an eigendecomposition of $\mathbf{C}$, where the eigenvalues and eigenvectors are given by

[TABLE]

The relative magnitudes of the eigenvalues describe how much variance the corresponding PC explains. The eigenvalues are sorted based on magnitude in descending order and collected in the eigenvalue matrix $\mathbf{\Lambda}$, and the eigenvectors are collected as column vectors in the transformation matrix $\mathbf{U}$ in the same order. Finally, this gives us the $M$ uncorrelated series $\mathbf{S}$, which can be referred to as the Scores of the PCs. If the input series $\mathbf{X}$ are highly correlated, then we can represent a high share of the variance using only a low number $M_{PC}<<M$ of the first PCs.

We can also invert the transformation in (4). Since it can be shown that $\mathbf{U}$ is an orthonormal matrix [10], the inverse of the matrix $\mathbf{U}$ is simply the transpose of the matrix:

[TABLE]

From the above relation we can state that PC $j$ appears in measurement $i$ with a magnitude given by the coefficient $u_{ij}$.

II-C Complex Principal Component Analysis

One limitation with conventional PCA is that a travelling wave, occuring in multiple measurements with different phase, can not be captured by only one PC. However, an extended version of PCA, often referred to as Complex Principal Component Analysis (CPCA), is better suited for this purpose. The method, described in detail in [4], uses the same equations as conventional PCA presented above, but some additional steps are required prior to the analysis.

CPCA is, in short, the same equations as presented above applied to complex time series, where the real part of the series is the input series, and the complex part is generated by applying the Hilbert Transform. To get sensible results when using the Hilbert Transform, it is important to make sure that the signal to be transformed contains only one frequency, and that any non-oscillatory trends are subtracted. This can be achieved using Empirical Mode Decomposition (EMD), or by applying PCA in two layers, discussed further in section II-D.

The Hilbert Transform is also strongly influenced by end effects [4]. When taking the transform of a measured series, this often results in very high values at the beginning and end of the series. To avoid that this influences consecutive steps in the algorithm, roughly $10\text{\,}\%$ of the series should be tapered in each end[4] after taking the Hilbert Transform.

After applying the necessary pre-processing to ensure a well-posed Hilbert Transform, the complex series can be written on the form

[TABLE]

where $H(\cdot)$ denotes the Hilbert Transform, and $\mathbf{x}_{i}$ is the real input series (with zero mean).

The complex series are collected in a matrix,

[TABLE]

The complex covariance matrix is given by

[TABLE]

where the superscript $(\cdot)^{H}$ denotes the conjugate transpose (not to be confused with the Hilbert Transform). Diagonalizing results in real eigenvalues[4] and complex eigenvectors. Analogously to (8), we have

[TABLE]

where $\mathbf{V}$ is the matrix of complex eigenvectors $\mathbf{v}_{i}$, with elements $v_{ij}$, and $\mathbf{Z}$ contains the Scores of the CPCs $\mathbf{z}_{j}$ as column vectors. The complex components of the eigenvectors describe the amplitude and phase of the contribution of each of the CPCs to each of the measurements. Here, the correlation matrix is Hermitian [4], the matrix $\mathbf{V}$ is unitary, and thus its inverse is the conjugate transpose.

II-D Two-layer combination of PCA and CPCA

Applying the Hilbert Transform requires the input signal to have only one frequency, and that any non-oscillatory trends are subtracted. This can be achieved by using EMD to decompose the input signal into Intrinsic Mode Functions (IMFs) and a Residual[7]. The Hilbert Transform can then be applied to the IMFs, before CPCA is applied on the resulting complex series.

An alternative to using EMD for decomposition into IMFs, is to apply a to a two-layer structure of PCA: In the first layer conventional PCA is applied on the input series, resulting in a number of PCs. In the second layer, the PCs from the first layer are detrended by subtracting the residual from EMD, the Hilbert Transform is applied, before finally CPCA is applied. Thus, EMD is used only for computing the non-oscillatory residual, not for decomposition into IMFs. Detrending by EMD is described in [11].

In equations we can write, for the first layer (where PCA is applied),

[TABLE]

The scores from the first layer are detrended (i.e. the residual computed using EMD is subtracted),

[TABLE]

the Hilbert Transform is applied,

[TABLE]

(Here, $H(\mathbf{S}^{\prime})$ denotes series-wise application of the Hilbert Transform, similarly as in (9). Also, although not shown with equations here, in a practical implementation roughly $10\text{\,}\%$ should be tapered in each end at this stage, to limit end effects.)

This gives the complex series which constitutes the input of the second layer (where CPCA is applied),

[TABLE]

Combining the above equations, we can write

[TABLE]

Neglecting the residual and defining the matrix $\mathbf{W}=\mathbf{UV}$, we get

[TABLE]

Finally, we can write

[TABLE]

This gives the contribution of the CPC $\mathbf{z}_{j}$, to the measurement $\mathbf{x}_{i}$ as the coefficient $w_{ij}$, i.e. element $(i,j)$ of the matrix $\mathbf{W}$.

Note that if no components are discarded in either of the layers, then the matrices $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{W}$ will be $M\times M$-matrices. If components are discarded, such that only $M_{PC}$ components are kept from the first layer and $M_{CPC}$ from the second layer, $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{W}$ will be of dimension $M\times M_{PC}$, $M_{PC}\times M_{CPC}$ and $M\times M_{CPC}$, respectively.

Comparing modal analysis with CPCA as described, we see that the coefficients $w_{ij}$ in (19) become very similar to the observability phasors $\phi_{ij}$ in (2). Thus, if the decomposition successfully produces CPCs with frequency and damping similar to that of the modes of the system, then the coefficients $w_{ij}$ can be considered estimates of observability mode shapes.

II-E k-Means Clustering

Clustering, in general, refers to the task of dividing a set of observations into groups, such that observations belonging to the same group are relatively similar, while the observations belonging to different groups are relatively different [12]. In the proposed method, clustering is used to determine which mode estimates, or observations, are similar enough to be included when computing averages of different modes and mode shapes. The k-Means clustering algorithm is widely used due to its computational simplicity. However, the method is prone to noise [13], and requires the number of clusters to be determined before the clustering.

To determine the number of clusters, a simple approach is to carry out the clustering multiple times with a range of different numbers of clusters, and then determining the most optimal solution by assessing some measure of the ”goodness” of each clustering. The Silhouette index can be adopted for this purpose. For a thorough description of the k-Means algorithm and formal definitions of the silhouette index, the interested reader is referred to [13].

II-F Exponential fitting by Linear Regression

Linear regression can be used for estimating the decay rate of an exponentially decaying series. We would like to describe the data $Y_{i}$ measured at $X_{i}$ using the exponential function $Y_{i}=\alpha e^{\beta X_{i}}+\varepsilon$, where $\varepsilon$ is the error in the estimate. The parameters $\alpha$ and $\beta$ are estimated by taking the logarithm of the values $Y_{i}$ in the series, and then using Simple Linear Regression to find the linear relationsship between $X_{i}$ and $\ln(Y_{i})$. Expressions for $\alpha$ and $\beta$ are given in [14].

Once we know the parameters, we can compute the average squared error by $R=\frac{1}{N}\sum_{i=1}^{N}\left(\alpha e^{\beta X_{i}}-Y_{i}\right)^{2}$.

III Proposed method

The two main parts of the proposed method are described in detail below; in short, the purpose of the first part of the method is to provide numerous estimates of modes using dimensionality reduction techniques, while in the second part averaged estimates of modes and mode shapes are computed using clustering. We assume $M$ frequency measurement streams from which we want to extract information about modes and mode shapes.

Part I - Mode estimation

The first part of the method acts on a sliding time window containing $N$ samples (for each of the $M$ series). The length of the sliding window should be sufficiently long to capture low frequency oscillations, i.e. in the range $5-$10\text{,}\mathrm{s}$$.

For each time window, we would like to decompose the oscillatory behaviour into a number of components, and find the amplitude and phase with which this behaviour happens in each measurement. Each time window is analyzed according to the following sequence:

The two-layer combination of PCA and CPCA, descibed in section II-D, is applied to the $M\times N$ matrix of input series:

1. (a)

Each of the time series are normalized such as to have zero mean. 2. (b)

PCA is applied, 3. (c)

The scores are detrended, 4. (d)

The Hilbert Transform is applied, and 10% of the ends of the resulting complex series are tapered 5. (e)

CPCA is applied. 2. 2.

The frequency and decay rate of the resulting CPCs are estimated. The frequency is computed using the built-in MATLAB-function meanfreq(). The decay rate is computed using exponential regression on the absolute value of the CPCs ($\mathrm{Abs}(\mathbf{z}_{i})=\sqrt{\mathrm{Re}(\mathbf{z}_{i})^{2}+\mathrm{Im}(\mathbf{z}_{i})^{2}}$), described in II-F. 3. 3.

To limit the amount of uninteresting observations and noise in the clustering, only the estimates with a frequency between $0.1\text{\,}\mathrm{H}\mathrm{z}$ and $2\text{\,}\mathrm{H}\mathrm{z}$ (in the range of electromechanical oscillations) and a regression error below $4\cdot 10^{-3}$, computed as described in section II-F, are used in the further analysis. The CPCs passing this criterion, along with the corresponding coefficients describing the amplitude and phase of their presence in the measurements, are expected to be relatively good estimates of modes and mode shapes. 4. 4.

Each of the accepted CPCs are rotated such that the largest coefficient, corresponding to the longest observability phasor in the mode shape, is at $0\text{\,}\mathrm{\SIUnitSymbolDegree}$. 5. 5.

Finally, each mode estimate and the corresponding mode shape is stored as a $2M+2$-dimensional point on the form

[TABLE]

where

[TABLE]

Once we have the mode estimates, it does not matter from which CPCs they were derived. Assuming that CPC $j$ fulfilled the criterions in step 3, then the observability phasor of measurement $i$ would be given by the coefficient $w_{ij}$ in (19), i.e., we have $G_{i}=w_{ij}$. In the following, the stored points on the form (20) are referred to as ”observations”, to avoid confusion with the final, averaged mode estimates resulting from Part II, which are referred to as ”mode estimates”.

Step 3 in the sequence above is included to filter out bad observations, due to the k-Means algorithm being susceptible to noise. If a CPC resulting from the two-layer PCA-decomposition successfully captures a mode, then the absolute value of the CPC should decay (or grow) exponentially. Applying regression as described in section II-F gives an estimate of the decay rate of the mode, along with the average error in the estimate. A large error implies that the CPC does not decay or grow exponentially, and therefore is probably not a good estimate of one of the modes of the system. The threshold value of $4\cdot 10^{-3}$ is found empirically, resulting in reasonable balance between acceptance and rejection of observations.

Part II - Computing mode averages

The PCA-decomposition in Part I is able to produce snapshots of modes and mode shapes, but with some uncertainty. To digest the numerous observations and provide more reliable estimates, we would like to compute averaged modes- and mode shapes.

To determine which of the observations to include when averaging each mode, we apply a clustering algorithm: Observations that are similar are be grouped together, and the different groups/clusters are assumed to be associated with different modes. The averaged mode shapes are then computed as the centroids of the clusters.

Assuming we have $Q$ point estimates resulting from Part I, these can be organized in a $Q\times(2M+2)$-matrix:

[TABLE]

Clustering of the points is carried out using the k-Means algorithm as described in II-E: The data is clustered multiple times with different numbers of clusters (for instance 1 to 10 clusters), and the final clustering is chosen to be the one with the highest silhouette index.

Finally, the averaged modes are derived from the centroids of the clusters.

IV Application

The method is tested on simulated data from The Kundur Two Area System, performed in DigSILENT PowerFactory, and on PMU data recorded during an oscillatory event in the Nordic Power System. In the case with simulated data, the analysed time series has very high modal content and no ambient noise, and should be regarded as a demonstrational example. The case with PMU data contains ambient noise, and demonstrates the applicability of the method to measurements from a real system.

IV-A Simulated data from the Kundur Two Area System

This system is widely used for studying local and interarea modes. The system parameters are described in [9], but the model is modified according to the following: Generators 1 & 2 emulate hydro power plants, with governor model HYGOV, inertia constants $H=$10\text{,}\mathrm{s}$$ and synchronous reactances $X_{d}=1.2$ & $X_{q}=0.9$. Generators 3 & 4 emulate thermal power plants, with governor model TGOV1. All generators are equipped with the excitation system model SEXS. Finally, the load flow is modified, in order to achieve lower damping of the interarea mode, such that Generators 1, 2, 3 and 4 produce $600\text{\,}\mathrm{M}\mathrm{W}$, $500\text{\,}\mathrm{M}\mathrm{W}$, $801\text{\,}\mathrm{M}\mathrm{W}$ and $900\text{\,}\mathrm{M}\mathrm{W}$, respectively. Other parameters are as described in [9].

Short circuits with a clearing time of $20\text{\,}\mathrm{m}\mathrm{s}$ are applied subsequently near each of the four generators, with $20\text{\,}\mathrm{s}$ between each short circuit. The resulting frequencies (measured at the generator buses), are shown in Fig. 1.

Figs. 2 and 3 show histograms resulting from Part I of the method, indicating the number of observations with a given frequency, damping and observability phasor for each generator. Considering the histogram for frequency in Fig. 2, we see that higher densities of observations form around frequencies $0.3,0.5$ and $0.7\text{\,}\mathrm{H}\mathrm{z}$.

Part II of the method performs the clustering of the densely populated areas, and results in 7 clusters. The lower plot in Fig. 1 shows the time instants where observations assigned to the different clusters are obtained.

The coloring of the columns of the histograms indicate the share of the observations being assigned to different clusters, for instance, most of the $\approx$0.5\text{,}\mathrm{H}\mathrm{z}-observations are assigned to the cluster shown in blue, and most of the $\approx$0.7\text{\,}\mathrm{H}\mathrm{z}-observations are assigned to the cluster shown in red, and so on. The same colors are used in the lower part of Fig. 1 and Figs. 2 to 4 to indicate the clustering division. Considering the observations around $0.5\text{\,}\mathrm{H}\mathrm{z}$ in Fig. 2, shown in blue, we can study Fig. 3 to find that in these observations, the observability phasors for generators 1 & 2 (shown in the two upper plots) in average appear to have a real part around -0.4, and an imaginary part of about zero. For generators 3 & 4, the corresponding real and imaginary parts appear to be situated around 0.7 and zero. This indicates that this is an interarea mode, where generators 1 & 2 are swinging against generators 3 & 4. Considering the lower plot in Fig. 2, the decay rate of the same mode appears to be less accurately estimated due to the high variance compared to the other histograms, but the average appears to be around $0.3\text{\,}\mathrm{1}\mathrm{/}\mathrm{s}$.

Once the clustering is performed, we know which observations to include when averaging each mode. The averaged modes from three of the clusters are shown in the upper part Fig. 4. For comparison, the corresponding modes computed using modal analysis in DigSILENT PowerFactory are shown below, indicating that the method is capable of providing accurate estimates of system modes. The most significant deviation is the decay rate of the third mode shown, where the proposed method produces an estimate of $-0.37\text{\,}\mathrm{1}\mathrm{/}\mathrm{s}$, while modal analysis gives the value $-0.939\text{\,}\mathrm{1}\mathrm{/}\mathrm{s}$. Note that this mode is nevertheless very well damped, and that it is a mode associated with turbine governors. The method is thus not limited to identification of electromechanical modes.

As mentioned, the method results in 7 clusters, representing 7 averaged modes, whereof three are shown in 4. By comparing with results from modal analysis, it is found that two of the averaged mode estimates, number 6 and 7, do not resemble any of the modes from modal analysis.

Further, it is found that mode estimate 2 and 3 are very similar, except that the mode shape is rotated $180\text{\,}\mathrm{\SIUnitSymbolDegree}$, which is probably due to the mode shapes from Part I being rotated such that the longest observability phasor lies at $0\text{\,}\mathrm{\SIUnitSymbolDegree}$. Due to some inaccuracy in the CPCA decomposition, the next longest observability phasor will in some estimates appear the longest, thus flipping the whole mode shape around. The same is the case for estimates 1 and 5, the averaged mode shapes are very similar, but appear at different angular positions.

IV-B PMU measurements from oscillatory event

The recorded PMU data is shown in the upper part of Fig. 6. Applying the method yields the three mode shapes shown in Fig. 8. Studying Fig. 7, which shows a zoomed view at the beginning of the oscillatory period, we see that the oscillations are most prominent in the two signals shown in green and yellow, corresponding to locations S and SW1 in Fig. 5. These locations are oscillating roughly in phase. Considering locations SW2 and SE, the amplitudes of the oscillations are a bit lower, and the respective phase shifts are about $135\text{\,}\mathrm{\SIUnitSymbolDegree}$ to the left, and $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ to the right. Mode Estimate 1 fits well with this description, indicating that the method is able to produce a reasonable estimate of the mode and corresponding mode shape.

From the lower plot in Fig. 6, it can be seen that the observations from which Mode estimate 1 is computed, are captured shortly after the oscillations started, Mode estimate 2 was captured shortly after that, and and finally Mode estimate 3 was captured at the end of the period with oscillations. The three mode shapes are similar, but in particular the angle of the phasor corresponding to location SE in Fig. 5 (shown with dashed, red lines) changes during the course of the oscillations. Also, the frequencies of Estimates 2 and 3 are a bit higher than the first. This could indicate that some remedial action was taken to mitigate the oscillations.

V Discussion

The above results are promising; first and foremost, the method is able to produce relatively accurate estimates of modes of the power system. In the simulated case, the mode estimates are very similar to analytically calculated modes. In the case with PMU data, the modes are as expected when considering amplitude and phase of oscillations in the time series.

The results from the case with PMU data reveals another potential benefit of the described empirical approach, namely that the method could facilitate tracking of how modes changes during a period. In the presented case, the oscillatory mode and the corresponding mode shape appears to change slightly during the course of the oscillations. This is information that would be difficult to obtain using modal analysis, which could potentially contribute to increased situational awareness.

When applying the method to the Kundur Two Area system, two of the resulting mode estimates are copies of other estimates, appearing at different angular positions. This is caused by the way the observations are stored, i.e. with the longest observability phasor at $0\text{\,}\mathrm{\SIUnitSymbolDegree}$. This could be avoided in a post processing step, checking if any of the mode estimates are replicates of each other, or by storing the observations in a format that is independent of the angle of the mode shape.

Further testing needs to be done to improve the method, and facilitate development into a form that is useful for grid operators. Imaginably, the dynamics decomposition in Part I would run continuously, analyzing numerous PMU data streams on the go, while the clustering in Part II could be performed on a longer-term moving window, for instance of lengths on the scale of minutes or hours. Based on some measure of certainty of each mode estimate (e.g. based on the degree of variance, or the number of observations in the cluster from which it is computed), the estimates could be stored, along with the time of the day when they were obtained, which could further be used to indicate how modes and mode shapes change during the course of load variations throughout a day. This could in turn be taken into account when designing future wide area control applications.

For the implementation, there is a lot of potential for improvement. Other input quantities than the frequency could be tested as input to the method, for instance the voltage angle. Further, in Part I of the method, the criterion for determining whether observations should be included in the clustering part could take many forms. If a clustering algorithm less susceptible to noise was used, the criterion could be less strict, or even be omitted. Other variants of the dynamics decomposition could also be tested, for instance as described in [7], where EMD is used to decompose the input signals into IMFs rather than the two-layer combination of PCA and CPCA implemented here.

VI Conclusion

Through testing on simulated data and recorded PMU data, the presented method appears to be a useful tool for digesting numerous PMU measurements into useful knowledge of the system, in the form of information about the oscillatory modes and their mode shapes. The main conclusion is that the method is feasible, and shows promising result on measured and simulated data. Further work is required to improve the different constituents of the method, including the dynamics decomposition using CPCA, the filtering/handling of noise/bad observations and the clustering part. Further work on measured data is necessary to verify the robustness and efficiency of the method for online use.