AGGREGATED MODEL OF LARGE-SCALE WIND FARMS FOR POWER SYSTEM SIMULATION SOFTWARE TOOLS

Abstract

A method of modeling an equivalent wind turbine generator (WTG) system for a wind farm having a plurality of WTG units includes determining an impact factor of each WTG unit of the plurality of WTG units, determining an equivalent single WTG unit model parameters of the wind farm based on the impact factor of each WTG unit, and determining an effective wind speed of the wind farm to use as the equivalent WTG input wind speed. The method produces a model of static and/or dynamic wind farm behavior. Additionally, a software configured to execute a method of modeling an equivalent wind turbine generator (WTG) system for a wind farm having a plurality of WTG units.

Claims

1. A method of modeling an equivalent wind turbine generator (WTG) system for a wind farm having a plurality of WTG units, the method comprising: determining an impact factor of each WTG unit of the plurality of WTG units; determining an equivalent single WTG unit model parameters of the wind farm based on the impact factor of each WTG unit; and determining an effective wind speed of the wind farm to use as the equivalent WTG input wind speed.

2. The method of claim 1 wherein the method produces a model of static wind farm behavior.

3. The method of claim 1 wherein the method produces a model of dynamic wind farm behavior.

4. The method of claim 1 wherein the plurality of WTG units are fixed speed units.

5. The method of claim 1 wherein a frequency response technique is employed to determine the impact factor of each WTG unit.

6. The method of claim 1 wherein the wind farm features a plurality of wind speed inputs.

7. The method of claim 6 wherein the plurality of wind speed inputs are in a plurality of locations throughout the wind farm.

8. The method of claim 1 wherein the plurality of WTG units feature a plurality of machine parameters.

9. The method of claim 8 further comprising determining an effect of plurality of machine parameters on each WTG unit of the plurality of WTG units.

10. The method of claim 1 wherein the impact factor of a WTG unit is the proportion of the WTG unit output current increment to a total wind farm output current increment.

11. The method of claim 1 wherein the mechanical input of the equivalent WTG unit is the sum of a total mechanical input of all WTG units in the wind farm.

12. The method of claim 1 wherein the mechanical input increment of the equivalent WTG unit is a sum of a total mechanical input increment of all WTG units in the wind farm.

13. The method of claim 1 further comprising determining an equivalent collector system model parameter of the wind farm based on the impact factor and an equilibrium point of each WTG unit.

14. The method of claim 1 wherein determining an equivalent single WTG unit model parameters of the wind farm based on the impact factor of each WTG unit includes first determining an equivalent electrical side of the windfarm including at least one generator or converter based on associated impact factors, and then determining an equivalent mechanical side of the windfarm including at least one wind turbine based on the equivalent electrical side and the impact factors.

15. A computer comprising machine readable medium having stored thereon one or more sequences of instructions configured to execute the method of claim 1.

16. A method of modeling an equivalent system for parallel systems having a plurality of units, the method comprising: determining an impact factor of each unit proportional to its operating point of the plurality of units; determining an equivalent single unit model parameters of the overall parallel system based on the impact factor of each unit; and determining an effective input of the overall parallel system to use as the equivalent input.

17. The method of claim 16 wherein the parallel systems include any of parallel converters, PV farms, battery banks and renewables.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The present invention is further detailed with respect to the following drawings that are intended to show certain aspects of the present of invention, but should not be construed as limit on the practice of the invention, wherein:

[0013] FIG. 1 shows a virtual wind farm based on an Asynchronous machine according to forms of the present disclosure;

[0014] FIG. 2 shows a virtual wind farm based on a Synchronous machine according to forms of the present disclosure;

[0015] FIG. 3A shows a schematic drawing of a two-zone wind farm according to embodiments of the present invention;

[0016] FIG. 3B shows a schematic diagram of a fixed-speed wind turbine generator having a collector system with capacitors to provide reactive power;

[0017] FIG. 4A shows a schematic of a two fixed-speed WTG wind farm having a collector system with capacitors to provide reactive power;

[0018] FIG. 4B shows a schematic of a semi aggregated two fixed-speed WTG wind farm having a collector system with capacitors to provide reactive power;

[0019] FIG. 5 is a schematic of a wind farm with 80 wind turbine generator units;

[0020] FIGS. 6A-6D are graphs showing active power dynamic behavior of system using Full aggregation, Zone aggregation, Semi aggregation and I.F aggregation models in four different scenarios covering the combinations of various wind speed inputs and different WTGs parameters in a wind farm system according to embodiments of the present invention;

[0021] FIGS. 7A-7D are graphs showing Full aggregation, Zone aggregation, Semi aggregation, and I.F. aggregation methods overall error and simulation time in four different scenarios covering the combinations of various wind speed inputs and different WTGs parameters in a wind farm system according to embodiments of the present invention;

[0022] FIG. 8 is an equivalent circuit diagram of the I.F. Agg. Model according to embodiments of the present invention;

[0023] FIG. 9 is a schematic drawing of a DFIG WTG according to embodiments of the present invention;

[0024] FIGS. 10A and 10B show a DFIG WTG and a Fixed-Speed WTG, respectively;

[0025] FIG. 11A shows windfarm d-q axis circuits;

[0026] FIG. 11B shows equivalent model d-q axis circuits of those shown in FIG. 11A;

[0027] FIG. 12 is a drawing showing rotating the d-q frame by .sub.0;

[0028] FIG. 13 shows a schematic of a 4-WTGs windfarm;

[0029] FIG. 14 shows a schematic of a 20-WTGs windfarm;

[0030] FIG. 15 is a graph showing PCC voltage, phase A current and active and reactive power for FIG. 13 DFIG windfarm with similar WTG parameters and its Full Agg, Zone Agg and WD Agg models;

[0031] FIG. 16 is a graph showing PCC voltage, phase A current and active and reactive power for FIG. 13 DFIG windfarm with similar WTG parameters, its Full Agg and Zone Agg models with modified controllers and its WD Agg model;

[0032] FIG. 17 is a graph showing PCC voltage, phase A current and active and reactive power for FIG. 13 DFIG windfarm with unequal WTG parameters, its WD Agg model and its modified Full Agg and Zone Agg models;

[0033] FIG. 18 is a graph showing wind speed, phase A current and active and reactive power for FIG. 14 DFIG windfarm with unequal WTG parameters, its Full Agg and Zone Agg models with modified controllers and its WD Agg model; and

[0034] FIG. 19 is a graph showing PCC voltage, phase A current and active and reactive power for FIG. 13 fixed-speed windfarm with unequal WTG parameters, its modified Full Agg, Zone Agg and WD Agg models.

DETAILED DESCRIPTION OF THE INVENTION

[0035] The present invention has utility as a method for aggregately modeling a wind farm capable of quantifying the contribution of each WTG in a large-scale wind farm and as a power system simulation software for large-scale wind farms for considering the impact factor and contribution of each of hundreds of wind turbine-generator units making up the large-scale wind farm as seen by a power grid from the point of common coupling. The method of the present invention is highly accurate while being efficient and low-cost for large-scale wind farms.

[0036] The present invention provides an Impact Factor aggregation (I.F. Agg.) method that includes the contribution of each WTG unit, based on its parameters and operating point, within the equivalent model of wind farm. The method provides a computationally efficient model for a wind farm that significantly improves the accuracy compared with Full Agg. method. The reason is that I.F. Agg. method includes the effects of WTGs with different parameters and/or operating points. The inventive I.F. Agg. method analytically calculates the contribution of each WTG unit as a weight function in frequency domain. This technique allows one to determine the best set of equivalent model parameters to improve the model accuracy over the frequency range of interest. Most of existing methods develop and test the performance of their aggregated models mainly for fixed-speed wind farms to explain the main concept of the methods for the simplest wind farm configuration. Furthermore, the inventive method includes the effect of wind farm collector system in the equivalent model that is less considered and discussed in the other existing methods.

[0037] The present invention additionally provides a Weighted Dynamic (WD) model for a windfarm which has several advantages over the existing methods. In this method, the contribution of each WTG in the aggregated model is quantified which significantly increases the model accuracy compared to the Full Agg model. The derived set of aggregated dynamic equations results in a simpler modeling approach compared to equivalent admittance methods. Also, the windfarm equivalent WTG is obtained through a one-time calculation leading to a much lower computation burden compared to the optimization algorithms. Furthermore, the whole windfarm is modeled with a single WTG resulting in a simpler model compared to the Zone Agg and Semi Agg models. It is worth mentioning that the proposed equivalent turbine for the whole windfarm has not been clearly addressed in the existing literature which provides a simpler model and better insight regarding the system behavior.

[0038] The performance of the inventive method is evaluated based on the comparisons of a 4-WTGs DFIG windfarm and the obtained equivalent model. Moreover, a time-domain simulation of 20-WTGs DFIG windfarm with a variable wind speed curve is studied to verify the applicability of the proposed model in a more realistic scenario. Finally, a 4-WTGs fixed speed windfarm are studied to demonstrate the generality of the inventive method. Simulations results conducted under identical and unequal WTGs operating conditions demonstrate a superior or at least similar performance of the proposed method compared to the existing approaches. The inventive model results are compared with the Full Agg and Zone Agg models because the computational burden of these methods are the superior of the existing ones and comparable to the inventive method.

