An attractive advantage of sparse component analysis (SCA) is that it can transform the time-domain signal into the frequency domain through appropriate linear transformation method, such as short-time Fourier transform, wavelet transform, etc., making the observed measurement signals sparse in the frequency domain. The sparsity of observation signal means that most time-frequency points are zero, and only a few time-frequency points have large values in the time-frequency domain (Ye et al., 2018). In this paper, one-dimensional time-domain signal is transformed into two-dimensional time-frequency-domain signal by short-time Fourier transform, and the rectangular window function is selected. And 1) could be further expressed as (Yi et al., 2017): [ x 1 ( t , f ) x 2 ( t , f ) ⋮ x p ( t , f ) ] = [ a 1 a 2 ⋯ a Q ] [ s 1 ( t , f ) s 2 ( t , f ) ⋮ s Q ( t , f ) ] = [ a 11 ⋯ a 1 Q a 22 ⋯ a 2 Q ⋮ ⋱ ⋮ a P 1 ⋯ a P Q ] [ s 1 ( t , f ) s 2 ( t , f ) ⋮ s Q ( t , f ) ] (3) where [ a 1 a 2 ⋯ a Q ] is the column vector of mixing matrix A . a P Q is the attenuation factor for the Q-th source signal reaching the P-th observation point.

When each source signal is sparse and disjoint, the observed signal has at most one source signal at a certain time-frequency point. Assuming that only the source signal s j ( t v , f v ) exists at the time-frequency point ( t v , f v ) , 3) can be expressed as: [ x 1 ( t v , f v ) x 2 ( t v , f v ) ⋮ x p ( t v , f v ) ] = [ a 1 a 2 ⋯ a Q ] [ 0 ⋮ s j ( t v , f v ) ⋮ 0 ] = [ a 11 ⋯ a 1 Q a 22 ⋯ a 2 Q ⋮ ⋱ ⋮ a p 1 ⋯ a P Q ] [ 0 ⋮ s j ( t v , f v ) ⋮ 0 ] (4)

Then, (4) can be simplified as: x 1 ( t v , f v ) a 1 j = x 2 ( t v , f v ) a 2 j = ⋯ = x P ( t v , f v ) a P j = s j ( t v , f v ) (5)

From a geometrical point of view, if the real part or imaginary part of x 1 ( t , f ) , x 2 ( t , f ) is used as the horizontal axis and vertical axis to draw a scatter plot, x 1 ( t v , f v ) / x 2 ( t v , f v ) = a 1 j / a 2 j is constant for the sampling times belonging to the source signal s j . That is, the time-frequency points belonging to the source signal s j will determine a straight line. Similarly, the time-frequency points of different source signals will be clustered into corresponding clustering lines. The direction vector of these clustering lines is the attenuation coefficient ratio of each source signal mixed to the observation point, which is the column vector a Q of the mixing matrix A . And the mixing matrix A can be estimated by the direction of each clustering line in the scatter diagram.

Based on the obtained mixing matrix A , the L1 norm is used to solve the underdetermined equation, then the system source signal is extracted.

The premise of estimating the mixing matrix A is that the time-frequency points of the mixed signal can form a clustering line with certain directionality in the real or imaginary part scatter diagram. In this paper, the fuzzy C-means method is used to estimate the mixing matrix A by calculating the clustering center of the clustering line, which is suitable for the blind source decomposition problem with any number of measurement signal channels. However, when there are too many time-frequency points, the clustering direction in the scatter diagram is complicated, resulting in much calculation time and poor identification accuracy.

The source signal propagates outward in the energy form, however, the energy will be attenuated during the propagation process, the final energy of the same source signal reaching each observation is a little different. Therefore, an algorithm combing the frequency energy peak point method (Jin et al., 2021) and the fuzzy C-means method is proposed to estimate the mixing matrix A .

The main principle of the proposed method to estimate the mixture matrix A is: since the energy of each source signal is the largest near its own frequency point, the frequency value of each source can be determined by using the frequency peak method. The clustering straight line direction of the corresponding time-frequency point at the maximum energy is the clustering direction of each source signal. Then, the real part or imaginary part time-frequency dispersion points corresponding to each peak frequency point are normalized and mapped to the unit circle. Finally, the clustering center is calculated by using the fuzzy C-means algorithm, so as to obtain the mixing matrix A .

