# 1. Introduction

In the recent decades, a few serious power system contingencies occurring around the world pertinent to voltage collapse have drawn wide attentions for electrical engineers [1-2]. Over the years, a variety of analyzing methods have been proposed to give insight into the mechanism and assessment of voltage stability. However, due to the complexity of this issue, it is difficult for these attempts to provide the comprehensive and time-saving solutions.

The existing analyzing methods can be classified into two categories: static and dynamic. So far, the static methods, which can give instructive information for power system planning, operation and control with ignorance of transient behavior, are widely applied in utilities. Especially, the conventional PV and QV curves based analysis methods [3-6], which are applied to investigate the relationships between voltage and active or reactive power in a specific load increasing direction, are very useful in determining voltage stability margin and developing remedial action strategies. However, such methods are only applicable for certain increasing modes and cannot provide a full view of voltage stability. On the other hand, the utilized continuation power flow technique is rather time-consuming. Furthermore, with the integration of large-scale renewable energy sources, especially for wind power, the new challenges and issues to the voltage stability will appear in the smart grid environment. Three factors, the variability and intermittency of wind energy, the wind turbine generator (WTG) models, and the pattern of grid-connection, should be taken into account for voltage stability studies incorporating wind power [7-9]. In a nutshell, how to effectively conduct the voltage stability analysis incorporating intermittent wind power and develop more powerful analysis tool become very imperative and necessary.

In this paper, an analytical method called as PQ Plane analysis, is proposed for voltage stability studies incurporating wind power. Two factors, the reactive power capability of WTG and the wake effect of wind farm, are paid special attention to during the modelling of wind farms of doubly-fed induction generator (DFIG) type. The stability boundary and voltage distribution are derived and shown on the PQ plane. Furthermore, a set of indices are defined to indicate the voltage security situation of a load bus. The proposed method is derived in a single-load-infinite-bus system, and can be easily expanded to complex power systems by applying Thevenin Equivalent Circuits [10] or the newly proposed Coupled Single-Port Circuits [11]. Some preliminary and meaningful conclusions are drawn.

The rest of the paper is organized as follows: In Section 2, the static modelling of wind farm involving the reactive power capability of DFIG and the wake effect is discussed. The PQ plane analysis method incorporating wind generation is proposed and detailed in Section 3. The case studies with different simulation scenarios are carried out in Section 4. The conclusions and research prospects can be found in Section 5.

# 2. Modelling of Wind Farms for Voltage Stability Studies

In this section, a static model of wind farm of DFIG type is proposed with consideration of reactive power capability and wake effect.

## 2.1 Active power output of wind farms with consideration of wake effect

The active power output of a DFIG Pwt with respect to the wind speed v can be expressed as [12]

where vci, vco, vr and Pwtr are cut-in speed, cut-out speed, rated speed and rated active power, respectively.

In conventional model of wind farm, the wind speed each wind turbine captures is assumed to be the same. Actually, in practical wind farms, the wind turbines located downwind gain a lower wind speed than those located upwind, in that the kinetic energy of wind is partly consumed by the latter ones. Such aforementioned situation is called as wake effect. A lot of wake effect models have been developed to date with specific applicability. In this paper, the Jensen Model [13], as illustrated in Fig. 1 is applied to describe the wake effect in a wind farm located on flat terrain. In accordance with the Jensen model, the wind speed behind a wind turbine vwake with respect to the distance from the wind turbine plane x can be presented as

**Fig. 1.**The wake effect of wind turbines

where v1 is the original wind speed; Rwt is the wind turbine radius; ξx is the damping coefficient; and ξT is the thrust coefficient of the wind turbine. ξx can be obtained by

where kw is the turbulence intensity, which can be estimated empirically; σwake and σ0 are mean-square deviation of natural turbulence and turbulence caused by WTGs, respectively.

For a wind farm with nrow rows and ncol columns in a chessboard-pattern layout, the distance between two adjacent rows of WTGs is drow, the wind speed in the ith row can be calculated by

