CN104836225A - Power system transient stability determination method based on first-order dimensionality reduction phase path - Google Patents

Abstract

The invention discloses a power system transient stability determination method based on first-order dimensionality reduction phase path, and relates to the safe stable operation technology of the power system. The method is based on the WAMS system, and the dimensionality reduction conversion of the power angle track in the n-dimension power angle space, and the operation track of the multi-computer system can be mapped to be the motion in one point, and the new phase plane can be formed by the motion parameters after the mapping conversion, and then the power system transient stability determination can be carried out by observing the geometrical characteristics of the tracks in the phase plane, and then the fast and accurate transient stability of the power system can be realized. The power system transient stability determination method is advantageous in that the useless elements having the unstability of the power angle in the generators can be removed by the dimensionality reduction conversion, and the grouping is not required, and therefore the wrong determination caused by the grouping can be prevented, and the transient stability of the system can be determined quickly and effectively.

Description

A kind of power system transient stability method of discrimination based on single order dimensionality reduction phase path

Technical field

The present invention relates to the safe and stable operation technology of electric power system, relate more specifically to a kind of power system transient stability based on single order dimensionality reduction phase path and differentiate.

Background technology

Operationally often can there is different faults or disturbance in electric power system, as short trouble, disconnection fault etc., when fault occurs, the automatic control equipment in electric power system can perceive the generation of fault, and automatically excises fault, by Fault Isolation; At failture evacuation or after disappearing, take as the measure such as reclosing, auto-parallel restores electricity; By existing supervisory control and data acquisition system SCADA, dispatcher can instruct the ruuning situation of whole electric power system, comprises the relevant information such as generation and disappearance of fault.

Electric power system transient stability refers to that electric power system is when normally running, after being subject to a large disturbance, can from original running status, and be not transitioned into lock-out a new running status, and stably run under new running status.The transient stability of electric power system is the significant problem that researcher pays close attention to always, the behavior of electric power system can characterize with a series of differential-algebraic equation, in theory, only require and solve these equations, the track of system state amount for the time can be obtained, thus judge the stability of this system; But it is very difficult for solving such equation, along with the complexity all the more of system, scale is huge all the more, and the difficulty of this task increases especially widely.

At present, conventional transient stability analysis method has time domain simulation method, based on the direct method of lyapunov energy function, take Extended Equal Area Criterion as the system equivalent method of representative, and take intelligent algorithm as the steady analytical method temporarily etc. instructed, but these methods all also exist respective shortcoming, as speed is comparatively slow, large etc. to dominant pattern dependence.

Along with WAMS (WideArea Measurement System, the WAMS) popularization and application in electric power system, occur much studying based on the transient stability analysis of power system of WAMS.Sun Wen etc. are at periodical " electric power network technique " 2009,33 (14): 16-20 " the electric power system on-line transient stability analytical methods " delivered are by being mapped as the simple motion on 1 dimension coordinate axle by the rocking process of multimachine system, and a kind of stability criterion of practicality is given according to the amount of exercise situation of change after mapping, but this method needs to observe a and v two curves simultaneously, it is very inconvenient to use; Xie Huan etc. are at periodical " Proceedings of the CSEE " 2006, " the power system transient stability identification based on phase path concavity and convexity " delivered in 26 (05): 38-42 proposes the new method that a kind of phase plane utilizing generator's power and angle δ and angular velocity omega to form carrys out analytical system transient stability, the phasor concavity and convexity transient stability evaluation demonstrating two-dimentional system is theoretical, and give new criterion, but this method is very large to the dependence of dominant pattern.And in transient process, leading Failure Model probably changes, at this moment erroneous judgement may be caused.

Summary of the invention

For prior art Problems existing, the invention provides a kind of power system transient stability method of discrimination based on single order dimensionality reduction phase path.The method is based on WAMS system, dimensionality reduction conversion is carried out to the phase-swing curves tieed up in space, merit angle at n, the running orbit of multimachine system is mapped as the motion of a point, the kinematic parameter after mapping transformation is utilized to form new phase plane, carried out the differentiation of electric power system transient stability by the geometric properties observing track in phase plane, achieve power system transient stability fast and accurately and differentiate.Advantage is: dimensionality reduction conversion eliminates in each generator and do not have contributive composition to merit angle unstability, and does not need to hive off, and avoids the erroneous judgement that mistake of hiving off causes, can the transient stability of fast and effeciently judgement system.

