CN106294993A - A kind of transient energy function analysis method considering that inverter current is saturated - Google Patents

Abstract

The present invention provides a kind of and considers that the transient energy function that inverter current is saturated analyzes method, based on inverter virtual synchronous generator control policy grounds, by setting up the virtual synchronous generator amature equation of motion and considering the current saturation factor of inverter, first integral method is utilized to build transient energy function, for in transient energy function and path-dependent can not product term use linear path method approximation, then use BCU method to obtain critical energy value, and obtain critical clearing time.The present invention considers that virtual rotation inertia and the saturated transient energy function of inverter current are analyzed the stability to the power system containing virtual synchronous electromotor that method can be quantitative and be estimated, the transient energy function constructed considers the saturated factor of inverter, closer to reality, solve traditional online evaluation conservative and cause the problem of system instability, be a kind of important supplement to time-domain simulation method.

Description

A kind of transient energy function analysis method considering that inverter current is saturated

Technical field

The present invention designs a kind of power system stability and control technology, a kind of considers saturated temporary of inverter current State energy function analyzes method.

Background technology

At present, along with the electric power electricity that is continuously increased and distributed power source is connected with electrical network of distributed power source permeability Sub-device uses the mode of Digital Circuit Control, and transient response speed is fast, and almost without inertia, when DG accounts for certain proportion Time, the disturbance that electrical network is little can cause safety and stability problem.In this context, Chinese scholars proposes virtual synchronous generating Machine technology, achieves the friendly access of distributed power source by controlling the operation logic of inverter simulation synchronous generator.

The research Chinese scholars focus of virtual synchronous electromotor is simulated synchronous generator inertia at combining inverter And damping, although some scholar considers virtual rotation inertia inverter transient energy function analysis, but all have ignored with void The impact on transient stability of the current saturation characteristic of the droop control inverter of plan inertia, the critical clearing time so obtained The most conservative, when system online evaluation can because the mute time long thus cause system this significant problem unstable, Directly influencing national economy, therefore, the instability problem caused for conservative needs to consider the saturated factor of inverter current Impact on transient stability.

Summary of the invention

It is an object of the invention to provide a kind of transient energy function considering that inverter current is saturated and analyze method, solution Traditional consider that virtual rotation inertia inverter transient energy function is analyzed the critical clearing time that method obtains and too protected Keep, when system online evaluation can because the mute time long thus cause system this significant problem unstable.

This method comprises the following steps:

Step 1, sets up the detail mathematic model of the power system containing virtual synchronous electromotor, this model bag to inverter Include and be similar to two order mode type equation of rotor motion of conventional electric generators and nonlinear load equation and network equation；

Step 2, utilizes first integral method to build transient energy function, equation of rotor motion is write as first-order system, profit By first integral method, both sides are integrated simultaneously, and saturated for inverter current factor is taken into account structure meter and the energy of damping Amount type liapunov function, in transient energy function and path-dependent can not product term use linear path method approximation；

Step 3, uses BCU method to obtain critical energy value: when first calculating fault, path is to the fault clearance moment and extremely Exit point, utilizes the parameter in fault clearance moment to calculate the energy in this moment, then calculates the minimal gradient point of post-fault system And obtain leading unstable uneven point according to solving power balance equation, and by this dotted state gain of parameter critical energy value；

Step 4, solves acquisition critical clearing time when transient state energy is equal to the energy at controlling unstable equilibrium point.

The present invention tries to achieve the critical fault mute time and carries out with traditional saturated transition energy obtained of inverter of ignoring Contrast is closer in reality.

Below in conjunction with Figure of description, the present invention is described further.

Accompanying drawing explanation

Fig. 1 is to consider that virtual rotation inertia inverter transient energy function analyzes method flow diagram.

Fig. 2 is three virtual synchronous electromotor Psychotria rubra (Lour.) Poir. dot system analogous diagram.

Fig. 3 is that BCU method obtains transition energy and CCT flow chart.

Fig. 4 is system gross energy and kinetic energy and potential energy simulation result figure.

Detailed description of the invention

