CN103248049B - Containing the acquisition methods in DFIG wind energy turbine set voltage stability domain of electric power system section - Google Patents

Abstract

The invention discloses a kind of acquisition methods calculated containing DFIG wind energy turbine set voltage stability domain of electric power system section, relate to field of power, comprise after mathematical modeling is carried out to DFIG wind energy turbine set, the state parameter of wind energy turbine set is retrained; The continuous power flow that have employed prediction-correction determines a static voltage stability boundary locus, when using Continuation Method, the Jacobian matrix of prediction step and correction step is all revised, the constraint function formula of DFIG wind energy turbine set is added Jacobian matrix; Finally by the replacement of the linear parameter of this two-dimensional parameter and the calculating of section point normal equation, construct the expression formula that this contains section, DFIG power system security territory, thus make DFIG wind energy turbine set entirely inject rear system voltage stabilizes analysis intuitively, reliably, provide Data support for analyzing electrical network static electric voltage stability with security domain.

Description

Method for acquiring voltage stability domain tangent plane of electric power system of wind power plant containing DFIG

Technical Field

The invention belongs to the field of analysis and planning of power systems, and relates to a method for acquiring a static voltage stability domain tangent plane of a DFIG wind power plant and a power system.

Background

Wind power generation is an important form of wind energy utilization, the wind energy is renewable, pollution-free, large in energy and wide in prospect, and the rapid development of clean energy is a strategic choice of countries in the world. Wind has randomness and uncertainty, which makes it difficult for the power grid to absorb wind power on a large scale. The method for applying the security domain is applied to the power system comprising the wind power plant, and the security margin and the optimal control information can be provided by tracking the boundary of the system operating point and the security domain, so that the online real-time security monitoring, defense and control of the power system are more scientific and effective. The document "research report on security domain method of power system" summarizes the classification and hierarchical relationship of security domains and the related application of the security domain method. The document "calibration step-station operating conditioning for double-fed induction generators windtbins" proposes the establishment of a model of a double-fed induction generator (DFIG). At present, a wind power plant based on DFIG starts grid-connected power generation, but the research on the voltage stability of a power system after the DFIG wind power plant is injected based on a security domain tangent plane method is very little, and a power grid has high requirements on the voltage stability of the power system after the DFIG wind power plant is injected, so that the research on the static voltage stability domain tangent plane of the power system of the DFIG wind power plant has practical significance on early warning of the safety of the power grid.

Disclosure of Invention

Technical problem to be solved

The invention aims to solve the problem of the research on the voltage stability of a power system after the injection of a large-scale DFIG wind turbine generator.

(II) technical scheme

In reality, the types of generators used in most wind power places are consistent, and if the wind power plant uses a DFIG (doubly Fed Induction Generator) wind turbine, the wind power plant uses the same type of wind turbine and cannot be mixed with other types of power generation devices. The invention provides a method for solving a tangent plane of a static voltage stability safety domain of a power system containing a DFIG (doubly Fed Induction Generator), which comprises the steps of firstly carrying out mathematical modeling on the DFIG and constraining a model by utilizing seven DFIG state variables; then, tracking a two-dimensional boundary of static voltage stability of the power system by using a continuous power flow prediction and correction method to obtain a static voltage stability boundary track, and laying a foundation for determining a static voltage stability safety domain tangent plane of the power system containing the DFIG; and finally, constructing an expression of the security domain tangent plane of the DFIG-containing power system by replacing the parameters which are in linear relation with the two-dimensional parameters and calculating a tangent plane point normal equation.

(III) advantageous effects

The invention solves the problem of the research on the voltage stability of the power system after the large-scale DFIG wind turbine generator is injected.

Drawings

Fig. 1 is a schematic diagram of a calculation flow of a stable domain tangent plane acquisition method of the present invention.

Detailed Description

The following describes the embodiments of the present invention with reference to the drawings and examples. The following examples are intended to illustrate the invention, but are not intended to limit the scope of the invention.

The calculation flow shown in fig. 1.

The invention comprises the following steps:

1) determining mathematical model of electric power system of DFIG wind power plant

