Identification of Low-Frequency Oscillation Modes Using PMU Based Data-Driven Dynamic Mode Decomposition Algorithm

Power system inter-area oscillations curtail the power transferring capabilities of the transmission lines in a large interconnected power system. Accurate identification of dominant modes and associated contributing generators is important to avoid power system failures by taking appropriate remedial measures. This paper proposes a multi-channel Improved Dynamic Mode Decomposition (IDMD) algorithm-based modal analysis technique using Synchrophasors measurement. First, a reduced-order dynamic power system model is estimated and using this model dominant oscillation modes, corresponding modes shapes, damping ratio, coherent group of generators, participation factors are determined. To improve the accuracy data stacking technique is used to capture detailed information of the system. An optimal hard threshold technique is utilized to select the most optimal model order to avoid uncertainties due to the presence of high level of measurement noise. The study results show that the proposed algorithm gives an accurate and robust solution even in systems having high level of noise in the measurement data. The performance of the proposed technique is tested on simulated data from two-area four-machine system and wNAPS 41-bus 16-generator system with PMU measurements corrupted with different levels of measurement noise. To further strengthen the viewpoint, the proposed method is validated on real-time PMU measurement from ISO New England data to validate the accuracy of the proposed work.


I. INTRODUCTION
Low-frequency electromechanical oscillations exist inherently in any large interconnected power system. These oscillations can be easily observed even with normal power system operation and may be excited due to switching operations, generation and load variations, faults, etc., [1], [2]. The accurate information about these low-frequency oscillation modes is crucial to determine the characteristics of the power system. These oscillations can be broadly classified as interarea, local, limit cycle and forced oscillations. Out of these, the damping of inter-area oscillation is most important as it affects the stability of large interconnected power system and limits the power-transferring capacity of the transmission lines exchanging power between two areas [3], [4].
The associate editor coordinating the review of this manuscript and approving it for publication was Amin Hajizadeh.
Conventionally, the analysis of electromechanical oscillations was conducted by developing the power system dynamic model from the set of differential equations and finding the modal analysis parameters by linearizing the system around the operating point [5]. This method may give inaccurate results, as the power system dynamics changes continuously, and it is almost impossible to update the algebraic model in real-time.
After various blackouts around the world, the Phasor Measurement Unit (PMU) based Wide Area Monitoring System (WAMS) is now widely accepted due to the numerous benefits it offers in real-time monitoring and control of highly stressed large interconnected power system. With the availability of these synchronized measurements, various signal-processing techniques were utilized to estimate the modal analysis parameters using synchronized PMU measurements such as voltage magnitude, voltage angle, frequency, active power reactive power etc. collected from multiple locations [6]- [8].
Various advanced signal-processing techniques are Fast Fourier Transform (FFT) [10], Matrix Pencil (MP) [11], [12], Prony Analysis (PA) [13], [14], Principal Component Analysis (PCA) [15], Empirical Mode Decomposition (EMD) [16], Multivariate Empirical Mode Decomposition (MEMD) [17], Hilbert Transform (HT) [18], wavelet transform(WT) [19], robust recursive least square (RRLS), regularized robust recursive least square (R3LS) algorithm [20], all these methods have their own merits and demerits. FFT is restricted to estimation of frequency components, while MP, PA and PCA estimate both frequency, damping, and the magnitude of the measurement data [21]. These methods are based on eigenvalue analysis of the system. They, however do not estimate mode shapes and participation factor, which are the most critical parameters to characterize low-frequency oscillations. In addition, the performance of these methods is sensitive to the presence of measurement noise and signal offset. When the PMU measurements are subjected to high level of measurement noise, which is unavoidable due to defects in measurement and communication devices or may be due to extreme operating conditions that includes random load fluctuations. These conditions affects the estimation accuracy in two ways. First, with noisy PMU measurements, results may contain significant errors due to the presence of spurious frequency components due to the presence of noise. Second, it is quite challenging to determine the accurate model order of the system for estimation of low frequency oscillation modes that may also introduce significant errors [22], [23].
Therefore, noise consideration should be also taken in to account for characterization of low-frequency oscillations. This problem is overcome by EMD technique that processes the raw PMU data using filter banks. However, it can lead to the phenomena of mode mixing, leading to many irrelevant artificial modes [16]. HT and WT analyze the oscillation features in both time and frequency domain. This magnifies the identification of the onset of any event. These algorithms, however, are only applicable to systems having a single dominant mode [19]. RRLS and R3LS are based on the identification of autoregressive moving average (ARMA) model to estimate dominant oscillation modes from ambient PMU data [20], [24]. The aforementioned methods also have one common limitation when directly applied, as denoising method is the estimation of an accurate reduced-order model of the system for estimation of dominant low-frequency oscillation modes. These methods are known as eigenvector analysis methods as they estimate mode shapes from right eigenvectors estimated from state matrix A of the estimated reduced power system model.
Recently eigensystem realization algorithms (ERA) and Dynamic mode decomposition (DMD) algorithms have proved very accurate in identification of eigenvalues of a dynamic system. In DMD, both state matrix A and state variables are known, unlike ERA where only state matrix A can be determined. This gives an additional advantage to DMD over ERA as it can provide additional information about the generators contributing to inter-area mode [25].
The DMD provides a spatiotemporal decomposition of a high-dimensional dataset of a given system in time. This advantage is not offered by any other model identification method. The system model can be estimated utilizing the dominant spatiotemporal modes obtained from given measurement matrix with reduced system model order obtained from singular value decomposition. The future state can also be predicted using initial system conditions and obtained a reduced-order system. These advantages extend the use of DMD in various power system applications [26].
Ref. [27] utilizes the DMD algorithm based on energy metric for analysis of power system disturbances. Ref. [22] utilizes an optimized DMD algorithm to extract low-frequency oscillation modes from wide-area measurement datasets using variable projection and finite difference style approximation method. In [28] DMD algorithm is utilized to estimate the frequency and amplitude of measurement data on different power system scenarios with different levels of measurement noise. Ref. [29] proposes a randomized DMD algorithm that reduces the computational cost by radically reducing the size of the measurement matrix to compute the low-frequency oscillation mode. In [30], DMD algorithm is implemented on frequency and voltage measurements to determine oscillation frequency, damping ratios and the amplitude of measurement. Ref. [31] discusses koopman based DMD algorithm for determination of modal analysis parameters of the power system. Ref. [18] utilizes the DMD framework to interpret the global dynamic behavior of a wide-area system in terms of spatial and temporal patterns containing a single tone frequency component. It highlights the superiority of DMD over koopman and prony methods in identifying dominant oscillation modes. Ref. [32] compares the performance of DMD algorithm with ERA considering three power system applications that include RLC circuit measurement, modal analysis using PMU measurements and power signal analysis corrupted with harmonics.
It can be noted that in all the above methods, the effect of measurement noise has not given sufficient consideration. Although measurements in real-time systems are often corrupted with, high level of measurement noise and WAMS has no exception. Therefore, the accuracy of all the methods available in the literature is affected due to the presence of measurement noise that makes the extraction of dominant frequency modes more challenging.
The contributions of this research work are highlighted as follows: 1. Utilization of data stacking technique to increases the dimensions of the measurement matrix. This increases the number of estimated eigenvalues that help in identifying the detailed dynamics of the low-frequency oscillation modes for better characterization of system dynamics. VOLUME 9, 2021 2. Application of hard threshold technique on singular values obtained from SVD to obtain judicious rank truncation and denoised DMD eigenvalues to accurately capture the system dominant modes. 3. Extend the application of DMD from dominant mode estimation to the determination of mode shapes, coherent group of generators and participation factor. Hence, develop a systematic approach to reveal the dynamic characteristics of the electromechanical oscillation modes in one holistic framework. This will prove to be very helpful in identifying the generators contributing to inter-area modes and take corrective actions to improve the damping of these modes. This paper is structured as follows. In section II, the mathematical background of DMD algorithm is briefly reviewed along with the concept of data stacking and hard threshold techniques for rank determination. Section III discusses various modal analysis parameters required to carry out the analysis. Results obtained from considered test systems are presented in Section IV and finally, conclusions are drawn in Section V.

