Analysis of subsynchronous current propagation path of subsynchronous oscillation induced by renewable energy integrated to the power grid

: This paper deals with operational parameter to study the propagation path of the subsynchronous current generated by the wind farm. First, the complex frequency domain operation parameter model of the transmission network, thermal power generator and other power system components are established. Then, the wind farm is taken as the subsynchronous current injection point and the subsynchronous current ’ s propagation path is analysed according to the frequency impedance characteristic of the transmission network. The in ﬂ uences of the current frequency, the position of the wind farms and the structure of the AC system on the harmonic current distribution are analysed and studied with this method. Then the digital simulation method is used to demonstrate the conclusion drawn from the theoretical method. On the basis of the method proposed, the subsynchronous component shunt coef ﬁ cient of the power system actual operation can be obtained initially.


Introduction
Today, centralised and grid-connected wind farms as well as photovoltaic (PV) power stations have drawn much attention in Chinese renewable energy utilisation. However, concentrated interconnection to the weak power grid of large-scale wind farms or PV power plants may contribute to severe risk of inducing subsynchronous oscillation (SSO) problem, which is a hindrance to the development of Chinese renewable energy accommodation. In recent years, a large-scale wind farm in the northwest of China faced frequent SSO problems. The wind farm produces a large number of subsynchronous and supersynchronous frequency current, which slides into electricity delivery channels through the wind farm and then widely spreads to the transmission network [1]. In further critical conditions, subsynchronous frequency components may induce large thermal power plants' torsional vibration protection to trip, which are at a distance of more than 100 km away from the wind farm. As a consequence, analysis of the propagation path of the subsynchronous current produced by wind farms has great significance for monitoring and analysing the risk of an SSO. This paper studied the use of operational parameter to obtain the subsynchronous current shunt coefficient, and then study the propagation path. The authors in references [2][3][4][5] focused on the distribution of harmonic current in shunt coefficient network, and obtained results based on three-phase stationary coordinate system. During the study of the SSO, the generator has to be considered, but the generator impedance is a time-dependent parameter in stationary coordinate [6], so the method presented in references [2][3][4][5] cannot accurately analyse the shunt coefficient of the subsynchronous current. The authors in references [7][8][9][10][11][12][13] show the effect of injecting mode current on the electrical damping coefficient, which is analysed qualitatively, but the shunt coefficient of the injection current is not given. The relevant idea can provide some reference for calculating the shunt coefficient. Peng et al. [14] qualitatively analyse the complementary subsynchronous and supersynchronous frequency components in synchronously rotating reference frame (SRF), which are the same components and have the same shunt coefficients, but the authors do not propose a method to analyse the shunt coefficient of the analytic calculation. Peng et al. [15] verified the shunt coefficient of the injection current under torsional modal frequencies in single-machine infinite-bus power system, but the shunt coefficient of subsynchronous current is not a complete analysis and the shunt coefficient in multi-machine system is not analysed.
In this paper, the complex frequency domain operation parameter model of the power system is established to analyse SSO. The point of common coupling (PCC) of renewable energy integrated to the power grid is assumed as the subsynchronous current injection point. The propagation path of the subsynchronous current in the transmission network is analysed by using the shunt coefficient. On the basis of the method proposed in this paper, the subsynchronous component shunt coefficient of the power system actual operation can be obtained initially.
The rest of this paper is organised as follows. Section 2 discusses the subsynchronous frequency current shunt coefficient in singlemachine infinite-bus power system. Section 3 deals with the subsynchronous frequency current shunt coefficient in multi-machine infinite-bus power system. Section 4 explains digital simulation method used to verify the conclusion drawn from the theoretical method. Section 5 analyses the subsynchronous component shunt coefficient of the actual power system in the northwest of China and verifies the actual data from the phasor measurement unit (PMU).
2 Theoretical analysis method of subsynchronous current shunt coefficient in single-machine infinite-bus power system

Power system component operational impedance in SRF
When an SSO occurs, the frequency components are neither fundamental frequency (50 Hz or 60 Hz) nor low frequency (<10 Hz), the power system component operational impedance is a timedependent parameter in the stationary coordinate. Therefore, the electromagnetic transient process of the power system is required to be considered when establishing a mathematical model of operational impedance for the analysis of SSOs [15].
2.1.1 Operation parameters of synchronous motor: The concept of operational parameters is normally used in the analytical and frequency domain analysis of the synchronous motor transition process. Shijie et al. [ 16] based on the d-q-axis equivalent circuit get its flux linkage equation and voltage equation, then derive the synchronous motor operation parameters. The conclusions are briefly reviewed in the following part.
Generator operation reactance at q-axis equivalent circuit: Generator operation reactance at d-axis equivalent circuit: As the measured method to obtain the generator operation reactance is more complex and difficult to obtain. Generator operation impedance at d-q-axis can be obtained from formulae (1) and (2) by using nameplate values such as X aq , X ad , X D and X Q1