The invention discloses and a kind of consider virtual rotation inertia and the saturated transient energy function analysis side of inverter current Method, based on inverter virtual synchronous generator control policy grounds, by setting up the virtual synchronous generator amature equation of motion And consider the current saturation factor of inverter, utilize first integral method to build transient energy function, in transient energy function With can not product term employing linear path method approximating of path-dependent, then use BCU method to obtain critical energy value, finally lead to The crash time crossing acquisition critical clearing time and time-domain-simulation acquisition is contrasted.The consideration virtual rotation inertia of the present invention The transient energy function saturated with inverter current analyze that method can be quantitative to the power train containing virtual synchronous electromotor The stability of system is estimated, and the transient energy function constructed considers the saturated factor of inverter, closer to reality, solves Traditional online evaluation conservative causes the problem that system is unstable, is a kind of important supplement to time-domain simulation method.Tool Body embodiment is following (as shown in Figure 1):

Step 1, sets up the detail mathematic model of the power system containing virtual synchronous electromotor, this model bag to inverter Include and be similar to two order mode type equation of rotor motion of conventional electric generators and nonlinear load equation and network equation；

Step 2, utilizes first integral method to build transient energy function, equation of rotor motion is write as first-order system, profit By first integral method, both sides are integrated simultaneously, and saturated for inverter current factor is taken into account structure meter and the energy of damping Amount type liapunov function, in transient energy function and path-dependent can not product term use linear path method approximation；

Step 3, uses BCU method to obtain critical energy value: when first calculating fault, path is to the fault clearance moment and extremely Exit point, utilizes the parameter in fault clearance moment to calculate the energy in this moment, then calculates the minimal gradient point of post-fault system And obtain leading unstable uneven point according to solving power balance equation, and by this dotted state gain of parameter critical energy value；

Step 4, solves acquisition critical clearing time when transient state energy is equal to the energy at controlling unstable equilibrium point.

Specifically

The first step, sets up the detail mathematic model containing power system, including:

1. the equation of rotor motion of virtual synchronous electromotor is set up.Inverter uses virtual synchronous generator control strategy, The two order mode type equation of rotor motion that foundation is similar to conventional electric generators are as follows；

$\begin{array}{c}\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\θ}}}_i={\widetilde{\mathrm{\ω}}}_i\\ M_i{\stackrel{\·}{\widetilde{\mathrm{\ω}}}}_i=P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}\end{array}\right.\\ P_{emi}=\frac{E_i{E}_j}{\sqrt{E_i^2-2E_i{E}_j\cos ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)+E_j^2}}\\ I_{ij}=\left\{\begin{array}{c}{I}_{ij},I_{ij}\≤I_{\max }\\ I_{\max },I_{ij}>I_{\max }\end{array}\right.\end{array}---\left(1\right)$

In formula, E represents node voltage, θ_i、Be respectively i-th electromotor relative to the rotor angle in the center of inertia and Angular frequency, unit is respectively rad, rad s-¹、M_iBeing the inertia constant of i-th electromotor, unit is s²·rad-¹, P_MiRepresent The mechanical output of i-th virtual synchronous electromotor, P_emiRepresenting the electromagnetic power of i-th virtual synchronous electromotor, inverter is maximum Electric current is I_maxAndθ_i=δ_i-δ₀,

In formula, δ_iBeing the rotor velocity of i-th electromotor, unit is rad, ω_iFor angular frequency, unit is rad s^-1, P_COIFor the accelerating power in the center of inertia, m is the quantity of virtual synchronous electromotor, δ₀、ω₀It is respectively the rotor angle speed in the center of inertia Degree and angular frequency.

2. nonlinear load model

$\begin{array}{c}{P}_{li}=f_{pi}\left(E_i\right)\\ Q_{li}=f_{qi}\left(E_i\right)\end{array}---\left(2\right)$

In formula, P_li, Q_liIt is respectively the idle and active power of load absorption.

3. network equation

Active power and reactive power from node i injection network be:

$\begin{array}{c}0=P_i+\underset{j=1}{\overset{n}{\Σ}}{B}_{ij}{E}_i{E}_j\sin ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)\\ 0=Q_i+\underset{j=1}{\overset{n}{\Σ}}{B}_{ij}{E}_i{E}_j\cos ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)\\ i=m+1,\mathrm{...},n\end{array}---\left(3\right)$

In formula, P_i, Q_iIt is respectively the active power in node i injection network and reactive power.

Second step, obtains transient energy function.

