Dynamic equivalent of wind farm model for power system stability studies

ABSTRACT A new dynamic equivalencing method for stability assessment of a grid-integrated wind farm is proposed in this article. The accuracy of the method is validated for a 34-bus system with 28-unit wind farm connected to Indian utility system. This wind farm consists of several wind turbines of two different ratings. The electrical parameters of the equivalent generator are derived from the mathematical model of the squirrel-cage induction generator. The parameters of the equivalent wind-turbine generator are optimized to yield minimum deviation from the detailed system response using genetic algorithm. The small-signal and transient stability responses of the study system with detail wind farm and equivalent model are simulated using MATLAB. Equivalent model eigenvalues are compared to the centre of inertia based detailed system eigenvalue. In addition, the computed eigenvalues and time-domain responses of the proposed equivalent model, detailed wind farm are compared against weighted model proposed earlier. In most of the investigated cases, the average error of dynamic responses between the proposed equivalent and detailed models are less when compared to weighted model. Thus, the large-signal responses of the proposed equivalent model show superior agreement with detailed system response.


Introduction
In developing countries like India, the steady state and dynamic operation of the power system are mainly influenced by increased installed capacity of wind generation. Moreover, the power evacuation from the large-scale wind farm (WF) is limited by the low short circuit power of the grid. In the immediate assessment of the stability simulation of such grid integrated WF, an accurate single-equivalent model of the WF is significantly needed. Such an equivalent model is useful in carrying out steady state and dynamic simulation studies that form an integral part of power system planning. Grid code standards and requirements for integrating the large-scale WF into bulk transmission system are discussed in [1].
The qualitative study of simple power system with large-scale wind power penetration is investigated in [2]. It was found that the impact of wind power mainly contributes damping. The power oscillation damping effect is improved particularly when squirrel-cage induction generators (SCIGs) are used.
Most of the WFs in India guise tree-like structure in which the wind-turbine generators (WTGs) of different ratings are mixed. In such case, simple aggregation of the WF yields inaccurate solution [3,4]. Identical WTGs with same operational points of WTG in a WF is assumed. The effect of equivalent network impedance is not included in the equivalent WTG. Multi-turbine equivalent model suggested for an irregular wind distribution over the large area of the WF. Aggregation technique used to develop an equivalent model for the fixed speed WTG-based WF developed in [5] includes the swing equations, the impedance of the generator transformer and generator impedances. In addition, the injected real power by the aggregated model has more error compared to detail model under transient condition. This is due to that fact the aggregated model does not accounting the feeder impedance between WTG to point of common coupling (PCC) point.
A method of aggregating fixed speed-wind-turbine generators (FSWTGs) based WF is studied in [6] and yields two equivalent WTGs. In the case of uniform wind distribution in the WF, a simple equivalent model is used. In addition a refined model is also suggested to take into account non-uniform wind speed. At PCC, the responses of variable single equivalent with compensating capacitor model approximately matches with the detailed system. Single-equivalent model of WF that consists of fixed speed WTGs is developed in [7] and the simulated results are validated against the field measurements. Short-and long-time frame simulations have been taken for comparison and the error between aggregated and detailed models becomes larger if all the WTGs are operated at different operating conditions. In this equivalent model, the mechanical input power is assumed to be constant. Also, mechanical parameters used in the simulation are inaccurate [7]. However, in the case of non-uniform wind velocity in large-scale WF, some of the WTGs are tripped due to high/low velocity of wind and the single-equivalent model cannot be used in such event. Analytical approach used in the article [8] is to derive the equivalent (collector) impedance by equating the losses within the detailed WF branch to the total losses of equivalent line/cable. However, all turbines are assumed to be operating in the same point. An aggregated model was developed for the mixed WF with SCIG and doubly fed induction generator (DFIG) as reported in [9]. However, the method for obtaining the equivalent model generator parameters is not straightforward. According to this, genetic algorithm-based optimization is used to minimize the error between the responses of equivalent single generator model and detailed WF. Also, the response of equivalent model is compared with weighted equivalent model reported in the earlier literature [10].
The dynamic equivalent model is developed in order to improve the operating efficiency of WF [11]. Using phasor measurement unit (PMU), the data are obtained from the WF and equivalent model of the WF is derived from parameter identification technique. However, the equivalent model parameters are optimized based on stochastic approximation. An equivalent WF model was developed from synchro-phasor measurement data by the use of parameter identification technique discussed in [12]. The improved genetic algorithm (GA) was employed to develop the equivalent WF for preserving the basic structure, characteristics and control patterns of actual WF. Using the deterministic data, the dynamic multi-machine equivalent modelling of WF is developed [13].
In summary, the methods in [3][4][5][6][7][8][9] suggest an aggregated model of similar type with identical rating. Also the equivalent generator parameter just used the re-scaled value of WTGs in the WF. The aggregation method proposed here yields the equivalent system dynamic parameters. It can reduce the system order and simulation time considerably. The method proposed here is based on aggregation technique adopted in [14]. In this reference, several induction motors in an industrial power plant are aggregated to yield a single equivalent which is valid for dynamic conditions. Here, this method is extended to WTG system using SCIG which involves Kron's reduction technique.
The major advantages of this proposed method are as follows: (1) WF of any network topology can be aggregated using this method. Most of the methods fail, if meshed topology or any complex network is opted for WF.
(2) WF consisting of different ratings of generators can be aggregated as single equivalent. In other aggregation methods, if the WF consists of different ratings, multiple clusters are formed and it is represented as multiple equivalents. But simple aggregation of WF is applicable only for identical WTGs. (3) Equivalent parameters are accurately obtained from the available parameters of detailed WF. (4) The electrical parameters are obtained such that it is operated in steady state; the generation produced by the aggregate model is identical to that of the detail WF model. (5) The mechanical parameters of the equivalent WF model are optimized to match the dynamic performance of the model with that of the detailed WF. GA is used in this optimization process.
The paper is organized as follows. Section 2 briefly describes the modelling of the synchronous generator (SG) and the WTG components. In addition, the modelling of two-mass model and wake effect are included. A brief description of the WF considered for investigation is presented in Section 3. The aggregation process to yield the single equivalent based on weighted model method and proposed method are explained in Section 4 and Section 5, respectively. Eigenvalues and time-domain responses of the study system with detail WF representation, proposed equivalent and weighted model reported in [10] are compared in Section 6. Finally, Section 7 concludes the paper.

