Modelling and Stability Analysis of Wind Power Plants Connected to Weak Grids

Featured Application: A uniﬁed frequency-domain modelling and stability analysis is presented to predict the impact of mutual interactions between a wind power plant and a weak grid on the stability margin. Abstract: It is important to develop modelling tools to predict unstable situations resulting from the interactions between the wind power plant and the weak power system. This paper presents a uniﬁed methodology to model and analyse a wind power plant connected to weak grids in the frequency-domain by considering the dynamics of the phase lock loop (PLL) and controller delays, which have been neglected in most of the previous research into modelling of wind power plants to simplify modelling. The presented approach combines both dq and positive / negative sequence domain modelling, where a single wind turbine is modelled in the dq domain but the whole wind power plant connected to the weak grid is analysed in the positive / negative sequence domain. As the proposed modelling of the wind power plant is systematic and modular and based on the decoupled positive / negative sequence impedances, the application of the proposed methodology is relevant for transmission system operators (TSOs) to assess stability easily with a very low compactional burden. In addition, as the analytical dq impedance models of the single wind turbine are provided, the proposed methodology is an optimization design tool permitting wind turbine manufacturers to tune their converter control. As a case study, a 108 MW wind power plant connected to a weak grid was used to study its sensitivity to variations in network short-circuit level, X / R ratio and line series capacitor compensation (Xc / Xg).


Introduction
With the increasing number of solar and wind power plants in the electrical grid, the installation ratio of traditional synchronous generators to the total production capacity is significantly decreasing, which is leading to weak transmission grids [1]. Worldwide, the transmission system operators are in the process of updating their grid connection requirements for generation units to ensure that transmission networks can continue to be operated in a secure manner when renewable energy resources replace traditional thermal generation units. It is important to have simulation models, which can reliably predict the stability margins under contemporary operational conditions. The grid-stability of such systems can be improved in multiple ways through combinations of properties of the generation units and an additional balance of plant equipment such as STATCOMS, battery storage, etc. is used to analyse the wind power plant stability based on the system poles and based on Nyquist criterion of the open-loop system. The presented method is scalable and adaptable to different scales of the power system along with a low computational burden.
In the following section, a wind turbine with full-scale converter configuration is considered to be modelled. The dynamics of the PLL and the time delay, which have been usually neglected in the literature on wind power plant stability analysis [24], are also included to have a more accurate model. Then, a 108 MW wind power plant connected to a weak grid is analysed under different conditions.

Full-Scale Converter-Based Wind Turbines
The full-scale converter configuration, which is shown in Figure 1, is one of the popular configurations for variable speed wind turbines. In this configuration, the grid side converter usually controls the dc-link voltage, and the machine side converter controls the active power [26]. If the small oscillations of the dc-link voltage are neglected, the dynamics of the wind turbine around and above the fundamental frequency can be predicted by modelling of the grid side converter [27]. Therefore, the dynamics of the machine side converter can be neglected in this paper.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 15 system. The presented method is scalable and adaptable to different scales of the power system along with a low computational burden. In the following section, a wind turbine with full-scale converter configuration is considered to be modelled. The dynamics of the PLL and the time delay, which have been usually neglected in the literature on wind power plant stability analysis [24], are also included to have a more accurate model. Then, a 108 MW wind power plant connected to a weak grid is analysed under different conditions.

