CN103248049A - Method of obtaining electrical power system voltage stability domain tangent plane containing DFIG (Doubly Fed Induction Generator) wind power plant - Google Patents

Abstract

The invention discloses a method of obtaining an electrical power system voltage stability domain tangent plane containing a DFIG (Doubly Fed Induction Generator) wind power plant and relates to the field of electrical power systems. The method comprises the following steps: after the DFIG wind power plant is subjected to mathematical modeling, the state parameters of the wind power plant are restrained; a quiescent voltage stable boundary track is determined by a prediction-correction continuous tide algorithm; when the continuous tide algorithm is applied, jacobian matrixes in the prediction step and the correction step are modified; the constraint functional expression of the DFIG wind power plant is added into the jacobian matrixes; and finally, the expression of security domain tangent plane of the electrical power system containing a DFIG is constructed through the replacement of a parameter in linear relation with a two dimensional parameter and the calculation of a point-normal equation of the tangent plane, so that the system voltage stability can be intuitively and reliably analyzed after the DFIG wind power plant is fully injected, and the data support is provided for the power grid quiescent voltage stability analyzed by a security domain.

Description

The acquisition methods that contains DFIG wind energy turbine set voltage stability domain of electric power system section

Technical field

The invention belongs to power system analysis and planning field, relate to the acquisition methods in DFIG wind energy turbine set, section, power system steady state voltage stability territory.

Background technology

Wind power generation is the important form of Wind Power Utilization, and wind energy is the energy renewable, pollution-free, that energy is big, prospect is wide, and greatly developing clean energy resource is the strategic choice of countries in the world.Wind has randomness and uncertainty, and this just makes the electrical network wind-powered electricity generation of dissolving on a large scale have certain difficulty.And the method in application safety territory is applied in the electric power system that contains wind energy turbine set, can provide margin of safety and optimal control information by tracing system operating point and security domain border, thereby make the online actual time safety supervision of electric power system, defend and control more science and more effective.Document " the method research commentary of power system security territory " has been summarized security domain classification and hierarchical relational, and the related application of security domain method.Document " Calculating Steady-State Operating Conditions for Doubly-Fed Induction Generator Wind Turbines " has proposed the foundation of a kind of double fed induction generators (Doubly Fed Induction Generator is hereinafter to be referred as DFIG) model.At present, wind energy turbine set based on DFIG has begun to generate electricity by way of merging two or more grid systems, but inject the very few of back power system voltage stabilization based on security domain section method research DFIG wind energy turbine set, and electrical network has high requirements to injecting the back power system voltage stabilization, so there is realistic meaning in research section, DFIG wind energy turbine set power system steady state voltage stability territory to the early warning power grid security.

Summary of the invention

(1) technical problem that will solve

After the objective of the invention is to solve extensive DFIG wind-powered electricity generation unit injection, the problem of power system voltage stabilization Journal of Sex Research.

(2) technical scheme

In the reality, the employed generator machine of most wind energy turbine set unanimity is used DFIG wind-powered electricity generation unit as this wind energy turbine set, and then this wind energy turbine set will be used same type wind-powered electricity generation unit, and can not use with other type Blast Furnace Top Gas Recovery Turbine Unit (TRT).The present invention proposes the method that find the solution in a kind of DFIG of containing power system steady state voltage stability security domain section, this method is at first carried out mathematical modeling to DFIG, utilizes seven DFIG state variables that model is retrained; Utilize prediction and the continuous tide method of proofreading and correct that the two-dimentional border of power system steady state voltage stability is followed the trail of then, obtain a static voltage stability boundary locus, for determining that containing DFIG power system steady state voltage stability security domain section lays the first stone; At last by to the replacement of the linear parameter of this two-dimensional parameter and the calculating of section point French equation, construct the expression formula that this contains section, DFIG power system security territory.

(3) beneficial effect

After the invention solves extensive DFIG wind-powered electricity generation unit injection, the problem of power system voltage stabilization Journal of Sex Research.

Description of drawings

Fig. 1 is stable region of the present invention section acquisition methods calculation process schematic diagram.

Embodiment

Below in conjunction with drawings and Examples, the specific embodiment of the present invention is described further.Following examples only are used for explanation the present invention, but are not used for limiting the scope of the invention.

Calculation process as shown in Figure 1.

The present invention includes following steps:

1) determines the mathematical models of power system of DFIG wind energy turbine set

1.1) at first think all double-fed fan motor units be same model and wind energy turbine set around wind speed evenly distribute, so can be with a double-fed blower fan model representation wind energy turbine set.Wherein, the mechanical energy of double-fed fan motor machine output expression formula is

$P_m=\frac{1}{2}\mathrm{\ρ}{\mathrm{AC}}_p{U}^3---\left(1\right)$

Wherein, ρ is atmospheric density; A is the swept area of wind energy conversion system; U is wind speed; C _pBe the power coefficient of wind energy conversion system, show the useful wind energy ratio that wind wheel obtains from wind, with tip-speed ratio T _TsrRelevant.

1.2) following contain and join the power flow equation formula for containing the electric power system of DFIG wind energy turbine set:

$f_1=(1+\mathrm{\λ})P_i+P_{\mathrm{im}}-V_i\underset{j\∈i}{\mathrm{\Σ}}{V}_j(G_{\mathrm{ij}}\mathrm{ocs}{\mathrm{\δ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\δ}}_{\mathrm{ij}})=0$