where v1 is the wind speed in the first row (the original wind speed), and ξd is the damping coefficient between two adjacent rows of WTGs, which is given by

Then the active power output of a single DFIG in the ith row Pwt,i can be calculated based on Eq. (1). Ignoring the losses in the collection systems, the total output of a wind farm equipped with DFIG is

Fig. 2 illustrates the active power output of a wind farm with and without considering the wake effect.

**Fig. 2.**Active power output of wind farms with respect to wind speed

## 2.2 Reactive power capability curve of DFIG

The equivalent circuit of a DFIG is given in Fig. 3, in which the stator voltage , stator current , rotor voltage and rotor current are restricted by

**Fig. 3.**Equivalent circuit of DFIG

where s is the slip ratio.

From (7) and (8), we have

where matrix Y, G and B can be derived from the parameters of the equivalent circuit.

The complex power outputs of the stator and rotor are expressed as

If the losses in the boosting transformer are negligible, the total power output can be given by

where Ps is the stator active power; Qs is the stator reactive power; Pr is the rotor active power; Qgc is the grid-side convertor reactive power.

Among the four variables, , , /s and , is determined by the grid voltage, and the other three variables have certain limits for security reasons, which would determine the reactive power capability. With a given , to analyze the effect of limit, the output should be expressed with respect to . By substituting (10) into (12) and (13), we have

where Zs = Rs + jXs; Zr = Rr /s + jXr; Zm = jXm. From (14), (15) and (16), the total power output can be derived and expressed as

From (17), with given Vs, Ir and s, it can be drawn as an ellipse on the PQ plane. As the absolute value of imaginary part of Zm+Zs is much larger than that of the real part, the center of the ellipse is very close to negative Q-axis.

Similarly, with a given , to analyze the effect of limit, the output should be expressed with respect to . By substituting (11) into (13), we have

From (12), (14) and (18), the total power output can be derived and expressed as

Fig. 4 shows the reactive power capability determined by limits of the stator current and rotor current. The solid line given in Fig. 4 denotes the reactive power capability with respect to the active power output.

**Fig. 4.**Reactive power capability determined by the stator current limit and the rotor current limit

As mentioned above, the slip ratio s is assumed to be constant. In actual operation of DFIGs, however, the rotational speed of the wind turbine changes with the wind speed to capture the maximum power. A commonly-applied speed control law is given as follows [15]:

where nmin is the minimum speed; n1 is the synchronous speed; nmax is the maximum speed; kopt is the power adjustment coefficient; P1, P2 and P3 are three demarcation points, respectively.

As the slip ratio changes, the shape of the reactive power capability curve will change. The corresponding solution procedures are given as follows:

1) Assign the initial active power P0=0 and the increment ΔP; Let i=1; 2) Pi=P0+iΔP; 3) Calculate rotational speed ni according to Eq. (20); 4) Calculate slip ratio si; 5) Calculate the reactive power upper limit Qmax,i and lower limit Qmin,i according to Eq. (17) and (19) respectively; 6) If Qmax,i≤Qmin,i or Pi≥1, go to Step 7; else, i=i+1, go to Step 2; 7) Plot the reactive power capability curve.

The reactive power capability curves with different s are given in Fig. 5. It can be seen that as the slip ratio decreases or the rotational speed increases, the operation region expands and the reactive power capability increases. The actual reactive power capability curve is made up of points in different curves. When P < P1 and P2 ≤ P < P3, the actual reactive power capability curve will overlap curves with s=0.3 and s=0 respectively. In other areas, it crosses a set of curves with constant s as the slip ratio keeps changing.

**Fig. 5.**Capability curves with different constant s and with the speed control law

## 2.3 Reactive power capability curve of wind farm and voltage control strategies

With the consideration of the wake effect, the total power output of wind farm of DFIG-type can be calculated according to the following procedures:

1) Assign the initial wind speed v0=0 and the increment Δv; Let k =1; 2) vk = v0 +kΔv; 3) Calculate the wind speeds vki according to Eqs. (4) and (5); 4) Calculate the active output of DFIGs in each row Pwt, ki according to Eq. (1); 5) Calculate the active output of the whole wind farm Pwf,k according to Eq. (6); 6) Calculate the rotational speeds of DFIGs in each row nki according to Eq. (20); 7) Calculate the slip ratios of DFIGs in each row ski; 8) Calculate the reactive power upper limit Qwtmax,ki and lower limit Qwtmin,ki of DFIGs in each row according to Eqs. (17) and (19), respectively; 9) Calculate the reactive power upper limit Qwfmax,ki and lower limit Qwfmin,ki of the whole wind farm according to the following Equations

10) If vi ≥ vco, go to Step 11; else, i =i +1, go to Step 2; 11) Plot the reactive power capability curve.

Fig. 6 shows the reactive power capability of a wind farm with/without consideration of the wake effect as well as the conventional constant power factor control strategy.

**Fig. 6.**Reactive power capability of a wind farm

From Fig. 6, with conventional control strategy, there is a large amount of reactive reserve unutilized, especially when the active power output is relatively low. To make full use of the reactive reserve for the sake of voltage security enhancement, the constant voltage control mode is applied in this paper. When there is an amount of reactive reserve, the voltage at the point of common coupling (PCC) is maintained at 1 p.u.; when the reactive reserve is depleted, the reactive output of the wind farm is maintained at the maximum value. The corresponding comparative analysis between the two control modes can be found in Section 4.

# 3. PQ Plane Analysis

As mentioned above, the conventional PV and QV curves based analysis methods only reflect the voltage stability in specific load-increasing directions, which are difficult to provide the full picture of voltage stability. In this section, an analytical method called the PQ plane analysis is proposed for voltage stability analysis and assessment incorporating wind generation.

## 3.1 Voltage stability region and voltage profiles

The single-load-infinite-bus system as shown in Fig. 7 is employed as benchmark to address the proposed method. The symbols , , , P , Q, R and X represents the infinite bus voltage, transmission line current, load bus voltage, active load, reactive load, transmission line resistance and reactance, respectively. According to circuit principles, we have

**Fig. 7.**Single-load-infinite-bus system

From the definition of complex power, the current can be expressed as

By substituting (24) into (23), we have

From (25), we can get

Further, we can derive

Eq. (28) can be converted into

With a constant V, the locus of Eq. (29) is a circle on the PQ plane, with the center

and radius

In this paper, the loci of the power system operating points with the specific constant V on the PQ plane are called the constant voltage lines, which shows the range of load that corresponds to secure voltage.

From (28), the voltage stability boundary can also be derived. Rewrite the equation with respect to V, we have

On the voltage stability boundary, there is only one solution for V. Thus, the discriminant of Eq. (32) should be equal to zero:

Accordingly, we can get the equation of voltage stability boundary:

The discriminant of Eq. (34) equals to zero, indicating that the voltage stability boundary is a parabola on the PQ plane.

Fig. 8 displays the voltage stability boundary and a set of constant voltage lines pertinent to the single-load-infinite-bus system in which X is much larger than R. It can be seen that the stability boundary is a parabola. The parts inside the parabola are the stability region. The constant voltage lines are circles tangent to the stability boundary, i.e. the stability boundary is the envelope of all constant voltage lines.

**Fig. 8.**Voltage stability boundary and constant voltage lines

For a load bus, only the first quadrant of PQ plane as shown in Fig. 9 needs to be elaborately studied. It can be seen from Fig. 9 that the voltage constant lines do not intersect in the first quadrant. As the load increases, the operating point moves closer to the stability boundary and the bus experiences a voltage dip. Furthermore, it can also be found that the information provided by PV and QV curves based analysis methods is easily included in PQ plane. For instance, the PV curve based analysis method assumes that the load increases with a given power factor cosφ, and calculates the active power margin to the stability boundary ΔPL. On the PQ plane, the load increasing direction is represented as the direction of , and ΔPL is the length of AC. The QV curve based analysis methods assumes that there is a synchronous condenser installed at the load bus and investigates the relationship between V and Q while maintaining P constant. On the PQ plane, the increasing direction is represented as the direction of , and the reactive power margin ΔQc is the length of AD.

