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

Abstract

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 firstly 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 find that instability can occur at high levels of power generations or low short-circuit ratio (SCR) conditions. In addition, the dominant mode is identified to be highly related to the PLL, and its modal damping is decreased when the bandwidths of PLLs in different generators are close. Detailed model-based time domain simulations verified the analysis above.

1. 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,6,7,8,9].

As indicated in [5], an active power oscillation of 15 Hz happens in offshore wind farms connected to VSC-HVDC in the North Sea. Furthermore, the authors in [6] prove that interactions between the wind farm and VSC-HVDC will contribute to generating electrical oscillations of around 10 Hz. Meanwhile, the oscillations of wind farms in DVC timescale have led to some serious accidents around the world. In December 2009, oscillations of around 20 Hz occurred in wind farms in Texas, which resulted in a great numbers of wind turbines tripping off the grid [7]. In December 2012, 6–8 Hz oscillations were observed in the northern area of China [8], causing many wind generators to trip. Additionally, as reported in [9], the wind farm in the Guyuan area, located in the northwest part of Hebei Province, China, also witnessed an oscillation of 8.1 Hz in March 2013, causing a loss of 66% of the wind power. These practical cases indicate that the DVC timescale stability of wind turbines deserves great attention.

Recently, some researchers have begun to focus on DVC timescale modeling and the stability analysis of various power converter-based devices, especially under the weak grid condition [10,11,12,13,14,15]. With respect to modeling, the authors in the literature [10] have proposed a small-signal model for a weak-grid-connected VSC, dedicated to DVC stability analysis. The model considers the DVC timescale dynamics, including the DVC, terminal voltage control (TVC) and phase-locked loop (PLL), while neglecting the fast inner current control dynamics. Likewise, a reduced small-signal model of a grid-connected DFIG, for DVC timescale stability analyses, is proposed in [11], and eigenvalue analysis implies that the PLL-related mode increases the risk of instability under a weak grid condition. Lately, a frequency-domain, small-signal model of multi-terminal high-voltage direct current (MTDC) systems in the DVC timescale has also been proposed in the literature [12]. In terms of dynamics analysis, the motion equation concept is extended to describe VSC DVC dynamics in [13,14], where the damping and synchronizing powers of a VSC are computed for quantified evaluation. Similarly, the DVC timescale stability for a DFIG is formally put forward in [15], and based on the concept of motion equation, the equivalent inertia, damping and synchronization coefficients are derived to characterize the DFIG DVC’s dynamic properties. However, for multi DFIGs, to the knowledge of the authors, modeling and stability analysis considering DVC timescale dynamics have not yet been deeply investigated.

The work presented herein is an extension of [16], where a system comprised of two DFIGs in parallel connection to a weak grid is investigated. The root locus analysis in [16] showed that the interactions of the two DFIGs can lead to instability when the AC grid strength is weak to a certain extent. Differently from [16], a comprehensive small-signal modeling process of multi DFIGs, covering DVC timescale dynamics, will be presented in detail herin. The model can count the dynamic interactions between the rotor side converter (RSC) active power control (APC), the TVC and the grid side converter (GSC) DVC, and the PLL of different DFIGs. Besides, in terms of analysis, except for eigenvalue locus sensitivity studies, the participation factor curves are also drawn for comparative study, which can reveal the relative impact of different controls on the system’s dominant mode in a transparent way, as well as providing guidance on controller parameter tunings for stability improvement.

This paper is organized as follows. The problem setup is given in Section 2, where the DVC timescale dynamics of grid-tied DFIGs are introduced, and a preliminarily analysis of the interaction relations of multi DFIGs in weak grids is also offered. Then in Section 3, a reduced-order, small-signal model of multi DFIGs for DVC timescale dynamics analysis is proposed. In Section 4, the eigenvalue locus and participation factor curves of multi control loops are illustrated to reveal the stability characteristics of multi DFIGs, considering different grid strengths, operating points and different bandwidths of PLL. In Section 5, detailed model-based time domain simulations are conducted to validate the analysis above. Conclusions are drawn in Section 6.

2. 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.

2.1. DC-Link Voltage Control Timescale Dynamics of DFIGs

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].

