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 determining an equivalent single model for a system of parallel systems having a plurality of units with given dynamic equations to feed power to an electrical grid through a point of common coupling, the method comprising: determining an impact factor of each unit of the plurality of units based on a proportion of a unit dynamic variable increment to the total parallel systems dynamic variable increment, wherein dynamic variables used in the determination of the impact factor of each unit of the plurality of units is selected from a group consisting of current, voltage, power, and mechanical speed; and determining coefficients of the dynamic variable in the equivalent single model based on a weighted average of dynamic variable coefficients for each unit of the plurality of units based on the dynamic variable impact factor of each unit of the plurality of units as the weighting; wherein the determining of coefficients of the dynamic variable reduces a computational burden of a computer and its memory usage for modeling the system; and wherein the method determines the equivalent single model for the system of parallel systems having the plurality of units to feed power to the electrical grid through the point of common coupling based on the weighted average of dynamic variable coefficients for each unit of the plurality of units.

2. The method of claim 1 wherein the method produces a static model of the parallel systems.

3. The method of claim 1 wherein the method produces a dynamic model of the parallel systems.

4. The method of claim 3 wherein the parallel systems are a wind farm and features a plurality of wind speed inputs.

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

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

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

8. The method of claim 1 wherein the plurality of 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 unit of the plurality of units.

10. The method of claim 1 wherein a mechanical input of the equivalent unit is the sum of a total mechanical input of all units.

11. The method of claim 1 wherein a mechanical input increment of the equivalent unit is a sum of a total mechanical input increment of all units of the plurality of units in the parallel systems.

12. The method of claim 1 further comprising determining an equivalent collector system model parameter of the parallel systems based on the impact factor and an equilibrium point of each unit of the plurality of units.

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

14. A non-transitory machine readable medium having stored thereon one or more sequences of instructions configured to execute the method of claim 1 based on mechanical inputs, electrical inputs, or a combination thereof.

15. The method of claim 1 wherein the parallel systems are one of wind turbines, parallel converters, photovoltaic farms, or battery banks.

16. A method of modeling a single 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 using a summation of weighted average of the impact factors for each of the plurality of units; and determining an effective input of the overall parallel system to use as the equivalent input by optimizing the modeling of the equivalent system for the parallel systems to feed power to an electrical grid through a point of common coupling; wherein the optimizing of the modeling of the equivalent system reduces the computational burden of a computer and its memory usage for modeling the equivalent system; and wherein the method determines the single equivalent system of the parallel systems having the plurality of units to feed the power to the electrical grid through the point of common coupling based on the weighted average of the dynamics variable coefficients for each unit of the plurality of units.

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

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) 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:

(2) FIG. 1 shows a virtual wind farm based on an Asynchronous machine according to forms of the present disclosure;

(3) FIG. 2 shows a virtual wind farm based on a Synchronous machine according to forms of the present disclosure;

(4) FIG. 3A shows a schematic drawing of a two-zone wind farm according to embodiments of the present invention;

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

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

(7) 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;

(8) FIG. 5 is a schematic of a wind farm with 80 wind turbine generator units;

(9) 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;

(10) 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;

(11) FIG. 8 is an equivalent circuit diagram of the I.F. Agg. Model according to embodiments of the present invention;

(12) FIG. 9 is Doubly-Fed Induction Generator (DFIG) WTG schematic;

(13) FIG. 10A is DFIG WTG, and FIG. 10B is Fixed-Speed WTG;

(14) FIG. llA are windfarm d-q axis circuits and FIG. 10B are equivalent model d-q axis circuits

(15) FIG. 12 depicts rotating the d-q frame by 0;

(16) FIG. 13 depicts a 4-WTGs windfarm;

(17) FIG. 14 depicts a 20-WTGs windfarm.

(18) FIG. 15 depicts PCC voltage, phase A current and active and reactive power for FIG. l3 DFIG windfarm with similar WTG parameters and its Full Agg, Zone Agg and WD Agg models;

(19) FIG. 16 depicts PCC voltage, phase A current and active and reactive power for FIG. l3 DFIG windfarm with similar WTG parameters, its Full Agg and Zone Agg models with modified controllers and its WD Agg model;

(20) FIG. 17 depicts PCC voltage, phase A current and active and reactive power for FIG. l3 DFIG windfarm with unequal WTG parameters, its WD Agg model and its modified Full Agg and Zone Agg models;

(21) FIG. 18 depicts 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

(22) FIG. 19 depicts PCC voltage, phase A current and active and reactive power for FIG. l3 fixed-speed windfarm with unequal WTG parameters, its modified Full Agg, Zone Agg and WD Agg models.

DETAILED DESCRIPTION OF THE INVENTION

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

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

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

(26) 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

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

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

(29) 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.5ρπr.sup.2V.sub.W.sup.3 is 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.

(30) $\begin{array}{cc}{C}_p\left(\lambda \right)=\frac{C_{\mathrm{pm}}}{{\lambda }_{\mathrm{opt}}^2}\left(2{\lambda }_{\mathrm{opt}}-\lambda \right)\lambda ,& \mathrm{Equation}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.

