Application of Synchrosqueezed Wavelet Transforms for Extraction of the Oscillatory Parameters of Subsynchronous Oscillation in Power Systems

The most classical subsynchronous oscillation (SSO) mode extraction methods have some shortcomings, such as lower mode identification and poor anti-noise properties. Thus, this paper proposes a new time-frequency analysis method, namely, synchrosqueezed wavelet transforms (SWT). SWT combines the advantages of empirical mode decomposition (EMD) and wavelet, which has the adaptability of EMD, and improve the ability of anti-mode mixing on EMD and wavelet. Thus, better anti-noise property and higher mode identification can be achieved. Firstly, the SSO signal is transformed by SWT and its time-frequency spectrum is obtained. Secondly, the attenuation characteristic of each intrinsic mode type (IMT) component in its time-frequency spectrum is analyzed by an automatic identification algorithm, and determine which IMT component needs reconstruction. After that, the selected IMT components with divergent characteristic are reconstructed. Thirdly, high-accuracy detection for mode parameter identification is achieved by the Hilbert transform (HT). Simulation and application examples prove the effectiveness of the proposed method.


Introduction
Subsynchronous oscillation (SSO) is a classic problem in power systems.A steam turbine and electric generator can have mechanical SSO modes between the various turbine stages and the generator.Now, with the increasing penetration of renewable energy sources, especially the wind power generation, and broad application of the high-power electronic technology, electric power system with renewable energy sources has evolved, which is constituted by complicated AC-DC systems with multiple energy resources and multiple conversions.As a result, new SSO problems are constantly prominent [1,2].SSO arising from grid-connection of the new energy first occurred in a wind farm in Texas in 2009 [3,4], and similar SSO have been reported in the Buffalo Ridge area of Canada [5] and in some parts of north China [6].It is a serious threat to the security and stability of power system operation.
The high electromechanical coupling characteristics of the power system and the strong non-linearity usually spawn multiple SSO modes simultaneously.Therefore, the SSO signal is a typical non-linear and non-stationary signal.The frequency of the SSO is generally between 5 Hz and 45 Hz, and each mode frequency is often close to each other.Moreover, the SSO signals measured Energies 2018, 11, 1525 2 of 18 through a phasor measurement unit (PMU) in actual electrical networks often contain some noise.This requires that the mode extraction method have a high mode resolution and some anti-noise ability.
In recent years, time-frequency analysis has become a hot field in signal processing.It is good at analyzing non-stationary signals.Synchrosqueezed wavelet transform (SWT) is a new time-frequency analysis method proposed by Daubechies et al. in 2011 [7].She called it: an empirical mode decomposition (EMD)-like tool.The anti-noise ability and time-frequency resolution of SWT are improved on the basis of wavelet transform (WT).The calculation principle of SWT is not complicated.After the WT, the corresponding frequency of each scale is given by the derivative of the wavelet coefficients with respect to time.Add the scale of the same frequency to obtain the higher resolution time-frequency curve.The above steps greatly improve the energy divergence of the WT and time-frequency resolution.In other words, it also improves the readability of the time-frequency [8].SWT retains the advantages of both EMD and WT.SWT is adaptive like EMD and does not depend on the mother wavelet; the mode mixing problem is improved greatly; moreover, SWT also has good anti-noise ability.At present, SWT has been applied to geophysics [9], civil engineering [10], low-frequency oscillation analysis [11] and other fields, and achieved good results.
In this paper, SWT is used to analyze the SSO signal firstly, and its time-frequency spectrum is obtained.Secondly, the attenuation characteristic of each intrinsic mode type (IMT) component in its time-frequency spectrum is analyzed, to determine which IMT component needs attention and reconstruction.After that, the selected IMT components with divergence characteristic are reconstructed.Thirdly, high-accuracy detection for mode parameter identification is achieved by the Hilbert transform (HT).

Forward Transformation
Given a signal s(t), the wavelet coefficients W s (a, b) are obtained by continuous wavelet transform (CWT): where variable a is the scale factor, it controls the scaling of wavelet and is inversely proportional to frequency.b is the shift factor, it controls the translation of the wavelet and associated with time.We can understand the relationship between scale factor and actual frequency in this way.Small scale factor represents a wavelet in the time compression, which measured details of the signal, and it also means that the corresponding analysis frequency is high.By contrast, large scale factor represents a wavelet in the time extension, which measured approximations of the signal, and it means a low analysis frequency.ψ is an appropriately chosen wavelet , and ψ( t−b a ) is the complex conjugate of ψ( t−b a ).The observation made in [12] is that, despite W s (a, b) disperses in a, no matter how a changes its oscillation in b refers to the original frequency ω.Taking partial derivative of instantaneous frequency W S (a, b) with respect to b and get ω s (a, b): In Equation (2), the frequency variable ω and the scale variable a are "binned", and denote In actual calculation, a, b and ω are discretized, i.e., W s (a,b) is computed only at a discrete a k , with a k − a k−1 = (∆a) k , and its synchrosqueezed transform T s (ω, b) is determined only at the