$f_2=(1+\mathrm{\λ})Q_i+Q_{\mathrm{im}}-V_i\underset{j\∈i}{\mathrm{\Σ}}{V}_j(G_{\mathrm{ij}}\sin {\mathrm{\δ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\cos {\mathrm{\δ}}_{\mathrm{ij}})=0$

f ₃=V _scos(θ _s)+R _sI _scos(φ _s)-(X _s+X _m)I _ssin(φ _s)+X _mI _rsin(φ _r)=0

f ₄=V _ssin(θ _s)-R _sI _ssin(φ _s)+(X _s+X _m)I _scos(φ _s)-X _mI _rcos(φ _r)=0

f ₅=V _rcos(θ _r)-sX _mI _ssin(φ _s)-R _rI _rcos(φ _r)+s(X _s+X _m)I _rsin(φ _r)=0 （2）

f ₆=V _rsin(θ _r)+sX _mI _scos(φ _s)-R _rI _rsin(φ _r)-s(X _s+X _m)I _rcos(φ _r)=0

f ₇=P _m-V _sI _scos(θ _s-φ _s)+V _rI _rcos(θ _r-φ _r)=0

f ₈=Q _m-V _sI _ssin(θ _s-φ _s)=0

$f_9=P_m+I_s^2{R}_s+I_r^2{R}_r+P_{\mathrm{GB},\mathrm{Loss}}-k{(1-s)}^3=0$

In the formula, V _iBe node voltage size, δ _iBe node voltage phase angle, V _sBe stator voltage size, θ _sBe stator voltage phase angle, V _rBe rotor voltage size, θ _rBe rotor voltage phase angle, I _rBe rotor current size, φ _rBe rotor current phase angle, I _sBe stator current size, φ _sBe the stator current phase angle, s is slippage.

To contain V _r, θ _r, I _S, φ _S, I _r, φ _r, seven formulas of seven unknown numbers of s, i.e. f in the formula (2) ₃-f ₉Be added to system contain in the ginseng power flow equation after, its Jacobian matrix is revised, the calculating of predicting-proofreading and correct is to obtain containing the static neutrality point of DFIG electric power system.

2) contain DFIG power system steady state voltage stability border initial point with the calculating of continuous tide method

After introducing the double-fed blower fan, system contains the ginseng tide model and can be designated as:

f(V _i,δ _i,V _r,θ _r,I _s,φ _s,I _r,φ _r,s,λ,P _m)=0 (3)

In the formula, V _i, δ _i, V _r, θ _r, I _s, φ _s, I _r, φ _r, s is the state vector of system; λ is that reflection system generator is meritorious, reactive power is injected and the vector of load variations; p _mFor the wind energy turbine set mechanical output is injected vector (being system's variable element vector).

Elder generation is p fixedly _mConstant, λ is free running parameter, depicts the P-U curve of an active power and voltage again by the continuous tide method, and draws a saddle joint bifurcation (Saddle Node Bifurcation is hereinafter to be referred as the SNB point), is the initial boundary point.Below be the continuous tide computational methods.

2.1) calculating in prediction step

At first to determine the prediction direction of variable.Be example with the tangential direction, formula (3) is differentiated,

$f_{V_i}{d}_{V_i}+f_{{\mathrm{\δ}}_i}{d}_{{\mathrm{\δ}}_i}+f_{V_r}{d}_{V_r}+f_{{\mathrm{\θ}}_r}{d}_{{\mathrm{\θ}}_r}+f_{I_s}{d}_{I_s}+f_{{\mathrm{\φ}}_s}{d}_{{\mathrm{\φ}}_s}$

（4）

$+f_{I_r}{d}_{I_r}+f_{{\mathrm{\φ}}_r}{d}_{{\mathrm{\φ}}_r}+f_s{d}_s+f_s{d}_s+f_{\mathrm{\λ}}{d}_{\mathrm{\λ}}=0$

If current solution is in the ordinary position of solution curve, namely Jacobian matrix is nonsingular, then with λ as parameterized variables, the change direction of λ is+1 before critical point, the change direction of λ is-1 after critical point.Can get following formula

$\left[\begin{array}{cccccccccc}{f}_{{1V}_i}& f_{{1\mathrm{\δ}}_i}& f_{{1V}_r}& f_{{1\mathrm{\θ}}_r}& f_{{1I}_S}& f_{{1\mathrm{\φ}}_S}& f_{{1I}_r}& f_{{1\mathrm{\φ}}_r}& f_{1s}& f_{1\mathrm{\λ}}\\ f_{{2V}_i}& f_{{2\mathrm{\δ}}_i}& f_{{2V}_r}& f_{{2\mathrm{\θ}}_r}& f_{{2I}_S}& f_{{2\mathrm{\φ}}_S}& f_{{2I}_r}& f_{{2\mathrm{\φ}}_r}& f_{2s}& f_{2\mathrm{\λ}}\\ f_{{3V}_i}& f_{{3\mathrm{\δ}}_i}& f_{{3V}_r}& f_{{3\mathrm{\θ}}_r}& f_{{3I}_S}& f_{{3\mathrm{\φ}}_S}& f_{{3I}_r}& f_{{3\mathrm{\φ}}_r}& f_{3s}& f_{3\mathrm{\λ}}\\ f_{{4V}_i}& f_{{4\mathrm{\δ}}_i}& f_{{4V}_r}& f_{{4\mathrm{\θ}}_r}& f_{{4I}_S}& f_{{4\mathrm{\φ}}_S}& f_{{4I}_r}& f_{{4\mathrm{\φ}}_r}& f_{4s}& f_{4\mathrm{\λ}}\\ f_{{5V}_i}& f_{{5\mathrm{\δ}}_i}& f_{{5V}_r}& f_{{5\mathrm{\θ}}_r}& f_{{5I}_S}& f_{{5\mathrm{\φ}}_S}& f_{{5I}_r}& f_{{5\mathrm{\φ}}_r}& f_{5s}& f_{5\mathrm{\λ}}\\ f_{{6V}_i}& f_{{6\mathrm{\δ}}_i}& f_{{6V}_r}& f_{{6\mathrm{\θ}}_r}& f_{{6I}_S}& f_{{6\mathrm{\φ}}_S}& f_{{6I}_r}& f_{{6\mathrm{\φ}}_r}& f_{6s}& f_{6\mathrm{\λ}}\\ f_{{7V}_i}& f_{{7\mathrm{\δ}}_i}& f_{{7V}_r}& f_{{7\mathrm{\θ}}_r}& f_{{7I}_S}& f_{{7\mathrm{\φ}}_S}& f_{{7I}_r}& f_{{7\mathrm{\φ}}_r}& f_{7s}& f_{7\mathrm{\λ}}\\ f_{{8V}_i}& f_{{8\mathrm{\δ}}_i}& f_{{8V}_r}& f_{{8\mathrm{\θ}}_r}& f_{{8I}_S}& f_{{8\mathrm{\φ}}_S}& f_{{8I}_r}& f_{{8\mathrm{\φ}}_r}& f_{8s}& f_{8\mathrm{\λ}}\\ f_{{9V}_i}& f_{{9\mathrm{\δ}}_i}& f_{{9V}_r}& f_{{9\mathrm{\θ}}_r}& f_{{9I}_S}& f_{9{\mathrm{\φ}}_S}& f_{{9I}_r}& f_{{9\mathrm{\φ}}_r}& f_{9s}& f_{9\mathrm{\λ}}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{d}_{V_i}\\ d_{{\mathrm{\δ}}_i}\\ d_{V_r}\\ d_{{\mathrm{\θ}}_r}\\ d_{I_s}\\ d_{{\mathrm{\φ}}_s}\\ d_{I_r}\\ d_{{\mathrm{\φ}}_r}\\ d_s\\ d_{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \±1\end{array}\right]---\left(5\right)$

Can determine prediction direction thus.

If current solution is in the position of solution curve near critical point, be that Jacobian matrix is near unusual, the coefficient matrix of formula (5) also will be near unusual, then answer selection mode variable (for example voltage of the node k of rate of change maximum) as the parametrization variable quantity, and with λ as the common variation amount, the slope of state variable is as tangential direction, and can determine prediction direction by following formula this moment:

$\left[\begin{array}{cccccccccc}{f}_{{1V}_i}& f_{{1\mathrm{\δ}}_i}& f_{{1V}_r}& f_{{1\mathrm{\θ}}_r}& f_{{1I}_S}& f_{{1\mathrm{\φ}}_S}& f_{{1I}_r}& f_{{1\mathrm{\φ}}_r}& f_{1s}& f_{1\mathrm{\λ}}\\ f_{{2V}_i}& f_{{2\mathrm{\δ}}_i}& f_{{2V}_r}& f_{{2\mathrm{\θ}}_r}& f_{{2I}_S}& f_{{2\mathrm{\φ}}_S}& f_{{2I}_r}& f_{{2\mathrm{\φ}}_r}& f_{2s}& f_{2\mathrm{\λ}}\\ f_{{3V}_i}& f_{{3\mathrm{\δ}}_i}& f_{{3V}_r}& f_{{3\mathrm{\θ}}_r}& f_{{3I}_S}& f_{{3\mathrm{\φ}}_S}& f_{{3I}_r}& f_{{3\mathrm{\φ}}_r}& f_{3s}& f_{3\mathrm{\λ}}\\ f_{{4V}_i}& f_{{4\mathrm{\δ}}_i}& f_{{4V}_r}& f_{{4\mathrm{\θ}}_r}& f_{{4I}_S}& f_{{4\mathrm{\φ}}_S}& f_{{4I}_r}& f_{{4\mathrm{\φ}}_r}& f_{4s}& f_{4\mathrm{\λ}}\\ f_{{5V}_i}& f_{{5\mathrm{\δ}}_i}& f_{{5V}_r}& f_{{5\mathrm{\θ}}_r}& f_{{5I}_S}& f_{{5\mathrm{\φ}}_S}& f_{{5I}_r}& f_{{5\mathrm{\φ}}_r}& f_{5s}& f_{5\mathrm{\λ}}\\ f_{{6V}_i}& f_{{6\mathrm{\δ}}_i}& f_{{6V}_r}& f_{{6\mathrm{\θ}}_r}& f_{{6I}_S}& f_{{6\mathrm{\φ}}_S}& f_{{6I}_r}& f_{{6\mathrm{\φ}}_r}& f_{6s}& f_{6\mathrm{\λ}}\\ f_{{7V}_i}& f_{{7\mathrm{\δ}}_i}& f_{{7V}_r}& f_{{7\mathrm{\θ}}_r}& f_{{7I}_S}& f_{{7\mathrm{\φ}}_S}& f_{{7I}_r}& f_{{7\mathrm{\φ}}_r}& f_{7s}& f_{7\mathrm{\λ}}\\ f_{{8V}_i}& f_{{8\mathrm{\δ}}_i}& f_{{8V}_r}& f_{{8\mathrm{\θ}}_r}& f_{{8I}_S}& f_{{8\mathrm{\φ}}_S}& f_{{8I}_r}& f_{{8\mathrm{\φ}}_r}& f_{8s}& f_{8\mathrm{\λ}}\\ f_{{9V}_i}& f_{{9\mathrm{\δ}}_i}& f_{{9V}_r}& f_{{9\mathrm{\θ}}_r}& f_{{9I}_S}& f_{9{\mathrm{\φ}}_S}& f_{{9I}_r}& f_{{9\mathrm{\φ}}_r}& f_{9s}& f_{9\mathrm{\λ}}\\ & & & & e_k^t& & & & & 0\end{array}\right]\left[\begin{array}{c}{d}_{V_i}\\ d_{{\mathrm{\δ}}_i}\\ d_{V_r}\\ d_{{\mathrm{\θ}}_r}\\ d_{I_s}\\ d_{{\mathrm{\φ}}_s}\\ d_{I_r}\\ d_{{\mathrm{\φ}}_r}\\ d_s\\ d_{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \±1\end{array}\right]---\left(6\right)$

Wherein,

Be that k element is+1, all the other elements are 0 row vector.Should notice that formula (6) is different with the coefficient matrix of formula (5).Owing to selected state variable as parameterized variables, though Jacobian matrix near unusual (supposing that its order is n-1), the coefficient matrix that also can prove (6) is nonsingular.

According to formula (5) or the definite prediction direction of formula (6), it is as follows to calculate future position:

$\left[\begin{array}{c}{\tilde{V}}_i\\ {\widetilde{\mathrm{\δ}}}_i\\ {\tilde{V}}_r\\ {\widetilde{\mathrm{\θ}}}_r\\ {\tilde{I}}_S\\ {\widetilde{\mathrm{\φ}}}_S\\ {\tilde{I}}_r\\ {\widetilde{\mathrm{\φ}}}_r\\ \tilde{s}\\ \widetilde{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{c}{V}_{i0}\\ {\mathrm{\δ}}_{i0}\\ V_{r0}\\ {\mathrm{\θ}}_{r0}\\ I_{S0}\\ {\mathrm{\φ}}_{S0}\\ I_{r0}\\ {\mathrm{\φ}}_{r0}\\ s_0\\ {\mathrm{\λ}}_0\end{array}\right]+\mathrm{\σ}\left[\begin{array}{c}{d}_{V_i}\\ d_{{\mathrm{\δ}}_i}\\ d_{V_r}\\ d_{{\mathrm{\θ}}_r}\\ d_{I_S}\\ d_{{\mathrm{\φ}}_S}\\ d_{I_r}\\ d_{{\mathrm{\φ}}_r}\\ d_s\\ d_{\mathrm{\λ}}\end{array}\right]---\left(7\right)$

Wherein, σ is the prediction step.

2.2) proofread and correct the calculating in step

Proofreading and correct in the step, if prediction direction is obtained by formula (5), should be earlier fixing λ, adopt the vertical correction method, with Find the solution the power flow equation of formula (1) for initial value.Be example with the inferior method of newton-pressgang, the iteration form is as follows:

$\left[\begin{array}{cccccccccc}{f}_{{1V}_i}& f_{{1\mathrm{\δ}}_i}& f_{{1V}_r}& f_{{1\mathrm{\θ}}_r}& f_{{1I}_S}& f_{{1\mathrm{\φ}}_S}& f_{{1I}_r}& f_{{1\mathrm{\φ}}_r}& f_{1s}& f_{1\mathrm{\λ}}\\ f_{{2V}_i}& f_{{2\mathrm{\δ}}_i}& f_{{2V}_r}& f_{{2\mathrm{\θ}}_r}& f_{{2I}_S}& f_{{2\mathrm{\φ}}_S}& f_{{2I}_r}& f_{{2\mathrm{\φ}}_r}& f_{2s}& f_{2\mathrm{\λ}}\\ f_{{3V}_i}& f_{{3\mathrm{\δ}}_i}& f_{{3V}_r}& f_{{3\mathrm{\θ}}_r}& f_{{3I}_S}& f_{{3\mathrm{\φ}}_S}& f_{{3I}_r}& f_{{3\mathrm{\φ}}_r}& f_{3s}& f_{3\mathrm{\λ}}\\ f_{{4V}_i}& f_{{4\mathrm{\δ}}_i}& f_{{4V}_r}& f_{{4\mathrm{\θ}}_r}& f_{{4I}_S}& f_{{4\mathrm{\φ}}_S}& f_{{4I}_r}& f_{{4\mathrm{\φ}}_r}& f_{4s}& f_{4\mathrm{\λ}}\\ f_{{5V}_i}& f_{{5\mathrm{\δ}}_i}& f_{{5V}_r}& f_{{5\mathrm{\θ}}_r}& f_{{5I}_S}& f_{{5\mathrm{\φ}}_S}& f_{{5I}_r}& f_{{5\mathrm{\φ}}_r}& f_{5s}& f_{5\mathrm{\λ}}\\ f_{{6V}_i}& f_{{6\mathrm{\δ}}_i}& f_{{6V}_r}& f_{{6\mathrm{\θ}}_r}& f_{{6I}_S}& f_{{6\mathrm{\φ}}_S}& f_{{6I}_r}& f_{{6\mathrm{\φ}}_r}& f_{6s}& f_{6\mathrm{\λ}}\\ f_{{7V}_i}& f_{{7\mathrm{\δ}}_i}& f_{{7V}_r}& f_{{7\mathrm{\θ}}_r}& f_{{7I}_S}& f_{{7\mathrm{\φ}}_S}& f_{{7I}_r}& f_{{7\mathrm{\φ}}_r}& f_{7s}& f_{7\mathrm{\λ}}\\ f_{{8V}_i}& f_{{8\mathrm{\δ}}_i}& f_{{8V}_r}& f_{{8\mathrm{\θ}}_r}& f_{{8I}_S}& f_{{8\mathrm{\φ}}_S}& f_{{8I}_r}& f_{{8\mathrm{\φ}}_r}& f_{8s}& f_{8\mathrm{\λ}}\\ f_{{9V}_i}& f_{{9\mathrm{\δ}}_i}& f_{{9V}_r}& f_{{9\mathrm{\θ}}_r}& f_{{9I}_S}& f_{9{\mathrm{\φ}}_S}& f_{{9I}_r}& f_{{9\mathrm{\φ}}_r}& f_{9s}& f_{9\mathrm{\λ}}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δ}}_{V_i}\\ {\mathrm{\Δ}}_{{\mathrm{\δ}}_i}\\ {\mathrm{\Δ}}_{V_r}\\ {\mathrm{\Δ}}_{{\mathrm{\θ}}_r}\\ {\mathrm{\Δ}}_{I_s}\\ {\mathrm{\Δ}}_{{\mathrm{\φ}}_s}\\ {\mathrm{\Δ}}_{I_r}\\ {\mathrm{\Δ}}_{{\mathrm{\φ}}_r}\\ {\mathrm{\Δ}}_s\\ {\mathrm{\Δ}}_{\mathrm{\λ}}\end{array}\right]=-\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6\\ f_7\\ f_8\\ f_9\\ 0\end{array}\right]---\left(8\right)$

