Analysis of the Propagation Characteristic of Subsynchronous Oscillation in Wind Integrated Power System

This paper proposes oscillation propagation factors to analyze power oscillations caused by the interharmonics of doubly fed induction generators (DFIG) at different points in the power system. First, a dynamic model of the DIFG is built, including the asynchronous generator, its transmission system, converters and the control systems. Then, the state space expression is formed by deducing the input and output matrices. From this, the oscillation propagation factor is proposed and denoted to exhibit the propagation mechanism of interharmonics in the view of frequency domain, by deducing the multi-input-multi-output transfer functions matrix. Along with this, the sensitivity of propagation is calculated for adjusting the parameters to block the oscillation propagating path. Finally, the modified four machine system with two DFIGs and the New-England 39 bus system with two DFIGs is used as a test system to verify the effectiveness of the oscillation propagation factor. From this the simulation results demonstrate that the subsynchronous interharmonics of DFIGs injected into the grid will propagate to the different points of the system and results in oscillation of the power. The oscillation propagation factor could quantize the oscillation magnitude propagating from one point to other point in the wind integrated power system.


Introduction
New energy sources such as wind power and photovoltaic cells have become an important part of the power system.With the continuous increase of the proportion of new energy sources, the stability of the power system caused by the access of new energy has attracted wide attention [1][2][3].Wind power is one of the best solutions for new energy generation [4].At present, doubly fed induction generator (DFIG) and permanent magnet synchronous generator (PMSG) are the mainstream machine types for wind power generation [5].In order to improve the stability and reliability of the wind integrated power system, it is necessary to make a much deeper study on the phenomenon of subsynchronous oscillation in the doubly-fed induction generation system.
The DFIG rotor-side converters and grid-side converters contain high power electronics.When the random variation of wind speed causes the frequency variation of the DFIG rotor excitation current, due to the switching effect during the conversion or regulation process, the rotor current frequency and the grid frequency will cause the subsynchronous frequency interharmonics in the DFIG current [6,7].Reference [7] shows that the DFIG can act as a source of significant subsynchronous interharmonic current emissions.When the interharmonic current is injected into the grid, it will cause voltage, current and power to oscillate, causing various faults and the instability of the power system [8].Subsynchronous oscillation caused by DFIG has been studied in terms of damping and impedance, such as eigenvalue analysis [9,10], frequency scanning method [11,12], and equivalent impedance [13][14][15].References [9,10] analyze and control the subsynchronous oscillation based on eigenvalue analysis.Frequency scanning method is performed to investigate subsynchronous control interaction [11,12].An impedance-based analysis approach is used to characterize subsynchronous oscillation [13][14][15].At present, there is no literature to quantitatively study the propagation law of DFIG subsynchronous interharmonics.
The rest of the paper is structured as follows: Section 2 gives the dynamic model of DFIG.In Section 3, a small signal state space model for the power system is derived from the DFIG model, the synchronous generator state equation and the network equation considering the load.Then we calculate the system transfer function matrix and define the oscillation propagation factor to analyze the system power oscillation.Section 4 validates the effectiveness of the oscillation propagation factor based on the 4-machine 2-area system and the New-England 39 bus system.The discussion and conclusion are provided in Sections 5 and 6 respectively.

Dynamic Model of Wind Integrated Power System
The doubly fed induction generator structure is shown in Figure 1.It mainly includes asynchronous generator, transmission system, rotor-side converter, grid-side converter and control system [16,17].
Energies 2018, 11, x 2 of 17 frequency and the grid frequency will cause the subsynchronous frequency interharmonics in the DFIG current [6,7].Reference [7] shows that the DFIG can act as a source of significant subsynchronous interharmonic current emissions.When the interharmonic current is injected into the grid, it will cause voltage, current and power to oscillate, causing various faults and the instability of the power system [8].Subsynchronous oscillation caused by DFIG has been studied in terms of damping and impedance, such as eigenvalue analysis [9,10], frequency scanning method [11,12], and equivalent impedance [13][14][15].References [9,10] analyze and control the subsynchronous oscillation based on eigenvalue analysis.Frequency scanning method is performed to investigate subsynchronous control interaction [11,12].An impedance-based analysis approach is used to characterize subsynchronous oscillation [13][14][15].At present, there is no literature to quantitatively study the propagation law of DFIG subsynchronous interharmonics.
The rest of the paper is structured as follows: Section 2 gives the dynamic model of DFIG.In Section 3, a small signal state space model for the power system is derived from the DFIG model, the synchronous generator state equation and the network equation considering the load.Then we calculate the system transfer function matrix and define the oscillation propagation factor to analyze the system power oscillation.Section 4 validates the effectiveness of the oscillation propagation factor based on the 4-machine 2-area system and the New-England 39 bus system.The discussion and conclusion are provided in Sections 5 and 6 respectively.