II. DYNAMIC MODE DECOMPOSITION
DMD is one of the potential algorithms to obtain a reducedorder model of any nonlinear system. It has two forms that can achieve the mathematical formulation of any system. First is by assuming the linear dependency beyond the critical number of snapshots and introducing a companion matrix to approximate the remaining infinite number of irrelevant snapshots using linear koopman operator. The second is obtaining the similarity matrix from the constructed observation matrix and components of the singular value decomposition obtained from snapshot matrix. The latter method can provide the best approximation of the system described by a large dataset given that two aspects are satisfied. First is expressing the few channel measurement data matrix into a highresolution row dimension matrix to obtain large eigenvalues from the measurement matrix to better approximate the system. The second is the accurate determination of optimal rank r of the system to construct the similarity matrix for reducedorder model of the system [26], [32].
A mathematical definition of DMD is stated as follows: For any linear time-invariant dynamic system, the present state z k can be correlated to next state z k+1 as where, z R n , with sampling interval t.
The time-domain expression of z(t) can be conveniently derived if the eigenvalues of A and its eigenvector matrix are known.
where is the diagonal matrix having continuous eigenvalues. In discrete form, it can be written as: The data measurement matrix is formed by data gathering into subsequent overlapping sets with other being timeshifted as The modal analysis procedure can be carried out using these two sequential time-shifted measurement data matrices. Based on the above background, the key steps of the DMD algorithm for processing PMU datasets are as follows: 1. Compute the SVD on the matrix as obtained in 4(4a) as where U and V are orthonormal 2. Express the matrix A in terms of SVD components and observation matrix Z 2 as Construct a similarity matrix based on the selected truncated rank r obtained from constructed singular values as where A R r×r Ã represent a similar feature as that of A but with dimensionally reduced matrix discarding the irrelevant features. 4. The lower rank dynamic model is defined as are the eigenvectors and eigenvalues respectively of the similarity matrix. 5. Compute the reconstruction matrix of DMD as i W i is a matrix consisting of ith DMD modes and b is a matrix consisting of the initial amplitude of each mode. The desired oscillation modes, damping ratio, mode shapes, and participation factors can be determined using the eigenvalues and eigenvectors obtained. The concept of data stacking and optimal hard threshold technique is explained in the following subsections: A. DATA STACKING Consider one-channel measurements Z represented by a row vector. This will give only one eigenvalue of the system that cannot give the best approximation of the dynamics of the system. To address this issue data stacking technique using the hankelization concept is adopted to create time-shifted copies of datasets. This increases the row dimensions of the measurement matrix to obtain the large eigenvalues of the system to capture the detailed system dynamics [28], [29].
The generalized form of the augmented matrix with stacking number s can be written as The overlapping set of the augmented matrix Z aug can be formed as The value of s depends on the number of frequency components present in the signal. A higher value of s should be selected for signal corrupted with a high level of measurement noise. In this study, the value of s varies from 200 to 450.

