Modeling and Stability Analysis of Weak-Grid Tied Multi-DFIGs in DC-Link Voltage Control Timescale

: The DC-link voltage control (DVC) timescale (i.e., the frequency dynamics covering converter outer controls) instabilities in wind generation have gained increased attention recently. This paper presents DVC timescale modeling and stability analysis for multi doubly-fed induction generators (DFIGs) connected to weak AC grids. A reduced-order, small-signal model of a grid-tied multi-DFIG system, designed for DVC dynamics analysis, is ﬁrstly proposed. The model allows for the dynamic interactions among the DC-link voltage control, active power control (APC), terminal voltage control (TVC) and phase-locked loop (PLL). Eigenvalue and participation factor analyses are conducted to explore the potential instabilities and correlated critical factors for such a multi-machine system. The sensitivity studies ﬁnd that instability can occur at high levels of power generations or low short-circuit ratio (SCR) conditions. In addition, the dominant mode is identiﬁed to be highly related to the PLL, and its modal damping is decreased when the bandwidths of PLLs in di ﬀ erent generators are close. Detailed model-based time domain simulations veriﬁed the analysis above.


Introduction
The last few decades have witnessed the rapid development of wind power generation around the world. As one of the state-of-art wind power generation facilities, doubly-fed induction generator (DFIG)-based wind turbines (WTs) have occupied a large market share, owing to the advantages of optimized power capture, variable speed operation with reduced mechanical stress, decreased rating grid-connected voltage source converters (VSCs), and lower costs [1]. With large numbers of DFIG-based WTs integrated into the power grid, the impact of wind generation on power system stability has drawn considerable attention from engineers and researchers. In particular, the stability of wind turbines during weak grid operation is one of the critical concerns.
In a WT, there exist different categories of energy storage elements with different time constants, such as the rotor mass of the generator, the DC-link capacitor and the AC filter inductor. Accordingly, hierarchical control is adopted to regulate the states of these storage elements within a certain level. Owing to this multi-timescale control feature, the control dynamics of a WT can be classified into three layers: rotor speed control timescale (1 Hz), DC-link voltage control (DVC) timescale (10 Hz), and AC current control timescale (100 Hz) [2]. Previous studies have primarily focused on the impact of DFIGs' dynamics on the stability of power systems in the rotor speed control timescale [3,4], considering that the rotor speed and drive chain control bandwidths of DFIG wind turbines are close to the synchronous generator's rotor swing mode. However, with the rapidly increasing penetration of wind power into the grid, more electromagnetic oscillations of power systems begin to occur in the range of the DVC timescale [5-9].

Dynamic Interactions for Multi-DFIGs in DC-Link Voltage Control Timescale
This section first gives an introduction to the multi timescale dynamics in DFIGs, in which DVC timescale dynamics will be clarified. Then, a preliminary analysis is offered to illustrate the interactions between control loops in a DVC timescale among multi DFIGs, under weak grid conditions. Figure 1 depicts the diagram of a grid-tied wind farm, where the DFIG-based wind turbines are radially connected to the common point bus. Though there exist many control strategies for power converters in different applications, such as the model predictive control [17], power synchronization control [18] and virtual inertia control [19], the vector control is still the most-commonly utilized strategy in the wind power industry [20], and will be considered herein.