(31) $\begin{array}{cc}{C}_p\left({\lambda }^{\prime }\right)=\frac{C_{\mathrm{pm}}}{{\lambda }_{\mathrm{opt}}^{{\prime }^2}}\left(2{\lambda }_{\mathrm{opt}}^{\prime }-{\lambda }^{\prime }\right){\lambda }^{\prime },& \mathrm{Equation}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.s)/Δ.sub.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.2−X.sub.ssX.sub.rr)].sup.2+[R.sub.rX.sub.ss+s.sub.0R.sub.sX.sub.rr].sup.2.

(32) 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.b−w.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.

(33) 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 I 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.

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

(35) $\begin{array}{cc}{V}_W=\frac{1}{n}{\left(\underset{k=1}{\overset{n}{.Math.}}{V}_{W_k}^3\right)}^{\frac{1}{3}},& \mathrm{Equation}5\end{array}$
where V.sub.W is the wind speed that provides a power equals to the total power of wind farm.

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

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

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

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

(40) 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

(41) $P_{m_{k_0}}=C_{p_k}\left({\lambda }_{k_0}^{\prime },{\beta }_{k_0}\right)P_{w_{k_0}}$

(42) 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.0≈rw.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=1−s.sub.0 can be substituted in Equation 4 and solved it for s.sub.0.

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

(44) $\begin{array}{cc}\Delta T_m=X_m{i}_{\mathrm{dr}0}\Delta i_{\mathrm{qs}}-X_m{i}_{\mathrm{qr}0}\Delta i_{\mathrm{ds}}-X_m{i}_{\mathrm{ds}0}\Delta i_{\mathrm{qr}}++X_m{i}_{\mathrm{qs}0}\Delta i_{\mathrm{dr}}-2H\frac{d{\mathrm{\Delta \omega }}_r}{\mathrm{dt}}.& \mathrm{Equation}6\end{array}$

(45) For a WTG and based on Equations 1-3, one obtains Equation 7.

(46) $\begin{array}{cc}{{\omega }_{r0}\Delta T_m+T_{m0}{\mathrm{\Delta \omega }}_r=P_{W_0}\frac{\partial C_p}{\partial {\lambda }^{\prime }}.Math.}_{{\lambda }^{\prime }={\lambda }_0^{\prime }}{\mathrm{\Delta \lambda }}^{\prime }+C_{p0}\Delta P_W,& \mathrm{Equation}7\end{array}$
where Δλ′=rΔW.sub.r/V.sub.W0. Solving Equation 7 for ΔT.sub.m yields Equation 8.

(47) $\begin{array}{cc}\Delta T_m=e_m{\mathrm{\Delta \omega }}_r,e_m=\frac{1}{{\omega }_{r_0}}{(\frac{\partial C_p}{\partial {\lambda }^{\prime }}.Math.}_{{\lambda }^{\prime }={\lambda }_0^{\prime }}\frac{{\mathrm{rP}}_{W_0}}{V_{W_0}}-T_{m_0}).\mathrm{where}.& \mathrm{Equation}8\end{array}$

(48) The linearized voltage equations of the machine in frequency domain are shown in Equations 9-12.

(49) $\begin{array}{cc}\Delta v_{\mathrm{qs}}=A_{\mathrm{qs}}\Delta i_{\mathrm{qs}}+\underset{-}{B_{\mathrm{qs}}}\Delta i_{\mathrm{ds}}+C_{\mathrm{qs}}\Delta i_{\mathrm{qr}}+D_{\mathrm{qs}}\Delta i_{\mathrm{dr}},& \mathrm{Equation}9\end{array}$$\begin{array}{cc}\Delta v_{\mathrm{ds}}=A_{\mathrm{ds}}\Delta i_{\mathrm{qs}}+B_{\mathrm{ds}}\Delta i_{\mathrm{ds}}+C_{\mathrm{ds}}\Delta i_{\mathrm{qr}}+D_{\mathrm{ds}}\Delta i_{\mathrm{dr}},& \mathrm{Equation}10\end{array}$$\begin{array}{cc}\Delta v_{\mathrm{qr}}=0=A_{\mathrm{qr}}\Delta i_{\mathrm{qs}}+B_{\mathrm{qr}}\Delta i_{\mathrm{ds}}+C_{\mathrm{qr}}\Delta i_{\mathrm{qr}}+D_{\mathrm{qr}}\Delta i_{\mathrm{dr}},& \mathrm{Equation}11\end{array}$$\begin{array}{cc}\Delta v_{\mathrm{dr}}=0=A_{\mathrm{dr}}\Delta i_{\mathrm{qs}}+B_{\mathrm{dr}}\Delta i_{\mathrm{ds}}+C_{\mathrm{dr}}\Delta i_{\mathrm{qr}}+D_{\mathrm{dr}}\Delta i_{\mathrm{dr}},& \mathrm{Equation}12\end{array}$$\begin{array}{cc}{A}_{\mathrm{qs}}=B_{\mathrm{ds}}=R_s+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{ss}},B_{\mathrm{qs}}=-A_{\mathrm{ds}}=X_{\mathrm{ss}},C_{\mathrm{qs}}=D_{\mathrm{ds}}=\frac{j\omega }{{\omega }_b}{X}_m,D_{\mathrm{qs}}=-C_{\mathrm{ds}}=X_m,& \end{array}$$\begin{array}{cc}\mathrm{where}\mathrm{and}& \end{array}$$\begin{array}{cc}{A}_{\mathrm{qr}}=\frac{j\omega }{{\omega }_b}{X}_m-X_d{i}_{\mathrm{dr}0},A_{\mathrm{dr}}=-s_0{X}_m+X_q{i}_{\mathrm{dr}0},B_{\mathrm{qr}}=s_0{X}_m+X_d{i}_{\mathrm{qr}0},B_{\mathrm{dr}}=\frac{j\omega }{{\omega }_b}{X}_m-X_q{i}_{\mathrm{qr}0},C_{\mathrm{qr}}=R_r+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{rr}}+X_d{i}_{\mathrm{ds}0},C_{\mathrm{dr}}=-s_0{X}_{\mathrm{rr}}-X_q{i}_{\mathrm{ds}0},D_{\mathrm{qr}}=s_0{X}_{\mathrm{rr}}-X_d{i}_{\mathrm{qs}0},D_{\mathrm{dr}}=R_r+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{rr}}+X_q{i}_{\mathrm{qs}0}.& \end{array}$

