A method for constructing power system stability domain considering randomness of new energy systems and power grid failures

Abstract

The present invention is based on the method for constructing the transient energy function of the power system based on Lyapunov transient stability theory, obtains the historical data of DFIG and obtains the wind speed sequence in the time domain, constructs a random transient energy function model of the power system containing wind power plants based on random processes, and realizes the intuitive embodiment of randomness by setting DFIGs with different grid-connected capacities. Then, on this basis, the stable domain and node stable nominal values ​​are constructed, and according to their distribution positions in the stable domain and the movement trend of the overall unstable equilibrium point of the system in the stable domain, the identification of different fault locations and fault types and the judgment of the overall operation of the system are realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of wind power generation, and more specifically, relates to a method for constructing a random transient energy function and a system stability domain of a new energy power system based on Lyapunov transient stability theory. Background Art

[0002] Due to the rapid development of new energy power systems, new energy systems represented by wind power injection are connected to the power grid on a large scale, and the resulting power system is more prone to transient processes than traditional power systems. Due to the grid-connected operation of wind turbines, the randomness of wind energy itself is inevitably introduced into the traditional power system. These changes are bound to change the operating mechanism of traditional power systems, and the transient stability analysis process will also present new characteristics. Compared with the determined transient energy function, this paper introduces the randomness of wind speed into it to obtain a random transient energy function, which makes it more in line with the operating characteristics of new energy power systems. Therefore, it can be applied to the transient stability analysis of various new energy systems.

[0003] Doubly fed induction generator (DFIG) is the most widely used wind turbine in the field of wind power generation. It consists of a wound asynchronous generator with a stator winding directly connected to the phase grid and a bidirectional back-to-back IGBT voltage source converter installed on the rotor winding. By analyzing the internal structure of DFIG, it can be seen that the key factors affecting the stability of the system are the randomness of the input evolution over time and the volatility of the system.

[0004] As for the randomness evolving over time, with the deepening of the study of random processes, Gaussian white noise, Wiener process and Weibull distribution are all used to describe the randomness of wind speed, but there is room for further improvement in terms of calculation time and calculation accuracy. The analysis of system stability based on the direct method has the advantage of being convenient and intuitive. With the deepening of the study of transient stability of power systems, the direct method is used in the calculation of critical clearing time (CCT) and the application of extended equal area criterion (EEAC). The direct method was then further developed into the transient energy function method. As for the volatility of the system, load changes, fault locations and fault holding times are more random. The general differential equations (ODEs) commonly used in traditional power systems are no longer applicable to the expanding power system containing DFIG. Although a random transient energy function suitable for wind turbines can be obtained from the DFIG itself, the grid connection of wind turbines requires the access of a large number of power electronic devices to control them. Characterizing the transient energy function of DFIG from a single perspective obviously cannot accurately represent its impact on the overall operating state. It can be seen that the establishment of a complete stochastic differential equation that can reflect the transient stability of the power system containing DFIG is very important for studying the transient stability analysis under the new situation when new energy is connected to the grid.

[0005] In terms of stable domain construction, the boundary of the stable domain of the power system composed of unstable equilibrium points and limit cycles is currently mainly obtained based on the theory of stable manifolds. With the improvement of computer computing power, stable manifolds have been further developed and applied to Gaussian random matrices. However, a common problem with this method is that this stable domain cannot be directly applied to high-order power systems. Therefore, it is necessary to group large power systems before solving them. The large number of simplified calculations generated in this process makes it difficult to ensure the credibility of the results. With the development of LaSalle's invariance principle, the application of Lyapunov stability theory in power systems has gradually been promoted. Based on the positive invariant set theory, the absorption domain of its nonlinear dynamic equilibrium points is analyzed, and an estimation algorithm that expands the absorption domain (region of attraction, ROA) in each iteration can be obtained. On this basis, the stable state trajectory of the system dynamics is used to construct the ROA of the system state of the power system containing new energy centered on stable equilibrium. However, due to the influence of the sampling algorithm, the above research ideas are difficult to obtain actual physical meanings in terms of the stable domain boundary and confidence. In summary, there are two main directions for the current research on the stability domain of power systems containing wind turbines. The first is a mechanism-based study, which studies and identifies the location of the system's stability domain and the movement trend of each state point based on the mechanism of the power system. The second is a method-based study, which uses the perspective of convex sets to obtain elliptical sets, budget constraint sets, interval sets, and convex hulls to represent the random characteristics of the power system. However, method-based research has limitations in practical applications and does not have a good physical explanation. Summary of the invention

[0006] The purpose of the present invention is to overcome the deficiencies of the prior art and provide a method for constructing a random transient energy function and a system stability domain of a new energy power system based on Lyapunov transient stability theory. Based on Lyapunov transient stability theory and the motion model of each subsystem in the new energy power system, random transient energy functions are constructed for the traditional power system and DFIG respectively, and the traditional transient energy function is improved by considering the existence of wind speed as a random source input. On this basis, a new method for constructing a stability domain for a new energy power system is proposed, so as to more intuitively determine the operating status of the power system and each node, as well as the location and type of fault.

[0007] In order to achieve the above-mentioned invention object, the method for constructing the random transient energy function and system stability domain of the new energy power system based on Lyapunov transient stability theory of the present invention comprises the following steps:

[0008] S1: Obtain the wind speed sequence of the historical wind farm and the historical wind turbine speed sequence of the DFIG, and then obtain the wind speed suitable for the system operation step length that develops over time according to the random wind speed model pre-set based on the basic conditions of the Weibull distribution to obtain the wind speed sequence;

[0009] S2: The following method is used to construct a random transient energy function model of the traditional power system based on the equation of motion and energy conservation, including the random transient energy function models of the synchronous generator side, the transmission network side and the load side.

[0010] Based on the relationship between the generator terminal current, voltage, voltage phase angle at its node and system power angle, and after all parameters are converted to per-unit values, the rotor motion equation can be directly used with mechanical power P m and electromagnetic power P e Instead of mechanical torque T m and electromagnetic torque T e , so the potential energy function W of the synchronous generator can be obtained SynL As shown below:

[0011]

[0012] Where V represents the generator terminal voltage, θ represents the voltage phase angle at the generator node, δ represents the power angle, ψ D , md , f , Q , mq and ψ g They represent the flux linkage of the damping winding D, d-axis, the excitation winding, the damping winding Q, q-axis and the damping winding g, respectively. Dl 、x fl 、x Ql 、x gl 、x l 、x md and x mq They represent the damping winding D, the excitation winding, the damping winding Q, g and the leakage inductance of the stator, as well as the mutual inductance between the dq axis windings, u f and R f They represent the voltage and impedance of the excitation winding respectively. When the subscript s is used in the text, it refers to the corresponding parameters when the system is in a stable state.

[0013] Assume that the number of synchronous generators in the power system is n g The overall energy function W of the synchronous generator part is obtained. Syn As shown below:

[0014]

[0015] Where M represents the inertia time constant, ωr Indicates the rotor speed.

[0016] For the dissipated energy on the transmission network side and the load side, when a fault occurs in the transmission network, the state variables such as the amplitude and phase angle of the node voltage, the active and reactive power of the line will change. From the energy function, these changes represent the conversion of the overall kinetic energy of the power grid into potential energy, and these potential energies will eventually be absorbed by the transmission network and the load. Therefore, the impact of the potential energy of the transmission network on the overall energy function of the system cannot be ignored. On high-voltage transmission lines, the conductivity is small enough to be ignored. Through flow calculation and related derivation and simplification, the energy function W of the transmission network is obtained. Line As shown below:

[0017]

[0018] Among them, V i represents the voltage on node i, θ i represents the voltage phase angle of node i, B ij represents the susceptance between nodes i and j. g The first node is a synchronous generator node, the nth g +1 node is a balance node, nth g +2 to nth g +m+1 nodes are load nodes.

[0019] On the load side, changes in the system as a whole will also affect the load's consumption of active and reactive power, which is still manifested as the conversion and absorption of potential energy from the perspective of energy function. Therefore, the potential energy absorbed in the load W Load The expression is as follows:

[0020]

[0021] Among them, Q Lml , U ml , P Lml ,θ ml They represent the reactive power consumption, voltage, active power consumption and voltage phase angle at the mlth load node respectively. Generally speaking, the active and reactive power injection and consumption of the power grid system are dynamically balanced.

[0022] S3: The following method is used to construct a random transient energy function model of the DFIG part based on the equation of motion and energy conservation, which consists of four parts: the random transient energy function represented by active output and reactive output, the random transient energy function stored in the inductance inside the DFIG, and the random transient energy function represented by kinetic energy.

[0023] For DFIG in the transient process, θ is constantly changing, but the magnetomotive force generated in different coordinate systems is consistent, which meets the basic requirements of the coordinate transformation method. Considering that the state quantity is transformed from time-varying to time-invariant, the coordinate transformation is divided into two steps as shown below.

[0024] Step 1: 3s / 2s transformation, assuming that the effective number of turns of each phase of the three-phase winding is N3, and the effective number of turns of each phase of the two-phase winding is N2, then according to the principle of equal instantaneous magnetomotive force, and considering that the power is unchanged, the conversion to matrix form is shown as follows:

[0025]

[0026] Step 2: 2s / 2r transformation, if the d and q axis currents are set to i d 、i q , and then C 2r / 2s Taking the inverse, we get the transformation matrix C from the two-phase stationary coordinate system to the two-phase rotating coordinate system 2s / 2r As shown below:

[0027]

[0028] Combining the transformation matrix obtained by 3s / 2s transformation and 2s / 2r transformation, the equivalent circuit diagram of the corresponding DFIG is obtained. Figure 1 It is the d and q axis equivalent circuit diagram of DFIG, such as Figure 1 It should be noted that for the rotor winding itself, the angle leading the stator winding is θ r .in accordance with Figure 1 , the voltage, current and flux equations of DFIG in the two-phase rotating coordinate system can be obtained, and then the differential equations of the flux in all directions are obtained by substituting them into the voltage equation and transforming them as shown in the following formula:

[0029]

[0030] Among them, ψ sd , sq , rd , rq Respectively represent the magnetic flux components of the stator and rotor windings on the d and q axes, R s and R r Respectively represent the stator and rotor winding resistance, u sd 、u sq 、u rd 、u rq Represent the voltage components of the stator and rotor windings on the d and q axes, ω1 and ω s Represent stator speed and slip speed respectively, L ls and L lr They represent the stator and rotor winding self-inductance and the mutual inductance between groups Lm The difference.