B. OPTIMAL TRUNCATION OF SINGULAR VALUES USING HARD THRESHOLD TECHNIQUE
In the DMD algorithm, nonsingular values are arranged in decreasing order while computing the SVD of the signal matrix Z. Thus keeping useful signal information represented by first k singular values and discarding the remaining (s-k) part representing the noise. The selection of k is quite trivial, and difficulty lies in finding the optimal value especially when the measurement matrix is corrupted with high level of measurement noise. To solve this issue, the authors have utilized hard singular value thresholding technique proposed in [33] to select the optimal rank r for measurement matrix Z containing useful signal information Z s plus additive white noise Z n .
This approach estimates the optimal threshold value τ of the singular values as given by τ = ω (β) σ median (8) where, σ median is the median singular value of the noisy measurement matrix and and µ B is the marcenko-pastur distribution. Thus, useful signal information represented by k singular values is kept while the remaining (s-k) part representing the noise is discarded. The performance of the proposed algorithm is also verified using ERA algorithm briefly described in Table 1.

III. MODAL ANALYSIS PARAMETERS
Modal analysis is a systematic methodology of determining the characteristic of the power system in terms of frequency, damping ratio, mode shape and participation factor. These parameters should be carefully examined to investigate the dynamic behavior of the power system under different characteristic frequencies called modes [35].