Inverse Transformation
Reconstruction can be obtained by inverse transformation of T s : where is the complex conjugate of the Fourier transform of ψ(ξ).Finally, extraction of the real part is the reconstruction.

The Procedure of Parameter Extraction and Calculation Based on Synchrosqueezed Wavelet Transforms
In summary, the method procedure is divided into the following steps: (1) Transform the SSO signal by SWT using Equation ( 3) and obtain it's a time-frequency spectrum.
(2) Search for the IMT components which need to be concerned in the time-frequency spectrum which is the divergent SSO mode.
Image observation of time-frequency picture is usual practical means to search for the divergent SSO mode components.In order to realize a quantitative indicator that can be automatically identified by computer and achieve automatic selection, this paper puts forward an automatic identification algorithm to solve it.
Suppose there is a signal with multifrequency modes.After obtaining the time-frequency spectrum by SWT, the frequency of the mode that we want to analyze is the target frequency f a , and the frequency of its adjacent mode is set to f b .
Suppose the upper limit is and n is the number of mode components in the spectrum.If the signal only contains a single frequency mode, f aH and f aL can be obtained by exceeding the upper and lower limits of the frequency in the spectrum, but its value should not exceed the limit by very much.
In the time-frequency spectrum, assume A j is the energy of the rectangular window.The length of the window is the number of sampling point K, and the width is the frequency range [f aL , f aH ].The equation is: where p is the location of the first sample point, and j is an integer, j = 1, 2, . . .Move rectangular window, and continuously calculate three to four values of A j .If A j keeps getting bigger, it shows that this mode is not attenuated with time.We need to pay close attention to this mode and analyze it further.
(3) Reconstruct the divergence mode components by inverse SWT.Identify the parameters of the reconstructed mode components by HT.With the HT method, the mode component's oscillation frequency and attenuation coefficients are obtained.Readers may refer to [13] for the details of the HT method.
Figure 1 shows the flow chart of the proposed method.Figure 1 shows the flow chart of the proposed method.

Simulation Signal Analysis
For verifying the performance of SWT for the SSO signals, in this section, a numerical simulation and IEEE first benchmark model are used to calculate and analyze.Suppose there is a signal with multifrequency modes, and the simulated signal can be approximated by: where Ai is the amplitude, fi is the oscillation frequency, σi is damping factor, Φi is the phase shift, λ(t) is white noise.A signal simulated with Equation ( 6) is described in Table 1.In addition, white noise is added to the signal, and the signal-to-noise ratio (SNR) is up to 20 dB.The simulation lasts 10 s, the result is shown in Figure 2. SWT does not depend on the mother wavelet selection [14,15], so this paper chooses Morlet wavelet as mother wavelet.Morlet wavelet is a non-orthogonal wavelet-based function which doesn't have a scaling function.It is composed of a complex exponential (carrier) multiplied by a Gaussian window: where