Operation parameters of transmission line:
When the transmission line is represented by a concentrated impedance parameter, the effect of the ground capacitance is ignored. In reference [6], the transmission line's electromagnetic transient model is obtained by Park transformation under the synchronous coordinates. The voltage equation is as follows: The magnetic chain equation is as follows: where c d is the flux linkage of the transmission line at the d-axis; c q the flux linkage of the transmission line at the q-axis; u d the voltage of the transmission line at the d-axis; u q the voltage of the transmission line at the q-axis; i d the current of the transmission line at the d-axis; i q the current of the transmission line at the q-axis; and r l is the resistance of the transmission line. As the transmission line is a stationary element, its d-q-axis parameters are equal, so the operation reactance of the transmission line at d-q-axis is as follows: The operation parameters of the transmission line can be obtained from formulae (4)-(6) as shown in (7) 2.2 Subsynchronous current shunt coefficient in single-machine infinite-bus power system The shunt coefficient of the subsynchronous current is analysed and calculated by using the first benchmark model (line complement degree is 0.74) for computer simulation of subsynchronous resonance [17]. The renewable energy is taken as the subsynchronous current injection point, and the influence of different positions of PCC is to be considered. The subsynchronous current is injected in different PCCs as shown in Fig. 1, such as the main transformer high-pressure side (A) and the end of the line (B). Then, the shunt coefficient at the generator side can be calculated by calculating the ratio of the subsynchronous current amplitude at the generator side with the injection current amplitude at the PCC. The magnitude and phase change of the injection current to the generator depends entirely on the frequency impedance characteristics of the network [14]. Taking the subsynchronous frequency component at point A as an example, the generator d-axis shunt coefficient is obtained. Considering the subsynchronous injection current at point A as an example, the generator d-axis shunt coefficient is as follows: where Z L p , Z T p and Z s p , respectively, represent the operational reactance of the transmission line, transformer and equivalent inductance of infinite-power supply, Z Gd p is the operational reactance of the generator d-axis.
The generator q-axis shunt coefficient is as follows: where Z Gq p is the operational reactance of the generator q-axis.
Assuming the frequency of injection current in a three-phase stationary coordinate as w i , converted to SRF the frequency then is (w 0 − w i ). Substituting p = j(w 0 − w i ) into formulae (8) and (9), the shunt coefficient can be derived.
The shunt coefficient can be obtained in a similar method when the injection current is at point B. Finally, the variation curves of the current shunt coefficient with frequency under different PCCs are shown in Fig. 2.
In Fig. 2, injecting different frequency currents shows different shunt coefficients in the same grid structure; and injecting subsynchronous current at different PCCs show different shunt coefficients in the same grid structure.  3 Theoretical analysis method of subsynchronous current shunt coefficient in multi-machine infinite-bus power system

Unified SRF
The two machines infinite-bus power system is shown in Fig. 3. When the output level of thermal power units, like G 1 , G 2 , is not inconsistent or the outlet parameters of different thermal power plants, like R L1 , X L1 , R L2 and X L2 , are not approximately the same, the SRF of G 1 may be different from the SRF of G 2 . In order to accurately analyse the diversion law of subsynchronous frequency components generated by the renewable energy, the SRF should be unified in the two machines system and multi-machine system.
Zhenghang et al. [18] introduce four different reference frames of multi-machine power system: absolute angle synchronised reference frame, relative angle synchronised reference frame, reference frame in which some machine is taken as a reference machine, and inertia centre reference frame. As the analysis of singlemachine infinite-system's subsynchronous current shunt coefficient in Section 2.2 is based on SRF of the thermal motor unit, so the analysis of two-machine system or multi-machine system can use the relative angle synchronised reference frame. Regarding the SRF of the generator G 1 as a reference, the SRF of any other generator is converted to the reference frame.
As shown in Fig. 4, the angle between the d-axis of the synchronous frame of the generator G 1 and the a-axis of the three-phase stationary frame is u 1 , the angle between the d-axis of synchronous frame of the generator G 2 and the a-axis of three-phase stationary frame is u 2 , the angle between the d-axis of G 1 and G 2 is Du.
The generator terminal current can be represented by (10), and based on Park transformation, the current can be represented by (11) Substitute v 1 = v into formula (11) When Du = u 1 − u 2 , the relationship of d-q-axis of the terminal current of G 2 between the SRF of G 2 and the SRF of G 1 can be shown in the below equation Assume The relationship of d-q-axis of the terminal current of G 2 between the SRF of G 2 and the SRF of G 1 can be obtained by a similar method as shown in the following equation: From the relation between voltage and current (19), the relationship of d-q-axis of operational reactance of G 2 between the SRF of G 2 and the SRF of G 1 can be shown in (20)