A. FREQUENCY AND DAMPING RATIO
After estimating the reduced-order similarity matrix, Ã, the frequency and damping ratio of the system can be calculated using eigenvalues of the system matrix Ã(λ i = a i + jb i ). The frequency f i and damping ratio ξ i of i th oscillation mode in the system and can be expressed as φ i is called the right eigenvector of the system matrix A associated with the eigenvalue λ i . It is also called a mode shape and determines the contribution of the system state (machine) in each dominant mode. It also gives information regarding the behavior of oscillations against each other, i.e., inter-area, local, forced etc. The length of the mode shape determines the contribution of the system state in that mode. For every eigenvalue λ i , there exists an eigenvector ψ i , which satisfies ψ i is called the left eigenvector of the system matrix A associated with the eigenvalue λ i . These along with the initial system states determine the amplitude of oscillations.

C. PARTICIPATION FACTOR
The estimation of the relationship between states and the modes utilizing right and left eigenvectors is dependent on units and scaling associated with the state variables. To solve this problem, a matrix called participation factor is proposed in [36] obtained by multiplying the right and left eigenvectors which is independent of the choice of the units and provides an improved alternative to associate the state and the modes of the system. Using the estimated φ i and ψ i , the participation factor P can be calculated as With The element p ki = φ ki ψ ik is called the participation factor, which gives an overview of the contribution of the system generators in each of the dominant modes.
This value is normalized by dividing the obtained values with the maximum participation factor as follows:

IV. TEST CASES
The performance and applicability of the proposed algorithm are tested utilizing simulated PMU data from two-area fourmachine system and wNAPS 16 generator 41-bus system. The obtained results are also validated on real-time PMU measurements from ISO New England dataset.

