A novel method for locating the source of sustained oscillation in power system using synchrophasors data

Large interconnected power systems are usually subjected to natural oscillation (NO) and forced oscillation (FO). NO occurs due to system transient response and is characterized by several oscillation modes, while FO occurs due to external perturbations driving generation sources. Compared to NO, FO is considered a more severe threat to the safe and reliable operation of power systems. Therefore, it is important to locate the source of FO so corrective actions can be taken to ensure stable power system operation. In this paper, a novel approach based on two-step signal processing is proposed to characterize FO in terms of its frequency components, duration, nature, and the location of the source. Data recorded by the Phasor Measurement Units (PMUs) in a Wide Area Monitoring System (WAMS) is utilized for analysis. As PMU data usually contains white noise and appears as multi-frequency oscillatory signal, the first step is to de-noise the raw PMU data by decomposing it into a series of intrinsic mode functions (IMF) using Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN) technique. The most appropriate IMF containing the vital information is selected using the correlation technique. The second step involves various signal processing and statistical analysis tools such as segmented Power Spectrum Density (PSD), excess kurtosis, cross PSD etc. to achieve the desired objectives. The analysis performed on the simulated two-area four-machine system, reduced WECC-179 bus 29 machine system, and the real-time power system PMU data set from ISO New England, demonstrates the accuracy of the proposed method. The proposed approach is independent of complex network topologies and their characteristics, and is also robust against measurement noise usually contained in PMU data.


Introduction
With the rapid deployment of Phasor Measurement Units (PMUs) in modern power systems, both natural oscillation (NO) and forced oscillation (FO) are now commonly observed. The presence of these oscillations introduces negative impacts such as reduced system reliability and equipment life, limited power-transferring capability and higher power losses, and sometimes can result in cascaded failure and even blackout [1][2][3].
The occurrence of NO is much more common compared to FO. However, in recent years many cases of FO have been observed, noticeable in Nordic power system, Chinese power grid, Indian power grid, western North American power system (wNAPS), US eastern interconnections (EI) [4][5][6][7]. NO can be efficiently suppressed by proper tuning of Power System Stabilizer (PSS) and generator rescheduling which substantially increase the damping ratio of the system to suppress poorly damped NO [8,9]. However, to counter FO, timely detection and isolation of the injecting source is required for safe and reliable grid operation. It is reported in literature that a poorly designed PSS can act as a source of FO [10]. In addition, with large-scale integration of renewable energy systems (RES) in the grid and continuous increase in variable loads with intermittent characteristics, the cases of FO have increased considerably. Among RES, wind energy systems are the major contributor of FO, primarily due to the fluctuation in wind speed and disturbances caused by doubly fed induction generator wind turbine (DFIG-WT) control strategies [11][12][13]. Compared to NO, FO initiates and spreads more rapidly, and has a much higher amplitude. Therefore, FO can produce more harmful effects than NO [14]. Several research works have been carried out to locate the source of FO in power system so corrective actions can be taken.
In [15], the idea based on a hybrid dynamic simulation approach using PMU data is presented. The simulation is run by replacing a portion of dynamic states and system variables with PMU measurements, and any deviations between simulation and PMU data locate the source of FO. Reference [16] proposes a method based on the flow of dissipation energy in the system. It demonstrates that network transient energy is equivalent to the energy dissipated by the damping torque, and the source dissipating abnormal transient energy is the source of sustained oscillation. Reference [17] utilizes Prony analysis to locate the source of oscillation using the damping torque coefficients of the generators, and the generator having a negative damping torque coefficient is considered as the source of oscillation. A two-step approach based on phasor and energy analysis is proposed in [18], and the realization of the total energy of the generator sets is obtained using a Hamiltonian function. The bus associated with the FO can then be detected according to the direction of the flow of injected energy obtained with phasor analysis. Reference [19] proposes the concept of oscillation phasors realized from oscillatory PMU dataset using Prony and first-wave algorithm, to identify the source of FO. The generator set having the leading phasor is identified as the source of FO. Reference [20] discusses a novel method based on the development of frequency response function for a particular generator set using PMU data and detection of oscillation source depends on the difference between the actual and estimated current spectra. Measurement noise is also taken into account to determine if the deviation between the actual and estimated current spectra is due to measurement noise or FO. An equivalent power system model in the frequency domain considering noisy PMU data and uncertain generator parameters is developed in [21] and a numerical procedure based on Bayesian framework is adopted to identify the source of FO. Reference [22] proposes a double stage mode decomposition technique, while [23] discusses an unknown input observer (UIO) based method to locate the source of FO. Initially, the linearized power system model considering the unknown oscillation sources from a prime mover, excitation system etc. is derived. The comparison between the model responses generated by UIOs and recorded responses establishes a strong relationship between the residual and oscillation sources. Hence, it can easily identify prime mover/excitation system of the generator set that gives rise to FO. A wave propagation-based technique is utilized in [24] to estimate the arrival time of various waves from different generation sources at reference bus using terminal voltages, and the wave associated with the leading generator is considered as the source of oscillation. Reference [25] proposes a distributed cooperative scheme utilizing phasor data concentrator (PDC) at the local control center to detect the source of FO. The energy exchange among local PDCs is realized using an average consensus algorithm, which determines the cut-set energy and locates the source of oscillation in a distributed manner.
Several other methods have also been reported in the literature, including graph theory-based methods [24], artificial intelligence-based methods [26], and equivalent circuit based method [27]. The summary of the performance of different methods for locating the source of FO is given in Table 1.
The contributions of this research work are highlighted as follows: Performance comparison of various methods in the literature for locating the source of FO. Utilization of ICEEMDAN algorithm to de-noise the raw PMU data to determine the actual low frequency oscillation modes presented in the system. Utilization of statistical analysis tools to accurately discriminate between NO and FO modes presented in the power system. Utilization of signal processing tools to characterize FO in terms of its frequency components, duration and location of source.
The rest of the paper is organized as follows. Section 2 describes various signal processing and statistical analysis tools used to formulate the proposed methodology, while Section 3 explains the systematic methodology adopted to characterize and predict the source of oscillation. Section 4 and 5 discuss the results obtained using simulated and real-time data sets, respectively. Finally, the conclusions from the presented work are drawn in Section 6.