Dynamic Model of Wind Integrated Power System
The doubly fed induction generator structure is shown in Figure 1.It mainly includes asynchronous generator, transmission system, rotor-side converter, grid-side converter and control system [16,17].

Rotor-side converter
Grid-side converter In Figure 1, s U is the stator voltage of DFIG, s I is the stator current; r U is the rotor voltage, r I is the rotor current; 1 U is the AC side voltage of the grid-side converter, and 3 r I is the AC side current of the grid-side converter.
3 r L is the filter between the AC side of the grid-side converter and the line.Ignoring the transient process of the motor stator flux, according to the DFIG flux and voltage equations [16,18], the dynamic model of the asynchronous generator is (1).In Figure 1, U s is the stator voltage of DFIG, I s is the stator current; U r is the rotor voltage, I r is the rotor current; U 1 is the AC side voltage of the grid-side converter, and I r3 is the AC side current of the grid-side converter.L r3 is the filter between the AC side of the grid-side converter and the line.X r3 is the filter reactance.C is the capacitance value of the intermediate capacitor, and U dc is the DC voltage of the capacitor.U W is the DFIG bus voltage; I W is the DFIG bus injection current.The positive direction of the physical quantity uses the generator convention.
Ignoring the transient process of the motor stator flux, according to the DFIG flux and voltage equations [16,18], the dynamic model of the asynchronous generator is (1).
Energies 2019, 12, 1081 3 of 16 where E d and E q are the d-axis and q-axis components of the stator transient voltage, respectively; U dr and U qr are the d-axis and q-axis components of the rotor winding voltage, respectively; I ds and I qs are the d-axis and q-axis components of the stator winding current, respectively; s l is the slip; R r is the rotor winding resistance; X m is the excitation reactance; X rr = X m + X r , where X r is the rotor leakage reactance; T J is the inertia time constant of the whole drive; T e is the electromagnetic torque; T m is the mechanical torque; D is the damping coefficient; s 0 is the slip at steady state, ω 0 is the synchronous angular frequency.
Considering the DC voltage dynamics, the dynamic model of the intermediate capacitor C in Figure 1 is: where P r is the active power from the rotor-side converter to the intermediate capacitor; P r3 is the active power from the intermediate capacitor to the grid-side converter.
The dynamic model of the filtered reactance is: where I dr3 and I qr3 , are the d-axis and q-axis components of the current on the AC side of the grid-side converter, respectively; U d1 and U q1 are the d-axis and q-axis components of the AC side voltage of the grid-side converter, respectively; The control objective of the rotor-side converter is to achieve decoupling independent control of the active power and reactive power on the stator side by controlling the excitation voltage [19].The control block diagram of the rotor-side converter is available from [17].
The dynamic model of the rotor-side converter is: where x 1 , x 2 x 3 and x 4 are the introduced state variables; K i1 , K i2 K i3 and K i4 , are the integral coefficients of the corresponding PI controller; P s and Q s are the active and reactive powers of the stator winding, respectively; I dr and I qr are the d-axis and q-axis components of the rotor winding current, respectively.The superscript * indicates the reference value of the corresponding physical quantity.The control objective of the grid-side converter is to maintain the intermediate capacitor voltage as stable, and control the reactive power exchange between the rotor side and the grid to zero [20].The control block diagram of the grid-side converter is available from [16].
The dynamic model of the grid-side converter is: where x 5 , x 6 and x 7 are the introduced state variables; K i5 , K i6 and K i7 are the integral coefficients of the corresponding PI controller.
According to (6), the mathematical model of all DFIGs in the wind integrated power system is (7). where T is the DFIG current vector; n is the number of DFIGs in wind integrated power system.The synchronous generator and its excitation system use a fourth order model.The linearized synchronous generator state equation is (8).
where ∆X G is the synchronous generator state vector; ∆U G is the bus voltage vector of synchronous generator; ∆I G is the synchronous generator injection current vector; A G , B G , C G and D G are state matrix, input matrix, output matrix, and direct transfer matrix, respectively.The load uses the constant impedance model.The power network equation considering load is (9) where  is the node admittance matrix considering the load.By combining (7), ( 8) and ( 9), ( 10) can be obtained. where Energies 2019, 12, 1081 5 of 16 The active power injected into the grid by DFIG is P W : P W = P s + P r3 P s = I ds U ds + I qs U qs P r3 = I dr3 U ds + I qr3 U qs (11) By linearizing (11), we get, According to the mathematical model of DFIG, we get (13): where X = X ss − X 2 m /X rr is the transient reactance of DFIG.X ss = X m + X s , X s is the stator leakage reactance.Substituting (13) into (12), the active power of DFIG can be written as (14).
By combining ( 9), ( 10) and ( 14), the small-signal state space model of the whole system can be obtained: where ∆y = [∆P W1 , ∆P W2 , • • • , ∆P Wk ] T is the output vector; k is the number of generators in wind integrated power system; ∆u is the input vector; A = A − BD −1 C is the state matrix; B is the input matrix; C is the output matrix; D is the direct transfer matrix.