A. TWO-AREA FOUR MACHINE SYSTEM
The Data for inter-area oscillation is generated using two-area four-machine model using Matlab-based dynamic simulation power system toolbox. The machines are modeled with static excitation system having conventional PSS with twocascaded lag-lead compensation block. As represented in Fig. 1, the network is divided into two Areas a 1 -(G 1 -G 2 ) and a 2 -(G 3 -G 4 ). The PMUs are installed at highly stressed buses i.e. bus 7 in Area 1 and bus 9 in Area 2. A detailed description of the system is given in [37]. The sampling rate of PMU is taken as 60 samples per second.
The inter-area oscillation cases are produced by increasing the stress on the system along with the introduction of the dynamic event. This is done by increasing the equivalent load at bus 7 to 1000 MW and at bus 9 to 2000 MW. The dynamic event considered in the present case study is a single-line to ground fault, applied to the line between bus 7 and bus 9.
Frequency signals, one from each Area (Bus 7 in Area 1 and Bus 9 in Area 2) are selected for modal analysis parameter estimation and signal reconstruction. A data window of 10 seconds is used for the construction of the augmented Hankel matrix, and the value of stacking factor s varies from 200 to 450 for all the cases considered.
The result of modal analysis for a different level of additive noise is given in Table 2. In case of higher level of additive noise, the higher model order r is required for obtaining high correlation between actual and DMD spectra. For noiseless data, the value of r=16 gives accurate results. With the additive noise of 2%, the optimal threshold value τ obtained from SVD of noise is 7.125. Therefore, rank r for the measurement matrix Z having singular values larger than the obtained τ is 23. With additive noise of 5%, the value of τ is 5.67 and the corresponding value of rank r is 26. For 8% of additive measurement noise, the value of τ is 5.13 and the corresponding value of rank r is 31 as shown in Fig. 2.
For all the test cases considered with different levels of measurement noise, the reconstructed denoised data obtained from the reduced-order system is shown in Fig. 3. It can be concluded that the reduced-order model obtained from the proposed method can accurately recover the actual PMU data from noisy PMU measurement.
From Table 2, it can be concluded that even with the increase in the level of measurement noise, the frequency and damping ratio error is within the acceptable limits. This proves the accuracy and robustness of the proposed methodology in identifying the dominant oscillation of the wide-area system.
Using the estimated left and right eigenvectors, the mode shapes associated with four inter-area modes having a damping ratio of less than 10% are plotted as shown in Fig. 4. Mode shape 1 reveals that the generator G 1 in Area 1 oscillates against generator G 4 in Area 2, while generator G 1 has a significant role in this oscillation. Mode shape 2 demonstrates that the generators G 1 and G 2 in Area 1 and generators G 3 and G 4 in Area 2 are coherent and oscillate against each other. While the generators G 1 and G 4 have a significant role in this oscillation. Mode shape 3 reveals that G 1 and G 3 oscillate against each other and exhibit higher participation in this oscillation. Mode 4 reveals that Generator G 1 in Area 1 oscillates against generator G 4 in Area 2 while G 1 exhibits high participation in this oscillation mode.
The accuracy of mode shapes estimated by IDMD algorithm is verified using ERA algorithm, as shown in Fig. 5. The phase displacement expressed in mode shapes for all the generators determined by both the algorithms may vary but    can accurately identify the contribution of generators in each mode.
The participation factors corresponding to dominant inter-area modes are estimated via (13) and shown in Fig. 6. It can be observed that Generators G 1 has a higher participation factor in mode 1 while G 4 dominates in oscillation mode 2. For mode 3, both generators G 1 and G 3 possess significant participation factor while  the generator G 1 has the highest participation factor in mode 4.
The participation factor estimated by the proposed improved DMD algorithm is also compared with ERA algorithms. It can be concluded from the obtained results that the participation factor obtained from both the algorithms is almost similar.
The result of modal analysis for a different level of additive noise is given in Table 3. For noiseless data, the value of r=27 gives accurate results. With the additive noise of 2%, the optimal threshold value τ obtained from SVD of noise is 12.436. Therefore, rank r for the measurement matrix Z having singular values larger than the obtained τ is 33. With additive noise of 5%, the value of τ is 9.213, and the corresponding value of rank r is 41. For 8% of additive measurement noise, the value of τ is 7.32 and the corresponding value of rank r is 47 as shown in Fig. 7.
For all the test cases considered with different levels of measurement noise, the reconstructed denoised data obtained from the reduced-order system is shown in Fig. 9.
It can be concluded that the reduced-order model obtained from the proposed system can accurately recover the actual PMU data from noisy PMU measurement. Table 3 shows that the percentage frequency and damping ratio errors for all the cases considered are within the acceptable limits and it proves the robustness of the proposed methodology against measurement noise.
The mode shapes of inter-area modes of all the generators estimated from the proposed method are shown in Fig. 10. As can be observed from Fig. 10, the mode shapes associated with mode 1 suggest that generators G 4 , G 6 , G 12 , G 14 in Area 2 and generators G 3 and G 11 in Area 3 are coherent and oscillate against generators G 7 , G 8 , G 15 and G 16 in Area 4 at a frequency of 0.23 Hz. For mode 2, generators comprising of G 1 , G 2 , G 3 , G 4 , G 5 , G 6 , G 12 and G 14 in area 1 and 2  are coherent and oscillates against generators G 3 , G 7 , G 8 , G 11 , G 15 , G 16 in Area 3 and 4 respectively at a frequency of 1.29 Hz. For mode shapes associated with mode 3, the generators G 3 and 11 in Area 3 oscillate against generators G 7 , G 15 , G 16 in Area 4 at a frequency of 0.89 Hz. For mode 4, generator G 1 , G 2 , G 4 , G 5 , G 6 , G 9 , G 10 , G 13 , G 12 and G 14 in Area 1 and 2 form a cluster of coherent group and swings against the cluster of other coherent groups of generators comprising of generators, G 3 , G 7 , G 8 , G 11 , G 15 and G 16 in Area 3 and 4 respectively at a frequency of 1.47 Hz.
The mode shapes estimated by the ERA, as shown in Fig. 11, are in close proximity to the mode shapes estimated by the proposed algorithm.
The participation factor of all the generators in each dominant mode is shown in Fig. 12. It can be observed from Fig.12 that generators G 14 and G 15 exhibit higher participation factor in mode 1 while G 15 plays a significant role in mode 2. For mode 3, the two largest participation factors in this inter-area mode are from generator G 3 and G 7, while G 6 and G 7 perform higher participation in inter-area oscillation mode 4.
The participation factor estimated by the proposed improved DMD algorithm is also compared with ERA algorithms. It can be concluded from the obtained results that the participation factor obtained from both the algorithm is almost similar.