Simulation Signal Analysis
For verifying the performance of SWT for the SSO signals, in this section, a numerical simulation and IEEE first benchmark model are used to calculate and analyze.Suppose there is a signal with multifrequency modes, and the simulated signal can be approximated by: where A i is the amplitude, f i is the oscillation frequency, σ i is damping factor, Φ i is the phase shift, λ(t) is white noise.A signal simulated with Equation ( 6) is described in Table 1.In addition, white noise is added to the signal, and the signal-to-noise ratio (SNR) is up to 20 dB.The simulation lasts 10 s, the result is shown in Figure 2. SWT does not depend on the mother wavelet selection [14,15], so this paper chooses Morlet wavelet as mother wavelet.Morlet wavelet is a non-orthogonal wavelet-based function which doesn't have a scaling function.It is composed of a complex exponential (carrier) multiplied by a Gaussian window: where . Its Fourier transform is shown in Equation (8).
Energies 2018, 11, 1525 where t is time and ω corresponds to frequency.The time-frequency resolution is regulated by changing parameter σ.Through a large amount of simulated examples, it indicates that, Morlet wavelet has good performance of resolution both in time and frequency domain for the SSO signals when σ = 12.
Energies 2018, 11, x 5 of 16  The simulated signal of Figure 2 is transformed by SWT, and the result is shown in Figure 3.In Figure 3, it is clear and intuitive to observe that the energy distribution of each oscillation mode varies with time and frequency.The mode frequencies are 15 Hz, 20 Hz, 25 Hz and 32 Hz, respectively.Among them, the energy of the two oscillation modes, 15 Hz and 20 Hz, decreases gradually with time and they are positive damping modes, while the energy of the other two oscillation modes, 25 Hz and 32 Hz, enhances gradually and they are negative damping modes.In particular, 32 Hz oscillation mode with negative damping property is much stronger, which mode needs to be paid attention to.However, image observation method requires the participation of the people, so the method is relatively inefficient.The proposed automatic recognition algorithm in Section 3 can solve this problem.Take the mode frequency of 15 Hz as an example, the frequency of its adjacent mode is 20 Hz, and the spectrum contains four mode components, so the target frequency fa = 15, the frequency of its adjacent mode fb = 20, and n = 4: Each decomposition amount could generally appear Gibbs phenomenon at its boundary while WT.To overcome negative impact of the end effect, we select the location of the first sample point p The simulated signal of Figure 2 is transformed by SWT, and the result is shown in Figure 3.In Figure 3, it is clear and intuitive to observe that the energy distribution of each oscillation mode varies with time and frequency.The mode frequencies are 15 Hz, 20 Hz, 25 Hz and 32 Hz, respectively.Among them, the energy of the two oscillation modes, 15 Hz and 20 Hz, decreases gradually with time and they are positive damping modes, while the energy of the other two oscillation modes, 25 Hz and 32 Hz, enhances gradually and they are negative damping modes.In particular, 32 Hz oscillation mode with negative damping property is much stronger, which mode needs to be paid attention to.
Energies 2018, 11, x 5 of 16 where t is time and ω corresponds to frequency.The time-frequency resolution is regulated by changing parameter σ.Through a large amount of simulated examples, it indicates that, Morlet wavelet has good performance of resolution both in time and frequency domain for the SSO signals when σ = 12.The simulated signal of Figure 2 is transformed by SWT, and the result is shown in Figure 3.In Figure 3, it is clear and intuitive to observe that the energy distribution of each oscillation mode varies with time and frequency.The mode frequencies are 15 Hz, 20 Hz, 25 Hz and 32 Hz, respectively.Among them, the energy of the two oscillation modes, 15 Hz and 20 Hz, decreases gradually with time and they are positive damping modes, while the energy of the other two oscillation modes, 25 Hz and 32 Hz, enhances gradually and they are negative damping modes.In particular, 32 Hz oscillation mode with negative damping property is much stronger, which mode needs to be paid attention to.However, image observation method requires the participation of the people, so the method is relatively inefficient.The proposed automatic recognition algorithm in Section 3 can solve this problem.Take the mode frequency of 15 Hz as an example, the frequency of its adjacent mode is 20 Hz, and the spectrum contains four mode components, so the target frequency fa = 15, the frequency of its adjacent mode fb = 20, and n = 4: Each decomposition amount could generally appear Gibbs phenomenon at its boundary while WT.To overcome negative impact of the end effect, we select the location of the first sample point p However, image observation method requires the participation of the people, so the method is relatively inefficient.The proposed automatic recognition algorithm in Section 3 can solve this problem.Take the mode frequency of 15 Hz as an example, the frequency of its adjacent mode is 20 Hz, and the spectrum contains four mode components, so the target frequency f a = 15, the frequency of its adjacent mode f b = 20, and n = 4: Frequency upper limit:

Hz
Each decomposition amount could generally appear Gibbs phenomenon at its boundary while WT.To overcome negative impact of the end effect, we select the location of the first sample point p equals 50, or 100, or more, which depends on the number of sample points.According to Equation ( 5), we let p = 50.The number of sampling point K = 100, and the frequency range is [13.75 Hz,16.25 Hz].Move rectangular window, and continuously calculate three values of A j .The results are shown in Table 2.In Table 2, A j of the first and second modes (15 Hz and 20 Hz) decreases with time.The damping ratio of these two modes is positive, while the others (25 Hz and 32 Hz) increase with time.Those two modes have negative damping ratios which need attention and further analyses.The judgment of the automatic algorithm is consistent with above image observation method.It is not necessary to reconstruct each mode and identify the parameters of each mode component by the proposed method which can greatly improve the efficiency and speed of identification.
In order to validate the accuracy of SWT, we reconstruct all four modes.The results are shown in Figure 4a.The full band signal reconstruction of SWT is shown in Figure 4b. Figure 5 is the error percentage of the reconstructed signal comparing with the original signal.As shown in the figure, the error percentage of the reconstructed signal can be controlled within 1.1%.The result is quite ideal.In Table 2, Aj of the first and second modes (15 Hz and 20 Hz) decreases with time.The damping ratio of these two modes is positive, while the others (25 Hz and 32 Hz) increase with time.Those two modes have negative damping ratios which need attention and further analyses.The judgment of the automatic algorithm is consistent with above image observation method.It is not necessary to reconstruct each mode and identify the parameters of each mode component by the proposed method which can greatly improve the efficiency and speed of identification.
In order to validate the accuracy of SWT, we reconstruct all four modes.The results are shown in Figure 4a.The full band signal reconstruction of SWT is shown in Figure 4b.In Table 2, Aj of the first and second modes (15 Hz and 20 Hz) decreases with time.The damping ratio of these two modes is positive, while the others (25 Hz and 32 Hz) increase with time.Those two modes have negative damping ratios which need attention and further analyses.The judgment of the automatic algorithm is consistent with above image observation method.It is not necessary to reconstruct each mode and identify the parameters of each mode component by the proposed method which can greatly improve the efficiency and speed of identification.
In order to validate the accuracy of SWT, we reconstruct all four modes.The results are shown in Figure 4a.The full band signal reconstruction of SWT is shown in Figure 4b.After that we identify the parameters of the IMTs by Hilbert transform.The results are shown in Table 3.The method is abbreviated as SWT-HT.From Table 3, SWT-HT has high accuracy in the identification of the frequency and attenuation factor.The error of the method is less than 2.5%.