2.2. 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²/(Z_gP_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.

3. 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;

(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

$$u_{\mathrm{tq}1}^{\mathrm{p}1}=U_{\mathrm{t}1}\sin ({\theta }_{\mathrm{t}0}-{\theta }_{\mathrm{pll}1})$$

(1)

$${\theta }_{\mathrm{pll}1}=\frac{u_{\mathrm{tq}1}^{\mathrm{p}1}}{U_{\mathrm{t}1}}(k_{\mathrm{p}1}^{\mathrm{pll}}+\frac{k_{\mathrm{i}1}^{\mathrm{pll}}}{s})\frac{1}{s}$$

(2)

Combining Equations (1) and (2), the linearized PLL model can be deduced as

$$\Delta {\theta }_{\mathrm{pll}1}=\frac{\Delta u_{\mathrm{tq}1}^{\mathrm{p}1}}{U_{\mathrm{t}1}}(k_{\mathrm{p}1}^{\mathrm{pll}}+\frac{k_{\mathrm{i}1}^{\mathrm{pll}}}{s})\frac{1}{s}$$

(3)

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

$$\{\begin{array}{l}{f}_{\mathrm{d}}^{\mathrm{sm}}=\cos {\theta }_{\mathrm{pllm}}{f}_{\mathrm{d}}^{\mathrm{pm}}-\sin {\theta }_{\mathrm{pllm}}{f}_{\mathrm{q}}^{\mathrm{pm}}\\ f_{\mathrm{q}}^{\mathrm{sm}}=\sin {\theta }_{\mathrm{pllm}}{f}_{\mathrm{d}}^{\mathrm{pm}}+\cos {\theta }_{\mathrm{pllm}}{f}_{\mathrm{q}}^{\mathrm{pm}}\end{array}$$

(4)

where m denotes the number of DFIG, which can be 1, 2, 3…n.

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

$$\{\begin{array}{l}{f}_{\mathrm{d}}^{\mathrm{s}1}=\cos {\theta }_{\mathrm{m}}{f}_{\mathrm{d}}^{\mathrm{sm}}-\sin {\theta }_m{f}_{\mathrm{q}}^{\mathrm{sm}}\\ f_{\mathrm{q}}^{\mathrm{s}1}=\sin {\theta }_{\mathrm{m}}{f}_{\mathrm{d}}^{\mathrm{sm}}+\cos {\theta }_m{f}_{\mathrm{q}}^{\mathrm{sm}}\end{array}$$

(5)

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

$${\overrightarrow{U}}_{\mathrm{tm}}-{\overrightarrow{U}}_{\mathrm{pcc}}=j{X}_{\mathrm{tm}}{\overrightarrow{I}}_{\mathrm{tm}}\begin{array}{cc},& m=1,2,\dots ,\end{array}n$$

(6)

$${\overrightarrow{I}}_{\mathrm{tm}}={\overrightarrow{I}}_{\mathrm{sm}}+{\overrightarrow{I}}_{\mathrm{cm}}\begin{array}{cc},& m=1,2,\dots ,\end{array}n$$

(7)

$${\overrightarrow{U}}_{\mathrm{pcc}}-{\overrightarrow{U}}_{\mathrm{g}}=j{X}_{\mathrm{g}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\overrightarrow{I}}_{\mathrm{ti}}}$$

(8)

Based on Equation (5), Equations (6)–(8) can be manipulated as

$$u_{\mathrm{tdm}}^{\mathrm{s}1}+j{u}_{\mathrm{tqm}}^{\mathrm{s}1}=j{X}_{\mathrm{tm}}(i_{\mathrm{tdm}}^{\mathrm{s}1}+j{i}_{\mathrm{tqm}}^{\mathrm{s}1})+j{X}_{\mathrm{g}}{\displaystyle \sum _{i=1}^n(i_{\mathrm{tdi}}^{\mathrm{s}1}+j{i}_{\mathrm{tqi}}^{\mathrm{s}1})}+u_{\mathrm{gd}}^{\mathrm{s}1}+j{u}_{\mathrm{gq}}^{\mathrm{s}1}$$

(9)

Dividing the real and imaginary components of (9), we get

$$u_{\mathrm{tdm}}^{\mathrm{s}1}=-X_{\mathrm{tm}}(i_{\mathrm{sqm}}^{\mathrm{s}1}+i_{\mathrm{cqm}}^{\mathrm{s}1})-X_{\mathrm{g}}{\displaystyle \sum _{i=1}^n(i_{\mathrm{sqi}}^{\mathrm{s}1}+i_{\mathrm{cqi}}^{\mathrm{s}1})}+u_{\mathrm{gd}}^{\mathrm{s}1}$$

(10)

$$u_{\mathrm{tqm}}^{\mathrm{s}1}=X_{\mathrm{tm}}(i_{\mathrm{sdm}}^{\mathrm{s}1}+i_{\mathrm{cdm}}^{\mathrm{s}1})+X_{\mathrm{g}}{\displaystyle \sum _{i=1}^n(i_{\mathrm{sdi}}^{\mathrm{s}1}+i_{\mathrm{cdi}}^{\mathrm{s}1}})+u_{\mathrm{gq}}^{\mathrm{s}1}$$

(11)

Linearizing Equations (10) and (11) gives

$$\Delta u_{\mathrm{tdm}}^{\mathrm{s}1}=-X_{\mathrm{tm}}(\Delta i_{\mathrm{sqm}}^{\mathrm{s}1}+\Delta i_{\mathrm{cqm}}^{\mathrm{s}1})-X_{\mathrm{g}}{\displaystyle \sum _{i=1}^n(\Delta i_{\mathrm{sqi}}^{\mathrm{s}1}+\Delta i_{\mathrm{cqi}}^{\mathrm{s}1})}$$

(12)

$$\Delta u_{\mathrm{tqm}}^{\mathrm{s}1}=X_{\mathrm{tm}}(\Delta i_{\mathrm{sdm}}^{\mathrm{s}1}+\Delta i_{\mathrm{cdm}}^{\mathrm{s}1})+X_{\mathrm{g}}{\displaystyle \sum _{i=1}^n(\Delta i_{\mathrm{sdi}}^{\mathrm{s}1}+\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}})$$

(13)

The stator and rotor circuit equations of the m-th DFIG can be expressed as

$${\overrightarrow{\mathit{U}}}_{\mathrm{tm}}=-j{\omega }_1{L}_{\mathrm{mm}}{\overrightarrow{\mathit{I}}}_{\mathrm{rm}}-j{\omega }_1{L}_{\mathrm{sm}}{\overrightarrow{\mathit{I}}}_{\mathrm{sm}}-R_{\mathrm{s}}{\overrightarrow{\mathit{I}}}_{\mathrm{sm}}$$

(14)

$${\overrightarrow{\mathit{U}}}_{\mathrm{rm}}=-R_{\mathrm{rm}}{\overrightarrow{\mathit{I}}}_{\mathrm{rm}}+j({\omega }_1-{\omega }_{\mathrm{rm}})(-L_{\mathrm{rm}}{\overrightarrow{\mathit{I}}}_{\mathrm{rm}}-L_{\mathrm{mm}}{\overrightarrow{\mathit{I}}}_{\mathrm{sm}})$$

(15)

Ignoring the stator resistance and giving the slip ratio of DFIGm the value sr_m = (ω₁ − ω_rm)/ω₁, the dq components of the stator and rotor voltages can be further deduced as

$$\Delta u_{\mathrm{tdm}}^{\mathrm{s}1}=X_{\mathrm{mm}}\Delta i_{\mathrm{rqm}}^{\mathrm{s}1}+X_{\mathrm{sm}}\Delta i_{\mathrm{sqm}}^{\mathrm{s}1}$$

(16)

$$\Delta u_{\mathrm{tqm}}^{\mathrm{s}1}=-X_{\mathrm{mm}}\Delta i_{\mathrm{rdm}}^{\mathrm{s}1}-X_{\mathrm{sm}}\Delta i_{\mathrm{sdm}}^{\mathrm{s}1}$$

(17)

$$\Delta u_{\mathrm{rdm}}^{\mathrm{s}1}=-R_{\mathrm{rm}}\Delta i_{\mathrm{rdm}}^{\mathrm{s}1}+s{r}_{\mathrm{m}}{X}_{\mathrm{mm}}\Delta i_{\mathrm{sqm}}^{\mathrm{s}1}+s{r}_{\mathrm{m}}{X}_{\mathrm{rm}}\Delta i_{\mathrm{rqm}}^{\mathrm{s}1}$$

(18)

$$\Delta u_{\mathrm{rqm}}^{\mathrm{s}1}=-R_{\mathrm{rm}}\Delta i_{\mathrm{rqm}}^{\mathrm{s}1}-s{r}_{\mathrm{m}}{X}_{\mathrm{mm}}\Delta i_{\mathrm{sdm}}^{\mathrm{s}1}-s{r}_{\mathrm{m}}{X}_{\mathrm{rm}}\Delta i_{\mathrm{rdm}}^{\mathrm{s}1}$$

(19)

Synthesizing Equations (12), (13), (16) and (17), it can be derived that

$$\Delta u_{\mathrm{tdm}}^{\mathrm{s}1}={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i-1}^{utd}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i}^{utd}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}}$$

