Synthetical Modal Parameters Identification Method of Damped Oscillation Signals in Power System

Abstract

It is vital to improve the stability of the power system by accurately identifying the modal parameters of damped low-frequency oscillations (DLFO) and controlling the oscillation in time. A new method based on empirical mode decomposition (EMD), stochastic subspace identification (SSI), and Prony algorithms, called synthetical modal parameters identification (SMPI) method, is developed by efficiently matching the modal parameters of DLFO which are acquired from the SSI and Prony algorithm. In this approach, EMD is used for denoising the raw oscillation signals thereby enhancing the noise resistance, and then using the SSI and Prony algorithms to identify the precise modal parameters assisted by parameter matching. It is demonstrated that the proposed SMPI method holds great accuracy in identifying full modal parameters including natural frequencies, damping ratios, amplitudes, and phase angles with simulated signals with known modal parameters and real-time signals from some power system case studies. The strategy of SMPI has effectively overcome the weakness of a single approach, and the identification results are promising to heighten the stabilization of power systems. Besides, SMPI shows the potential to troubleshoot in different fields, such as construction, aeronautics, and marine, for its satisfactory robustness and generalization ability.

1. Introduction

The safety of the electric power supply affects national economy and social development, while the key to the safe operation of the power system lies in the stability of the system [1]. During the operation of the interconnected power system, oscillations frequently occur, where the low-frequency oscillation (LFO) is a detrimental fault threatening the stability of the system, which refers to the relative swing between the rotors of generators running in parallel in power systems under disturbance, and leads to an incessant oscillation within 0.1~2.5 Hz for lack of damping [2,3]. The 1996 system blackout in the Western Electricity Coordinating Council (WECC) system makes people aware of the importance of monitoring and controlling the LFO [4]. The first step to analyze LFO is modal identification, which acquires the modal parameters of the oscillations for taking suitable measures to enhance the stability of the system.

LFO mode identification methodology can be mainly divided into the analysis based on numerical solutions and that based on measured signals. The traditional linearization method is to identify these low-frequency modes through eigenvalue analysis, which requires a lot of information about the modeling of the system, lacking satisfactory precision. With the help of the wide area measurement system (WAMS) and phase measurement unit (PMU), the operation data of the power grid with a high-precision time scale can be transmitted to the control center, realizing the real-time data collection and monitoring of wide-area power grid [5]. The methods of analyzing the measured data and identifying the LFO mode primarily include fast Fourier transform (FFT) [6], wavelet transform [7], Prony [8], stochastic subspace identification (SSI) [9], auto-regressive and moving average (ARMA) [10], Hilbert–Huang Transform (HHT) [11], and so forth. FFT cannot analyze the damping characteristics and local characteristics, so it is not suitable for nonlinear and non-stationary signals. Wavelet transform has the problems of frequency overlap and adaptive basis selection and is only suitable for transient and non-stationary signals. Both Prony and SSI have problems with estimation of order and sensitivity to noise which bring in false modes when analyzing the oscillation signal of the non-stationary power system. The ARMA model method establishes a mathematical model for ordered discrete random data to obtain the inherent parameters of the system, whereas it is easily affected by noise and still has the problem of order determination as Prony and SSI. Although the HHT algorithm is suitable for non-stationary and nonlinear signals, it is affected by the end effect. These methods mentioned above have a certain ability to identify partial modal parameters, but the accuracy of identification is unacceptable yet.

The damped low-frequency oscillation (DLFO) signal is one of the most common transient signal forms, which is often related to mechanical and circuit faults, especially in the power system. To solve the problems of insufficient precision and parameters in modal identification of DLFO, the synthetical modal parameters identification (SMPI) method, a comprehensive method, is presented by integrating the advantages of empirical mode decomposition (EMD), SSI, and Prony. In the SMPI method, the DLFO signals are, firstly, decomposed into several intrinsic mode functions (IMF) to filter high-frequency noise by EMD. Then, the Prony and SSI approaches are adopted to identify the modal parameters of the filtered signals after EMD, respectively. Lastly, the proposed parameters matching method is developed to match the accurate modal parameters of DLFO in line with the similar frequencies recognized by SSI and Prony and avoid the difficulty of estimation of model order. The proposed SMPI method is demonstrated by implementing case studies.

In what follows, Section 2 introduces the theory and methods of SMPI including the EMD, SSI, and Prony algorithm. In Section 3, SMPI is employed to identify the modal parameters of an ideal simulated DLFO signal to verify the effectiveness of the proposed algorithm. To validate the robustness of the method, the simulated signal under lower SNR is identified by SMPI in Section 4. In Section 5, to verify the generalization ability of SMPI, a real-time DLFO signal in the power system is taken as the subject of modal identification. Section 6 summarizes the conclusions of this study.