Anti-Noise Ability of Synchrosqueezed Wavelet Transforms
We change the noise level of Equation ( 6).The time-frequency analysis with two different noise levels by SWT is shown in Figure 6.
Energies 2018, 11, x 7 of 16 After that we identify the parameters of the IMTs by Hilbert transform.The results are shown in Table 3.The method is abbreviated as SWT-HT.From Table 3, SWT-HT has high accuracy in the identification of the frequency and attenuation factor.The error of the method is less than 2.5%.

Anti-Noise Ability of Synchrosqueezed Wavelet Transforms
We change the noise level of Equation ( 6).The time-frequency analysis with two different noise levels by SWT is shown in Figure 6.As shown in Figure 6, noise do not influence the precision of time-frequency analysis severely.The behavior is still clear when SNR = 10.The mode parameters of identification in Figure 6b are shown in Table 4.
This proves the SWT method presents good performance against noise to a certain degree, especially for frequency performance.However, the error of the attenuation factor is larger when SNR reaches 10 dB.It is thus indicated that if the noise is strong (SNR is lower than 10 dB), denoising pretreatment is necessary.This proves the SWT method presents good performance against noise to a certain degree, especially for frequency performance.However, the error of the attenuation factor is larger when SNR reaches 10 dB.It is thus indicated that if the noise is strong (SNR is lower than 10 dB), denoising pretreatment is necessary.The simulated signal in Equation ( 6) is decomposed by EMD, the intrinsic mode function (IMF) components and residue are shown in Figure 7.According to the reference [16], EMD fails to extract mode components, in the condition that the frequencies of two adjacent mode components are very close.That is f 1 < f 2 and f 1 /f 2 > 0.5.The simulated signal in Equation ( 6) is decomposed by EMD, the intrinsic mode function (IMF) components and residue are shown in Figure 7.According to the reference [16], EMD fails to extract mode components, in the condition that the frequencies of two adjacent mode components are very close.That is f1 < f2 and f1/f2 > 0.5.The adjacent mode frequency ratios of the above simulated signal are 0.75, 0.8 and 0.78 respectively which are all in the double frequency range.The Fourier spectrum of the IMF1 component in Figure 8 shows that IMF1 still contains four frequency modes.In other words, EMD fails to decompose the signal effectively and modes are mixing.The adjacent mode frequency ratios of the above simulated signal are 0.75, 0.8 and 0.78 respectively which are all in the double frequency range.The Fourier spectrum of the IMF1 component in Figure 8 shows that IMF1 still contains four frequency modes.In other words, EMD fails to decompose the signal effectively and modes are mixing.
Reference [17] proposed an improved method based on EMD, but it is only valid for two components and fails for multiple components.Moreover, the improved EMD will undoubtedly increase the complexity of the signal analysis, which is not conducive to the rapid and accurate extraction and analysis of the signals.The adjacent mode frequency ratios of the above simulated signal are 0.75, 0.8 and 0.78 respectively which are all in the double frequency range.The Fourier spectrum of the IMF1 component in Figure 8 shows that IMF1 still contains four frequency modes.In other words, EMD fails to decompose the signal effectively and modes are mixing.To further compare the anti-mode mixing ability of SWT with EMD, we simulate a signal, which contains two kinds of mode frequencies without damping.Based on Section 3, EMD can't identify the modes of 15 Hz and 20 Hz.T, let one mode frequency have a fixed value of 15 Hz, and the other mode frequency constantly increases to keep away from 20 Hz.In this way, when the other mode frequency reaches 30 Hz, EMD can identify two modes.The result of EMD and their FFT are shown in Figure 9.When the second mode frequency is less than 30 Hz, that is f 1 /f 2 > 0.5, EMD emerge mode-mixing.According to the above, SWT can identify the modes of 15 Hz and 20 Hz, which does improve the ability of anti-mode mixing.Reference [17] proposed an improved method based on EMD, but it is only valid for two components and fails for multiple components.Moreover, the improved EMD will undoubtedly increase the complexity of the signal analysis, which is not conducive to the rapid and accurate extraction and analysis of the signals.
To further compare the anti-mode mixing ability of SWT with EMD, we simulate a signal, which contains two kinds of mode frequencies without damping.Based on Section 3, EMD can't identify the modes of 15 Hz and 20 Hz.T, let one mode frequency have a fixed value of 15 Hz, and the other mode frequency constantly increases to keep away from 20 Hz.In this way, when the other mode frequency reaches 30 Hz, EMD can identify two modes.The result of EMD and their FFT are shown in Figure 9.When the second mode frequency is less than 30 Hz, that is f1/f2 > 0.5, EMD emerge mode-mixing.According to the above, SWT can identify the modes of 15 Hz and 20 Hz, which does improve the ability of anti-mode mixing.