In the time-frequency domain, the energy distribution of a single measurement signal is firstly calculated, and then the energy of multiple channels at the same frequency point is added, which is: E ( f ) = ∑ i = 1 m ∫ t 1 t 2 ( ( Re ( x i ( t , f ) ) ) 2 + ( Im ( x i ( t , f ) ) ) 2 ) d t (6) where E(f) is the sum of the energy of all observation signals at each frequency point; x _i is the i-th time-frequency point; Re ( x i ( t , f ) ) and Im ( x i ( t , f ) ) are the real part and imaginary part of the time-frequency point respectively; m is the number of the observation signals.

SCA algorithm utilizes the sparsity of the source signal in time-frequency domain to recover the source signal. The sparse solution is obtained by establishing and solving the optimization problem (P ₁ ) in (7): ( P 1 ) : { s ∼ ( t , f ) = min | s ( t , f ) | l 1 s u b j e c t t o A s ( t , f ) = x ( t , f ) (7) where s ∼ ( t , f ) is the estimation of the source signal s ( t , f ) ; ( P 1 ) is a convex optimization problem. The solution of the optimization problem is the L1 norm minimum solution (Donoho, 2006).

The power system requires accurate generator parameters and mathematical models when performing small-signal stability calculation and analysis. The kinetic energy of the generator is an important parameter of the generator, which is determined by the rotational inertia of the generator rotor. The rotational inertia refers to the inertia of the generator rotor during the rotation process, which affects the electromechanical transient process of each generator and is of great significance to the stability and safety of the power system. When the system is disturbed, the generator rotor will accelerate or decelerate due to the unbalanced torque, and the rotational inertia can be expressed as the ability of the system to return to the original stable operation. The smaller the rotational inertia, the worse the generator’s resistance to disturbance. When the proposed method is applied to the actual system, it is necessary to use the rotational inertia of the generator in each region as a quantitative index, and select the measurement signals of the generator with the smallest rotational inertias as the input signals of SCA algorithm, which can greatly reduce the number of the input signals and effectively improve the calculation speed.

The low-frequency oscillation signal of power system is essentially the multimodal coupling response, which can be expressed by (8). The oscillation frequency has a large deviation in a period of time, which shows obvious nonlinear characteristics: x ( t ) = Φ q ( t ) = ∑ i = 1 n φ i q i ( t ) (8) where Φ =[ φ 1 , φ 2 ,⋯ φ n ] is the mode shape matrix, the mode shape of the i-th column φ i is related to the i-th modal response vector q i ( t ) , and can also be transformed into the form of (9) (Yang and Satish, 2013): q ( t ) = Φ − 1 x ( t ) (9)

Comparing (1) and (8), the viewpoints of characteristic parameter identification of electromechanical oscillation mode and blind source separation are consistent, that is, the perturbed dynamic response of multi-modal coupling is regarded as a linear group of multiple single-modal responses. The single-modal response q(t) is equivalent to a special form of the source signal s(t), and the mixing matrix A contains the information of the modal matrix, satisfying Φ = A .

In this paper, Hilbert transform (HT) is used to extract the oscillation frequency and damping ratio (Gibbard and Vowles, 2010). Performing the HT on a single-mode oscillation component q _i (t), we can obtain: h ( t ) = 1 π ∫ − ∞ + ∞ q ( t ) t − τ d t (10) Where τ is the integral variable.

Introduce an analytic function consisting of q(t) and h(t): Z ( t ) = q ( t ) + j h ( t ) = A ( t ) e − j φ ( t ) (11) Where A(t) is the envelope of the signal after Hilbert transform and φ ( t ) is the instantaneous phase, namely: { A ( t ) = [ q 2 ( t ) + h 2 ( t ) ] 1 / 2 φ ( t ) = tan − 1 [ h ( t ) q ( t ) ] (12)

The instantaneous frequency can be obtained: f ( t ) = 1 2 π × d A ( t ) d t (13)

The system oscillation signal can be deduced from the generator swing equation, which can be expressed as: x ( t ) = = ∑ i = 1 n φ i q i ( t ) = ∑ i = 1 n φ i u i e − λ i t c o s ( ω d i t + θ i ) (14) Where q _i (t) is a cosine curve that decays exponentially; u i and θ i are the amplitude and phase respectively, which are determined by the initial conditions; λ i is the decay coefficient; ω d i is the oscillation frequency.

For a damped oscillaion, the single-modal response signal q _i (t) can be expressed as: q i ( t ) = u i e − ξ ω 0 t c o s ( ω 0 1 − ξ 2 t + θ i ) (15)

By comparing 14) and (15), we can obtain: { ξ ω 0 = λ i ω d i = ω 0 1 − ξ 2 (16)

The damping ratio of the signal is calculated as: ξ = − λ i ω d i 2 + λ i 2 (17)

The modal parameter extraction of an electromechanical oscillation based on the SCA algorithm is shown in Figure 1:

The five single-mode oscillation signals shown in Eq. 18 are constructed with frequencies of 0.6 Hz, 1.1 Hz, 0.4 Hz, 0.8 Hz, and 0.2 Hz, in that order: { S 1 = e − 0.1 t ⁡ sin ( 2 π × 0.6 t + π 6 ) S 2 = e − 0.05 t ⁡ sin ( 2 π × 1.1 t + 0 ) S 3 = e − 0.15 t ⁡ cos ( 2 π × 0.4 t + π 3 ) S 4 = e − 0.25 t ⁡ sin ( 2 π × 0.8 t + π 4 ) S 5 = e − 0.08 t ⁡ sin ( 2 π × 0.2 t + π 5 ) (18)

Using a three-dimensional random mixing matrix A = [ 0.6 0.4 0.53 0.45 0.41 0.44 0.65 0.5 0.68 0.48 0.5 0.45 0.35 0.67 0.61 ] , it is mixed into three coupled oscillating signals with non-smooth characteristics as the test signals. The three test signals are sampled separately with the sampling frequency of 0.01Hz, and the sampling time is 15s. The spectrum analysis of the three test signals is performed using window Fourier transform (WFT), as shown in Figure 2.

Short time Fourier transform (STFT) is used to transform the three groups of test signals into time-frequency domain, and the window function is chosen to be a rectangular window of length 10s, and the window was shifted by four points each time. Figure 3 shows the scatter plot of the real part of the test signal x. From Figure 3, we can know that there are five obvious clustered straight lines, indicating that there are five source signals in the test signal x.

According to (6), the total energy of the three groups of test signals at each frequency point is calculated, and the frequency corresponding to each energy peak is detected using the frequency energy peak point method, as shown in Figure 4. The normalized clustering center of the five straight lines is calculated using the fuzzy C-mean algorithm, and the mixing matrix A can finally be accurately estimated as: A = [ 0.4826 0.6458 0.6662 - 0.4350 0 . 4478 0 . 5492 0 . 6030 0 . 4770 - 0 . 6257 0 . 7329 0 . 6777 0 . 4243 0 . 5484 - 0 . 6230 0 . 5040 ] .

The L1 norm minimization method was used to separate the five source signal components, and the time domain and frequency domain distribution of the five source signals is shown in Figure 5. As shown in Figure5, each component has the characteristics of periodic oscillation, and the calculated frequencies are basically the same as the actual frequencies. Meanwhile, the proposed algorithm decomposes the multi-modal coupled oscillation signal into multiple source signals with the single frequency, without modal aliasing. The detected frequencies are arranged in the order from low to high, which effectively solves the problem of uncertainty in the order of source signals.

Hilbert transform was utilized to extract modal parameters of the separated components, and the oscillation frequencies and damping ratios are obtained, as shown in Table 1. It can be seen from Table 1 that the results of the extracted modal parameters are quite close to their true values, indicating the effectiveness of the proposed method.

In this section, the modified IEEE 14-generater 5-area simplified system is used as an example (Ding et al., 2019), as shown in Figure 6. There are five areas in this system, but since area one is strongly coupled with area 2, the system can be considered to have four main areas, which are defined as Area 1, Area 2, Area 3, and Area 4, respectively. Area one contains five generators, namely G101, G201, G202, G203, and G204; Area2 contains two generators, G301, G302; Area3 contains four generators, namely G401, G402, G403, G404; Area4 contains G501, G502, G503. And two aggregated PMSG-based wind farms are connected to Buses 504 and 201, respectively in Figure 6. Furthermore, the parameters of the two wind farms are configured in accordance with Literature (Dejian et al., 2022a; Dejian et al., 2022b). Moreover, the active outputs of the SGs in the same area are reduced accordingly to keep the operating point the same as that before the SG reduction.

Inter-area mode is more likely to excite a poorly damped oscillation, which involves more generators and has a wide range of influence, so it is usually considered the mode of most interest. The 14-generator test system was linearized around an operating point, and small-signal stability analysis (SSSA) results show that the system has three inter-area modes. as shown in Table 2.

In this paper, the rotational inertia of the generators in each area of the system is used as a quantitative index, and the angular frequency of the generator with the smallest rotational inertia in each area is selected as the input signal. Figure 7 shows the rotational inertia of each generator. According to the signal selection principle and the inertia results in Figure 7, generator G203 is selected as the representative of Area1, in which its inertia is lower than the other generators in the same area. Similarly, the signal of generator G301 is selected as representative of Area2, generator G401 is se-lected as representative of Area3, and generator G501 is selected as representative of Area4.

In order to verify the effectiveness of the proposed method for identifying the interarea modes, a representative two-phase grounded short circuit is applied at bus 506, the sampling frequency of 100 Hz and the sampling time is 50 s. And the simulation was carried out using power system toolbox (PST). The angular frequency of each generator within 10s after the fault is used as the input signal in this paper. The input signals of generator G203, generator G301, generator G401, and generator G501 were first detrended, and the time domain angular frequencies of these four generators were transformed into the frequency domain using WFT. Three interarea oscillation modes could be identified, with frequencies of 0.34 Hz, 0.51 Hz, and 0.68 Hz, respectively. Figure 8 shows the time-domain distributions and the corresponding frequency-domain distributions of the angular frequencies of the four generators.

The four input signals are decomposed using SCA algorithm, and three single-mode oscillation components can be obtained. The time-domain and frequency-domain distributions of each single-mode component are shown in Figure 9. HT was used to calculate the oscillation frequency and damping ratio, and the obtained modal parameter results were compared with WFT results and SSSA results, as shown in Table 3. From Figure3, we can know that all the three inter-area modal parameters can be extracted by using only four input signals of generator G203, G301, G401 and G501, which effectively saves calculation time.

While identifying the oscillation frequency and damping ratio of the system, SCA algorithm can also estimate the modal matrix utilizing the clustering characteristics.

The modal matrix is estimated as: [ 0.1902 0.4766 - 0.6595 - 0.0357 - 0.1971 0.451 0.7167 - 0.0124 0.5477 - 0.4067 - 0.0325 - 0.0364 - 0.7905 - 0.6355 - 0.2221 - 0.9985 ]

The mode shape corresponding to the three inter-area oscillation modes are shown in Figure 10. For Mode 1, generator 401 oscillates against generator 501. In this paper, the input signal is the representative of each region, so it can also be considered that mode one is the inter-area oscillation mode, with generator 401 in Area3 oscillates against generator 501 in Area4. By the same token, the oscillation area clusters of mode two is Area two and Area three against Area 4, the oscillation area clusters of mode three is Area two against Area 3. The estimated oscillation frequencies, damping ratios and mode shapes using the proposed method in this paper are summarized in Table 4.

In this section, a set of real PMU measurements from north power grid of China was used to further illustrate the adaptability of the algorithm in practical system. The system comprises two major areas connected through four 500-kV transmission lines. When the maintenance schedule requires, heavy power is delivered between the two areas, the electromechanical oscillation is one of the major threats to the power grid stability. The data are collected from four PMUs located in the two areas. The input signals were recorded for a duration of 15 s with the sampling frequency of 30 Hz. The time-domain oscillation waveforms and the corresponding frequency-domain distributions are shown in Figure 11. As can be seen from Figure 11, the practical power grid contains an inter-area oscillation mode with the frequency of 0.28 Hz.

The SCA algorithm was used to extract oscillation features from PMU measurement data, the decomposed time-domain source signals and their spectrum are shown in Figure 12. Then the HT algorithm was applied to the decomposed source signals to extract oscillation frequencies and damping ratios, as shown in Table 5. As can be seen from Table 5, the proposed method identified the oscillation mode with the frequency of 0.28 Hz, which is consistent with the spectrum analysis, and the corresponding mode shape is shown in Figure 13. From Figure 13, it can be seen that PMU2 is mainly involved in the oscillation of this mode.

To sum up, it can be seen that the proposed algorithm in this paper can not only effectively extract the modal parameters from the test signal and the simulated signal, but also has good applicability to the real measured signal.