If above-mentioned trend is calculated convergence, then can obtain on the P-U curve a bit, begin new prediction step calculating then.If trend is dispersed, the measure of two kinds of correspondences is arranged: the first reduces step-length σ and predicts new point with formula (7), and uses formula (8) iteration again with the vertical correction method; Another kind of way be selection mode variable (for example voltage of the node k of rate of change maximum) as parameterized variables, λ as common variables, is adopted level correction method solution power flow equation, the iteration form is as follows:

$\left[\begin{array}{cccccccccc}{f}_{{1V}_i}& f_{{1\mathrm{\δ}}_i}& f_{{1V}_r}& f_{{1\mathrm{\θ}}_r}& f_{{1I}_S}& f_{{1\mathrm{\φ}}_S}& f_{{1I}_r}& f_{{1\mathrm{\φ}}_r}& f_{1s}& f_{1\mathrm{\λ}}\\ f_{{2V}_i}& f_{{2\mathrm{\δ}}_i}& f_{{2V}_r}& f_{{2\mathrm{\θ}}_r}& f_{{2I}_S}& f_{{2\mathrm{\φ}}_S}& f_{{2I}_r}& f_{{2\mathrm{\φ}}_r}& f_{2s}& f_{2\mathrm{\λ}}\\ f_{{3V}_i}& f_{{3\mathrm{\δ}}_i}& f_{{3V}_r}& f_{{3\mathrm{\θ}}_r}& f_{{3I}_S}& f_{{3\mathrm{\φ}}_S}& f_{{3I}_r}& f_{{3\mathrm{\φ}}_r}& f_{3s}& f_{3\mathrm{\λ}}\\ f_{{4V}_i}& f_{{4\mathrm{\δ}}_i}& f_{{4V}_r}& f_{{4\mathrm{\θ}}_r}& f_{{4I}_S}& f_{{4\mathrm{\φ}}_S}& f_{{4I}_r}& f_{{4\mathrm{\φ}}_r}& f_{4s}& f_{4\mathrm{\λ}}\\ f_{{5V}_i}& f_{{5\mathrm{\δ}}_i}& f_{{5V}_r}& f_{{5\mathrm{\θ}}_r}& f_{{5I}_S}& f_{{5\mathrm{\φ}}_S}& f_{{5I}_r}& f_{{5\mathrm{\φ}}_r}& f_{5s}& f_{5\mathrm{\λ}}\\ f_{{6V}_i}& f_{{6\mathrm{\δ}}_i}& f_{{6V}_r}& f_{{6\mathrm{\θ}}_r}& f_{{6I}_S}& f_{{6\mathrm{\φ}}_S}& f_{{6I}_r}& f_{{6\mathrm{\φ}}_r}& f_{6s}& f_{6\mathrm{\λ}}\\ f_{{7V}_i}& f_{{7\mathrm{\δ}}_i}& f_{{7V}_r}& f_{{7\mathrm{\θ}}_r}& f_{{7I}_S}& f_{{7\mathrm{\φ}}_S}& f_{{7I}_r}& f_{{7\mathrm{\φ}}_r}& f_{7s}& f_{7\mathrm{\λ}}\\ f_{{8V}_i}& f_{{8\mathrm{\δ}}_i}& f_{{8V}_r}& f_{{8\mathrm{\θ}}_r}& f_{{8I}_S}& f_{{8\mathrm{\φ}}_S}& f_{{8I}_r}& f_{{8\mathrm{\φ}}_r}& f_{8s}& f_{8\mathrm{\λ}}\\ f_{{9V}_i}& f_{{9\mathrm{\δ}}_i}& f_{{9V}_r}& f_{{9\mathrm{\θ}}_r}& f_{{9I}_S}& f_{9{\mathrm{\φ}}_S}& f_{{9I}_r}& f_{{9\mathrm{\φ}}_r}& f_{9s}& f_{9\mathrm{\λ}}\\ & & & & e_k^t& & & & & 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δ}}_{V_i}\\ {\mathrm{\Δ}}_{{\mathrm{\δ}}_i}\\ {\mathrm{\Δ}}_{V_r}\\ {\mathrm{\Δ}}_{{\mathrm{\θ}}_r}\\ {\mathrm{\Δ}}_{I_s}\\ {\mathrm{\Δ}}_{{\mathrm{\φ}}_s}\\ {\mathrm{\Δ}}_{I_r}\\ {\mathrm{\Δ}}_{{\mathrm{\φ}}_r}\\ {\mathrm{\Δ}}_s\\ {\mathrm{\Δ}}_{\mathrm{\λ}}\end{array}\right]=-\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6\\ f_7\\ f_8\\ f_9\\ 0\end{array}\right]---\left(9\right)$

