CN105449665A - Time lag electric power system stability discrimination method based on SOD-PS - Google Patents

Abstract

The invention discloses a time lag electric power system stability discrimination method based on pseudospectral discretization of solution operators (SOD-PS). The method comprises the following steps of: (1) establishing a time lag electric power system model; according to a relation between characteristic values of the time lag electric power system model and solution operator characteristic values of the time lag electric power system model, converting the calculation of the characteristic value of the time lag electric power system model into the calculation of the solution operator characteristic value, and converting a time lag electric power system stability judging problem into a problem of calculating a maximum characteristic value of module values of the solution operators; (2) utilizing a pseudospectral method to discretize the solution operators, and obtaining discretized matrixes of the solution operators; (3), utilizing a sequential method or a sub-space method to calculate the maximum characteristic value [mu] of the module values of the discretized matrixes; and (4) judging the stability of the time lag electric power system according to the magnitude of the characteristic value [mu]. The method has the advantages that the calculation amount is small, and the judgment is rapid and accurate.

Description

Based on the time-lag power system stability method of discrimination of SOD-PS

Technical field

The present invention relates to a kind of time-lag power system stability method of discrimination based on SOD-PS.

Background technology

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 [[1] WuHX of system unstability, TsakalisKS, HeydtGT.Evaluationoftimedelayeffectstowide-areapowersyst emstabilizerdesign.IEEETrans.PowerSyst., 2004,19 (4): 1935-1941.].Therefore, when utilizing wide area measurement information to carry out electric power system closed-loop control, the impact of time lag must be taken into account.

[[2] ox is newborn in invention, Ye Hua, Wang Chunyi, Deng. the time-lag power system characteristic value approximate based on Pad é calculates and Convenient stable criterion .201210271783.8:[P] .] 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.

Invention [[3] Ye Hua, Wang Yanyan, Liu Yutian. based on the extensive time-lag power system characteristic value computational methods .201510055743.3. China of EIGD, 201510055743.3 [P] .] propose a kind of extensive time-lag power system characteristic value calculating based on display IGD (ExplicitIGD, EIGD).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, there is provided a kind of based on SOD-PS (pseudo-spectrum discretization Solution operator (Pseudospectraldiscretizationofsolutionoperator, SOD-PS) time-lag power system stability method of discrimination), the method only needs to calculate a maximum characteristic value of Solution operator discretization matrix norm value or a pair Con-eigenvalue, just can judge the stability of system under fixed time lag.There is amount of calculation little, differentiate feature accurately.

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

Based on the Power System Delay Convenient stable criterion of SOD-PS, comprise the steps:

Step (1): set up time-lag power system model; Relation between the characteristic value of foundation time-lag power system model and time-lag power system solution to model operator characteristic value, changes into the characteristic value calculating Solution operator by the characteristic value calculating time-lag power system model; Thus calculate the maximum eigenvalue problem of the modulus value of Solution operator by judging that the problem of time-lag power system stability is converted into;

Step (2): adopt pseudo-spectral method to carry out discretization to Solution operator, obtains the discretization matrix of Solution operator;

Step (3): adopt sequential method or subspace method to carry out the maximum characteristic value μ of the discretization matrix norm value of the Solution operator that calculation procedure (2) obtains;

Step (4): the stability judging time-lag power system according to the size of characteristic value μ.

The step of described step (4) is:

If the modulus value of characteristic value μ is greater than 1, then time-lag power system is in minor interference unsure state;

If the modulus value of characteristic value μ equals 1, then time-lag power system is in the state of neutrality;

If the modulus value of characteristic value μ is less than 1, then time-lag power system is in the state of asymptotically stability.

The time-lag power system model of described step (1) is as follows:

In formula: for the state variable vector of electric power 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 mode 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. initial condition) that Δ x (0) is system state variables, and be abbreviated as

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

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

In formula: λ is characteristic value, v is right characteristic vector corresponding to characteristic value.

Described Solution operator is defined as follows:

Solution operator T (h) is defined as initial condition transfer to h (transfer step-length, 0≤h≤τ _m) linear operator of time-lag power system solution segmentation after the moment.

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 characteristic value of described time-lag power system model and Solution operator characteristic value:

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

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

In formula: σ (T (h)) represents the spectrum of Solution operator.

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

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

T _M,N＝Τ _M+Τ′ _M,N*(I _Nn-U _N) ^-1*U _M,N(1.5)

In formula:

$T_M=\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\&CircleTimes;I_n---\left(1.6\right)$

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

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

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

In formula (1.5), M and N is given positive integer, I _nnfor Nn rank unit matrix (diagonal entry is 1, and other elements are 0), subscript-1 representing matrix inversion operation.

In formula (1.6), Q is 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, I _nfor n rank unit matrix, 0 ^{(M+1) × M}for (M+1) × M rank null matrix, for Kronecker amasss computing.Matrix Τ _mfor height sparse matrix, and with time-lag power system state matrix ${\tilde{A}}_i(i=0,...,m)$ Irrelevant.

In formula (1.7), 0 ^{qM × N}for QM × N rank null matrix,

In formula (1.8), i=0 ..., m, E _ielement determined by Lagrange coefficient completely.

In formula (1.9), i=0 ..., m, F _ielement determined by Lagrange coefficient completely.

The step of described step (3) 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 (3-1): 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 (3-2): calculate

Step (3-3): calculate

Step (3-4): calculate w ^k=Τ _m* v ^k+ Τ ' _m,N* q ^k.

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

Formula (1.8) is substituted into, can obtain:

$p^k=U_{M,N}*v^k=\mathrm{vec}(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\tilde{A}}_i*V^k*E_i^T)$

In formula: subscript T representing matrix transposition.K is kth time iteration; M is the number of time lag.

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.

It should be noted that sparsely can realize, thus reduce amount of calculation, improve computational efficiency.

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

After formula (1.9) is substituted 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:

$b=(I_{\mathrm{Nn}}-U_N)*y=y-\mathrm{vec}(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\tilde{A}}_i*Y*{F_i}^T)$

It should be noted that sparsely can realize, thus reduce amount of calculation, improve computational efficiency.

Beneficial effect of the present invention:

The method only needs to calculate a maximum characteristic value of Solution operator discretization matrix norm value or a pair Con-eigenvalue, just can judge the stability of system under fixed time lag.There is amount of calculation little, differentiate feature fast and accurately.

Accompanying drawing explanation

Fig. 1 is method flow diagram of the present invention;

Fig. 2 (a) is time-lag system characteristic value;

Fig. 2 (b) is Solution operator characteristic value;

Fig. 3 is four machine two district systems;

Fig. 4 is matrix T _m,Nspectrum, i.e. σ (T _m,N);

Fig. 5 is matrix spectrum, namely

Embodiment

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

As shown in Figure 1, based on the Power System Delay Convenient stable criterion of SOD-PS, comprise the steps:

Step (1): set up time-lag power system model;

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

Step (3): the characteristic value μ that adopt sequential method or subspace method (as implicit restarted Arnoldi algorithm) Solution operator discretization matrix norm value that calculation procedure (2) obtains is maximum;

Step (4):

If the modulus value of characteristic value μ is greater than 1, then time-lag power system is in minor interference unsure state;

If the modulus value of characteristic value μ equals 1, then time-lag power system is in the state of neutrality;

If the modulus value of characteristic value μ is less than 1, then time-lag power system is in the state of asymptotically stability.

1. time-lag power system model

After considering wide-area communication time-delay, one group of time-delayed differential equations that electric power system can be following describes:

In formula: for the state variable vector of electric power 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 mode matrix.

The characteristic equation of the linearized system that formula (1) represents is:

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

In formula: λ and v is respectively characteristic value and corresponding right characteristic vector.

2. Solution operator