DC-Link Voltage Control Timescale Dynamics of DFIGs
Energies 2020, 13, x FOR PEER REVIEW 3 of 18 This section first gives an introduction to the multi timescale dynamics in DFIGs, in which DVC timescale dynamics will be clarified. Then, a preliminary analysis is offered to illustrate the interactions between control loops in a DVC timescale among multi DFIGs, under weak grid conditions. Figure 1 depicts the diagram of a grid-tied wind farm, where the DFIG-based wind turbines are radially connected to the common point bus. Though there exist many control strategies for power converters in different applications, such as the model predictive control [17], power synchronization control [18] and virtual inertia control [19], the vector control is still the most-commonly utilized strategy in the wind power industry [20], and will be considered herein. As depicted in Figure 1, the vector control of a DFIG is hierarchically formulated, and referring to [15], the associated dynamics can be classified into three layers, namely, the rotor speed control timescale (1 s), the DC-link voltage control timescale (0.1 s) and the AC current control timescale (0.01 s). The dynamics of the rotor speed control timescale are close to the electromechanical dynamics of the power systems, since its time constant is close to the synchronous generator's inertia constant. For the typical vector-controlled VSCs in wind turbine applications, the control bandwidths of the inner current control loops of the RSC and GSC are around 100 Hz, while the outer control loops are usually designed to a decade slower [21], around 10 Hz. Generally, the RSC controls the total active power sent to the grid, as well as regulating the reactive power to maintain the terminal voltage around the nominal value. By contrast, the GSC is responsible for controlling the DC-link voltage, keeping it around the rated value. Besides, unity power factor control is normally adopted by GSC. The frequency range of dynamics covering the outer control loops, such as the GSC DVC, the PLL, the RSC active power control and the terminal voltage control, are defined as the DC-link voltage control timescale dynamics [15]. As depicted in Figure 1, the vector control of a DFIG is hierarchically formulated, and referring to [15], the associated dynamics can be classified into three layers, namely, the rotor speed control timescale (1 s), the DC-link voltage control timescale (0.1 s) and the AC current control timescale (0.01 s). The dynamics of the rotor speed control timescale are close to the electromechanical dynamics of the power systems, since its time constant is close to the synchronous generator's inertia constant. For the typical vector-controlled VSCs in wind turbine applications, the control bandwidths of the inner current control loops of the RSC and GSC are around 100 Hz, while the outer control loops are usually designed to a decade slower [21], around 10 Hz. Generally, the RSC controls the total active power sent to the grid, as well as regulating the reactive power to maintain the terminal voltage around the nominal value. By contrast, the GSC is responsible for controlling the DC-link voltage, keeping it around the rated value. Besides, unity power factor control is normally adopted by GSC. The frequency range of dynamics covering the outer control loops, such as the GSC DVC, the PLL, the RSC active power control and the terminal voltage control, are defined as the DC-link voltage control timescale dynamics [15].

DC-Link Voltage Control Timescale Interactions of Multi DFIGs in Weak Grid Condition
In China, as well as many other places in the world, large numbers of wind plants are located far from the load center, resulting in power delivery through a long AC transmission line [3]. Meanwhile, a wind plant usually contains many identical machines parallel-connected to the point of common connection (PCC) [9,22]. Therefore, the stability study of the system configured as depicted in Figure 2 is of practical significance.

DC-Link Voltage Control Timescale Interactions of Multi DFIGs in Weak Grid Condition
In China, as well as many other places in the world, large numbers of wind plants are located far from the load center, resulting in power delivery through a long AC transmission line [3]. Meanwhile, a wind plant usually contains many identical machines parallel-connected to the point of common connection (PCC) [9,22]. Therefore, the stability study of the system configured as depicted in Figure 2 is of practical significance. The IEEE Task Force suggests the use of a short-circuit ratio (SCR) to prescribe the relative strength of the AC grid to a DC terminal [23]. Quantitatively, this index is calculated as SCR = V 2 /(ZgPrated), in which Zg is the grid impedance and Prated is the rated power of the DC terminal. So, the strength of the AC grid tends to be weaker for longer transmission lines. In the field, some wind plants can be confronted with very weak grid situations (SCR < 2) [24].
The interactions of multi DFIGs in the DVC timescale under weak grid conditions can be complicated. As shown in Figure 2, when the AC transmission line is long, the bus voltage at the PCC, as well as the DFIGs' terminal voltages, will become sensitive to disturbances from both the wind turbine side and the grid side. On the one hand, the DC-link voltage control, the PLL, the active power control and the terminal voltage control will interact with each other through the coupling of terminal voltage dynamics in a single DFIG, due to their close control bandwidths. On the other hand, the variations of the output powers of each DFIG will affect the PCC voltage and the terminal voltage of other DFIGs, considering the circuit power flows; in turn, a change of terminal voltage can drive the DVC control actions of other DFIGs, as well as altering their output powers. Consequently, multi DFIGs also interact with each other in the DVC timescale. Under weak grid conditions, the terminal voltage of a DFIG and the PCC's voltage are more sensitive to the wind turbine's power variations; therefore, the weaker grid condition will aggravate the DVC timescale interactions among multiple DFIGs.