**Fig. 9.**Voltage stability boundary and constant voltage lines in and nearby the first quadrant

## 3.2 Discussions of PV curve based analysis method

According to Fig. 9, for the conventional PV curve based analysis, when there is an error between the actual load-increasing direction and the assumed direction, an error of the load margin will occur. Moreover, even if the calculated load margin ΔPL of two systems are same, they may have different performance in other given load-increasing directions. Fig. 10 shows the stability boundaries of two systems with different parameters: SLIB1 and SLIB2. Their stability boundaries intersect in the first quadrant. Assume that the two systems have the same initial operating point A. With a certain assigned load-increasing direction, they all reach B point, the point of intersection of their stability boundaries as shown in Fig. 10. Obviously, they share the same load margin ΔP1, 2, AB. If the load increases along , the margin of SLIB1 ΔP1, AC is smaller than that of SLIB2 ΔP2, AC; if the load increases along , however, the results are exactly the opposite. Hence, it can be found that the load margin of SLIB1 varies more intensely with respect to the load-increasing direction variations. And it also can be concluded that SLIB2 is more stable than SLIB1 from a voltage stability viewpoint.

**Fig. 10.**Conventional PV curve based analysis

## 3.3 Voltage stability indices

In this section, two voltage stability indices based on PQ plane analysis incorporating the uncertainties of load-increasing direction and wind generation are proposed. Assume that the power factor angle of the load increment φ is a random variable subject to a certain probability distribution φ~Dφ with a mean φμ and a variance φσ. The load margin in the direction of φ can be calculated by solving the following equations:

Two solutions can be obtained from (35), and only the solution named Pmargin in the first quadrant is retained. As Pmargin is a function of φ:

It must subject to a certain probability distribution Pmargin~ DPmargin. In our work, the mean and variance of Pmargin are calculated approximately to avoid the heavy computational burden. According to probability theories, with given φμ and φσ, Pmargin μ and Pmargin, σ can be estimated respectively by

where f’(φμ) and f’’(φμ) can be estimated by difference quotients. Assign a small increment, and substitute φ=φμ-Δφ and φ=φμ+Δφ into Eq. (36) to get f(φμ-Δφ) and f(φμ+Δφ). Then f’(φμ) and f’’(φμ) can be calculated by

The coefficient of variation of Pmargin can be calculated by

We define Pmargin,μ and Pmargin,CV as voltage stability indices, called the Mean Margin to Voltage Instability (MMVI) and the Variation Margin to Voltage Instability (VMVI) respectively. MMVI represents the expectation of the load margin. The larger the index, the more secure the system. VMVI represents the variance of the load margin. The larger the index, the less robust the system is to the randomness of the load.

Conventionally, the solutions of voltage stability analysis mainly focus on how the current state is to the stability boundary. In fact, voltage violation is also an important issue that needs to be considered in practical operation. As the initial operating point and the load-increasing direction vary, there does not exist a fixed relationship between the distance to voltage instability and the distance to voltage violation. It is necessary, therefore, to define separate indices to indicate the risk of voltage violation. Fortunately, it is convenient to implement by using voltage constant lines on the PQ plane. By setting a voltage lower bound Vmin, and plotting the corresponding constant voltage line on the PQ plane, the secure region can be determined, as given in Fig. 11. By solving the following equations:

**Fig. 11.**Margin to voltage instability and margin to voltage violation

The margin to voltage violation with an initial operating point A and load-increasing direction can be obtained.

Similarly, the solution P’margin is a function of φ:

Replacing function f by function g in Eq. (37) through (41), the mean and variance coefficient of P’margin μ, P’margin, CV can be derived. We define P’margin μ as the Mean Margin to Voltage Violation (MMVV) and P’margin, CV the Variation Margin to Voltage Violation (VMVV) to conduct the voltage stability evaluation.