[0039] It is to be understood that in instances where a range of values are provided that the range is intended to encompass not only the end point values of the range but also intermediate values of the range as explicitly being included within the range and varying by the last significant figure of the range. By way of example, a recited range of from 1 to 4 is intended to include 1-2, 1-3, 2-4, 3-4, and 1-4.

[0040] According to embodiments, the software element can be used for both steady-state and dynamic analyses of power systems including large-scale wind farms. The system significantly reduces the computational burden of a computer and its memory usage for power system simulation by replacing a large-scale wind farm including hundreds of wind turbine-generator units with an equivalent functional model. Compared to the existing aggregated models, the inventive functional model according to embodiments of the present invention supports a wind farm with different wind speed zones. It also supports modeling of a wind farm with different ratings of wind turbine generator units

[0041] Embodiments of present invention utilize Impact Factors (I.F.) of a WTG in a wind farm to quantify the contribution of each WTG in an aggregation model. According to embodiments, the I.F. aggregation method uses the frequency response technique to find the best match between the aggregated model parameters and wind farm based on d-q reference frame model of the wind farm generators. Using the I.F. concept results in the model having least amount of error and simulation time overall compared to Full aggregation, Zone aggregation and Semi aggregation methods for both steady-state and transient analyses. According to embodiments, the performance of the method is evaluated based on time-domain simulation of fixed-speed wind farm including 80 WTGs. The time-domain investigation compares the simulation results for the aggregation of the wind farm by Full aggregation, Zone aggregation, Semi aggregation and I.F. aggregation methods under four different scenarios. These test scenarios cover the combinations of various wind speed inputs and different WTGs parameters in the wind farm test system.

[0042] Furthermore, the inventive method includes the effect of wind farm collector system in the equivalent model that is less considered and discussed in the other existing methods. According to embodiments of the I.F. Agg. method, a wind farm including 80 WTGs is fully modeled using MATLAB/SIMULINK software tool. The time domain dynamic and steady state behavior of this wind farm is obtained and used as a reference to evaluate and compare different aggregation methods. Various test scenarios are defined including a combination WTGs with similar/different parameters and a wind farm with uniform/nonuniform wind speed distributions. The present invention also defines a normalized index to quantify computational burden, and accuracy to present superior features of the inventive I.F. Agg. method compared with the other methods.

[0043] FIG. 3A shows the diagram of a wind farm including tens of fixed-speed WTG units that are connected to a grid at the point of common coupling (PCC) through a collector system. In a large-scale wind farm, the wind speed distribution at WTG units is not necessarily uniform. To improve the accuracy, the wind farm can be partitioned in zones as shown in FIG. 3A. In partitioning of a large-scale wind farm, in addition to wind speed, several operating point parameters such as slip ratio, and active power must be used in a clustering algorithm. A fixed speed wind farm includes several WTG units with a given schematic as shown in FIG. 3B. The aerodynamic model of the turbine is shown in Equation 1.

P.sub.m=C.sub.p(,)P.sub.W=T.sub.m.sub.r,Equation 1:

where P.sub.W=0.5r.sup.2V.sub.W.sup.3 the wind power. The parameters , r, and V.sub.W denote the air density, turbine radius, and the wind speed, respectively. C.sub.p(,) is the turbine coefficient that is a function of the pitch angle of the turbine blades, , and the tip speed ratio, r.sub.l/V.sub.W, where w.sub.l is the turbine shaft speed. P.sub.m, T.sub.m, and w.sub.r are the generator mechanical power, torque, and speed, respectively. For a given pitch angle, C.sub.p can be estimated with a quadratic function as Equation 2.

[00001]$\begin{array}{cc}{C}_p\left(\right)=\frac{C_{\mathrm{pm}}}{{}_{\mathrm{opt}}^2}.Math.\left(2.Math.{}_{\mathrm{opt}}-\right).Math.,& \mathrm{Equation}.Math..Math.2\end{array}$

where C.sub.pm is the maximum of Cp that occurs at =.sub.opt. The turbine coefficient in Equation 2 can be referred to generator shaft in Equation 3.

[00002]$\begin{array}{cc}{C}_p\left({}^{}\right)=\frac{C_{\mathrm{pm}}}{{}_{\mathrm{opt}}^{{}^2}}.Math.\left(2.Math.{}_{\mathrm{opt}}^{}-{}^{}\right).Math.{}^{},& \mathrm{Equation}.Math..Math.3\end{array}$

where the referred parameters are =rw.sub.r/V.sub.W=G and .sub.opt=G.sub.opt, and G_=w.sub.r/w.sub.l is gear box turns ratio in a WTG. Using this notation, the steady-state model of WTG at the generator side is shown in Equation 4.

T.sub.m+T.sub.e=D.sub.r/.sub.b,Equation 4:

where T.sub.e=(X.sub.m.sup.2R.sub.rs.sub.0V.sub.e)/T.sub.e is the generator electric torque, w.sub.b is the based frequency, D is mechanical damping coefficient, and T.sub.e=[R.sub.sR.sub.r+s.sub.0(X.sub.m.sup.2X.sub.ssX.sub.rr)].sup.2+[R.sub.rX.sub.ss+s.sub.0R.sub.sX.sub.rr].sup.2.

[0044] Vs is the effective voltage at PCC, R.sub.s and R.sub.r are the stator and rotor resistances, and X.sub.m is the magnetizing reactance, respectively. The slip of induction generator at the operating point is s.sub.0=(w.sub.bw.sub.r)/w.sub.b and the machine reactances are X.sub.ss=X.sub.m+X.sub.ls and X.sub.rr=X.sub.m+X.sub.lr, where X.sub.ls and X.sub.lr are leakage reactances of stator and rotor, respectively.

[0045] The schematics of equivalent system corresponding to Full Agg. method is depicted on FIG. 3B, in which the wind farm is replaced with a single unit WTG assuming that the WTG parameters and wind speed of all zones (zone 1 and 2 in FIG. 3A) are the same for all n WTG units within the farm. Based on these assumptions, the per unitized parameters of the generator in Full Agg. model is the same as that of calculated for a single WTG unit in section A, provided that S.sub.beq=nS.sub.b is used as the base power where S.sub.b is the rated power of WTGs and S.sub.beq is the rated power of the wind farm.

[0046] Two limitations of the Full Agg. model are ambiguities in the definition of mechanical parameters and modeling of real wind farms including machines with different ratings and various wind speeds within the zones. It has previously been proposed that the total mechanical power of the wind farm is calculated and applied to the equivalent generator without considering a model for wind turbine. The variable wind speed at different zones causes steady-state and dynamic errors when aggregated model for a large-scale wind farm is used. To mitigate the error of wind speed mismatch at different zone, the concept of equivalent effective wind speed of wind farm with unison WTGs is defined as Equation 5.

[00003]$\begin{array}{cc}.Math.V_W=\frac{1}{n}.Math.{\left(\underset{k=1}{\overset{n}{.Math.}}.Math.V_{W_k}^3\right)}^{\text{?}},.Math.\text{?}.Math.\text{indicates text missing or illegible when filed}& \mathrm{Equation}.Math..Math.5\end{array}$

where V.sub.W is the wind speed that provides a power equals to the total power of wind farm.

[0047] Zone Agg. method partitions the wind farm into a few zones with respect to wind speed variations, and other operating point parameters of WTG units. Then, a Full Agg. model is associated to each zone to represent WTGs within the zone with an equivalent system as shown in FIG. 2(a). This method significantly improves the accuracy of equivalent models since each Full Agg. model is assigned to a zone with an almost fixed wind speed. However, improving the accuracy in this method is obtained at the cost of model complexity since the number of zones may be high in large-scale wind farms. Furthermore, the other limitations of Full Agg. method, including the inability to model WTGs with different parameters, will be directly induced to the Zone Agg. method.

[0048] To further improve the aggregated model accuracy, Semi Agg. method is proposed, in which the wind farm generators are represented with a single per unitized generator similar to the one in Full Agg. method. However, the wind turbines are individually modeled to calculate mechanical torque separately, as depicted on FIG. 4B. This method keeps the mechanical part of wind farm intact and only the generators are aggregated. Therefore, the method will become computationally inefficient for large-scale wind farms although it provides more accurate results compared with Full and Zone Agg. methods.

[0049] The inventive I.F. Agg. method quantifies the impact of each WTG within a wind farm to develop a more accurate equivalent system for the wind farm. This method introduces an equivalent WTG unit for the wind farm and determines its parameters based on weighted average of the WTGs within the farm. The weighting function is defined as the incremental ratio of WTG current to the wind farm current at PCC.

[0050] The weighted averaging technique can be analytically realized in frequency domain that needs the full model of WTG to be linearized about its operating point. Then, the technique will be used for WTG models in frequency domain to obtain the equivalent model parameters. An advantage of using I.F. Agg. method is that it can also define an equivalent RC model for the collector system of the wind farm. It will be shown that the weighted averaging technique significantly improves the accuracy of aggregation model while it remains computationally efficient and addresses the limitation of existing methods.

