CN105468909A - Time delay power system electromechanical oscillation mode computing method based on SOD-PS-R R algorithm - Google Patents

Abstract

The invention discloses a time delay power system electromechanical oscillation mode computing method based on the SOD-PS-R (Solution Operator Discretization-Pesudo Spectrum and Rotation) algorithm. According to the relation between an eigenvalue of a time delay power system model and an operator eigenvalue of the time delay power system model, computing the eigenvalue of the time delay power system model is converted into computing the operator eigenvalue, so that a problem about computing an electromechanical oscillation mode of a time delay power system is converted into a problem about computing the operator eigenvalue; the method is used for the electromechanical oscillation mode of a large-scale time delay power system. The SOD-PS-R requires computing the maximal setting number eigenvalue of an operator discretization approximate matrix norm only, therefore one computing is enough for obtaining the electromechanical oscillation mode of the large-scale time delay power system.

Description

Based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm

Technical field

The present invention relates to the time-lag power system electromechanic oscillation mode computing method based on SOD-PS-R (SolutionOperatorDiscretization-PesudoSpectrumandRotation, Solution operator puppet spectrum discretize and rotation) algorithm.

Background technology

Along with the rise of global energy internet, the scale of interconnected electric power system increases gradually, and interval low-frequency oscillation problem is more remarkable.Traditional solution installs power system stabilizer, PSS (PowerSystemStabilizer, PSS), but derive from locality due to its feedback control signal, can not the inter-area oscillations of effective damping interconnected electric power system.

The appearance of WAMS (Wide-AreaMeasurementSystem, WAMS) brings new opportunity to the development of extensive interconnected electric power system stability analysis and control.Based on the interconnected network low-frequency oscillation control of the Wide-area Measurement Information that WAMS provides, by introducing the wide area feedback signal of effectively reflection inter-area oscillation mode, good damping control performance can be obtained, it is for solving the inter-area low-frequency oscillation problem in interconnected network, and then improve the ability to transmit electricity of system and provide new control device, have well and application prospect widely.

When wide area signal transmits and processes in the WAMS communication network be made up of different communication medium (as optical fiber, telephone wire, digital microwave, satellite etc.), there is the communication delay changed between tens to hundreds of millisecond.Time lag be cause Systematical control to be restrained losing efficacy, operation conditions worsens and a kind of major incentive of system unstability.Therefore, when utilizing wide area measurement information to carry out electric system closed-loop control, the impact of time lag must be taken into account.

In modern power systems, what small interference stability was mainly paid close attention to is electromechanical oscillations problem.Eigenvalue Method based on state-space model is the powerful tool of research aircraft electric oscillation.At present, researchist has proposed the effective ways of many calculating large-scale electrical power system Critical eigenvalues, mainly comprise the preference pattern analytic approach based on reduced order system, AESOPS algorithm and S-matrix method, the subspace iteration methods such as the sequential methods such as power method, Newton method, Rayleigh Rayleigh quotient iteration and simultaneous iterative, Arnoldi algorithm, the double iterative again decomposed and Jacobi-Davidson method.These methods all make use of the openness of augmented state matrix when calculating section eigenwert, majority method is all by carrying out spectral transformation to original system thus the distribution of change characteristic spectrum, then ask for the eigenwert of system, then obtain the critical eigenvalue of original system by inverse transformation.But above-mentioned method does not all consider the impact of time lag.Chinese invention patent calculates and Convenient stable criterion .201210271783.8:[P based on the time-lag power system eigenwert that Pad é is approximate]. utilize Pade approximation polynomial to approach Time Delay, and then the critical eigenvalue of the computing system rightmost side, and judge the time lag stability of system.Chinese invention patent is based on the extensive time-lag power system eigenwert computing method .201510055743.3. China of EIGD, 201510055743.3 [P]. propose a kind of extensive time-lag power system eigenwert based on display IGD (ExplicitIGD, EIGD) and calculate.Utilize the critical eigenvalue of the system rightmost side calculated, the stability of system under fixed time lag can be judged.These time lag Convenient stable criterion, all need the critical eigenvalue of interior by Multiple-Scan [0.1,2.5] Hz low-frequency oscillation frequency range, the close imaginary axis, could judge the time lag stability of system.

Summary of the invention

Object of the present invention is exactly to solve the problem, and provides a kind of computing method of the time-lag power system electromechanic oscillation mode character pair value based on SOD-PS-R algorithm, in order to obtain the electromechanic oscillation mode of extensive time-lag power system.SOD-PS-R algorithm only need calculate the maximum setting of the modulus value of Solution operator discretize approximate matrix several eigenwert, by once calculating, just can obtain the electromechanic oscillation mode of extensive time-lag power system.

To achieve these goals, the present invention adopts following technical scheme:

Communication delay is introduced in the modeling process of actual large-scale electrical power system.Time-lag power system comprises without time-lag power system, wide area Feedback Delays, wide area damping control and wide area output time lag four part.And time-lag power system model is carried out linearization near its steady-state operation point, obtain the inearized model for time-lag power system Small signal stability analysis.

By spectral mapping theorem, described time-lag power system is positioned near the imaginary axis, the eigenvalue λ of frequency within the scope of [0,2] Hz, is mapped as the eigenwert of Solution operator T (h), and is distributed in thick and fast near unit circle.

Solution operator T (h) of described time-lag power system is carried out puppet spectrum discretize, by calculating finite dimension discretize matrix T _m,Nsetting several eigenwert that modulus value is maximum, can be converted into the critical eigenvalue that time-lag power system electromechanic oscillation mode is corresponding.

In order to improve described time-lag power system eigenwert convergence, need to increase T _m,Neigenwert between relative distance.The T that the critical eigenvalue λ-conversion utilizing rotation of coordinate preprocess method the damping ratio of time-lag power system can be less than given constant ζ is correspondence _m,Nthe eigenwert that modulus value is greater than 1.