Full-Scale Converter-Based Wind Turbines
The full-scale converter configuration, which is shown in Figure 1, is one of the popular configurations for variable speed wind turbines. In this configuration, the grid side converter usually controls the dc-link voltage, and the machine side converter controls the active power [26]. If the small oscillations of the dc-link voltage are neglected, the dynamics of the wind turbine around and above the fundamental frequency can be predicted by modelling of the grid side converter [27]. Therefore, the dynamics of the machine side converter can be neglected in this paper.  Figure 2 shows the grid side converter of a wind turbine with its control in the dq frame. The grid side converter is a voltage source converter, and its control consists of four main parts: grid synchronization, dc-link voltage controller, current controller and pulse width modulation (PWM) [28]. The grid synchronization loop estimates the phase angle of the point-of-connection (PoC) voltage, and this angle is used to transform the signals between the abc frame and the dq frame [29]. The current references (Id ref and Iq ref ) are generated from the dc-link voltage and reactive power controller. The error of these references is compared to the measured grid currents (Id and Iq) and fed to proportional integral (PI) controllers (Gcc).    Figure 2 shows the grid side converter of a wind turbine with its control in the dq frame. The grid side converter is a voltage source converter, and its control consists of four main parts: grid synchronization, dc-link voltage controller, current controller and pulse width modulation (PWM) [28]. The grid synchronization loop estimates the phase angle of the point-of-connection (PoC) voltage, and this angle is used to transform the signals between the abc frame and the dq frame [29]. The current references (I d ref and I q ref ) are generated from the dc-link voltage and reactive power controller. The error of these references is compared to the measured grid currents (I d and I q ) and fed to proportional integral (PI) controllers (G cc ).
The outputs of the current controllers are fed to the PWM modulator to generate the final switching signals of the insulated-gate bipolar transistors (IGBTs). The time-delay of the PWM and the discrete control can be modelled by using the Padé approximation [30]. synchronization, dc-link voltage controller, current controller and pulse width modulation (PWM) [28]. The grid synchronization loop estimates the phase angle of the point-of-connection (PoC) voltage, and this angle is used to transform the signals between the abc frame and the dq frame [29]. The current references (Id ref and Iq ref ) are generated from the dc-link voltage and reactive power controller. The error of these references is compared to the measured grid currents (Id and Iq) and fed to proportional integral (PI) controllers (Gcc).  Figure 3 shows the model of the current controller, the feedforward terms, the PWM modulator and the filter choke of the grid side converter in the dq frame. The measurement delay is very small compared to the PWM delay and can be neglected. An outer dc-link voltage control loop is implemented in the GSC and is responsible to keep the dc-link voltage constant. The bandwidth of this outer control loop is very low compared to the fundamental frequency. This paper focuses on the oscillations around and above the fundamental frequency. Therefore, such low-bandwidth outer loops can be neglected [13]. As it is assumed that the dc-link voltage is constant, the dynamics of the dc-link controller have been neglected. Because of the dynamics of the grid synchronization loop, the dq domain model comprises a fixed global frame and a local controller frame. In the model, the measurement delay is neglected, but the PWM delay (G d ) is included as defined in [30]. The grid voltage is added to the output current controller (G cc ) through a high-pass/low-pass filter (G f ). The outputs of the current controllers are fed to the PWM modulator to generate the final switching signals of the insulated-gate bipolar transistors (IGBTs). The time-delay of the PWM and the discrete control can be modelled by using the Padé approximation [30]. Figure 3 shows the model of the current controller, the feedforward terms, the PWM modulator and the filter choke of the grid side converter in the dq frame. The measurement delay is very small compared to the PWM delay and can be neglected. An outer dc-link voltage control loop is implemented in the GSC and is responsible to keep the dc-link voltage constant. The bandwidth of this outer control loop is very low compared to the fundamental frequency. This paper focuses on the oscillations around and above the fundamental frequency. Therefore, such low-bandwidth outer loops can be neglected [13]. As it is assumed that the dc-link voltage is constant, the dynamics of the dc-link controller have been neglected. Because of the dynamics of the grid synchronization loop, the dq domain model comprises a fixed global frame and a local controller frame. In the model, the measurement delay is neglected, but the PWM delay (Gd) is included as defined in [30]. The grid voltage is added to the output current controller (Gcc) through a high-pass/low-pass filter (Gf).

Small-Signal Admittance Model in dq-Domain
Current controller Plant PWM and Delay  Figure 3. Dq small signal model of the current controller, PLL and the filter inductance (Lf). According to Figure 3, the small signal output currents and dq output voltage can be obtained by the following matrices: According to Figure 3, the small signal output currents and dq output voltage can be obtained by the following matrices: By neglecting the dynamics of the grid synchronization loop, e.g., PLL, the dq frame of the plant and the dq frame of the controller would be the same, i.e., By substituting (2) and (3) in (1), the admittance matrix can be obtained as where the dynamics of the grid synchronization loop are ignored.
For the grid synchronization loop, the synchronous reference frame phase-locked loop (SRF-PLL) can be used and the small signal model is shown in Figure 4.
where the dynamics of the grid synchronization loop are ignored.
For the grid synchronization loop, the synchronous reference frame phase-locked loop (SRF-PLL) can be used and the small signal model is shown in Figure 4. When small-signal perturbations are added to the Vg voltage, the estimated angle will have a small error (Δθ) between the local (controller) and global (plant) frames because of the PLL dynamics. In Figure 4b, fdq c corresponds to the estimated angle and fdq g corresponds to the actual angle. Based on Figure 4, the voltage and current in the controller frame and plant frame under small-signal perturbations can be derived as Using (1), (2), (5) and (6) and some simplifications, the small signal admittance matrix can be obtained by When small-signal perturbations are added to the V g voltage, the estimated angle will have a small error (∆θ) between the local (controller) and global (plant) frames because of the PLL dynamics.
In Figure 4b, f dq c corresponds to the estimated angle and f dq g corresponds to the actual angle. Based on Figure 4, the voltage and current in the controller frame and plant frame under small-signal perturbations can be derived as Using (1), (2), (5) and (6) and some simplifications, the small signal admittance matrix can be obtained by