For achieving the above object, the present invention adopts following technical scheme:

Based on a power system transient stability method of discrimination for single order dimensionality reduction phase path, the method comprises the following step;

A, electric power system fault are at t=t _jmoment disappears, and gets the merit angle δ of i-th generator in the electric power system of t WAMS system log (SYSLOG) _iand angular velocity omega _i, i=1,2 ..., n, n be the number of units of all generators in electric power system;

B, by t merit angle δ _iand angular velocity omega _ibe transformed to the kinematic parameter relative to system inertia center COI:

$\left\{\begin{array}{c}{\mathrm{\θ}}_i={\mathrm{\δ}}_i-{\mathrm{\δ}}_{\mathrm{COI}}\\ {\widetilde{\mathrm{\ω}}}_i={\mathrm{\ω}}_i-{\mathrm{\ω}}_{\mathrm{COI}}\end{array}\right.i=\mathrm{1,2},\·\·\·,n$

Wherein, ${\mathrm{\δ}}_{\mathrm{COI}}=\frac{1}{M_T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{\mathrm{\δ}}_i,{\mathrm{\ω}}_{\mathrm{COU}}=\frac{1}{M_T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{\mathrm{\ω}}_i,M_T=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i,$ M _ibe the inertia time constant of i-th generator, θ _ibe the merit angle of i-th generator relative to COI, be the angular speed of i-th generator relative to COI;

C, calculate t through one dimension angle R, the one dimension angular speed v of single order dimensionality reduction and one dimension angular acceleration a:

$\left\{\begin{array}{c}R=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\θ}}_i^2}\\ v=\frac{\mathrm{dR}}{\mathrm{dt}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left({\mathrm{\θ}}_i{\widetilde{\mathrm{\ω}}}_i\right)}{R}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\θ}}_i}{R}{\widetilde{\mathrm{\ω}}}_i\\ a=\frac{\mathrm{dv}}{\mathrm{dt}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\widetilde{\mathrm{\ω}}}_i^2+{\mathrm{\θ}}_i\frac{d{\widetilde{\mathrm{\ω}}}_i}{\mathrm{dt}})-v^2}{R}\end{array}\right.$

D, when v is less than zero, make t=t+ Δ t, Δ t is the sampling interval duration of WAMS, returns step a; Otherwise, carry out step e;

E, calculating $\frac{d^{2}v}{{\mathrm{dR}}^2}=\frac{\∂\left(\frac{a}{v}\right)}{\∂a}\·\frac{\mathrm{da}}{\mathrm{dR}}+\frac{\∂\left(\frac{a}{v}\right)}{\∂v}\·\frac{\mathrm{dv}}{\mathrm{dR}}=\frac{v\frac{\mathrm{da}}{\mathrm{dt}}-a^2}{v^3}$

When make t=t+ Δ t, return step a; Otherwise judgement system unstability, and record current time t.

Compared with prior art, the present invention has the following advantages and beneficial effect:

(1) dimensionality reduction mapping method of the present invention not only remains in system the useful information showing merit angle unstability, but also eliminates in each generator angular speed and do not have contributive composition to merit angle unstability.

(2) the present invention proposes and utilize the concavity and convexity of track to judge the method for power system transient stability in R-v phase plane, can the transient stability of fast and effeciently judgement system.

(3) method that the present invention proposes does not need to hive off, and avoids the erroneous judgement that mistake of hiving off causes, and judges that the time of unstability is faster than general needs equivalent method of hiving off.

Accompanying drawing explanation

Fig. 1 is flow chart of the present invention.

Fig. 2 is the movement locus of system in space, two-dimentional merit angle.

Fig. 3 is stable R-v phase path.

Fig. 4 is unstable R-v phase plane.

Fig. 5 is New England 39 node system.

The R-v phase path of system stability when Fig. 6 is circuit 1-39 fault.

The power-angle curve of system stability when Fig. 7 is circuit 1-39 fault.