Discretize due to described time-lag power system Solution operator T (h) is counted and is necessary for integer, requires that the maximum time lag after rotation of coordinate is real number.When the anglec of rotation is not 0, obtain postrotational approximation characteristic equation and corresponding Solution operator discretize approximate matrix T " _m,N.

For described time-lag power system discretize approximate matrix T " _m,Nexponent number be (QM+1) n, when system scale is larger, need to adopt sparse features value-based algorithm, as implicit restarted Arnoldi algorithm, calculate the maximum setting of its modulus value several eigenwert μ ".Calculating μ " afterwards, first obtaining the eigenvalue λ of corresponding time-lag power system after rotation of coordinate with it according to spectrum mapping relations ".Then, the eigenvalue λ " e of time-lag power system is obtained by coordinate despining ^{-j θ}.

Based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, comprise the steps:

S1: set up time-lag power system model; Relation between the eigenwert of foundation time-lag power system model and time-lag power system solution to model operator eigenwert, changes into the eigenwert calculating Solution operator by the eigenwert calculating time-lag power system model; Thus the problem calculating time-lag power system electromechanic oscillation mode is converted into the eigenvalue problem calculating Solution operator;

S2: carry out discretize to Solution operator by pseudo-spectral method, obtains the discretize matrix of Solution operator;

S3: utilize rotation of coordinate preprocess method, the discretize matrix rotation of Solution operator is carried out approximate obtaining postrotational Solution operator discretize approximate matrix, the damping ratio of time-lag power system model is less than the critical eigenvalue λ of setting constant ζ, is transformed to the eigenwert μ that corresponding postrotational Solution operator discretize approximate matrix modulus value is greater than 1 ";

S4: adopt sequential method or subspace method, from the eigenwert μ of step S3 " in find the modulus value of Solution operator discretize approximate matrix maximum setting several eigenwert μ ";

S5: obtain eigenwert μ from step S4 " afterwards, successively by μ ", through spectrum mapping, coordinate despining and newton's checking treatment, finally obtains the eigenvalue λ of time-lag power system model, the electromechanic oscillation mode of eigenvalue λ and corresponding time-lag power system.

Time-lag power system model in described step S1 is as follows:

In formula: for the state variable vector of electric system, n is system state variables sum.T is current time.0< τ ₁< τ ₂< ... < τ _i< τ _mfor the time lag constant of Time Delay, wherein maximum time lag is τ _m. for system state matrix, it is dense matrix. for system time lags state matrix, it is sparse matrix.Δ x (t) is the increment of t system state variables, Δ x (t-τ _i) be t-τ _ithe increment of moment system state variables, for the increment of t system state variables derivative.The initial value (i.e. starting condition) that Δ x (0) is system state variables, and be abbreviated as

The secular equation of the linearized system that formula (1.1) represents is:

$({\tilde{A}}_0+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{e}^{-{\mathrm{\&lambda;\&tau;}}_i})v=\mathrm{\&lambda;}v---\left(1.2\right)$

In formula: λ is eigenwert, v is right proper vector corresponding to eigenwert.

In described step S2, Solution operator T (h): X → X is defined as the initial condition in space X the linear operator of time-lag power system solution segmentation after transferring to the h moment; Wherein, h is transfer step-length, 0≤h≤τ _m;

Wherein, s is integration variable, and θ is variable, with be respectively the state of 0 and h+ θ moment time-lag power system.

Relation between the eigenwert of described time-lag power system model and time-lag power system solution to model operator eigenwert:

From spectral mapping theorem, between the eigenwert μ of Solution operator T (h) and the eigenvalue λ of time lag system, there is following relation:

$\mathrm{\&lambda;}=\frac{1}{t}ln\mathrm{\&mu;},\mathrm{\&mu;}\&Element;\mathrm{\&sigma;}\left(T\right(t\left)\right)\backslash \left\{0\right\}---\left(1.4\right)$

In formula: σ (T (h)) represents the spectrum of Solution operator, represent set difference computing.

The eigenwert that time lag system is positioned at left half complex plane is mapped within Solution operator unit circle, and the eigenwert that time lag system is positioned at right half complex plane is mapped as the eigenwert that Solution operator modulus value is greater than 1, and is positioned at outside unit circle.Therefore, utilize the eigenwert of Solution operator, just can judge the stability of former time lag system.If Solution operator at least exists the eigenwert that a modulus value is greater than 1, then can judge that former time lag system is unstable, if the modulus value of all eigenwerts of Solution operator is all less than 1, then former time lag system is asymptotically stability.

Solution operator T (h) is the Infinite Dimension Linear operator describing X → X mapping.In order to calculate the eigenwert of Solution operator and obtain former time-lag power system electromechanic oscillation mode characteristic of correspondence value according to this, first pseudo-spectral method (Pesudospectral is adopted, PS) discretize is carried out to Solution operator T (h), obtain corresponding with Solution operator, a finite dimensional approximate matrix, and then calculate the eigenwert of approximate matrix and obtain former time-lag power system electromechanic oscillation mode characteristic of correspondence value.

The step of described step S2 is as follows:

Discretize matrix T that is corresponding with Solution operator T (h), Standard basis form _m,Nbe expressed as follows:

$T_{M,N}=T_M^{\left(1\right)}+T_{M,N}^{\left(2\right)}{(I_{Nn}-U_N^{\left(2\right)})}^{-1}{U}_{M,N}^{\left(1\right)}---\left(1.5\right)$

In formula:

$T_M^{\left(1\right)}=\left[\begin{array}{ccc}{1}^{M\&times;1}& & \\ I^{(Q-2)M}& & \\ & T_M^Q& 0^{(M+1)\&times;M}\end{array}\right]\&CircleTimes;I_n={\tilde{T}}_M^{\left(1\right)}\&CircleTimes;I_n---\left(1.6\right)$

$T_{M,N}^{\left(2\right)}=\left[\begin{array}{c}{T}_{M,N}^1\\ 0^{QM\&times;N}\end{array}\right]\&CircleTimes;I_n={\tilde{T}}_{M,N}^{\left(2\right)}\&CircleTimes;I_n---\left(1.7\right)$