[0051] The first step of I.F. method needs to determine the operating point of WTG units. Based on Equation 1, the input mechanical power corresponding to the k-th WTG unit is

[00004]$P_{m_{k_0}}=C_{p_k}\left({}_{k_0}^{},{}_{k_0}\right).Math.P_{w_{k_0}}$

where the subscript 0 signifies quantities at the operating point. For the sake of brevity, subscript k is removed within the rest this section till it is needed for merging the equations. The rated slip of high power induction generators is small (e.g. 0.005 for MW-scale generators), therefore, the mechanical speed of generator shaft can be approximated as r.sub.r0=w.sub.b/p where p is the number of pole pair of generator in a WTG. As the dynamic of pitch control system is slow compared with the power system dynamics, the pitch angle .sub.0 can be assumed constant corresponding to a fixed wind speed. Thus, P.sub.m0 is given as P.sub.m0=Cp(.sub.0)P.sub.W0 where .sub.0rw.sub.b/(p.Math.V.sub.W). Hence, to obtain the slop at WTG operating point, the per unitized P.sub.m0 and w.sub.r0=1s.sub.0 can be substituted in Equation 4 and solved it for so.

[0052] The next step in I.F. Agg. method is to obtain the impact factors of WTGs that are used for weighted averaging of machine parameters to obtain an equivalent WTG for the wind farm. This averaging can be appropriately performed in frequency domain to cover the frequency range of interest for power system studies. The linearized mechanical model of an squirrel cage induction generator is shown in Equation 6.

[00005]$\begin{array}{cc}.Math..Math.T_m=X_m.Math.i_{\mathrm{dr}.Math..Math.0}.Math..Math..Math.i_{\mathrm{qs}}-X_m.Math.i_{\mathrm{qr}.Math..Math.0}.Math..Math..Math.i_{\mathrm{ds}}-X_m.Math.i_{\mathrm{ds}.Math..Math.0}.Math..Math..Math.i_{\mathrm{qr}}++.Math.X_m.Math.i_{\mathrm{qs}.Math..Math.0}.Math..Math..Math.i_{\mathrm{dr}}-2.Math.H.Math.\frac{d.Math..Math.{}_r}{\mathrm{dt}}.& \mathrm{Equation}.Math..Math.6\end{array}$

[0053] For a WTG and based on Equations 1-3, one obtains Equation 7.

[00006]$\begin{array}{cc}{{}_{r.Math..Math.0}.Math..Math..Math.T_m+T_{m.Math..Math.0}.Math.{}_r=P_{W_0}.Math.\frac{C_p}{{}^{}}.Math.}_{{}^{}={}_0^{}}.Math.{}^{}+C_{p.Math..Math.0}.Math..Math..Math.P_W,& \mathrm{Equation}.Math..Math.7\end{array}$

where =rW.sub.r/V.sub.W0. Solving Equation 7 for T.sub.m yields Equation 8.

[00007]$\begin{array}{cc}.Math..Math.T_m=e_m.Math.{}_r,.Math.e_m=\frac{1}{{}_{r_0}}.Math.{(\frac{C_p}{{}^{}}.Math.}_{{}^{}={}_0^{}}.Math.\frac{{\mathrm{rP}}_{W_0}}{V_{W_0}}-T_{m_0})..Math.\mathrm{where}.& \mathrm{Equation}.Math..Math.8\end{array}$

[0054] The linearized voltage equations of the machine in frequency domain are shown in Equations 9-12.

[00008]$\begin{array}{cc}.Math..Math.v_{\mathrm{qs}}=A_{\mathrm{qs}}.Math..Math..Math.i_{\mathrm{qs}}+\underset{-}{B_{\mathrm{qs}}}.Math..Math..Math.i_{\mathrm{ds}}+C_{\mathrm{qs}}.Math..Math..Math.i_{\mathrm{qr}}+D_{\mathrm{qs}}.Math..Math..Math.i_{\mathrm{dr}},& \mathrm{Equation}.Math..Math.9\\ .Math..Math.v_{\mathrm{ds}}=A_{\mathrm{ds}}.Math..Math..Math.i_{\mathrm{qs}}+B_{\mathrm{ds}}.Math..Math..Math.i_{\mathrm{ds}}+C_{\mathrm{ds}}.Math..Math..Math.i_{\mathrm{qr}}+D_{\mathrm{ds}}.Math..Math..Math.i_{\mathrm{dr}},& \mathrm{Equation}.Math..Math.10\\ .Math..Math.v_{\mathrm{qr}}=0=A_{\mathrm{qr}}.Math..Math..Math.i_{\mathrm{qs}}+B_{\mathrm{qr}}.Math..Math..Math.i_{\mathrm{ds}}+C_{\mathrm{qr}}.Math..Math..Math.i_{\mathrm{qr}}+D_{\mathrm{qr}}.Math..Math..Math.i_{\mathrm{dr}},& \mathrm{Equation}.Math..Math.11\\ .Math..Math.v_{\mathrm{dr}}=0=A_{\mathrm{dr}}.Math..Math..Math.i_{\mathrm{qs}}+B_{\mathrm{dr}}.Math..Math..Math.i_{\mathrm{ds}}+C_{\mathrm{dr}}.Math..Math..Math.i_{\mathrm{qr}}+D_{\mathrm{dr}}.Math..Math..Math.i_{\mathrm{dr}},& \mathrm{Equation}.Math..Math.12\\ A_{\mathrm{qs}}=B_{\mathrm{ds}}=R_s+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{ss}},B_{\mathrm{qs}}=-A_{\mathrm{ds}}=X_{\mathrm{ss}},.Math.C_{\mathrm{qs}}=D_{\mathrm{ds}}=\frac{j.Math..Math.}{{}_b}.Math.X_m,D_{\mathrm{qs}}=-C_{\mathrm{ds}}=X_m,& \\ \mathrm{where}.Math..Math.\mathrm{and}& \\ A_{\mathrm{qr}}=\frac{j.Math..Math.}{{}_b}.Math.X_m-X_d.Math.i_{\mathrm{dr}.Math..Math.0},A_{\mathrm{dr}}=-s_0.Math.X_m+X_q.Math.i_{\mathrm{dr}.Math..Math.0},.Math.B_{\mathrm{qr}}=s_0.Math.X_m+X_d.Math.i_{\mathrm{qr}.Math..Math.0},B_{\mathrm{dr}}=\frac{j.Math..Math.}{{}_b}.Math.X_m-X_q.Math.i_{\mathrm{qr}.Math..Math.0},.Math.C_{\mathrm{qr}}=R_r+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{rr}}+X_d.Math.i_{\mathrm{ds}.Math..Math.0},C_{\mathrm{dr}}=-s_0.Math.X_{\mathrm{rr}}-X_q.Math.i_{\mathrm{ds}.Math..Math.0},.Math.D_{\mathrm{qr}}=s_0.Math.X_{\mathrm{rr}}-X_d.Math.i_{\mathrm{qs}.Math..Math.0},D_{\mathrm{dr}}=R_r+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{rr}}+X_q.Math.i_{\mathrm{qs}.Math..Math.0}.& \end{array}$

[0055] Solving Equations 11 and 12 for i.sub.qr and i.sub.dr and substituting the solutions in Equations 9 and 10 yield Equations 13 and 14.

v.sub.qs=.sub.q(j)i.sub.qs+.sub.q(j)i.sub.ds,Equation 13:

v.sub.ds=.sub.d(j)i.sub.qs+.sub.d(j)i.sub.ds,Equation 14:

[0056] Finally, by solving Equation 14 for i.sub.ds and substituting the solution in Equation 13, v.sub.qsk for the k-th unit can be expressed in terms of two transfer functions K.sub.k(jw) and G.sub.k(jw) as given by Equation 15.

v.sub.qs.sub.k=K.sub.k(j)i.sub.qs.sub.k+G.sub.k(j)v.sub.ds.sub.k,Equation 15:

where

K.sub.k(j)=.sub.q.sub.k(j).sub.q.sub.k(j).sub.d.sub.k(j)/.sub.d.sub.k(j),

G.sub.k(j)=.sub.q.sub.k(j)/.sub.d.sub.k(j).

[0057] The WTG units are connected in parallel through a collector system that is often design for negligible power losses (e.g. less than 2%) at rated power of wind farm. Thus, to develop equivalent system of a wind farm including n WTGs, it can be assumed that v.sub.qskv.sub.qs and v.sub.dskv.sub.ds for k=1, 2, . . . , n where v.sub.qs and v.sub.ds are dq-components of the wind farm at the point of common coupling. Thus, by applying summation over k=1, 2, . . . , n in Equation 15, the wind farm model in frequency domain can be expressed as Equation 16.

[00009]$\begin{array}{cc}n.Math..Math..Math..Math.v_{\mathrm{qs}}=.Math..Math.i_{\mathrm{qs}}.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.\frac{K_k\left(j.Math..Math.\right).Math..Math..Math.i_{{\mathrm{qs}}_k}}{.Math..Math.i_{\mathrm{qs}}}+.Math..Math.v_{\mathrm{ds}}.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.G_k\left(j.Math..Math.\right),& \mathrm{Equation}.Math..Math.16\end{array}$