C. ANALYSIS OF OSCILLATION ON ACTUAL PMU DATA
The performance of the proposed algorithm is also tested utilizing real-time PMU measurements from ISO New England. The system is divided into 3 areas as shown in Fig. 13.  The dataset contains PMU measurements from 12 substations of area 1, while line 7 and line 21 connect the external areas 2 and 3 [39]. In this work, the authors have utilized datasets reported on July 20, 2017, to test the effectiveness of the proposed algorithm.
Frequency signals from substations 3, 6, 7 and 9 are selected for modal analysis parameter estimation and signal reconstruction. A data window of 10 seconds is used to construct the augmented Hankel matrix, and the value of stacking factor s varies from 200 to 450.
The results of modal analysis is given in Table 4. It is found that even for actual PMU data, the proposed algorithm is capable of recovering the actual PMU data from noisy PMU measurements. The value of most optimal rank r is 21 for obtaining high correlation between actual and DMD spectra and the corresponding value of τ obtained from SVD of noise is 11.35 as shown in Fig. 14.
The reconstructed denoised data obtained from the reduced-order system is shown in Fig. 15. It can be concluded that the reduced order model obtained from the proposed   system can accurately recover the actual PMU data from noisy PMU measurement.
From Table 4, it can be concluded that the proposed method is capable in identifying the dominant oscillation modes of the wide-area system with measurements from a real-time system. Using the estimated left and right eigenvectors, the mode shapes associated with four inter-area modes having damping ratio of less than 10% are plotted as shown in Fig. 16. The mode shape 1 reveals that the generator G 1 in Area 1 oscillates against generators in Area 2 and Area 3 respectively, while the generator G 1 has a significant role in this oscillation. Mode shape 2 demonstrates that the generators in Area 1 and Area 2 are coherent and oscillate against generators in Area 3, while the generator G 1 and generators in Area 3 have a significant role in this oscillation. Mode shape 3 reveals that generators in Area 1 and Area 3 are coherent and oscillate against generators in Area 2, while generators in area 2 exhibit higher participation in this oscillation. Mode 4 reveals that Generator G 1 and G 2 in Area 1 oscillate against Area 2 and Area 3 respectively, while G 1 exhibits high participation in this oscillation mode.
The accuracy of mode shapes estimated by IDMD algorithm is verified using ERA algorithm, as shown in Fig. 17. The phase displacement expressed in mode shapes for all the generators determined by both the algorithms may vary but can accurately identify the contribution of generators in each mode.
The participation factors corresponding to dominant interarea modes are estimated via (13) and shown in Fig. 18. It can be observed that Generators G 1 in Area 1 has a higher participation factor in mode 1 while generator G 2 in Area 1 and generators in Area 3 dominates in oscillation mode 2. For mode 3, generators in Area 2 possess significant participation factor while the generator G 1 in Area 1 has the highest participation factor in mode 4.
The participation factor estimated by the proposed improved DMD algorithm is also compared with ERA algorithms. It can be concluded from the obtained results that    the participation factor obtained from both the algorithms is almost similar.

V. CONCLUSION
This work presented an Improved DMD algorithm to estimate the entire modal analysis parameters like dominant modes, damping ratio, mode shapes, coherent group of generators, and participation factor in one holistic framework. The introduction of data stacking technique and optimal hard threshold technique for determination of reduced model order substantially improves the performance of standard DMD algorithm in identifying the system dominant modes even in the presence of high level of measurement noise.
The proposed approach is tested on two-area four-machine system and wNAPS 41 bus 16 generator system with PMU measurements corrupted with different levels of measurement noise. The obtained results are also validated on real-time PMU datasets from ISO New England. It can be observed from the results obtained that even in the presence of high level of measurement noise the proposed algorithm can accurately recover the actual PMU measurements and estimate the dominant modes of the system with high accuracy.