$U_{M,N}^{\left(1\right)}=\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i---\left(1.8\right)$

$U_N^{\left(2\right)}=\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i---\left(1.9\right)$

In formula (1.5) and formula (1.6), Q, M and N are given positive integer, 1 ^{m × 1}for the M dimensional vector that element is 1 entirely, I ^{(Q-2) M}for (Q-2) M rank unit matrix, 0 ^{(M+1) × M}for (M+1) × M rank null matrix, for Kronecker amasss computing.Matrix for height sparse matrix, and with system state matrix ${\tilde{A}}_i(i$ $=0,\mathrm{...},m)$ Irrelevant.

In formula (1.7), 0 ^{qM × N}for QM × N rank null matrix, matrix for height sparse matrix, and with system state matrix irrelevant.

In formula (1.8), its element is determined by Lagrange coefficient completely.

In formula (1.9), its element is determined by Lagrange coefficient completely.

In described step S3, utilize rotation of coordinate preprocess method, coordinate axis is rotated counterclockwise θ angle, the damping ratio of time-lag power system model is less than the critical eigenvalue λ of setting constant ζ (ζ=sin θ), is transformed to T _m,Nthe eigenwert μ that modulus value is greater than 1.It is the equal damping ratio of ζ that the postrotational imaginary axis correspond to damping ratio in former coordinate system.

If the eigenwert of time-lag power system model is λ after rotation of coordinate ^'so, the λ λ ' e in formula (1.2) ^{-j θ}replace, the secular equation after rotation of coordinate can be obtained:

$({\tilde{A}}_0^{\&prime;}+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{e}^{-{\mathrm{\&lambda;}}^{\&prime;}{\mathrm{\&tau;}}_i^{\&prime;}})v={\mathrm{\&lambda;}}^{\&prime;}v---\left(1.10\right)$

In formula:

${\tilde{A}}_0^{\&prime;}={\tilde{A}}_0{e}^{j\mathrm{\&theta;}},{\tilde{A}}_i^{\&prime;}={\tilde{A}}_i{e}^{j\mathrm{\&theta;}},i=1,\mathrm{...},m$

λ′＝λe ^jθ

τ′ _i＝τ _ie ^-jθ,i＝1,...,m

Due to h be real number, Q is integer, by relation known, the maximum time lag τ ' after rotation of coordinate _maxmust be real number.But, when θ ≠ 0, by τ ' _ithe known τ ' of expression formula _i(i=1 ..., m) must be plural number.Obviously, this and τ ' _maxfor real number is conflicting.Therefore, in the secular equation after rotation of coordinate, get:

τ′ _i＝τ _ie ^-jθ≈τ _i,i＝1,...,m

Thus, obtain approximation characteristic equation as follows:

$({\tilde{A}}_0^{\&prime;}+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{e}^{-{\mathrm{\&lambda;}}^{\&prime;\&prime;}{\mathrm{\&tau;}}_i})v={\mathrm{\&lambda;}}^{\&prime;\&prime;}v---\left(1.11\right)$

" be the approximate value of λ ' in formula: λ.

The approximate matrix T of discretize matrix after the Solution operator rotation that formula (1.11) is corresponding " _m,Ncan be expressed as:

$T_{M,N}^{\&prime;\&prime;}=T_M^{\left(1\right)}+T_{M,N}^{\left(2\right)}{(I_{Nn}-U_N^{\&prime;\&prime;})}^{-1}{U}_{M,N}^{\&prime;\&prime;}$

In formula:

$U_{M,N}^{\&prime;\&prime;}=\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}$

$U_N^{\&prime;\&prime;}=\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}$

Obviously, when anglec of rotation θ=0 °, T " _m,Njust deteriorate to T _m,N.

If T " _m,Neigenwert to be μ meet following relation between ", itself and λ ":

${\mathrm{\&lambda;}}^{\&prime;\&prime;}=\frac{1}{h}{\mathrm{ln\&mu;}}^{\&prime;\&prime;},{\mathrm{\&mu;}}^{\&prime;\&prime;}\&Element;\mathrm{\&sigma;}\left(T_{M,N}^{\&prime;\&prime;}\right)\backslash \left\{0\right\}.$

In described step S4, matrix T " _m,Nexponent number be (QM+1) n.For large-scale electrical power system, matrix T " _m,Nexponent number will be very huge.Therefore, in the discretize matrix T that application Solution operator is corresponding " _m,Nwhen solving the time lag eigenwert of large-scale electrical power system, iterative characteristic value-based algorithm (sequential method or subspace method) must be adopted to calculate the maximum setting of its modulus value several eigenwert.

The step of described step S4 is as follows:

Suppose, when the secondary iteration of kth, to need to calculate T " _m,Nwith vector product, w ^k=T " _m, concrete steps are as follows:

Step (41): by vector v ^kby row boil down to one matrix correspondingly, have: v ^k=vec (V ^k), wherein vec () is for being the computing of column vector by matrix compression.

Step (42): calculate

Step (43): calculate

Step (44): calculate $w^k=T_M^{\left(1\right)}{v}^k+T_{M,N}^{\left(2\right)}{q}^k.$

The step of described step (42) is as follows:

By U " _m,Nexpression formula substitutes into, and can obtain:

$\begin{array}{c}{p}^k-U_{M,N}^{\&prime;\&prime;}{v}^k\\ =(\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;})vec\left(V^k\right)\\ =vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{V}^k{E_i}^T\right)\\ =e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{V}^k{E_i}^T\right)\\ =e^{j\mathrm{\&theta;}}vec(\underset{i=0}{\overset{m}{\&Sigma;}}\&lsqb;{\tilde{A}}_i{v}_1^k,\mathrm{...},{\tilde{A}}_i{v}_{QM+1}^k\&rsqb;{E_i}^T)\end{array}$