1.1) firstly, all the doubly-fed wind turbines are considered to be of the same type and the wind speed around the wind power plant is uniformly distributed, so that the wind power plant can be represented by a doubly-fed wind turbine model. Wherein, the mechanical energy output expression of the double-fed wind motor is

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

Wherein ρ is the air density; a is the swept area of the wind turbine; u is the wind speed; c_pThe wind energy utilization coefficient of the wind turbine shows the useful wind energy ratio obtained by the wind turbine from the wind and the tip speed ratio T_tsrIt is related.

1.2) following equation of power flow containing parameters of the electric power system of the wind power plant containing DFIG:

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

$f_2=(1+\mathrm{\λ})Q_i+Q_{im}-V_i\underset{j\∈i}{\Σ}{V}_j(G_{ij}{\mathrm{sin\δ}}_{ij}+B_{ij}{\mathrm{cos\δ}}_{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_{GB,Loss}-k{(1-s)}^3=0$

in the formula, V_iIs the magnitude of the voltage at the node,_iis the phase angle of the node voltage, V_sIs the magnitude of the stator voltage, theta_sIs the phase angle of the stator voltage, V_rIs the magnitude of the rotor voltage, theta_rIs the phase angle of the rotor voltage, I_rIs the magnitude of the rotor current phi_rIs the phase angle of the rotor current, I_sIs the magnitude of the stator current phi_sIs the stator current phase angle and s is the slip.

Will contain V_r,θ_r,I_S,φ_S,I_r,φ_rS seven equations of seven unknowns, i.e. f in equation (2)₃-f₉And after the power flow is added into a parameter-containing power flow equation of the system, the Jacobian matrix of the system is modified, and then prediction-correction calculation is carried out to obtain a static critical stable point of the DFIG-containing power system.

2) Method for calculating static voltage stable boundary initial point of DFIG-containing power system by using continuous power flow method

After the doubly-fed wind turbine is introduced, a parameter-containing power flow model of the system can be recorded as follows:

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,φ_rS is the state vector of the system; lambda is a vector reflecting active power, reactive power injection and load change of a system generator; p is a radical of_mA vector (i.e., a system variable parameter vector) is injected for wind farm mechanical power.

Firstly fix p_mConstant, lambda is free variation parameter, and continuous power flow methodA P-U curve of the active power and the voltage is plotted, and a saddle node bifurcation point (SNB point) is obtained, which is the initial boundary point. The following is a continuous power flow calculation method.

2.1) calculation of the prediction step

The prediction direction of the variable is first determined. Taking the tangential direction as an example, the equation (3) is differentiated to obtain

$\begin{array}{c}{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}\\ +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\end{array}---\left(4\right)$

If the current solution is at the trivial position of the solution curve, i.e. the jacobian matrix is not singular, then lambda is taken as the parameterized variable, the direction of change of lambda before the critical point is +1 and the direction of change of lambda after the critical point is-1. The following formula can be obtained

$\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)$

From which the prediction direction can be determined.

If the current solution is located at a position of the solution curve close to the critical point, that is, the jacobian matrix is close to singular, and the coefficient matrix of equation (5) is also close to singular, then the state variable (for example, the voltage of the node k with the largest change rate) should be selected as the parametric variation, and λ should be taken as the normal variation, and the slope of the state variable is taken as the tangential direction, at which time the prediction direction can be determined by the following equation:

$\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,is a row vector with the k-th element being +1 and the remaining elements being 0. Note that the coefficient matrix of equation (6) is different from that of equation (5). Since the state variables are selected as parameterized variables,even if the jacobian matrix is close to singular (assuming its rank is n-1), the coefficient matrix of (6) can be proven to be non-singular.

From the prediction direction determined by equation (5) or equation (6), the prediction point can be calculated as follows:

$\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) calculation of the correction step