Oscillation Propagation Factor
The Laplace transform is performed on the (15) and the initial condition is considered to be zero: Eliminate ∆X(s) in ( 16), then the transfer function matrix between the output variable ∆y and the input variable ∆u is: where W(s) is a matrix function.The matrix element W ij (s) represents the sinusoidal response from the jth element ∆u j of the input variable to the ith element ∆y i of the output variable.For the jth element ∆u j of the input variable with angular frequency ω, the amplitude of the ith element ∆y i of the output variable is W ij (jω) , the phase shift is ∠W ij (jω), and the total response ∆y i is equal to the linear sum of the individual input responses: After the DFIG interharmonics are injected into the grid, the resulting disturbances will respond at different points in the power system.In the multi-input-multi-output systems, we define oscillation propagation factors to characterize the effects of the input variable disturbances with different frequencies on the oscillations amplitude of the output variables.The oscillation propagation factor is the ratio of an element of the output vector to the input vector.The oscillation propagation factor k i (s) of the ith element of the output variable is denoted and calculated by (19) [21].
In ( 19), the unit of k i (s) is dB; |∆y i | is the magnitude of the ith element of the output vector and ∆u 2 is the 2 norm of the input vector.The larger the oscillation propagation factor, the stronger is the influence of the system disturbance at this point, and the larger the amplitude of the output vector element oscillation.Equation ( 17) can be expressed as: It can be seen that when the elements W ij (s) of the transfer function matrix W(s) do not have zero pole cancellation, the pole of the transfer function W ij (s) is the same as the eigenvalue of the system state matrix.Consider the zero-polar format of the transfer function W ij (s): Substitute s = jω into (21) and replace pole p u with the eigenvalue of the state matrix A. Then, we obtain the magnification W ij (jω) of the transfer function W ij (jω), that is the amplitude gain of output.
When the angular frequency ω of the input variable approaches the angular frequency ω u of the system, the factor |jω − σ u − jω u | in the denominator of (21) decreases, so the amplification factor W ij (s) increases and the oscillation propagation factor becomes larger.And when ω = ω i , the smaller σ i , that is, the smaller the damping ratio, the larger the oscillation propagation factor.
The sensitivity D(ω) of oscillation propagation factor to system parameters is calculated as follows: where k i (s) is the oscillation propagation factor of generators; ω is the angular frequency; p is the system parameters, such as the active power of the DFIGs and the PI controller parameters of these DFIGs.The sensitivity can be solved by perturbation theory.Sensitivity indicates the influence of system parameters changes on the system oscillation propagation factor.The positive and negative sensitivity reflects the increase and decrease of the oscillation propagation factor when the system parameters fluctuate.The magnitude of the sensitivity reflects the degree of influence of system parameters on the oscillation propagation factor.

Simulation Results and Analysis
In order to verify the validity of the oscillation propagation factor proposed in this paper, the propagation of the subsynchronous oscillation in the test system is analyzed by the oscillation propagation factor.When the interharmonic frequency is changed, and when the system damping ratio is changed, the oscillation propagation factor is calculated to compare the strength of the subsynchronous oscillation of each generator.

Two Area Four Machine System
The modified 4-machine 2-area system [22] with two DFIGs is shown in Figure 2. The DFIG parameters are listed below.Rotor winding resistance R r = 0.0013 p.u. Excitation reactance X m = 2.6 p.u. Rotor leakage reactance X r = 2.9 p.u. Filter reactance X r3 = 5 p.u. Capacitance C = 13.29 p.u. Inertia time constant T J = 3.4 p.u.The small-signal stability analysis was performed on the system of Figure 2, and seven oscillation modes were obtained, as shown in Table 1.The active power of DFIGs and the input vector are P DFIG1 = 0.4p.u., P DFIG2 = 0.4p.u.,u = 0 in Table 1.Eigenvalues for each mode are given in the form λ u = σ u ± jω u , from which the frequencies f u = ω u 2π and the damping ratios The real part of the system eigenvalue is negative, and the system is small-signal stable.
sensitivity reflects the increase and decrease of the oscillation propagation factor when the system parameters fluctuate.The magnitude of the sensitivity reflects the degree of influence of system parameters on the oscillation propagation factor.