Modeling of Multi DFIGs Connected to Grid
In this section, a detailed modeling process for multi DFIGs parallel-connected to an AC grid is presented. Since this paper mainly studies the dynamic characteristics in the DVC timescale, referring to the modeling consideration for a single DFIG in [15], the following assumptions are made: (1) Both the stator and rotor windings of the generator adopt the motor convention. The fast flux dynamics are ignored, and an algebraic generator model is utilized; (2) The fast dynamic processes of the RSC and GSC's inner current loops are neglected; (3) As the mechanical rotor speed control of wind turbines is much slower than the converter control dynamics, its output power reference (see Figure 1) is assumed to be constant; The IEEE Task Force suggests the use of a short-circuit ratio (SCR) to prescribe the relative strength of the AC grid to a DC terminal [23]. Quantitatively, this index is calculated as SCR = V 2 /(Z g P rated ), in which Z g is the grid impedance and P rated is the rated power of the DC terminal. So, the strength of the AC grid tends to be weaker for longer transmission lines. In the field, some wind plants can be confronted with very weak grid situations (SCR < 2) [24].
The interactions of multi DFIGs in the DVC timescale under weak grid conditions can be complicated. As shown in Figure 2, when the AC transmission line is long, the bus voltage at the PCC, as well as the DFIGs' terminal voltages, will become sensitive to disturbances from both the wind turbine side and the grid side. On the one hand, the DC-link voltage control, the PLL, the active power control and the terminal voltage control will interact with each other through the coupling of terminal voltage dynamics in a single DFIG, due to their close control bandwidths. On the other hand, the variations of the output powers of each DFIG will affect the PCC voltage and the terminal voltage of other DFIGs, considering the circuit power flows; in turn, a change of terminal voltage can drive the DVC control actions of other DFIGs, as well as altering their output powers. Consequently, multi DFIGs also interact with each other in the DVC timescale. Under weak grid conditions, the terminal voltage of a DFIG and the PCC's voltage are more sensitive to the wind turbine's power variations; therefore, the weaker grid condition will aggravate the DVC timescale interactions among multiple DFIGs.