When Fig. 8 is circuit 1-39 fault mistake hive off under δ-ω phase path.

The R-v phase path of system instability when Fig. 9 is circuit 26-29 fault.

The power-angle curve of system instability when Figure 10 is circuit 26-29 fault.

When Figure 11 is circuit 26-29 fault, unit 39 is the δ-ω phase path of delayed group.

When Figure 12 is circuit 26-29 fault, unit 30-37 is the δ-ω phase path of advanced group.

Embodiment

Below in conjunction with accompanying drawing, the present invention is further illustrated.

Embodiment one

Based on a power system transient stability method of discrimination for single order dimensionality reduction phase path, as shown in Figure 1, the method comprises the following step;

A, electric power system fault are at t=t _jmoment disappears, and gets the merit angle δ of i-th generator in the electric power system of t WAMS system log (SYSLOG) _iand angular velocity omega _i, i=1,2 ..., n, n be the number of units of all generators in electric power system;

B, by t merit angle δ _iand angular velocity omega _ibe transformed to the kinematic parameter relative to system inertia center COI: $\left\{\begin{array}{c}{\mathrm{\θ}}_i={\mathrm{\δ}}_i-{\mathrm{\δ}}_{\mathrm{COI}}\\ {\widetilde{\mathrm{\ω}}}_i={\mathrm{\ω}}_i-{\mathrm{\ω}}_{\mathrm{COI}}\end{array}\right.i=\mathrm{1,2},\·\·\·,n\backslash *\mathrm{MERGEFORMAT}---\left(1\right)$

Now, the equation of rotor motion of system becomes

$\left\{\begin{array}{c}\frac{d{\mathrm{\θ}}_i}{\mathrm{dt}}={\widetilde{\mathrm{\ω}}}_i\\ \frac{d{\widetilde{\mathrm{\ω}}}_i}{\mathrm{dt}}=\frac{{\mathrm{\ω}}_N}{M_i}(P_{\mathrm{mi}}-P_{\mathrm{ei}})-\frac{{\mathrm{\ω}}_N}{M_T}\end{array}\right.i=\mathrm{1,2},\·\·\·,n\backslash *\mathrm{MERGERORMAT}---\left(2\right)$

In formula, θ _i=δ _i-δ _cOIbe the merit angle of i-th generator relative to COI, be the angular speed of i-th generator relative to COI, ${\mathrm{\ω}}_{\mathrm{COI}}=\frac{1}{M_T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{\mathrm{\ω}}_i;P_{\mathrm{COI}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}).$

C, calculate t through one dimension angle R, the one dimension angular speed v of single order dimensionality reduction and one dimension angular acceleration a:

$\left\{\begin{array}{c}R=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\θ}}_i^2}\\ v=\frac{\mathrm{dR}}{\mathrm{dt}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left({\mathrm{\θ}}_i{\widetilde{\mathrm{\ω}}}_i\right)}{R}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\θ}}_i}{R}{\widetilde{\mathrm{\ω}}}_i\\ a=\frac{\mathrm{dv}}{\mathrm{dt}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\widetilde{\mathrm{\ω}}}_i^2+{\mathrm{\θ}}_i\frac{d{\widetilde{\mathrm{\ω}}}_i}{\mathrm{dt}})-v^2}{R}\end{array}\right.\backslash *\mathrm{MERGERORMAT}---\left(3\right)$

The performance of 3.1 system merit angle unstabilitys

R in formula characterizes the degree that the angle swing of original system merit is opened, if formula is set up, then and original system merit angle unstability; If original system merit angle unstability, then formula is also had to set up.

$R\left(t\right){|}_{t\→\∞}\→\∞\backslash *\mathrm{MERGEFORMAT}---\left(4\right)$

The relation of parameter and original system parameter after 3.2 mappings

With inertia center for initial point, θ _ifor reference axis and the n that launches ties up in space, merit angle, each generator's power and angle θ _ichange can be expressed as a track in this space.And the R representative of formula definition is the Euclidean distance at point on a certain moment track and inertia center.Accompanying drawing 2, to represent 2 dimension spaces, merit angle of two machine systems, analyzes R, v and θ _i, relation.