Simulation Results and Analysis
In order to verify the validity of the oscillation propagation factor proposed in this paper, the propagation of the subsynchronous oscillation in the test system is analyzed by the oscillation propagation factor.When the interharmonic frequency is changed, and when the system damping ratio is changed, the oscillation propagation factor is calculated to compare the strength of the subsynchronous oscillation of each generator.

Two Area Four Machine System
The modified 4-machine 2-area system [22] with two DFIGs is shown in Figure 2 p.u.The small-signal stability analysis was performed on the system of Figure 2, and seven oscillation modes were obtained, as shown in Table 1.The active power of DFIGs and the input vector are are derived.The real part of the system eigenvalue is negative, and the system is small-signal stable.The interharmonic d and q axis components obtained by Park's transformation [23] are sinusoidal, and the phase of the d-axis component lags behind the phase of the q-axis component π 2 .Therefore, the input ∆u is as shown in (24).
where A i is the amplitude of interharmonic; ω is the angular frequency of interharmonic.Assume interharmonic currents from DFIG1 is with the angular frequency of ω rad/s and amplitude of 0.001 p.u. is injected to the power grid.The output power of each DFIG and synchronous generator is taken as the output variable, and the system transfer function is established.The oscillation propagation factor is shown in Figure 3. Figure 4 shows the oscillation propagation factors for different generators near the oscillation frequency of system mode 1-4.