Modeling of Multi DFIGs Connected to Grid
In this section, a detailed modeling process for multi DFIGs parallel-connected to an AC grid is presented. Since this paper mainly studies the dynamic characteristics in the DVC timescale, referring to the modeling consideration for a single DFIG in [15], the following assumptions are made: (1) Both the stator and rotor windings of the generator adopt the motor convention. The fast flux dynamics are ignored, and an algebraic generator model is utilized; (2) The fast dynamic processes of the RSC and GSC's inner current loops are neglected; (3) As the mechanical rotor speed control of wind turbines is much slower than the converter control dynamics, its output power reference (see Figure 1) is assumed to be constant; Energies 2020, 13, 3689 5 of 18 (4) The system is lossless.
There exist various types of PLLs for the synchronization of power converters in different applications [25,26], and the synchronous reference frame (SRF) PLL is the most widely used technique [20] in three-phase power systems. Assuming the terminal voltage angle and the PLL output angle of the DFIG at the synchronous speed coordinate are θ t and θ pll respectively, referring to [15], the SRF PLL model for DFIG1 can be expressed as Combining Equations (1) and (2), the linearized PLL model can be deduced as The voltage phase angle relations between the PLL coordinate and the synchronous rotating coordinate are shown in Figure 3. In a steady state, the PLL coordinate coincides with the synchronous rotating coordinate, i.e., θ pll1 = 0, and the d-axis of the PLL coordinate is oriented to the terminal voltage U t . When subjected to a disturbance, the PLL cannot immediately locate the terminal voltage phase angle, and the angle error between the terminal voltage and the d-axis of the phase-locked loop coordinate system becomes θ t0 -θ pll1 . Therefore, in the dynamic process, the relationship of the dq components of a variable projected in the two coordinates can be expressed as where m denotes the number of DFIG, which can be 1, 2, 3 . . . n.
Energies 2020, 13, x FOR PEER REVIEW 5 of 18 (4) The system is lossless. There exist various types of PLLs for the synchronization of power converters in different applications [25,26], and the synchronous reference frame (SRF) PLL is the most widely used technique [20] in three-phase power systems. Assuming the terminal voltage angle and the PLL output angle of the DFIG at the synchronous speed coordinate are θt and θpll respectively, referring to [15], the SRF PLL model for DFIG1 can be expressed as Combining Equations (1) and (2), the linearized PLL model can be deduced as The voltage phase angle relations between the PLL coordinate and the synchronous rotating coordinate are shown in Figure 3. In a steady state, the PLL coordinate coincides with the synchronous rotating coordinate, i.e., θpll1 = 0, and the d-axis of the PLL coordinate is oriented to the terminal voltage Ut. When subjected to a disturbance, the PLL cannot immediately locate the terminal voltage phase angle, and the angle error between the terminal voltage and the d-axis of the phaselocked loop coordinate system becomes θt0-θpll1. Therefore, in the dynamic process, the relationship of the dq components of a variable projected in the two coordinates can be expressed as where m denotes the number of DFIG, which can be 1, 2, 3…n. Figure 3. Vector diagram representing DFIG1 control dynamics in (a) a steady state and (b) a dynamic process.
Taking the synchronous rotation coordinate of DFIG1 as the global common reference, and assuming that the steady-state phase angle difference between the terminal bus voltage of the DFIGm and that of DFIG1 is θm, the dq components of the variables of the DFIGm projected to this common coordinate can be written as  Taking the synchronous rotation coordinate of DFIG1 as the global common reference, and assuming that the steady-state phase angle difference between the terminal bus voltage of the DFIGm and that of DFIG1 is θ m , the dq components of the variables of the DFIGm projected to this common coordinate can be written as in which f can represent the components of voltages and currents, and the superscripts s1 and sm denote the local synchronous coordinates of DFIG1 and DFIGm, respectively. The circuit equations of the power network (See Figure 2) are Based on Equation (5), Equations (6)-(8) can be manipulated as Dividing the real and imaginary components of (9), we get Linearizing Equations (10) and (11) gives The stator and rotor circuit equations of the m-th DFIG can be expressed as Ignoring the stator resistance and giving the slip ratio of DFIGm the value sr m = (ω 1 − ω rm )/ω 1 , the dq components of the stator and rotor voltages can be further deduced as Synthesizing Equations (12), (13), (16) and (17), it can be derived that Energies 2020, 13, 3689 The output variables of the m-th DFIG can be obtained as The equations related to DC-link voltage control dynamics are The RSC active power control-and terminal voltage control-related equations are Linearizing Equations (29)-(32) gives Energies 2020, 13, 3689 Combining all the linearized equations regarding the DFIG and grid subsystems, the overall small-signal model can be visually formulated by a diagram representation, as given in Figure 4. Here, each DFIG's control system is independently described as a block, containing the DC-link voltage control, terminal voltage control, active power control and phase-locked loop. The power flow equations of the power grid act as the linkage for the interactions of multi-DFIGs. As is shown, each DFIG outputs its currents to the grid block, while the grid block feeds back its power and voltage information for each DFIG's control systems.
Energies 2020, 13, x FOR PEER REVIEW 8 of 18 The small-signal model in Figure 4 can also be expressed in the canonical state-space form as (39) in which ∆x represents the state variables associated with the control dynamics of each DFIG, ∆u are the reference signals and ∆y is the output vector. The eigenvalues and eigenvectors can be calculated from the state matrix A, which will be used to study the small-signal stability of this system.