[0031] There is dissipated energy of magnetic field in the inductance of DFIG, which is also potential energy converted from kinetic energy in transient process. Incorporating this potential energy into the energy function model can make the model more applicable. According to the characteristics of energy, adding the right sides of the four equations after integration of the above formula can obtain the potential energy W converted in the DFIG inductance represented by magnetic flux and voltage. L As shown below:

[0032]

[0033] According to the actual situation, the voltage change and transmission network loss caused by stator transient, stator resistance, and stator speed change can be ignored. Substitute the generator terminal voltage, the voltage phase angle and power angle at the node into the stator dq axis current expression to obtain the active and reactive output of DFIG. Due to the power balance of the system, for the active part ΔP=P l +P k -P out , for the reactive part ΔQ=Q l +Q k -Q out , where P l and Q l Represent the active and reactive losses of the system, P k and Q k Respectively represent the active power and reactive power consumed by the load. In fact, ∫P out dδ represents the potential energy converted from kinetic energy, part of which is expressed as W L is stored in the inductor in the form of ∫P out Therefore, when the loss is ignored and the per-unit value is taken, the active output P out Output energy W corresponding to the system active power P It can be expressed as follows:

[0034]

[0035] Among them, s r is the generator slip, P s =u sd i sd +u sq i sq ,i sd and i sq Represent the current components of the stator and dq axes respectively. For reactive output Q out There is no simplified formula like the above, so it is necessary to calculate the rotor side reactive output Q rFirst, the reactive output of the rotor side is converted to the stator side. At the same time, the stator d and q axis currents are represented by vectors and substituted into the reactive output Q represented by the stator side. s After taking the per-unit value, we can get the output energy Q corresponding to the overall reactive output and system reactive power out It can be expressed as follows:

[0036]

[0037] Since the grid directly consumes reactive power, and the change of reactive power leads to voltage change, which affects the voltage stability of the system, that is, the potential energy corresponding to the reactive power output of the generator is directly related to the voltage. Therefore, this part of transient energy can be expressed by ∫Q out / VdV is used to express it.

[0038] Based on the operating mechanism of DFIG, its active power can be expressed as P m =P S +P R , where P m Represents mechanical input power, P S , P R denote the power output on the stator and rotor sides respectively. And P S =T e ω s , P m =T m ω r , where T e represents the electromagnetic torque on the rotor side, ω s Indicates the stator speed, T m represents the mechanical torque on the rotor side, ω r Indicates the rotor speed.

[0039] When DFIG is running stably, T e =T m , the power output on the rotor side can be expressed as When it enters the transient process, the motion equations of the wind turbine, drive shaft and doubly fed generator rotor side are as follows:

[0040]

[0041] Where J is the moment of inertia of the DFIG, D is the damping coefficient on the rotor side, θ r represents the rotation angle, n p represents the number of pole pairs of the doubly-fed generator, L m represents the mutual inductance between the stator and rotor windings, i sd 、i sq 、i rd 、i rqrepresents the current components of the stator and rotor windings on the d-axis and q-axis respectively. Therefore, set M = J / n p , D T =D / n P Based on the first and third formulas of the above formula, the transient process expression of the rotor motion equation can be obtained as shown below:

[0042] -D T Δω r dΔω r =MΔω r dΔω r -T m dθ r +T e dθ r

[0043] By integrating both sides of the equation above, we can obtain the transient energy function W represented by the rotor motion equation: R , as shown below:

[0044]

[0045] Among them, W ke It represents the kinetic energy of DFIG in transient process, W pe It represents the potential energy converted from kinetic energy of DFIG in transient process. Therefore, by substituting the motion equations of wind turbine, transmission shaft and rotor side of doubly fed generator into the above equation, it can be transformed into the following equation:

[0046]

[0047] Figure 2 is the relationship between DFIG rotor speed and output power. Figure 2 It can be seen that the input power and rotor speed of DFIG are linear in the AB and CD parts. Therefore, the ratios set in the grid-connected startup area and the constant speed area are expressed as K B and K N In the BC part, based on the coordinates of point B (P outB ,ω rB ) and the coordinates of point C (P outC ,ω rC ), the present invention uses a power exponential function to fit it. Therefore, the output power of DFIG can be defined as shown in the following formula:

[0048]

[0049] Assuming that the energy loss ΔP of the DFIG itself is approximately a constant, the relationship between the input power and the output power of the power system can be expressed as P m =P out+ΔP. Therefore, combining the motion equations of the wind turbine, drive shaft and DFIG rotor side and its output power, the DFIG rotor speed can be defined as follows:

[0050]

[0051] Among them, the input power of DFIG can be expressed as ρ represents the ambient gas density, R represents the radius of the windward surface of the DFIG, v represents the wind speed, C p (λ, β) represents the wind energy utilization coefficient, where λ and β represent the blade tip speed ratio and pitch angle, respectively. Therefore, when the DFIG operates in the grid-connected startup area and the constant speed area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows:

[0052]

[0053] Among them, K i To represent K B and K N When the DFIG operates in the maximum wind energy tracking area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows:

[0054]

[0055] set up is the proportional coefficient of the DFIG input part, so the above two equations can be simplified and uniformly expressed as the random transient energy function represented by kinetic energy inside the DFIG as shown below.

[0056]

[0057] In summary, through the above analysis, all transient energy functions are integrated. For the transient energy function of the entire power grid system, assuming that there are s wind turbines in the power grid, n g synchronous generators, a total of s+n g generators, of which one is a synchronous generator as a balancing generator. The overall random transient energy function can be expressed as:

[0058]

[0059] S4: The following method is used to construct a transient stability assessment model of the DFIG-based power system containing wind power plants and the stable nominal values ​​of each node.

[0060] In order to characterize the transient stability of the power system, the present invention constructs a transient stability assessment model for the power system containing wind farms based on DFIG. The model is based on the above-mentioned stochastic transient energy function (STEF) and random matrix theory, and consists of a basic region (BR), a stochastic region (SR) and a hazard region (HR). The center of the stable region is defined as a trap point (TP). In addition, the present invention proposes a new bus stable nominal value (BSNV) as an evaluation criterion. Figure 3 The flowchart of the specific implementation method of the present invention for constructing the stability domain of the power system containing DFIG based on the random transient energy function and random matrix theory is as follows. Figure 3 shown.

[0061] exist Figure 3 Based on the data obtained from the wind farm and the derived STEF, the mathematical model of STEF can be obtained as W k , arranged according to the node number as W kΓ As the input layer of the calculation stability domain, k represents the time series, Γ represents the node number, the eigenvalue matrix corresponding to the input layer is λ and tr(λ) is the trace of the eigenvalue matrix. Therefore, the distance between the HR boundary and TP can be defined as follows:

[0062]

[0063] n represents the number of nodes in the DFIG-based power system with wind farms. For the boundary of SR, a new random matrix can be obtained by substituting the parameters obtained from the stable operation of the system into STEF. Then repeat the above steps to calculate the boundary of SR, and finally set the boundary of BR to 1;

[0064] In the BSNV calculation process, the active power and reactive power of each node are used to calculate the STEF of the corresponding node, so the corresponding random matrix W can be obtained. Γ and the corresponding maximum nominal value Load Γ , which can be used to describe the fluctuation of various parameters in the node. The degree of eccentricity of BSNV can be obtained through the following mathematical transformation:

[0065]

[0066] Among them, Λ represents a unitary matrix, Load Γ is the complex domain representation of the node stable nominal value, WΓ represents the STEF of the corresponding node in the time domain, represents its conjugate transposed matrix, W λ represents its eigenvalue matrix; the maximum value in the eigenvalue matrix is ​​defined as the BSNV of the corresponding node. After normalization, its horizontal and vertical coordinates in the power system stability domain composed of BR / SR / HR are represented as BSNV x =Re[Load Γ ] and BSNV y =Im[Load Γ ];

[0067] It is assumed that when more than 60% of BSNV enters HR, the DFIG-based power system with wind farm will be judged to be in an unstable state. Therefore, the position of the unstable equilibrium point (UEP) of the entire power system in the stable domain can be defined as shown in the following formula. It should be noted that it has the same characteristics as BSNV. Since BSNV is a standardized value, the following formula can be used to reflect the overall operation of the power system with wind farm:

[0068]

[0069] Among them, R UEP Indicates the distance between UEP and TP in the power system, Load Γ Represents the BSNV of node Γ. When UEP is in BR, it means that the power system is stable. When UEP enters SR, it can be determined that the randomness of the power system caused by the input tends to be obvious, and the various parameters of the system during operation will also fluctuate significantly, especially the parameters corresponding to the wind farm grid connection point. When UEP enters HR, it can be determined that the power system has entered an unstable operation state.

[0070] The present invention is based on the method for constructing the transient energy function of the power system based on Lyapunov transient stability theory, obtains the historical data of DFIG and obtains the wind speed sequence in the time domain, constructs a random transient energy function model of the power system containing wind power plants based on random processes, and realizes the intuitive embodiment of randomness by setting DFIGs with different grid-connected capacities. Then, on this basis, the stable domain and node stable nominal values ​​are constructed, and according to their distribution positions in the stable domain and the movement trend of the overall unstable equilibrium point of the system in the stable domain, the identification of different fault locations and fault types and the judgment of the overall operation of the system are realized.

[0071] The present invention has the following beneficial effects:

[0072] (1) Different from the traditional transient stability analysis method of power system, the present invention is based on Lyapunov stability theory and uses the motion equations and energy distribution of each part of the system to construct a random transient energy function, introduces randomness into the transient stability analysis of the power system containing DFIG, and overcomes the conservatism of the traditional energy function.

[0073] (2) Traditional methods often ignore the randomness of wind speed when dealing with transient stability problems, and cannot accurately simulate the actual working scene of wind turbines. Therefore, they are conservative in describing the system operation process. The present invention analyzes the entire system based on the consideration of random wind speed information, thereby enhancing robustness.

[0074] (3) Based on the random transient energy function, the present invention applies it to construct the stability domain of the power system, and displays the stability trend of the entire system and the parameter fluctuation intensity of the node itself by characterizing the distribution position of each node in the stability domain, so as to more intuitively judge the position of the fluctuation node and the unstable node. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 is the d-axis and q-axis equivalent circuit diagram of the DFIG in the present invention;