First integral method is utilized to build transient energy function.Equation of rotor motion is write as first-order system, is utilized first Both sides are integrated by integration method simultaneously, and saturated for inverter current factor is taken into account structure meter and energy type Lee of damping Ya Punuofu function such as formula (4)～(7), in transient energy function and path-dependent can not product term use linear path method Approximation, detailed process is as follows:

$V({\widetilde{\mathrm{\ω}}}_{gi},\mathrm{\θ},E)=V_K+V_P---\left(4\right)$

V_P=V_P1+V_P2+V_P3+V_P4+V_damping (5)

V_P4=V_P41+V_P42+V_P43 (6)

$\left\{\begin{array}{c}{V}_K=\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{M}_i{\widetilde{\mathrm{\ω}}}_{gi}^2\\ V_{P1}=-\underset{i=1}{\overset{m}{\Σ}}{P}_{M_i}({\mathrm{\θ}}_i-{\mathrm{\θ}}_i^s)\\ V_{P2}=\underset{i=1}{\overset{n+m}{\Σ}}{P}_i({\mathrm{\θ}}_i-{\mathrm{\θ}}_i^s)\\ V_{P3}=\underset{i=1}{\overset{n+m}{\Σ}}\frac{Qi}{a{\left(E_i^s\right)}^a}({\left(E_i\right)}^a-{\left(E_i^s\right)}^a)\\ V_{P4}=\underset{i=1}{\overset{m}{\Σ}}\underset{j=i+1}{\overset{m}{\Σ}}\frac{E_i{E}_j}{\sqrt{E_i^2-2E_i{E}_j\cos ({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)+E_j^2}}{I}_{ij}\\ -\underset{i=1}{\overset{m}{\Σ}}\underset{j=i+1}{\overset{m}{\Σ}}\frac{E_i^s{E}_j^s}{\sqrt{{\left(E_i^s\right)}^2-2E_i^s{E}_j^s\cos ({\mathrm{\θ}}_i^s-{\mathrm{\θ}}_j^s)+{\left(E_j^s\right)}^2}}{I}_{ij}\\ -\underset{i=m}{\overset{n+m-1}{\Σ}}\underset{j=i+1}{\overset{n+m}{\Σ}}{B}_{ij}\left(E_i{E}_j\cos \right({\mathrm{\θ}}_i-{\mathrm{\θ}}_j)-E_i^s{E}_j^s\cos ({\mathrm{\θ}}_i^s-{\mathrm{\θ}}_j^s\left)\right)\\ -\frac{1}{2}\underset{i=1}{\overset{n+m}{\Σ}}{B}_{ii}(E_i^2-{\left(E_i^s\right)}^2)\\ V_{damping}=\underset{i=1}{\overset{m}{\Σ}}{\∫}_{t_s}^t{D}_i{\mathrm{\ω}}_i\frac{{\mathrm{d\θ}}_i}{dt}=\underset{i=1}{\overset{m}{\Σ}}{\∫}_{{\mathrm{\θ}}_{i,s}}^{{\mathrm{\θ}}_i}{D}_i{\mathrm{\ω}}_i{\mathrm{d\θ}}_i\end{array}\right.---\left(7\right)$

Wherein, V is described energy function, V_KFor the kinetic energy of virtual synchronous electromotor, V_PFor the potential energy that system is total, V_P1For All virtual synchronous generator mechanical power inputs the rotor potential energy caused, V_P2The potential energy caused for whole active loads, V_P3For The potential energy that all reactive load causes, V_P4For storage and the potential energy in network, behalf stable equilibrium point, i, j are the index of node Value, n is node number, D_iIt is the damping of i-th electromotor, B_ij、B_iiBe respectively transadmittance between node i, j and node i from Admittance, E_i、E_jIt is the voltage at node i, j respectively,It is the angular frequency of i-th electromotor, θ_iBe i-th electromotor relative to The rotor angle in the center of inertia, wherein a is that constant value is generally 2.

3rd step, obtains critical energy value.

In conjunction with Fig. 3, BCU method is used to obtain critical energy value step as follows

Step 3-1, for formula (1) electric power system model, writes out its pinch system such as formula (8)

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\θ}}}_i=P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}=f_i\left(\mathrm{\θ}\right)\\ {\mathrm{\θ}}_m=\frac{\underset{i=1}{\overset{m-1}{\Σ}}{M}_i{\mathrm{\θ}}_i}{M_m}\end{array}\right.---\left(8\right)$