Theoretical framework
This section presents a brief overview of various signal processing and statistical analysis tools utilized in analyzing the simulated and real-time FO cases.

ICEEMDAN
The raw PMU data is de-noised using ICEEMDAN algorithm [28]. The basic EMD algorithm cannot be utilized due to the problem of mode mixing [29], which refers to the presence of multiple frequency components contained in the IMFs after EMD decomposition. Numerous improved versions of EMDs have been reported in the literature. Ensemble EMD Table 1 Summary of performance of various methods available in the literature

Methods Advantages Disadvantages
Hybrid dynamic simulation [15] Independent of the type, model and complexity of the system Only applicable to post-fault analysis Dissipation energy method [16] Applicable with PMU measurements to determine oscillation energy flow.
Requirement of prior information about the model restricts its use in large interconnected system Damping torque method [17] Clear physical meaning, with PMU measurement applicable for implementation in wide area systems Often requires rotor speed, load angle and electromagnetic torque of generator, which are not directly available Phasor and energy analysis method [18] Robust to data quality issues, can detect oscillations from multiple sources Construction of energy function for a large interconnected power system is quite complex Oscillation phasors based method [19] Strong physical meaning as there is no oscillation phasor in case of stationary active power With the increasing trend of renewable energy integration into the grid, the detection of oscillation source becomes more challenging Frequency response function method [20], Bayesian Method [21] Free from load modeling assumptions, robust against measurement noise as well as uncertainties in generator parameters Detailed dynamic modeling is required, which affects the accuracy in case accurate models are not available Double Stage Mode Decomposition Method [22] Independent of the type, model and complexity of the system Relatively complex and requires high computational time due to signal decomposition in two stages, criteria for selection of most suitable IMF for analysis is missing, consideration of measurement noise is not taken into account Unknown input observer (UIO) based method [23] Easier to estimate state variables related to the source Detailed and accurate system modeling is required lacks a guide in choosing an appropriate amplitude of added noise to realize the signal, as it generates different mode numbers with different noise components [30]. CEEMDAN overcomes this shortcoming by adding a pair of positive and negative noise components to the actual signal, and hence, it tackles the problem of inconsistent mode number in the decomposition process. However, its decomposed modes consist of persistent residual noise that makes it difficult to denoise an already noisy signal [28]. ICEEMDAN algorithm overcomes the problem of mode mixing, inconsistent IMF number under variable noise, and requires lower number of sifting iteration per IMF with reduced computational burden. Therefore, it is the most suitable algorithm for decomposing nonstationary dynamic signals and is utilized in this work to de-noise the PMU data. The flow chart of ICEEMDAN algorithm is shown in Fig. 1. E k (.) is the k th IMF operator obtained by ICEEMD AN and M(.) computes the mean of the signal x [n].
For instance, is a realization of unit variance white Gaussian noise, and r k is the residue obtained at the kth realization. The algorithm is terminated when the residue fails to decompose as an IMF.