(20)

$$\Delta u_{\mathrm{tqm}}^{\mathrm{s}1}={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i-1}^{utq}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i}^{utq}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}}$$

(21)

$$\Delta i_{\mathrm{sdm}}^{\mathrm{s}1}={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i-1}^{is}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i}^{is}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}}$$

(22)

$$\Delta i_{\mathrm{sqm}}^{\mathrm{s}1}={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i-1}^{is}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,2i}^{is}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}}$$

(23)

The output variables of the m-th DFIG can be obtained as

$$\begin{array}{l}\Delta P_{\mathrm{em}}=1.5(\Delta u_{\mathrm{tdm}}^{\mathrm{s}1}{i}_{\mathrm{tdm}0}^{\mathrm{s}1}+u_{\mathrm{tdm}0}^{\mathrm{s}1}\Delta i_{\mathrm{tdm}}^{\mathrm{s}1}+\Delta u_{\mathrm{tqm}}^{\mathrm{s}1}{i}_{\mathrm{tqm}0}^{\mathrm{s}1}+u_{\mathrm{tqm}0}^{\mathrm{s}1}\Delta i_{\mathrm{tqm}}^{\mathrm{s}1})\\ ={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-3}^{\mathrm{pe}}}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-2}^{\mathrm{pe}}}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-1}^{\mathrm{pe}}}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i}^{\mathrm{pe}}}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}\end{array}$$