Then path during fault is usedAsk for exit point θ^EP, it is that to there is pinch system steady for projection path Fixed the most borderline, formula (1) path when can obtain its fault.Detection exit point θ^EPIt is to be arrived first by projection path Individual local potential energy maximum.The θ value namely obtained by formula (1) brings power deviation amount equation after fault into:

$P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}=f_i\left(\mathrm{\θ}\right)---\left(9\right)$

When meeting condition f_i* d θ=0, namelyRear acquisition θ^EP.In order to correctly obtain exit point, detection Precision chooses 10^-5。

Step 3-2, with exit point θ EP as initial point, is integrated, along integration the pinch system of formula (8) gained Curve looks for first minima shown in formula (10):

$F\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{m}{\Σ}}{f}_i^2\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{m}{\Σ}}{\[P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}\]}^2---\left(10\right)$

First minima obtained is referred to as minimal gradient point θ^MGP。

Step 3-3, with minimal gradient point θ^EPFor initial value, by Newton-Raphson method iterative (m-1) individual fault Rear power deviation amount equation:

$\left\{\begin{array}{c}{f}_i\left(\mathrm{\θ}\right)=P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}=0\\ {\mathrm{\θ}}_m=\frac{\underset{i=1}{\overset{m-1}{\Σ}}{M}_i{\mathrm{\θ}}_i}{M_m}\end{array}\right.---\left(11\right)$

Obtain the controlling unstable equilibrium point CUEP of pinch system.The state parameter of CUEP is substituted into formula (7) transient state energy Function calculates critical energy value V_cr。

4th step, obtains critical clearing time and is estimated transient stability.

Make the transient state energy value of formula (7) gained equal to the critical energy value V obtained in the 3rd step_cr, thus try to achieve critical Fault clearing time go forward side by side line stabilization assessment.

The consideration virtual rotation inertia inverter transient energy function analysis method that the present invention proposes is analyzed by emulation Effectiveness.Electromotor based on traditional IEEE 3 machine 9 node system changes virtual synchronous electromotor into and containing current limit On device, respectively three phase short circuit fault is set at different circuits and is analyzed, obtain different critical clearing times and system Energy includes that the kinetic energy of system and potential energy are shown in accompanying drawing 4, and the CCT result obtained with time-domain-simulation contrasts.Analogue system figure Seeing accompanying drawing 2, the parameter of virtual synchronous electromotor is as shown in table 1, and the critical clearing time obtained is as shown in table 2.

The different virtual synchronous generator parameter of table 1

Parameter	Size
		Filter capacitor	35uF
Dead resistance	0.5′Ω
		Load	50′Ω
DC voltage	650V
		Line voltage	380V
Line impedance	0.01+j0.377′Ω
		PR proportionality constant kp	150
PR integration time constant ki	4
		Mains frequency f0	50Hz
Reactive-power control coefficient kq	1×10-4
		Voltage regulation coefficient kv	3.5×10-2
Rated capacity S1	247.5MVA
		Rated capacity S2	192.0MVA
Rated capacity S3	128.0MVA

Table 2 considers that the saturated different circuit three phase short circuit fault of inverter current obtains CCT

Malfunctioning node	Faulty line	CCT transient energy function	CCT time-domain-simulation
				4	4-5	0.160	0.162
4	4-6	0.155	0.160
				5	5-7	0.165	0.168
6	6-9	0.174	0.179
				7	7-8	0.175	0.18
8	8-9	0.130	0.135

Table 3 does not consider that the saturated different circuit three phase short circuit fault of inverter current obtains CCT

Malfunctioning node	Faulty line	CCT transient energy function	CCT time-domain-simulation
				4	4-5	0.164	0.168
4	4-6	0.159	0.164
				5	5-7	0.169	0.173
6	6-9	0.178	0.182
				7	7-8	0.180	0.184
8	8-9	0.136	0.140

Contrast according to table 2 and table 3 and understand, traditional transient energy function analysis not considering the saturated factor of inverter current The critical clearing time that method obtains is longer, if fault clearing time length is just than the time considering the saturated factor of inverter current Can cause the instability of system, and utilize critical clearing time (CCT) that transient energy function method obtains and time-domain-simulation to obtain The error arrived is the least, and also the checking present invention considers virtual rotation inertia and the saturated transient energy function analysis side of inverter current The effectiveness of method.Therefore, the inventive method, closer to reality, solves traditional online evaluation conservative and causes system unstable Problem, be a kind of important supplement to time-domain simulation method.