[0076] Figure 2 is a relationship diagram between the DFIG rotor speed and the output power in the present invention;

[0077] Figure 3 It is a flowchart of a specific implementation method of constructing a stable domain of a power system containing DFIG derived from the random transient energy function and random matrix theory of the present invention;

[0078] Figure 4 It is a flowchart of a specific implementation method of a method for constructing a random transient energy function and a system stability domain of a new energy power system based on Lyapunov transient stability theory of the present invention;

[0079] Figure 5 It is a curve diagram of the distribution trend and envelope of the random transient energy function of the power system containing the wind power plant under different DFIG grid-connected capacities after the occurrence of wind speed fluctuation and grid-side three-phase fault in the present invention;

[0080] Figure 6 It is a distribution diagram of the stability domain and node stability nominal value of the power system containing wind power plants after wind speed fluctuation and grid-side three-phase fault under different DFIG grid-connected capacities obtained in the present invention based on random transient energy function and random matrix theory;

[0081] Figure 7 It is a comparative curve diagram of voltage change trend at the grid connection point after wind speed fluctuation and grid-side three-phase fault occur under different DFIG grid connection capacities in the present invention;

[0082] Figure 8 It is the motion trajectory of the stable domain of the power system containing DFIG and its dominant unstable equilibrium point obtained based on the potential energy boundary method in the present invention. DETAILED DESCRIPTION

[0083] The specific implementation of the present invention is described below in conjunction with the accompanying drawings so that those skilled in the art can better understand the present invention. It should be noted that in the following description, when the detailed description of known functions and designs may dilute the main content of the present invention, these descriptions will be omitted here.

[0084] Example

[0085] In order to better illustrate the technical solution of the present invention, the derivation process of the wind speed model involved in the present invention is first briefly described.

[0086] Because the wind in a wind farm is random and volatile, the wind speed model generally uses a four-component superposition method, where the four components are basic wind, gust, gradual wind and random wind.

[0087] For basic wind, the specific value of basic wind can be obtained by analyzing the measured wind data and further determining the relevant parameters of Weibull distribution. It should be noted that basic wind directly affects the power delivered by wind turbines to the grid. Therefore, basic wind can be defined as V B =AΓ(1+1 / K), where A represents the proportionality coefficient, Γ(1+1 / K) represents the gamma function, and K represents the scale parameter of the Weibull distribution.

[0088] For gusts, the dramatic change in wind speed will inevitably affect the output power of the wind turbine, thereby affecting voltage fluctuations. According to its physical characteristics, gusts can be defined as V G =V S (t1≤t≤t1+T g ),in t、T G , t1 represent the actual time, gust duration and gust occurrence time respectively, and maxG represents the maximum value of the gust.

[0089] For the gradually changing wind, since it has a linear property within a certain period of time, it can be defined as V R =V γ (T 1R ≤t≤T 2R ), where V γ =maxR[1-(tT 2R ) / (T 1R -T 2R )], maxR is the maximum value of the gradient wind, T 1R , T2R Respectively represent the occurrence time and end time.

[0090] For random wind, since the random wind speed changes at different heights are not the same, this paper uses Weibull distribution as the initial value of random wind, and then uses the Wiener process as the development trend of random wind. The mathematical expression of the change of random wind speed can be defined as. Therefore, the input wind speed of DFIG can be expressed as V N .

[0091] Since the overall model of the new energy power system is complex and the internal mechanical and electrical structures are different, it is very difficult to directly use the traditional transient energy function, and it cannot reflect the impact of random sources on the changes in the system energy function. Therefore, the present invention has made adaptive improvements to the transient energy function and proposed a new random transient energy function. This random transient energy function divides the power system containing a wind power plant into seven parts, namely, the random transient energy functions of the synchronous generator, transmission network and load of the traditional power system, and the random transient energy function represented by the kinetic energy and potential energy of the DFIG inside the wind farm, wherein the potential energy includes the potential energy random transient energy function represented by active output and reactive power, and the random transient energy function stored in the internal inductance. The above seven parts cover all parts of the power system containing a wind power plant, so they can be used to calculate the transient energy function of the system as a whole and each node. At the same time, by making full use of the random information of the wind speed, the impact of the randomness of the wind speed on the energy changes of the system as a whole and each part can be more intuitively reflected.

[0092] In the power system, the stability domain is a relatively intuitive and convenient model for judging the operating state of the transient process of the power system. By depicting the operating trajectory of the unstable equilibrium point of the system, the transient operation process of the system can be tracked and determined. Based on the above-mentioned random transient energy function, the present invention adopts the basic idea of ​​the single-loop theorem on the basis of the stability domain theory to construct a stability domain for the power system containing a wind power plant. At the same time, in order to adapt to the possible existence of different locations and different types of random faults, the location of each node is marked in the stability domain to achieve the following goals:

[0093] (1) Comprehensively analyze and judge the overall operation of the system to determine whether it is in the stable zone, random zone or risk zone;

[0094] (2) Each node is individually labeled and its parameter fluctuation is determined by the degree of eccentricity;

[0095] (3) Identify and judge different possible fault types and faults occurring in different locations.

[0096] In constructing the random transient energy function, it is necessary to introduce a wind speed model for the wind turbine model and actual operating conditions. In order to solve the input-output volatility problem caused by the random fluctuation of wind speed, the present invention introduces the kinetic energy of the random source through the DFIG into the random transient energy function on the basis of the motion equations of each part. This can avoid the conservatism caused by the relaxation of the overall constraints and protect the relevant mechanical components in the wind turbine.

[0097] With the improvement of computer performance and the need to analyze practical problems, the classical state space model and the network structure retention model have been proposed one after another. In the present invention, the network structure retention model will be used as the basis for analysis. Therefore, for a synchronous generator, it can be divided into the potential energy converted in the internal inductance of the generator and the rotor kinetic energy and the power output energy converted therefrom. For a typical 6th-order synchronous generator equipped with a PSS controller, the d-axis contains an excitation winding f and a damping winding D, and the q-axis contains damping windings Q and g. On the premise of ignoring the voltage changes caused by stator transients, resistance and rotor speed changes, the voltage equation and flux equation in the dq coordinate system are substituted into the current equation, and the leakage reactance and leakage inductance of each part in the dq axis are integrated respectively. Finally, the right side of the obtained equation is added together, which is the potential energy converted in the inductance of the synchronous generator expressed in flux and voltage.

[0098] Based on the relationship between the generator terminal current, voltage, voltage phase angle at its node and system power angle, and after all parameters are converted to per-unit values, the rotor motion equation can be directly used with mechanical power P m and electromagnetic power P m Instead of mechanical torque T m and electromagnetic torque T e , so the potential energy function W of the synchronous generator can be obtained SynL As shown below:

[0099]

[0100] Where V represents the generator terminal voltage, θ represents the voltage phase angle at the generator node, δ represents the power angle, ψ D , md , f , Q , mq and ψ g They represent the flux linkage of the damping winding D, d-axis, the excitation winding, the damping winding Q, q-axis and the damping winding g, respectively. Dl 、x fl 、x Ql 、x gl 、x l 、x md and x mqThey represent the damping winding D, the excitation winding, the damping winding Q, g and the stator leakage inductance, as well as the mutual inductance between the dq axis windings, u f and R f They represent the voltage and impedance of the excitation winding respectively. When the subscript s is used in the text, it refers to the corresponding parameters when the system is in a stable state.

[0101] Because the randomness of synchronous generators is not obvious, they obey the equation of motion and kinetic energy formula both in stable operation and instability, and the synchronous generator is in the output line, and the rotor side is not connected to the transmission network, so there is no need to calculate its electromagnetic power. At the same time, the calculation process takes the per-unit value. Based on its equation of motion without considering energy dissipation, both sides are integrated to obtain the energy function represented by the kinetic energy of a single synchronous generator and then add it to the potential energy function of the synchronous generator. Assume that the number of synchronous generators in the power system is n g The overall energy function W of the synchronous generator part is obtained. Syn As shown below:

[0102]

[0103] Where M represents the inertia time constant, ω r represents the rotor speed, and D represents the damping coefficient.

[0104] For the dissipated energy on the transmission network side and the load side, when a fault occurs in the transmission network, the amplitude and phase angle of the node voltage, and the active and reactive power state variables of the line will change. From the energy function, these changes represent the conversion of the overall kinetic energy of the power grid into potential energy, and these potential energies will eventually be absorbed by the transmission network and the load. Therefore, the impact of the potential energy of the transmission network on the overall energy function of the system cannot be ignored. On high-voltage transmission lines, the conductivity is small enough to be ignored. Through flow calculation and related derivation and simplification, the energy function W of the transmission network is obtained. Line As shown below:

[0105]

[0106] Among them, V i represents the voltage on node i, θ i represents the voltage phase angle of node i, B ij represents the susceptance between nodes i and j. g The first node is a synchronous generator node, the nth g +1 node is a balance node, nth g +2 to nth g +m+1 nodes are load nodes.

[0107] On the load side, changes in the system as a whole will also affect the load's consumption of active and reactive power, which is still manifested as the conversion and absorption of potential energy from the perspective of energy function. Therefore, the potential energy absorbed in the load W Load The expression is as follows:

[0108]

[0109] Among them, Q Lml , U ml , P Lml ,θ ml They represent the reactive power consumption, voltage, active power consumption and voltage phase angle at the mlth load node respectively. Generally speaking, the active and reactive power injection and consumption of the power grid system are dynamically balanced.

[0110] For DFIG, the time-varying differential equations obtained through the three-phase stationary coordinate system are very difficult to analyze the operating state of the generator. And as for the transient energy function, in the transient process, θ is constantly changing, but the magnetomotive force generated in different coordinate systems is consistent, so the coordinate transformation method can be used.

[0111] The three-phase symmetrical stationary winding rotates at a synchronous speed. When connected to a three-phase power supply, a sinusoidally distributed rotating magnetic motive force F will be generated. Under the condition that F remains unchanged, the three-phase stationary coordinate system can be transformed into a two-phase stationary coordinate system. After precise mathematical processing, the time-varying state quantity in the stationary coordinate system and the time-invariant state quantity in the rotating coordinate system can be considered equivalent. Only for the dq two-phase coordinate system, the same rotating magnetic motive force can be obtained by making it rotate at a synchronous speed. In this coordinate system, the state quantity is converted from time-varying to time-invariant, so the coordinate transformation is divided into two steps.