(50) 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:

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

(52) $\begin{array}{cc}\Delta v_{q{s}_k}=K_k\left(j\omega \right)\Delta i_{q{s}_k}+G_k\left(j\omega \right)\Delta v_{d{s}_k},& \mathrm{Equation}15\end{array}$$\mathrm{where}$$K_k\left(j\omega \right)={\alpha }_{q_k}\left(j\omega \right)-\frac{{\beta }_{q_k}\left(j\omega \right){\alpha }_{d_k}\left(j\omega \right)}{{\beta }_{d_k}\left(j\omega \right)},G_k\left(j\omega \right)={\beta }_{q_k}\left(j\omega \right)/{\beta }_{d_k}\left(j\omega \right).$

(53) 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.qsk≅v.sub.qs and v.sub.dsk≅v.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.

(54) 0$\begin{array}{cc}n\Delta v_{\mathrm{qs}}=\Delta i_{\mathrm{qs}}\underset{k=1}{\overset{n}{.Math.}}\frac{K_k\left(j\omega \right)\Delta i_{{\mathrm{qs}}_k}}{\Delta i_{\mathrm{qs}}}+\Delta v_{\mathrm{ds}}\underset{k=1}{\overset{n}{.Math.}}{G}_k\left(j\omega \right),& \mathrm{Equation}16\end{array}$

(55) The impact factor is defined as u.sub.k ≙Δi.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.

(56) $\begin{array}{cc}\Delta v_{\mathrm{qs}}=\frac{K^{\prime }\left(j\omega \right)}{n}\Delta i_{\mathrm{qs}}+\frac{G^{\prime }\left(j\omega \right)}{n}\Delta v_{\mathrm{ds}},\mathrm{where}{K}^{\prime }\left(j\omega \right)\stackrel{\Delta }{=}\underset{k=1}{\overset{n}{.Math.}}{K}_k{u}_k,G^{\prime }\left(j\omega \right)\stackrel{\Delta }{=}\underset{k=1}{\overset{n}{.Math.}}{G}_k.& \mathrm{Equation}17\end{array}$

(57) 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:

(58) 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, 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,

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

(60) 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

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

(62) $\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}-{\omega }_s{C}_{C_k}\right).& \mathrm{Equation}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.

(63) $\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}-{\omega }_s{C}_C\right)==\underset{k=1}{\overset{n}{.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}-{\omega }_s{C}_{C_k}\right).& \mathrm{Equation}20\end{array}$

(64) Considering a low loss collector system, V.sub.sk≅V.sub.s for k=1, 2, . . . , n. Thus, Equation 20 yields Equation 21.

(65) $\begin{array}{cc}{C}_C=\frac{1}{{\omega }_s{X}_m}+\frac{\left(X_{\mathrm{ls}}+X_{\mathrm{lr}}\right)/{\omega }_s}{{\left(\frac{R_r}{s}\right)}^2+{\left(X_{\mathrm{ls}}+X_{\mathrm{lr}}\right)}^2}--\underset{k=1}{\overset{n}{.Math.}}\left(\frac{1}{{\omega }_s{X}_{m_k}}+\frac{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)/{\omega }_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}21\end{array}$
To determine the parameters of equivalent wind turbine for a wind farm, assuming,

(66) $P_w=\underset{k=1}{\overset{n}{.Math.}}{P}_{w_k}$

(67) Equation 22 is obtained.

(68) $\begin{array}{cc}{P}_W=\frac{1}{2}\rho {\mathrm{AV}}_W^3=\underset{k=1}{\overset{n}{.Math.}}{P}_{W_k}=\underset{k=1}{\overset{n}{.Math.}}\frac{1}{2}\rho A_k{V}_{W_k}^3,& \mathrm{Equation}22\end{array}$
where

(69) $A=\stackrel{n}{\underset{k=1}{.Math.}}{A}_k$
is the equivalent surface of all WTGs. Thus, the equivalent radius, r, and wind speed, Vw, can be expressed as Equation 23.