Claims (5)

1. consider that the transient energy function that inverter current is saturated analyzes a method, comprise the following steps:

Step 1, sets up the detail mathematic model of the power system containing virtual synchronous electromotor to inverter, and this model includes class It is similar to two order mode type equation of rotor motion of conventional electric generators and nonlinear load equation and network equation；

Step 2, utilizes first integral method to build transient energy function, equation of rotor motion is write as first-order system, utilize head Both sides are integrated by secondary integration method simultaneously, and saturated for inverter current factor is taken into account structure meter and the energy type of damping Liapunov function, in transient energy function and path-dependent can not product term use linear path method approximation；

Step 3, uses BCU method to obtain critical energy value: when first calculating fault, path is to the fault clearance moment and to outlet Point, utilizes the parameter in fault clearance moment to calculate the energy in this moment, then calculates minimal gradient point the root of post-fault system Leading unstable uneven point is obtained according to solving power balance equation, and by this dotted state gain of parameter critical energy value；

Step 4, solves acquisition critical clearing time when transient state energy is equal to the energy at controlling unstable equilibrium point.

Method the most according to claim 1, it is characterised in that the mathematical model in step 1 is

$\begin{array}{c}\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\θ}}}_i={\widetilde{\mathrm{\ω}}}_i\\ M_i{\stackrel{\·}{\widetilde{\mathrm{\ω}}}}_i=P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}\end{array}\right.\\ P_{emi}=\frac{E_i{E}_j}{\sqrt{E_i^2-2E_i{E}_j\cos \left({\mathrm{\θ}}_i-{\mathrm{\θ}}_j\right)+E_j^2}}\\ I_{ij}=\left\{\begin{array}{c}{I}_{ij},I_{ij}\≤I_{\max }\\ I_{\max },I_{ij}>I_{\max }\end{array}\right.\end{array}---\left(1\right)$

In formula, E represents node voltage, θ_i、It is respectively i-th electromotor relative to the rotor angle in the center of inertia and angular frequency, M_iIt is the inertia constant of i-th electromotor, P_MiRepresent the mechanical output of i-th virtual synchronous electromotor, P_emiRepresent i-th void Intend the electromagnetic power of synchronous generator, I_ijFor line current between node i and node j, inverter maximum current is I_max；

θ_i=δ_i-δ₀,In formula, δ_i It is the rotor velocity of i-th electromotor, ω_iFor angular frequency, P_COIFor the accelerating power in the center of inertia, m is virtual synchronous generating The quantity of machine, δ₀、ω₀It is respectively rotor velocity and the angular frequency in the center of inertia.

Method the most according to claim 2, it is characterised in that in step 2, saturated for inverter current factor is taken into account The energy type liapunov function of structure meter and damping is:

$V({\widetilde{\mathrm{\ω}}}_{gi},\mathrm{\θ},E)=V_K+V_P---\left(2\right)$

V_P=V_P1+V_P2+V_P3+V_P4+V_damping (3)

$\left\{\begin{array}{c}{V}_K=\frac{1}{2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{M}_i{\widetilde{\mathrm{\ω}}}_{gi}^2\\ V_{P1}=-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{P}_{M_i}\left({\mathrm{\θ}}_i-{\mathrm{\θ}}_i^s\right)\\ V_{P2}=\underset{i=1}{\overset{n+m}{\mathrm{\Σ}}}{P}_i\left({\mathrm{\θ}}_i-{\mathrm{\θ}}_i^s\right)\\ V_{P3}=\underset{i=1}{\overset{n+m}{\mathrm{\Σ}}}\frac{Qi}{a{\left(E_i^s\right)}^a}\left({\left(E_i\right)}^a-{\left(E_i^s\right)}^a\right)\\ V_{P4}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{m}{\mathrm{\Σ}}}\frac{E_i{E}_j}{\sqrt{E_i^2-2E_i{E}_j\cos \left({\mathrm{\θ}}_i-{\mathrm{\θ}}_j\right)+E_j^2}}{I}_{ij}\\ -\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{m}{\mathrm{\Σ}}}\frac{E_i^s{E}_j^s}{\sqrt{{\left(E_i^s\right)}^2-2E_i^s{E}_j^s\cos \left({\mathrm{\θ}}_i^s-{\mathrm{\θ}}_j^s\right)+{\left(E_j^s\right)}^2}}{I}_{ij}\\ -\underset{i=1}{\overset{n+m-1}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n+m}{\mathrm{\Σ}}}{B}_{ij}\left(E_i{E}_j\cos \left({\mathrm{\θ}}_i-{\mathrm{\θ}}_j\right)-E_i^s{E}_j^s\cos \left({\mathrm{\θ}}_i^s-{\mathrm{\θ}}_j^s\right)\right)\\ -\frac{1}{2}\underset{i=1}{\overset{n+m}{\mathrm{\Σ}}}{B}_{ii}\left(E_i^2-{\left(E_i^s\right)}^2\right)\\ V_{damping}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\∫}_{t_s}^t{D}_i{\mathrm{\ω}}_i\frac{{\mathrm{d\θ}}_i}{dt}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\∫}_{{\mathrm{\θ}}_{i,s}}^{{\mathrm{\θ}}_i}{D}_i{\mathrm{\ω}}_i{\mathrm{d\θ}}_i\end{array}\right.---\left(4\right)$