In the correction step, if the predicted direction is obtained from equation (5), λ should be fixed first, and a vertical correction method is employed, so thatSolving the power flow equation of the formula (1) for the initial value. Taking the newton-raphson method as an example, the iteration format 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 the power flow calculation converges, a point on the P-U curve can be obtained, and then a new prediction step calculation is started. If the power flow diverges, there are two corresponding measures: one is to reduce the step size sigma, predict a new point by using the formula (7), and iterate by using the formula (8) again by using a vertical correction method; another method is to select a state variable (for example, the voltage of the node k with the largest change rate) as a parameterized variable, use λ as a common variable, solve the current equation by using a horizontal correction method, and have the following iteration format:

$\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)$

in the correction step, if the predicted direction itself is obtained from equation (6), the tidal flow equation should also be solved as in equation (9).

3) Calculating the SNB point accurate value

The constraint equation for accurately solving the SNB initial point can be described as

The boundary surface of the voltage stability region formed by the SNB points can be described as

$\left\{\begin{array}{c}f\left(V_i,{\mathrm{\δ}}_i,V_r,{\mathrm{\θ}}_r,I_s,{\mathrm{\φ}}_s,I_r,{\mathrm{\φ}}_r,s,\mathrm{\λ},P_m\right)=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_id_idV_rdθ_rdI_sdφ_sdI_rdφ_rds]

In the formula, V_i,_i,V_r,θ_r,I_S,φ_S,I_r,φ_rS is the n-dimensional state vector of the system; lambda is a vector reflecting active power, reactive power injection and load change of a system generator; p_mInjecting a vector (namely a system variable parameter vector) into the mechanical power of the wind power plant, wherein P is firstly injected_mFixing the value; and y is a right eigenvector of the jacobian matrix of the power flow equation. And taking the solution obtained by the continuous power flow as an initial value, and solving equation variables by using a Newton-Raphson method. The initial point is accurately derived.

4) Calculating SNB point track containing DFIG power system stable region

The algorithm of the stable domain SNB point trajectory also applies a prediction-correction method, and equation (10) can be abbreviated as

φ(z)＝0(11)

At this time

z^t＝[V_{i i}V_rθ_rI_Sφ_SI_rφ_rsλP_m]∈R²ⁿ⁺²,

$\mathrm{\φ}=\left[\begin{array}{c}f\left(V_i,{\mathrm{\δ}}_i,V_r,{\mathrm{\θ}}_r,I_s,{\mathrm{\φ}}_s,I_r,{\mathrm{\φ}}_r,s,\mathrm{\λ},P_m\right)\\ 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}$

The dimension of the Jacobian matrix determined by the above formula is (2n +1) × (2n +2), and a parameterized equation is added therefor

(z-z₁)^tv＝τ(12)

Wherein v is z₁Tangent vector of point, z₀For the initial point, τ is the step size. v is determined by the following formula:

4.1) calculation of the prediction step

First assume that a point z has been calculated from step 2)₀On the power system stability domain SNB point trajectory. From this point on, the prediction direction is determined by:

$\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{\λ}& {\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]\mathrm{---}\left(\mathrm{16}\right)$

from the prediction direction determined in (14), the prediction points can be calculated as follows:

$\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,φ_rS is the n-dimensional state vector of the system; y is a right eigenvector of a Jacobian matrix of the power flow equation; lambda is a vector reflecting active power, reactive power injection and load change of a system generator; p is a radical of_mInjecting a vector (namely a system variable parameter vector) for the mechanical power of the wind power plant; τ is the prediction step size.

4.2) calculation of the correction step

In the correction format, a one-dimensional correction equation is added to be combined with the correction equation expressed by the formula (10)

$\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{\λ}& {\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-1\\ 0\end{array}\right]\mathrm{---}\left(\mathrm{16}\right)$

The corrected variable is