Stability Analysis of Multi DFIGs Connected to Weak Grids
Based on the proposed small-signal model, this section will investigate the stability of the multi-DFIG system. The commonly-used eigenvalue locus and associated participation factor analysis will be utilized as the analysis tools. Various operating conditions, such as different grid strengths, operating points and controller tunings, will be examined to identify the potential unstable conditions, and explore the impact of these parameters. The small-signal model in Figure 4 can also be expressed in the canonical state-space form as in which ∆x represents the state variables associated with the control dynamics of each DFIG, ∆u are the reference signals and ∆y is the output vector. The eigenvalues and eigenvectors can be calculated from the state matrix A, which will be used to study the small-signal stability of this system.

Stability Analysis of Multi DFIGs Connected to Weak Grids
Based on the proposed small-signal model, this section will investigate the stability of the multi-DFIG system. The commonly-used eigenvalue locus and associated participation factor analysis Energies 2020, 13, 3689 9 of 18 will be utilized as the analysis tools. Various operating conditions, such as different grid strengths, operating points and controller tunings, will be examined to identify the potential unstable conditions, and explore the impact of these parameters.

Eigenvalue and Participation Factor Analysis
Eigenvalues and eigenvectors can be readily calculated from the state matrix as given in Equation (37), and are commonly used to determine system modal characteristics and stability. The participation factor is one of the useful tools in eigenvalue analysis, which is defined as [27] in which φ ki and ψ ik are the kth entry of the right and left eigenvectors, respectively. The participation factor, as expressed in Equation (40), is an index which can weigh the relative participation level of the kth state variable in the ith mode. For a normalized eigenvector, the sum of the participation factors to a mode for each state variable is equal to unity. Besides, the dynamics of a state variable corresponding to a mode are usually linked to the associated controls, for instance, DC-link voltage greatly affects the dynamic of DC-link voltage. Therefore, different controls' relative impacts on one mode can also be quantified by the participation factor. In addition, as the participation factor is calculated from eigenvalues and eigenvectors, its value also changes with system operating condition variations. So, in the following, along with the eigenvalue locus, the participation factor curves under different conditions will be drawn to determine the closely-related controls of the system dominant mode.