In accompanying drawing 2, at a time t=t _a, system motion has arrived A point, now the coordinate of A point represent generator 1 and generator 2 merit angle now respectively, then size be one dimension angle R now ^a.After a very little time period Δ t, along with the change at merit angle, system motion is to B point, same, and the coordinate of B point represents generator merit angle now, size be R ^b.

And in process from point A to point B, can be seen by accompanying drawing 2 do obviously as can be seen from formula, one dimension angular speed v be by weighted sum and obtaining, and weight is θ _i/ R, can be regarded as the θ of i-th generator _ito " percentage contribution " of R; In the some A of accompanying drawing 2, then have according to formula

$v^{\mathrm{AB}}=\cos {\mathrm{\α}}_1^A{\widetilde{\mathrm{\ω}}}_1^{\mathrm{AB}}+\cos {\mathrm{\α}}_2^A{\widetilde{\mathrm{\ω}}}_2^{\mathrm{AB}}\backslash *\mathrm{MERGEFORMAT}---\left(5\right)$

This illustrates, v ^aBbe with about α ^aone projection synthesis.Tie up in space, merit angle at n, also have similar projection relation, such as formula shown.

$v=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\cos {\mathrm{\α}}_i{\widetilde{\mathrm{\ω}}}_i\backslash *\mathrm{MERGEFORMAT}---\left(6\right)$

In formula, α _ifor R and θ _ithe space angle of reference axis, α _ibecome in time.

This relation can be explained like this: in process from point A to point B, mathematically can think track first from A point to D point, and then to B point; Due to namely do not have influence on the stability of system from A point to D point, and ensuing D → B just the real increase to R contribute to some extent.That is to say, with in containing not having contributive composition to merit angle unstability, only select and assign to synthesize v to the merit angle contributive one-tenth of unstability, Here it is weight cos α _i=θ _ithe meaning of/R.

As a same reason, when track is from B point to C point, also can analyze by similar method.Can find from the analysis to accompanying drawing 2, the transfer pair system of through type is carried out single order dimensionality reduction and is not only remained in system the useful information showing merit angle unstability, but also filtering does not have contributive composition to merit angle unstability.

D, when v is less than zero, make t=t+ Δ t, Δ t is the sampling interval of WAMS, returns step a; Otherwise, carry out step e;

E, calculating $\frac{d^{2}v}{{\mathrm{dR}}^2}=\frac{\∂\left(\frac{a}{v}\right)}{\∂a}\·\frac{\mathrm{da}}{\mathrm{dR}}+\frac{\∂\left(\frac{a}{v}\right)}{\∂v}\·\frac{\mathrm{dv}}{\mathrm{dR}}=\frac{v\frac{\mathrm{da}}{\mathrm{dt}}-a^2}{v^3}$

When make t=t+ Δ t, return step a; Otherwise judgement system unstability, and record current time t.

The phase plane of system after 5.1 single order dimensionality reductions

Using R as abscissa, v, as ordinate, gets final product the R-v phase plane of construction system.By formula, because R >=0 perseverance is set up, thus the R-v track of system be only positioned at phase plane first, fourth quadrant or tangent with v axle, there will not be second and third quadrant; And due to fourth quadrant v < 0, therefore track only can in first quartile unstability, if phase path passes through R axle, owing to having v=0 on R axle, so track is inherently orthogonal with R axle when passing through.

Accompanying drawing 3 is the R-v phase plane trajectory of certain systems stabilisation, and 4, accompanying drawing represents the R-v phase plane trajectory of certain unstability system.

Accompanying drawing 3 and accompanying drawing 4 show, if system stability, then track can penetrate fourth quadrant from first quartile, enter into backswing state (v<0), and now v starts oppositely to increase and R constantly reduces; And if system is unstable; then track will start to disperse at certain point of first quartile; due to generally all certain margin of safety can be had when electric power system normally runs; so generally all certain decelerating phase can be experienced after Failure elimination; and system unstability is not enough just because of slowing down power(SDP); cause system again to enter boost phase, namely track there will be one " lower salient point ", and this indicates the unstability of system.

The mathematical feature of 5.2 unstability tracks

In R-v phase plane, the first derivative of track and second dervative are