Cross correlation
The cross-correlation function is widely used to estimate the degree of similarity between two signals.
For finite energy systems, the cross-correlation of discrete-time signals x(n) and y(n) can be expressed as where l = 0, ± 1, ± 2…, is the window length varying from − ∞ ≤ l ≤ ∞. It is the same as the cross-covariance function for the zero-mean random process.

Segmented PSD
For a discrete random signal x(n), segmented PSD can be defined as a DTFT of an auto-correlation sequence r xx and is given by where and l = 0, ± 1, ± 2…, is the length of the window segment used.
The factors affecting PSD are overlapping segment, frequency resolution, and window function. The data is divided into N small segments of length ls, and the total number of segments has 50% overlapping between consecutive segments. There are various methods to estimate PSD such as yule walker, Welch, Burgs etc., while in this work, welch transform is used.

Cross PSD
For discrete random signals x(n) and y(n), cross PSD can be defined as a DTFT of a cross-correlation sequence r xy and is given by: where r xy is the cross correlation sequence as defined by (1).

Kurtosis
Kurtosis is the measure of deviation of the probability distribution from a normal distribution and can be  mathematically expressed as the fourth moment around the mean of the distribution as [31]: where X is the length of the input sequence, μ and σ are the mean value and variance of X, respectively. β is usually referred to as Pearson kurtosis. Low kurtosis indicates that data has flattop near mean while high kurtosis indicates distinct peaks near the mean.
As the value of kurtosis in case of normal distribution function is 3, to make the value of kurtosis zero, it is re-defined (known as excess kurtosis) as γ = β-3 as: Excess kurtosis is used in this study by Gaussian distribution as a reference. A window of 60 samples is considered and excess kurtosis is computed for each window.

Methodology
Based on the technical terms described in Section 2, the flow chart of the systematic procedure adopted to carry out the desired objectives is shown in Fig. 2. The entire analysis is divided into two stages. The first stage involves de-noising of raw PMU data using mode decomposition based on ICEEMDAN algorithm that transforms the signal into a series of IMFs known as decomposed modes after preprocessing (DMAPs). However, all the IMFs obtained may or may not represent the actual information contained in the signal. Therefore, the IMF selection methodology based on the correlation coefficient ( xy ) is utilized to select an optimal  IMF to carry out the analysis. The higher the absolute value of xy (τ) is, the higher its relative degree of similarity with the actual PMU data will be. The second stage includes determination of oscillation frequency, type of oscillation, its time duration, and source location (in case FO is detected). Segmented PSD using Welch transform is used to estimate the frequency of oscillation, while the time duration of the oscillation is determined by expressing these segmented PSDs in colormaps. Once low-frequency oscillation is detected, identification of the type of sustained oscillation is carried out using properties of segmented PSD and excess kurtosis as described in Table 2. The excess kurtosis value of the selected DMAP is calculated over a window size of one or two cycles depending on the length of the signal. If FO is detected, the source location can be determined by calculating CPSD between the selected DMAP and the raw PMU data from buses associated with generation sources. If the oscillations are originated from one source, the peak of the CPSD for that particular source will be the highest.

