Enhanced modelling of doubly fed induction generator in load flow analysis of distribution systems

Funding information Department of Science and Technology, Government of India, Grant/Award Numbers: DST/INSPIRE/04/2016/002381, ECR/2018/002565 Abstract This paper presents an accurate framework for incorporating the doubly fed induction generator based wind farms in load flow analysis of distribution systems. The proposed model for doubly fed induction generator is based on an accurate equivalent circuit and also takes into consideration the voltage dependent reactive power limits associated with the doubly fed induction generator. In addition, the developed approach is capable of handling a doubly fed induction generator operating in either maximum power point tracking or nonmaximum power point tracking modes. The proposed PQ model of doubly fed induction generator easily fits into the forward–backward technique that is typically employed for obtaining the load flow solution of distribution networks. Test results on 19-bus and 133bus systems are reported to indicate the effectiveness of the proposed approach.


INTRODUCTION
In recent years, wind energy-based power generating systems [referred to as wind turbine generating (WTG) units ] have evolved significantly due to technical advancements in power electronics and control of electric drives. During early installations, WTG were predominantly fixed speed units. On the other hand, present day installations are dominated by variable speed wind generators (VSWG). Two of the most popular configurations of VSWG are the doubly fed induction generator (DFIG) [1] and the generator with front-end converter. With the penetration levels reaching significant levels, it is essential to develop accurate models (for various WTG) and computationally efficient power system tools to understand their impact on various aspects related to power grid. The paper focuses on developing an accurate methodology for performing the load flow analysis of a distribution system with grid connected DFIG. Several models are proposed in the literature for incorporating wind generators in load flow solvers for both transmission and distribution systems. Most of the approaches proposed in the early days focused on fixed speed WTG and the heart of fixed speed WTG is the induction generator. corresponding bus (i.e. the node to which the WTG is connected) either as PX or PQ or PV bus. In PX-based modelling [2,3], the active power at the node corresponding to which WTG is connected is computed using wind velocity and is used as one of the characterising equations for solving the load flow problem. In addition, the magnetising reactance of induction generator is incorporated as an equivalent shunt reactance. Another approach that is commonly employed for modelling WTG is based on PQ modelling [4][5][6][7][8][9]. In this choice of modelling, turbine characteristics determine active power output and the reactive power is expressed as function of active power (which is obtained based on steady-state equivalent circuit of induction generator). In [4,[10][11][12], the authors have reported RX-based model for incorporating induction generators in load flow analysis wherein, an iterative algorithm is adopted.
Various models have also been reported in the literature for incorporating DFIG (variable speed type) based wind farms in load flow analysis. Similar to fixed speed WTG, DFIGbased WTG's can also be incorporated in load flow framework using PQ or PV or RX models [5,[13][14][15][16][17][18]. Divya et al. [5] proposed a model for incorporating DFIG-based WTG in distribution system load flow considering the limits associated with the converter. However, this model requires the machine to be modelled in DQ reference frame and also does not account for stator resistance (though its value is relatively small). Another approach for determining the steady-state performance of DFIG without using the DQ transformation is reported in [14] . For both these models, power is determined from turbine characteristics and depending upon the control specification, the bus to which the DFIG is connected is modelled as PQ or PV bus. An RX-based framework for handling DFIG in load flow studies of transmission networks is proposed in [18].

Shortcomings of the existing approaches
In general, most of the approaches adopt the steady-state equivalent circuit for obtaining the load flow solution with DFIG. However, in [19,20], the authors show that the existing equivalent circuit of the DFIG is inaccurate (when the DFIG is operating in super-synchronous mode) and subsequently propose an accurate model for DFIG. Based on the accurate model, an approach for solving the load flow problem with DFIGs is proposed in [17]. It is to be noted that the algorithm proposed in [17] is applicable for transmission systems and does not consider non-MPPT mode of operation. Furthermore, most of the approaches do not take into consideration the voltage dependent reactive power limits of the DFIG. This is of prime importance in distribution networks, since the voltage profile at the point of common coupling can be low under certain operating conditions (which would impact both the active and reactive power injection).

Contributions of this paper
This paper presents an approach for incorporating DFIGbased wind farms in load flow analysis of distribution systems. The proposed framework for obtaining the load flow solution is based on an accurate equivalent circuit and also takes into consideration the voltage dependent reactive power limits associated with the DFIG. The voltage dependent reactive power limits are developed by taking the effect of control action (i.e. due to back to back converter in the rotor) into consideration. In addition, the developed approach is capable of handling a DFIG operating in either maximum power point tracking (MPPT) or non-MPPT modes. The paper is organised as follows. Section 2 of the paper describes the generic framework for incorporating DFIG in distribution system load flow problem. Section 3 reports an approach to compute the state variables of the DFIG and Section 4 presents a way to handle the limits associated with the DFIG. The performance of the proposed approach is reported in Section 5 followed by conclusions in Section 6.

METHODOLOGY
This section describes the proposed methodology for solving the load flow problem of a distribution network in the presence of DFIG. Of the several approaches proposed in the literature for obtaining the load flow solution of distribution network, this work adopts the forward backward approach [21]. In the forward backward approach, the load flow solution is obtained in an iterative manner in which, each iteration involves two steps.
One of the pre-processing steps needed to apply the forward-backward approach is the network renumbering. The entire distribution network must be renumbered based on tree like parent chain relationship [21]. In the backward sweep, the current in each upstream branch is updated (starting from the remote load end) based on the load currents (which are computed by assuming the voltage at the node to be fixed during each iterative step or from on-measurem) and the downstream branch currents (i.e. application of KCL). In the forward sweep, starting from the primary substation/grid, the voltage at the downstream nodes (or child nodes) is updated based on the branch current obtained in backward sweep (i.e. application of KVL). This iterative process is repeated till the deviation in voltage (at nodes) obtained between consecutive iterations is very small.
The process of obtaining the load flow solution of a distribution network with DFIGs using forward backward approach involves four major steps. Initially, the process is started by assuming a flat start (with voltage magnitude angle as 1∠0 • ).

1.
Step-1: Node Renumbering: The first step in this method is node renumbering. This is done to generate a parent-child node relationship which is useful for backward and forward sweep [21]. The node renumbering procedure remains same even for distribution systems integrated with DFIG's. 2.
Step-2: Computation of load currents: In the first iteration, the load current is computed using a flat voltage profile. In all the subsequent iterations, the load current is computed based on the voltage profile obtained in the forward sweep phase of previous iteration. (a) Nodes at which load is present: The load current depends on the type of load existing at that particular node(or bus). Typically, the loads in distribution network are modelled either as constant impedance loads (or) constant current loads (or) constant power loads. Accordingly, the load current can be computed as  where Z L is the load impedance, I spe L is the specified load current, S L is the specified load power and V is the node voltage at the particular ith node. (b) Nodes at which DFIG is present: The net current injected by the DFIG to the grid can be computed as where S WG is the complex power injected by the DFIG, V PCC is the voltage on the grid side of the coupling transformer (i.e. the voltage at the node to which the DFIG is integrated) and the negative sign indicates that the current flows from the DFIG to the grid. The power injected by the DFIG depends on its mode of operation. Depending on the operating slip (s; which is determined by the wind velocity), the DFIG can operate in either sub-synchronous or super-synchronous mode. The power flow in the DFIG during the sub and super synchronous modes of operation is shown in Fig Neglecting losses (which is not adopted in this paper), P m (= P s − P r ) is the mechanical power input from turbine and is given by where is air density (kg/m 3 ), R is the radius swept by the rotor blades (m), v w is the wind velocity (m/s) and C p is performance coefficient which is a function of tip speed ratio ( ) and pitch angle( ) of the turbine. For MPPT mode of operation, the performance coefficient is set to C p,max which depends on the turbine characteristics. The negative sign (−P r ) in front of the rotor active power is because of the current convention employed in modelling the DFIG (i.e. the rotor current is assumed to be positive if its flowing into the rotor). During MPPT mode of operation, the control of rotor side converter (RSC) enables us to extract the maximum available mechanical power (P m ) and also allows us to inject the desired amount of reactive power from the stator (Q s ) [22]. The grid side converter (GSC) processes the same amount of active power as that of RSC in order to maintain a constant dc-link voltage. However, an additional amount of reactive power can be injected to the grid using GSC (Q GSC ) [22]. The network operator can choose (based on some optimization technique [23]) the amount of reactive power that DFIG needs to inject to the grid. However, the amount of reactive power that the DFIG can inject is limited based on the machine constraints which is discussed in detail in Section 3. Further, in certain scenarios, the DFIG is operated in non-MPPT mode [24], during which, the amount of current injected by the DFIG can be computed as follows: where S s WG is specified power at the DFIG integrated bus. It must be noted that the amount of current that is injected to the grid is strongly dependent on the machine and converter and as a result, the injected current must be updated taking the limits into consideration. 3.
Step 3: Backward sweep: After computing the load currents, the branch currents are computed using parent-child node relationship (obtained using node renumbering) and the process starts from the branch in downstream and terminates with the branch at the source. The current in the upstream branch can be computed using the current in the downstream branches and the load currents as where I i j is the branch current between ith (parent) and j th (child) node, I L j is the load current at the child node, Y j is the shunt admittance of the child node and I bk indicates the current from downstream branch. 4.
Step 4: Forward sweep: In this step, the node voltages are computed using the branch currents obtained using backward sweep and the process moves from upstream to downstream (source to load end). The general expression 4: Check for violations of current limits and adjust the power output of DFIG accordingly (as discussed in Section 4).
5: Recompute DFIG currents with the adjusted outputs (as discussed in Section 4).
6: Compute the branch currents using backward approach.
7: Compute node voltage using forward approach.

8: end while
for computing the node voltage is where V j is voltage of the child node, V i is the parent node voltage and Z i j is the impedance between parent and child nodes. The steps involved in forward-backward approach with the inclusion of DFIG in the distribution systems is presented in Algorithm 1. In order to avoid convergence issues, the limits checking of DFIG currents is typically avoided in the first few iterations.

COMPUTATION OF STATE VARIABLES OF DFIG
The node to which the DFIG is integrated can usually be modelled as PV bus or PQ bus depending on the control objective of the back to back converter [17] . Typically, in DFIGs integrated to distribution systems, the control specification to the back to back converter is the desired reactive power i.e. the DFIG is operated in constant PQ mode. As a result, the node to which the DFIG is connected is modelled as PQ node, thus making it feasible to handle it in the manner specified in Section 2.
The amount of active and reactive power injected by the DFIG is limited by the constraints associated with the converter and the induction machine. Hence, it is necessary to check if these limits are violated. To achieve this, it is necessary to compute the state variables of the DFIG.
The state variables of the DFIG can be computed using the steady-state equivalent circuit shown in Figure 2. In this paper, the state variables of the DFIG when operating in supersynchronous mode are computed using an accurate equivalent circuit as shown in Figure 2(b) [19]. The only approach to be proposed in the literature based on accurate equivalent circuit is reported in [17]. However, as opposed to the modelling carried out in [17], this work considers the DFIG to operate in non-MPPT mode as well. As a result, the state variables are computed based on a set of 7 non-linear equations as opposed to 6 in [17]. The complete modelling is outlined in this paper for the sake of completeness.
For a given operating condition, the state variables of DFIG are obtained by solving the loop equations of the DFIG (i.e. application of KVL to stator and rotor) along with the power balance equations.
The set of equations that describe the behaviour of DFIG during sub [based on Figure 2 where ⃗ v s = V s ∠ s , ⃗ i s = I s ∠ s represent the stator voltages and currents in phasor form, respectively, while ⃗ v r = V r ∠ r , ⃗ i r = I r ∠ r represent the rotor voltages and currents phasor form, respectively. The power balance equations of the DFIG (based on the chosen current convention) are In MPPT mode, the wind velocity determines the P m (under this condition c p = c max p ) and in non-MPPT mode of operation or in scenarios where the DFIG needs to inject more reactive power (discussed in Section 4), the grid operator determines the amount of P m that must be extracted from the turbine (under this scenario, c p = P spe m P max m ). Under non-MPPT mode, the extracted power from the turbine is contained by changing the pitch angle of the turbine. In either case, P m forms the input in the power balance Equation (9) and Q wg is the reactive power specified by the grid operator.
For a grid connected DFIG, the stator voltage is determined by the grid and as a result, the set of unknown variables (or state variables) are i s ∠ s , i r ∠ r and v r ∠ r . The set of equations given by Equation (8) are essentially complex equations and the non-linear equations required to obtain the state variables can be obtained by equating the real and imaginary parts. The resulting set of equations that need to be solved for determining the state variables of DFIG operating in sub-synchronous mode are V s cos( s ) + I s R s cos( s ) − I s X s sin( s ) + I r X m sin( r ) = 0 V s sin( s ) + I s R s sin( s ) + I s X s cos( s ) − I r X m cos( r ) = 0 V r cos( r ) − I r R r cos( r ) + sI r X r sin( r ) − sI s X m sin( s ) = 0 V r sin( r ) − I r R r sin( r ) − sI r X r cos( r ) + sI s X m cos( s ) = 0 Similarly, the set of non-linear equations that need to be solved for computing state variables of DFIG operating in super-synchronous mode are V s cos( s ) + I s R s cos( s ) − I s X s sin( s ) − I r X m sin( r ) = 0 V s sin( s ) + I s R s sin( s ) + I s X s cos( s ) − I r X m cos( r ) = 0 V r cos( r ) − I r R r cos( r ) − sI r X r sin( r ) − sI s X m sin( s ) = 0 V r sin( r ) − I r R r sin( r ) + sI r X r cos( r ) − sI s X m cos( s ) = 0 where where k indicates the present iteration count and J DFIG is Jacobin matrix. This is repeated till the convergence criterion is satisfied. In the power balance equations Q wg is the reactive power processed by the RSC. In addition, some amount of reactive power can be processed by the DFIG. However, typically, the GSC is operated in UPF mode and as a result, in this work, we assume the GSC of the DFIG to be operated in UPF mode.

HANDLING OF DFIG LIMITS
In general, the parameters that impact the amount of complex power that can be injected by the DFIG are (a) maximum stator current (I max s ), (b) maximum rotor current (I max r ), (c) mechanical power and (d) the rotor voltage [25]. For the typical operating conditions, the limit due to rotor voltage may be neglected. The typical PQ limits of a 1.5 MW DFIG due to these constraints is shown in Figure 3. As highlighted in [17], if the DFIG is injecting reactive power to the grid, the maximum amount of reactive power, the maximum amount of complex power injected into the grid is limited by the rotor current limit. On the other hand, if the DFIG is absorbing the reactive power from the grid, then the maximum amount of complex power is limited by the stator current limit.
Since the PQ limits of DFIG are dependent on stator voltage, the active and reactive power scheduling of DFIG carried out for one operating point may not be within the operating limits for a different operating condition. Hence, it is necessary to account for these limits and adjust the output of the DFIG (either active or reactive power) in case of any violation.
Reduction of reactive power output: If the DFIG is injecting reactive power to the grid and the rotor current at any point in the iterative process exceeds the maximum value (I max r ), then the reactive power output of the DFIG is set to (14) On the other hand, if the DFIG is absorbing reactive power from the grid and the stator current at any point in the iterative process exceeds the maximum value (I max s ), then the reactive power output of the DFIG is set to (15) Reduction of active power output: In certain scenarios, the active power output of DFIG is changed keeping the reactive power output unaltered. If the DFIG is injecting reactive power to the grid and if the rotor current at any point in the iterative process exceeds the maximum value (I max r ), then the active power output of the DFIG (by maintaining the same reactive power) is set to Equation (16) If the DFIG is absorbing reactive power from the grid and if the stator current at any point in the iterative process exceeds the maximum value (I max s ), then the active power output of the DFIG is set to Equation (17) .

RESULTS AND DISCUSSIONS
In this section, the performance of the proposed algorithm is demonstrated on two practical distribution systems of Indian Southern Grid namely 19-bus system and the 133-bus system [26]. A program based on the proposed algorithm is developed in Python environment. For the purpose of analysis, a DFIG of 1.5 MW rating is assumed to be integrated to the two systems. As mentioned in Section 3, Q spec wg is typically specified by the grid operator. In this paper, the case studies are chosen such that the effectiveness of the proposed algorithm (in handling various limits of DFIG) can be demonstrated.

19-bus system
To start with, the results on 19-bus system is presented [26].
The single line diagram of the system is shown in Figure 4. The DFIG is placed at 19th bus (integrated via a transformer) and the behaviour of system in sub-synchronous and supersynchronous modes of operation is presented.
1. Case A (v w = 9 m/s): Under this scenario, the mechanical power that can be extracted from the turbine is 0.582 p.u.
(on machine base of 1.5 MVA). The DFIG is assumed to be operating in MPPT mode with the reactive power specified to be 0.3 p.u. The typical voltage profile of the system under this operating condition is shown in Figure 5(a). For this operating scenario, it can be observed that the stator and rotor currents are well within the limits [see Figure 5(b)], as a result of which the DFIG is capable of processing the specified reactive power with the given mechanical power input without any reduction. 2. Case B (v w = 10 m/s): When the wind velocity is at 10 m/s, the DFIG continues to operate in sub-synchronous mode with the maximum possible power that can be extracted from the turbine is 0.79 p.u. on machine base. The reactive power specification is set to 0.3 p.u.. Unlike the scenario in Case A, it is observed that the rotor current limit is hit during the iterative process due to which either the active and reactive power injected into the grid needs to be reduced. Results corresponding to both the scenarios are reported.  analysis is carried out with the existing approach [14]. As highlighted earlier, existing approaches do not account for voltage dependent reactive power limits (arising due to current limits of the stator and rotor) which will result in an inaccurate load flow solution.
To illustrate this, the stator and rotor currents obtained at the end of load flow solution (corresponding to the current scenario) obtained using the existing approach is shown in Figure 6(c). Observation into Figure 6(c) indicates that although the stator currents are well within the limits, the rotor currents are much higher than the rated value. In practice, this will not be possible, since the current limits of the converter control loop do not allow the rotor current limits to be violated. As a result, it can be concluded that the results obtained using the existing approach are inaccurate.  Figure 7(c) indicates that the voltage profile is slightly better in the operating scenario where reactive power reduction is employed. In transmission system there is strong coupling between Q and | V | (due to large x r ) giving an impetus for reducing the active power output for a gain in system voltage profile (in most of the scenarios). However, in distribution systems, due to a low x r ratio, it will not be possible to conclude the superiority of one strategy (i.e. active or reactive power reduction) over other.  power that can be extracted from the turbine is 1 p.u. and as a result the reactive power specification is set to 0 p.u. (i.e. UPF mode of operation). The resulting voltage profile of the network is shown in Figure 8(a). It is interesting to note that, even with zero reactive power injection to the grid, the stator and rotor current limits are violated [as shown in Figure 8

133-bus system
The second system considered for case study is a 133-node distribution network [26]. The system comprises of two DFIGs placed at 72nd(1st) and 95th(2nd) bus (integrated via a transformer). Due to space constraints, results pertaining to a single case study (where the active or reactive power injected to the grid reduces) is presented. For this case study, the wind veloci-ties at 1st and 2nd DFIG locations are considered to be 9 m/s and 9.5 m/s, respectively, (wherein the DFIG's are operating in sub-synchronous mode). The maximum possible active power outputs from DFIGs corresponding to the chosen wind velocities are P m1 = 0.58 p.u. and P m2 = 0.68 p.u. respectively. The reactive power specifications to the control loops of the DFIG are chosen to be Q wg1 = 0.3 p.u. and Q wg2 = 0.4 p.u., respectively.
1. Scenario 1 (Reactive Power reduction): For the operating point under consideration, the rotor current limit is reached in generator connected at node 95 [shown in Figure 9(a)]. As a result, the reactive power injected by the DFIG reduces from 0.4 p.u. to 0.36 p.u. [shown in Figure 9(b)]. On the other hand, the DFIG at bus 72 is able to inject the desired reactive power. Figure 9(c) indicates the stator and rotor currents in both the DFIGs at the load flow when the solution is obtained using the existing approach [14]. It can be observed from Figure 9(c), that the rotor current of DFIG located at node 95 [indicated by I r2 in Figure 9(c)] is well above the limits. This indicates the inability of existing approach to handle  Figure 10(c). Even for this system, it can be observed that voltage profile is better when the DFIG operates in MPPT mode (i.e. active power is given priority in case of any limit violation).

CONCLUSION
In this paper, an accurate framework to incorporate a DFIG in load flow model of distribution system is presented. The pro-posed model is capable of taking into consideration various limits associated with the converter and can also handle non-MPPT mode of operation. Case studies reported on 19-bus system and 133-bus system indicate the effectiveness of the proposed algorithm. Owing to a relatively large r x ratio of distribution systems, it is observed that a better voltage profile is obtained by giving priority to active power output of DFIG instead of reactive power in case of any current limit violations. The proposed approach considers a balanced distribution network making it feasible to analyse the network as a single phase equivalent. Our future work in this direction is to propose a model for DFIG for unbalanced distribution networks.