To assess the impact of local wind integration on voltage stability, the single-load-infinite-bus system with wind power integration as shown in Fig. 12 is employed during analysis. The wind generation is connected to the grid at the load bus with Pwind+jQwind power injection. In this test system, the total active and reactive loads are

**Fig. 12.**Single-load-infinite-bus system with wind power integration

**Fig. 13.**The way to incorporate wind power in PQ Plane Analysis

To plot voltage stability boundary and constant voltage lines on the PQ plane based on the net load P and Q, a shift with value of (Pwind, Qwind) is added to implement the required solutions.

# 4. Case Studies

Based on the aforementioned single-load-infinite-bus test system shown in Fig. 7, two simulation scenarios involving the studies on the relationship between power system parameters and voltage stability and the impacts of wind farm of DFIG type on voltage stability are designed for case studies. Some meaningful conclusions are drawn based on the studies.

## 4.1 Impacts of key parameters of power system on voltage stability

Impacts of transmission line ratio of reactance X and resistance R Keeping the infinite bus voltage E and the transmission line modulus of impedance |Z| constant, the stability boundaries with different X/R are plotted in Fig. 14. It can be seen that the opening direction of the parabola rotates clockwise as X/R reduces.

**Fig. 14.**Voltage stability boundaries with different X/R

Fig. 15 shows the influence of X/R on stability boundaries in and nearby the first quadrant. Because of the rotation of the parabola, the part in the first quadrant becomes steeper, and the boundaries with different X/R intersect. This situation may cause the problem as discussed in Section 3.2. Table 1 shows the results of MMVI, VMVI and ΔQc with different X/R. The initial load is set to be 0, and cosφ=cos(π/6)=0.866. Other factors of systems are kept constant. It can be seen that the MMVI gradually reduces as increase of X/R. The VMVI and X/R have a positive correlation, indicating that the variation of load margin reduces as X/R becomes larger. Furthermore, the larger ΔQc does not necessarily result in larger load margin, and that is what we need to be careful with when applying the conventional QV analysis.

**Fig. 15.**Voltage stability boundaries with different X/R in and nearby the first quadrant

**Table 1.**MMVI, VMVI and ΔQc with different X/R

To understand the impact of X/R on voltage profiles, the changes of constant voltage lines must be studied. Fig. 16 plots the constant voltage lines with different X/R, from which we can see that the center of the circle rotates clockwise as decrease of X/R.

**Fig. 16.**Constant voltage lines with different X/R

Fig. 17 shows the influence of X/R on constant voltage lines in and nearby the first quadrant. All of the lines intersect at a common point. Obviously, it can be inferred that the tendencies of MMVV and VMVV are similar to those of the MMVI and VMVI.

**Fig. 17.**Constant voltage lines with different X/R in and nearby the first quadrant

Impacts of infinite bus voltage E and transmission line modulus of impedance |Z|

Fig. 18 shows the voltage stability boundaries with different E or |Z|. As increase of E or decrease of |Z|, the stability region expands evenly in all directions, and no intersection occurs. The MMVI, VMVI and ΔQc values with different E are given in Table 2, from which we can see that MMVI increases steadily with the change of E, while VMVI remains unchanged.

**Fig. 18.**Voltage stability boundaries with different E or |Z|

**Table 2.**MMVI, VMVI and ΔQc with different E

Table 3 and Table 4 show the results of MMVV, VMVV and ΔQc, VV (reactive margin to voltage violations) with different E and |Z| respectively. It can be seen that as increase of E or decrease of |Z|, the MMVV shows an increasing trend. The VMVV decreases as increase of E, and it remains unchanged as decrease of |Z|. This means that the high E value can reduce the risk of serious voltage violations, but low |Z | value would not have similar influence.

**Table 3.**MMVV, VMVV and ΔQc, VV, VV with different E

**Table 4.**MMVV, VMVV and ΔQc, VV with different |Z|