$\left[\begin{array}{c}{V_i}^{\left(n+1\right)}\\ {\mathrm{\δ}}_i^{\left(n+1\right)}\\ {V_r}^{\left(n+1\right)}\\ {\mathrm{\θ}}_r^{\left(n+1\right)}\\ I_S^{\left(n+1\right)}\\ {\mathrm{\φ}}_S^{\left(n+1\right)}\\ I_r^{\left(n+1\right)}\\ {\mathrm{\φ}}_r^{\left(n+1\right)}\\ s^{\left(n+1\right)}\\ y^{\left(n+1\right)}\\ {\mathrm{\λ}}^{\left(n+1\right)}\\ P_m^{\left(n+1\right)}\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^{\left(n\right)}\\ 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{\Δ\δ}}_i^{\left(n\right)}\\ {{\mathrm{\ΔV}}_r}^{\left(n\right)}\\ {\mathrm{\Δ\θ}}_r^{\left(n\right)}\\ {\mathrm{\ΔI}}_S^{\left(n\right)}\\ {\mathrm{\Δ\φ}}_S^{\left(n\right)}\\ {\mathrm{\ΔI}}_r^{\left(n\right)}\\ {\mathrm{\Δ\φ}}_r^{\left(n\right)}\\ {\mathrm{\Δs}}^{\left(n\right)}\\ {\mathrm{\Δy}}^{\left(n\right)}\\ {\mathrm{\Δ\λ}}^{\left(n\right)}\\ {\mathrm{\ΔP}}_m^{\left(n\right)}\end{array}\right]---\left(17\right)$

The definition and the related calculation of the function g in the above prediction-correction method are as follows:

$\begin{array}{c}g=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\end{array}---\left(18\right)$

the jacobian matrix is a sparse matrix because it contains more zero elements.

Under the condition that an SNB point on the boundary of the voltage stability domain is known as an initial point, an SNB point track containing the stability domain of the DFIG power system can be obtained by utilizing the prediction-correction method.

5) Calculating tangent plane containing DFIG power system stable region

Linearizing equation (10) at the SNB point to obtain

$\begin{array}{c}{f}_{V_i}{|}_{*}{\mathrm{\ΔV}}_i+f_{{\mathrm{\δ}}_i}{|}_{*}{\mathrm{\Δ\δ}}_i+f_{V_r}{|}_{*}{\mathrm{\ΔV}}_r+f_{{\mathrm{\θ}}_r}{|}_{*}{\mathrm{\Δ\θ}}_r+f_{I_S}{|}_{*}{\mathrm{\ΔI}}_S+f_{{\mathrm{\φ}}_S}{|}_{*}{\mathrm{\Δ\φ}}_S\\ +f_{I_r}{|}_{*}{\mathrm{\ΔI}}_r+f_{{\mathrm{\φ}}_r}{|}_{*}{\mathrm{\Δ\φ}}_r+f_s{|}_{*}\mathrm{\Δ}s+f_{\mathrm{\λ}}{|}_{*}\mathrm{\Δ}\mathrm{\λ}+f_{P_m}{|}_{*}{\mathrm{\ΔP}}_m=0\end{array}---\left(19\right)$

In the formula,a Jacobian matrix at the corresponding SNB point of the system; f. of_λ|_*Andat this point, the pairs of λ and P are respectively of the formula (10)_mThe derivative matrix of (a).

The left eigenvector omega corresponding to the zero eigenvalue of the Jacobian matrix at the SNB point is multiplied by the above formula (19) to obtain

$\begin{array}{c}{\mathrm{\ωf}}_{V_i}{|}_{*}{\mathrm{\ΔV}}_i+{\mathrm{\ωf}}_{{\mathrm{\δ}}_i}{|}_{*}{\mathrm{\Δ\δ}}_i+{\mathrm{\ωf}}_{V_r}{|}_{*}{\mathrm{\ΔV}}_r+{\mathrm{\ωf}}_{{\mathrm{\θ}}_r}{|}_{*}{\mathrm{\Δ\θ}}_r+{\mathrm{\ωf}}_{I_S}{|}_{*}{\mathrm{\ΔI}}_S+{\mathrm{\ωf}}_{{\mathrm{\φ}}_S}{|}_{*}{\mathrm{\Δ\φ}}_S\\ +{\mathrm{\ωf}}_{I_r}{|}_{*}{\mathrm{\ΔI}}_r+{\mathrm{\ωf}}_{{\mathrm{\φ}}_r}{|}_{*}{\mathrm{\Δ\φ}}_r+{\mathrm{\ωf}}_s{|}_{*}\mathrm{\Δ}s+{\mathrm{\ωf}}_{\mathrm{\λ}}{|}_{*}\mathrm{\Δ}\mathrm{\λ}+{\mathrm{\ωf}}_{P_m}{|}_{*}{\mathrm{\ΔP}}_m=0\end{array}---\left(20\right)$