2. Synthetical Modal Parameters Identification (SMPI) Method

This paper focuses on the theory and methods of SMPI proposed for modal parameter identification of DLFO in power systems. According to previous research, the SSI method is weak in identifying amplitude and phase angle of non-stationary signals [5], and the computation of damping ratio is inaccurate by the Prony approach [12]. Besides sensitivity to noise, estimation of model order is a difficulty for SSI and Prony. Under-estimating model order would lead to the omission of dominant modes of the original signal, while over-estimation of the model order would bring in fictitious modes, especially with the amount of noise.

Accordingly, it is impracticable to utilize a single algorithm to obtain the full modal parameter of DLFO. Addressing this problem, this paper has presented a comprehensive method, i.e., the SMPI method, based on the parameter matching approach which fully makes use of the strengths of SSI approach and Prony method, to find out full modal parameters. In addition, SMPI can handily avoid the difficulty of estimating model orders. The described synthetical modal parameters identification method is shown in Figure 1.

As Figure 1 shows, in the developed SMPI method, the raw signals are initially denoising by EMD. And the filtered signals are put into modal identification with Prony and SSI, respectively. Prony calculates natural frequencies, amplitudes, and phase angles of the signals, while SSI extracts frequencies and damping ratios. In respect of the same natural frequencies, all true modal parameters of DLFO, involving natural frequencies, damping ratios, amplitudes, and phase angles, are matched finally. The detailed implementations of EMD, SSI, and Prony approaches will be presented in the following parts.

2.1. Empirical Mode Decomposition (EMD) Method

In essence, the EMD method [13] is widely used in signal pre-processing, which eliminates high-frequency noise from the raw signal by filtering. The procedure of EMD functioned in denoising the oscillation signal is described as follows:

Step 1: Derivative the input signal x(t) and extract its extreme value (minimum value and maximum value).

Step 2: Fit the upper and lower envelopes of x(t) with cubic spline interpolation function, and get the maximum value e_max(t) and minimum value e_min(t) of the envelope curve.

Step 3: Calculate the mean value m(t) of the upper and lower envelope curves as Equation (1), and then subtract m(t) from x(t) to obtain the residual r(t) as Equation (2).

$$m\left(t\right)=\frac{1}{2}\left(e_{max}\left(t\right)+e_{min}\left(t\right)\right)$$

(1)

$$r\left(t\right)=x\left(t\right)-m\left(t\right)$$

(2)

Step 4: Determine whether r(t) meets IMF conditions: (a) The number of local extrema, that is, the total number of local minima and local maxima, and the number of zero crossings differ by at most one; (b) The mean value of the upper and lower envelopes constructed from the local extrema is zero [13].

Step 5: If the IMF condition is met, take r(t) as one of the IMF components decomposed from x(t), and go to Step 6; If the condition is not met, take r(t) as the new input signal and repeat Step 1 to Step 5 until the condition is met.

Step 6: Determine whether the stop screening condition is met, that is, the residual r(t) is monotonous or a constant. If the condition is met, go to Step 8. If the condition is not met, go to Step 7.

Step 7: Take the difference between x(t) and r(t) as the new input signal, and repeat Step 1 to Step 6.

$$x_{new}\left(t\right)=x\left(t\right)-r\left(t\right)$$

(3)

Step 8: Reshape the original signal into a new signal c_EMD(t) which consisted of total n IMFs c_i(t) and the last residual r_n(t) as Equation (4).

$$c_{EMD}\left(t\right)={\displaystyle \sum _{i=1}^n{c}_i}\left(t\right)+r_n\left(t\right)$$

(4)

2.2. Stochastic Subspace Identification (SSI) Approach

The SSI is an identification method based on discrete state space equations of linear systems, which calculates partial modal parameters [14]. The SSI is split into covariance-driven SSI and data-driven SSI [15]. The covariance-driven SSI uses the covariance as a statistic to describe the data correlation [16], while the data-driven SSI projects the line space of the “future” onto the line space of the “past” to describe the data correlation [17,18]. In this paper, we utilize the covariance-driven SSI for its strength in noise resistance and data processing accuracy [19], which is based on a discrete-time state space model, takes the covariance as the statistic, and uses the singular value decomposition (SVD) to obtain the state matrix and its eigenvalues, and, finally, gains the modal parameters of the system [20]. The subspace involved in this method refers to the state space model of a multi-degree-of-freedom system [19,21]. The state model of SSI is [22]