Impact Morlet Mother Wavelet Compact Support on Anti-Mode Mixing Ability of Synchrosqueezed Wavelet Transforms
Based on Section 4.1, SWT can identify the modes of 15 Hz and 20 Hz, so in this section the other mode frequency constantly decreases to keep away from 20 Hz.When the other mode frequency reaches 17 Hz (f 1 /f 2 = 0.88), two modes are almost overlapping in time-frequency spectrum.And there can be no doubt that two modes can't also be recognized effectively when two mode frequencies are too closer.The result of SWT is shown in Figure 10a.
According to the analyzing results of Wu et al. [14] and You et al. [15], the resolution of SWT is determined by the Fourier transform of the mother wavelet compact support ψ.The smaller support is, the higher resolution is.When σ > 5, σ is approximately equal to the central frequency.The central frequency is inversely proportional to the support, so in theory, the value of σ is bigger and resolution of SWT is higher.Changed the value of σ to 25 in the mother wavelet, the time-frequency spectrum of the signal is shown in Figure 10b.The original mixing modes can be accurately identified.It is thus clear that SWT can sift the mode-mixing signal into good mode by increasing σ in mother wavelet.Further analysis indicates that increasing σ is conducive to the extraction of high frequency signals.While σ shouldn't be too big.It will reduce the energy of mode in time-frequency spectrum and influence the readability of the spectrum.Decreasing σ is conducive to the extraction of low-frequency signals.But too small σ will let the mode float in time-frequency spectrum and also influence the readability of the spectrum.Therefore, the appropriate value of σ selected is especially important, and the value of σ needs to be adjusted on the basis of different signal characteristics to obtain the best time-frequency characteristic.In this paper, when σ equals 12, SWT can have a good effect on the SSO signal between 5 Hz and 50 Hz.The original mixing modes can be accurately identified.It is thus clear that SWT can sift the mode-mixing signal into good mode by increasing σ in mother wavelet.Further analysis indicates that increasing σ is conducive to the extraction of high frequency signals.While σ shouldn't be too big.It will reduce the energy of mode in time-frequency spectrum and influence the readability of the spectrum.Decreasing σ is conducive to the extraction of low-frequency signals.But too small σ will let the mode float in time-frequency spectrum and also influence the readability of the spectrum.Therefore, the appropriate value of σ selected is especially important, and the value of σ needs to be adjusted on the basis of different signal characteristics to obtain the best time-frequency characteristic.In this paper, when σ equals 12, SWT can have a good effect on the SSO signal between 5 Hz and 50 Hz.the spectrum.Decreasing σ is conducive to the extraction of low-frequency signals.But too small σ will let the mode float in time-frequency spectrum and also influence the readability of the spectrum.Therefore, the appropriate value of σ selected is especially important, and the value of σ needs to be adjusted on the basis of different signal characteristics to obtain the best time-frequency characteristic.In this paper, when σ equals 12, SWT can have a good effect on the SSO signal between 5 Hz and 50 Hz.

Performance of Synchrosqueezed Wavelet Transforms in IEEE First Benchmark Model
IEEE First benchmark model for computer simulation is a classic modal [18].It divides the steam turbine shafting into six mass blocks.They are: high pressure cylinder (HP), intermediate pressure cylinder (IP), two low pressure cylinders (LPA and LPB), generator mass (GEN) and exciter mass (EXC).Each mass block sequences in an axis.The oscillation frequencies are: 15.7 Hz, 20.2 Hz, 25.6 Hz and 32.3 Hz.The wiring diagram is showed in Figure 11.The generator capacity is 892.4MVA, the transmission system voltage is 500 kV, and the system frequency is selected as 60 Hz.The three-phase short-circuit fault occurs in bus B. The fault is The generator capacity is 892.4MVA, the transmission system voltage is 500 kV, and the system frequency is selected as 60 Hz.The three-phase short-circuit fault occurs in bus B. The fault is triggered at 1.5 s, and it lasts 0.075 s.The simulation time is set to 7 s. Figure 12 shows the generator speed deviation signal that starts form 1.5 s. White noise is added to the signal, and SNR is up to 20 dB.The signal is analyzed by SWT and EMD, respectively.Results are shown in Figure 13.
Energies 2018, 11, x 11 of 16 triggered at 1.5 s, and it lasts 0.075 s.The simulation time is set to 7 s. Figure 12 shows the generator speed deviation signal that starts form 1.5 s. White noise is added to the signal, and SNR is up to 20 dB.The signal is analyzed by SWT and EMD, respectively.Results are shown in Figure 13.In Figure 13c,d, the IMF1 still contains two frequency modes.It means that EMD fails to decompose the signal, and modes mixing problem occurs.However, Figure 13a,b shows that SWT can clearly identify the two modes.
Identification results by SWT-HT are listed in Table 5.They prove that two mode components have negative damping which will seriously endanger the safe operation of the turbine.In Figure 13c,d, the IMF1 still contains two frequency modes.It means that EMD fails to decompose the signal, and modes mixing problem occurs.However, Figure 13a,b shows that SWT can clearly identify the two modes.
Identification results by SWT-HT are listed in Table 5.They prove that two mode components have negative damping which will seriously endanger the safe operation of the turbine.To further discuss the anti-noise ability of SWT, SNR is up to 10 dB.Its time-frequency spectrum is shown in Figure 14 and identification results are shown in Table 6.This proves that SWT can still accurately analyze SSO signal when SNR is up to 10 dB, especially for frequency.To further discuss the anti-noise ability of SWT, SNR is up to 10 dB.Its time-frequency spectrum is shown in Figure 14 and identification results are shown in Table 6.This proves that SWT can still accurately analyze SSO signal when SNR is up to 10 dB, especially for frequency.