Effect of Grid Strengths
For brevity, a system composed of two DFIGs in parallel with an AC grid is considered. The parameters of the two DFIGs are assigned as identical, and given in Appendix A, whereinthe proportional-integral (PI) control parameters primarily refer to the model in Matlab [28].
The effect of grid strengths on stability of this system is studied first. Grid strength is measured by SCR. The active power outputs of DFIG1 and DFIG2 are set to 1.5 MW and 2 MW, respectively. The eigenvalue locus for SCR changing from 4 to 1.1 is given in Figure 5a, and the participation factor curves on system dominant mode for different controls in DFIG1 and DFIG2 are depicted in Figure 5b,c. It can be observed that, when SCR > 1.4, the real part of λ 1, 2 lies closest to the imaginary axis; therefore, λ 1, 2 is dominant when grid strength is relatively stiff. This mode hardly changes with SCR variations. The participation factor curves for SCR > 1.4 reveal that this mode is mainly participated in by the DVC-related state variables. Since the output power of DFIG2 is larger than that of DFIG1, the participation level of DFIG2 is also relatively higher.
However, with the SCR reduced to a certain low level, i.e., SCR < 1.4, another pair of complex modes, λ 5, 6 , sharply moves towards the right half of the plane and becomes the dominant mode, even crossing the imaginary axis when the SCR approaches 1.1. The participation factor curves in Figure 5b,c show that this mode is mainly participated in by the PLLs of the two DFIGs. Therefore, it can be concluded that the interactions of the PLLs of the two generators play a leading role in the stability of weak-grid tied multi-DFIG systems. Additionally, when the SCR is around 1.4, the locations of λ 1, 2 and λ 5, 6 are very close in the complex plane, and it can be found that both the DVC and PLL have high levels of participation in the two modes; by contrast, when λ 1, 2 and λ 5, 6 are separated in comparably high or low SCR conditions, the participation degree of the PLL in λ 1, 2 and that of the DVC in λ 5, 6 become smaller. This result implies that the close interaction between the PLL and the DVC comes under the condition that their related modes are locate close together in the complex plane. Indeed, this phenomenon is known as near modal resonance [29,30], where two eigenvalues located close together in the complex plane will significantly interact with each other, and a slight change of the system parameters could possibly induce a large damping variation in one of the two modes, as well as enlarging the risk of instability. The results of the eigenvalue locus and related participation factor curves in Figure 5 confirm this general phenomenon.
The effect of grid strengths on stability of this system is studied first. Grid strength is measured by SCR. The active power outputs of DFIG1 and DFIG2 are set to 1.5 MW and 2 MW, respectively. The eigenvalue locus for SCR changing from 4 to 1.1 is given in Figure 5a, and the participation factor curves on system dominant mode for different controls in DFIG1 and DFIG2 are depicted in Figure  5b,c. It can be observed that, when SCR > 1.4, the real part of λ1, 2 lies closest to the imaginary axis; therefore, λ1, 2 is dominant when grid strength is relatively stiff. This mode hardly changes with SCR variations. The participation factor curves for SCR > 1.4 reveal that this mode is mainly participated in by the DVC-related state variables. Since the output power of DFIG2 is larger than that of DFIG1, the participation level of DFIG2 is also relatively higher.

Effect of Operating Points
In this subsection, the effect of the operating points on the stability of multi DFIGs is studied. In this case, active power outputs of DFIG2 are set to 2 MW, and the SCR is fixed at 1.1. The eigenvalue locus for DFIG1 s output powers varying from 0.5 p.u. to 1.0 p.u. is depicted in Figure 6a. As the active power output of DFIG1 increases, similarly, the dominant eigenvalues λ 5, 6 move toward the right half of the plane, and migrate into the unstable region, until P e1 = 1.0 p.u. Likewise, participation factor analyses in Figure 6b,c imply that the PLLs occupy the largest participation share of this mode. Besides, the participation level of the RSC APC and the TVC also grow with P e1 's increase. The analysis shows that heavy loading is detrimental to the stability of a multi-DFIG system with weak grid connection.
active power output of DFIG1 increases, similarly, the dominant eigenvalues λ5, 6 move toward the right half of the plane, and migrate into the unstable region, until Pe1 = 1.0 p.u. Likewise, participation factor analyses in Figure 6b,c imply that the PLLs occupy the largest participation share of this mode. Besides, the participation level of the RSC APC and the TVC also grow with Pe1's increase. The analysis shows that heavy loading is detrimental to the stability of a multi-DFIG system with weak grid connection.