$\frac{\mathrm{dv}}{\mathrm{dR}}=\frac{\mathrm{dv}}{\mathrm{vdt}}=\frac{a}{v}\backslash *\mathrm{MERGEFORMAT}---\left(7\right)$

$\frac{d^{2}v}{{\mathrm{dR}}^2}=\frac{\&PartialD;\left(\frac{a}{v}\right)}{\&PartialD;a}\&CenterDot;\frac{\mathrm{da}}{\mathrm{dR}}+\frac{\&PartialD;\left(\frac{a}{v}\right)}{\&PartialD;v}\&CenterDot;\frac{\mathrm{dv}}{\mathrm{dR}}=\frac{v\frac{\mathrm{da}}{\mathrm{dt}}-a^2}{v^3}\backslash *\mathrm{MERGEFORMAT}---\left(8\right)$

Mathematically known, curve occurs that flex point illustrates that the second dervative of this curve equals zero, namely when time, curve is in this interval epirelief (concave function); When time, curve convex under this interval (convex function).

Can infer from accompanying drawing 3 and accompanying drawing 4, if phase path has occurred becoming convex function from concave function at first quartile, then system likely can unstability.If phase path enters fourth quadrant from first quartile through R axle, then system enters the decelerating phase, and known is stable at this pendulum subsystem.

5.3 based on the electric power system transient stability criterion of R-v phase plane

Real-time Transient Stability Criterion based on R-v phase plane is: after failure vanishes, if system is concave function at the track of R-v phase plane first quartile always, namely has the track of one quadrant changes convex function into by concave function, namely by become then system unstability.

Embodiment two

Based on a power system transient stability method of discrimination for single order dimensionality reduction phase path, for IEEE10 machine 39 node modular system, system as shown in Figure 5.

Utilize the method that the present invention proposes, the merit angle of each generator after register system fault and angular speed, through type and formula carry out the conversion of single order dimensionality reduction, form R-v phase plane, calculate the second dervative of track in phase plane, by judging its symbol to differentiate the stability of now system.

Fault one: circuit 1-39, near bus 1 place, three-phase shortcircuit occurs when t=0, failure vanishes after 0.1s.

The method utilizing the present invention to propose is analyzed system now, and the R-v phase-plane diagram of system as shown in Figure 6.

As can be seen from accompanying drawing 6, after failure vanishes, the phase path of system is in recessed state always, and passes perpendicularly through R axle, enters reversal phase, and criterion shows now system at first quartile perseverance is less than zero, now can judge that system is stable.What accompanying drawing 7 showed is the power-angle curve of each generator, shows that now system is stable equally.

As a comparison, the concavity and convexity of track in equivalent δ-ω phase plane is adopted to carry out the method for judgement of stability, but now to hive off to this system, be divided into advanced group S and delayed group A two groups respectively, carry out merit angle and angular speed equivalence again, finally utilize the concavity and convexity of phase path to carry out the judgement of transient stability.If now unit 30-32 is divided into advanced group, all the other are delayed group, and δ-ω phase path equivalent so as shown in Figure 8.

δ-ω phase path in accompanying drawing 8 has occurred by the convex situation of concave change when fourth quadrant, and when two-shipper equivalence, likely there will be reverse unstability, if when adopting concavity and convexity to judge, this situation can be mistaken for reverse unstability, and in fact system is stable, hiving off of the mistake that has its source in judged by accident in this example.Adopt method of the present invention then can avoid hiving off, thus avoid erroneous judgement.

Fault two: near bus 26 place, three-phase shortcircuit occurs at 0 moment circuit 26-29, failure vanishes after 0.135s.

The method utilizing the present invention to propose is analyzed system now, and the R-v phase-plane diagram of system as shown in Figure 9.

In accompanying drawing 9, R-v phase plane has occurred by recessed turn of convex situation at first quartile, illustrates that now system is unstable, according to carried criterion, occur by become moment be approximately 0.32s, namely after failure vanishes, 0.185s judges system unstability.

The power-angle curve of each generator as shown in Figure 10.

Accompanying drawing 10 shows now system merit angle unstability, if judge system unstability according to engineering experience, namely maximum work angular difference exceedes certain larger value, and continues the regular hour, and threshold value is taken as 180 ° here, then approximately will arrive 1.35s and just can judge system unstability.