## 4.2 Impacts of wind power integration on voltage stability

### 4.2.1 Impacts of the wake effect

Fig. 19 displays the changes of MMVI and VMVI with respect to the wind speed. The initial load is set to be 0, and power factor of load interment is set to be 0.9. From Fig. 19 we can see that the shape of the MMVI curve is similar to the output curve shown in Fig. 2: the higher wind speed results in higher voltage stability margin, but the wake effect slows this tendency. In the range of medium wind speed (10-15m/s), the calculated MMVI with the wake effect taken into account is much lower than that without consideration of the wake effect. On the other hand, it can be found that the VMVI is not sensitive to wind speed variation no matter whether the wake effect is considered. The curves of MMVV and VMVV are highly similar to those of MMVI and VMVI and are not given here.

**Fig. 19.**Impacts of the wake effect on voltage stability

### 4.2.2 Impacts of the voltage control strategies based on reactive power capability of wind farm

In this section, the conventional constant power factor control strategy and the constant voltage control strategy based on the reactive power capability discussed in Section 2 are applied for comparisons. In the constant-voltage control mode, the power factor cosφwind=0.955.

Fig. 20 shows the voltage stability boundaries in different control modes and different wind speeds. When v=19m/s, the stability boundaries of the two control modes overlap each other, because the active power output reaches the maximum and there is no reactive power reserve left. As a result, the reactive output stays at the maximum no matter which control strategy is implemented. However, in a relatively low wind speed, the studies pertinent to different control strategies are significant. When v=8m/s, the stability boundary corresponding to constant voltage control is shifted towards the positive Q-axis direction compared with the one corresponding to constant power factor control. In the first quadrant, the shift considerably expands the stability region, and increases the stability margin. In fact, in the region nearby the load-increasing direction, the stability margin is close to that with the wind speed of 19m/s, i.e. the stability margin variation resulting from the wind speed variation is largely mitigated. The load margin gap between the high and low wind speeds situations is greatly narrowed by implementing constant voltage control strategy, so that the MMVI can maintain a relatively high value in a wide wind speed range. Fig. 22 demonstrates another advantage of applying the constant voltage control strategy: the VMVI drops greatly in the range of low wind speed, indicating that the risk of serious voltage instability problem is reduced in a large extent.

**Fig. 20.**Voltage stability boundaries in different control modes and different wind speeds

**Fig. 21.**MMVI with respect to wind speed in different control modes

**Fig. 22.**VMVI with respect to wind speed in different control modes

The voltage violation issue can be studied in a similar way. The constant voltage lines in different control modes and different wind speeds are shown in Fig. 23. The impact of constant voltage control strategy on constant voltage lines is even greater: the security region with constant voltage mode and a wind speed of 8m/s entirely covers the security region with constant power factor mode the speed of 19m/s.

**Fig. 23.**Constant voltage lines in different control modes and different wind speeds

This phenomenon is shown more clearly in Fig. 24. With the constant voltage control, the value of MMVV in the low wind speed even surpasses that in the high wind speed. Moreover, Fig. 25 shows that the system performs better in low wind speeds. In summary, compared with the conventional constant power factor control strategy, the constant voltage control strategy significantly improves power system’s voltage security in the range of low wind speed.

**Fig. 24.**MMVV with respect to wind speed in different control modes

**Fig. 25.**VMVV with respect to wind speed in different control modes

# 5. Conclusion

In this paper, an analytical approach for voltage stability evaluation incorporating wind power on PQ plane is proposed and studied in a single-load-infinite-bus system. A set of indices are defined with consideration of the load uncertainties. Based on the proposed method, the case studies are performed to investigate the impacts of power system key parameters and wind power integration on voltage stability. Furthermore, the impact analysis of wind farms with different control strategies is conducted. Some important and meaningful conclusions are obtained from the simulation results.

Future work involving the application to large-scale complex power system with remedy action strategies taken into account is under way by the authors to further enhance the model and method proposed.