Proofreading and correct in the step, if prediction direction itself is obtained by formula (6), also should separate power flow equation by formula (9).

3) calculate SNB point exact value

Accurately find the solution the SNB initial point, the constraint equation that SNB is ordered can be described as

The voltage stability domain boundary face that is made of the SNB point can be described as

$\left\{\begin{array}{c}f(V_i,{\mathrm{\δ}}_i,V_r,{\mathrm{\θ}}_r,I_s,{\mathrm{\φ}}_s,I_r,{\mathrm{\φ}}_r,s,\mathrm{\λ},P_m)=0\\ f_{V_i}{|}_{*}y+f_{{\mathrm{\δ}}_i}{|}_{*}y+f_{V_r}{|}_{*}y+f_{{\mathrm{\θ}}_r}{|}_{*}y+f_{I_S}{|}_{*}y+f_{{\mathrm{\φ}}_S}{|}_{*}y+f_{I_r}{|}_{*}y+f_{{\mathrm{\φ}}_r}{|}_{*}y+f_s{|}_{*}y=0\\ y^{t}y=1\end{array}\right.---\left(10\right)$

y ^t=[dV _i dδ _i dV _r dθ _r dI _s dφ _s dI _r dφ _r ds]

In the formula, V _i, δ _i, V _r, θ _r, I _S, φ _S, I _r, φ _r, s is the n dimension state vector of system; λ is that reflection system generator is meritorious, reactive power is injected and the vector of load variations; P _mInject vector (being system's variable element vector) for the wind energy turbine set mechanical output, this moment is earlier with P _mValue is fixing; Y is the right characteristic vector of power flow equation Jacobian matrix.The solution that obtains with continuous tide is initial value, carries out the solving equation variable with the inferior method of newton-pressgang.Accurately draw initial point.

4) calculate the SNB locus of points that contains DFIG power system stability territory

The algorithm of this stable region SNB locus of points is applied forecasting-bearing calibration equally, formula (10) can be noted by abridging be

φ(z)=0 （11）

At this moment

z ^t=[V _i δ _i V _r θ _r I _S φ _S I _r φ _r s λ P _m]∈R ²ⁿ⁺²,

$\mathrm{\φ}=\left[\begin{array}{c}f(V_i,{\mathrm{\δ}}_i,V_r,{\mathrm{\θ}}_r,I_s,{\mathrm{\φ}}_s,I_r,{\mathrm{\φ}}_r,s,\mathrm{\λ},P_m)\\ f_{V_i}{|}_{*}y+f_{{\mathrm{\δ}}_i}{|}_{*}y+f_{V_r}{|}_{*}y+f_{{\mathrm{\θ}}_r}{|}_{*}y+f_{I_S}{|}_{*}y+f_{{\mathrm{\φ}}_S}{|}_{*}y+f_{I_r}{|}_{*}y+f_{{\mathrm{\φ}}_r}{|}_{*}y+f_s{|}_{*}y\\ y^{t}y-1\end{array}\right]:R^{2n+2}\→R^{2n+1}$

Dimension by the determined Jacobian matrix of following formula is (2n+1) * (2n+2), therefore adds a parametrization equation

(z-z ₁) ^tv=τ （12）

In the formula, v is z ₁The tangent vector of point, z ₀Be initial point, τ is step-length.V is determined by following formula:

${\mathrm{\φ}}_z{|}_{z=z_0}v=0---\left(13\right)$

4.1) calculating in prediction step