5-Machine Power System
Figure 15 shows the wiring diagram of a multi-machine system.S represents infinite power system.There are four turbine generators in this system, and they are called G1-G4 for short.The capacity of generator is 600 MVA, the transmission system voltage is 220 kV and the system frequency is 50 Hz.RT1 + jXT1 = RT4 + jXT4 = 0.0002 + j0.02, RT2 + jXT2 = RT3 + jXT3 = 0.0004 + j0.04, RL +

5-Machine Power System
Figure 15 shows the wiring diagram of a multi-machine system.S represents infinite power system.There are four turbine generators in this system, and they are called G1-G4 for short.The capacity of generator is 600 MVA, the transmission system voltage is 220 kV and the system frequency is 50 Hz.RT1 + jXT1 = RT4 + jXT4 = 0.0002 + j0.02, RT2 + jXT2 = RT3 + jXT3 = 0.0004 + j0.04, RL + jXL = 0.0052 + j0.054, XS = 0.03, Xc = 0.044.This SSO modeling process simplifies the stream turbine shafting to a four-mass model: the exciter mass, generator mass, low-pressure cylinder mass and high-pressure cylinder mass.The simulation condition are as follows: a three-phase short circuit fault occurs in bus A. The total simulation time is 10 s.The start time of fault is at 1 s and fault lasts 16.9 ms. Figure 16 shows the torque diagram of G1.
In Figure 16, the oscillation phenomenon was observed, and the torque curve is divergent.The G1 has three natural frequencies of shafting: 24.65 Hz, 32.39 Hz and 51.10 Hz. Figure 17 shows time-frequency spectrum by SWT.According to the automatic recognition algorithm proposed in Section 3, the general trend of mode components is determined and the results are shown in Table 7.This SSO modeling process simplifies the stream turbine shafting to a four-mass model: the exciter mass, generator mass, low-pressure cylinder mass and high-pressure cylinder mass.The simulation condition are as follows: a three-phase short circuit fault occurs in bus A. The total simulation time is 10 s.The start time of fault is at 1 s and fault lasts 16.9 ms. Figure 16 shows the torque diagram of G1.
In Figure 16, the oscillation phenomenon was observed, and the torque curve is divergent.The G1 has three natural frequencies of shafting: 24.65 Hz, 32.39 Hz and 51.10 Hz. Figure 17 shows time-frequency spectrum by SWT.According to the automatic recognition algorithm proposed in Section 3, the general trend of mode components is determined and the results are shown in Table 7.In Table 7, the values Aj of the first mode (24 Hz) increases with time, and the trend of first mode components is divergent, which needs attention and further analyses, while Aj of the second mode (32 Hz) decreases with time, and does not need its parameters to be identified.
Identification results by SWT-HT in Table 8 prove that the above judgement is correct, and the first mode component has negative damping which will seriously endanger the safe operation of turbine.In order to verify the effectiveness of the proposed method, stochastic subspace identification (SSI) [19,20] is chosen to be compared with it.Because SSI has good precision and noise immunity, and is suitable for its control group.The results are shown in Table 9.
In the Table 9, SSI gets three modes.The third mode component is identified by SSI, but not detected by SWT-HT.Because the order of the initial selection system is different from the system itself, the selected order is greater than the order of the system, the third mode component may be the false component.Even if the third mode component is existent, the damping ratio is positive and relatively high.Thus, the energy of the third mode component is very small and has negligible impact on the system.Furthermore, SWT-HT and SSI have the same high accuracy in the identification of the frequency and damping ratio, and moreover, SSI needs to identify all the oscillation parameters to determine whether the SSO occurs.It is not as fast and flexible as SWT.In Table 7, the values A j of the first mode (24 Hz) increases with time, and the trend of first mode components is divergent, which needs attention and further analyses, while A j of the second mode (32 Hz) decreases with time, and does not need its parameters to be identified.
Identification results by SWT-HT in Table 8 prove that the above judgement is correct, and the first mode component has negative damping which will seriously endanger the safe operation of turbine.In order to verify the effectiveness of the proposed method, stochastic subspace identification (SSI) [19,20] is chosen to be compared with it.Because SSI has good precision and noise immunity, and is suitable for its control group.The results are shown in Table 9.
In the Table 9, SSI gets three modes.The third mode component is identified by SSI, but not detected by SWT-HT.Because the order of the initial selection system is different from the system itself, the selected order is greater than the order of the system, the third mode component may be the false component.Even if the third mode component is existent, the damping ratio is positive and relatively high.Thus, the energy of the third mode component is very small and has negligible impact on the system.Furthermore, SWT-HT and SSI have the same high accuracy in the identification of the frequency and damping ratio, and moreover, SSI needs to identify all the oscillation parameters to determine whether the SSO occurs.It is not as fast and flexible as SWT.