(24)

$$\begin{array}{l}\Delta P_{\mathrm{cm}}=1.5(\Delta u_{\mathrm{tdm}}^{\mathrm{s}1}{i}_{\mathrm{cdm}0}^{\mathrm{s}1}+u_{\mathrm{tdm}0}^{\mathrm{s}1}\Delta i_{\mathrm{cdm}}^{\mathrm{s}1}+\Delta u_{\mathrm{tqm}}^{\mathrm{s}1}{i}_{\mathrm{cqm}0}^{\mathrm{s}1}+u_{\mathrm{tqm}0}^{\mathrm{s}1}\Delta i_{\mathrm{cqm}}^{\mathrm{s}1})\\ ={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-3}^{\mathrm{pc}}}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-2}^{\mathrm{pc}}}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-1}^{\mathrm{pc}}}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i}^{\mathrm{pc}}}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}\end{array}$$

(25)

$$\begin{array}{l}\Delta P_{\mathrm{rm}}=1.5(\Delta u_{\mathrm{rdm}}^{\mathrm{s}1}{i}_{\mathrm{rdm}0}^{\mathrm{s}1}+u_{\mathrm{rdm}0}^{\mathrm{s}1}\Delta i_{\mathrm{rdm}}^{\mathrm{s}1}+\Delta u_{\mathrm{rqm}}^{\mathrm{s}1}{i}_{\mathrm{rqm}0}^{\mathrm{s}1}+u_{\mathrm{rqm}0}^{\mathrm{s}1}\Delta i_{\mathrm{rqm}}^{\mathrm{s}1})\\ ={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-3}^{\mathrm{pr}}}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-2}^{\mathrm{pr}}}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-1}^{\mathrm{pr}}}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i}^{\mathrm{pr}}}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}\end{array}$$

(26)

$$\Delta U_{\mathrm{tm}}={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-3}^{Ut}}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-2}^{Ut}}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-1}^{Ut}}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i}^{Ut}}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}$$

(27)

$$\begin{array}{ll}\Delta u_{\mathrm{tqm}}^{\mathrm{pm}}& ={\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-3}^{utq}}\Delta i_{\mathrm{rdi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-2}^{utq}}\Delta i_{\mathrm{rqi}}^{\mathrm{s}1}+{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i-1}^{utq}}\Delta i_{\mathrm{cdi}}^{\mathrm{s}1}\\ & +{\displaystyle \sum _{\mathrm{i}=1}^n{k}_{m,4i}^{utq}}\Delta i_{\mathrm{cqi}}^{\mathrm{s}1}-u_{\mathrm{tdm}0}^{\mathrm{pm}}\Delta {\theta }_{\mathrm{pllm}}\end{array}$$

(28)

The equations related to DC-link voltage control dynamics are