Solution operator is defined as initial condition 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 characteristic value μ of Solution operator T (h) and the eigenvalue λ of time-lag system, there is following relation.Its graph-based, as shown in Fig. 2 (a) He Fig. 2 (b).

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

In formula: σ (T (h)) represents the spectrum of Solution operator.

The characteristic value that the middle time-lag system of Fig. 2 (a) is arranged in left half complex plane is mapped within Fig. 2 (b) unit circle, and the characteristic value that time-lag system is positioned at right half complex plane is mapped as the characteristic value that Solution operator modulus value is greater than 1, and be positioned at outside unit circle.Therefore, utilize the characteristic value of Solution operator, just can judge the stability of former time-lag system.

If Solution operator at least exists the characteristic value that a modulus value is greater than 1, then can judge that former time-lag system is unstable,

If the modulus value of all characteristic values of Solution operator is all less than 1, then former time-lag system is asymptotically stability.

Solution operator T (h) is Infinite Dimension Linear operator.In order to the characteristic value calculating Solution operator also judges the stability of time-lag system according to this, first pseudo-spectral method (Pesudospectral is adopted, PS) discretization is carried out to T (h), obtain corresponding with Solution operator, a finite dimensional approximate matrix, and then calculate the characteristic value of approximate matrix and judge the stability of former time-lag system.

3. based on the Solution operator discretization of puppet spectrum

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

T _M,N＝Τ _M+Τ′ _M,N*(I _Nn-U _N) ^-1*U _M,N(1.5)

In formula:

$T_M=\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\&CircleTimes;I_n---\left(1.6\right)$

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

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

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

In 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 Τ _mfor height sparse matrix, and with system mode matrix ${\tilde{A}}_i(i$ $=0,...,m)$ Irrelevant.

In formula (1.7), 0 ^{qM × N}for QM × N rank null matrix, i _nfor n rank unit matrix.

In formula (1.8), i=0 ..., m, E _ielement determined by Lagrange coefficient completely.

In formula (1.9), i=0 ..., m, F _ielement determined by Lagrange coefficient completely.

4. sparse realization

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 discretization matrix T that application Solution operator is corresponding _m,Nwhen judging the time lag stability of large-scale electrical power system, iterative characteristic value-based algorithm (sequential method or subspace method) must be adopted to calculate several maximum characteristic values of its modulus value.

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

(1) 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.

(2) calculate formula (1.8) is substituted into, can obtain:

$p^k=U_{M,N}^{\left(1\right)}*v^k=\mathrm{vec}(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\tilde{A}}_i*V^k*E_i^T)$

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.It should be noted that sparsely can realize, thus reduce amount of calculation, improve computational efficiency.

(3) calculate

After formula (1.9) is substituted into, known there is no Explicit Expression.Thus, iterative algorithm is adopted to calculate q here ^k=(I _nn-U _n) ^-1* p ^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:

$b=(I_{\mathrm{Nn}}-U_N)*y=y-\mathrm{vec}(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\tilde{A}}_i*Y*{F_i}^T)$

It should be noted that sparsely can realize, thus reduce amount of calculation, improve computational efficiency.

(4) w is calculated ^k=Τ _m* v ^k+ Τ ' _m,N* q ^k.

Utilize four Ji Liang district example systems, verify the validity of the sparse features value calculating method of the extensive time-lag power system based on display infinitesimal generator discretization that the present invention proposes.All analyses are all carried out in Matlab and on InterCorei54 × 3.4GHz8GBRAM personal computer.

Four Ji Liang district example systems as shown in Figure 3.All generators adopt high-gain Thyristor Excitation System, and installing is with the power system stabilizer, PSS (PSS) of rotating speed deviation for input.On this basis, at G ₁upper installing is with G ₁and G ₃relative rotation speed deviation delta ω ₁₃for the wide area PSS of feedback signal, to improve the damping to inter-area low-frequency oscillation pattern further.Suppose that Feedback Delays and the input delay of wide area PSS are respectively τ ₁=150ms and τ ₂=100ms.The dimension of system state variables and algebraic variable is respectively n=56 and l=22.