Subsynchronous current shunt coefficient in two machines infinite-bus power system
After the operational reactance of two-machine systems and multimachine system is obtained under unified SRF, combine the calculation method of the shunt coefficient proposed in Section 2.2, the shunt coefficient of the system shown in Fig. 3 can be calculated. Considering subsynchronous injection current at the bus bar of G 1 and G 2 as an example, the shunt coefficient of G 1 is shown as follows: (see (21)) (see (22)) Substitute p = j(w 0 − w i ) into formulae (21) and (22), the d-q-axis current shunt coefficient of G 1 at the unified SRF can be derived.
The shunt coefficient when the injection current is at different PCCs can be obtained in a similar manner. By combining the calculated results K d1 and K q1 , the shunt coefficient of generators G 1 and G 2 in a three-phase stationary frame can be derived by inverse Park transformation. Set the parameters of the system shown in Fig. 2 according to Table 1 and set the values of R L1 and X L1 to 0, the angle Du between the d-axis of G 1 and G 2 can be calculated through power flow calculation. Substituting Du into formulae (17)-(25), the variation curves of G 1 and G 2 current shunt coefficient with frequency in different operating modes are shown in Figs. 5 and 6.
Comparing Figs. 5 and 6 shows that the shunt coefficient is highly influenced by the system component parameters and the change of the grid structure; the shunt coefficient is not so much influenced by the generator output level, and different generators have different shunt coefficients.

Simulation analysis
To verify the correctness of the above theoretical analysis, the simulation model was established in PSCAD (power systems computeraided design), and the subsynchronous current from 1 to 50 Hz (fundamental frequency of this model is 60 Hz) was injected at different positions, and the corresponding shunt curve was obtained and compared with the theoretical calculation result.

Comparison of single-machine infinite-bus power system analytical results
The simulation model as shown in Fig. 1 is established in PSCAD, the shunt coefficient curves obtained from simulation are compared with the curves obtained from theoretical calculations mentioned in Section 2.2. The compared graphs based on the subsynchronous current injection position are shown in Figs. 7 and 8.
In Figs. 7 and 8, the trend of PSCAD simulation curves of the machine-side current shunt coefficient in d-q-axis and the theoretical analysis curve is consistent, only a small part does not match. Also, the conclusions derived from Section 2.2 are validated: different frequency injection current shows different shunt coefficients in the same grid structure; and subsynchronous injection current at different PCCs shows different shunt coefficients in the same grid structure.

Comparison of two-machine infinite-bus power system analytical results
The simulation model as shown in Fig. 3 is established in PSCAD, the operating mode is adjusted according to the component parameters and generator output level shown in Table 1. The shunt coefficient curves obtained from simulation are compared with the curves obtained from theoretical calculation in Section 3.2. The compared graphs according to different operating modes are shown in Figs. 9 and 10.
In Figs. 9 and 10, the trend of PSCAD simulation curves of the G 1 and G 2 current shunt coefficient and the theoretical analysis curve is consistent, only a small part does not match. Also, the conclusions derived from Section 3.2 are validated: the shunt coefficient is highly influenced by system component parameters and the change of the grid structure; the shunt coefficient is little influenced by the generator output level, and different generators have different shunt coefficients.   The simplified diagram of the large thermal power units and the system nearby where the PMU data is recorded from is shown in Fig. 11.
When the number of grid-connected generators is different, the method mentioned in this paper is used to calculate the shunt coefficient, and combined with the analysis results of PMU recording data, the contrast analysis diagram is shown in Fig. 12.
In Fig. 12, the number of grid-connected generators has a greater impact on the shunt coefficient. PMU recording data analysis results are consistent with the theoretical calculation trends, and there is a little difference. The possible reason is that the theoretical equivalence of the network at the time of theoretical calculations is somewhat simplistic, but the actual system is much more complex.

Conclusion
This paper used operational parameter to get the subsynchronous current diversion rule, and then study the propagation path. The study shows that the different frequency injection current exhibits different shunt coefficients in the same grid structure, and injection subsynchronous current at different PCCs show different shunt coefficients in the same grid structure. Then the analysis method is extended from single-machine infinite-bus power system to multi-machine infinite-bus power system. The conclusions derived from theoretical calculation validated by simulation that the shunt coefficients are highly influenced by system component parameters and the change of the grid structure, and the shunt coefficients are not so much influenced by the generator output level, and different generators have different shunt coefficients. Finally, the analysis of the propagation path about subsynchronous current produced by wind farms has great significance for monitoring and analysing the risk of SSO.