Effect of PLL Bandwidths
According to the aforementioned analysis, the dominant mode under weak grid conditions is largely related to the PLLs of the two DFIGs. This section checks the impact of PLL controller tuning. In this case, the active power outputs of DFIG1 and DFIG2 are both set to 2 MW, and the SCR is given as 1.1. The bandwidth of the PLL in DFIG2 is assigned as 14 Hz, and the eigenvalue locus for the variations of the PLL bandwidth of DFIG1 ranging from 2 Hz to 20 Hz is depicted in Figure 7a. The result shows that with the bandwidth of the PLL in DFIG1 approaching that of DFIG2, the dominant mode moves to the right and becomes unstable; however, when the bandwidths of the two PLLs diverge, the system resumes stability. Similar to the case mentioned in Section 4.2, the close bandwidths of the PLLs in the two DFIGs will generate a near modal resonance condition [29,30], which is detrimental to closed-loop stability. Besides, with the increase of the PLL bandwidth in DFIG1, participation factor analysis shows that the participation degree of the PLL in DFIG1 gradually decreases, while the participation of the PLL in DFIG2 gets enlarged. This means that the effectiveness of enhancing the damping of the dominant mode by increasing the PLL bandwidth of one generator is limited, because the modal damping will be governed by the other generator. The case study implies that separating the bandwidths of the PLLs in the two DFIGs is beneficial for system stability.
DFIG1, participation factor analysis shows that the participation degree of the PLL in DFIG1 gradually decreases, while the participation of the PLL in DFIG2 gets enlarged. This means that the effectiveness of enhancing the damping of the dominant mode by increasing the PLL bandwidth of one generator is limited, because the modal damping will be governed by the other generator. The case study implies that separating the bandwidths of the PLLs in the two DFIGs is beneficial for system stability.

Simulation Validations
In this section, detailed model-based time domain simulations are carried out to verify the analysis given in Section 4. A two-DFIG system in connection with a weak grid system model is built in Matlab/Simulink. The main structure of the DFIG model is inherited from the demo model in Simulink [28], and as this paper focuses on converter electromagnetic dynamics, the mechanical part, including the turbine and related controls, in the demo model is abandoned, and constant power/voltage references as converter control inputs are assumed (See Figure 1). Additionally, compared to the proposed reduced-order model, the detailed simulation model includes the dynamics of current controls, the generator and the network circuit. The model's parameters and

Simulation Validations
In this section, detailed model-based time domain simulations are carried out to verify the analysis given in Section 4. A two-DFIG system in connection with a weak grid system model is built in Matlab/Simulink. The main structure of the DFIG model is inherited from the demo model in Simulink [28], and as this paper focuses on converter electromagnetic dynamics, the mechanical part, including the turbine and related controls, in the demo model is abandoned, and constant power/voltage references as converter control inputs are assumed (See Figure 1). Additionally, compared to the proposed reduced-order model, the detailed simulation model includes the dynamics of current controls, the generator and the network circuit. The model's parameters and operating conditions are set to be the same as those in the case studies given in Section 4 for comparison. For all the simulation scenarios, the system reaches a steady state before t = 13 s; then, a small disturbance, i.e., a 0.5% step-up increase of the active power reference in DFIG1, is enacted at t = 13 s, and the responses of DC-link voltages and the active powers of the two DFIGs are captured. Figure 8 shows the simulation results for different grid strength conditions. The DC voltages and active powers of the two generators oscillate with the decrease of SCR, which indicates that the system stability deteriorates with the reduction of grid strength. This tendency is in accord with the eigenvalue locus given in Figure 5a. In addition, Figure 8 shows that sustained oscillation happens at SCR = 1.1, and this also strongly agrees with the analysis given in Figure 5a, where the poles move to the right-half of the plane at SCR = 1.1. Figure 9 gives the simulation results for different output powers of DFIG1. Similarly, it reveals that system stability is reduced with the increase of the loading level, and a loss of stability even happens when the power outputs of DFIG reach 1 p.u. This result is also in good agreement with the analysis given in Figure 6a. Figure 10 shows the simulation results for different PLL bandwidths. As analyzed in Section 4.4, the worst stability condition comes when the bandwidths of the PLLs in the two DFIGs are close. This conclusion is also validated in the simulation results given in Figure 10; with the increase of the PLL bandwidth in DFIG1, the system first becomes unstable after DFIG1 s PLL bandwidth reaches 8 Hz, then it restores stability when DFIG1 s PLL Energies 2020, 13, 3689 13 of 18 bandwidth further increases to 20 Hz. This is because the difference between the bandwidths of the two PLLs first decreased, then grew, during this process. Overall, the simulation studies verify the correctness of the eigenvalue analysis given in Section 4. powers of DFIG1. Similarly, it reveals that system stability is reduced with the increase of the loading level, and a loss of stability even happens when the power outputs of DFIG reach 1 p.u.. This result is also in good agreement with the analysis given in Figure 6a. Figure 10 shows the simulation results for different PLL bandwidths. As analyzed in Section 4.4, the worst stability condition comes when the bandwidths of the PLLs in the two DFIGs are close. This conclusion is also validated in the simulation results given in Figure 10; with the increase of the PLL bandwidth in DFIG1, the system first becomes unstable after DFIG1′s PLL bandwidth reaches 8 Hz, then it restores stability when DFIG1′s PLL bandwidth further increases to 20 Hz. This is because the difference between the bandwidths of the two PLLs first decreased, then grew, during this process. Overall, the simulation studies verify the correctness of the eigenvalue analysis given in Section 4.   level, and a loss of stability even happens when the power outputs of DFIG reach 1 p.u.. This result is also in good agreement with the analysis given in Figure 6a. Figure 10 shows the simulation results for different PLL bandwidths. As analyzed in Section 4.4, the worst stability condition comes when the bandwidths of the PLLs in the two DFIGs are close. This conclusion is also validated in the simulation results given in Figure 10; with the increase of the PLL bandwidth in DFIG1, the system first becomes unstable after DFIG1′s PLL bandwidth reaches 8 Hz, then it restores stability when DFIG1′s PLL bandwidth further increases to 20 Hz. This is because the difference between the bandwidths of the two PLLs first decreased, then grew, during this process. Overall, the simulation studies verify the correctness of the eigenvalue analysis given in Section 4.