Two-area four-machine system
The two-area four-machine model [32] is utilized to generate cases for analysis of FO in the system using the Matlab power system toolbox. The classical generator model in the dq rotor reference frame is used in the system and all the loads have been modulated with 2% of white Gaussian noise to simulate ambient noise. For the generation of FO, a sinusoidal signal with a peak value of 0.1 is used to modulate the shaft torque of G 1 for 2 minutes as shown in Fig. 3.  The raw PMU data at bus 7 is selected for mode decomposition, which has multi-frequency components due to presence of superimposed noise as shown in Fig. 4a. The signal x(t) has 12 modes. However, only three most relevant modes having the highest correlation coefficients with the raw PMU data have been demonstrated as shown in Fig. 4b. Here IMF 6 has the highest correlation coefficient of 0.78, so it is selected as the most relevant DMAP for analysis. The segmented PSD of the selected DMAP removes the frequency components due to ambient noise and shows the existence of the dominant mode characterized by a unique peak at 0.79 Hz. The representation of segmented PSD in colormap shown in Fig. 5a confirms the duration of the sustained oscillation between 20 and 120 s. The variation of excess kurtosis between 0 to − 1.5 between 20 and 120 s as shown in Fig. 5b confirms that sustained oscillation is a FO between the said duration. To locate the source of the oscillation, the voltage magnitudes of all the buses associated with the generation sources are considered and their CPSD with the selected DMAP are estimated. As shown in Fig. 5c, the magnitude of CPSD is highest for bus 1 at 0.79 Hz, which is associated with G 1 . Therefore, it can be concluded that G 1 is the source of the FO.

WECC 179 bus 29 machine system
The proposed scheme is also validated on the test case library specifically created for locating the source of FO in a wide-area system. The library includes simulated phasor data set on the reduced WECC 179-bus, 29machine system. The damping ratio of each generator is selected so that all the natural modes of oscillation, which are of no interest in the present study, are damped out. All the generator sets have a damping D of 4 with all the loads in the system modeled as constant MVA loads. PMU measurements are created at a sampling rate of 30 Hz by the TSAT software from Powertech Labs Inc. for a duration of 40 s, which is sufficient to measure oscillations up to 15 Hz. To reduce the complexity, the portion of the network as shown in Fig. 6 is utilized to test the effectiveness of the proposed scheme based on a summary of FO cases provided in [33].
The raw PMU data at bus 34 is selected for mode decomposition and its PSD shown in Fig. 7a has multifrequency components. The 3 most relevant DMAPs based on the correlation coefficients with the raw PMU data are shown in Fig. 7b. Here IMF 7 is selected as the most relevant DMAP for analysis as it has the highest correlation coefficient of 0.87 with the actual PMU data.
The segmented PSD of the selected DMAP removes the frequency components due to ambient noise and shows the existence of the dominant mode characterized by a unique peak at 0.47 Hz. The representation of segmented PSD in colormap shown in Fig. 8a suggests the existence of a narrow frequency band between 0.4 Hz to 0.6 Hz that exists almost for the entire length of the signal. The variation of excess kurtosis between 0 to − 1.3 as shown in Fig. 8b further confirms that the detected sustained oscillation is a case of FO. The CPSD estimated from all the bus voltages associated with the generator and selected DMAP are shown in Fig. 8c. As can be seen, bus 79 has the highest magnitude and therefore, it can be concluded that the generator associated with bus 79 is the source of FO.