Analysis of the Guohua Jinjie Plant Transmission System
The Guohua Jinjie plant transmission system with series compensation [21] is a typical point to power network, long-distance and high-capacity transmission mode.Because of the higher series compensation, SSO is a practical problem which affects the safe and stable operation of the power grid.The transmission system is shown in Figure 18.Jinxin line and Fuxin line of Xinzhou switching station are installed with 35% series compensation capacitors.

Analysis of the Guohua Jinjie Plant Transmission System
The Guohua Jinjie plant transmission system with series compensation [21] is a typical point to power network, long-distance and high-capacity transmission mode.Because of the higher series compensation, SSO is a practical problem which affects the safe and stable operation of the power grid.The transmission system is shown in Figure 18.Jinxin line and Fuxin line of Xinzhou switching station are installed with 35% series compensation capacitors.

Fugu power plant 2×600MW
Jinjie power plant 2×600MW The turbine shafts of the power plant adopt four-mass model.This model has three torsion frequencies: 13.19 Hz, 22.82 Hz and 28.19 Hz respectively.Among them, the damping of the third mode is weakest and in many ways of operation, there is a risk of SSO.We selected the rotor speed deviation of the generator as the analysis signal.At 0 s, a three-phase short circuit occurs in the Xinzhou-Shibei line and the fault time of duration is 75 ms.The waveform is shown in Figure 19 where the amplitude is per-unit value.After that, the signal is analyzed by SWT, time-frequency spectrum is shown in Figure 20.
There are three oscillation modes of the signal, as shown in Figure 20, and its frequency distribution is around 13 Hz, 22 Hz and 28 Hz, respectively.Using the automatic identification method in Section 3 to analyze its attenuation, let p = 100.It is found that the third mode component (28 Hz) only increases with time, that is A1 = 0.0186, A2 = 0.0193 and A3 = 0.0212, so this mode needs reconstruction.Reconstruction of the third mode is shown in Figure 21.The turbine shafts of the power plant adopt four-mass model.This model has three torsion frequencies: 13.19 Hz, 22.82 Hz and 28.19 Hz respectively.Among them, the damping of the third mode is weakest and in many ways of operation, there is a risk of SSO.We selected the rotor speed deviation of the generator as the analysis signal.At 0 s, a three-phase short circuit occurs in the Xinzhou-Shibei line and the fault time of duration is 75 ms.The waveform is shown in Figure 19 where the amplitude is per-unit value.After that, the signal is analyzed by SWT, time-frequency spectrum is shown in Figure 20.
There are three oscillation modes of the signal, as shown in Figure 20, and its frequency distribution is around 13 Hz, 22 Hz and 28 Hz, respectively.Using the automatic identification method in Section 3 to analyze its attenuation, let p = 100.It is found that the third mode component (28 Hz) only increases with time, that is A 1 = 0.0186, A 2 = 0.0193 and A 3 = 0.0212, so this mode needs reconstruction.Reconstruction of the third mode is shown in Figure 21.

Analysis of the Guohua Jinjie Plant Transmission System
The Guohua Jinjie plant transmission system with series compensation [21] is a typical point to power network, long-distance and high-capacity transmission mode.Because of the higher series compensation, SSO is a practical problem which affects the safe and stable operation of the power grid.The transmission system is shown in Figure 18.Jinxin line and Fuxin line of Xinzhou switching station are installed with 35% series compensation capacitors.

Fugu power plant 2×600MW
Jinjie power plant 2×600MW The turbine shafts of the power plant adopt four-mass model.This model has three torsion frequencies: 13.19 Hz, 22.82 Hz and 28.19 Hz respectively.Among them, the damping of the third mode is weakest and in many ways of operation, there is a risk of SSO.We selected the rotor speed deviation of the generator as the analysis signal.At 0 s, a three-phase short circuit occurs in the Xinzhou-Shibei line and the fault time of duration is 75 ms.The waveform is shown in Figure 19 where the amplitude is per-unit value.After that, the signal is analyzed by SWT, time-frequency spectrum is shown in Figure 20.
There are three oscillation modes of the signal, as shown in Figure 20, and its frequency distribution is around 13 Hz, 22 Hz and 28 Hz, respectively.Using the automatic identification method in Section 3 to analyze its attenuation, let p = 100.It is found that the third mode component (28 Hz) only increases with time, that is A1 = 0.0186, A2 = 0.0193 and A3 = 0.0212, so this mode needs reconstruction.Reconstruction of the third mode is shown in Figure 21.

Figure 1 .
Figure 1.Parameter extraction and calculation flow chart of the subsynchronous oscillation (SSO) mode based on synchrosqueezed wavelet transforms (SWT).
transform is shown in Equation(8).