Equivalent Small-Signal Admittance Model in Sequence Domain
In the previous section, the small-signal admittance model of the converter in the dq domain has been obtained. However, the admittance model can also be obtained in the positive-negative sequence domain [18,20]. The small-signal admittance matrix in the sequence domain can also be written by a 2 × 2 admittance matrix, where its off-diagonal elements show a coupling between the positive and negative components. This coupling is very small and can be neglected [7]. The relationship between the sequence domain admittance definition and dq domain admittance definition [20] can be obtained by the following equation as derived in Appendix A.
The detailed form of Equation (8) is obtained by substituting Equation (7). In this paper, the dq domain admittance is obtained, and then the equivalent sequence admittance is derived based on (8). For the linear stability analysis of the wind power plant in the sequence domain, the coupling terms may be neglected, and the conventional positive impedance is calculated as:

Total Wind Turbine Thévenin Impedance
As shown in Figure 5, there is a filter capacitor and a transformer in the output circuit of the wind turbine. Therefore, the total Thévenin impedance of the wind turbine can be equated to where Z p is the equivalent small-signal positive impedance of the wind turbine, Z cf is the filter capacitor impedance, Z lg is the leakage impedance transformed in the primary side, and K is the transformer voltage ratio. As can be seen, Z p is in parallel with Z cf and they are in series with Z lg . Then this equivalent impedance is transferred to the high-voltage side by the transformer ratio to the power of two (K 2 ).
where Zp is the equivalent small-signal positive impedance of the wind turbine, Zcf is the filter capacitor impedance, Zlg is the leakage impedance transformed in the primary side, and K is the transformer voltage ratio. As can be seen, Zp is in parallel with Zcf and they are in series with Zlg. Then this equivalent impedance is transferred to the high-voltage side by the transformer ratio to the power of two (K 2 ).

Case Study of a Wind Power Plant Connected to a Weak Grid
In order to perform an impedance-based stability analysis, an aggregated model of the wind power plant is required. For example, Figure 6 shows a 108 MW wind power plant with three strings and ten 3.6 MW wind turbines installed on each string. All wind turbines in the wind power plant are assumed to have the same structure. Therefore, the impedance equivalent of the aggregated wind turbines can be obtained based on per-unit values.

Case Study of a Wind Power Plant Connected to a Weak Grid
In order to perform an impedance-based stability analysis, an aggregated model of the wind power plant is required. For example, Figure 6 shows a 108 MW wind power plant with three strings and ten 3.6 MW wind turbines installed on each string. All wind turbines in the wind power plant are assumed to have the same structure. Therefore, the impedance equivalent of the aggregated wind turbines can be obtained based on per-unit values. If it is assumed that the voltages and currents of the wind turbines are the same, ten 3.6 MW wind turbines of each feeder can be aggregated into one 36 MW wind turbine (see Figure 7), where its equivalent series cable parameter can be obtained [31,32]   If it is assumed that the voltages and currents of the wind turbines are the same, ten 3.6 MW wind turbines of each feeder can be aggregated into one 36 MW wind turbine (see Figure 7), where its equivalent series cable parameter can be obtained [31,32] by  ZAT is the equivalent series impedance, and BAT is the equivalent shunt susceptance, respectively. Finally, three 36 MW parallel aggregated wind turbines can be aggregated as one 108 MW wind turbine (shown in Figure 7), where its equivalent series parameter can be calculated by