At first suppose by step 2) calculated 1 z ₀On this power system stability territory SNB locus of points.Begin to predict from this point, can determine prediction direction by following formula:

$\left[\begin{array}{cccccccccccc}{f}_{{1V}_I}& f_{{1\mathrm{\δ}}_i}& f_{{1V}_r}& f_{{1\mathrm{\θ}}_r}& f_{{1I}_S}& f_{{1\mathrm{\φ}}_S}& f_{{1I}_r}& f_{{1\mathrm{\φ}}_r}& f_{1s}& 0& f_{1\mathrm{\λ}}& f_{{1P}_m}\\ f_{{2V}_i}& f_{{2\mathrm{\δ}}_i}& f_{{2V}_r}& f_{{2\mathrm{\θ}}_r}& f_{{2I}_S}& f_{{2\mathrm{\φ}}_S}& f_{{2I}_r}& f_{{2\mathrm{\φ}}_r}& f_{2s}& 0& f_{2\mathrm{\λ}}& f_{{2P}_m}\\ f_{{3V}_i}& f_{{3\mathrm{\δ}}_i}& f_{{3V}_r}& f_{{3\mathrm{\θ}}_r}& f_{{3I}_S}& f_{{3\mathrm{\φ}}_S}& f_{{3I}_r}& f_{{3\mathrm{\φ}}_r}& f_{3s}& 0& f_{3\mathrm{\λ}}& f_{{3P}_m}\\ f_{{4V}_i}& f_{{4\mathrm{\δ}}_i}& f_{{4V}_r}& f_{{4\mathrm{\θ}}_r}& f_{{4I}_S}& f_{{4\mathrm{\φ}}_S}& f_{{4I}_r}& f_{{4\mathrm{\φ}}_r}& f_{4s}& 0& f_{4\mathrm{\λ}}& f_{{4P}_m}\\ f_{{5V}_i}& f_{{5\mathrm{\δ}}_i}& f_{{5V}_r}& f_{{5\mathrm{\θ}}_r}& f_{{5I}_S}& f_{{5\mathrm{\φ}}_S}& f_{{5I}_r}& f_{{5\mathrm{\φ}}_r}& f_{5s}& 0& f_{5\mathrm{\λ}}& f_{{5P}_m}\\ f_{{6V}_i}& f_{{6\mathrm{\δ}}_i}& f_{{6V}_r}& f_{6{\mathrm{\θ}}_r}& f_{{6I}_S}& f_{{6\mathrm{\φ}}_S}& f_{{6I}_r}& f_{{6\mathrm{\φ}}_r}& f_{6s}& 0& f_{6\mathrm{\λ}}& f_{{6P}_m}\\ f_{{7V}_i}& f_{{7\mathrm{\δ}}_i}& f_{{7V}_r}& f_{{7\mathrm{\θ}}_r}& f_{{7I}_S}& f_{{7\mathrm{\φ}}_S}& f_{{7I}_r}& f_{{7\mathrm{\φ}}_r}& f_{7s}& 0& f_{7\mathrm{\λ}}& f_{{7P}_m}\\ f_{{8V}_i}& f_{{8\mathrm{\δ}}_i}& f_{{8V}_r}& f_{{8\mathrm{\θ}}_r}& f_{{8I}_S}& f_{{8\mathrm{\φ}}_S}& f_{{8I}_r}& f_{{8\mathrm{\φ}}_r}& f_{8s}& 0& f_{8\mathrm{\λ}}& f_{{8P}_m}\\ f_{{9V}_i}& f_{{9\mathrm{\δ}}_i}& f_{{9V}_r}& f_{{9\mathrm{\θ}}_r}& f_{{9I}_S}& f_{{9\mathrm{\φ}}_S}& f_{{9I}_r}& f_{{9\mathrm{\φ}}_r}& f_{9s}& 0& f_{9\mathrm{\λ}}& f_{{9P}_m}\\ \frac{\∂g}{{\∂V}_i}& \frac{\∂g}{{\∂\mathrm{\δ}}_i}& \frac{\∂g}{{\∂V}_r}& \frac{\∂g}{{\∂\mathrm{\θ}}_r}& \frac{\∂g}{{\∂I}_S}& \frac{\∂g}{{\∂\mathrm{\φ}}_S}& \frac{\∂g}{{\∂I}_r}& \frac{\∂g}{{\∂\mathrm{\φ}}_r}& \frac{\∂g}{\∂s}& \frac{\∂g}{\∂y}& \frac{\∂g}{\∂\mathrm{\λ}}& \frac{\∂g}{{\∂P}_m}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {2y}^t& 0& 0\\ {\mathrm{\ΔV}}_i& {\mathrm{\Δ\δ}}_i& {\mathrm{\ΔV}}_r& {\mathrm{\Δ\θ}}_r& {\mathrm{\ΔI}}_S& {\mathrm{\Δ\φ}}_S& {\mathrm{\ΔI}}_r& {\mathrm{\Δ\φ}}_r& \mathrm{\Δs}& {\mathrm{\Δy}}^t& \mathrm{\Δ\λ}& {\mathrm{\ΔP}}_m\end{array}\right]\left[\begin{array}{c}{d}_{V_i}\\ d_{{\mathrm{\δ}}_i}\\ d_{V_r}\\ d_{{\mathrm{\θ}}_r}\\ d_{I_S}\\ d_{{\mathrm{\φ}}_S}\\ d_{I_r}\\ d_{{\mathrm{\φ}}_r}\\ d_s\\ d_y\\ d_{\mathrm{\λ}}\\ d_{P_m}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \±1\end{array}\right]---\left(14\right)$