As the above analysis, p be calculated ^k, first to calculate and then sue for peace, finally recompression is the column vector of a Nn dimension and is multiplied by coefficient e ^{j θ}.

It should be noted that sparsely can realize, thus reduce calculated amount, improve counting yield.

The step of described step (43) is as follows:

By U " _nafter expression formula substitutes into, known there is no Explicit Expression.Thus, iterative algorithm is adopted to calculate q here ^k=(I _nn-U " _n) ^-1p ^k.In solution procedure, relate to matrix-vector multiplication computing b=(I _nn-U " _n) y, wherein

First, vectorial y is pressed row boil down to one matrix and then, b=(I can be obtained _nn-U " _n) the sparse performing step of y is as follows:

$\begin{array}{c}b=(I_{Nn}-U_N^{\&prime;\&prime;})y\\ =(I_{Nn}-\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;})vec\left(Y\right)\\ =y-vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{{\mathrm{YF}}_i}^T\right)\\ =y-e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{{\mathrm{YF}}_i}^T\right)\\ =y-e^{j\mathrm{\&theta;}}vec(\underset{i=0}{\overset{m}{\&Sigma;}}\&lsqb;{\tilde{A}}_i{y}_1,\mathrm{...},{\tilde{A}}_i{y}_N\&rsqb;F_i^T)\end{array}$

It should be noted that sparsely can realize, thus reduce calculated amount, improve counting yield.

In described step S5, calculating μ " afterwards, after spectrum mapping, coordinate despining and newton's verification, obtain the eigenvalue λ of time-lag power system successively, before newton's verification, the computing formula of λ is as follows:

$\mathrm{\&lambda;}={\mathrm{\&lambda;}}^{\&prime;\&prime;}{e}^{-j\mathrm{\&theta;}}=\frac{1}{h}{\mathrm{ln\&mu;}}^{\&prime;\&prime;}{e}^{-j\mathrm{\&theta;}}.$

Beneficial effect of the present invention:

The first, when the SOD-PS-R algorithm that proposes of the present invention is for calculating critical eigenvalue corresponding to real system electromechanic oscillation mode, the scale of real system has been taken into full account, and the impact of communication delay.

The second, the SOD-PS-R algorithm that the present invention proposes only need calculate the maximum setting of the modulus value of Solution operator discretize approximate matrix several eigenwert, just can obtain time-lag power system electromechanic oscillation mode characteristic of correspondence value.

Three, the SOD-PS-R algorithm that the present invention proposes obtains the eigenwert of time-lag power system according to the discretize approximate matrix of Solution operator, the degree of accuracy of time-lag power system eigenwert increases with its sensitivity to time lag and reduces.

Four, method of the present invention is the time-lag power system electromechanic oscillation mode characteristic of correspondence value calculating method based on SOD-PS-R, core innovative point is following various means to combine, infinite dimensional time-lag power system is converted to finite dimensional discretize matrix, and it is rotated, the eigenwert eigenwert of discretize approximate matrix being converted to time-lag power system is then mapped by spectrum.

Five, the SOD-PS-R algorithm that the present invention proposes improves convergence of algorithm effect by rotation on the basis of SOD-PS algorithm.

Accompanying drawing explanation

Fig. 1 is time-lag power system schematic diagram.

Fig. 2 is electric system electromechanical oscillations pattern characteristic of correspondence value region.

The graph-based that Fig. 3 (a) and Fig. 3 (b) is spectral mapping theorem.

Fig. 4 (a) and Fig. 4 (b) is the spectrum mapping relations after X-axis rotate.

Fig. 5 is the process flow diagram based on SOD-PS-R time-lag power system electromechanic oscillation mode computing method.

Embodiment

Below in conjunction with accompanying drawing and embodiment, the invention will be further described.

As shown in Figure 1: in the modeling process of actual large-scale electrical power system, introduce communication delay.Time-lag power system comprises without time-lag power system, wide area Feedback Delays, wide area damping control and wide area output time lag four part, and the annexation between each several part as shown in the figure.In Fig. 1, the output without time-lag power system is y _f, y _dfdelayed wide area feedback signal during for considering as the input of damping controller, y _cfor the output of wide area damping control, y _dcfor the control inputs without time-lag power system.

As shown in Figure 2: the operational mode that the eigenwert of time-lag power system is corresponding mainly can be divided three classes, electromechanic oscillation mode, slow pattern, rapid pattern.Wherein electromechanic oscillation mode characteristic of correspondence value region is the object of this Algorithm for Solving.

The eigenwert that the middle time lag system of Fig. 3 (a) is arranged in left half complex plane is mapped within Fig. 3 (b) unit circle, and the eigenwert that time lag system is positioned at right half complex plane is mapped as the eigenwert that Solution operator modulus value is greater than 1, and be positioned at outside unit circle.Therefore, utilize the eigenwert of Solution operator, just can judge the stability of former time lag system.If Solution operator at least exists the eigenwert that a modulus value is greater than 1, then can judge that former time lag system is unstable, if the modulus value of all eigenwerts of Solution operator is all less than 1, then former time lag system is asymptotically stability.

As shown in Fig. 4 (a) He Fig. 4 (b): utilize rotation of coordinate preprocess method, coordinate axis is rotated counterclockwise θ angle, the damping ratio of time-lag power system is less than the critical eigenvalue λ of given constant ζ (=sin θ), is transformed to T _m,Nthe partial feature value μ that modulus value is greater than 1.As shown in the figure, the postrotational imaginary axis correspond to damping ratio in former coordinate system is the dotted line of ζ.

As shown in Figure 5: based on the time-lag power system electromechanic oscillation mode characteristic of correspondence value calculating method of SOD-PS-R, comprise the steps:

S1: set up time-lag power system model;

S2: compose discretize by puppet, obtains the finite dimension discretize matrix T of Solution operator T (h) _m,N;