[0112] Step 1, 3s / 2s transformation: Assume that the effective number of turns of each phase of the three-phase winding is N3, and the effective number of turns of each phase of the two-phase winding is N2. According to the principle of equal instantaneous magnetomotive force, and considering that the power is unchanged, the conversion into matrix form is shown as follows:

[0113]

[0114] Step 2, 2s / 2r transformation: If the d-axis and q-axis currents are set to i d 、i q , and then C 2r / 2s Taking the inverse, we get the transformation matrix C from the two-phase stationary coordinate system to the two-phase rotating coordinate system 2s / 2r As shown below:

[0115]

[0116] Combining the transformation matrix obtained by 3s / 2s transformation and 2s / 2r transformation, the equivalent circuit diagram of the corresponding DFIG is obtained. Figure 1 It is the d and q axis equivalent circuit diagram of DFIG, such as Figure 1 It should be noted that for the rotor winding itself, the angle leading the stator winding is θ r .in accordance with Figure 1 , the voltage, current and flux equations of DFIG in the two-phase rotating coordinate system can be obtained, and then the differential equations of the flux in all directions are obtained by substituting them into the voltage equation and transforming them as shown in the following formula:

[0117]

[0118] Among them, ψ sd , sq , rd , rq Respectively represent the magnetic flux components of the stator and rotor windings on the d and q axes, R s and R r Respectively represent the stator and rotor winding resistance, u sd 、u sq 、u rd 、u rq Represent the voltage components of the stator and rotor windings on the d and q axes, ω1 and ω s Represent stator speed and slip speed respectively, L ls and L lr They represent the stator and rotor winding self-inductance and the mutual inductance between groups L m The difference.

[0119] There is dissipated energy of magnetic field in the inductance of DFIG, which is also potential energy converted from kinetic energy in transient process. Incorporating this potential energy into the energy function model can make the model more applicable. According to the characteristics of energy, adding the right sides of the four equations after integration of the above formula can obtain the potential energy W converted in the DFIG inductance represented by magnetic flux and voltage. L As shown below:

[0120]

[0121] According to the actual situation, the voltage change and transmission network loss caused by stator transient, stator resistance, and stator speed change can be ignored. Substitute the generator terminal voltage, the voltage phase angle and power angle at the node into the stator dq axis current expression to obtain the active and reactive output of DFIG. Due to the power balance of the system, for the active part ΔP=P l +P k -P out , for the reactive part ΔQ=Q l +Q k -Qout , where P l and Q l Represent the active and reactive losses of the system, P k and Q k Respectively represent the active power and reactive power consumed by the load. In fact, ∫P out dδ represents the potential energy converted from kinetic energy, part of which is expressed as W L is stored in the inductor in the form of ∫P out Therefore, when the loss is ignored and the per-unit value is taken, the active output P out The output energy W corresponding to the system active power P It can be expressed as follows:

[0122]

[0123] Among them, s r is the generator slip, P s =u sd i sd +u sq i sq ,i sd and i sq Represent the current components of the stator and dq axes respectively. For reactive output Q out There is no simplified formula like the above, so it is necessary to calculate the rotor side reactive output Q r First, the reactive output of the rotor side is converted to the stator side. At the same time, the stator d and q axis currents are represented by vectors and substituted into the reactive output Q represented by the stator side. s After taking the per-unit value, we can get the output energy Q corresponding to the overall reactive output and system reactive power out It can be expressed as follows:

[0124]

[0125] Since the grid directly consumes reactive power, and the change of reactive power leads to voltage change, which affects the voltage stability of the system, that is, the potential energy corresponding to the reactive power output of the generator is directly related to the voltage. Therefore, this part of transient energy can be expressed by ∫Q out / VdV is used to express it.

[0126] Based on the operating mechanism of DFIG, its active power can be expressed as P m =P S +P R , where P m Represents mechanical input power, P S , P R denote the power output on the stator and rotor sides respectively. And P S =Te ω s , P m =T m ω r , where T e represents the electromagnetic torque on the rotor side, ω s Indicates the stator speed, T m represents the mechanical torque on the rotor side, ω r Indicates the rotor speed.

[0127] When DFIG is running stably, T e =T m , the power output on the rotor side can be expressed as When it enters the transient process, the motion equations of the wind turbine, drive shaft and doubly fed generator rotor side are as follows:

[0128]

[0129] Where J is the moment of inertia of the DFIG, D is the damping coefficient on the rotor side, θ r represents the rotation angle, ω m represents the speed of the wind turbine, n p represents the number of pole pairs of the doubly-fed generator, L m represents the mutual inductance between the stator and rotor windings, i sd 、i sq 、i rd 、i rq represents the current components of the stator and rotor windings on the d-axis and q-axis respectively. Therefore, set M = J / n p , D T =D / n P Based on the first and third formulas of the above formula, the transient process expression of the rotor motion equation can be obtained as shown below:

[0130] -D T Δω r dΔω r =MΔω r dΔω r -T m dθ r +T e dθ r

[0131] By integrating both sides of the equation above, we can obtain the transient energy function W represented by the rotor motion equation: R , as shown below:

[0132]

[0133] Among them, W keIt represents the kinetic energy of DFIG in transient process, W pe It represents the potential energy converted from kinetic energy of DFIG in transient process. Therefore, by substituting the motion equations of wind turbine, transmission shaft and rotor side of doubly fed generator into the above equation, it can be transformed into the following equation:

[0134]

[0135] Figure 2 is the relationship between DFIG rotor speed and output power. Figure 2 It can be seen that the input power and rotor speed of DFIG are linear in the AB and CD parts. Therefore, the ratios set in the grid-connected startup area and the constant speed area are expressed as K B and K N In the BC part, based on the coordinates of point B (P outB ,ω rB ) and the coordinates of point C (P outC ,ω rC ), the present invention uses a power exponential function to fit it. Therefore, the output power of DFIG can be defined as shown in the following formula:

[0136]

[0137] Assuming that the energy loss ΔP of the DFIG itself is approximately a constant, the relationship between the input power and the output power of the power system can be expressed as P m =P out +ΔP. Therefore, combining the motion equations of the wind turbine, drive shaft and DFIG rotor side and its output power, the DFIG rotor speed can be defined as follows:

[0138]

[0139] Among them, the input power of DFIG can be expressed as ρ represents the ambient gas density, R represents the radius of the windward surface of the DFIG, v represents the wind speed, C p (λ, β) represents the wind energy utilization coefficient, where λ and β represent the blade tip speed ratio and pitch angle, respectively. Therefore, when the DFIG operates in the grid-connected startup area and the constant speed area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows:

[0140]

[0141] Among them, K i To represent K B and K NWhen the DFIG operates in the maximum wind energy tracking area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows:

[0142]

[0143] set up is the proportional coefficient of the DFIG input part, so the above two equations can be simplified and uniformly expressed as the random transient energy function represented by kinetic energy inside the DFIG as shown below.

[0144]

[0145] Next, the random transient energy function represented by kinetic energy inside the DFIG derived above is verified. The verification method of the random transient energy function of the rest is similar. Based on Hamilton's theorem, the rotor angle θ and the rotor side speed ω r The relationship is shown as follows:

[0146]

[0147] in, is the derivative of the rotor angle, is the derivative of the rotor speed. Based on Formula 1 and Formula 3 in the random transient energy function represented by kinetic energy inside the DFIG and the above formula, its state space expression can be defined as shown below:

[0148]

[0149] STEF could have been converted to From the above formula, we can see that R(x)≥0 always holds true. Therefore, its derivative with respect to time is It is consistent with Lyapunov stability theory. When DFIG is in the maximum power point tracking area, its Hamilton criterion can be expressed as shown below:

[0150]

[0151] Through the same steps as above, it can be obtained that the random transient energy function represented by kinetic energy inside the DFIG still obeys the Lyapunov stability theory when the DFIG operates in the maximum wind energy tracking area. Therefore, it can be used as a part of the random transient energy function of the power system containing the DFIG.

[0152] In summary, through the above analysis, all transient energy functions are integrated. For the transient energy function of the entire power grid system, assuming that there are s wind turbines in the power grid, n g synchronous generators, a total of s+n ggenerators, of which one is a synchronous generator as a balancing generator. The overall random transient energy function can be expressed as:

[0153]

[0154] In order to characterize the transient stability of the power system, the present invention constructs a transient stability assessment model of the power system containing a wind farm based on DFIG. The model is based on the above-mentioned stochastic transient energy function (STEF) and random matrix theory, and consists of a basic region (BR), a stochastic region (SR) and a hazard region (HR). The boundary between BR and SR is defined as the boundary between the stable operation of the system within an acceptable range and the obvious random fluctuations of the node parameters. The boundary between SR and HR is defined as the boundary between the obvious random fluctuations of the node parameters and the large-scale fluctuations of the node parameters leading to system instability. The center of the stable region is defined as the trap point (TP). In addition, the present invention proposes a new bus stable nominal value (BSNV) as an evaluation criterion. Figure 3 The flowchart of the specific implementation method of the present invention for constructing the stability domain of the power system containing DFIG based on the random transient energy function and random matrix theory is as follows. Figure 3 shown.

[0155] exist Figure 3 Based on the data obtained from the wind farm and the derived STEF, the mathematical model of STEF can be obtained as W k , arranged according to the node number as W kΓ As the input layer of the calculation stability domain, k represents the time series, Γ represents the node number, the eigenvalue matrix corresponding to the input layer is λ and tr(λ) is the trace of the eigenvalue matrix. Therefore, the distance between the HR boundary and TP can be defined as follows:

[0156]

[0157] n represents the number of nodes in the DFIG-based power system with wind farms. For the boundary of SR, a new random matrix can be obtained by substituting the parameters obtained from the stable operation of the system into STEF. Then repeat the above steps to calculate the boundary of SR, and finally set the boundary of BR to 1;

[0158] In the BSNV calculation process, the active power and reactive power of each node are used to calculate the STEF of the corresponding node, so the corresponding random matrix W can be obtained. Γ and the corresponding maximum nominal value Load Γ , which can be used to describe the fluctuation of various parameters in the node. The degree of eccentricity of BSNV can be obtained through the following mathematical transformation:

[0159]

[0160] Among them, Λ represents a unitary matrix, Load Γ is the complex domain representation of the node stable nominal value, W Γ represents the STEF of the corresponding node in the time domain, represents its conjugate transposed matrix, W λ represents its eigenvalue matrix; the maximum value in the eigenvalue matrix is ​​defined as the BSNV of the corresponding node. After normalization, its horizontal and vertical coordinates in the power system stability domain composed of BR / SR / HR are represented as BSNV x =Re[Load Γ ] and BSNV y =Im[Load Γ ];

[0161] It is assumed that when more than 60% of BSNV enters HR, the DFIG-based power system with wind farm will be judged to be in an unstable state. Therefore, the position of the unstable equilibrium point (UEP) of the entire power system in the stable domain can be defined as shown in the following formula. It should be noted that it has the same characteristics as BSNV. Since BSNV is a standardized value, the following formula can be used to reflect the overall operation of the power system with wind farm:

[0162]

[0163] Among them, R UEP Indicates the distance between UEP and TP in the power system, Load Γ Represents the BSNV of node Γ; when UEP is in BR, it means that the power system is stable. When UEP enters SR, it can be determined that the randomness of the power system caused by the input tends to be obvious, and the various parameters of the system during operation will also fluctuate significantly, especially the parameters corresponding to the wind farm grid connection point. When UEP enters HR, it can be determined that the power system enters an unstable operation state.

[0164] Based on the above derivation process, the present invention proposes a method for constructing a random transient energy function and a system stability domain of a new energy power system based on Lyapunov transient stability theory. Figure 4 This is a flowchart of a specific implementation method of the method for constructing a random transient energy function and a system stability domain of a new energy power system based on Lyapunov transient stability theory. Figure 4 As shown in FIG. 1 , the specific steps of the method for constructing the random transient energy function and system stability domain of the new energy power system based on Lyapunov transient stability theory include:

[0165] S201: Get random source input wind speed sequence:

[0166] The wind speed sequence of the historical wind farm and the historical wind turbine speed sequence of the DFIG are then used to obtain the wind speed that is suitable for the system operation step size that develops over time based on the random wind speed model pre-set based on the basic conditions of the Weibull distribution.

[0167] Because the wind in a wind farm is random and volatile, the wind speed model generally uses a four-component superposition method, where the four components are basic wind, gust, gradual wind and random wind.

[0168] For basic wind, the specific value of basic wind can be obtained by analyzing the measured wind data and further determining the relevant parameters of Weibull distribution. It should be noted that basic wind directly affects the power delivered by wind turbines to the grid. Therefore, basic wind can be defined as V B =AΓ(1+1 / K), where A represents the proportionality coefficient, Γ(1+1 / K) represents the gamma function, and K represents the scale parameter of the Weibull distribution.

[0169] For gusts, the dramatic change in wind speed will inevitably affect the output power of the wind turbine, thereby affecting voltage fluctuations. According to its physical characteristics, gusts can be defined as V G =V S (t1≤t≤t1+T g ),in t、T G , t1 represent the actual time, gust duration and gust occurrence time respectively, and maxG represents the maximum value of the gust.

[0170] For the gradually changing wind, since it has a linear property within a certain period of time, it can be defined as V R =V γ (T 1R ≤t≤T 2R ), where V γ =maxR[1-(tT 2R ) / (T 1R -T 2R)], maxR is the maximum value of the gradient wind, T 1R , T 2R Respectively represent the occurrence time and end time.

[0171] For random wind, since the random wind speed changes at different heights are not the same, this paper uses Weibull distribution as the initial value of random wind, and then uses the Wiener process as the development trend of random wind. The mathematical expression of the change of random wind speed can be defined as. Therefore, the input wind speed of DFIG can be expressed as V N .

[0172] In order to improve the randomness of wind speed, the wind speed model in this embodiment adopts the synthesis of the above four wind speed types, namely basic wind, gust wind, gradual wind and random wind, and the expression is V = V B +V G +V R +V N

[0173] S202: The following method is used to construct a random transient energy function model of the traditional power system based on the equation of motion and energy conservation, including the random transient energy function models of the synchronous generator side, the transmission network side and the load side:

[0174] Based on the relationship between the generator terminal current, voltage, voltage phase angle at its node and system power angle, and after all parameters are converted to per-unit values, the rotor motion equation can be directly used with mechanical power P m and electromagnetic power P m Instead of mechanical torque T m and electromagnetic torque T e , so the potential energy function W of the synchronous generator can be obtained SynL As shown below:

[0175]

[0176] Where V represents the generator terminal voltage, θ represents the voltage phase angle at the generator node, δ represents the power angle, ψ D , md , f , Q , mq and ψ g They represent the flux linkage of the damping winding D, d-axis, the excitation winding, the damping winding Q, q-axis and the damping winding g, respectively. Dl 、x fl 、x Ql 、x gl 、x l 、x md and x mqThey represent the damping winding D, the excitation winding, the damping winding Q, g and the stator leakage inductance, as well as the mutual inductance between the dq axis windings, u f and R f They represent the voltage and impedance of the excitation winding respectively. When the subscript s is used in the text, it refers to the corresponding parameters when the system is in a stable state.

[0177] Assume that the number of synchronous generators in the power system is n g The overall energy function W of the synchronous generator part is obtained. Syn As shown below:

[0178]

[0179] Where M represents the inertia time constant, ω r represents the rotor speed, and D represents the damping coefficient.

[0180] For the dissipated energy on the transmission network side and the load side, when a fault occurs in the transmission network, the amplitude and phase angle of the node voltage, and the active and reactive power state variables of the line will change. From the energy function, these changes represent the conversion of the overall kinetic energy of the power grid into potential energy, and these potential energies will eventually be absorbed by the transmission network and the load. Therefore, the impact of the potential energy of the transmission network on the overall energy function of the system cannot be ignored. On high-voltage transmission lines, the conductivity is small enough to be ignored. Through flow calculation and related derivation and simplification, the energy function W of the transmission network is obtained. Line As shown below:

[0181]

[0182] Among them, V i represents the voltage on node i, θ i represents the voltage phase angle of node i, B ij represents the susceptance between nodes i and j. g The first node is a synchronous generator node, the nth g +1 node is a balance node, nth g +2 to nth g +m+1 nodes are load nodes.

[0183] On the load side, changes in the system as a whole will also affect the load's consumption of active and reactive power, which is still manifested as the conversion and absorption of potential energy from the perspective of energy function. Therefore, the potential energy absorbed in the load W Load The expression is as follows:

[0184]

[0185] Among them, Q Lml , U ml , P Lml,θ ml They represent the reactive power consumption, voltage, active power consumption and voltage phase angle at the mlth load node respectively. Generally speaking, the active and reactive power injection and consumption of the power grid system are dynamically balanced.

[0186] S203: The random transient energy function model of the DFIG part based on the equation of motion and energy conservation is constructed by the following method, which consists of four parts: the random transient energy function represented by the active output and the reactive output, the random transient energy function stored in the inductance inside the DFIG, and the random transient energy function represented by kinetic energy:

[0187] For DFIG in the transient process, θ is constantly changing, but the magnetomotive force generated in different coordinate systems is consistent, which meets the basic requirements of the coordinate transformation method. Considering that the state quantity is transformed from time-varying to time-invariant, the coordinate transformation is divided into two steps as shown below.

[0188] Step 1, 3s / 2s transformation: Assume that the effective number of turns of each phase of the three-phase winding is N3, and the effective number of turns of each phase of the two-phase winding is N2. According to the principle of equal instantaneous magnetomotive force, and considering that the power is unchanged, the conversion into matrix form is shown as follows:

[0189]

[0190] Step 2, 2s / 2r transformation: If the d-axis and q-axis currents are set to i d 、i q , and then C 2r / 2s Taking the inverse, we get the transformation matrix C from the two-phase stationary coordinate system to the two-phase rotating coordinate system 2s / 2r As shown below:

[0191]

[0192] Combining the transformation matrix obtained by 3s / 2s transformation and 2s / 2r transformation, the equivalent circuit diagram of the corresponding DFIG is obtained. Figure 1 It is the d and q axis equivalent circuit diagram of DFIG, such as Figure 1 It should be noted that for the rotor winding itself, the angle leading the stator winding is θ r .in accordance with Figure 1 , the voltage, current and flux equations of DFIG in the two-phase rotating coordinate system can be obtained, and then the differential equations of the flux in all directions are obtained by substituting them into the voltage equation and transforming them as shown in the following formula:

[0193]

[0194] Among them, ψ sd , sq , rd ,rq Respectively represent the magnetic flux components of the stator and rotor windings on the d and q axes, R s and R r Respectively represent the stator and rotor winding resistance, u sd 、u sq 、u rd 、u rq Represent the voltage components of the stator and rotor windings on the d and q axes, ω1 and ω s Represent stator speed and slip speed respectively, L ls and L lr They represent the stator and rotor winding self-inductance and the mutual inductance between groups L m The difference.

[0195] There is dissipated energy of magnetic field in the inductance of DFIG, which is also potential energy converted from kinetic energy in transient process. Incorporating this potential energy into the energy function model can make the model more applicable. According to the characteristics of energy, adding the right sides of the four equations after integration of the above formula can obtain the potential energy W converted in the DFIG inductance represented by magnetic flux and voltage. L As shown below:

[0196]

[0197] According to the actual situation, the voltage change and transmission network loss caused by stator transient, stator resistance, and stator speed change can be ignored. Substitute the generator terminal voltage, the voltage phase angle and power angle at the node into the stator dq axis current expression to obtain the active and reactive output of DFIG. Due to the power balance of the system, for the active part ΔP=P l +P k -P out , for the reactive part ΔQ=Q l +Q k -Q out , where P l and Q l Represent the active and reactive losses of the system, P k and Q k Respectively represent the active power and reactive power consumed by the load. In fact, ∫P out dδ represents the potential energy converted from kinetic energy, part of which is expressed as W L is stored in the inductor in the form of ∫P out Therefore, when the loss is ignored and the per-unit value is taken, the active output P out The output energy W corresponding to the system active power P It can be expressed as follows:

[0198]