$$P_{\mathrm{rm}}-P_{\mathrm{cm}}=U_{\mathrm{dcm}0}{C}_{\mathrm{m}}\frac{d{U}_{\mathrm{d}\mathrm{cm}}}{dt}$$

(29)

$$i_{\mathrm{cdm}}^{\mathrm{pm}}=(U_{\mathrm{dcm}}^{*}-U_{\mathrm{dcm}})(k_{\mathrm{pm}}^{\mathrm{Udc}}+\frac{k_{\mathrm{im}}^{\mathrm{Udc}}}{s})$$

(30)

The RSC active power control- and terminal voltage control-related equations are

$$i_{\mathrm{rdm}}^{\mathrm{pm}}=(P_{em}^{*}-P_{em})(k_{\mathrm{pm}}^{Pe}+\frac{k_{\mathrm{im}}^{Pe}}{s})$$

(31)

$$i_{\mathrm{rqm}}^{\mathrm{pm}}=(U_{\mathrm{tm}}^{*}-U_{\mathrm{tm}})(k_{\mathrm{pm}}^{\mathrm{Ut}}+\frac{k_{\mathrm{im}}^{\mathrm{Ut}}}{s})$$

(32)

Linearizing Equations (29)–(32) gives

$$\Delta P_{\mathrm{rm}}-\Delta P_{\mathrm{cm}}=U_{\mathrm{dcm}0}{C}_{\mathrm{m}}\frac{d\Delta U_{\mathrm{d}\mathrm{cm}}}{dt}$$

(33)

$$\Delta i_{\mathrm{cdm}}^{\mathrm{pm}}=(\Delta U_{\mathrm{dcm}}^{*}-\Delta U_{\mathrm{dcm}})(k_{\mathrm{pm}}^{\mathrm{Udc}}+\frac{k_{\mathrm{im}}^{\mathrm{Udc}}}{s})$$

(34)

$$\Delta i_{\mathrm{rdm}}^{\mathrm{pm}}=(\Delta P_{\mathrm{em}}^{*}-\Delta P_{em})(k_{\mathrm{pm}}^{Pe}+\frac{k_{\mathrm{im}}^{Pe}}{s})$$

(35)

$$\Delta i_{\mathrm{rqm}}^{\mathrm{pm}}=(\Delta U_{\mathrm{tm}}^{*}-\Delta U_{\mathrm{tm}})(k_{\mathrm{pm}}^{\mathrm{Ut}}+\frac{k_{\mathrm{im}}^{\mathrm{Ut}}}{s})$$

(36)

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.

The small-signal model in Figure 4 can also be expressed in the canonical state-space form as

$$\begin{array}{l}\Delta \dot{x}=\mathbf{A}\Delta \mathbf{x}+\mathbf{B}\Delta \mathbf{u}\\ \Delta \mathbf{y}=\mathbf{C}\Delta \mathbf{x}+\mathbf{D}\Delta \mathbf{u}\end{array}$$

(37)

$$\begin{array}{ll}\Delta \mathbf{x}=& [\Delta i_{\mathrm{rd}1}^{\mathrm{p}1}\Delta i_{\mathrm{rq}1}^{\mathrm{p}1}\Delta i_{\mathrm{cd}1}^{\mathrm{p}1}\Delta U_{\mathrm{dc}1}\Delta {\omega }_{\mathrm{pll}1}\Delta {\theta }_{\mathrm{pll}1}\cdot \cdot \cdot \Delta {\omega }_{\mathrm{plln}}\Delta {\theta }_{\mathrm{plln}}{]}^T\end{array}$$

(38)

$$\Delta \mathbf{u}=\begin{array}{ccc}\begin{array}{cc}[\Delta P_{\mathrm{e}1}^{\ast }& \Delta U_{\mathrm{t}1}^{\ast }\end{array}& \Delta U_{\mathrm{dc}1}^{*}& \begin{array}{cc}\Delta u_{\mathrm{tq}1}^{\mathrm{p}1\ast }\cdot \cdot \cdot \Delta U_{\mathrm{dcn}}^{*}& \Delta u_{\mathrm{tqn}}^{\mathrm{pn}\ast }\end{array}\end{array}{]}^T$$

(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.

4. 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.

4.1. 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]

$$p_{ki}={\varphi }_{ki}{\psi }_{ik}$$

(40)

in which ${\varphi }_{ki}$ and ${\psi }_{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.

4.2. 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, wherein the 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.

4.3. 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.

4.4. 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.

5. 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 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.

6. 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 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.