S3: utilize rotation of coordinate preprocess method, is less than the critical eigenvalue λ of given constant ζ, is transformed to T by the damping ratio of time-lag power system _m,Npostrotational approximate matrix T " _m,Nthe partial feature value μ that modulus value is greater than 1 ";

S4: employing sequential method or subspace method (as implicit restarted Arnoldi algorithm) carry out the Solution operator discretize approximate matrix T that calculation procedure S3 obtains " _m,Nthe maximum setting of modulus value several eigenwert μ ";

S5: calculating μ " afterwards, obtains the eigenvalue λ of time-lag power system successively through spectrum mapping, coordinate despining and newton's verification.

So far, setting several critical eigenvalue that the electromechanic oscillation mode of time-lag power system is corresponding has been calculated.

In described step S1, after considering wide-area communication time-delay, electric system can describe with following one group of time-delayed differential equations:

In formula: for the state variable vector of electric system, n is system state variables sum.T is current time.0< τ ₁< τ ₂< ... < τ _i< τ _mfor the time lag constant of Time Delay, wherein maximum time lag is τ _m. for system state matrix, it is dense matrix. for system time lags state matrix, it is sparse matrix.

The secular equation of the linearized system that above formula represents is:

$({\tilde{A}}_0+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{e}^{-{\mathrm{\&lambda;\&tau;}}_i})v=\mathrm{\&lambda;}v$

In formula: λ and v is respectively eigenwert and corresponding right proper vector.

In described step S2, Solution operator T (h): X → X is defined as the initial condition in space X transfer to h (transfer step-length, 0≤h≤τ _m) linear operator of For Solutions of Systems segmentation after the moment.

From spectral mapping theorem, between the eigenwert μ of Solution operator T (h) and the eigenvalue λ of time lag system, there is following relation.

$\mathrm{\&lambda;}=\frac{1}{t}ln\mathrm{\&mu;},\mathrm{\&mu;}\&Element;\mathrm{\&sigma;}\left(T\right(t\left)\right)\backslash \left\{0\right\}$

In formula: σ (T (h)) represents the spectrum of Solution operator, represent set difference computing.

The eigenwert that time lag system is positioned at left half complex plane is mapped within Solution operator unit circle, and the eigenwert that time lag system is positioned at right half complex plane is mapped as the eigenwert that Solution operator modulus value is greater than 1, and is positioned at outside unit circle.Therefore, utilize the eigenwert of Solution operator, just can judge the stability of former time lag system.If Solution operator at least exists the eigenwert that a modulus value is greater than 1, then can judge that former time lag system is unstable, if the modulus value of all eigenwerts of Solution operator is all less than 1, then former time lag system is asymptotically stability.

Solution operator T (h) is the Infinite Dimension Linear operator describing X → X mapping.In order to calculate the eigenwert of Solution operator and obtain former time lag system electromechanic oscillation mode characteristic of correspondence value according to this, first pseudo-spectral method (Pesudospectral is adopted, PS) discretize is carried out to T (h), obtain corresponding with Solution operator, a finite dimensional approximate matrix, and then calculate the eigenwert of approximate matrix and obtain former time lag system electromechanic oscillation mode characteristic of correspondence value.Discretize matrix T that is corresponding with Solution operator T (h), Standard basis form _m,Ncan be expressed as follows:

$T_{M,N}=T_M^{\left(1\right)}+T_{M,N}^{\left(2\right)}{(I_{Nn}-U_N^{\left(2\right)})}^{-1}{U}_{M,N}^{\left(1\right)}$

In formula:

$T_M^{\left(1\right)}=\left[\begin{array}{ccc}{1}^{M\&times;1}& & \\ I^{(Q-2)M}& & \\ & T_M^Q& 0^{(M+1)\&times;M}\end{array}\right]\&CircleTimes;I_n={\tilde{T}}_M^{\left(1\right)}\&CircleTimes;I_n$

$T_{M,N}^{\left(2\right)}=\left[\begin{array}{c}{T}_{M,N}^1\\ 0^{QM\&times;N}\end{array}\right]\&CircleTimes;I_n={\tilde{T}}_{M,N}^{\left(2\right)}\&CircleTimes;I_n$

In formula, Q, M and N are given positive integer, 1 ^{m × 1}for the M dimensional vector that element is 1 entirely, I ^{(Q-2) M}for (Q-2) M rank unit matrix, 0 ^{(M+1) × M}for (M+1) × M rank null matrix, for Kronecker amasss computing.Matrix for height sparse matrix, and with system state matrix irrelevant. 0 ^{qM × N}for QM × N rank null matrix, matrix for height sparse matrix, and with system state matrix irrelevant.

$U_{M,N}^{\left(1\right)}=\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i$

In formula, its element is determined by Lagrange coefficient completely.

$U_N^{\left(2\right)}=\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i$

In formula, its element is determined by Lagrange coefficient completely.

In described step S3, utilize rotation of coordinate preprocess method, coordinate axis is rotated counterclockwise θ angle, the damping ratio of time-lag power system is less than the critical eigenvalue λ of given constant ζ (=sin θ), is transformed to T _m,Nthe partial feature value μ that modulus value is greater than 1.It is the equal damping ratio of ζ that the postrotational imaginary axis correspond to damping ratio in former coordinate system.

If the eigenwert of time-lag power system is λ ' after rotation of coordinate, so by the λ λ ' e in spectrum mapping relations ^{-j θ}replace, the secular equation after rotation of coordinate can be obtained:

$({\tilde{A}}_0^{\&prime;}+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{e}^{-{\mathrm{\&lambda;}}^{\&prime;}{\mathrm{\&tau;}}_i^{\&prime;}})v={\mathrm{\&lambda;}}^{\&prime;}v$

In formula:

${\tilde{A}}_0^{\&prime;}={\tilde{A}}_0{e}^{j\mathrm{\&theta;}},{\tilde{A}}_i^{\&prime;}={\tilde{A}}_i{e}^{j\mathrm{\&theta;}},i=1,\mathrm{...},m$