In the formula, "[ mu ] n_*"indicates that the value of the corresponding coefficient is from the SNB point.

Let L be the load margin, i.e. lambda-_*＝λ₀+ L. And the SNB position in the formula can be simplified into

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

It is composed ofAre all zero, formula (21) can be regarded as L-P_mLocal approximation condition of voltage stability domain in space if one SNB point z is known₁The corresponding number of the point is substituted into the formula (21) to obtain a coefficient vectorThe coefficient vector is L-P_mVoltage stability domain boundary in space at point z₁The normal vector of (c). Let z₁The coordinate in this space is (L)₁,P_m1) Then, the boundary of the voltage stability region at z can be obtained according to the equation of the point-normal equation₁The tangent line expression of (A) is

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

L-P given by formula (22)_mAnd mapping a voltage stability domain boundary tangent plane equation in the space to a generator, a load and a DFIG wind power plant full injection space, wherein the following relation is considered:

(1+k_iL)P_i ⁰＝P_i|_*(23)

in the formula k_iThe active energy injection direction is the increasing direction of the generators and the loads of all nodes of the system. By combining the formula (23) and the formula (22), the boundary of the stable safety domain of the static voltage in the full injection space of the generator, the load and the wind farm is obtained at z₁Tangent plane equation of point

$\underset{i=1}{\overset{n}{\Σ}}{\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)$

And the equation (24) is a local tangent plane analytical equation of the static voltage stable security domain of the electric power system of the DFIG-containing wind power plant.

The above embodiments are merely illustrative, and not restrictive, and those skilled in the relevant art can make various changes and modifications without departing from the spirit and scope of the invention, and therefore all equivalent technical solutions also fall within the scope of the invention.

Claims (2)

1. A method for obtaining a voltage stability domain tangent plane of a power system of a wind power plant containing DFIG is characterized by comprising the following steps:

1) performing mathematical modeling on a DFIG wind power plant to obtain a static critical stability point of a DFIG-containing power system;

2) calculating an approximate value of a saddle node bifurcation point of a static voltage stability boundary of a DFIG-containing power system by using a continuous power flow method;

3) calculating a saddle node bifurcation point accurate value by using an accurate algorithm by taking the saddle node bifurcation point approximate value calculated by the continuous power flow as an initial value;

4) calculating a boundary track containing a stable region of the DFIG power system by taking the calculated saddle node bifurcation point accurate value as an initial point;

5) in the full injection space of a generator, a load and a wind power plant, determining a tangent plane equation of the boundary of a static voltage stability safety domain at a known point of a saddle node bifurcation by adopting a point-normal equation mode, wherein the tangent plane equation is a local tangent plane analytical formula of the static voltage stability safety domain of a power system of the wind power plant containing the DFIG;

the mathematical model of the DFIG wind power plant in the step 1) needs to be constrained by an equation containing state parameters of the wind power plant, a power flow equation containing the parameters is established, and in the power flow equation, the modified Jacobian matrix is predicted and corrected by modifying the Jacobian matrix so as to obtain a static critical stability point containing the DFIG power system;

the boundary equation of the voltage stabilization domain formed by the saddle node bifurcation points in the step 4) needs to call the established power flow equation;

adding a parameterized equation into a voltage stability domain boundary surface equation formed by the saddle node bifurcation points to form a new equation set;

in an equation set for calculating the saddle node bifurcation point locus in the stable domain of the DFIG power system, tracking a static voltage stable boundary by using a prediction-correction method again to form a new Jacobian matrix in a prediction step and a correction step.

2. The method for obtaining the tangent plane of the voltage stabilization domain of the electric power system of the DFIG-containing wind power plant as recited in claim 1, wherein in the step 2), a state parameter power flow model of the DFIG wind power plant of the system is used to modify a Jacobian matrix of a prediction step and a correction step of a continuous power flow method, and the state parameters of the DFIG wind power plant are added into corresponding elements of the Jacobian matrix to solve a power flow equation.