(1) approximate matrix T _m,Napproximation capability analysis

First, by Solution operator discretization Spectral radius σ (T _m,N) be converted into the approximate eigenvalue of time-lag power system then utilize newton to verify, obtain time-lag power system characteristic value accurately finally will be converted into σ (T), and by comparing σ (T _m,N) and σ (T), verify discretization matrix T _m,Nthe accuracy of Approximating Solutions operator T.

If the Feedback Delays of wide area damping control is τ=0.2045s, get step-length h=0.05s, then get M=N=20, T _m,Ndimension be (QM+1) n=5656.In order to anti-leak-stopping root, QR algorithm is utilized to calculate T _m,Nall Eigenvalues σ (T _m,N), wherein modulus value is greater than the partial feature value of 1e-3 as shown in Figure 4.And then, be converted to the characteristic value of time-lag system and carry out newton's verification, as shown in Figure 5.By the accurate characteristic value of time-lag power system the accurate profile value σ (T) of the time-lag power system Solution operator be converted to, adds in Fig. 4.

As shown in Figure 4, SOD-PS method accurately can calculate the partial feature value that Solution operator modulus value is greater than 0.45.Wherein, the maximum characteristic value of modulus value is μ=0.87495+0.48415i, and corresponding modulus value is 0.99997.Thus, can decision-making system neutrality, the time lag upper limit that system can be born is τ _m=0.2045s.Now, the characteristic value of the time-lag system corresponding with μ is the low frequency oscillation mode λ=-0.0005+10.1082i close to zero damping.

Above-mentioned analysis shows: the approximate matrix of the Solution operator that SOD-PS method obtains, can more adequately Approximating Solutions operator modulus value the best part characteristic value, can the time lag stability of judgement system exactly.

(2) analysis of the accuracy of delay margin calculating

Below by with the contrast utilizing LMI method to calculate the maximum time lag that four machine two district systems can be born, the accuracy of SOD-PS method is described.

(open loop) POWER SYSTEM STATE order of matrix number not comprising wide area damping control is 52 rank.LMI (LinearMatrixInequality in the robust control tool box (RobustControlToolbox) provided in order to Matlab can be utilized, LMI) solver, calculate delay margin [[5] YaoW of (closed loop) electric power system comprising wide area damping control, JiangL, WuQH, Deng .Delay-dependentstabilityanalysisofthepowersystemwithawi de-areadampingcontrollerembedded.IEEETrans.PowerSyst., 2011, 26 (1): 233-240. [6] WuM, HeY, SheJ-H, Deng .Delay-dependentcriteriaforrobuststabilityoftime-varying delaysystems.Automatica, 2004, 40 (8): 1435-1439.], depression of order must be carried out to divided ring system.Utilize the Schur depression of order, be 7 rank by open cycle system depression of order, its transfer function is as follows:

$G\left(s\right)=\frac{-0.55684\&times;{10}^{-4}{s}^6+0.1489s^5-73.45s^4-871.7s^3-3832s^2-1.593\&times;{10}^{4}s-387.1}{s^7+47.64s^6+1055s^5+9986s^4+8.864\&times;{10}^4{s}^3+3.489\&times;{10}^5{s}^2+1.121\&times;{10}^{6}s+2.6\&times;{10}^6}$

Before and after contrast depression of order, the frequency characteristics of open cycle system can find, in 0.2-2.5Hz low-frequency oscillation frequency range, it is very little that reduced-order model and full rank model frequency respond difference.And then utilize the gevp function that robust control tool box provides, the delay margin of trying to achieve system is 187.9ms.The analysis analyzed by the approximation capability of approximate matrix TM, N is known, and the system time lags stability margin utilizing SOD-PS method to calculate is 204.5ms.Thus known, there is larger conservative in LMI method, is (204.5-187.9)/204.5*100%=8.12%.

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 (9)

1., based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that, comprise the steps:

Step (1): set up time-lag power system model; Relation between the characteristic value of foundation time-lag power system model and time-lag power system solution to model operator characteristic value, changes into the characteristic value calculating Solution operator by the characteristic value calculating time-lag power system model; Thus calculate the maximum eigenvalue problem of the modulus value of Solution operator by judging that the problem of time-lag power system stability is converted into;

Step (2): adopt pseudo-spectral method to carry out discretization to Solution operator, obtains the discretization matrix of Solution operator;

Step (3): adopt sequential method or subspace method to carry out the maximum characteristic value μ of the discretization matrix norm value of the Solution operator that calculation procedure (2) obtains;

Step (4): the stability judging time-lag power system according to the size of characteristic value μ.

2., as claimed in claim 1 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that, the step of described step (4) is:

If the modulus value of characteristic value μ is greater than 1, then time-lag power system is in minor interference unsure state;

If the modulus value of characteristic value μ equals 1, then time-lag power system is in the state of neutrality;

If the modulus value of characteristic value μ is less than 1, then time-lag power system is in the state of asymptotically stability.

3., as claimed in claim 1 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

The time-lag power system model of described step (1) is as follows:

In formula: for the state variable vector of electric power 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; (i=0,1 ..., m) be system mode 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 characteristic equation of the linearized system that formula (1.1) represents is:

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

In formula: λ is characteristic value, v is right characteristic vector corresponding to characteristic value.

4., as claimed in claim 1 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

Described Solution operator is defined as follows:

Solution operator T (h) is defined as initial condition transfer to h, transfer step-length, 0≤h≤τ _mthe linear operator of time-lag power system solution segmentation after the moment;

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

5., as claimed in claim 1 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

Relation between the characteristic value of described time-lag power system model and Solution operator characteristic value:

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

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

In formula: σ (T (h)) represents the spectrum of Solution operator.

6., as claimed in claim 1 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

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

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

T _M,N＝Τ _M+Τ′ _M,N*(I _Nn-U _N) ^-1*U _M,N(1.5)

In formula:

$T_M=\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\&CircleTimes;I_n---\left(1.6\right)$

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

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

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

In formula (1.5), M and N is given positive integer, I _nnfor Nn rank unit matrix, subscript-1 representing matrix inversion operation;

In formula (1.6), Q is 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, I _nfor n rank unit matrix, 0 ^{(M+1) × M}for (M+1) × M rank null matrix, for Kronecker amasss computing; Matrix Τ _mfor height sparse matrix, and with time-lag power system state matrix (i=0 ..., m) irrelevant;

In formula (1.7), 0 ^{qM × N}for QM × N rank null matrix, in formula (1.8), i=0 ..., m, E _ielement determined by Lagrange coefficient completely;

In formula (1.9), i=0 ..., m, F _ielement determined by Lagrange coefficient completely.

7., as claimed in claim 1 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

The step of described step (3) 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 (3-1): 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 (3-2): calculate

Step (3-3): calculate

Step (3-4): calculate w ^k=Τ _m* v ^k+ Τ ' _m,N* q ^k.

8., as claimed in claim 7 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

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

Formula (1.8) is substituted into, can obtain:

$p^k=U_{M,N}*v^k=\mathrm{vec}(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\tilde{A}}_i*V^k*E_i^T)$

In formula: subscript T representing matrix transposition; K is kth time iteration; M is the number of time lag

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;

It should be noted that sparsely can realize, thus reduce amount of calculation, improve computational efficiency.

9., as claimed in claim 7 based on the time-lag power system stability method of discrimination of SOD-PS, it is characterized in that,

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

After formula (1.9) is substituted into, known there is no Explicit Expression;

Thus, iterative algorithm is adopted to calculate q here ^k=(I _nn-U _n) ^-1* p ^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:

$b=(I_{\mathrm{Nn}}-U_N)*y=y-\mathrm{vec}(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\tilde{A}}_i*Y*F_i^T)$

It should be noted that sparsely can realize, thus reduce amount of calculation, improve computational efficiency.