λ′＝λe ^jθ

τ _i′＝τ _ie ^-jθ,i＝1,...,m

Due to h be real number, Q is integer, by relation known, the maximum time lag τ ' after rotation of coordinate _maxmust be real number.But, when θ ≠ 0, by τ ' _ithe known τ ' of expression formula _i(i=1 ..., m) must be plural number.Obviously, this and τ ' _maxfor real number is conflicting.Therefore, in the secular equation after rotation of coordinate, get:

τ′ _i＝τ _ie ^-jθ≈τ _i,i＝1,...,m

Thus, obtain approximation characteristic equation as follows:

$({\tilde{A}}_0^{\&prime;}+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{e}^{-{\mathrm{\&lambda;}}^{\&prime;\&prime;}{\mathrm{\&tau;}}_i})v={\mathrm{\&lambda;}}^{\&prime;\&prime;}v$

" be the approximate value of λ ' in formula: λ.

The Solution operator discretize matrix T that above formula is corresponding " _m,Ncan be expressed as:

$T_{M,N}^{\&prime;\&prime;}=T_M^{\left(1\right)}+T_{M,N}^{\left(2\right)}{(I_{Nn}-U_N^{\&prime;\&prime;})}^{-1}{U}_{M,N}^{\&prime;\&prime;}$

In formula:

$U_{M,N}^{\&prime;\&prime;}=\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}$

$U_N^{\&prime;\&prime;}=\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}$

Obviously, when anglec of rotation θ=0 °, T " _m,Njust deteriorate to T _m,N.

If T " _m,Neigenwert to be μ meet following relation between ", itself and λ ":

${\mathrm{\&lambda;}}^{\&prime;\&prime;}=\frac{1}{h}{\mathrm{ln\&mu;}}^{\&prime;\&prime;},{\mathrm{\&mu;}}^{\&prime;\&prime;}\&Element;\mathrm{\&sigma;}\left(T_{M,N}^{\&prime;\&prime;}\right)\backslash \left\{0\right\}$

In described step S4, matrix T " _m,Nexponent number be (QM+1) n.For large-scale electrical power system, matrix T " _m,Nexponent number will be very huge.Therefore, in the discretize matrix T that application Solution operator is corresponding " _m,Nwhen solving the time lag eigenwert of large-scale electrical power system, iterative characteristic value-based algorithm (sequential method or subspace method) must be adopted to calculate the maximum setting of its modulus value several eigenwert.

Suppose, when the secondary iteration of kth, to need to calculate T " _m,Nwith vector product, w ^k=T " _m, concrete steps are as follows:

(1): by vector v ^kby row boil down to one matrix correspondingly, have: v ^k=vec (V ^k), wherein vec () is for being the computing of column vector by matrix compression.

(2): calculate

(3): calculate

(4): calculate $w^k=T_M^{\left(1\right)}{v}^k+T_{M,N}^{\left(2\right)}{q}^k.$

The step of described step (2) is as follows:

By U " _m,Nexpression formula substitutes into, and can obtain:

$\begin{array}{c}{p}^k-U_{M,N}^{\&prime;\&prime;}{v}^k\\ =(\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;})vec\left(V^k\right)\\ =vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{V}^k{E_i}^T\right)\\ =e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{V}^k{E_i}^T\right)\\ =e^{j\mathrm{\&theta;}}vec(\underset{i=0}{\overset{m}{\&Sigma;}}\&lsqb;{\tilde{A}}_i{v}_1^k,\mathrm{...},{\tilde{A}}_i{v}_{QM+1}^k\&rsqb;{E_i}^T)\end{array}$

As the above analysis, p be calculated ^k, first to calculate and then sue for peace, finally recompression is the column vector of a Nn dimension and is multiplied by coefficient e ^{j θ}.

It should be noted that sparsely can realize, thus reduce calculated amount, improve counting yield.

The step of described step (3) is as follows:

By U " _nafter expression formula substitutes into, known there is no Explicit Expression.Thus, iterative algorithm is adopted to calculate q here ^k=(I _nn-U " _n) ^-1p ^k.In solution procedure, relate to matrix-vector multiplication computing b=(I _nn-U " _n) y, wherein

First, vectorial y is pressed row boil down to one matrix and then, b=(I can be obtained _nn-U " _n) the sparse performing step of y is as follows:

$\begin{array}{c}b=(I_{Nn}-U_N^{\&prime;\&prime;})y\\ =(I_{Nn}-\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;})vec\left(Y\right)\\ =y-vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{{\mathrm{YF}}_i}^T\right)\\ =y-e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{{\mathrm{YF}}_i}^T\right)\\ =y-e^{j\mathrm{\&theta;}}vec(\underset{i=0}{\overset{m}{\&Sigma;}}\&lsqb;{\tilde{A}}_i{y}_1,\mathrm{...},{\tilde{A}}_i{y}_N\&rsqb;F_i^T)\end{array}$

It should be noted that sparsely can realize, thus reduce calculated amount, improve counting yield.

In described step S5, calculating μ " afterwards, after spectrum mapping, coordinate despining and newton's verification, obtain the eigenvalue λ of time-lag power system successively, before newton's verification, the computing formula of λ is as follows:

$\mathrm{\&lambda;}={\mathrm{\&lambda;}}^{\&prime;\&prime;}{e}^{-j\mathrm{\&theta;}}=\frac{1}{h}{\mathrm{ln\&mu;}}^{\&prime;\&prime;}{e}^{-j\mathrm{\&theta;}}.$

By reference to the accompanying drawings the specific embodiment of the present invention is described although above-mentioned; but not limiting the scope of the invention; one of ordinary skill in the art should be understood that; on the basis of technical scheme of the present invention, those skilled in the art do not need to pay various amendment or distortion that creative work can make still within protection scope of the present invention.

Claims (10)

1., based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, comprise the steps:

S1: set up time-lag power system model; Relation between the eigenwert of foundation time-lag power system model and time-lag power system solution to model operator eigenwert, changes into the eigenwert calculating Solution operator by the eigenwert calculating time-lag power system model; Thus the problem calculating time-lag power system electromechanic oscillation mode is converted into the eigenvalue problem calculating Solution operator;

S2: carry out discretize to Solution operator by pseudo-spectral method, obtains the discretize matrix of Solution operator;

S3: utilize rotation of coordinate preprocess method, the discretize matrix rotation of Solution operator is carried out approximate obtaining postrotational Solution operator discretize approximate matrix, the damping ratio of time-lag power system model is less than the critical eigenvalue λ of setting constant ζ, is transformed to the eigenwert μ that corresponding postrotational Solution operator discretize approximate matrix modulus value is greater than 1 ";

S4: adopt sequential method or subspace method, from the eigenwert μ of step S3 " in find the modulus value of Solution operator discretize approximate matrix maximum setting several eigenwert μ ";

S5: obtain eigenwert μ from step S4 " afterwards, successively by μ ", through spectrum mapping, coordinate despining and newton's checking treatment, finally obtains the eigenvalue λ of time-lag power system model, the electromechanic oscillation mode of eigenvalue λ and corresponding time-lag power system.

2., as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, the time-lag power system model in described step S1 is as follows:

In formula: for the state variable vector of electric system, n is system state variables sum; T is current time; 0< τ ₁< τ ₂< ... < τ _i< τ _mfor the time lag constant of Time Delay, wherein maximum time lag is τ _m; for system state matrix, it is dense matrix; for system time lags state matrix, it is sparse matrix; Δ x (t) is the increment of t system state variables, Δ x (t-τ _i) be t-τ _ithe increment of moment system state variables, for the increment of t system state variables derivative; The initial value that Δ x (0) is system state variables, and be abbreviated as

The secular equation of the linearized system that formula (1.1) represents is:

$({\tilde{A}}_0+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{e}^{-{\mathrm{\&lambda;\&tau;}}_i})v=\mathrm{\&lambda;}v---\left(1.2\right)$

In formula: λ is eigenwert, v is right proper vector corresponding to eigenwert.

3., as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, in described step S2, Solution operator T (h): X → X is defined as the initial condition in space X the linear operator of time-lag power system solution segmentation after transferring to the h moment; Wherein, h is transfer step-length, 0≤h≤τ _m;

Wherein, s is integration variable, and θ is variable, with be respectively the state of 0 and h+ θ moment time-lag power system.

4., as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, the relation between the eigenwert of described time-lag power system model and time-lag power system solution to model operator eigenwert:

From spectral mapping theorem, between the eigenwert μ of Solution operator T (h) and the eigenvalue λ of time lag system, there is following relation:

$\mathrm{\&lambda;}=\frac{1}{t}ln\mathrm{\&mu;},\mathrm{\&mu;}\&Element;\mathrm{\&sigma;}\left(T\right(t\left)\right)\backslash \left\{0\right\}---\left(1.4\right)$

In formula: σ (T (h)) represents the spectrum of Solution operator, represent set difference computing.

5., as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, the step of described step S2 is as follows:

Discretize matrix T that is corresponding with Solution operator T (h), Standard basis form _m,Nbe expressed as follows:

$T_{M,N}=T_M^{\left(1\right)}+T_{M,N}^{\left(2\right)}{(I_{Nn}-U_N^{\left(2\right)})}^{-1}{U}_{M,N}^{\left(1\right)}---\left(1.5\right)$

In formula:

$T_M^{\left(1\right)}=\left[\begin{array}{cccc}{1}^{M\&times;1}& & & \\ I^{(Q-2)M}& & & \\ & & & \\ & & T_M^Q& 0^{(M+1)\&times;M}\end{array}\right]\&CircleTimes;I_n={\tilde{T}}_M^{\left(1\right)}\&CircleTimes;I_n---\left(1.6\right)$

$T_{M,N}^{\left(2\right)}=\left[\begin{array}{c}{T}_{M,N}^1\\ 0^{QM\&times;N}\end{array}\right]\&CircleTimes;I_n={\tilde{T}}_{M,N}^{\left(2\right)}\&CircleTimes;I_n---\left(1.7\right)$

$U_{M,N}^{\left(1\right)}=\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i---\left(1.8\right)$

$U_N^{\left(2\right)}=\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i---\left(1.9\right)$

In formula (1.5) and formula (1.6), Q, M and N are given positive integer, 1 ^{m × 1}for the M dimensional vector that element is 1 entirely, I ^{(Q-2) M}for (Q-2) M rank unit matrix, 0 ^{(M+1) × M}for (M+1) × M rank null matrix, for Kronecker amasss computing; Matrix for height sparse matrix, and with system state matrix irrelevant;

In formula (1.7), 0 ^{qM × N}for QM × N rank null matrix, matrix for height sparse matrix, and with system state matrix irrelevant;

In formula (1.8), i=0 ..., m, its element is determined by Lagrange coefficient completely;

In formula (1.9), i=0 ..., m, its element is determined by Lagrange coefficient completely.

6. as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, in described step S3, utilize rotation of coordinate preprocess method, coordinate axis is rotated counterclockwise θ angle, the damping ratio of time-lag power system model is less than the critical eigenvalue λ of setting constant ζ, is transformed to T _m,Nthe eigenwert μ that modulus value is greater than 1; It is the equal damping ratio of ζ that the postrotational imaginary axis correspond to damping ratio in former coordinate system;

If the eigenwert of time-lag power system model is λ ' after rotation of coordinate, so the λ λ ' e in formula (1.2) ^{-j θ}replace, the secular equation after rotation of coordinate can be obtained:

$\left({\tilde{A}}_0^{\&prime;}+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{e}^{-{\mathrm{\&lambda;}}^{\&prime;}{\mathrm{\&tau;}}_i^{\&prime;}}\right)v={\mathrm{\&lambda;}}^{\&prime;}v---\left(1.10\right)$