[0058] The impact factor is defined as u.sub.ki.sub.qsk/i.sub.qs at w=0, i.e. the dc gain of incremental current ratios since w=0 in dq frame corresponding to the fundamental frequency of the generator in time domain. Then, Equation 16 can be expressed as Equation 17.

[00010]$\begin{array}{cc}.Math..Math.v_{\mathrm{qs}}=\frac{K^{}\left(j.Math..Math.\right)}{n}.Math..Math..Math.i_{\mathrm{qs}}+\frac{G^{}\left(j.Math..Math.\right)}{n}.Math..Math..Math.v_{\mathrm{ds}},.Math.\mathrm{where}.Math..Math.K^{}\left(j.Math..Math.\right).Math.\stackrel{}{=}.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.K_k.Math.u_k,G^{}\left(j.Math..Math.\right).Math.\stackrel{}{=}.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.G_k.& \mathrm{Equation}.Math..Math.17\end{array}$

[0059] Finally, by updating the base apparent power to the rating of wind farm, i.e. SWF=nSWTG, Equation 17 can be rearranged in wind farm per unit system as Equation 18.

v.sub.qs=K(j)i.sub.qs+G(j)v.sub.ds,Equation 18:

[0060] Based on Equation 18, the frequency domain model of wind farm is formulated similar to a single unit WTG as given in Equation 15. The schematic diagram of the equivalent circuit for this model is shown in FIG. 8. As the wind farm model is developed based on linearized model of WTGs, the parameters of equivalent model of wind farm in FIG. 8 can be readily determined using the impact factors u.sub.k, 1, 2, . . . , n, e.g. X.sub.rr.sup.[pu]=.sub.k=1.sup.nX.sub.rr.sub.k.sup.[pu]u.sub.k, X.sub.ss.sup.[pu]=.sub.k=1.sup.nX.sub.ss.sub.k.sup.[pu]u.sub.k, etc. The main difference between the inventive I.F. Agg. and Full Agg. methods and in fact the key contribution of the I.F. Agg. is to calculate the parameters of equivalent model based on weighted averaging using impact factors u.sub.k to provide more accurate results compared with existing methods.

[0061] An equivalent collector system and shunt capacitors can be defined with equivalent R.sub.C and C.sub.C as shown in FIG. 3B. R.sub.C can be determined from the steady state power balance or it can be obtained based on weighted averaging method as R.sub.C.sup.[pu]=.sub.k=1.sup.nR.sub.C.sub.k.sup.[pu]u.sub.k.

[0062] Using FIG. 3B and considering the steady state equivalent circuit of induction generator, the per phase reactive demand from grid for k-th unit is shown in Equation 19.

[00011]$\begin{array}{cc}{Q}_{d_k}=V_{s_k}^2\left(\frac{1}{X_{m_k}}+\frac{X_{\mathrm{ls}}+X_{\mathrm{lr}}}{{\left(\frac{R_r}{s}\right)}^2+{\left(X_{\mathrm{ls}}+X_{\mathrm{lr}}\right)}^2}-{}_s.Math.C_{C_k}\right).& \mathrm{Equation}.Math..Math.19\end{array}$

where V.sub.sk is the effective voltage at the terminal of the k-th WTG unit. The equivalent capacitor C.sub.C is obtained based on reactive power balance, as given by Equation 20.

[00012]$\begin{array}{cc}{Q}_d=V_s^2\left(\frac{1}{X_m}+\frac{X_{\mathrm{ls}}+X_{\mathrm{lr}}}{{\left(\frac{R_r}{s}\right)}^2+{\left(X_{\mathrm{ls}}+X_{\mathrm{lr}}\right)}^2}-{}_s.Math.C_C\right)==\underset{k=1}{\overset{n}{.Math.}}.Math.V_{s_k}^2\left(\frac{1}{X_{m_k}}+\frac{X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}}{{\left(\frac{R_{r_k}}{s_k}\right)}^2+{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)}^2}-{}_s.Math.C_{C_k}\right).& \mathrm{Equation}.Math..Math.20\end{array}$

[0063] Considering a low loss collector system, V.sub.skV.sub.s for k=1, 2, . . . , n. Thus, Equation 20 yields Equation 21.

[00013]$\begin{array}{cc}{C}_C=\frac{1}{{}_s.Math.X_m}+\frac{\left(X_{\mathrm{ls}}+X_{\mathrm{lr}}\right)/{}_s}{{\left(\frac{R_r}{s}\right)}^2+{\left(X_{\mathrm{ls}}+X_{\mathrm{lr}}\right)}^2}--.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.\left(\frac{1}{{}_s.Math.X_{m_k}}+\frac{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)/{}_s}{{\left(\frac{R_{r_k}}{s_k}\right)}^2+{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)}^2}-C_{C_k}\right).& \mathrm{Equation}.Math..Math.21\end{array}$

[0064] To determine the parameters of equivalent wind turbine for a wind farm, assuming P.sub.W=.sub.k=1.sup.nP.sub.W.sub.k, Equation 22 is obtained.

[00014]$\begin{array}{cc}{P}_W=\frac{1}{2}.Math..Math..Math.{\mathrm{AV}}_W^3=\underset{k=1}{\overset{n}{.Math.}}.Math.P_{W_k}=\underset{k=1}{\overset{n}{.Math.}}.Math.\frac{1}{2}.Math..Math..Math.A_k.Math.V_{W_k}^3,& \mathrm{Equation}.Math..Math.22\end{array}$

where A=.sub.k=1.sup.nA.sub.k is the equivalent surface of all WTGs. Thus, the equivalent radius, r, and wind speed, V.sub.W, can be expressed as Equation 23.

[00015]$\begin{array}{cc}.Math.r={\left(\underset{k=1}{\overset{n}{.Math.}}.Math.r_k^2\right)}^{\frac{1}{2}},V_W={\left(\underset{k=1}{\overset{n}{.Math.}}.Math.A_k.Math.V_{W_k}^3/\underset{k=1}{\overset{n}{.Math.}}.Math.A_k\right)}^{\text{?}}..Math.\text{?}.Math.\text{indicates text missing or illegible when filed}& \mathrm{Equation}.Math..Math.23\end{array}$

[0065] The equivalent mechanical power, P.sub.m=.sub.k=1.sup.nP.sub.m.sub.k, and incremental power P.sub.m={dot over ()}.sub.k=1.sup.nP.sub.m.sub.k, equations yield Equation 24 and Equation 25.

[00016]$\begin{array}{cc}{C}_p\left({}^{}\right).Math.P_W=\underset{k=1}{\overset{n}{.Math.}}.Math.C_{p_k}\left({}_k^{}\right).Math.P_{W_k}& \mathrm{Equation}.Math..Math.24\\ P_W.Math.\frac{C_p}{{}^{}}.Math.\frac{r}{V_W}.Math..Math..Math.{}_r=\underset{k=1}{\overset{n}{.Math.}}.Math.P_{W_k}.Math.\frac{C_{p_k}}{{}_k^{}}.Math.\frac{r_k}{V_{W_k}}.Math..Math..Math.{}_{r_k}.& \mathrm{Equation}.Math..Math.25\end{array}$

[0066] Therefore, equivalent C.sub.pm and .sub.opt in Equation 3 can be obtained from the simultaneous solutions of Equation 24 and Equation 25. Finally, .sub.opt and gear-box ratio, G, for the equivalent wind turbine generator can be defined based on weighted average of .sub.optk with respect to radii r.sub.k for k=1, 2, . . . , n as Equation 26.

[00017]$\begin{array}{cc}\frac{{}_{\mathrm{opt}}}{r}.Math.\stackrel{}{=}.Math.\frac{1}{n}.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.\frac{{{}_{\mathrm{opt}}}_k}{r_k},G.Math.\stackrel{}{=}.Math.\frac{{}_{\mathrm{opt}}^{}}{{}_{\mathrm{opt}}}.& \mathrm{Equation}.Math..Math.26\end{array}$

[0067] Next, the performance and accuracy of the inventive I.F. Agg. method in comparison with Full, Zone, and Semi Agg. methods are compared using a fixed speed wind farm study system. According to embodiments, the system includes 80 WTG units as shown in FIG. 5 that are interconnected using a collector system with a given schematic in FIG. 3B. The wind farm is fully simulated using detailed parameters given in Table 1 and the simulation test results are used to generate the reference curves and points for calculation of the accuracy.

[0068] Furthermore, the simulation time of the 80-WTG wind farm is considered as the reference to compare the computational efficiency of different methods. Two types of generators with different parameters and ratings (Types I and II) are used in four different test scenarios A, B, C, and D, in which WTGs can have various wind speeds. The parameters of wind turbine generators Type I and II are listed in Table I and the details of test scenarios are as follows: [0069] A. All WTG units are Type I with the same wind speeds, V.sub.Wk=20 m/s for 1k80; [0070] B. All WTG units are Type I with different wind speeds at