Stability Analysis
By obtaining the equivalent wind power plant impedance and the equivalent grid impedance, Figure 8 can be used for stability analysis. Based on this figure, the transfer function for the current can be obtained as For the stability analysis, the poles of the system can be obtained by solving the following equation: As the wind power plant is stable when connected to a strong grid, i.e., Vwt/Zwt has no right-half plane poles, an alternative option for the stability analysis is to assess the Nyquist stability criterion for Zgrid/Zwt. A weak grid model is considered here including inductive reactance (Xg), resistance (R), and series capacitive reactance (Xc) for compensation. To illustrate sensitivity to variation in the network, short-circuit level, Xg/R ratio and line series compensation, three scenarios, are studied, all based on a fictitious wind turbine converter controller, the parameters of which are listed in [21]. The wind turbine has been tuned for a strong grid. In Appendix B, all the variables and parameters used in the figures are summarized.  Z AT is the equivalent series impedance, and B AT is the equivalent shunt susceptance, respectively. Finally, three 36 MW parallel aggregated wind turbines can be aggregated as one 108 MW wind turbine (shown in Figure 7), where its equivalent series parameter can be calculated by

Stability Analysis
By obtaining the equivalent wind power plant impedance and the equivalent grid impedance, Figure 8 can be used for stability analysis. Based on this figure, the transfer function for the current can be obtained as For the stability analysis, the poles of the system can be obtained by solving the following equation: As the wind power plant is stable when connected to a strong grid, i.e., V wt /Z wt has no right-half plane poles, an alternative option for the stability analysis is to assess the Nyquist stability criterion for Z grid /Z wt . A weak grid model is considered here including inductive reactance (X g ), resistance (R), and series capacitive reactance (X c ) for compensation. To illustrate sensitivity to variation in the network, short-circuit level, X g /R ratio and line series compensation, three scenarios, are studied, all based on a fictitious wind turbine converter controller, the parameters of which are listed in [21]. The wind turbine has been tuned for a strong grid. In Appendix B, all the variables and parameters used in the figures are summarized. 6.1.1. Short Circuit Ratio (SCR) Variations for X g /R = 20 and X c /X g = 0 In this scenario, the SCR is changed from 1.25 to 1.5, X g /R = 20 and X c /X g = 0 (without series compensation). Some poles of the system are shown in Table 1. As it can be seen, there is an unstable pole for SCR = 1.25. Table 1. Poles for short circuit ratio (SCR) variations for X g /R = 20 and X c /X g = 0. Nyquist criterion-based stability analyses for SCR = 1.25 and SCR = 1.5 are shown in Figure 9. When the SCR = 1.25, there is an encirclement of the point −1, which shows that the system is unstable. Nyquist criterion-based stability analyses for SCR = 1.25 and SCR = 1.5 are shown in Figure 9. When the SCR = 1.25, there is an encirclement of the point −1, which shows that the system is unstable. In this scenario, Xg/R is decreased from 20 to 5 and to 2.5, when the SCR = 1.25. As it can be seen from the pole analysis in Table 2 and the Nyquist analysis in Figure 10, the system becomes more stable by decreasing the Xg/R ratio.  6.1.2. X g /R variations for SCR = 1.25 and X c /X g = 0