Δu cos( ) cos(
) where i A is the amplitude of interharmonic; ω is the angular frequency of interharmonic.
Assume interharmonic currents from DFIG1 is with the angular frequency of rad/s and amplitude of 0.001 p.u. is injected to the power grid.The output power of each DFIG and synchronous generator is taken as the output variable, and the system transfer function is established.The oscillation propagation factor is shown in Figure 3. Figure 4 shows the oscillation propagation factors for different generators near the oscillation frequency of system mode 1-4.It can be seen from Figures 3 and 4 that the oscillation propagation factors of the output power of the same generator are different for input variables of different angular frequency.When the angular frequency of the input variable is the same, the oscillation propagation factors of the output power of each generator are also different.When the input variable angular frequency is close to the oscillation angular frequency of system mode 1-4, the oscillation propagation factors of each generator exceed 0 dB, and the interharmonics of the DFIG are amplified by the system propagation.When the input variable frequency is close to the frequencies of mode 3 and 4, the amplitude of the output power oscillation of the DFIG is larger than the amplitude of the synchronous generator output power oscillation.When the input variable frequency is close to the frequencies of mode 1 and 2, the amplitude of the output power oscillation of the DFIG is smaller than the amplitude of the synchronous generator output power oscillation.The oscillation propagation factor of DFIG2 is given in Figure 5.It can be seen from Figures 3 and 4 that the oscillation propagation factors of the output power of the same generator are different for input variables of different angular frequency.When the angular frequency of the input variable is the same, the oscillation propagation factors of the output power of each generator are also different.When the input variable angular frequency is close to the oscillation angular frequency of system mode 1-4, the oscillation propagation factors of each generator exceed 0 dB, and the interharmonics of the DFIG are amplified by the system propagation.When the input variable frequency is close to the frequencies of mode 3 and 4, the amplitude of the output power oscillation of the DFIG is larger than the amplitude of the synchronous generator Energies 2019, 12, 1081 9 of 16 output power oscillation.When the input variable frequency is close to the frequencies of mode 1 and 2, the amplitude of the output power oscillation of the DFIG is smaller than the amplitude of the synchronous generator output power oscillation.The oscillation propagation factor of DFIG2 is given in Figure 5.
angular frequency of the input variable is the same, the oscillation propagation factors of the output power of each generator are also different.When the input variable angular frequency is close to the oscillation angular frequency of system mode 1-4, the oscillation propagation factors of each generator exceed 0 dB, and the interharmonics of the DFIG are amplified by the system propagation.When the input variable frequency is close to the frequencies of mode 3 and 4, the amplitude of the output power oscillation of the DFIG is larger than the amplitude of the synchronous generator output power oscillation.When the input variable frequency is close to the frequencies of mode 1 and 2, the amplitude of the output power oscillation of the DFIG is smaller than the amplitude of the synchronous generator output power oscillation.The oscillation propagation factor of DFIG2 is given in Figure 5.The abscissa of the red dotted line in Figure 5 is the system oscillation angular frequency.The abscissa of the three maximum points of the oscillation propagation factor of DFIG2 is 90.79 rad/s, 89.61 rad/s and 75.12 rad/s, respectively.These three angular frequency are almost equal to the system oscillation angular frequencies of mode 1, 2 and 3 in Table 1, respectively.When the inter-harmonics injected into the power network of the DFIG approaches these frequencies, resonance occurs, and a large amplitude power oscillation takes place.When the amplitude of the fluctuation of the inter-harmonic current becomes larger or smaller, the amplitude of the power oscillation also increases or decreases.The proportional relationship between the two amplitudes does not change, and corresponds to the ordinate value of the oscillation propagation factor characteristic curve at each frequency point.
The effect of the output power of DFIGs on the oscillation propagation factor is analyzed below.When the DFIG1 output power reduces by 25%, the small-signal analysis result is shown in Table 2.The active power of DFIGs and the input vector are P DFIG1 = 0.3p.u., P DFIG2 = 0.4p.u.,u = 0.The sensitivity with system subsynchronous oscillation frequency respect to DFIG active power and the PI controller parameters of DFIGs is shown in Figure 7. From Figure 7, it can be seen that the sensitivity of the oscillation propagation factor to the active power of DFIG2 is positive, indicating that as the active power of DFIG2 decreases, the oscillation propagation factor will decrease.The active power of DFIG2 is more sensitive to mode 4, which demonstrates that regulation of the active power of DFIG2 could block the oscillation propagation path.
When the interharmonic currents from DFIG1 are injected to the power grid, the oscillation propagation factor of DFIG2 output power near the frequency of modes 1-4 is shown in Figure 6.The abscissa of the dotted line in Figure 6 is the system oscillation angular frequency.Comparing Tables 1 and 2 after reducing the DFIG power, the damping ratio of system mode 1 decreases, and the damping ratios of system mode 2 and 3 increase.As can be seen from Figure 6, the oscillation propagation factor at the system angular frequency of mode 1 is increased.The oscillation propagation factor at the system angular frequency of mode 2 and 3 are decreased.The larger the system damping ratio, the smaller the oscillation propagation factor.In order to verify the accuracy of the oscillation propagation factor analysis, the time domain simulations are carried out.The interharmonic currents from DFIG1 are injected at bus 6.The interharmonic current angular frequency is 89.64 rad/s, which is simulated under the condition that the original system and DFIG1 power are reduced respectively.The simulation results are shown in Figure 8.
When the amplitude of the active power oscillation is almost constant, the amplitude of the power oscillation can be obtained from the ordinate of the curve in the figure.The oscillation propagation factor can then be calculated based on (19), the definition of the oscillation propagation factor.In Figure 8, the interharmonic current at the original system bus 6 causes DFIG2 active power oscillation, the oscillation amplitude is 0.0124, and the oscillation propagation factor is 8.77, which is basically the same as the oscillation propagation factor 18.75 / 20 10 8.66 = at 89.61 rad/s in Figure 5.After the power of DFIG1 is reduced, the system damping increases, the oscillation propagation factor decreases, and the power oscillation amplitude decreases.The dynamic frequency domain analysis method based on transfer function can quantitatively describe the proportional relationship between the amplitude of power oscillation and the amplitude of interharmonic current.The abscissa of the red dotted line in Figure 5 is the system oscillation angular frequency.The abscissa of the three maximum points of the oscillation propagation factor of DFIG2 is 90.79 rad/s, 89.61 rad/s and 75.12 rad/s, respectively.These three angular frequency are almost equal to the system oscillation angular frequencies of mode 1, 2 and 3 in Table 1, respectively.When the inter-harmonics injected into the power network of the DFIG approaches these frequencies, resonance occurs, and a large amplitude power oscillation takes place.When the amplitude of the fluctuation of the interharmonic current becomes larger or smaller, the amplitude of the power oscillation also increases or decreases.The proportional relationship between the two amplitudes does not change, and corresponds to the ordinate value of the oscillation propagation factor characteristic curve at each frequency point.
The effect of the output power of DFIGs on the oscillation propagation factor is analyzed below.When the DFIG1 output power reduces by 25%, the small-signal analysis result is shown in Table 2.The active power of DFIGs and the input vector are .The sensitivity with system subsynchronous oscillation frequency respect to DFIG active power and the PI controller parameters of DFIGs is shown in Figure 6.From Figure 6, it can be seen that the sensitivity of the oscillation propagation factor to the active power of DFIG2 is positive, indicating that as the active power of DFIG2 decreases, the oscillation propagation factor will decrease.The active power of DFIG2 is more sensitive to mode 4, which demonstrates that regulation of the active power of DFIG2 could block the oscillation propagation path.When the interharmonic currents from DFIG1 are injected to the power grid, the oscillation propagation factor of DFIG2 output power near the frequency of modes 1-4 is shown in Figure 7.The abscissa of the dotted line in Figure 7 is the system oscillation angular frequency.Comparing Table 1 and Table 2 after reducing the DFIG power, the damping ratio of system mode 1 decreases, and the damping ratios of system mode 2 and 3 increase.As can be seen from Figure 7, the oscillation propagation factor at the system angular frequency of mode 1 is increased.The oscillation propagation factor at the system angular frequency of mode 2 and 3 are decreased.The larger the system damping ratio, the smaller the oscillation propagation factor.