[0199] Among them, s r is the generator slip, P s =u sd i sd +u sq i sq ,i sd and i sq Represent the current components of the stator and dq axes respectively. For reactive output Q out There is no simplified formula like the above, so it is necessary to calculate the rotor side reactive output Q r First, the reactive output of the rotor side is converted to the stator side. At the same time, the stator d and q axis currents are represented by vectors and substituted into the reactive output Q represented by the stator side. s After taking the per-unit value, we can get the output energy Q corresponding to the overall reactive output and system reactive power out It can be expressed as follows:

[0200]

[0201] Since the grid directly consumes reactive power, and the change of reactive power leads to voltage change, which affects the voltage stability of the system, that is, the potential energy corresponding to the reactive power output of the generator is directly related to the voltage. Therefore, this part of transient energy can be expressed by ∫Q out / VdV is used to express it.

[0202] Based on the operating mechanism of DFIG, its active power can be expressed as P m =P S +P R , where P m Represents mechanical input power, P S , P R denote the power output on the stator and rotor sides respectively. And P S =T e ω s , P m =T m ω r , where T e represents the electromagnetic torque on the rotor side, ω s Indicates the stator speed, T m represents the mechanical torque on the rotor side, ω r Indicates the rotor speed.

[0203] When DFIG is running stably, T e =T m , the power output on the rotor side can be expressed as When it enters the transient process, the motion equations of the wind turbine, drive shaft and doubly fed generator rotor side are as follows:

[0204]

[0205] Where J is the moment of inertia of the DFIG, D is the damping coefficient on the rotor side, θ r represents the rotation angle, ω m represents the speed of the wind turbine, n p represents the number of pole pairs of the doubly-fed generator, L m represents the mutual inductance between the stator and rotor windings, i sd 、i sq 、i rd 、i rq represents the current components of the stator and rotor windings on the d-axis and q-axis respectively. Therefore, set M = J / n p , D T =D / n P Based on the first and third formulas of the above formula, the transient process expression of the rotor motion equation can be obtained as shown below:

[0206] -D T Δω r dΔω r =MΔω r dΔω r -T m dθ r +T e dθ r

[0207] By integrating both sides of the equation above, we can obtain the transient energy function W represented by the rotor motion equation: R , as shown below:

[0208]

[0209] Among them, W ke It represents the kinetic energy of DFIG in transient process, W pe It represents the potential energy converted from kinetic energy of DFIG in transient process. Therefore, by substituting the motion equations of wind turbine, transmission shaft and rotor side of doubly fed generator into the above equation, it can be transformed into the following equation:

[0210]

[0211] Figure 2 is the relationship between DFIG rotor speed and output power. Figure 2 It can be seen that the input power and rotor speed of DFIG are linear in the AB and CD parts. Therefore, the ratios set in the grid-connected startup area and the constant speed area are expressed as K B and K N In the BC part, based on the coordinates of point B (P outB ,ωrB ) and the coordinates of point C (P outC ,ω rC ), the present invention uses a power exponential function to fit it. Therefore, the output power of DFIG can be defined as shown in the following formula:

[0212]

[0213] Assuming that the energy loss ΔP of the DFIG itself is approximately a constant, the relationship between the input power and the output power of the power system can be expressed as P m =P out +ΔP. Therefore, combining the motion equations of the wind turbine, drive shaft and DFIG rotor side and its output power, the DFIG rotor speed can be defined as follows:

[0214]

[0215] Among them, the input power of DFIG can be expressed as ρ represents the ambient gas density, R represents the radius of the windward surface of the DFIG, v represents the wind speed, C p (λ, β) represents the wind energy utilization coefficient, where λ and β represent the blade tip speed ratio and pitch angle, respectively. Therefore, when the DFIG operates in the grid-connected startup area and the constant speed area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows:

[0216]

[0217] Among them, K i To represent K B and K N When the DFIG operates in the maximum wind energy tracking area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows:

[0218]

[0219] set up is the proportional coefficient of the DFIG input part, so the above two equations can be simplified and uniformly expressed as the random transient energy function represented by kinetic energy inside the DFIG as shown below.

[0220]

[0221] In summary, through the above analysis, all transient energy functions are integrated. For the transient energy function of the entire power grid system, assuming that there are s wind turbines in the power grid, n g synchronous generators, a total of s+n ggenerators, of which one is a synchronous generator as a balancing generator. The overall random transient energy function can be expressed as:

[0222]

[0223] S204: The transient stability assessment model of the power system containing wind power plants based on DFIG and the stable nominal value of each node are constructed by the following method:

[0224] In order to characterize the transient stability of the power system, the present invention constructs a transient stability assessment model for the power system containing wind farms based on DFIG. The model is based on the above-mentioned stochastic transient energy function (STEF) and random matrix theory, and consists of a basic region (BR), a stochastic region (SR) and a hazard region (HR). The center of the stable region is defined as a trap point (TP). In addition, the present invention proposes a new bus stable nominal value (BSNV) as an evaluation criterion. Figure 3 The flowchart of the specific implementation method of the present invention for constructing the stability domain of the power system containing DFIG based on the random transient energy function and random matrix theory is as follows. Figure 3 shown.

[0225] exist Figure 3 Based on the data obtained from the wind farm and the derived STEF, the mathematical model of STEF can be obtained as W k , arranged according to the node number as W kΓ As the input layer of the calculation stability domain, k represents the time series, Γ represents the node number, the eigenvalue matrix corresponding to the input layer is λ and tr(λ) is the trace of the eigenvalue matrix. Therefore, the distance between the HR boundary and TP can be defined as follows:

[0226]

[0227] n represents the number of nodes in the DFIG-based power system with wind farms. For the boundary of SR, a new random matrix can be obtained by substituting the parameters obtained from the stable operation of the system into STEF. Then repeat the above steps to calculate the boundary of SR, and finally set the boundary of BR to 1;

[0228] In the BSNV calculation process, the active power and reactive power of each node are used to calculate the STEF of the corresponding node, so the corresponding random matrix W can be obtained.Γ and the corresponding maximum nominal value Load Γ , which can be used to describe the fluctuation of various parameters in the node. The degree of eccentricity of BSNV can be obtained through the following mathematical transformation:

[0229]

[0230] Among them, Λ represents a unitary matrix, Load Γ is the complex domain representation of the node stable nominal value, W Γ represents the STEF of the corresponding node in the time domain, represents its conjugate transposed matrix, W λ represents its eigenvalue matrix; the maximum value in the eigenvalue matrix is ​​defined as the BSNV of the corresponding node. After normalization, its horizontal and vertical coordinates in the power system stability domain composed of BR / SR / HR are represented as BSNV x =Re[Load Γ ] and BSNV y =Im[Load Γ ];

[0231] It is assumed that when more than 60% of BSNV enters HR, the DFIG-based power system with wind farm will be judged to be in an unstable state. Therefore, the position of the unstable equilibrium point (UEP) of the entire power system in the stable domain can be defined as shown in the following formula. It should be noted that it has the same characteristics as BSNV. Since BSNV is a standardized value, the following formula can be used to reflect the overall operation of the power system with wind farm:

[0232]

[0233] Among them, R UEP Indicates the distance between UEP and TP in the power system, Load Γ Represents the BSNV of node Γ. When UEP is in BR, it means that the power system is stable. When UEP enters SR, it can be determined that the randomness of the power system caused by the input tends to be obvious, and the various parameters of the system during operation will also fluctuate significantly, especially the parameters corresponding to the wind farm grid connection point. When UEP enters HR, it can be determined that the power system has entered an unstable operation state.

[0234] In order to illustrate the technical effect of the present invention, a specific example is used to experimentally verify the present invention. In this embodiment, the power grid model adopts a 3-machine 9-node standard power grid, and 50 2MW wind turbines with a capacity of 100MW are integrated into the power grid model at the 4th node through the bus mode. Since the embodiment involves different DFIG grid-connected capacities, a grid-connected switch is set at the external connection of each DFIG for control. Four types of faults are set on the grid side, namely: single-phase short circuit grounding fault, two-phase short circuit fault, two-phase short circuit grounding fault and three-phase short circuit grounding fault. When determining the type of fault, the fault point is uniformly set between nodes 7 / 8, and then a three-phase short circuit grounding fault is set. The positions of the fault points are: fault point 1 is placed between nodes 5 / 7 / 8, fault point 2 is placed between nodes 6 / 8 / 9, and fault point 3 is placed between nodes 4 / 5 / 6. The fault occurrence time and fault duration are both set to 2.0s and 0.3s. The potential energy boundary method is used as a comparison method to compare with the method for constructing the random transient energy function and system stability domain of the new energy power system based on Lyapunov transient stability theory of the present invention.

[0235] Figure 5 It is a curve diagram of the distribution trend of the random transient energy function and its envelope and energy spectrum of the power system containing a wind power plant under different DFIG grid-connected capacities after wind speed fluctuation and grid-side three-phase fault occur in the present invention. Figure 5 (a) shows the transient energy function envelope of the system containing only synchronous generators, (b), (c), (d) respectively show the random transient energy function envelope and energy spectrum when the grid-connected capacity of wind turbines accounts for 3%, 25%, and 50% of the grid-connected power generation of synchronous generators in the system. By comparison, it can be clearly seen that, first of all, for the traditional power system containing only synchronous generators, due to the lack of the influence of random sources, the random transient energy function will degenerate into the traditional transient energy function. At this time, the system has no randomness, such as Figure 5 As shown in (a) in the figure, as the proportion of wind turbine grid-connected capacity to system capacity increases, the changes in the overall random transient energy function of the system become more complex and dispersed, and the spacing between the upper and lower envelopes also increases. Figure 5It can be seen that whether it is a traditional power system or a power system after the wind turbine is connected to the grid, the impact of this fault on the entire system is almost the same. Random fluctuations will not exceed the impact of the power grid failure on the power system, indicating that although the random factors are expanded, they are still within the controllable range. When the wind turbine grid-connected capacity is 3% of the system, the lower envelope of the transient energy function is parallel to the envelope of the traditional energy function, indicating that under stable conditions, the lower limit of the transient energy function is not affected by random wind speed. With the increase of the wind turbine grid-connected capacity, the random changes of the upper and lower envelopes are more obvious. At the same time, the color in the energy spectrum from light to dark represents the increase in energy intensity. It can be more intuitively seen that with the increase of the wind turbine grid-connected capacity, the overall energy fluctuation of the system is also increasing, and the possible distribution positions of the random transient energy function at the same time are also more dispersed.