SCR
In this scenario, X g /R is decreased from 20 to 5 and to 2.5, when the SCR = 1.25. As it can be seen from the pole analysis in Table 2 and the Nyquist analysis in Figure 10, the system becomes more stable by decreasing the X g /R ratio.   Table 3 and Figure 11 show the corresponding poles and Nyquist criteria under Xc/Xg variations for SCR = 1.25 and Xg/R = 20. It can be seen that the real part of the pole around 70 Hz becomes more negative (more stable) along with more compensation. 6.1.3. X c /X g Variations for SCR = 1.25 and X g /R = 20 Table 3 and Figure 11 show the corresponding poles and Nyquist criteria under X c /X g variations for SCR = 1.25 and X g /R = 20. It can be seen that the real part of the pole around 70 Hz becomes more negative (more stable) along with more compensation. Table 3. Poles for X c /X g variations for SCR = 1.25 and X g /R = 20.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 10 of 15  In order to validate the frequency-domain analysis, a time-domain simulation was performed for SCR = 1.25, Xg/R = 20 and Xc/Xg = 0. The current injected to the grid and its fast Fourier transform (FFT) are shown in Figure 12. As it can be seen, some oscillations around 73 Hz (459 rad) are increased until the wind turbine goes to trip. This result is matched with the frequency-domain result in Table  1, which indicates there is one unstable mode as 4.61 + j449. In order to validate the frequency-domain analysis, a time-domain simulation was performed for SCR = 1.25, Xg/R = 20 and X c /X g = 0. The current injected to the grid and its fast Fourier transform (FFT) are shown in Figure 12. As it can be seen, some oscillations around 73 Hz (459 rad) are increased until the wind turbine goes to trip. This result is matched with the frequency-domain result in Table 1, which indicates there is one unstable mode as 4.61 + j449.
for SCR = 1.25, Xg/R = 20 and Xc/Xg = 0. The current injected to the grid and its fast Fourier transform (FFT) are shown in Figure 12. As it can be seen, some oscillations around 73 Hz (459 rad) are increased until the wind turbine goes to trip. This result is matched with the frequency-domain result in Table  1, which indicates there is one unstable mode as 4.61 + j449.

Conclusions
This paper presents a unified impedance-based modelling and analysis tool to predict unstable conditions resulting from the interactions between the wind power plant and the weak grid. The proposed approach allows to have advantages of both the dq domain modelling and sequence domain modelling. This enables transmission system operators (TSOs) to do stability analyses of large power systems more efficiently, and it enables wind turbine manufacturers to optimise the controller parameters to increase the stability margin under a weak grid connection. A 108 MW wind power plant has been studied and the short circuit ratio (SCR) of the grid is considered to be 1.25 and 1.5. The poles of the closed-loop system in the frequency-domain are calculated to find the frequency and damping of the oscillatory modes. Nyquist criteria of the open-loop system are plotted to predict easily the stability margin of the wind power plant connected to the weak grid. The results of a simulated wind power plant have shown that a wind power plant which is tuned for a strong grid can become unstable when connected to a weak grid. The proposed method is a valuable tool when adapting and retuning the control system to actual grid conditions. Author Contributions: E.E. and F.B. proposed the initial methodology, T.L. proposed how to apply the method for weak grids, J.G.N. provided the data, and P.C.K reviewed and edited the paper.

Conclusions
This paper presents a unified impedance-based modelling and analysis tool to predict unstable conditions resulting from the interactions between the wind power plant and the weak grid. The proposed approach allows to have advantages of both the dq domain modelling and sequence domain modelling. This enables transmission system operators (TSOs) to do stability analyses of large power systems more efficiently, and it enables wind turbine manufacturers to optimise the controller parameters to increase the stability margin under a weak grid connection. A 108 MW wind power plant has been studied and the short circuit ratio (SCR) of the grid is considered to be 1.25 and 1.5. The poles of the closed-loop system in the frequency-domain are calculated to find the frequency and damping of the oscillatory modes. Nyquist criteria of the open-loop system are plotted to predict easily the stability margin of the wind power plant connected to the weak grid. The results of a simulated wind power plant have shown that a wind power plant which is tuned for a strong grid can become unstable when connected to a weak grid. The proposed method is a valuable tool when adapting and retuning the control system to actual grid conditions. Author Contributions: E.E. and F.B. proposed the initial methodology, T.L. proposed how to apply the method for weak grids, J.G.N. provided the data, and P.C.K. reviewed and edited the paper.
Funding: This research was funded by "Vestas Wind System A/S" and "Innovation Fund Denmark, grant number 8054-00023B".

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
A Relationship between dq-domain admittance/impedance and sequence-domain admittance/impedance. In this part, the relationship between dq-domain admittance/impedance and sequence-domain admittance/impedance is derived.
∆i d ∆i q = 1 1 −j j ∆i p ∆i n (A5) As already presented in this paper, the small-signal behaviour of a three-phase converter can be modelled in dq domain by ∆i d ∆i q = y dd (s) y dq (s) y qd (s) y qq (s) Based on (A8), the relationship between sequence-domain admittance definition and dq-domain admittance definition can be obtained by y pp (s) y pn (s) y np (s) y nn (s) = 1 1 −j j −1 y dd (s) y dq (s) y qd (s) y qq (s)

Appendix B
In the following table, all the variables and parameters used in the figures are defined.