In formula:

${\tilde{A}}_0^{\&prime;}={\tilde{A}}_0{e}^{j\mathrm{\&theta;}},{\tilde{A}}_i^{\&prime;}={\tilde{A}}_i{e}^{j\mathrm{\&theta;}},i=1,...,m$

λ′＝λe ^jθ

τ′ _i＝τ _ie ^-jθ,i＝1,...,m

Due to h be real number, Q is integer, by relation known, the maximum time lag τ ' after rotation of coordinate _maxmust be real number; But, when θ ≠ 0, by τ ' _ithe known τ ' of expression formula _i(i=1 ..., m) must be plural number; Obviously, this and τ ' _maxfor real number is conflicting; Therefore, in the secular equation after rotation of coordinate, get:

τ′ _i＝τ _ie ^-jθ≈τ _i,i＝1,...,m

Thus, obtain approximation characteristic equation as follows:

$({\tilde{A}}_0^{\&prime;}+\underset{i=1}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{e}^{-{\mathrm{\&lambda;}}^{\&prime;\&prime;}{\mathrm{\&tau;}}_i})v={\mathrm{\&lambda;}}^{\&prime;\&prime;}v---\left(1.11\right)$

" be the approximate value of λ ' in formula: λ;

The approximate matrix T of discretize matrix after the Solution operator rotation that formula (1.11) is corresponding " _m,Nbe expressed as:

$T_{M,N}^{\&prime;\&prime;}=T_M^{\left(1\right)}+T_{M,N}^{\left(2\right)}{(I_{Nn}-U_N^{\&prime;\&prime;})}^{-1}{U}_{M,N}^{\&prime;\&prime;}$

In formula:

$U_{M,N}^{\&prime;\&prime;}=\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}$

$U_N^{\&prime;\&prime;}=\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}$

Obviously, when anglec of rotation θ=0 °, T " _m,Njust deteriorate to T _m,N;

If T " _m,Neigenwert to be μ meet following relation between ", itself and λ ":

${\mathrm{\&lambda;}}^{\&prime;\&prime;}=\frac{1}{h}{\mathrm{ln\&mu;}}^{\&prime;\&prime;},{\mathrm{\&mu;}}^{\&prime;\&prime;}\&Element;\mathrm{\&sigma;}\left(T_{M,N}^{\&prime;\&prime;}\right)\backslash \left\{0\right\}.$

7., as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, the step of described step S4 is as follows:

Suppose, when the secondary iteration of kth, to need to calculate T " _m,Nwith vector product, concrete steps are as follows:

Step (41): by vector v ^kby row boil down to one matrix i=1 ..., QM+1; Correspondingly, have: v ^k=vec (V ^k), wherein vec () is for being the computing of column vector by matrix compression;

Step (42): calculate

Step (43): calculate

Step (44): calculate $w^k=T_M^{\left(1\right)}{v}^k+T_{M,N}^{\left(2\right)}{q}^k.$

8., as claimed in claim 7 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, the step of described step (42) is as follows:

By U " _m,Nexpression formula substitutes into, and can obtain:

$\begin{array}{c}{p}^k=U_{M,N}^{\&prime;\&prime;}{v}^k\\ =\left(\underset{i=0}{\overset{m}{\&Sigma;}}{E}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}\right)vec\left(V^k\right)\\ =vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{V}^k{E_i}^T\right)\\ =e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{V}^k{E_i}^T\right)\\ =e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}\&lsqb;{\tilde{A}}_i{v}_1^k,...,{\tilde{A}}_i{v}_{QM+1}^k\&rsqb;{E_i}^T\right)\end{array}$

As the above analysis, p be calculated ^k, first to calculate i=0 ..., m, and then sue for peace, finally recompression is the column vector of a Nn dimension and is multiplied by coefficient e ^{j θ};

It should be noted that sparsely can realize, thus reduce calculated amount, improve counting yield.

9., as claimed in claim 7 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, the step of described step (43) is as follows:

By U " _nafter expression formula substitutes into, known there is no Explicit Expression; Thus, iterative algorithm is adopted to calculate q here ^k=(I _nn-U " _n) ^-1p ^k; In solution procedure, relate to matrix-vector multiplication computing b=(I _nn-U " _n) y, wherein

First, vectorial y is pressed row boil down to one matrix i=1 ..., N; And then, b=(I can be obtained _nn-U " _n) the sparse performing step of y is as follows:

$\begin{array}{c}b=\left(I_{Nn}-U_N^{\&prime;\&prime;}\right)y\\ =\left(I_{Nn}-\underset{i=0}{\overset{m}{\&Sigma;}}{F}_i\&CircleTimes;{\tilde{A}}_i^{\&prime;}\right)vec\left(Y\right)\\ =y-vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i^{\&prime;}{{\mathrm{YF}}_i}^T\right)\\ =y-e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}{\tilde{A}}_i{{\mathrm{YF}}_i}^T\right)\\ =y-e^{j\mathrm{\&theta;}}vec\left(\underset{i=0}{\overset{m}{\&Sigma;}}\&lsqb;{\tilde{A}}_i{y}_1,...,{\tilde{A}}_i{y}_N\&rsqb;{F_i}^T\right)\end{array}$

It should be noted that sparsely can realize, thus reduce calculated amount, improve counting yield.

10. as claimed in claim 1 based on the time-lag power system electromechanic oscillation mode computing method of SOD-PS-R algorithm, it is characterized in that, in described step S5, calculating μ " afterwards; after spectrum mapping, coordinate despining and newton's verification, obtain the eigenvalue λ of time-lag power system successively; before newton's verification, the computing formula of λ is as follows:

$\mathrm{\&lambda;}={\mathrm{\&lambda;}}^{\&prime;\&prime;}{e}^{-j\mathrm{\&theta;}}=\frac{1}{h}{\mathrm{ln\&mu;}}^{\&prime;\&prime;}{e}^{-j\mathrm{\&theta;}}.$