[0236] Figure 6 It is a distribution diagram of the stability domain and node stability nominal value of the power system containing wind power plants after wind speed fluctuation and grid-side three-phase fault under different DFIG grid-connected capacities obtained in the present invention based on random transient energy function and random matrix theory. Figure 6 (a), (b), (c) and (d) show that when the grid-connected capacity of wind turbines increases, the penetration rate of wind turbines into the original traditional power system as a whole becomes higher and higher. By depicting different circles with different colors, the blue circle represents BR, the green circle represents SR, and the red circle represents HR. Therefore, it can be seen that when the grid-connected capacity of wind turbines is 0%, the stability domain of the entire power system degenerates into the traditional stability domain. At this time, the width of SR is 0, indicating that the power system at this time is not random. The larger the capacity of wind turbines connected to the grid, the wider the width of SR, indicating that the randomness and volatility of wind speed penetrate the traditional power system through wind turbines, and the greater the impact on the overall stability of the power system. Therefore, the simultaneous existence of BR, SR and HR in the stability domain is more in line with the construction of the stability domain of the power system under the large-scale grid connection of new energy. And with the increase of wind turbine grid-connected capacity, the increase of SR width indicates that the grid-connected wind turbine increases the randomness of the power system. Even if there are nodes distributed in SR due to randomness, the system will not enter the adaptive state. The increase of HR width indicates that the area where faults may occur and cause the transient stability of the system to decrease increases, and the power system becomes more vulnerable. The decrease of BR width indicates that the connection of wind turbines to the grid will eventually reduce the range of instability in the stable domain of the power system, and at this time, the stable nominal values ​​of more nodes will be distributed in RC and HC.

[0237] Figure 7This is a comparative curve diagram of the voltage change trend of the grid connection point after wind speed fluctuation and grid-side three-phase fault under different DFIG grid connection capacities in the present invention. When DFIGs with different capacities are connected to the grid, once a power system fault occurs, the voltage fluctuation of the wind farm grid connection point is greater than the voltage fluctuation of other nodes, such as Figure 7 As shown by the red line in Figure 1, this means that the penetration of wind energy into the power system is intuitively represented by the voltage change at the wind farm connection point. Similarly, after the fault is eliminated, the system with higher DFIG capacity requires longer recovery time.

[0238] Figure 8 It is the motion trajectory of the stable domain of the DFIG-containing power system and its dominant unstable equilibrium point obtained based on the potential energy boundary method in the present invention. The maximum potential energy represents the maximum energy that can be converted into potential energy along the fault trajectory. Therefore, the fault can be cleared before the fault trajectory leaves the area of ​​the green line envelope, and all kinetic energy can be converted into potential energy, and the system can be determined to be stable. However, there are two problems that arise at the same time: (1) The construction of this stable domain is based on the reduced-order model. Although it can represent the reality of the original model to a certain extent, it is certainly not the stable domain directly obtained from the original model. Therefore, if the trajectory in the projection crosses the stable domain range of the green line envelope, it does not necessarily guarantee that the fault trajectory has left the stable boundary of the original system. On the contrary, if the trajectory in the projection moves within the stable domain range of the green line envelope, it does not necessarily guarantee that the fault trajectory is still in the stable domain of the original system. (2) For high-order systems, the characterization is more based on an "intuition" and lacks relevant theoretical support. In addition, describing each node in the entire system as a node itself involves a lot of approximations and simplifications, and the description of the details of each node is lost. The BC / RC / HC stability domain proposed in this paper is not only based on actual theoretical proof, but also through the characterization of each node in the stability domain, it can more intuitively determine the influence relationship between each node and the degree of influence on the overall power system.

[0239] Then the comparison design of different fault points is as described above. And the wind turbine generators in each group are set to 3% of the grid-connected capacity. Therefore, Table 2 shows the comparison table of the spectrum radius of each node in different fault locations and in the case of no fault.

[0240] Table 1 Summary of comparison of spectrum radius of each node at different fault locations

[0241]

[0242] By comparison, it can be clearly seen that when the fault point moves, the spectral radius of the node closest to the fault point will fluctuate more violently than the spectral radius of other nodes. For example, when the fault occurs at position 1, the node 7 closest to it fluctuates the most, and the inward contraction amplitude reaches 0.1575. At the same time, the fluctuation of node 5 connected by a 150km transmission line is the second largest. Therefore, the location of the three-phase short circuit grounding fault can be located between nodes 5 and 7. Similarly, when the fault occurs at position 2, the node 9 closest to it fluctuates the most, followed by node 6 connected by a 150km transmission line. Since the wind speed input of the grid-connected wind turbine is a random input, at the node position where the fault will not affect or has a small impact, the random fluctuations that affect the wind turbine dominate the dominant position. Therefore, the spectral radius of these nodes will fluctuate in a small range, resulting in a situation where the spectral radius increases even though the fault occurs. At the same time, when the fault occurs, the average spectral radius of all generator nodes and other nodes except the generator nodes will decrease, indicating that the overall stability of the power system containing wind turbines is reduced due to the fault. It can be noticed that when a three-phase grounding fault occurs at the node where the wind turbine is connected to the grid, there are both violent fluctuations and random fluctuations of the wind turbine. At this time, the overall grid system is extremely fragile, that is, at the fault point 3 in the table, both the spectral radius of other nodes and the spectral radius of the generator node are greatly reduced.

[0243] In the simulation analysis of different fault types, the fault point position is controlled to be the same, and fault type 1 is defined to represent single-phase short circuit grounding, fault type 2 represents two-phase short circuit, fault type 3 represents two-phase short circuit grounding, and fault type 4 represents three-phase short circuit grounding. Similarly, the spectral radius of other nodes in Table 1 represents line nodes and load nodes other than generator nodes. Through the embodiment, Table 2 is obtained to obtain a comparison summary table of the spectral radius of each node of different fault types.

[0244] Table 2 Summary of comparison of spectrum radius of each node for different fault types

[0245]

[0246] It can be clearly seen that as the fault becomes more serious, that is, from single-phase short circuit to grounding, two-phase short circuit to two-phase short circuit to grounding and the most serious three-phase short circuit to grounding, the average spectral radius of each node in the power system when the above four faults occur when the fault point is determined is smaller than the average spectral radius in the case of no fault, and the average spectral radius of the transmission circuit node and the load node shows a trend of gradually decreasing with the severity of the fault. At the same time, it can be noted that because the fault point is designed between nodes 5 and 7 in the simulation, the maximum fluctuation of the spectral radius of other nodes occurs at node 7, and the fluctuation is also more severe as the severity of the fault increases. For the generator node, the influence of two-phase short circuit fault and two-phase short circuit grounding fault on the maximum fluctuation of its spectral radius shows that the tolerance of the generator itself to two-phase short circuit grounding fault is higher than that of two-phase short circuit fault. In summary, the PEBS method lacks reasonable mathematical analysis and theoretical basis for power systems other than MIMO. It can only obtain a local approximation of the stability boundary related to the original model by describing its method steps, but it cannot consistently give a conservative stability estimate. Therefore, the stability domain characteristic evaluation model of the power system containing BC / HC / RC constructed in this paper is more in line with reality. Because the entire system is not abstracted as a motion trajectory, it is more intuitive and accurate in reflecting randomness.

[0247] Although the above describes the illustrative specific embodiments of the present invention to facilitate those skilled in the art to understand the present invention, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations using the concept of the present invention are protected.

Claims