Modelling of synchronous generator and WTG components
This section presents the modelling of SG and WTG. The modelling of the SG can be found in detail in [15][16][17][18]. The major components of the FSWTG are wind turbine, gear box and SCIG. These are briefly explained in this section. Furthermore, the modelling of wind speed and wake effect is also presented.

Modelling of synchronous generator
The dynamic modelling of SG is cast in direct and quadrature coordinates rotating at synchronous speed [15,16]. The classical model of the SGs is taken for study, so rotor angle and rotor speed of SG are state variables: The electric torque of the synchronous machine is expressed by The stator algebraic equations are represented as follows: For i = 1,… nsg, where nsg is number of SGs

Modelling of wind turbine
An aerodynamic rotor model gives the relationship between the mechanical power of the wind turbine to the rotational speed. The power coefficient (C p ) is decided by the tip speed ratio (λ) and pitch angle (f). In this paper, the detailed WF has 200 kW and 250 kW FSWTGs with passive stall technique for controlling the speed of the turbine. Hence, the pitch angle control is not included in this model. The tip-speed ratio is given by The wind turbine captures the wind energy through the blades, and delivers mechanical power (P wt ) to the shaft. The power extracted from the wind varies as the cube of wind speed and is expressed as The mechanical torque (T m ) of wind turbine is given by The power coefficient C p is The C p versus λ and power curve of a wind turbine are provided by the manufacturer. The values of wind wheel coefficients, C 1 to C 6 , the exponent x and pitch angle f are given in the Appendix.