[00018]$V_{W_k}=\{\begin{array}{cc}20& 1k20\\ 16& 21k40\\ 12& 41k60\\ 8& 61k80\end{array};$ [0071] C. The first half of WTGs are Type I, the second half are Type II, and they all operate at same wind speeds, V.sub.Wk=20 m/s for 1k80; [0072] D. The first half of WTGs are Type I, the second half are Type II, and the wind speed of WTGs are:

[00019]$V_{W_k}=\{\begin{array}{cc}20& 1k10.Math..Math.\mathrm{or}.Math..Math.41k50\\ 16& 11k20.Math..Math.\mathrm{or}.Math..Math.51k60\\ 12& 21k30.Math..Math.\mathrm{or}.Math..Math.61k12\\ 8& 31k40.Math..Math.\mathrm{or}.Math..Math.71k80\end{array}.$

TABLE-US-00001 TABLE I PARAMETERS OF TYPE I AND II WTG UNITS Name of the parameter Type I Type II Units GENERATOR S.sub.b 150 110 [kVA] V.sub.s 460 460 [V] f.sub.s 60 60 [Hz] H 0.3096 0.3175 [s] D 0.0114 0.006839 [pu] R.sub.s 0.01282 0.01597 [pu] R.sub.r 0.00702 0.009103 [pu] X.sub.m 2.503 2.183 [pu] X.sub.ss 2.55351 2.23942 [pu] X.sub.rr 2.55351 2.23942 [pu] P (poles) 4 4 COLLECTOR & GRID R.sub.c 0.05 0.05 [pu] C.sub.c 2.63 2.63 [mF] V.sub.th 460 460 [V] R.sub.th 0.01 0.01 [pu] TURBINE V.sub.W.sub.b 12 11.33 [m/s] r 10.6 10 [m] G 23.55 23.55 C.sub.pm 0.4 0.4 .sub.opt 8 8

[0073] These four scenarios cover all events that can occur for a wind farm including different machine types and various wind speeds. The tests start at t.sub.0=3 s by applying a small signal disturbance, that is a limited 3-phase connection to ground via resistances R.sub.f=0.09 pu at PCC for 3 cycles. After removing this small signal disturbance, the wind farm operating point will be back to its prior operating point at t=3.

[0074] To investigate the performance and accuracy of the aggregation methods the Total Normalized Simulation Time and Error (TNSTE) criterion is defined and used to evaluate the proposed and existing methods. TNSTE consists of three main components as: [0075] 1) Normalized simulation time, STE, that is defined as the ratio of detailed wind farm study system simulation time to the one for an aggregation method; [0076] 2) Steady state error of the aggregated methods at their operating points, given by Equation 27.

[00020]$\begin{array}{cc}{e}_{\mathrm{ss}}=.Math.\frac{{\hat{f}}_{\mathrm{ss}}-f_{\mathrm{ss}}}{f_{\mathrm{ss}}}.Math..& \mathrm{Equation}.Math..Math.27\end{array}$ fss in (27) denotes the steady state of f{V, I, P, Q} where V, I, P, and Q are the voltage, current, real, and reactive power at the point of common coupling of wind farm and grid. {circumflex over ()}f represents the corresponding quantity to f.sub.ss that is obtained from an aggregated model for the study state wind farm; [0077] 3) Transient error of the aggregated methods at a specific operating point following to a small signal disturbance that is defined as Equation 28.

[00021]$\begin{array}{cc}{e}_{\mathrm{Trans}}=.Math.\frac{\frac{1}{T}.Math.{}_{t_0}^{T+t_0}.Math.\left(\hat{f}\left(t\right)-f\left(t\right)\right).Math.\mathrm{dt}}{f_{\mathrm{ss}}\left(t_0^{-}\right)}.Math.,& \mathrm{Equation}.Math..Math.28\end{array}$ where t0 is the fault time and f.sub.ss (t.sub.0.sup.) is the steady state quantity prior to the small signal disturbance.

[0078] TNSTE is defined as the summation of these three normalized components as Equation 29.

TNSTE=STE+e.sub.ss+e.sub.Trans.Equation 29:

[0079] The study system is simulated in MATLAB/SIMULINK software tool and the test results for active power following the small disturbance are depicted in FIGS. 6A-6d corresponding to the four test scenarios A, B, C, and D. For each scenario, the accuracy of the inventive I.F. Agg. method and one of the Full, Zone, and Semi Agg. methods are compared with the detailed 80-unit wind farm test results. The conclusions of active power simulations are: [0080] 1) if the WTG parameters and wind speed of the units in a wind farm are similar (an unrealistic assumption), all equivalent aggregated systems accurately describe the wind farm behavior (FIG. 6A); [0081] 2) The accuracy of Full Agg. method is highly sensitive to wind speed variations between the zones (FIGS. 6B and 6D); [0082] 3) The Full Agg. method is also slightly sensitive to non-unison WTG units with different parameter. The Zone and Semi Agg. methods are more accurate compared with Full Agg. method in all test scenarios; [0083] 4) The I.F. Agg. method provides more accurate results compared with Full and Semi Agg. methods in all test scenarios. However, Zone Agg. method in terms of active power matching is slightly more accurate in test scenarios (b) and (d).

[0084] The last conclusion is expected since in scenarios (b) and (d) the wind speed is different in four zones, thus, the Zone Agg. method that uses four equivalent WTGs corresponding to each zone provides better matching with reference in terms of active power. However, it will be shown in the next analysis (FIG. 6) that the advantages of I.F. Agg. method outweigh the Zone Agg. method when the computational efficiency and accuracy of matching other quantities (e.g. reactive power) are taken into account.

[0085] To further evaluating the performance and accuracy of the methods, simulation results are also studied for reactive power, current, and voltage of the wind farm and the results are compared based on TNSTE index as elaborated in Equations 27 to 29.

[0086] FIGS. 7A-7D show details of TNSTE components using color coded bars for each test scenario. The conclusion of test results are as follows: [0087] 1) in all test scenarios the normalized simulation times (STE) for Full Agg. and I.F. methods (shown with dark gray color) are significantly smaller than Zone and Semi Agg. methods. Thus, I.F. and Full Agg. methods are the most computationally efficient ways for developing equivalent systems for wind farms; [0088] 2) in case of uniform wind speed distribution (that is an unrealistic assumption) the normalized simulation time, TSE, is the dominant component in TNSTE index compared with the transient and steady state errors. Under this condition, I.F. and Full Agg. methods provide the best performance (FIGS. 7A and 7C); [0089] 3) in case of nonuniform wind speed distribution the steady state and transient errors comprise the most part of TNSTE. Except for the Zone Agg. method, the steady state and transient errors are dominant components of TNSTE compared with normalized simulation time. Thus, when the wind speeds at the zones are different quantities, I.F. and Zone. Agg. method provide the best performances for aggregation of the wind farm (FIGS. 7B and 7D); [0090] 4) The nonuniform wind speed distribution significantly increases the reactive power errors, as shown by light green and blue bars in FIGS. 7B and 7D. Thus, using only real power indicator is insufficient to evaluate the accuracy of aggregation methods; [0091] 5) Comparing the scenarios A and C for steady state and transient errors show that change of generator parameters slightly increase the errors, however, these are more sensitive to wind speed variations as shown in FIGS. 7B and 7D.