Time Domain Simulation Results
In order to verify the accuracy of the oscillation propagation factor analysis, the time domain simulations are carried out.The interharmonic currents from DFIG1 are injected at bus 6.The interharmonic current angular frequency is 89.64 rad/s, which is simulated under the condition that the original system and DFIG1 power are reduced respectively.The simulation results are shown in Figure 8.
When the amplitude of the active power oscillation is almost constant, the amplitude of the power oscillation can be obtained from the ordinate of the curve in the figure.The oscillation propagation factor can then be calculated based on (19), the definition of the oscillation propagation factor.In Figure 8, the interharmonic current at the original system bus 6 causes DFIG2 active power oscillation, the oscillation amplitude is 0.0124, and the oscillation propagation factor is 8.77, which is basically the same as the oscillation propagation factor 10 18.75/20 = 8.66 at 89.61 rad/s in Figure 5.After the power of DFIG1 is reduced, the system damping increases, the oscillation propagation factor decreases, and the power oscillation amplitude decreases.The dynamic frequency domain analysis method based on transfer function can quantitatively describe the proportional relationship between the amplitude of power oscillation and the amplitude of interharmonic current.

New England 39 Bus System
The New-England 39 bus system [24] with two DFIGs is shown in Figure 9.The parameters of DFIGs in New-England 39 bus system are same as that in 4-machine 2-area system.The small-signal stability analysis was performed on the system shown in Figure 9, and the subsynchronous oscillation modes were obtained as shown in Table 3.The active power of DFIGs and the input vector are  3.

New England 39 Bus System
The New-England 39 bus system [24] with two DFIGs is shown in Figure 9.The parameters of DFIGs in New-England 39 bus system are same as that in 4-machine 2-area system.The small-signal stability analysis was performed on the system shown in Figure 9, and the subsynchronous oscillation modes were obtained as shown in Table 3.The active power of DFIGs and the input vector are P DFIG1 = 0.4p.u., P DFIG2 = 0.4p.u.,u = 0 in Table 3.The New-England 39 bus system [24] with two DFIGs is shown in Figure 9.The parameters of DFIGs in New-England 39 bus system are same as that in 4-machine 2-area system.The small-signal stability analysis was performed on the system shown in Figure 9, and the subsynchronous oscillation modes were obtained as shown in Table 3.The active power of DFIGs and the input vector are  3.