According to (14) determined prediction direction, it is as follows to calculate future position:

$\left[\begin{array}{c}{\tilde{V}}_i\\ {\widetilde{\mathrm{\δ}}}_i\\ {\tilde{V}}_r\\ {\widetilde{\mathrm{\θ}}}_r\\ {\tilde{I}}_S\\ {\mathrm{\φ}}_S\\ {\tilde{I}}_r\\ {\widetilde{\mathrm{\φ}}}_r\\ \tilde{s}\\ \tilde{y}\\ \widetilde{\mathrm{\λ}}\\ {\tilde{P}}_m\end{array}\right]=\left[\begin{array}{c}{V}_{i0}\\ {\mathrm{\δ}}_{i0}\\ V_{r0}\\ {\mathrm{\θ}}_{r0}\\ I_{S0}\\ {\mathrm{\φ}}_{S0}\\ I_{r0}\\ {\mathrm{\φ}}_{r0}\\ s_0\\ y_0\\ {\mathrm{\λ}}_0\\ P_{m0}\end{array}\right]+\mathrm{\τ}\left[\begin{array}{c}{d}_{V_i}\\ d_{{\mathrm{\δ}}_i}\\ d_{V_r}\\ d_{{\mathrm{\θ}}_r}\\ d_{I_S}\\ d_{{\mathrm{\φ}}_S}\\ d_{I_r}\\ d_{{\mathrm{\φ}}_r}\\ d_s\\ d_y\\ d_{\mathrm{\λ}}\\ d_{P_m}\end{array}\right]---\left(15\right)$

In the formula, V _i, δ _i, V _r, θ _r, I _S, φ _S, I _r, φ _r, s is the n dimension state vector of system; Y is the right characteristic vector of power flow equation Jacobian matrix; λ is that reflection system generator is meritorious, reactive power is injected and the vector of load variations; p _mFor the wind energy turbine set mechanical output is injected vector (being system's variable element vector); τ is prediction step.

4.2) proofread and correct the calculating in step

In correction format, increased 1DC equation and formula (10) simultaneous and got correction equation