Figure 1 .
Figure 1.Parameter extraction and calculation flow chart of the subsynchronous oscillation (SSO) mode based on synchrosqueezed wavelet transforms (SWT).
t is time and ω corresponds to frequency.The time-frequency resolution is regulated by changing parameter σ.Through a large amount of simulated examples, it indicates that, Morlet wavelet has good performance of resolution both in time and frequency domain for the SSO signals when σ = 12.
Figure 5 is the error percentage of the reconstructed signal comparing with the original signal.As shown in the figure, the error percentage of the reconstructed signal can be controlled within 1.1%.The result is quite ideal.(a) (b)

Figure 4 .
Figure 4. IMT and full band signal reconstruction: (a) reconstruction of IMT; and (b) full band signal reconstruction.

Figure 4 .
Figure 4. IMT and full band signal reconstruction: (a) reconstruction of IMT; and (b) full band signal reconstruction.
Figure 5 is the error percentage of the reconstructed signal comparing with the original signal.As shown in the figure, the error percentage of the reconstructed signal can be controlled within 1.1%.The result is quite ideal.(a) (b)

Figure 4 .
Figure 4. IMT and full band signal reconstruction: (a) reconstruction of IMT; and (b) full band signal reconstruction.

Figure 6 .
Figure 6.Time-frequency analysis of signal with different white noise levels by SWT: (a) SNR is 15 dB; and (b) SNR is 10 dB.

Figure 6 .
Figure 6.Time-frequency analysis of signal with different white noise levels by SWT: (a) SNR is 15 dB; and (b) SNR is 10 dB.
Benchmark Model IEEE First benchmark model for computer simulation is a classic modal [18].It divides the steam turbine shafting into six mass blocks.They are: high pressure cylinder (HP), intermediate pressure cylinder (IP), two low pressure cylinders (LPA and LPB), generator mass (GEN) and exciter mass (EXC).Each mass block sequences in an axis.The oscillation frequencies are: 15.7 Hz, 20.2 Hz, 25.6 Hz and 32.3 Hz.The wiring diagram is showed in Figure 11.

4. 4 .
Performance of Synchrosqueezed Wavelet Transforms in IEEE First Benchmark Model IEEE First benchmark model for computer simulation is a classic modal [18].It divides the steam turbine shafting into six mass blocks.They are: high pressure cylinder (HP), intermediate pressure cylinder (IP), two low pressure cylinders (LPA and LPB), generator mass (GEN) and exciter mass (EXC).Each mass block sequences in an axis.The oscillation frequencies are: 15.7 Hz, 20.2 Hz, 25.6 Hz and 32.3 Hz.The wiring diagram is showed in Figure 11.

Figure 17 .
Figure 17.Time-frequency distribution of torque of G1.

Figure 17 .
Figure 17.Time-frequency distribution of torque of G1.

Figure 18 .
Figure 18.Jinjie plant transmission system with series compensation.

Figure 19 .
Figure 19.Signal of rotor speed deviation of the Jinjie plant.

Figure 18 .
Figure 18.Jinjie plant transmission system with series compensation.

Figure 18 .
Figure 18.Jinjie plant transmission system with series compensation.

Figure 19 .
Figure 19.Signal of rotor speed deviation of the Jinjie plant.Figure 19.Signal of rotor speed deviation of the Jinjie plant.

Figure 19 .
Figure 19.Signal of rotor speed deviation of the Jinjie plant.Figure 19.Signal of rotor speed deviation of the Jinjie plant.

Table 2 .
Automatic recognition algorithm of intrinsic mode types (IMTs).

Table 2 .
or 100, or more, which depends on the number of sample points.According to Equation (5), we let p = 50.The number of sampling point K = 100, and the frequency range is[13.75Hz,16.25Hz].Move rectangular window, and continuously calculate three values of Aj.The results are shown in Table2.Automatic recognition algorithm of intrinsic mode types (IMTs).

Table 2 .
or 100, or more, which depends on the number of sample points.According to Equation (5), we let p = 50.The number of sampling point K = 100, and the frequency range is[13.75Hz,16.25Hz].Move rectangular window, and continuously calculate three values of Aj.The results are shown in Table2.Automatic recognition algorithm of intrinsic mode types (IMTs).

Table 3 .
The parameters of the IMTs.

Table 3 .
The parameters of the IMTs.

Table 4 .
Parameters of IMTs with 10 dB noise.

Table 4 .
Parameters of IMTs with 10 dB noise.

Table 5 .
Results of identification by SWT-HT.

Table 5 .
Results of identification by SWT-HT.

Table 6 .
SWT results with two different level noise.

Table 6 .
SWT results with two different level noise.

Table 7 .
Distinguish attenuation of IMT by the automatic recognition algorithm.

Table 8 .
Results of identification by SWT-HT.

Table 7 .
Distinguish attenuation of IMT by the automatic recognition algorithm.

Table 8 .
Results of identification by SWT-HT.

Table 9 .
Mode parameters of the simulated signal base on SSI.

Table 9 .
Mode parameters of the simulated signal base on SSI.

Table 9 .
Mode parameters of the simulated signal base on SSI.