Frequency Domain Simulation Results
The subsynchronous harmonics injected into the grid by DFIG1 are taken as input variables, and the active power of each generator is the output variable.The frequency response of each generator is as shown in Figure 10.The subsynchronous harmonics injected into the grid by DFIG1 are taken as input variables, and the active power of each generator is the output variable.The frequency response of each generator is as shown in Figure 10.It can be seen from Figure 10 that when the interharmonic frequency is close to the oscillation frequency of the system, the oscillation propagation factors of each generator increase.Near the system subsynchronous frequency, the oscillation propagation factors of the two DFIG and part of the generator exceed 0 dB, and the output power is greatly affected by the subsynchronous oscillation caused by the interharmonics.When the input variable frequency is close to the frequencies of mode 3 and 4, the amplitude of the output power oscillation of the DFIG is larger than the amplitude of the synchronous generator output power oscillation.When the input variable frequency is close to the frequencies of mode 1 and 2, the amplitude of the output power oscillation of the DFIG is smaller than the amplitude of the synchronous generator output power oscillation.The oscillation It can be seen from Figure 10 that when the interharmonic frequency is close to the oscillation frequency of the system, the oscillation propagation factors of each generator increase.Near the system subsynchronous frequency, the oscillation propagation factors of the two DFIG and part of the generator exceed 0 dB, and the output power is greatly affected by the subsynchronous oscillation caused by the interharmonics.When the input variable frequency is close to the frequencies of mode 3 and 4, the amplitude of the output power oscillation of the DFIG is larger than the amplitude of the synchronous generator output power oscillation.When the input variable frequency is close to the frequencies of mode 1 and 2, the amplitude of the output power oscillation of the DFIG is smaller than the amplitude of the synchronous generator output power oscillation.The oscillation propagation factor of DFIG2 is given in Figure 11.The abscissa of the red dotted line in Figure 11 is the system oscillation angular frequency.In Figure 11, the two maximum value points of the DFIG2 oscillation propagation factor are 90.77rad/s and 74.91 rad/s, respectively, which are close to the system oscillation angular frequency.When the interharmonic frequency is in the low frequency band, the power oscillation of DFIG2 is hardly caused.When the subsynchronous frequency interharmonics of the DFIG1 are injected into the power network, resonance is induced, and a large amplitude power oscillation of DFIG2 occurs.
The sensitivity with subsynchronous oscillation frequency respect to DFIG active power and the PI controller parameters of DFIGs are shown in Figure 12.Changes in controller parameters can cause changes in the propagation factor, and the amplitude of the oscillation of DFIGs changes.From Figure 12 it can be seen that K is highly sensitive to the oscillation mode 3 and mode 4, which means the regulation of K is effective for reducing the propagation aptitude.For the mode 1 and mode 2, K is not sensitive, and the sensitivity is almost equal to zero, which means the regulation of 1 i K could not help for reducing the propagation aptitude of mode 1 and mode 2. By this way, the parameters could be optimized to satisfy the engineering requirements.The abscissa of the red dotted line in Figure 11 is the system oscillation angular frequency.In Figure 11, the two maximum value points of the DFIG2 oscillation propagation factor are 90.77rad/s and 74.91 rad/s, respectively, which are close to the system oscillation angular frequency.When the interharmonic frequency is in the low frequency band, the power oscillation of DFIG2 is hardly caused.When the subsynchronous frequency interharmonics of the DFIG1 are injected into the power network, resonance is induced, and a large amplitude power oscillation of DFIG2 occurs.
The sensitivity with subsynchronous oscillation frequency respect to DFIG active power and the PI controller parameters of DFIGs are shown in Figure 12.Changes in controller parameters can cause changes in the propagation factor, and the amplitude of the oscillation of DFIGs changes.From Figure 12 it can be seen that K i1 is highly sensitive to the oscillation mode 3 and mode 4, which means the regulation of K i1 is effective for reducing the propagation aptitude.For the mode 1 and mode 2, K i1 is not sensitive, and the sensitivity is almost equal to zero, which means the regulation of K i1 could not help for reducing the propagation aptitude of mode 1 and mode 2. By this way, the parameters could be optimized to satisfy the engineering requirements.

Time Domain Simulation Results
In the MATLAB, a time domain simulation model of New-England 39 bus System with two DFIGs is established, and an inter-harmonic current with an angular frequency of 90.77 rad/s and an amplitude of 0.001 is injected at bus 6.The DFIG2 active power simulation results are shown in Figure 13.The blue curve is the time domain simulation result of DFIG2 active power when the DFIG1 power is increased by a factor of 2. After the power of DFIG1 is increased, the system damping increases, the oscillation propagation factor decreases, and the power oscillation amplitude decreases.In Figure 13, the interharmonic current at bus 6 causes DFIG2 active power oscillation, the oscillation amplitude is 0.0039, and the oscillation propagation factor is 2.76, which is basically the same as the oscillation propagation factor 10 8.866/20 = 2.78 at 90.77 rad/s in Figure 11.12 it can be seen that K is highly sensitive to the oscillation mode 3 and mode 4, which means the regulation of K is effective for reducing the propagation aptitude.For the mode 1 and mode 2, K is not sensitive, and the sensitivity is almost equal to zero, which means the regulation of 1 i K could not help for reducing the propagation aptitude of mode 1 and mode 2. By this way, the parameters could be optimized to satisfy the engineering requirements.In the MATLAB, a time domain simulation model of New-England 39 bus System with two DFIGs is established, and an inter-harmonic current with an angular frequency of 90.77 rad/s and an amplitude of 0.001 is injected at bus 6.The DFIG2 active power simulation results are shown in Figure 13.The blue curve is the time domain simulation result of DFIG2 active power when the DFIG1 power is increased by a factor of 2. After the power of DFIG1 is increased, the system damping increases, the oscillation propagation factor decreases, and the power oscillation amplitude decreases.In Figure 13, the interharmonic current at bus 6 causes DFIG2 active power oscillation, the oscillation amplitude is 0.0039, and the oscillation propagation factor is 2.76, which is basically the same as the oscillation propagation factor 8.866 / 20 10 2.78 = at 90.77 rad/s in Figure 11.

Figure 13
The output active power of DFIG2