$$\left\{\begin{array}{l}{x}_{s(}{}_{k+1)}=A{x}_{s(}{}_{k)}+w_k\hfill \\ y_k=C{x}_{s(}{}_{k)}+v_k\hfill \end{array}\right.$$

(5)

where k is the time instant, i.e., the sampling point number of the discrete signal; vector x_s is the state vector of the discrete-time system; vector y_k is the system output vector at time instant k; matrix A is the state matrix of the discrete-time state space equation; matrix C is the system output matrix, also called the observation matrix; and w_k is the process noise and v_k is the measurement noise at time instant k with zero mean.

The measured data y_k, i.e., the system output vector, are constructed into the Hankel matrix, as shown in the following equation [14]:

$$H=Y_{0/2M-1}=\frac{1}{\sqrt{N}}\left(\begin{array}{cccc}{y}_0& y_1& \cdots & y_{N-1}\\ y_1& y_2& \cdots & y_N\\ \cdots & \cdots & \cdots & \cdots \\ y_{M-1}& y_M& \cdots & y_{M+N-2}\\ y_M& y_{M+1}& \cdots & y_{M+N-1}\\ y_{M+1}& y_{M+2}& \cdots & y_{M+N}\\ \cdots & \cdots & \cdots & \cdots \\ y_{2M-1}& y_{2M}& \cdots & y_{2M+N-2}\end{array}\right)=\frac{Y_{0/M-1}}{Y_{M/2M-1}}=\left(\frac{Y_p}{Y_f}\right)$$

(6)

where the number of row blocks is 2M; the number of columns is N which is the number of measured data as well; and Y_p = Y_0/M−1 and Y_f = Y_M/2M−1 denote the past and future parts, respectively, of the block Hankel matrix. Since the number of rows in each row block y_k is l (the number of measured points), and the number of columns is 1, the matrix Y_0/2M−1 consists of 2M × l rows and N columns.

The amount of measured data is limited, i.e., N is not infinite in actual measurement conditions. According to the Hankel matrix obtained above, the covariance of the output vector can be constructed as [19]

$$R_M=\frac{1}{N}{\displaystyle \sum _{k=0}^{N-1}{y}_{k+M}}{y}_k^T$$

(7)

The Hankel matrix in Equation (6) can be converted to the Toeplitz matrix as [23].

$$T_{1/M}=Y_f{Y}_p^T=\left(\begin{array}{cccc}{R}_M& R_{M-1}& \cdots & R_1\\ R_{M+1}& R_M& \cdots & R_2\\ \cdots & \cdots & \cdots & \cdots \\ R_{2M-1}& R_{2M-2}& \cdots & R_M\end{array}\right)$$

(8)

To obtain the state matrix A, the observable matrix O_M should be determined. Singular value decomposition (SVD) is used to perform the above factorization, i.e.,

$$T_{1/M}=US{V}^T=\left(\begin{array}{ll}{U}_1\hfill & U_2\hfill \end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& S_2=0\end{array}\right)\left(\begin{array}{c}{V}_1^T\\ V_2^T\end{array}\right)=U_1{S}_1{V}_1^T$$

(9)

where matrix U₁ and V₁ are orthogonal matrices; matrix S₁ is a diagonal matrix composed of positive singular matrices. The number of singular values of the S₁ matrix is 2m, which is the rank of the S₁ matrix; m is the order of the system.

T_1/M can be divided as

$$T_{1/M}={\mathit{O}}_M{\mathsf{\Gamma }}_M$$

(10)

where O_M is the observable matrix and Γ_M is the extended observability matrix.

The matrices O_M and Γ_M are expressed by

$$\left\{\begin{array}{l}{O}_M=U_1{S}_1^{1/2}\\ {\mathsf{\Gamma }}_M=I^{-1}{S}_1^{1/2}{V}_1^T\end{array}\right.$$

(11)

where I is the unit matrix.

The output matrix C of the system can be obtained by the front l row of matrix O_M.

$$C=O_M(1:l)$$

(12)

The state matrix A can be calculated as follows:

$$A=S_1^{-1/2}{U}_1^T{T}_{2/M+1}{V}_1{S}_1^{-1/2}$$

(13)

The state matrix A and output matrix C are both identified from all measured signals rather than one single signal. Hence, the modal parameters can be attained through the eigen decomposition of the state matrix A as follows:

$$A=\Psi \Lambda {\Psi }^{-1}$$

(14)

where Ψ is the complex eigenvector matrix and Λ = diag(μ₁, μ₂, _…, μ_m) is a diagonal matrix composed of the eigenvalues of the system μ_m.

The eigenvalues of the system λ_m can be calculated from the eigenvalues of the state matrix A:

$${\lambda }_m=\frac{\ln {\mu }_m}{\Delta t}$$

(15)

where Δt is the time interval.

According to the above principle, the modal frequency f_s and damping ratio ξ_s corresponding to each eigenvalue can be calculated by

$$\left\{\begin{array}{l}{f}_s=\frac{\sqrt{\left|{\lambda }_m\right|}}{2\pi }\\ {\xi }_s=-\frac{{\lambda }_m+{\lambda }_m^{\ast }}{2\sqrt{{\lambda }_m{\lambda }_m^{\ast }}}\end{array}\right.$$

(16)

2.3. Prony Algorithm

The Prony algorithm is a mathematical model in which equally spaced sampled data is represented by a linear combination of complex exponential functions [24]. The LFO mode can also be represented by a complex exponential function, so Prony can be applied to LFO. The Prony algorithm is fast and reliable, which does not depend on the mathematical model of the system [25,26].

Given a complex-valued data sequence x(n), n = 0, 1, …, N − 1, based on N original data, the traditional Prony method fits an exponential model to the data in the least-squares sense, and estimates $\widehat{x}(n)$ [27,28,29], as

$$\widehat{x}(n)={\displaystyle \sum _{i=1}^p}{A}_i{e}^{j{\theta }_i}{e}^{\left({\alpha }_i+j2\pi f_i\right)\Delta t}={\displaystyle \sum _{i=1}^p}{b}_i{z}_i^n n=0,1,\cdots ,N-1$$

(17)

where p is equal to the order of system; j is the imaginary unit; the subscript i denotes the i-th mode; A_i, f_i, θ_i and α_i signify amplitude, frequency, initial phase, and damping factor, respectively; b_i and z_i are the i-th residue and polynomial root, correspondingly; and Δt is the time interval.

Considering the linear prediction model is

$$\left[\begin{array}{c}{x}_{(p)}\\ x_{(p+1)}\\ ⋮\\ x_{(N-1)}\end{array}\right]=\left[\begin{array}{cccc}{x}_{(p-1)}& x_{(p-2)}& \cdots & x_{(0)}\\ x_{(p)}& x_{(p-1)}& \cdots & x_{(1)}\\ ⋮& ⋮& \ddots & ⋮\\ x_{(N-2)}& x_{(N-3)}& \cdots & x_{(N-p-1)}\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_p\end{array}\right]$$

(18)

which can be denoted as d = Da, where a is the vector of coefficients, and p is the polynomial order assumed to be known a priori.

Estimation of a is obtained by solving Equation (18) in the least-squares sense as

$$\widehat{a}={\left(D^{\mathrm{H}}D\right)}^{-1}{D}^{\mathrm{H}}d$$

(19)

where D^H is the complex-conjugate transpose of D.

Using the previously estimated coefficients, the roots of the characteristic polynomial expressed as

$$z^p-\left(a_1{z}^{p-1}+a_2{z}^{p-2}+\dots +a_p{z}^0\right)=0$$

(20)

where each root zⁱ, i = 1, …, p, is a discrete-time approximation of its respective continuous-time eigenvalue in the Z-domain.

The i-th mode frequency is calculated as

$$f_{p(i)}=\arctan \left[\mathrm{Im}\left(z_i\right)/\mathrm{Re}\left(z_i\right)\right]/(2\pi \Delta t) i=1,\dots ,p$$

(21)

Using the polynomial roots, z_i, the linear regression model given by

$$\left[\begin{array}{c}{x}_{(0)}\\ x_{(1)}\\ ⋮\\ x_{(N-1)}\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ z_1^1& z_2^1& \cdots & z_p^1\\ ⋮& ⋮& \ddots & ⋮\\ z_1^{N-1}& z_2^{N-1}& \cdots & z_p^{N-1}\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_p\end{array}\right]$$

(22)

which can be denoted as x = Vb, where V is a Vandermonde matrix.

Using the least square method, vector b is estimated as

$$\widehat{b}={\left(V^{\mathrm{H}}V\right)}^{-1}{V}^{\mathrm{H}}x$$

(23)

The i-th mode frequency, amplitude, and initial phase are, respectively, calculated as

$$\left\{\begin{array}{l}{A}_{p(i)}=\left|b_i\right|\\ {\theta }_{p(i)}=\arctan \left[\mathrm{Im}\left(b_i\right)/\mathrm{Re}\left(b_i\right)\right]\\ f_{p(i)}=\arctan \left[\mathrm{Im}\left(z_i\right)/\mathrm{Re}\left(z_i\right)\right]/(2\pi \Delta t)\end{array}\right.i=1,\cdots ,p$$

(24)

2.4. Parameter Matching Approach

The Prony algorithm only extracts precisely natural frequencies, amplitudes, and phase angles of DLFO, while SSI accurately obtains natural frequencies and damping ratios. Due to noise interference and overestimation of order, algorithms often identify modes containing false information, i.e., false modes (or fictitious modes). In line with the features of SSI and Prony, it is needed to execute parameter matching to search for the exact amplitudes, phase angles, and damping ratios according to the similar frequencies from SSI and Prony. In respect of the previous investigation, there is no such method that integrates SSI with Prony to identify the full modal parameters of DLFO accurately.

The steps for this technique are as follows:

Step 1: Calculate the dimensions d_s and d_p of the frequency f_s acquired by SSI and the frequency f_p computed by Prony, respectively;

Step 2: Determine whether the number of iterations si is less than or equal to d_s. If not met the condition, the modal parameter result, i.e., the successful matched frequencies f_mi, damped ratios ξ_mi, phase angles θ_mi and A_mi will be output; if the condition is met, jump to Step 3;

Step 3: Calculate the absolute error Δf, the minimum absolute error Δf_min, and the relative error δ_f between f_s and f_p;

Step 4: Determine whether the number of iterations pi is less than or equal to d_p. If the condition is not met, go to Step 2, carrying out the next iteration; if the condition is met, determine whether the parameter matching conditions are fulfilled in Step 5;

Step 5: Determine whether f_p(i) meets the matching condition, that is, the number of successful matching js is less than d_s, and the relative error δ_f(pi) is less than the allowable error value ε_a, and the absolute error Δf_pi is just the minimum value Δf_min. If the condition is not met, go to Step 4 to carry out the next iteration; if the condition is met, go to Step 6 to carry out the parameter matching; and

Step 6: Matching parameters, that is, f_mi = f_s(si), ξ_mi = ξ_s(si), θ_mi = θ_p(pi), A_{mi =} A_p(pi). After completing a round of pairing, return to Step 4 to carry out the next iteration.

The flowchart of parameter matching of DLFO through the SMPI method describes the detailed process of the steps above as Figure 2.

3. Validation of SMPI with Simulated Signals

To verify the performance of the method in modal identification, a simulated signal with known modal parameters is identified with SMPI in this section.

3.1. Test of the Simulated Signal

The ideal simulated signal close to the DLFO of the power system is constructed for verifying the effectiveness of SMPI, as

$$\begin{array}{c}y=2e^{-0.12t}\cdot \cos (2\pi \cdot 0.33t-1.5\pi )\\ +4e^{-0.13t}\cdot \cos (2\pi \cdot 0.78t+0.5\pi )\\ +6e^{-0.71t}\cdot \cos (2\pi \cdot 1.00t+1.5\pi )\end{array}$$

(25)

which is a synthetic signal consisting of three dominant modes, and the three oscillation frequency ranges are all within the frequency range (0.1–2.5 Hz) of the LFO of the power system introduced previously. And its true modal parameters are shown in

Table 1. Modal parameters of the ideal simulated signal.

True Modes	Frequency (Hz)	Amplitude (μm)	Damping Ratio	Phase Angle (Rad)
1	0.33	2	0.0579	1.5708
2	0.78	4	0.0265	1.5708
3	1.00	6	0.1130	−1.5708

.

To further simulate the real oscillation signal, Gaussian white noise at 30 dB SNR is added into the simulated signal as Figure 3, and it will be referred to as signal A.

3.2. Oscillation Signal Denoising

In this subsection, the EMD method is utilized to filter and smooth the DLFO signal with noise. The simulated signal with Gaussian white noise is decomposed into many IMFs with different mode orders and a residual r (i.e., the last IMF in Figure 4) by EMD. After eliminating higher-order noise and trend term components, five effective IMFs and a residual signal are exacted as Figure 5.

After denoising, the new signal c_EMD(t) is reconstructed by IMFs and r. The comparison of the synthesized signal by EMD c_EMD(t) with the original signal with noise (i.e., signal A) is displayed in Figure 6.

The curve of c_EMD(t) is smooth and extremely coincidence with that of signal A, which indicates that EMD has filtered high-frequency noise and retained dominant modes of the original signals. c_EMD(t) is conducive to modal identification by SSI and Prony in the next subsection.

3.3. Modal Parameters Identification with SSI and Prony

In this section, the modal parameters of DLFO signal A are identified in respect of SSI and Prony approaches. Firstly, SSI is adopted to calculate the frequencies and damping ratios of c_EMD(t) and the results are displayed in

Table 2. Modal parameters of signal A identified with SSI approach.

Modes	Frequency (Hz)	Damping Ratio
1	0.3318	0.0589
2	0.7807	0.0258
3	1.0072	0.1119
4	2.0134	0.0104

.

As shown in Table 2, there are four model orders identified by SSI, including one false mode (the last one). And only frequencies and damping ratios have been recognized from each model order, lacking other modal parameters. Subsequently, it is necessary to identify more modal parameters of c_EMD(t), i.e., amplitudes and phase angles in another way. After Prony analysis, the frequencies, amplitudes, and phase angles of the first 99 modes are attained in

Table 3. Modal parameters of signal A identified with Prony approach.

Modes	Frequency (Hz)	Amplitude (μm)	Phase Angle (Rad)
1	0.3305	2.0785	1.5942
2	0.7808	3.9976	1.6011
3	1.0069	5.7695	−1.5118
4	1.9077	0.0179	−0.1935
5	2.4386	0.0354	−0.6599
⁝	⁝	⁝	⁝
99	49.4514	0.0014	−2.7861

.

Although Prony has identified other more modal parameters, amplitude, and phase angle in addition to frequency, compared with SSI analysis, there were still a large number of false modes.

To avoid missing dominant modes, the number of model orders is set relatively higher in this study, so the false modes are generated owing to the over-estimation of model order, which is difficult to determine in SSI and Prony. Therefore, it is needed to extract the true modes from the trivial modes through a novel method.

With respect of Table 2 and Table 3, SSI has recognized 4 modes while 99 modes have been found by Prony. Although false modes exist in Table 2 and Table 3, the true modes of DLFO exist as well. Besides, identified frequencies are both displayed in these two tables except for damping ratios, amplitudes, and phase angles. Therefore, the key to seek for each true model order is to match two sets of the results of identified modal parameters in Table 2 and Table 3 according to the principle of the similar frequency. This work will be completed in the next subsection.

3.4. Parameter Matching

According to the modal parameter identification results obtained by the SSI and Prony algorithm, they are paired with each other according to the principle of the similar frequency, whose detailed procedure has been illustrated in Section 2.4.

If the relative error between f_s and f_p is less than 0.005, it is accepted that f_s is one of the natural frequencies of DLFO, and then the corresponding modal parameters can be found with respect to the natural frequencies. Using parameter matching, there are three sets of similar frequencies have been matched successfully, that is, 0.3318 Hz, 0.7807 Hz, and 1.0072 Hz which have been determined as natural frequencies of DLFO. And other corresponding modal parameters are the modal parameters of true modes in these three groups, respectively, summarized in

Table 4. Modal parameters matching results of signal A.

True Modes	Frequency (Hz)	Amplitude (μm)	Damping Ratio	Phase Angle (Rad)
1	0.3318	2.0785	0.0589	1.5942
2	0.7807	3.9976	0.0258	1.6011
3	1.0072	5.7695	0.1119	−1.5118

.

3.5. Validation of Identified Results

Substituting the matching results in Table 4 into the classic decaying oscillation signal mathematical equation [30,31], the governing equation of three components in DLFO is

$$X_i=A_i{e}^{(-2\pi f_i{\xi }_i\Delta t)}\cos (\sqrt{1-{\xi }_i^2}\cdot 2\pi f_i\Delta t+{\theta }_i)i=1, 2, 3$$

(26)

In which the subscript i indicates the i-th dominant mode and X_i is the i-th dominant component.

In light of Equation (26), the three components of DLFO signal A are reconstructed as Figure 7.

The three signal components are superimposed to obtain a synthesized signal after SMPI as Equation (27).

$$Y={\displaystyle \sum X_{mi}}$$

(27)

where Y is the synthesized signal with SMPI.

Subsequently, the comparison between Y and c_EMD(t) is shown in Figure 8. And the comparison of estimated modal parameters of SMPI with true modal parameters is described in

Table 5. Comparison of estimated modal parameters of SMPI with true modal parameters.

Modal Parameters	True Modes	True Value	Estimated Value	Relative Error
Frequency (Hz)	1	0.33	0.3318	0.0055
	2	0.78	0.7807	0.0009
	3	1.00	1.0072	0.0072
Amplitude (μm)	1	2.00	2.0785	0.0392
	2	4.00	3.9976	0.0006
	3	6.00	5.7695	0.0384
Damping ratio	1	0.0579	0.0589	0.0167
	2	0.0265	0.0258	0.0252
	3	0.1130	0.1119	0.0098
Phase angle (rad)	1	1.5708	1.5942	0.0149
	2	1.5708	1.6011	0.0193
	3	−1.5708	−1.5118	0.0375

.

In light of Table 5, the value of estimated modal parameters is very close to that of true modal parameters, and the relative errors between them are comparatively small. The frequencies calculated by SMPI are all within 1% relative error. And the relative error of damping ratios and phase angles are less than 3% and 4%, respectively. As for the relative error of amplitude, the smallest is less than 0.1%, and the biggest one is not more than 4%.

As Figure 8, the curve of signal Y synthesized with the identified results of SMPI is found in good coincidence with that of c_EMD(t), whose dominant modal parameters have been extracted successfully.

4. Robustness Test of SMPI under Lower SNR

To verify the robustness of SMPI, the Gaussian white noise in the ideal simulated signal increased to 20 dB SNR. The new signal will be referred to as signal B following and its curve is illustrated in Figure 9.

The five effective IMFs and a residual signal are extracted by EMD to eliminate high-frequency noise as Figure 10, and the comparison of the synthesized signal c_EMD(t) with signal B in Figure 11.

The curve of c_EMD(t) is smooth and extremely coincidence with that of signal B, which indicates that c_EMD(t) has retained the dominant modes of the original DLFO signal. And c_EMD(t) is conducive to modal identification by SSI and Prony following.

Using the SSI method and Prony algorithm to identify c_EMD(t), the estimated modal parameters are shown in

Table 6. Modal parameters of signal B identified with SSI approach.

Modes	Frequency (Hz)	Damping Ratio
1	0.3306	0.0607
2	0.7813	0.0270
3	1.0083	0.1124
4	2.1442	0.0214

and

Table 7. Modal parameters of signal B identified with Prony approach of signal B.

Modes	Frequency (Hz)	Amplitude (μm)	Phase Angle (Rad)
1	0.3275	2.0790	1.7038
2	0.7804	4.3000	1.6164
3	1.0129	6.1616	−1.5975
4	1.6755	0.0614	−2.9024
5	2.3345	0.0218	2.8943
⁝	⁝	⁝	⁝
99	49.4658	0.0018	−2.7895

, respectively.

To extract the true modes from the trivial modes, the parameter matching method is applied to match the full modal parameters in the similar frequency of Table 6 and Table 7, and the matching results are displayed in

Table 8. Modal parameters matching results of signal B.

True Modes	Frequency (Hz)	Amplitude (μm)	Damping Ratio	Phase Angle (Rad)
1	0.3306	2.0790	0.0607	1.7038
2	0.7813	4.3000	0.0270	1.6164
3	1.0083	6.1616	0.1124	−1.5975

.

As shown in Table 8, there are three dominant components. In light of Equation (26), the three components of signal B are reconstructed as Figure 12.

The three signal components are superimposed to obtain a synthesized signal after SMPI as Equation (27). Subsequently, the comparison between Y and c_EMD(t) is shown in Figure 13. And the comparison of estimated modal parameters of SMPI with true modal parameters is described in

Table 9. Comparison of estimated modal parameters of SMPI with true modal parameters.

Modal Parameters	True Mode	True Value	Estimated Value	Relative Error
Frequency (Hz)	1	0.33	0.3306	0.0017
	2	0.78	0.7813	0.0017
	3	1.00	1.0083	0.0083
Amplitude (μm)	1	2.00	2.0790	0.0395
	2	4.00	4.3000	0.750
	3	6.00	6.1616	0.269
Damping ratio	1	0.0579	0.0607	0.0492
	2	0.0265	0.0270	0.0192
	3	0.1130	0.1124	0.0050
Phase angle (rad)	1	1.5708	1.7038	0.0846
	2	1.5708	1.6164	0.0290
	3	−1.5708	−1.5975	0.0170

.

In light of Table 9, the value of estimated modal parameters is close to that of true modal parameters, and the relative errors between them are comparatively small as well. The frequencies calculated by SMPI are still all within 1% relative error. While the relative errors of amplitude, damping ratio, and phase angle are less than 8%, 5%, and 8.5%, respectively.

Although the relative errors have increased in a sense, due to the introduction of more noise, the curve fitted by the identified results of SMPI is consistent with the curve of the original simulated signal whose main oscillating features are obtained without omission in Figure 13. Hence, the proposed SMPI method has shown satisfactory robustness under the condition of noise with different SNRs.

5. Generalization Ability Test of SMPI with Real-Time Signals

To test the effectiveness in modal identification of real-time DLFO signals, SMPI is adopted to identify a set of sampling signals from a power system, whose sampling frequency is 100 Hz, the sampling time is 20 s and sampling points is 2000, as shown in Figure 14.

Totally nine IMFs and a residual have been decomposed from the real-time signal in the power system polluted by noise in the process of the EMD method, as shown in Figure 15. After EMD denoising, the four effective IMFs and a residual signal are extracted as Figure 16.

The comparison of the synthesized signal by EMD with the real-time signals is in Figure 17. The curve of c_EMD(t) fits the oscillating features of the polluted real-time signals smoothly, which illustrates that the EMD method is workable in the signal pre-process.

Afterward, the SSI and Prony algorithms are utilized to purchase the modal parameters of the synthesized signal after EMD denoising, and their identification results are presented in

Table 10. Modal parameters of real-time signal identified with SSI approach.

Modes	Frequency (Hz)	Damping Ratio
1	0.5390	0.0089
2	0.9811	0.0891
3	14.7210	0.0583
4	26.8585	0.0607
5	35.5469	0.0573
6	38.2654	0.0120
7	45.5538	0.0124

and

Table 11. Modal parameters of real-time signal identified with Prony approach.

Modes	Frequency (Hz)	Amplitude (μm)	Phase Angle (Rad)
1	0.5397	0.0579	2.7909
2	0.9778	0.3053	1.8003
3	1.6909	0.3829	−3.1258
4	2.2648	0.0343	1.7157
5	2.6177	0.0217	2.0594
⁝	⁝	⁝	⁝
98	49.6196	0.0339	1.7125

.

There are many high-frequency fictitious modes so the parameter matching method is carried out to filter out the true modes, and the matching results are displayed in

Table 12. Modal parameters matching results of real-time signal.

True Modes	Frequency (Hz)	Amplitude (μm)	Damping Ratio	Phase Angle (Rad)
1	0.5390	0.0579	0.0089	2.7909
2	0.9811	0.3053	0.0891	1.8003

. There are two dominant components of the real-time signal, and they can be reconstructed as shown in Figure 18.

These two dominant components of real-time signal are superimposed to obtain a synthesized signal Y after SMPI as Equation (27). Subsequently, the comparison between Y and c_EMD(t) is shown in Figure 19.

Given that the signal sampling is unstable and the algorithm needs a certain convergence process at the beginning, there are a few errors Y and c_EMD(t) of the real-time signal. After 0.6 s, the relative errors are controlled well and gradually decrease, owing to the ability of SMPI to fast convergence and the stability of signal sampling. In Figure 19, the curve of signal Y synthesized with the identified results of SMPI is found in good coincidence with that of c_EMD(t), whose dominant modal parameters have been extracted successfully.

In other words, the developed SMPI is a practicable tool, which holds the generalization ability to identify the full modal parameters including frequency, damping ratio, amplitude, and phase angle of DLFO real-time signal in the power system.

6. Conclusions

To solve the problems of insufficient precision and parameters in modal identification of damped low-frequency oscillation (DLFO), a comprehensive method, called synthetical modal parameters identification (SMPI) method, is presented by integrating the advantages of empirical mode decomposition (EMD), stochastic subspace identification (SSI) and Prony algorithm assisted by parameter matching. In the SMPI method, the DLFO signals are, firstly, denoised by EMD. Afterward, the Prony and SSI approaches are adopted to identify the modal parameters of the filtered signals after EMD, respectively. Lastly, the proposed parameters matching method is developed to match the accurate modal parameters of DLFO in line with the similar frequencies recognized by SSI and Prony, and avoid the difficulty of estimation of model order. The proposed SMPI method is demonstrated by ideal simulated signals with known modal parameters and real-time signals from power system case studies. The primary conclusions are as follows:

(1)

To solve the problem that SSI and Prony are both sensitive to noise, EMD shows the potential to denoise the original DLFO signals and cut down the occurrence of fictitious modes to some extent, enhancing the accuracy of modal identification;

(2)

Integrating the strengths of Prony and SSI by the parameter matching method is capable of purchasing precise and full modal parameters of DLFO, and handily avoiding the difficulty of estimating model orders in Prony and SSI; and

(3)

Through the case studies on simulated signals with known modal parameters and real-time signals from some power systems, it is demonstrated that SMPI holds satisfactory precision, robustness, and generalization ability.

This study on the proposed SMPI method is conducive to improve the stability of the power system by controlling and even avoiding the oscillation in time. Furthermore, SMPI shows the potential for prognostics and health management in different conditions and fields, such as construction, aeronautics and marine for its robustness and generalization ability.