Sentence surely according to the equivalent δ-ω phase path based on two groups, first hive off, accompanying drawing 11 gives using unit 39 as delayed group, δ-ω phase path when all the other units are leading group.Accompanying drawing 11 has also occurred by recessed turn convex, and describe the unstability of its system that also can correctly judge, it occurs that the time of flex point is approximately 0.55s, although faster than engineering experience, but slower than method herein.

If now unit 30-37 is divided into advanced group, unit 38-39 is delayed group, and accompanying drawing 12 gives δ-ω phase path equivalent in this case.Can find, the track concavity and convexity with reference to the accompanying drawings in 12 cannot correctly judge system unstability, and this is that mistake is hived off and caused.

Except being presented above two fault examples, table 1 also to list under other faults method herein and equivalent δ-ω phase plane method to the recognition time of system instability.The time showing method identification herein in table wants Zao than δ-ω phase path, this is because this method remains by mapping the composition contributed to some extent system unstability, eliminate and influential part is not had to unstability, so unstability shows more early, earlier can identify the feature of unstability.

The comparison of the lower three kinds of method recognition times of table 1 different unstability fault

Claims (1)

1., based on a power system transient stability method of discrimination for single order dimensionality reduction phase path, it is characterized in that, the method comprises the following step;

A, electric power system fault are at t=t _jmoment disappears, and gets the merit angle δ of i-th generator in the electric power system of t WAMS system log (SYSLOG) _iand angular velocity omega _i, i=1,2 ..., n, n be the number of units of all generators in electric power system;

B, by t merit angle δ _iand angular velocity omega _ibe transformed to the kinematic parameter relative to system inertia center COI:

$\left\{\begin{array}{c}{\mathrm{\&theta;}}_i={\mathrm{\&delta;}}_i-{\mathrm{\&delta;}}_{\mathrm{COI}}\\ {\widetilde{\mathrm{\&omega;}}}_i={\mathrm{\&omega;}}_i-{\mathrm{\&omega;}}_{\mathrm{COI}}\end{array}\right.I=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n$

Wherein, ${\mathrm{\&delta;}}_{\mathrm{COI}}=\frac{1}{M_T}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{\mathrm{\&delta;}}_i,{\mathrm{\&omega;}}_{\mathrm{COI}}=\frac{1}{M_T}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{\mathrm{\&omega;}}_i,M_T=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i,$ M _ibe the inertia time constant of i-th generator, θ _ibe the merit angle of i-th generator relative to COI, be the angular speed of i-th generator relative to COI;

C, calculate t through one dimension angle R, the one dimension angular speed v of single order dimensionality reduction and one dimension angular acceleration a:

$\left\{\begin{array}{c}R=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_i^2}\\ v=\frac{\mathrm{dR}}{\mathrm{dt}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left({\mathrm{\&theta;}}_i{\widetilde{\mathrm{\&omega;}}}_i\right)}{R}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&theta;}}_i}{R}{\widetilde{\mathrm{\&omega;}}}_i\\ a=\frac{\mathrm{dv}}{\mathrm{dt}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\widetilde{\mathrm{\&omega;}}}_i^2+{\mathrm{\&theta;}}_i\frac{d{\widetilde{\mathrm{\&omega;}}}_i}{\mathrm{dt}})-v^2}{R}\end{array}\right.$

D, when v is less than zero, make t=t+ Δ t, Δ t is the sampling interval duration of WAMS, returns step a; Otherwise, carry out step e;

E, calculating $\frac{d^{2}v}{{\mathrm{dR}}^2}=\frac{\&PartialD;\left(\frac{a}{v}\right)}{\&PartialD;a}\&CenterDot;\frac{\mathrm{da}}{\mathrm{dR}}+\frac{\&PartialD;\left(\frac{a}{v}\right)}{\&PartialD;v}\&CenterDot;\frac{\mathrm{dv}}{\mathrm{dR}}=\frac{v\frac{\mathrm{da}}{\mathrm{dt}}-a^2}{v^3}$

When time, make t=t+ Δ t, return step a; Otherwise judgement system unstability, and record current time t.