Discussion
The evaluation of hazard degree of the subsynchronous oscillation needs to identify the propagation range of the subsynchronous oscillation.In order to maintain the stability of the power system, different methods are required according to different subsynchronous oscillation propagation ranges.The method for dealing with subsynchronous oscillation of power system is: When the oscillation propagation factor of the synchronous generator is larger than a given value, the subsynchronous oscillation propagation range includes the synchronous generators.The subsynchronous oscillation caused by the interharmonics is serious, and the DFIG that causes the subsynchronous harmonics needs to be shed [25].When the synchronous generator is not included in the subsynchronous oscillation propagation range, the subsynchronous oscillation caused by the interharmonics is less harmful.It is not necessary to shed the DFIG that causes the subsynchronous harmonics.We can adjust the operating parameters of the DFIGs or take other measures.
Meanwhile, if a detailed electromagnetic transient model is established for a complete actual power system, it will take a lot of time to determine the subsynchronous oscillation propagation range by time domain simulation analysis.The method can quickly determine the subsynchronous oscillation propagation range by linearizing the transient model and using the oscillation propagation factor.When the DFIGs output is different, the oscillation propagation factor is different, then the oscillation amplitude at a different frequency is different.According to the participation factor [22], only the state variables with the highest participation in the oscillation mode can be obtained, but according to the oscillation propagation factor, the correlation degree between each generator and the oscillation with any frequency can be obtained.Through the sensitivity curve, it can be known that the change of the controller parameters will cause the change of the oscillation

Discussion
The evaluation of hazard degree of the subsynchronous oscillation needs to identify the propagation range of the subsynchronous oscillation.In order to maintain the stability of the power system, different methods are required according to different subsynchronous oscillation propagation ranges.The method for dealing with subsynchronous oscillation of power system is: When the oscillation propagation factor of the synchronous generator is larger than a given value, the subsynchronous oscillation propagation range includes the synchronous generators.The subsynchronous oscillation caused by the interharmonics is serious, and the DFIG that causes the subsynchronous harmonics needs to be shed [25].When the synchronous generator is not included in the subsynchronous oscillation propagation range, the subsynchronous oscillation caused by the interharmonics is less harmful.It is not necessary to shed the DFIG that causes the subsynchronous harmonics.We can adjust the operating parameters of the DFIGs or take other measures.
Meanwhile, if a detailed electromagnetic transient model is established for a complete actual power system, it will take a lot of time to determine the subsynchronous oscillation propagation range by time domain simulation analysis.The method can quickly determine the subsynchronous oscillation propagation range by linearizing the transient model and using the oscillation propagation factor.When the DFIGs output is different, the oscillation propagation factor is different, then the oscillation amplitude at a different frequency is different.According to the participation factor [22], only the state variables with the highest participation in the oscillation mode can be obtained, but according to the oscillation propagation factor, the correlation degree between each generator and the oscillation with any frequency can be obtained.Through the sensitivity curve, it can be known that the change of the controller parameters will cause the change of the oscillation propagation factor, and the amplitude of the oscillation of DFIG2 will change.

Conclusions
The oscillation propagation factor which could characterize the subsynchronous oscillation magnitude propagation, is proposed in this paper.The interharmonics injected into the grid by DFIG will cause continuous oscillation of power.The oscillation propagation factor based on frequency domain method is used to study the propagation characteristics of DFIG interharmonic injected into the grid.The main conclusions are as follows: The oscillation propagation factor can quantify the strength of the subsynchronous oscillation, and could give the propagation intensity to different directions or different devices.The oscillation propagation factor generally has a maximum point at the oscillation frequency of the system, and the amplitude of the power oscillation increases significantly.The sensitivities of the operating parameters and system parameters to the oscillation propagation factor are deduced to give the guidance for blocking the oscillation propagation path.

3 rX
is the filter reactance.C is the capacitance value of the intermediate capacitor, and dc U is the DC voltage of the capacitor.UW is the DFIG bus voltage; IW is the DFIG bus injection current.The positive direction of the physical quantity uses the generator convention.

X
. The DFIG parameters are listed below.Rotor winding resistance 0= p.u. Capacitance C = 13.29 p.u. Inertia time constant 3.4 J T =

Figure 3 .
Figure 3. Oscillation propagation factor of each generator.

Figure 4 .
Figure 4. Oscillation propagation factor near the oscillation frequency.

Figure 4 .
Figure 4. Oscillation propagation factor near the oscillation frequency.

Figure 8 .
Figure 8.The output active power of DFIG2.

Figure 10 .
Figure 10.Oscillation propagation factor of each generator

Figure 10 .
Figure 10.Oscillation propagation factor of each generator.

Figure 12
Figure 12The sensitivity of DFIG2 oscillation propagation factor.

Figure 13 .
Figure 13.The output active power of DFIG2.