1. A method for constructing random transient energy function and system stability domain of new energy power system based on Lyapunov transient stability theory. The following steps are involved: S1: Obtain the wind speed sequence of the historical wind farm and the historical wind turbine speed sequence of the DFIG, and then obtain the wind speed suitable for the system operation step length that develops over time according to the random wind speed model pre-set based on the basic conditions of the Weibull distribution to obtain the wind speed sequence; S2: The following method is used to construct the random transient energy function model of the traditional power system based on the equation of motion and energy conservation, including the random transient energy function models of the synchronous generator side, the transmission network side and the load side; Based on the relationship between the generator terminal current, voltage, voltage phase angle at its node and system power angle, and after all parameters are converted to per-unit values, the rotor motion equation can directly use the mechanical power P m and electromagnetic power P e Instead of mechanical torque T m and electromagnetic torque T e , so the potential energy function W of the synchronous generator can be obtained SynL As shown below: Where V represents the generator terminal voltage, θ represents the voltage phase angle at the generator node, δ represents the power angle, ψ D , md , f , Q , mq and ψ g They represent the flux linkage of the damping winding D, d-axis, the excitation winding, the damping winding Q, q-axis and the damping winding g, respectively. Dl 、x fl 、x Ql 、x gl 、x l 、x md and x mq They represent the damping winding D, the excitation winding, the damping winding Q, g and the stator leakage inductance, as well as the mutual inductance between the dq axis windings, u f and R f They represent the voltage and impedance of the excitation winding respectively. Meanwhile, the subscript s in the text refers to the corresponding parameters in the stable state of the system. Assume that the number of synchronous generators in the power system is n g The overall energy function W of the synchronous generator part is obtained. Syn As shown below: Where M represents the inertia time constant, ω r Indicates the rotor speed; For the dissipated energy on the transmission network side and the load side, when a fault occurs in the transmission network, the amplitude and phase angle of the node voltage, the active and reactive power state variables of the line will change. From the energy function, these changes represent the conversion of the overall kinetic energy of the power grid into potential energy, and these potential energies will eventually be absorbed by the transmission network and the load. Therefore, the influence of the potential energy of the transmission network on the overall energy function of the system cannot be ignored. On high-voltage transmission lines, the conductivity is small enough to be ignored. Through power flow calculation and related derivation and simplification, the energy function W of the transmission network is obtained. Line As shown below: Among them, V i represents the voltage on node i, θ i represents the voltage phase angle of node i, B ij represents the susceptance between nodes i and j; the first n g The first node is a generator node, the nth g +1 node is a balance node, nth g +2 to nth g +m+1 nodes are load nodes; On the load side, changes in the system as a whole will also affect the load's consumption of active power and reactive power, which is still manifested as the conversion and absorption of potential energy from the perspective of energy function. Therefore, the expression of the potential energy absorbed in the load is shown as follows: Among them, Q Lml , U ml , P Lml ,θ ml They represent the reactive power consumption, voltage amplitude, active power consumption and voltage phase angle at the mlth load node respectively. Generally speaking, the active and reactive power injection and consumption of the power grid system are dynamically balanced. S3: The random transient energy function model of the DFIG part based on the equation of motion and energy conservation is constructed by the following method, which consists of four parts: the random transient energy function represented by active output and reactive output, the random transient energy function stored in the inductance inside the DFIG, and the random transient energy function represented by kinetic energy; For DFIG, θ is constantly changing during the transient process, but the magnetomotive force generated in different coordinate systems is consistent, which meets the basic requirements of the coordinate transformation method. Considering that the state quantity is transformed from time-varying to time-invariant, the coordinate transformation is divided into the following two steps: Step 1, 3s / 2s transformation: Assume that the effective number of turns of each phase of the three-phase winding is N3, and the effective number of turns of each phase of the two-phase winding is N2. According to the principle of equal instantaneous magnetomotive force, and considering that the power is unchanged, the conversion into matrix form is shown as follows: Step 2, 2s / 2r transformation: If the d-axis and q-axis currents are set to i d 、i q , and then C 2r / 2s Taking the inverse, we get the transformation matrix C from the two-phase stationary coordinate system to the two-phase rotating coordinate system 2s / 2r As shown below: Combining the transformation matrix obtained by 3s / 2s transformation and 2s / 2r transformation, the equivalent circuit of the corresponding DFIG is obtained. It should be noted that for the rotor winding itself, the angle leading the stator winding is θ r ; Then, the voltage, current and flux equations of DFIG in the two-phase rotating coordinate system are established, and the differential equations of the flux in all directions are obtained by substituting the voltage equation and transforming it as shown in the following formula: Among them, ψ sd , sq , rd , rq Respectively represent the magnetic flux components of the stator and rotor windings on the d and q axes, R s and R r Respectively represent the stator and rotor winding resistance, u sd 、u sq 、u rd 、u rq Represent the voltage components of the stator and rotor windings on the d and q axes, ω 1 and ω s Represent stator speed and slip speed respectively, L ls and L lr They represent the stator and rotor winding self-inductance and the mutual inductance between groups L m The difference between There is dissipated energy of magnetic field in the inductance of DFIG, which is also potential energy converted from kinetic energy in transient process. Incorporating this potential energy into the energy function model can make the model more applicable. According to the characteristics of energy, the right sides of the four equations after integration are added together to obtain the potential energy W converted in the DFIG inductance expressed by magnetic flux and voltage. L As shown below: According to the actual situation, the voltage change and transmission network loss caused by stator transient, stator resistance, and stator speed change are ignored, and the generator terminal voltage, the voltage phase angle and power angle at the node are substituted into the stator dq axis current expression to obtain the active and reactive output of DFIG; due to the power balance of the system, for the active part ΔP=P l +P k -P out , for the reactive part ΔQ=Q l +Q k -Q out , where P l and Q l Represent the active and reactive losses of the system, P k and Q k Respectively represent the active power and reactive power consumed by the load; in fact, ∫P out dδ represents the potential energy converted from kinetic energy, part of which is expressed as W L is stored in the inductor in the form of ∫P out dθ is injected into the grid; therefore, when the loss is ignored and the per-unit value is taken, the active output P out The output energy W corresponding to the system active power P It can be expressed as follows: Among them, s r is the generator slip, P s =u sd i sd +u sq i sq ,i sd and i sq Represent the current components of the stator and dq axes respectively. For the reactive output Q out There is no simplified formula like the above, so it is necessary to calculate the rotor side reactive output Q r First, the reactive output of the rotor side is converted to the stator side. At the same time, the stator d and q axis currents are represented by vectors and substituted into the reactive output Q represented by the stator side. s After taking the per-unit value, we can get the output energy Q corresponding to the overall reactive output and system reactive power out It can be expressed as follows: Since the grid directly consumes reactive power, and the change of reactive power leads to voltage change, which affects the voltage stability of the system, that is, the potential energy corresponding to the reactive power output of the generator is directly related to the voltage. Therefore, this part of transient energy can be expressed by ∫Q out It is expressed as VdV; Based on the operating mechanism of DFIG, its active power can be expressed as P m =P S +P R , where P m Represents mechanical input power, P S , P R denote the power output on the stator and rotor sides respectively; and P S =T e ω s , P m =T m ω r , where T e represents the electromagnetic torque on the rotor side, ω s Indicates the stator speed, T m represents the mechanical torque on the rotor side, ω r Indicates the rotor speed; When DFIG is running stably, T e =T m , the power output on the rotor side can be expressed as When it enters the transient process, the motion equations of the wind turbine, drive shaft and doubly fed generator rotor side are as follows: Where J is the moment of inertia of the DFIG, D is the damping coefficient on the rotor side, θ r represents the rotation angle, n p represents the number of pole pairs of the doubly-fed generator, L m represents the mutual inductance between the stator and rotor windings, i sd 、i sq 、i rd 、i rq Represent the current components of the stator and rotor windings on the d-axis and q-axis respectively, so set M = J / n p , D T =D / n P Based on the first and third formulas of the above formula, the transient process expression of the rotor motion equation can be obtained as shown below: -D T Give r dSee r =MDo r dSee r -T m dth r +T e dth r By integrating both sides of the equal sign in the above equation, we can obtain the transient energy function W represented by the rotor motion equation: R , as shown below: Among them, W ke It represents the kinetic energy of DFIG in transient process, W pe It represents the potential energy converted from kinetic energy of DFIG in transient process; therefore, by substituting the motion equations of wind turbine, transmission shaft and rotor side of doubly fed generator into the above equation, it can be transformed into the following equation: Since the input power and rotor speed of DFIG are linear in the AB and CD parts, the ratios set in the grid-connected startup area and the constant speed area are expressed as K B and K N ; In the BC part, based on the coordinates of point B (P outB ,ω rB ) and the coordinates of point C (P outC ,ω rC ), and a power exponential function is used to fit it, so the output power of DFIG can be defined as follows: Assuming that the energy loss ΔP of the DFIG itself is approximately a constant, the relationship between the input power and the output power of the power system can be expressed as P m =P out +ΔP, therefore, combining the equations of motion of the wind turbine, drive shaft and DFIG rotor side and its output power, the DFIG rotor speed can be defined as follows: Among them, the input power of DFIG can be expressed as ρ represents the ambient gas density, R represents the radius of the windward surface of the DFIG, v represents the wind speed, C p (λ, β) represents the wind energy utilization coefficient, where λ and β represent the tip speed ratio and pitch angle, respectively. Therefore, when the DFIG operates in the grid-connected startup area and the constant speed area, the random transient energy function represented by kinetic energy inside the DFIG can be defined as follows: Among them, K i To represent K B and K N , and when DFIG operates in the maximum wind energy tracking area, the random transient energy function represented by kinetic energy inside DFIG can be defined as follows: set up is the proportional coefficient of the DFIG input part, so the above two equations can be simplified and uniformly expressed as the random transient energy function represented by kinetic energy inside the DFIG as shown below: In summary, through the above analysis, all transient energy functions are integrated. For the transient energy function of the entire power grid system, assuming that there are s wind turbines in the power grid, n g synchronous generators, a total of s+n g generators, and a synchronous generator as a balance generator; thus the overall random transient energy function can be expressed as: S4: The transient stability assessment model of the power system containing wind power plants based on DFIG and the stable nominal value of each node are constructed by the following method; In order to characterize the transient stability of the power system, a transient stability assessment model of the power system containing wind farms based on DFIG is constructed. The model is based on the above-mentioned stochastic transient energy function (STEF) and random matrix theory, and consists of a basic region (BR), a stochastic region (SR) and a hazard region (HR). The center of the stable region is defined as the trap point (TP). In addition, a new bus stable nominal value (BSNV) is proposed as an evaluation criterion. The stability domain of the power system containing DFIG is constructed based on random transient energy function and random matrix theory: Based on the data obtained from the wind farm and the derived STEF, the mathematical model that can be used to obtain the STEF is W k , arranged according to the node number as W kΓ As the input layer of the calculation stability domain, k represents the time series, Γ represents the node number, the eigenvalue matrix corresponding to the input layer is λ and tr(λ) is the trace of the eigenvalue matrix. Therefore, the distance between the HR boundary and TP can be defined as follows: n represents the number of nodes in the DFIG-based power system with wind farms; for the boundary of SR, a new random matrix can be obtained by substituting the parameters obtained by the stable operation of the system into STEF; then repeat the above steps to calculate the boundary of SR, and finally set the boundary of BR to 1; In the BSNV calculation process, the active power and reactive power of each node are used to calculate the STEF of the corresponding node, so the corresponding random matrix W can be obtained. Γ and the corresponding maximum nominal value Load Γ , which can be used to describe the fluctuation of various parameters in the node. The degree of eccentricity of BSNV can be obtained through the following mathematical transformation: Among them, Λ represents a unitary matrix, Load Γ is the complex domain representation of the node stable nominal value, W Γ represents the STEF of the corresponding node in the time domain, represents its conjugate transposed matrix, W λ represents its eigenvalue matrix; the maximum value in the eigenvalue matrix is ​​defined as the BSNV of the corresponding node. After normalization, its horizontal and vertical coordinates in the power system stability domain composed of BR / SR / HR are represented as BSNV x =Re[Load Γ ] and BSNV y =Im[Load Γ ]; It is assumed that when more than 60% of BSNV enters HR, the DFIG-based power system with wind farm will be judged to be in an unstable state. Therefore, the position of the unstable equilibrium point (UEP) of the entire power system in the stable domain can be defined as shown in the following formula. It should be noted that it has the same characteristics as BSNV. Since BSNV is a standardized value, the following formula can be used to reflect the overall operation of the power system with wind farm: Among them, R UEP Indicates the distance between UEP and TP in the power system, Load Γ Represents the BSNV of node Γ; when UEP is in BR, it means that the power system is stable. When UEP enters SR, it can be determined that the randomness of the power system caused by the input tends to be obvious, and the various parameters of the system during operation will also fluctuate significantly. When UEP enters HR, it can be determined that the power system enters an unstable operation state.

Images