(70) $\begin{array}{cc}r={\left(\underset{k=1}{\overset{n}{.Math.}}{r}_k^2\right)}^{\frac{1}{2}},V_W={\left(\underset{k=1}{\overset{n}{.Math.}}{A}_k{V}_{W_k}^3/\underset{k=1}{\overset{n}{.Math.}}{A}_k\right)}^{\frac{1}{3}}.& \mathrm{Equation}23\end{array}$

(71) The equivalent mechanical power,
P.sub.m=Σ.sub.k=1.sup.nP.sub.m.sub.k,and
incremental power
ΔP.sub.m=Σ.sub.k=1.sup.nΔP.sub.m.sub.k,

(72) , equations yield Equation 24 and Equation 25.

(73) $\begin{array}{cc}{C}_p\left({\lambda }^{\prime }\right)P_W=\underset{k=1}{\overset{n}{.Math.}}{C}_{p_k}\left({\lambda }_k^{\prime }\right)P_{W_k}& \mathrm{Equation}24\end{array}$$\begin{array}{cc}{P}_W\frac{\partial C_p}{\partial {\lambda }^{\prime }}\frac{r}{V_W}\Delta {\omega }_r=\underset{k=1}{\overset{n}{.Math.}}{P}_{W_k}\frac{\partial C_{p_k}}{\partial {\lambda }_k^{\prime }}\frac{r_k}{V_{W_k}}\Delta {\omega }_{r_k}.& \mathrm{Equation}25\end{array}$

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

(75) 0$\begin{array}{cc}\frac{{\lambda }_{\mathrm{opt}}}{r}\stackrel{\Delta }{=}\frac{1}{n}\underset{k=1}{\overset{n}{.Math.}}\frac{{{\lambda }_{\mathrm{opt}}}_k}{r_k},G\stackrel{\Delta }{=}\frac{{\lambda }_{\mathrm{opt}}^{\prime }}{{\lambda }_{\mathrm{opt}}}.& \mathrm{Equation}26\end{array}$

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

(77) 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 1 and the details of test scenarios are as follows: A. All WTG units are Type I with the same wind speeds, V.sub.Wk=20 m/s for 1≤k≤80; B. All WTG units are Type I with different wind speeds at

(78) $V_{W_k}=\{\begin{array}{cc}20& 1\le k\le 20\\ 16& 21\le k\le 40\\ 12& 41\le k\le 60\\ 8& 61\le k\le 80\end{array};$ 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 1≤k≤80; D. The first half of WTGs are Type I, the second half are Type II, and the wind speed of WTGs are:

(79) $V_{W_k}=\{\begin{array}{cc}20& 1\le k\le 10\mathrm{or}41\le k\le 50\\ 16& 11\le k\le 20\mathrm{or}51\le k\le 60\\ 12& 21\le k\le 30\mathrm{or}61\le k\le 12\\ 8& 31\le k\le 40\mathrm{or}71\le k\le 80\end{array}.$

(80) 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 —

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

(82) 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: 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; 2) Steady state error of the aggregated methods at their operating points, given by Equation 27.

(83) $\begin{array}{cc}{e}_{\mathrm{ss}}=.Math.\frac{{\hat{f}}_{\mathrm{ss}}-f_{\mathrm{ss}}}{f_{\mathrm{ss}}}.Math..& \mathrm{Equation}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; 3) Transient error of the aggregated methods at a specific operating point following to a small signal disturbance that is defined as Equation 28.

(84) $\begin{array}{cc}{e}_{\mathrm{Trans}}=.Math.\frac{\frac{1}{T}{\int }_{t_0}^{T+t_0}\left(\hat{f}\left(t\right)-f\left(t\right)\right)\mathrm{dt}}{f_{\mathrm{ss}}\left(t_0^{-}\right)}.Math.,& \mathrm{Equation}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.

(85) TNSTE is defined as the summation of these three normalized components as Equation 29.
TNSTE=STE+e.sub.ss+e.sub.Trans.  Equation 29:

(86) 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: 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); 2) The accuracy of Full Agg. method is highly sensitive to wind speed variations between the zones (FIGS. 6B and 6D); 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; 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).

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

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

(89) FIGS. 7A-7D show details of TNSTE components using color coded bars for each test scenario. The conclusion of test results are as follows: 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; 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); 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); 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; 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.

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

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

(92) Furthermore, Xq and Xd in Equations 9-12 are:

(93) $\begin{array}{cc}{X}_q=\left(X_m{i}_{\mathrm{qs}0}+X_{\mathrm{rr}}{i}_{\mathrm{qr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right),& \mathrm{Equation}30\end{array}$$\begin{array}{cc}{X}_d=\left(X_m{i}_{\mathrm{ds}0}+X_{\mathrm{rr}}{i}_{\mathrm{dr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right),& \mathrm{Equation}31\end{array}$$\begin{array}{cc}{\mathrm{where}{\Delta }_{m0}=2\mathrm{Hj}\omega +\frac{\partial C_p}{\partial {\lambda }^{\prime }}.Math.}_{{\lambda }_0^{\prime }}\frac{{\mathrm{rP}}_{W_0}}{V_{W0}\left(1-s_0\right)}-\frac{P_{m0}}{{\left(1-s_0\right)}^2}.& \end{array}$

(94) α.sub.q,d and β.sub.q,d in Equations 13 and 14 are:

(95) $\begin{array}{cc}{\alpha }_{q,d}=A_{\mathrm{qs},\mathrm{ds}}+{\gamma }_q{C}_{\mathrm{qs},\mathrm{ds}}+{\gamma }_d{D}_{\mathrm{qs},\mathrm{ds}},& \mathrm{Equation}32\end{array}$$\begin{array}{cc}{\beta }_{q,d}=A_{\mathrm{qs},\mathrm{ds}}+{\gamma }_q{C}_{\mathrm{qs},\mathrm{ds}}+{\gamma }_d{D}_{\mathrm{qs},\mathrm{ds}},& \mathrm{Equation}33\end{array}$$\begin{array}{cc}\mathrm{where}& \end{array}$$\begin{array}{cc}{\gamma }_q=\left(\frac{A_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{A}_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{C}_{\mathrm{dr}}}\right),{\gamma }_d=\left(\frac{A_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{A}_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{D}_{\mathrm{qr}}}\right),& \end{array}$$\begin{array}{cc}{\lambda }_q=\left(\frac{B_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{B}_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{C}_{\mathrm{dr}}}\right),{\lambda }_d=\left(\frac{B_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{B}_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{D}_{\mathrm{qr}}}\right)& \end{array}$

(96) The incremental speed is Δw.sub.r=C(jw)Δi.sub.qs, and C(jw) is:

(97) $C\left(j\omega \right)=\frac{X_m\left(i_{\mathrm{dr}0}-i_{\mathrm{qr}0}{C}_1+i_{\mathrm{qs}0}{C}_3-i_{\mathrm{ds}0}{C}_2\right)}{e_m+2\mathrm{Hj}\omega },\mathrm{where}{C}_1=-{\alpha }_d/{\beta }_d,C_2={\gamma }_q+{\lambda }_q{C}_1,\mathrm{and}{C}_3={\gamma }_d+{\lambda }_d{C}_1.$

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

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

(100) 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.rm.sup.[pu], and the base apparent power is obtained as S.sub.b=Σ.sub.j=1.sup.nS.sub.bj.

(101) The electrical torque of induction machine can be found by Equation 35. Equation 35:

(102) $T_e=\frac{X_m^2{R}_r{s}_0.Math.v_s.Math.}{{\Delta }_{T_e}}\mathrm{where}{\Delta }_{T_e}={\left[R_s{R}_r+s_0\left(X_m^2-X_{\mathrm{ss}}{X}_{\mathrm{rr}}\right)\right]}^2+{\left[R_r{X}_{\mathrm{ss}}+s_0{R}_s{X}_{\mathrm{rr}}\right]}^2.$

(103) Steady-state stator d-q currents are shown by Equations 36 and 37.

(104) $\begin{array}{cc}{i}_{\mathrm{qr}0}=-\frac{s_0^2{X}_m{X}_{\mathrm{rr}}}{s_0^2{X}_{\mathrm{rr}}^2+R_r^2}{i}_{\mathrm{qs}0}-\frac{s_0{X}_m{R}_r}{s_0^2{X}_{\mathrm{rr}}^2+R_r^2}{i}_{\mathrm{ds}0},i_{\mathrm{dr}0}=\frac{s_0{X}_m{R}_r}{s_0^2{X}_{\mathrm{rr}}^2+R_r^2}{i}_{\mathrm{qs}0}-\frac{s_0^2{X}_m{X}_{\mathrm{rr}}}{s_0^2{X}_{\mathrm{rr}}^2+R_r^2}{i}_{\mathrm{ds}0}& \mathrm{Equation}36\end{array}$$\begin{array}{cc}{i}_{\mathrm{qs}0}=\frac{v_{\mathrm{qs}0}-\frac{B_{q0}}{B_{d0}}{v}_{\mathrm{ds}0}}{A_{q0}-A_{d0}},i_{\mathrm{ds}0}=\frac{v_{\mathrm{ds}0}-\frac{A_{d0}}{A_{q0}}{v}_{\mathrm{qs}0}}{B_{d0}-B_{q0}}& \mathrm{Equation}37\end{array}$
where A.sub.q0, B.sub.q0, A.sub.d0 and B.sub.d0 are shown by Equation 38.

(105) Equation 38:

(106) 0$\begin{array}{cc}{A}_{q0}=R_s+\frac{s_0{X}_m^2{R}_r}{R_r^2+s_0^2{X}_{\mathrm{rr}}^2},B_{q0}=X_{\mathrm{ss}}+\frac{s_0^2{X}_m^2{X}_{\mathrm{rr}}}{R_r^2+s_0^2{X}_{\mathrm{rr}}^2}& \mathrm{Equation}38\end{array}$$\begin{array}{cc}{A}_{d0}=-X_{\mathrm{ss}}+\frac{s_0^2{X}_m^2{X}_{\mathrm{rr}}}{R_r^2+s_0^2{X}_{\mathrm{rr}}^2},B_{d0}=R_s+\frac{s_0{X}_m^2{R}_r}{R_r^2+s_0^2{X}_{\mathrm{rr}}^2}& \end{array}$

(107) Second-order equation of s.sub.0 is

(108) ${\alpha }_{2}s\frac{2}{0}-{\alpha }_1{s}_0+{\alpha }_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.r−P.sub.mR.sub.r  equation 39:

(109) Mechanical linearized equation of induction machine is shown by Equation 40.
ΔT.sub.m=X.sub.mi.sub.dr0Δi.sub.qs−X.sub.mi.sub.qr0Δi.sub.ds−X.sub.mi.sub.ds0Δi.sub.qr+X.sub.mi.sub.qs0Δi.sub.dr−2HpΔω.sub.r  Equation 40:

(110) 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.dr  Equation 41:
Δv.sub.ds=Δ.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.dr  Equation 42:
Δv.sub.qr=Δ.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.dr=0  Equation 43:
Δv.sub.dr=Δ.sub.dr(jω)Δi.sub.qs+B.sub.dr(jω)Δi.sub.ds+C.sub.dr(jω)Δi.sub.qr+D.sub.dr(jω)Δ.sub.dr=0  Equation 44:

(111) Where Δ.sub.qs(jw), B.sub.qs(jw), C.sub.qs(jw), . . . are shown by Equation 45.

(112) $\begin{array}{cc}{A}_{\mathrm{qs}}=R_s+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{ss}},B_{\mathrm{qs}}=X_{\mathrm{ss}},C_{\mathrm{qs}}=\frac{j\omega }{{\omega }_b}{X}_m,D_{\mathrm{qs}}=X_m& \mathrm{Equation}45\end{array}$$\begin{array}{cc}{A}_{\mathrm{ds}}=-X_{\mathrm{ss}},B_{\mathrm{ds}}=R_s+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{ss}},C_{\mathrm{ds}}=-X_m,D_{\mathrm{ds}}=\frac{j\omega }{{\omega }_b}{X}_m& \end{array}$$\begin{array}{cc}{A}_{\mathrm{qr}}=\frac{j\omega }{{\omega }_b}{X}_m-\left(X_m{i}_{\mathrm{ds}0}+X_{\mathrm{rr}}{i}_{\mathrm{dr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{dr}0}& \end{array}$$\begin{array}{cc}{B}_{\mathrm{qr}}=s_0{X}_m+\left(X_m{i}_{\mathrm{ds}0}+X_{\mathrm{rr}}{i}_{\mathrm{dr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{qr}0}& \end{array}$$\begin{array}{cc}{C}_{\mathrm{qr}}=R_r+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{rr}}+\left(X_m{i}_{\mathrm{ds}0}+X_{\mathrm{rr}}{i}_{\mathrm{dr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{ds}0}& \end{array}$$\begin{array}{cc}{D}_{\mathrm{qr}}=s_0{X}_{\mathrm{rr}}-\left(X_m{i}_{\mathrm{ds}0}+X_{\mathrm{rr}}{i}_{\mathrm{dr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{qs}0}& \end{array}$$\begin{array}{cc}{A}_{\mathrm{dr}}=-s_0{X}_m+\left(X_m{i}_{\mathrm{qs}0}+X_{\mathrm{rr}}{i}_{\mathrm{qr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{dr}0}& \end{array}$$\begin{array}{cc}{B}_{\mathrm{dr}}=\frac{j\omega }{{\omega }_b}{X}_m-\left(X_m{i}_{\mathrm{qs}0}+X_{\mathrm{rr}}{i}_{\mathrm{qr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{qr}0}& \end{array}$$\begin{array}{cc}{C}_{\mathrm{dr}}=-s_0{X}_{\mathrm{rr}}-\left(X_m{i}_{\mathrm{qs}0}+X_{\mathrm{rr}}{i}_{\mathrm{qr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{ds}0}& \end{array}$$\begin{array}{cc}{D}_{\mathrm{dr}}=R_r+\frac{j\omega }{{\omega }_b}{X}_{\mathrm{rr}}+\left(X_m{i}_{\mathrm{qs}0}+X_{\mathrm{rr}}{i}_{\mathrm{dr}0}\right)\left(\frac{X_m}{{\Delta }_{m0}}\right)i_{\mathrm{qs}0}& \end{array}$$\begin{array}{cc}\mathrm{where}& \end{array}$$\begin{array}{cc}{\Delta }_{m0}=2\mathrm{Hj}\omega +\frac{\partial C_p}{\partial \lambda }\frac{{\mathrm{rP}}_W}{\left(G\right)V_W\left(1-s_0\right)}-\frac{P_{m0}}{{\left(1-s_0\right)}^2}& \end{array}$

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

(114) $\begin{array}{cc}{\alpha }_q=A_{\mathrm{qs}}+C_{\mathrm{qs}}\left(\frac{A_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{A}_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{C}_{\mathrm{dr}}}\right)+D_{\mathrm{qs}}\left(\frac{A_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{A}_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{D}_{\mathrm{qr}}}\right)& \mathrm{Equation}47\end{array}$$\begin{array}{cc}{\beta }_q=B_{\mathrm{qs}}+C_{\mathrm{qs}}\left(\frac{B_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{B}_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{C}_{\mathrm{dr}}}\right)+D_{\mathrm{qs}}\left(\frac{B_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{B}_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{D}_{\mathrm{qr}}}\right)& \end{array}$$\begin{array}{cc}{\alpha }_d=A_{\mathrm{ds}}+C_{\mathrm{ds}}\left(\frac{A_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{A}_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{C}_{\mathrm{dr}}}\right)+D_{\mathrm{ds}}\left(\frac{A_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{A}_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{D}_{\mathrm{qr}}}\right)& \end{array}$$\begin{array}{cc}{\beta }_d=B_{\mathrm{ds}}+C_{\mathrm{ds}}\left(\frac{B_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{B}_{\mathrm{dr}}}{C_{\mathrm{qr}}-\frac{D_{\mathrm{qr}}}{D_{\mathrm{dr}}}{C}_{\mathrm{dr}}}\right)+D_{\mathrm{ds}}\left(\frac{B_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{B}_{\mathrm{qr}}}{D_{\mathrm{dr}}-\frac{C_{\mathrm{dr}}}{C_{\mathrm{qr}}}{D}_{\mathrm{qr}}}\right)& \end{array}$

(115) Therefore the electrical impact factors can be found using Equation 48.

(116) $\begin{array}{cc}{u}_m=\frac{\Delta i_{\mathrm{qsm}}}{\Delta i_{\mathrm{qs}}}& \mathrm{Equation}48\end{array}$$\begin{array}{cc}\mathrm{where}:& \end{array}$$\begin{array}{cc}\Delta v_{\mathrm{qs}}=K\left(j\omega \right)\Delta i_{\mathrm{qs}}+G\left(j\omega \right)\Delta v_{\mathrm{ds}},& \mathrm{Equation}49\end{array}$$\begin{array}{cc}K\left(j\omega \right)={\alpha }_q\left(j\omega \right)-\frac{{\beta }_q\left(j\omega \right){\alpha }_d\left(j\omega \right)}{{\beta }_d\left(j\omega \right)}& \end{array}$$\begin{array}{cc}G\left(j\omega \right)=\frac{{\beta }_q\left(j\omega \right)}{{\beta }_d\left(j\omega \right)}& \end{array}$

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

(118) $\begin{array}{cc}C\left(j\omega \right)=\frac{X_m{i}_{\mathrm{dr}0}-X_m{i}_{\mathrm{qr}0}{C}_1\left(j\omega \right)}{C_0+2\mathrm{Hj}\omega }+\frac{-X_m{i}_{\mathrm{ds}0}{C}_2\left(j\omega \right)+X_m{i}_{\mathrm{qs}0}{C}_3\left(j\omega \right)}{C_0+2\mathrm{Hj}\omega }& \mathrm{Equation}50\end{array}$$\begin{array}{cc}\mathrm{where},& \end{array}$$\begin{array}{cc}{C}_0=\left(\frac{\partial C_p\left(\lambda ,\beta \right)}{\partial \lambda }\frac{{\mathrm{rP}}_W}{\left(G\right)V_W}-T_{m0}\right)\frac{1}{{\omega }_{r0}},C_1\left(j\omega \right)=-\frac{{\alpha }_d}{{\beta }_d}& \mathrm{Equation}51\end{array}$$\begin{array}{cc}{C}_2\left(j\omega \right)=\left(\frac{A_{\mathrm{qr}}\left(j\omega \right)-\frac{D_{\mathrm{qr}}\left(j\omega \right)}{D_{\mathrm{dr}}\left(j\omega \right)}{A}_{\mathrm{dr}}\left(j\omega \right)}{{\Delta }_q}\right)+\left(\frac{B_{\mathrm{qr}}\left(j\omega \right)-\frac{D_{\mathrm{qr}}\left(j\omega \right)}{D_{\mathrm{dr}}\left(j\omega \right)}{\beta }_{\mathrm{dr}}\left(j\omega \right)}{{\Delta }_q}\right)C_1\left(j\omega \right)& \end{array}$$\begin{array}{cc}{C}_3\left(j\omega \right)=\left(\frac{A_{\mathrm{dr}}\left(j\omega \right)-\frac{C_{\mathrm{dr}}\left(j\omega \right)}{C_{\mathrm{qr}}\left(j\omega \right)}{A}_{\mathrm{qr}}\left(j\omega \right)}{{\Delta }_d}\right)+\left(\frac{B_{\mathrm{dr}}\left(j\omega \right)-\frac{C_{\mathrm{dr}}\left(j\omega \right)}{C_{\mathrm{qr}}\left(j\omega \right)}{\beta }_{\mathrm{qr}}\left(j\omega \right)}{{\Delta }_d}\right)C_1\left(j\omega \right)& \end{array}$

(119) Therefore, the mechanical impact factors can be found as Equation 52.

(120) $\begin{array}{cc}{y}_k=\frac{\Delta {\omega }_{r_k}}{\Delta {\omega }_{r_{\mathrm{eq}}}}=\frac{C_k\left(j\omega \right)K\left(j\omega \right)}{K_k\left(j\omega \right)C\left(j\omega \right)}& \mathrm{Equation}52\end{array}$

(121) Equivalent system mechanical relations are shown using Equations 53 and 54.

(122) $\begin{array}{cc}\underset{i=1}{\overset{n}{.Math.}}{C}_{\mathrm{pi}}{P}_{\mathrm{Wi}}=C_{\mathrm{peg}}{P}_{\mathrm{Weq}},\mathrm{where}{P}_{W_{\mathrm{eq}}}=\frac{1}{2}\rho A_{\mathrm{eq}}{V}_{\mathrm{Weq}}^3& \mathrm{Equation}53\end{array}$$\begin{array}{cc}\underset{i=1}{\overset{n}{.Math.}}{P}_{\mathrm{Wi}}\frac{\partial C_{\mathrm{pi}}}{\partial {\lambda }_i^{\prime }}\frac{r_i{y}_i}{V_{\mathrm{Wi}}}=P_{W_{\mathrm{eq}}}\frac{\partial C_{\mathrm{peq}}}{\partial {\lambda }_{\mathrm{eq}}^{\prime }}\frac{r_{\mathrm{eq}}}{V_{\mathrm{Weq}}}& \mathrm{Equation}54\end{array}$$\begin{array}{cc}\mathrm{where}:& \end{array}$$\begin{array}{cc}{\lambda }^{\prime }=G\lambda ,A_{\mathrm{eq}}=\underset{k=1}{\overset{n}{.Math.}}{A}_k,r_{\mathrm{eq}}=\sqrt{A_{\mathrm{eq}}/\pi }& \mathrm{Equation}55\end{array}$$\begin{array}{cc}{C}_{p_{\mathrm{eq}}}=C_{{\mathrm{pm}}_{\mathrm{eq}}}f\left({\lambda }^{\prime }\right),G_{\mathrm{eq}}=\frac{{\lambda }_{{\mathrm{opt}}_{\mathrm{eq}}}^{\prime }}{{\lambda }_{{\mathrm{opt}}_{\mathrm{eq}}}},& \end{array}$$\begin{array}{cc}{\lambda }_{{\mathrm{opt}}_{\mathrm{eq}}}=\frac{r}{n}\underset{k=1}{\overset{n}{.Math.}}\frac{{\lambda }_{{\mathrm{opt}}_k}}{r_k}& \end{array}$

(123) λ′.sub.opt.sub.ca and C.sub.pm.sub.cq can be found by help of Equations 53 and 54. Equivalent collector parameters can be found by Equations 56 and 57.

(124) $\begin{array}{cc}{R}_{C_{\mathrm{eq}}}^{\left[\mathrm{pu}\right]}=\underset{k=1}{\overset{n}{.Math.}}{R}_{C_k}^{\left[\mathrm{pu}\right]}{u}_k& \mathrm{Equation}56\end{array}$$\begin{array}{cc}{C}_{C_{\mathrm{eq}}}=\frac{1}{{\omega }_s{X}_{m_{\mathrm{eq}}}}+\frac{\left(X_{{\mathrm{ls}}_{\mathrm{eq}}}+X_{{\mathrm{lr}}_{\mathrm{eq}}}\right)/{\omega }_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.}}\left(\frac{1}{{\omega }_s{X}_{m_k}}+\frac{\left(X_{{\mathrm{ls}}_k}+X_{{\mathrm{lr}}_k}\right)/{\omega }_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}57\end{array}$

(125) Other required equations are:

(126) $\begin{array}{cc}{V}_s=V_b\sqrt{\frac{R_s^2+{\omega }_e^2{L}_{\mathrm{ss}}^2}{R_s^2+{\omega }_b^2{L}_{\mathrm{ss}}^2}}& \mathrm{Equation}58\end{array}$$\begin{array}{cc}{\omega }_e=\frac{{\omega }_m+\sqrt{\max \left(0,{\omega }_m^2+X_{\mathrm{corr}}\right)}}{2}& \mathrm{Equation}59\end{array}$$\begin{array}{cc}{X}_{\mathrm{corr}}=3P\left(v_{\mathrm{qs}}{i}_{\mathrm{qs}}-2R_s{I}_s^2\right)/K_{\mathrm{tv}}& \mathrm{Equation}60\end{array}$$\begin{array}{cc}{K}_{\mathrm{tv}}=\frac{3\left(\frac{P}{2}\right)L_M^2{R}_r{V}_b^2}{R_r^2\left(R_s^2\right)+{\omega }_b^2{L}_{\mathrm{ss}}^2}& \mathrm{Equation}61\end{array}$

(127) While at least one exemplary embodiment has been presented in the foregoing description and attached appendix, it should be appreciated that a vast number of variations exist. It should also be appreciated that the exemplary embodiment or exemplary embodiments are only examples, and are not intended to limit the scope, applicability, or configuration of the described embodiments in any way. Rather, the foregoing description and incorporated references will provide those skilled in the art with a convenient roadmap for implementing the exemplary embodiment or exemplary embodiments. It should be understood that various changes may be made in the function and arrangement of elements without departing from the scope as set forth in the appended claims and the legal equivalents thereof.