Analysis of oscillation on actual PMU data
The proposed scheme is also tested using real time PMU dataset from ISO New England, which reports five cases of FO. The data is collected at a sampling frequency of 30 Hz and the system is divided into 3 areas. The dataset contains PMU measurements from 12 substations of area 1 while line 7 and line 21 connect the external areas 2 and 3, respectively as shown in Fig. 9 (http://web.eecs.utk.edu/~kaisun/Oscillation/index.html). Two sets of FO data reported on 17 June 2016 and 29 January 2018 are utilized to test the effectiveness of the proposed algorithm.

Forced oscillation case 1
This section analyze the FO event occurred on 17 June 2016 and the raw PMU data at substation 4 is selected for mode decomposition. The signal and its PSD illustrated in Fig. 10a show multi-frequency components due to the presence of noise. The 3 most relevant DMAPs based on the correlation coefficients with the raw PMU data are shown in Fig. 10b. Here IMF 7 is selected as the most relevant DMAP for analysis as it has the highest correlation coefficient of 0.67.  The segmented PSD of the selected DMAP removes the frequency components due to ambient noise and shows the existence of the dominant mode characterized by a unique peak at 0.27 Hz. The representation of segmented PSD in colormap shown in Fig. 11a, suggests the existence of 0.27 Hz components from 30 s to the entire length of the signal. The variation of excess kurtosis between 0 to − 1.5 as shown in Fig. 11b further confirms that the detected sustained oscillation is a case of FO.
CPSD are estimated for PMU measurements from G 1 , G 2 , line 7, line 21, and the selected DMAP as shown in Fig. 11c. The results show the highest magnitude for line 7, which is connected to external area 2. Therefore, it can be concluded that the source of oscillation is from a generator located in the external area 2.

Forced oscillation case 2
This section analyzes the FO event occurred on 29 January 2018. The raw PMU data at substation 4 is selected for mode decomposition and its PSD shown in Fig. 12a has multi-frequency components. The 3 most relevant DMAPs based on the correlation coefficients with the raw PMU data are shown in Fig. 12b. Here IMF 4 is selected as the most relevant DMAP for analysis as it has the highest correlation coefficient of 0.63 with the actual PMU data.
The segmented PSD of the selected DMAP removes the frequency components due to ambient noise and shows the existence of the dominant mode characterized by a unique peak at 1.57 Hz. The representation of segmented PSD in colormap shown in Fig. 13a confirms the duration of the sustained oscillation from 50 s to the entire length of the signal. The variation of excess kurtosis between 0 to − 1.5 from 50 to 240 s as shown in Fig. 13b  confirms that the sustained oscillation is a FO between the said duration. CPSD are estimated for PMU measurements from G 1 , G 2 , line 7, line 21, and the selected DMAP as shown in Fig. 13c. The results show the highest magnitude for generator G 2 and therefore, it can be concluded that the source of oscillation is from generator G 2 located in the internal area 1.

Conclusion
This paper has proposed a systematic two-step approach to locate the source of FO in power system using PMU measurements. It has shown that ICEEMDAN algorithm utilized in the first stage can effectively de-noise the raw PMU data to precisely capture the FO components presented in the system.
The approach is demonstrated using simulated PMU data and further validated using the real-time PMU dataset from ISO New England power system. The first real FO event analyzed has FO of 0.27 Hz from 30 to 180 s, originated from the generation source located outside the considered area. For the second event of 1.57 Hz FO, it is originated from the generator G 2 located in the internal area 1 during 50-240 s. The results obtained on the actual system show the robustness of the proposed system under practical influences.
Compared to other methods, the performance of the proposed method is independent of system dynamics, load modeling assumptions as well as size and complexity of the power system. It can accurately pinpoint the source of FO and provide accurate results even in systems having a low signal to noise ratio.