Wherein, V is described energy function, V_KFor the kinetic energy of virtual synchronous electromotor, V_PFor the potential energy that system is total, V_P1For all Virtual synchronous generator mechanical power inputs the rotor potential energy caused, V_P2The potential energy caused for whole active loads, V_P3For all The potential energy that reactive load causes, V_P4For storage and the potential energy in network, behalf stable equilibrium point, i, j are the index value of node, n For node number, D_iIt is the damping of i-th electromotor, B_ij、B_iiIt is respectively the transadmittance between node i, j and the self-admittance of node i, E_i、E_jIt is the voltage at node i, j respectively,It is the angular frequency of i-th electromotor, θ_iIt is that i-th electromotor is relative in inertia The rotor angle of the heart, wherein a is that constant value is generally 2.

Method the most according to claim 3, it is characterised in that the detailed process of described step 3 is:

Step 3.1, obtains its pinch system such as following formula to formula (1)

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\θ}}}_i=P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}=f_i\left(\mathrm{\θ}\right)\\ {\mathrm{\θ}}_m=\frac{\underset{i=1}{\overset{m-1}{\mathrm{\Σ}}}{M}_i{\mathrm{\θ}}_i}{M_m}\end{array}\right.---\left(5\right)$

Step 3.2, uses path during faultAsk for exit point θ^EP, θ^EPIt is that to there is pinch system steady for projection path Fixed more borderline, it is specially

θ is obtained by the detail mathematic model of power system^EPPath during fault, formula (1) after the θ value obtained brings fault into, power is inclined Residual quantity equation, detects exit point θ^EPIt is to be arrived first local potential energy maximum by projection path；Wherein departure equation is

$P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}=f_i\left(\mathrm{\θ}\right)---\left(6\right)$

When meeting condition f_i* θ is obtained behind d θ=0^EP；

Step 3.3, with exit point θ^EPFor initial point, being integrated the pinch system of formula (6) gained, the curve along integration goes First minima shown in seeking (7)

$F\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{m}{\Σ}}{f_i}^2\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{m}{\Σ}}{\[P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}\]}^2---\left(7\right)$

First minima obtained is minimal gradient point θ^MGP。

Step 3.4, with minimal gradient point θ^MGPFor initial value, by merit after Newton-Raphson method iterative (m-1) individual fault Rate departure equation

$\left\{\begin{array}{c}{f}_i\left(\mathrm{\θ}\right)=P_{M_i}-P_{emi}-\frac{M_i}{M_T}{P}_{COI}=0\\ {\mathrm{\θ}}_m=\frac{\underset{i=1}{\overset{m-1}{\Σ}}{M}_i{\mathrm{\θ}}_i}{M_m}\end{array}\right.---\left(8\right)$

Obtain the controlling unstable equilibrium point CUEP of pinch system；

Step 3.4, substitutes into formula (2) by state parameter θ, ω of CUEP and obtains critical energy value V_cr。

Method the most according to claim 4, it is characterised in that in step 4, if the critical energy value V that step 3.4 obtains_cr Transient state energy value equal to formula (2) gained, it is thus achieved that corresponding θ, ω, find in time-domain-simulation obtained θ, ω corresponding time Between value for critical clearing time.