Conclusions
In this paper, small-signal modeling and stability analysis for multi DFIGs tied to weak grids in the DC-link voltage control timescale are studied. The proposed model neglects the dynamics of the converter current's controls, which significantly reduced the order of the studied multi-DFIG system. Besides, the model can suitably allow for the possible interactions of multi DFIGs, covering the DVC timescale dynamics. In terms of sensitivity analysis, except for eigenvalue locus, the participation factor curves are also drawn to reveal the weight of the effect of different controls on the system dominant modes for various conditions. The impact of grid strengths, operating points and the bandwidths of PLL are mainly examined, and the major conclusions are drawn, as follows:  For a relatively strong grid, the system's dominant mode is highly participated in by the DVC, and this mode is hardly affected by other controls and grid strength variations. However, in very

Conclusions
In this paper, small-signal modeling and stability analysis for multi DFIGs tied to weak grids in the DC-link voltage control timescale are studied. The proposed model neglects the dynamics of the converter current's controls, which significantly reduced the order of the studied multi-DFIG system. Besides, the model can suitably allow for the possible interactions of multi DFIGs, covering the DVC timescale dynamics. In terms of sensitivity analysis, except for eigenvalue locus, the participation factor Energies 2020, 13, 3689 14 of 18 curves are also drawn to reveal the weight of the effect of different controls on the system dominant modes for various conditions. The impact of grid strengths, operating points and the bandwidths of PLL are mainly examined, and the major conclusions are drawn, as follows: • For a relatively strong grid, the system's dominant mode is highly participated in by the DVC, and this mode is hardly affected by other controls and grid strength variations. However, in very weak grids, the PLL-related mode becomes dominant, and has the risk of instability; • The damping of the dominant mode for a multi-DFIG system under weak grid condition is reduced with the decline of the grid's strength or the increase of the DFIG's power outputs.

•
In ultra-weak grid conditions, the stability of a two-DFIG system gets worse when the bandwidths of the PLLs of the two DFIGs are close, reaching a condition of near modal resonance. In contrast, to improve system stability, the bandwidths of the PLLs in different wind generators should be separated.
The current work is limited to investigating possible instabilities and identify critical factors; the countermeasure regarding the stability enhancement will be carried out in the near future. In addition, the interactions among DFIGs with different control strategies, e.g., virtual synchronous control and vector control, also deserve attention, since these controls are foreseen to coexist in the power systems of the future.

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