Modelling of squirrel-cage induction generator
SCIG is employed for wind power generation because of its simple construction and ruggedness. The dynamic modelling of induction generator is cast in direct and quadrature coordinates rotating at synchronous speed [15,16]. Three-phase stator and rotor voltage equations are referred to synchronously rotating reference frame with quadrature axis (q-axis) leading direct axis (d-axis) by 90 . Generator convention is used, i.e. stator current and rotor current are considered to be positive when they are leaving and entering the machine, respectively.
The stator dq voltage equations are shown in Figure 1, the two-axis diagram [17,18] forms relationship between terminal voltage, speed emf and resistive drop. The diagram for the rotor is obtained by replacing the stator quantities by corresponding rotor quantities and v s by v sv r.
The per unit stator and rotor dq-axis voltage equations of SCIG are as follows: v ds ¼ ÀR s i ds þ v s c qs À pc ds (10) v qs ¼ ÀR s i qs À v s c ds À pc qs (11) v dr ¼ ÀR r i dr þ sv s c qr À pc dr (12) v qr ¼ ÀR r i qr À sv s c dr À pc qr: The per unit stator and rotor dq-axis flux linkage equations are as follows: The third-order model for the induction generator is considered with the following assumptions.
(1) The stator transients are negligible, which implies pc ds and pc qs terms are zero in Equations (10) and (11). (2) The rotor is short circuited which implies that v dr and v qr are zero in Equations (12) and (13).
The transient reactance (X 0 s ) is given by The differential equations of the third-order model of SCIG are given by where the rotor open circuit time constant is T o 0 = L r /R r . The expression for electromagnetic torque is

Model of drive train
The drive train of the wind turbine is represented by a two-mass model as shown in Figure 2 and is described by the following equations: (for more details see [19]): where

Modelling of wind speed
Generally, wind speed is intermittent and stochastic in nature. The mechanical torque produced by WTG is directly proportional to the wind speed. Hence, the suitable wind speed model must be taken into account to simulate WTG dynamics. The four-component wind speed model is chosen in this paper [4] and defined as

Modelling of wake effect
A simple wake effect model is illustrated as shown in Figure 3 [20,21]. According to this model given by Equation (29) the wind speed seen by the turbine rotor 2 in row-2 is decreased, thereby affecting the power generation of row-2 WTGs installed within the wake region as shown in Figure 3:

Detailed WF system
The WF considered for investigation is located in Kayathar, Tamil Nadu state, India, with a total installed capacity of 5.85 MW. The WF is connected at bus 8, PCC of 34-bus system shown in Figure 4. It can be seen that the detail WF is given in Figure 5 which shows 28 numbers of WTGs arranged in 3 rows. The generators in the first row encounter higher wind speed and the speed decreases as we go down the rows, taking wake effect into consideration. The 200-kW wind generator is considered in second row except one which is rated of 250 kW. The rating of each SG is 250 MVA and load 1 is connected at bus 5 and load 2 at bus 6 with each capacity of real power 100 MW and reactive power of 40 MVAr. Grid is represented as SG in dynamic simulation. SG data are found in [18]. The data of study system are given in the Appendix (Tables A.1 The data for the equivalent system model are derived from the detailed representation of the WF. At bus 8, the entire WF is represented by the single-equivalent model as shown in Figure 6.

Weighted model equivalent of wind farm
Weighted-model-based equivalent method gives the simplest representation of aggregated WF model with acceptable precision in power system stability studies reported in [10]. Wind turbine facing same wind velocity can be grouped together and it can be represented as a single aggregate wind turbine model or an equivalent generator model.
The salient points that are discussed in the weighted model in reference [10] are as follows: (1) The first-, third-and fifth-order WTG models in the WF are compared (Electromagnetic transient) for the EMT type simulation using software packages PSS/E and PSCAD/EMTDC. (2) In this method, all the WTGs in the WF are considered as identically rated. According to this method, the equivalent generator parameters are obtained based on scaled down values of WTG in a WF. to aggregation of type-4 WTG (PMSG-based fully rated converter WTG) based WF investigated in [22] and the model used includes both stability, EMT type simulation studies.
However, the voltage drop across the line is depending on the real power flow as well as reactive power flow. Also, the equivalent weighted model's line impedance calculation is not straightforward because the WTG connection in detailed WF should be known prior whether it is connected in series or parallel. The wind speed variation in the WF is not considered. Also, wake effect model is not included in the equivalent WTG. In the proposed method the equivalent line impedance of the WF retained in the equivalent generator parameters through Kron's reduction. Furthermore, in this work, the stability investigation of detailed and equivalent WTG models is performed with wake effect. The methodology for obtaining the proposed dynamic equivalent model of WF is discussed in the next section.

Proposed dynamic equivalent wind farm model
Dynamic equivalent of WF is constructed in two-fold. First, the equivalent wind turbine model is developed using Equations (30)-(33). Second, the equivalent generator is yielded through modelling of SCIG.

Wind turbine aggregation
The power rating of the equivalent wind turbine is obtained by adding individual machine wind turbine ratings. The equivalent wind speed (U we ) is given by The equivalent generator mechanical power P WTe is given by where K WTe is equivalent generator power constant given by The equivalent WTG's compensating capacitor (C e ), turbine inertia constant (H te ), generator inertia constant (H ge ), stiffness constant (K she ) and damping constant (D she ) are calculated as follows: The compensating capacitor of the equivalent generator C e is included in the admittance matrix.

Equivalent generator model electrical parameters
WTGs are represented by fifth-order model using Equations (21), (22) and Equations (24)-(27) for detailed WF system, as shown in Figure 5. Bus admittance matrix is constructed for the study system shown in Figure 4. Kron's reduction is applied to admittance matrix ½Y to get the equivalent generator transient reactance.
The individual generator transformers and interconnecting lines of the detailed WF are included in the admittance matrix ½Y given by The transient admittance of the individual generators is absorbed in the diagonals of ½Y which is reduced to yield the admittance 1 R s þjX 0 s of the equivalent generator. Thus, we have The SCIG model is used to derive the other parameters of the equivalent generator.
At the PCC (bus 8) in Figure 4, the real power generated and reactive power absorbed by the detail WF at the rated wind speed are represented as P and Q. Using P and Q, the equivalent generator power factor angle, u can be obtained as For the equivalent SCIG, voltage behind transient reactance is given by 5.3 Procedure to obtain stator reactance of equivalent generator (X s ) The steps for arriving at the stator reactance of equivalent generator are given as follows: Step 1: Get the phasor form of the voltage pE 0 m using Equations (21) and (22). It is given by Step 2: For steady-state condition pE 0 m = 0. AssumeX r ffi X s and setting magnitude of stator current to unity in Equation (37), we can show that E 0 qs E 0 ds ¼ R r cos u À sX r sin u R r sin u þ sX r cos u : Step 3: Setting V s = 1 p.u. and I s = cos u À j sin u in Equation (38) yields Step 4: Solving Equations (40) and (41) yield Equation (42). The equivalent SCIG's stator reactance (X s ) is given by To calculate the stator reactance of the equivalent generator the value of (R r =s) has to be known. The procedure for obtaining (R r =s) is described below.

Finding the ratio of rotor resistance to full load slip of equivalent generator
The major steps in arriving at the ratio of rotor resistance to full load slip of equivalent generator are as follows: Step 1: Stator current of induction generator from Equation (37) is given by Step 2: Substituting Equation (43) in Equation (38) and considering steady-state condition with 1.0 p.u voltage yield .
Step 3: Rationalising Equation (44) and using Equation (39), we can obtain Step 4: It is assumed that the equivalent generator is operating at full load with rated speed delivering rated power at unity current magnitude. Solving Equation (45) for R r =s yields Equation (46).
The deduced expression of ratio of rotor resistance to full load slip of the equivalent generator is given by When the equivalent generator is operated at full load, the magnitude of slip j s j is close to R r and hence j R r s j becomes unity. Equations (42) and (46) are directly solved for X s and R r =s .

Determination of operating slip and finding
the magnetizing and leakage reactance of the equivalent generator The aggregated mechanical power is obtained by using Equation (31). The terminal voltage at bus-8 for the equivalent generator is calculated from the steady-state operating conditions of the detailed WF. Therefore, by knowing the mechanical power supplied to the equivalent generator and terminal voltage of the equivalent generator, the slip (s) is calculated iteratively from the simple quadratic equation as discussed in [23]. The magnetizing reactance (X m ) for the equivalent generator is given by The leakage reactance (X l ) of the equivalent generator is Hence, the known equivalent generator electrical parameters R s , R r =s, X s , X m and X l are sufficient to mimic the response of the detailed system. The optimization of the equivalent generator parameters is explained in the next section.

Optimization of single-equivalent generator parameters
Development of single-equivalent WF model and its parameters optimization were reported in the earlier literature [24]. The rotor circuit time constant, a moment of inertia and transient reactance of the equivalent generator are taken as control variables for optimization. However, simple WF model is investigated with lumped representation turbine and generator of FSWTG. Also, the SG dynamic influence on WTG is not included in the model. An optimization method has been used to minimize the error between responses of the detailed WF and single-equivalent model for a dynamic conditions and normal operating condition. Thus, the single-equivalent generator parameters are optimized. The optimization problem can be solved by using the GA [25]. In this paper, the following control variables are taken for equivalent generator, rotor resistance, inertia of wind turbine, inertia of generator, damping constant of WTG and shaft stiffness. Sum-squared error deviation (SSED) is taken as the main objective function (f) with input variables of voltage, real power and reactive power of detailed WF and equivalent generator of WF at the PCC (bus 8) as shown in Figure 4 and Figure 6, respectively.
The optimization problem can be stated as Subject to equality constraints: For the equivalent generator internal bus, the electrical output mismatch equation is For the equivalent generator electrical power output and detail system WTG power output mismatch equation is P WTe À X nt k¼1 P WTk ðU w ; C p Þ ¼ 0; for equivalent generator terminal bus, the real and reactive power mismatch equations are D shemin < D she < D shemax ; where i is time step in the non-linear time-domain simulation and N is the maximum number of time steps. The V L D ; P L D ; Q L D and V E ; P E ; Q E are variables of the detailed system and equivalent system, respectively. For the each time step, the difference between arrived at detail and equivalent system variables error are squared and added to the fitness function. The objective function (49) is minimized by satisfying the equality constraints given by Equations (50)-(55) for non-linear timedomain simulation and the steady-state operating equality constraints given by Equations (53) and (56) of an equivalent generator. The GA-based optimization algorithm is applied to the equivalent systems. The maximum and minimum values of control variables and the parameters of GA are listed in Table A.4. The equivalent generator's rotor resistance (R r ), inertia of wind turbine (H te ), generator (H ge ), damping co-efficient (D she ) and stiffness coefficient (K she ) are chosen as decision variables of the optimization problem. The optimized solution yields the best equivalent model parameters. Optimized equivalent generator parameters are summarized in Table 1.

System study and results
The proposed equivalent model is validated for small and large disturbances by conducting small-signal and large-signal analysis for the study system shown in Figure 4.

Small-signal model of the study system and equivalent models
To perform the small-signal stability analysis of the study system shown in Figure 4, the small-signal linear model is derived. More details about the development of state space model of a system are found in reference [26]. In this paper, the classical model of SG and fifthorder model of FSWTG are used for study. The state variables of the complete system are given by where subscript "sg" denotes SGs variables and superscript "i" refers to induction generator variables.
The state matrix of the study system is obtained by linearizing differential equations and eliminating the algebraic variables as discussed in [26] from Equations (1), (2), Equations (21), (22) and Equations (24)(25)(26). Generalized model has been derived analytically The eigenvalues of system are determined from the system matrix A S and presented in Tables 2 and 3. The eigenvalues of the study system with detail WF representation, weighted model and proposed equivalent model are presented in Tables 2 and 3.
From Table 2, it can be observed that rotor modes of synchronous machines with WF represented by detailed model are close agreement with proposed equivalent mode than weighted average model. The rotor modes, mechanical mode and electrical modes of the WF models represented in Table 3, eigenvalues of proposed equivalent are matching closely with the detailed WF than weighted average-based equivalent WF model. Therefore, the dynamic equivalent of the WF by proposed equivalent WF system preserves accurately for the given detail system.

Large-signal analysis of the study system and equivalent models
The study system with detailed WF, the WF replaced by weighted model and equivalent systems programmed for transient stability analysis using MAT-LAB [27] and their dynamic responses are compared.  magnitude at bus 8 for the study system with detailed WF representation and WF replaced by equivalent models are shown in Figure 7.

Case 2: Rotor speed instability
In this case, the damping of all three SGs is taken as zero. A three phase to ground fault is simulated at bus 8 at 1 sec and removed at 1.1 sec. Due to insufficient damping of SGs and severe fault leads to growing of oscillations in rotor speed of WTG shown in Figure 8(g)-(h). The impact of the growing of oscillations is reflected on real power, reactive power, and voltage at PCC as shown in Figure 8(a)-(f). This is referred as "rotor speed instability" of WF reported in [28]. The detailed, equivalent model dynamic responses are obtained and compared.

Case 3: Normal operation of the wind farm
Under normal operation, the dynamic resulting responses of detailed WF, equivalent by weighted model and proposed are shown in Figure 9. The time series of measured wind speed is considered for this simulation case as shown in Figure 10, applied to detailed WF. The wind speed for proposed equivalent model is calculated as per Equation (30).

Average error between study system with detailed WF and equivalent models
The average error of the responses has been computed for the all the above investigated cases for the detailed system and equivalent models. The average error is expressed as From Table 4, it can be seen that the average deviations in the active power, voltage, rotor angle of SG-2, wind turbine speed, SCIG's rotor speed and it is torsional angle responses between the proposed equivalent and detailed models are less than 1% except for rotor speed instability case, whereas for weighted model the deviations are more. Also proposed equivalent model responses at the PCC are in close agreement with the detailed model responses as shown in  It is inferred from Table 5 that time required to simulate equivalent WF model is very less compared to the study system with detailed WF representation. Hence, the proposed and weighted model equivalent of WF system has significantly reduced the simulation time compared to detailed WF for the time-domain analysis of power system. Processor -Intel ® Pentium ® -CPU B950 @ 2.10 GHz, 2 GB RAM is used for this simulation.

Comments on equivalent wind farm model
The WTG in a WF may be arranged to be in series or parallel or mixed [29]. In any case, the proposed method can be applied. The impact of the wake effect is observed to be less for the equivalent WF model. This is expected since the single-equivalent generator is compact, whereas the WTGs in the detailed WF are distributed. The eigenvalues of the study system with detail WF representation, WF replaced by proposed equivalent model and weighted model are presented in Tables 2 and 3. The damping and natural frequency of oscillation of proposed equivalent system are very close to detailed WF system when compare to weighted model. In most of the investigated cases, the average error of dynamic responses between the proposed equivalent and detailed models are less when compared to weighted model as shown in Table 4. The equivalent system responses exactly match with detailed system pre-fault conditions, but the response of the weighted model has more deviation during post-fault conditions. But the average error of the dynamic responses of the detailed system and proposed equivalent system are less with acceptable tolerance.

Conclusion
A method of aggregation of a complex WF that yields a single equivalent is proposed in this paper. The performance of the proposed equivalent is superior to the weighted model equivalent proposed earlier. The dynamic performance of the aggregated equivalent is tested against that of the complex WF. The test involves small-and large-signal disturbances that trigger the electromechanical dynamics. The equivalent WF is obtained by proposed method is also investigated for rotor speed instability case. In addition, the computed eigenvalues and time-domain responses are compared against weighted model. The eigenvalues obtained from the small-signal study of the proposed equivalent system are very close to detailed WF system than weighted model. Also the timedomain response from the large-signal study reveals that there is a good agreement between the proposed equivalent model and the detailed WF. In some cases, the average error between the detailed and proposed models response is less than 1%. The large-signal responses of the proposed equivalent show significantly superior agreement with detailed system dynamic response compared to weighted model in most of the investigated cases. Future scope of this proposed technique extended to deduce dynamic equivalent of WF with different WTG technologies with suitable assumptions.

Disclosure statement
No potential conflict of interest was reported by the authors.  wind speed at first and second wind turbine in m/sec, respectively U we, U wi wind speed at equivalent generator and i th WTG in m/sec, respectively v ds , i ds d-axis stator voltage and current in p. u., respectively v qs , i qs q-axis stator voltage and current in p. u., respectively v dr , i dr d-axis rotor voltage and current in p. u., respectively v qr , i qr q-axis rotor voltage and current in p. u., respectively V, V s , I s bus voltage and stator voltage phasor, stator current phasor in p.û V D , V E voltage of detailed wind farm system and equivalent model, respectively. v sg , v s rotor speed and synchronous speed of SG in p.u. v t ; v r angular speed of the wind turbine and SCIG in p.u. X distance between two wind turbine rotors in meters. X WTD distance between two wind turbine rotors in meters. X 0 q , X 0 d quadrature and direct axis transient reactance of SG, respectively X 0 s , X s , X r stator transient reactance, stator and rotor reactance of IG, respectively X m , X l magnetizing reactance and leakage reactance of IG, respectively X D , X E detailed and equivalent wind farm (WF) parameters. ½Y admittance matrix of the wind farm