[0092] The overall test results in FIGS. 7A-7D show that when the wind speed is changing between the zones (e.g. in large-scale wind farms) the I.F. and Zone Agg. methods are the top two effective aggregation methods. However, if the wind speed is almost the same for all WTG units (e.g. in small-scale wind farms (then I.S. and Full Agg. methods can be considered as the superior aggregation methods. Furthermore, clustering algorithms for large-scale wind farms need revising in selection their indicators, i.e. reactive power and computational efficiency must be taken into account.

[0093] Based on a newly defined impact factor of WTG units within a wind farm, this paper presents a systematic analytical method to develop an aggregated system for large-scale fixed-speed wind farms. The aggregated model is established based on linearized dq dynamic model of WTG in frequency domain. It also encompasses an equivalent circuit for the collector system of the wind farm that significantly improves the accuracy of the model specially in terms of reactive power balance. Conventional aggregation methods become highly inaccurate when the wind speed at different zones of a large-scale wind farm are unequal. The advantage of the proposed impact factor method is to improve the accuracy of the aggregated model by considering the contribution of each WTG in the equivalent system based on its operating point current. A study system including 80 WTG units is used for performance evaluation and verification of the method. The test results of the different test scenarios show the superior performance and accuracy of the proposed impact factor aggregation method specially for large-scale wind farms with different wind speed zones.

[0094] Furthermore, Xq and Xd in Equations 9-12 are:

[00022]$\begin{array}{cc}{X}_q=\left(X_m.Math.i_{\mathrm{qs}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{qr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right),& \mathrm{Equation}.Math..Math.30\\ X_d=\left(X_m.Math.i_{\mathrm{ds}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{dr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right),& \mathrm{Equation}.Math..Math.31\\ {{}_{m.Math..Math.0}=2.Math.\mathrm{Hj}.Math..Math.+\frac{C_p}{{}^{}}.Math.}_{{}_0^{}}.Math.\frac{{\mathrm{rP}}_{W_0}}{V_{W.Math..Math.0}\left(1-s_0\right)}-\frac{P_{m.Math..Math.0}}{{\left(1-s_0\right)}^2}..Math.\mathrm{where}& \end{array}$

[0095] .sub.q,d and .sub.q,d in Equations 13 and 14 are:

[00023]$\begin{array}{cc}{}_{q,d}=A_{\mathrm{qs},\mathrm{ds}}+{}_q.Math.C_{\mathrm{qs},\mathrm{ds}}+{}_d.Math.D_{\mathrm{qs},\mathrm{ds}},& \mathrm{Equation}.Math..Math.32\\ {}_{q,d}=A_{\mathrm{qs},\mathrm{ds}}+{}_q.Math.C_{\mathrm{qs},\mathrm{ds}}+{}_d.Math.D_{\mathrm{qs},\mathrm{ds}},& \mathrm{Equation}.Math..Math.33\\ \mathrm{where}& \\ {}_q=\left(\frac{A_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.A_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.C_{\mathrm{dr}}}\right),{}_d=\left(\frac{A_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.A_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.D_{\mathrm{qr}}}\right),& \\ {}_q=\left(\frac{B_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.B_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.C_{\mathrm{dr}}}\right),{}_d=\left(\frac{B_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.B_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.D_{\mathrm{qr}}}\right)& \end{array}$

[0096] The incremental speed is w.sub.r=C(jw)i.sub.qs, and C(jw) is:

[00024]$C\left(j.Math..Math.\right)=\frac{X_m\left(i_{\mathrm{dr}.Math..Math.0}-i_{\mathrm{qr}.Math..Math.0}.Math.C_1+i_{\mathrm{qs}.Math..Math.0}.Math.C_3-i_{\mathrm{ds}.Math..Math.0}.Math.C_2\right)}{e_m+2.Math.\mathrm{Hj}.Math..Math.},.Math.\mathrm{where}$$C_1=-{}_d/{}_d,C_2={}_q+{}_q.Math.C_1,\mathrm{and}.Math..Math.C_3={}_d+{}_d.Math.C_1.$

[0097] Virtual Synchronous Machine (VSM) is an inverter connected to the grid which is controlled by a new control method. This new control method helps the inverter acts similar to a synchronous generator in aspect of delivering active and reactive power to the grid. Therefore, VSM contributes to the frequency stability of the grid by an virtual inertia provided in the control loop. Conventional VSM control approach use a fixed value for virtual inertia. But more advanced control techniques change the virtual inertia value accordingly to achieve desired behavior. But changing just one parameter of synchronous machine will result in moving the operation point from the nominal point. Moreover, it may move the system eigenvalues away from the realistic values. To prevent the following issues extra control loops and protection should be added to the system.

[0098] Impact Factor Aggregation method obtains the equivalent d-q model of the wind farm by considering the contribution of each wind turbine generator (WTG) in the model, and use these equation to control an inverter to act like a wind farm. By changing of the operation point, the virtual wind farm will be modeled by connecting or disconnecting of some WTGs while the remaining connected WTGs working near to their operation point. This method resolve the issues mentioned above automatically and let the system work without extra protection.

[0099] To calculate virtual wind farm parameters, every parameter obtained by electrical impact factors in per unit. For example, the calculated X.sub.rr can be obtained by X.sub.rr.sup.[pu]=.sub.m=1.sup.nu.sub.rmX.sub.rrm.sup.[pu], and the base apparent power is obtained as S.sub.b=.sub.j=1.sup.nS.sub.bj.

[0100] The electrical torque of induction machine can be found by Equation 35.

[00025]$T_e=\frac{X_m^2.Math.R_r.Math.s_0.Math..Math.v_s.Math.}{{}_{T_e}}$${}_{T_s}={\left[R_s.Math.R_r+s_0\left(X_m^2-X_{\mathrm{ss}}.Math.X_{\mathrm{rr}}\right)\right]}^2+{\left[R_r.Math.X_{\mathrm{ss}}+s_0.Math.R_s.Math.X_{\mathrm{rr}}\right]}^2$$\mathrm{where}.$

[0101] Steady-state stator d-q currents are shown by Equations 36 and 37.

[00026]$\begin{array}{cc}{i}_{\mathrm{qr}.Math..Math.0}=-\frac{s_0^2.Math.X_m.Math.X_{\mathrm{rr}}}{s_0^2.Math.X_{\mathrm{rr}}^2+R_r^2}.Math.i_{\mathrm{qs}.Math..Math.0}-\frac{s_0.Math.X_m.Math.R_r}{s_0^2.Math.X_{\mathrm{rr}}^2+R_r^2}.Math.i_{\mathrm{ds}.Math..Math.0},.Math.i_{\mathrm{dr}.Math..Math.0}=\frac{s_0.Math.X_m.Math.R_r}{s_0^2.Math.X_{\mathrm{rr}}^2+R_r^2}.Math.i_{\mathrm{qs}.Math..Math.0}-\frac{s_0^2.Math.X_m.Math.X_{\mathrm{rr}}}{s_0^2.Math.X_{\mathrm{rr}}^2+R_r^2}.Math.i_{\mathrm{ds}.Math..Math.0}& \mathrm{Equation}.Math..Math.36\\ i_{\mathrm{qs}.Math..Math.0}=\frac{v_{\mathrm{qs}.Math..Math.0}-\frac{B_{q.Math..Math.0}}{B_{d.Math..Math.0}}.Math.v_{\mathrm{ds}.Math..Math.0}}{A_{q.Math..Math.0}-A_{d.Math..Math.0}},i_{\mathrm{ds}.Math..Math.0}=\frac{v_{\mathrm{ds}.Math..Math.0}-\frac{A_{d.Math..Math.0}}{A_{q.Math..Math.0}}.Math.v_{\mathrm{qs}.Math..Math.0}}{B_{d.Math..Math.0}-B_{q.Math..Math.0}}& \mathrm{Equation}.Math..Math.37\end{array}$

where A.sub.q0, B.sub.q0, A.sub.d0 and B.sub.d0 are shown by Equation 38.

[00027]$\begin{array}{cc}{A}_{q.Math..Math.0}=R_s+\frac{s_0.Math.X_m^2.Math.R_r}{R_r^2+s_0^2.Math.X_{\mathrm{rr}}^2},B_{q.Math..Math.0}=X_{\mathrm{ss}}+\frac{s_0^2.Math.X_m^2.Math.X_{\mathrm{rr}}}{R_r^2+s_0^2.Math.X_{\mathrm{rr}}^2}& \mathrm{Equation}.Math..Math.38\\ A_{d.Math..Math.0}=-X_{\mathrm{ss}}+\frac{s_0^2.Math.X_m^2.Math.X_{\mathrm{rr}}}{R_r^2+s_0^2.Math.X_{\mathrm{rr}}^2},B_{d.Math..Math.0}=R_s+\frac{s_0.Math.X_m^2.Math.R_r}{R_r^2+s_0^2.Math.X_{\mathrm{rr}}^2}& \end{array}$

[0102] Second-order equation of s.sub.0 is

[00028]${}_2.Math.s.Math.\frac{2}{0}-{}_1.Math.s_0+{}_0=0,$

where .sub.2, .sub.1 and .sub.0 are shown by Equation 39.

.sub.2=DR.sub.r+v.sub.s.sup.2, .sub.1=2DR.sub.r+v.sub.s.sup.2, .sub.0=DR.sub.rP.sub.mR.sub.rEquation 39:

[0103] Mechanical linearized equation of induction machine is shown by Equation 40.

T.sub.m=X.sub.mi.sub.dr0i.sub.qsX.sub.mi.sub.qr0i.sub.dsX.sub.mi.sub.ds0i.sub.qr+X.sub.mi.sub.qs0i.sub.dr2Hp.sub.rEquation 40:

[0104] The resulted four electrical linearized equations by reducing w.sub.r from d-q equations of Equations 41-44.

v.sub.qs=A.sub.qs(j)i.sub.qs+B.sub.qs(j)i.sub.ds+C.sub.qs(j)i.sub.qr+D.sub.qs(j)i.sub.drEquation 41:

v.sub.ds=A.sub.ds(j)i.sub.qs+B.sub.ds(j)i.sub.ds+C.sub.ds(j)i.sub.qr+D.sub.ds(j)i.sub.drEquation 42:

v.sub.qr=A.sub.qr(j)i.sub.qs+B.sub.qr(j)i.sub.ds+C.sub.qr(j)i.sub.qr+D.sub.qr(j)i.sub.drEquation 43:

v.sub.dr=A.sub.dr(j)i.sub.qs+B.sub.dr(j)i.sub.ds+C.sub.dr(j)i.sub.qr+D.sub.dr(j)i.sub.drEquation 44:

[0105] Where A.sub.qs(jw), B.sub.qs(jw), C.sub.qs(jw), . . . are shown by Equation 45.

[00029]$\begin{array}{cc}{A}_{\mathrm{qs}}=R_s+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{ss}},B_{\mathrm{qs}}=X_{\mathrm{ss}},C_{\mathrm{qs}}=\frac{j.Math..Math.}{{}_b}.Math.X_m,D_{\mathrm{qs}}=X_m& \mathrm{Equation}.Math..Math.45\\ A_{\mathrm{ds}}=-X_{\mathrm{ss}},B_{\mathrm{ds}}=R_s+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{ss}},C_{\mathrm{ds}}=-X_m,D_{\mathrm{ds}}=\frac{j.Math..Math.}{{}_b}.Math.X_m& \\ .Math.A_{\mathrm{qr}}=\frac{j.Math..Math.}{{}_b}.Math.X_m-\left(X_m.Math.i_{\mathrm{ds}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{dr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{dr}.Math..Math.0}& \\ .Math.B_{\mathrm{qr}}=s_0.Math.X_m+\left(X_m.Math.i_{\mathrm{ds}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{dr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{qr}.Math..Math.0}& \\ .Math.C_{\mathrm{qr}}=R_r+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{rr}}+\left(X_m.Math.i_{\mathrm{ds}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{dr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{ds}.Math..Math.0}& \\ .Math.D_{\mathrm{qr}}=s_0.Math.X_{\mathrm{rr}}-\left(X_m.Math.i_{\mathrm{ds}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{dr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{qs}.Math..Math.0}& \\ .Math.A_{\mathrm{dr}}=-s_0.Math.X_m+\left(X_m.Math.i_{\mathrm{qs}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{qr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{dr}.Math..Math.0}& \\ .Math.B_{\mathrm{dr}}=\frac{j.Math..Math.}{{}_b}.Math.X_m-\left(X_m.Math.i_{\mathrm{qs}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{qr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{qr}.Math..Math.0}& \\ .Math.C_{\mathrm{dr}}=-s_0.Math.X_{\mathrm{rr}}-\left(X_m.Math.i_{\mathrm{qs}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{qr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{ds}.Math..Math.0}& \\ .Math.D_{\mathrm{dr}}=R_r+\frac{j.Math..Math.}{{}_b}.Math.X_{\mathrm{rr}}+\left(X_m.Math.i_{\mathrm{qs}.Math..Math.0}+X_{\mathrm{rr}}.Math.i_{\mathrm{dr}.Math..Math.0}\right).Math.\left(\frac{X_m}{{}_{m.Math..Math.0}}\right).Math.i_{\mathrm{qs}.Math..Math.0}& \\ .Math.{}_{m.Math..Math.0}=2.Math.\mathrm{Hj}.Math..Math.+\frac{C_p}{}.Math.\frac{{\mathrm{rP}}_W}{\left(G\right).Math.V_W\left(1-s_0\right)}-\frac{P_{m.Math..Math.0}}{{\left(1-s_0\right)}^2}& \\ .Math.\mathrm{where}& \end{array}$

[0106] The resulted two d-q linearized equation by reducing rotor and mechanical linearized equations is shown by Equation 46.

v.sub.qs=.sub.q(j)i.sub.qs+.sub.q(j)i.sub.ds, v.sub.ds=.sub.d(j)i.sub.qs+.sub.d(j)i.sub.ds Equation 46:

where .sub.q(jw), .sub.q(jw), .sub.d(jw) and .sub.d(jw) are shown by Equation 47.

[00030]$\begin{array}{cc}{}_q=A_{\mathrm{qs}}+C_{\mathrm{qs}}\left(\frac{A_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.A_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.C_{\mathrm{dr}}}\right)+D_{\mathrm{qs}}\left(\frac{A_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.A_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.D_{\mathrm{qr}}}\right)& \mathrm{Equation}.Math..Math.47\\ {}_q=B_{\mathrm{qs}}+C_{\mathrm{qs}}\left(\frac{B_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.B_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.C_{\mathrm{dr}}}\right)+D_{\mathrm{qs}}\left(\frac{B_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.B_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.D_{\mathrm{qr}}}\right)& \\ {}_d=A_{\mathrm{ds}}+C_{\mathrm{ds}}\left(\frac{A_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.A_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.C_{\mathrm{dr}}}\right)+D_{\mathrm{ds}}\left(\frac{A_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.A_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.D_{\mathrm{qr}}}\right)& \\ {}_d=B_{\mathrm{ds}}+C_{\mathrm{ds}}\left(\frac{B_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.B_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}.Math.C_{\mathrm{dr}}}\right)+D_{\mathrm{ds}}\left(\frac{B_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.B_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}.Math.D_{\mathrm{qr}}}\right)& \end{array}$

[0107] Therefore the electrical impact factors can be found using Equation 48.

[00031]$\begin{array}{cc}{u}_m=\frac{.Math..Math.i_{\mathrm{qsm}}}{.Math..Math.i_{\mathrm{qs}}}& \\ \mathrm{where}.Math.\text{:}& \\ .Math..Math.v_{\mathrm{qs}}=K\left(j.Math..Math.\right).Math..Math..Math..Math.i_{\mathrm{qs}}+G\left(j.Math..Math.\right).Math..Math..Math.v_{\mathrm{ds}},& \mathrm{Equation}.Math..Math.49\\ K\left(j.Math..Math.\right)={}_q\left(j.Math..Math.\right)-\frac{{}_q\left(j.Math..Math.\right).Math.{}_d\left(j.Math..Math.\right)}{{}_d\left(j.Math..Math.\right)}& \\ G\left(j.Math..Math.\right)=\frac{{}_q\left(j.Math..Math.\right)}{{}_d\left(j.Math..Math.\right)}& \end{array}$

[0108] Relation between w.sub.r and i.sub.qs is w.sub.r=C(jw) i.sub.qs, where C(jw) is shown by Equation 50.

[00032]$\begin{array}{cc}C\left(j.Math..Math.\right)=\frac{X_m.Math.i_{\mathrm{dr}.Math..Math.0}-X_m.Math.i_{\mathrm{qr}.Math..Math.0}.Math.C_1\left(j.Math..Math.\right)}{C_0+2.Math..Math.\mathrm{Hj}.Math..Math.}+\frac{-X_m.Math.i_{\mathrm{ds}.Math..Math.0}.Math.C_2\left(j.Math..Math.\right)+X_m.Math.i_{\mathrm{qs}.Math..Math.0}.Math.C_3\left(j.Math..Math.\right)}{C_0+2.Math..Math.\mathrm{Hj}.Math..Math.}& \mathrm{Equation}.Math..Math.50\\ .Math.\mathrm{where},& \\ .Math.C_0=\left(\frac{C_p\left(,\right)}{}.Math.\frac{{\mathrm{rP}}_W}{\left(G\right).Math.V_W}-T_{m.Math..Math.0}\right).Math.\frac{1}{{}_{r.Math..Math.0}},C_1\left(j.Math..Math.\right)=-\frac{{}_d}{{}_d}& \mathrm{Equation}.Math..Math.51\\ C_2\left(j.Math..Math.\right)=\left(\frac{A_{\mathrm{qr}}\left(j.Math..Math.\right)-\frac{D_{\mathrm{qr}}\left(j.Math..Math.\right)}{D_{\mathrm{dr}}\left(j.Math..Math.\right)}.Math.A_{\mathrm{dr}}\left(j.Math..Math.\right)}{{}_q}\right)+\left(\frac{B_{\mathrm{qr}}\left(j.Math..Math.\right)-\frac{D_{\mathrm{qr}}\left(j.Math..Math.\right)}{D_{\mathrm{dr}}\left(j.Math..Math.\right)}.Math.{}_{\mathrm{dr}}\left(j.Math..Math.\right)}{{}_q}\right).Math.C_1\left(j.Math..Math.\right)& \\ C_3\left(j.Math..Math.\right)=\left(\frac{A_{\mathrm{dr}}\left(j.Math..Math.\right)-\frac{C_{\mathrm{dr}}\left(j.Math..Math.\right)}{C_{\mathrm{qr}}\left(j.Math..Math.\right)}.Math.A_{\mathrm{qr}}\left(j.Math..Math.\right)}{{}_d}\right)+\left(\frac{B_{\mathrm{dr}}\left(j.Math..Math.\right)-\frac{C_{\mathrm{dr}}\left(j.Math..Math.\right)}{C_{\mathrm{qr}}\left(j.Math..Math.\right)}.Math.{}_{\mathrm{qr}}\left(j.Math..Math.\right)}{{}_d}\right).Math.C_1\left(j.Math..Math.\right)& \end{array}$

[0109] Therefore, the mechanical impact factors can be found as Equation 52.

[00033]$\begin{array}{cc}{y}_k=\frac{.Math..Math.{}_{r_k}}{.Math..Math.{}_{r_{\mathrm{eq}}}}=\frac{C_k\left(j.Math..Math.\right).Math.K\left(j.Math..Math.\right)}{K_k\left(j.Math..Math.\right).Math.C\left(j.Math..Math.\right)}& \mathrm{Equation}.Math..Math.52\end{array}$

[0110] Equivalent system mechanical relations are shown using Equations 53 and 54.

[00034]$\begin{array}{cc}\underset{i=1}{\overset{n}{.Math.}}.Math.C_{\mathrm{pi}}.Math.P_{\mathrm{Wi}}=C_{\mathrm{peg}}.Math.P_{\mathrm{Weq}},\mathrm{where}.Math..Math.P_{W_{\mathrm{eq}}}=\frac{1}{2}.Math..Math..Math.A_{\mathrm{eq}}.Math.V_{\mathrm{Weq}}^3& \mathrm{Equation}.Math..Math.53\\ \underset{i=1}{\overset{n}{.Math.}}.Math.P_{\mathrm{Wi}}.Math.\frac{C_{\mathrm{pi}}}{{}_i^{}}.Math.\frac{r_i.Math.y_i}{V_{\mathrm{Wi}}}=P_{W_{\mathrm{eq}}}.Math.\frac{C_{\mathrm{peq}}}{{}_{\mathrm{eq}}^{}}.Math.\frac{r_{\mathrm{eq}}}{V_{\mathrm{Weq}}}& \mathrm{Equation}.Math..Math.54\\ \mathrm{where}.Math.\text{:}& \\ {}^{}=G.Math..Math.,A_{\mathrm{eq}}=\underset{k=1}{\overset{n}{.Math.}}.Math.A_k,r_{\mathrm{eq}}=\sqrt{A_{\mathrm{eq}}/}& \mathrm{Equation}.Math..Math.55\\ C_{p_{\mathrm{eq}}}=C_{{\mathrm{pm}}_{\mathrm{eq}}}.Math.f\left({}^{}\right),G_{\mathrm{eq}}=\frac{{}_{{\mathrm{opt}}_{\mathrm{eq}}}^{}}{{}_{{\mathrm{opt}}_{\mathrm{eq}}}},& \\ {}_{{\mathrm{opt}}_{\mathrm{eq}}}=\frac{r}{n}.Math.\underset{k=1}{\overset{n}{.Math.}}.Math.\frac{{}_{{\mathrm{opt}}_k}}{r_k}& \end{array}$

[0111] .sub.opt.sub.ea and C.sub.pm.sub.eq can be found by help of Equations 53 and 54. Equivalent collector parameters can be found by Equations 56 and 57.

[00035]$\begin{array}{cc}.Math.R_{C_{\mathrm{eq}}}^{\left[\mathrm{pu}\right]}=\underset{k=1}{\overset{n}{.Math.}}.Math.R_{C_k}^{\left[\mathrm{pu}\right]}.Math.u_k& \mathrm{Equation}.Math..Math.56\\ C_{C_{\mathrm{eq}}}=\frac{1}{{}_s.Math.X_{m_{\mathrm{eq}}}}+\frac{\left(X_{{\mathrm{ls}}_{\mathrm{eq}}}+X_{{\mathrm{lr}}_{\mathrm{eq}}}\right)/{}_s}{{\left(\frac{R_{r_{\mathrm{eq}}}}{s_{\mathrm{eq}}}\right)}^2+{\left(X_{{\mathrm{ls}}_{\mathrm{eq}}}+X_{{\mathrm{lr}}_{\mathrm{eq}}}\right)}^2}-\underset{k=1}{\overset{n}{.Math.}}.Math.\left(\frac{1}{{}_s.Math.X_{m_k}}+\frac{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)/{}_s}{{\left(\frac{R_{r_k}}{s_k}\right)}^2+{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)}^2}-C_{C_k}\right)& \mathrm{Equation}.Math..Math.57\end{array}$

[0112] Other required equations are:

[00036]$\begin{array}{cc}{V}_s=V_b.Math.\sqrt{\frac{R_s^2+{}_e^2.Math.L_{\mathrm{ss}}^2}{R_s^2+{}_b^2.Math.L_{\mathrm{ss}}^2}}& \mathrm{Equation}.Math..Math.58\\ {}_e=\frac{{}_m+\sqrt{\max \left(0,{}_m^2+X_{\mathrm{corr}}\right)}}{2}& \mathrm{Equation}.Math..Math.59\\ X_{\mathrm{corr}}=3.Math.P\left(v_{\mathrm{qs}}.Math.i_{\mathrm{qs}}-2.Math..Math.R_s.Math.I_s^2\right)/K_{\mathrm{tv}}& \mathrm{Equation}.Math..Math.60\\ K_{\mathrm{tv}}=\frac{3.Math.\left(\frac{P}{2}\right).Math.L_M^2.Math.R_r.Math.V_b^2}{R_r^2\left(R_s^2\right)+{}_b^2.Math.L_{\mathrm{ss}}^2}& \mathrm{Equation}.Math..Math.61\end{array}$

[0113] A large number of WTGs use induction generators with the stator directly connected to the grid. Due to the wide wind speed range, such induction machines operate at high slip away from their nominal speed. The high slip results in high rotor loss, low efficiency and heated rotor in WTGs with a squirrel cage rotor limiting the operating speed range and output power. Hence, the WTGs with a squirrel cage rotor that can efficiently operate close to the nominal speed are called fixed-speed WTGs. To expand the speed range of such induction machines, the rotor can be connect to the grid through an AC/DC/AC variable frequency converter forming a DFIG shown in FIGS. 9 and 10A. This converter delivers a portion of the rotor energy back to the grid resulting in a rotor loss reduction and a larger wind speed operating range.

[0114] FIG. 10B and FIG. 10A show a diagram of a fixed-speed and a DFIG WTG, respectively. Each WTG includes a wind turbine as the main mechanical part, generating mechanical power P.sub.m and mechanical torque T.sub.m according to Equation 62.

P.sub.m=T.sub.m.sub.m=C.sub.p(,)AV.sub.W.sup.3,Equation 62:

[0115] where w.sub.m is the mechanical speed of generator. The mechanical power is related to wind power by turbine coefficient C.sub.p(, ). This factor depends on the structure of the wind turbine. is the blade angle, =rw.sub.l=V.sub.W, C.sub.p(, ) expressed by Equation 63, V.sub.W is wind speed, is air density, A is area covered by the blades, r is blades radius and w.sub.l is blades rotational speed.

[00037]$\begin{array}{cc}{C}_p\left(,\right)=\frac{C_{\mathrm{pmax}}}{{}_{\mathrm{opt}}^2}.Math.\left(2.Math..Math.{}_{\mathrm{opt}}-\right).Math..& \mathrm{Equation}.Math..Math.63\end{array}$

[0116] While another control system can be applied, the present invention uses DFIG control approach, as shown in FIG. 9. In this control system, the Grid-Side Converter (GSC) controls the DC bus voltage to ensure the DC link voltage stability. A PI compensator uses the DC voltage error to generate the GSC d-axis reference current i*.sub.dg. Also, GSC can provide a desired amount of reactive power Q*.sub.g using the q-axis reference current i*.sub.qg. An inner current control loop generates the gating signals for the GSC by means of i*.sub.dg and i*.sub.qg. The Rotor-Side Converter (RSC) controls the induction machine operating under a large range of wind speeds. RSC measures the rotor speed w.sub.m and generates the torque reference T*.sub.e through T*.sub.e=K.sub.optw.sup.2.sub.m/G.sup.3, where K.sub.opt is defined as Equation 64.

K.sub.opt=r.sup.5C.sub.pmax/.sub.opt.sup.3,Equation 64:

and G is the gear-box ratio. The T*.sub.e signal is also used to form the rotor q-axis reference current i*.sub.qr for the RSC current controllers by Equation 65.

[00038]$\begin{array}{cc}{i}_{\mathrm{qr}}^{*}=-\frac{T_e^{*}}{\frac{3}{2}.Math.p.Math..Math.\frac{L_m}{L_s}.Math..Math.{}_s.Math.},& \mathrm{Equation}.Math..Math.65\end{array}$

where p is the number of machine poles, L.sub.m is the magnetizing inductance, L.sub.s is the stator self inductance and |.sub.s| is stator linkage flux and estimated with

[00039]$.Math.{}_s.Math.=\sqrt{\frac{2}{3}}.Math.\frac{V_s}{{}_s},$

where V.sub.s is the rms value of the stator voltage and w.sub.s is the synchronous speed. Also, the rotor d-axis reference current i*.sub.dr is set to zero to use all of the RSC capacity for active power delivery from the rotor windings as the required reactive power for the induction machine is provided by the GSC. It is worth mentioning that, the generality of the inventive method is not limited by the i*.sub.dr set value.

[0117] To find the dynamic model of the wind farm, a small-signal model of the WTG and windfarm should be derived. Small signal model of a WTG requires the steady-state calculations at the operating condition. The steady-state electro-mechanical relationship between the mechanical side and electrical side of a WTG are expressed as Equation 66.

T.sub.m+T.sub.e=D.sub.m.Equation 66:

[0118] T.sub.e is the electrical torque and D is the mechanical damping of the induction machine. The steady state speed, w.sub.m0, can be found with the assumption that the control system is stable so that T*.sub.e=T.sub.e under the steady-state and by substituting Equation 62 and T*.sub.e=K.sub.optw.sup.2.sub.m into Equation 66 and solving it for w.sub.m. Therefore, all other DFIG steady-state parameters can be expressed using Equation 67.

[00040]$\begin{array}{cc}{T}_{e_0}=K_{\mathrm{opt}}.Math.{}_{m_0}^2,I_{\mathrm{qr}}=-\frac{T_{e_0}}{\frac{3}{2}.Math.p.Math..Math.\frac{L_m}{L_s}.Math..Math.{}_s.Math.},I_{\mathrm{dr}}=0,& \mathrm{Equation}.Math..Math.67\\ T_{m_0}=T_{e_0}+D.Math..Math.{}_{m_0},I_{\mathrm{qs}}=-\frac{L_m}{L_s}.Math.I_{\mathrm{qr}},I_{\mathrm{ds}}=\frac{.Math.{}_s.Math.}{L_s}& \end{array}$