$\left[\begin{array}{cccccccccccc}{f}_{{1V}_I}& f_{{1\mathrm{\δ}}_i}& f_{{1V}_r}& f_{{1\mathrm{\θ}}_r}& f_{{1I}_S}& f_{{1\mathrm{\φ}}_S}& f_{{1I}_r}& f_{{1\mathrm{\φ}}_r}& f_{1s}& 0& f_{1\mathrm{\λ}}& f_{{1P}_m}\\ f_{{2V}_i}& f_{{2\mathrm{\δ}}_i}& f_{{2V}_r}& f_{{2\mathrm{\θ}}_r}& f_{{2I}_S}& f_{{2\mathrm{\φ}}_S}& f_{{2I}_r}& f_{{2\mathrm{\φ}}_r}& f_{2s}& 0& f_{2\mathrm{\λ}}& f_{{2P}_m}\\ f_{{3V}_i}& f_{{3\mathrm{\δ}}_i}& f_{{3V}_r}& f_{{3\mathrm{\θ}}_r}& f_{{3I}_S}& f_{{3\mathrm{\φ}}_S}& f_{{3I}_r}& f_{{3\mathrm{\φ}}_r}& f_{3s}& 0& f_{3\mathrm{\λ}}& f_{{3P}_m}\\ f_{{4V}_i}& f_{{4\mathrm{\δ}}_i}& f_{{4V}_r}& f_{{4\mathrm{\θ}}_r}& f_{{4I}_S}& f_{{4\mathrm{\φ}}_S}& f_{{4I}_r}& f_{{4\mathrm{\φ}}_r}& f_{4s}& 0& f_{4\mathrm{\λ}}& f_{{4P}_m}\\ f_{{5V}_i}& f_{{5\mathrm{\δ}}_i}& f_{{5V}_r}& f_{{5\mathrm{\θ}}_r}& f_{{5I}_S}& f_{{5\mathrm{\φ}}_S}& f_{{5I}_r}& f_{{5\mathrm{\φ}}_r}& f_{5s}& 0& f_{5\mathrm{\λ}}& f_{{5P}_m}\\ f_{{6V}_i}& f_{{6\mathrm{\δ}}_i}& f_{{6V}_r}& f_{6{\mathrm{\θ}}_r}& f_{{6I}_S}& f_{{6\mathrm{\φ}}_S}& f_{{6I}_r}& f_{{6\mathrm{\φ}}_r}& f_{6s}& 0& f_{6\mathrm{\λ}}& f_{{6P}_m}\\ f_{{7V}_i}& f_{{7\mathrm{\δ}}_i}& f_{{7V}_r}& f_{{7\mathrm{\θ}}_r}& f_{{7I}_S}& f_{{7\mathrm{\φ}}_S}& f_{{7I}_r}& f_{{7\mathrm{\φ}}_r}& f_{7s}& 0& f_{7\mathrm{\λ}}& f_{{7P}_m}\\ f_{{8V}_i}& f_{{8\mathrm{\δ}}_i}& f_{{8V}_r}& f_{{8\mathrm{\θ}}_r}& f_{{8I}_S}& f_{{8\mathrm{\φ}}_S}& f_{{8I}_r}& f_{{8\mathrm{\φ}}_r}& f_{8s}& 0& f_{8\mathrm{\λ}}& f_{{8P}_m}\\ f_{{9V}_i}& f_{{9\mathrm{\δ}}_i}& f_{{9V}_r}& f_{{9\mathrm{\θ}}_r}& f_{{9I}_S}& f_{{9\mathrm{\φ}}_S}& f_{{9I}_r}& f_{{9\mathrm{\φ}}_r}& f_{9s}& 0& f_{9\mathrm{\λ}}& f_{{9P}_m}\\ \frac{\∂g}{{\∂V}_i}& \frac{\∂g}{{\∂\mathrm{\δ}}_i}& \frac{\∂g}{{\∂V}_r}& \frac{\∂g}{{\∂\mathrm{\θ}}_r}& \frac{\∂g}{{\∂I}_S}& \frac{\∂g}{{\∂\mathrm{\φ}}_S}& \frac{\∂g}{{\∂I}_r}& \frac{\∂g}{{\∂\mathrm{\φ}}_r}& \frac{\∂g}{\∂s}& \frac{\∂g}{\∂y}& \frac{\∂g}{\∂\mathrm{\λ}}& \frac{\∂g}{{\∂P}_m}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {2y}^t& 0& 0\\ {\mathrm{\ΔV}}_i& {\mathrm{\Δ\δ}}_i& {\mathrm{\ΔV}}_r& {\mathrm{\Δ\θ}}_r& {\mathrm{\ΔI}}_S& {\mathrm{\Δ\φ}}_S& {\mathrm{\ΔI}}_r& {\mathrm{\Δ\φ}}_r& \mathrm{\Δs}& {\mathrm{\Δy}}^t& \mathrm{\Δ\λ}& {\mathrm{\ΔP}}_m\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δ}}_{V_i}\\ {\mathrm{\Δ}}_{{\mathrm{\δ}}_i}\\ {\mathrm{\Δ}}_{V_r}\\ {\mathrm{\Δ}}_{{\mathrm{\θ}}_r}\\ {\mathrm{\Δ}}_{I_S}\\ {\mathrm{\Δ}}_{{\mathrm{\φ}}_S}\\ {\mathrm{\Δ}}_{I_r}\\ {\mathrm{\Δ}}_{{\mathrm{\φ}}_r}\\ {\mathrm{\Δ}}_s\\ {\mathrm{\Δ}}_y\\ {\mathrm{\Δ}}_{\mathrm{\λ}}\\ {\mathrm{\Δ}}_{P_m}\end{array}\right]=-\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6\\ f_7\\ f_8\\ f_9\\ g\\ y^{t}y-4\\ 0\end{array}\right]---\left(16\right)$

Revised variable is

$\left[\begin{array}{c}{V}_i^{(n+1)}\\ {\mathrm{\δ}}_i^{(n+1)}\\ V_r^{(n+1)}\\ {\mathrm{\θ}}_r^{(n+1)}\\ I_S^{(n+1)}\\ {\mathrm{\φ}}_S^{(n+1)}\\ I_r^{(n+1)}\\ {\mathrm{\φ}}_r^{(n+1)}\\ s^{(n+1)}\\ y^{(n+1)}\\ {\mathrm{\λ}}^{(n+1)}\\ P_m^{(n+1)}\end{array}\right]=\left[\begin{array}{c}{V}_i^{\left(n\right)}\\ {\mathrm{\δ}}_i^{\left(n\right)}\\ V_r^{\left(n\right)}\\ {\mathrm{\θ}}_r^{\left(n\right)}\\ I_S^{\left(n\right)}\\ {\mathrm{\φ}}_s^{\left(n\right)}\\ I_r^{\left(n\right)}\\ {\mathrm{\φ}}_r^n\\ s^{\left(n\right)}\\ y^{\left(n\right)}\\ {\mathrm{\λ}}^{\left(n\right)}\\ P_m^{\left(n\right)}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\ΔV}}_i^{\left(n\right)}\\ \mathrm{\Δ}{\mathrm{\δ}}_i^{\left(n\right)}\\ \mathrm{\Δ}{V}_r^{\left(n\right)}\\ \mathrm{\Δ}{\mathrm{\θ}}_r^{\left(n\right)}\\ \mathrm{\Δ}{I}_S^{\left(n\right)}\\ \mathrm{\Δ}{\mathrm{\φ}}_S^{\left(n\right)}\\ \mathrm{\Δ}{I}_r^{\left(n\right)}\\ \mathrm{\Δ}{\mathrm{\φ}}_r^{\left(n\right)}\\ \mathrm{\Δ}{s}^{\left(n\right)}\\ \mathrm{\Δ}{y}^{\left(n\right)}\\ \mathrm{\Δ}{\mathrm{\λ}}^{\left(n\right)}\\ \mathrm{\Δ}{P}_m^{\left(n\right)}\end{array}\right]---\left(17\right)$

The definition of function g and correlation computations are as follows in above-mentioned prediction-bearing calibration:

$g=f_{V_i}{|}_{*}y+f_{{\mathrm{\δ}}_i}{|}_{*}y+f_{V_r}{|}_{*}y+f_{{\mathrm{\θ}}_r}{|}_{*}y$

（18）

$+f_{I_S}{|}_{*}y+f_{{\mathrm{\φ}}_S}{|}_{*}y+f_{I_r}{|}_{*}y+f_{{\mathrm{\φ}}_r}{|}_{*}y+f_s{|}_{*}y=0$

Because it is more that above-mentioned Jacobian matrix contains neutral element, institute thinks sparse matrix.

Be under the situation of initial point at the borderline known SNB point of voltage stability domain, utilize this prediction-correction method can obtain a SNB locus of points that contains DFIG power system stability territory.

5) calculate the section of containing DFIG power system stability territory

At SNB point place formula (10) linearisation is obtained

$f_{V_i}{|}_{*}\mathrm{\Δ}{V}_i+f_{{\mathrm{\δ}}_i}{|}_{*}{\mathrm{\Δ\δ}}_i+f_{V_r}{|}_{*}\mathrm{\Δ}{V}_r+f_{{\mathrm{\θ}}_r}{|}_{*}\mathrm{\Δ}{\mathrm{\θ}}_r+f_{I_S}{|}_{*}\mathrm{\Δ}{I}_S+f_{{\mathrm{\φ}}_S}{|}_{*}\mathrm{\Δ}{\mathrm{\φ}}_S$

（19）

$+f_{I_r}{|}_{*}\mathrm{\Δ}{I}_r+f_{{\mathrm{\φ}}_r}{|}_{*}\mathrm{\Δ}{\mathrm{\φ}}_r+f_s{|}_{*}\mathrm{\Δs}+f_{\mathrm{\λ}}{|}_{*}\mathrm{\Δ\λ}+f_{P_m}{|}_{*}\mathrm{\Δ}{P}_m=0$

In the formula, Be the Jacobian matrix of system at corresponding SNB point place; f _λ| _*With

Be respectively formula (10) at this some place to λ and P _mJacobian matrix.

To the left eigenvector ω of this SNB point place Jacobian matrix zero eigenvalue correspondence of following formula (19) premultiplication,

$\mathrm{\ω}{f}_{V_i}{|}_{*}\mathrm{\Δ}{V}_i+\mathrm{\ω}{f}_{{\mathrm{\δ}}_i}{|}_{*}\mathrm{\Δ}{\mathrm{\δ}}_i+\mathrm{\ω}{f}_{V_r}{|}_{*}\mathrm{\Δ}{V}_r+\mathrm{\ω}{f}_{{\mathrm{\θ}}_r}{|}_{*}\mathrm{\Δ}{\mathrm{\θ}}_r+\mathrm{\ω}{f}_{I_s}{|}_{*}\mathrm{\Δ}{I}_s+\mathrm{\ω}{f}_{{\mathrm{\φ}}_s}{|}_{*}\mathrm{\Δ}{\mathrm{\φ}}_s$

（20）

$+{\mathrm{\ωf}}_{I_r}{|}_{*}{\mathrm{\ΔI}}_r+{\mathrm{\ωf}}_{{\mathrm{\φ}}_r}{|}_{*}{\mathrm{\Δ\φ}}_r+{\mathrm{\ωf}}_s{|}_{*}\mathrm{\Δs}+{\mathrm{\ωf}}_{\mathrm{\λ}}{|}_{*}\mathrm{\Δ\λ}+{\mathrm{\ωf}}_{P_m}{|}_{*}{\mathrm{\ΔP}}_m=0$

In the formula, " | _*" represent that the value of coefficient of correspondence is from the SNB point.

If L is load margin, i.e. λ | _*=λ ₀+ L.Again owing to but SNB place in the formula with the following formula abbreviation is

${\mathrm{\ωf}}_{\mathrm{\λ}}{|}_{*}\mathrm{\ΔL}+{\mathrm{\ωf}}_{P_m}{|}_{*}{\mathrm{\ΔP}}_m=0---\left(21\right)$

Its Be zero, formula (21) can be considered L-P _mIf the local approximate condition of voltage stability domain in the space is a known SNB point z ₁, in the corresponding number substitution formula (21) with this point, can obtain coefficient vector

This coefficient vector is exactly L-P _mVoltage stability domain border in the space is at a z ₁The normal vector at place.If z ₁Coordinate in this space is (L ₁, P _M1), then can obtain the voltage stability domain border at z according to a French equation ₁The tangent line expression formula at place is

${\mathrm{\ωf}}_{\mathrm{\λ}}{|}_{*}^1(L-L_1)+{\mathrm{\ωf}}_{P_m}{|}_{*}^1(P_m-P_{m1})=0---\left(22\right)$

The L-P that formula (22) is provided _mSection, voltage stability domain border equation is mapped to generator, load and DFIG wind energy turbine set and injects the space entirely in the space, considers following relation:

$(1+k_{i}L){P_i}^0=P_i{|}_{*}---\left(23\right)$

K in the formula _iBe each node generator of system and the meritorious growing direction that injects of load.With formula (23) and formula (22) simultaneous, can obtain injecting static voltage stability security domain border, space entirely at z in generator, load and wind energy turbine set ₁The section equation of point

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ωf}}_{\mathrm{\λ}}{|}_{*}^1\frac{(P_i{|}_{*}-P_i{|}_{*}^1)}{k_i{P_i}^0}+{\mathrm{\ωf}}_{P_m}{|}_{*}^1(P_m{|}_{*}-P_m{|}_{*}^1)=0---\left(24\right)$

Formula (24) is exactly the power system steady state voltage stability security domain local tangential plane analytic expression that contains the DFIG wind energy turbine set of finding the solution.

Above execution mode only is used for explanation the present invention; and be not limitation of the present invention; the those of ordinary skill in relevant technologies field; under the situation that does not break away from the spirit and scope of the present invention; can also make a variety of changes and modification, so all technical schemes that are equal to also belong to protection category of the present invention.

Claims (8)

1. an acquisition methods that contains the voltage stability domain of electric power system section of DFIG wind energy turbine set is characterized in that comprising the steps:

1) the DFIG wind energy turbine set is carried out mathematical modeling;

2) contain the approximation that DFIG power system steady state voltage stability border SNB is ordered with the calculating of continuous tide method;

3) the SNB point approximation of calculating with continuous tide is that initial value calculates SNB point exact value with exact algorithm;

4) be initial point with the SNB point exact value that calculates, calculate the boundary locus that contains DFIG power system stability territory;

5) calculate the section of containing DFIG power system stability territory.

2. the acquisition methods in the voltage stability domain of electric power system section of the wind energy turbine set that contains DFIG according to claim 1, it is characterized in that, the Mathematical Modeling of DFIG wind energy turbine set need retrain with the equation that contains the wind energy turbine set state parameter in the described step 1), sets up the power flow equation that contains this parameter.

3. the acquisition methods in the voltage stability domain of electric power system section of the wind energy turbine set that contains DFIG according to claim 1, it is characterized in that, described step 2) in, the state parameter tide model of the DFIG wind energy turbine set of use system, the Jacobian matrix that prediction goes on foot and correction goes on foot to the continuous tide method is made amendment, the state parameter of DFIG wind energy turbine set is added in this Jacobian matrix corresponding element, find the solution power flow equation.

4. the acquisition methods in the voltage stability domain of electric power system section of the wind energy turbine set that contains DFIG according to claim 1, it is characterized in that: the voltage stability domain absorbing boundary equation that is made of the SNB point in the described step 4) need call the power flow equation of having set up.

5. the acquisition methods in the voltage stability domain of electric power system section of the wind energy turbine set that contains DFIG according to claim 4, it is characterized in that: the voltage stability domain bounding surface equation that the SNB point constitutes adds the parametrization equation, constitutes new equation group.

6. the acquisition methods that contains the voltage stability domain of electric power system section of DFIG wind energy turbine set according to claim 5, it is characterized in that, in the equation group of the SNB locus of points that calculates DFIG power system stability territory, the method of applied forecasting-correction is followed the trail of the static voltage stability border again, constitutes new prediction step and the Jacobian matrix of proofreading and correct in the step.

7. the acquisition methods that contains the voltage stability domain of electric power system section of DFIG wind energy turbine set according to claim 1, it is characterized in that, in the described step 5), according to a French equation, try to achieve known SNB point on the voltage stability domain border in the section expression formula at this some place.

8. the acquisition methods that contains the voltage stability domain of electric power system section of DFIG wind energy turbine set according to claim 7, it is characterized in that: determine to inject static voltage stability security domain border, space entirely at the section of SNB known point equation in generator, load and wind energy turbine set, this equation is the power system steady state voltage stability security domain local tangential plane analytic expression that contains the DFIG wind energy turbine set.