CN108280593A - Time-lag power system stability method of discrimination based on IGD-IRK - Google Patents
Time-lag power system stability method of discrimination based on IGD-IRK Download PDFInfo
- Publication number
- CN108280593A CN108280593A CN201810157511.2A CN201810157511A CN108280593A CN 108280593 A CN108280593 A CN 108280593A CN 201810157511 A CN201810157511 A CN 201810157511A CN 108280593 A CN108280593 A CN 108280593A
- Authority
- CN
- China
- Prior art keywords
- formula
- time
- power system
- matrix
- lag power
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Withdrawn
Links
- 238000000034 method Methods 0.000 title claims abstract description 82
- 239000011159 matrix material Substances 0.000 claims abstract description 74
- 239000013598 vector Substances 0.000 claims abstract description 62
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 27
- 238000013507 mapping Methods 0.000 claims abstract description 9
- 238000001228 spectrum Methods 0.000 claims abstract description 9
- 230000006698 induction Effects 0.000 claims abstract description 7
- 238000013016 damping Methods 0.000 claims description 10
- 238000006073 displacement reaction Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 230000009467 reduction Effects 0.000 claims description 5
- 230000009897 systematic effect Effects 0.000 claims description 5
- 238000013459 approach Methods 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 3
- 238000004519 manufacturing process Methods 0.000 claims description 2
- 230000010355 oscillation Effects 0.000 description 10
- 238000004891 communication Methods 0.000 description 4
- 238000005259 measurement Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000003190 augmentative effect Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 230000005611 electricity Effects 0.000 description 2
- 230000005662 electromechanics Effects 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000003381 stabilizer Substances 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000001934 delay Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
- G06Q10/063—Operations research, analysis or management
- G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/06—Energy or water supply
Landscapes
- Business, Economics & Management (AREA)
- Human Resources & Organizations (AREA)
- Engineering & Computer Science (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- General Physics & Mathematics (AREA)
- Development Economics (AREA)
- Health & Medical Sciences (AREA)
- Educational Administration (AREA)
- Marketing (AREA)
- Entrepreneurship & Innovation (AREA)
- Theoretical Computer Science (AREA)
- Tourism & Hospitality (AREA)
- Physics & Mathematics (AREA)
- General Business, Economics & Management (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Game Theory and Decision Science (AREA)
- Public Health (AREA)
- Water Supply & Treatment (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses the time-lag power system stability method of discrimination based on IGD IRK, including：Establish time-lag power system model；Discretization is carried out to infinitesimal generator using implicit Runge Kutta method, obtains the discretization matrix of infinitesimal generator；Convert infinite dimensional eigenvalue problem to the eigenvalue problem of finite dimension；The approximation of the maximum characteristic value of discretization matrix modulus value for the infinitesimal generator being calculated using implicit Arnoldi algorithm；Sparse realization is carried out using the property of Kronecker products, using the product problem of induction dimension-reduction algorithm iterative solution matrix inversion and vector；According to spectrum mapping relations, it converts the approximation of the maximum characteristic value of infinitesimal generator discretization matrix modulus value to the approximate eigenvalue of time-lag power system model；It is modified using Newton iteration method pairing approximation characteristic value, obtains accurate profile value；The stability of time-lag power system is judged according to the size of accurate profile value.
Description
Technical field
The present invention relates to Runge-Kutta discretization (the Infinitesimal Generator based on infinitesimal generator
Discretization Method with Implicit Runge-Kutta, IGD-IRK) time-lag power system stability
Method of discrimination.
Background technology
With the rise of global energy internet, the scale of interconnected electric power system gradually increases, section low-frequency oscillation problem
More significantly.Traditional solution is installation power system stabilizer, PSS (Power System Stabilizer, PSS), still
Since its feedback control signal is from locality, it is unable to the inter-area oscillations of effective damping interconnected electric power system.
Extensive interconnecting electric power is given in the appearance of Wide Area Measurement System (Wide-Area Measurement System, WAMS)
The development of system stability analysis and control brings new opportunity.Interconnected network low-frequency oscillation based on the WAMS Wide-area Measurement Informations provided
Control effectively reflects the wide area feedback signal of inter-area oscillation mode by introducing, can obtain preferable damping control performance,
To solve the problems, such as the inter-area low-frequency oscillation in interconnected network, and then the ability to transmit electricity for improving system provides new control hand
Section has good and is widely applied foreground.
Wide area signal is in the WAMS communications being made of different communication medium (such as optical fiber, telephone wire, digital microwave, satellite)
When transmitting and handle in network, there are the communication delays changed between tens to hundreds of milliseconds.Time lag is that system control law is caused to lose
Effect, operation conditions deteriorate and a kind of major incentive of system unstability.Therefore, electric system closed loop is carried out using wide area measurement information
When control, it is necessary to the influence of meter and time lag.
In modern power systems, concern is primarily with electromechanical oscillations problems for small interference stability.It is with state-space model
The eigenvalue Method on basis is to study the powerful tool of electromechanical oscillations.Currently, many calculating have been proposed in researcher
The effective ways of large-scale electrical power system Critical eigenvalues include mainly the selection mode analytic approach based on reduced order system,
AESOPS algorithms and S- matrix methods, the sequential methods such as power method, Newton method, Rayleigh Rayleigh quotient iterations and simultaneous iterative,
The subspace iteration methods such as Arnoldi algorithm, the double iterative decomposed again and Jacobi-Davidson methods.These methods exist
The sparsity of augmented state matrix is all utilized when calculating section characteristic value, most methods are all by carrying out spectrum change to original system
The distribution for changing to change characteristic spectrum, then seeks the characteristic value of system, then obtain the key feature of original system by inverse transformation
Value.But method mentioned above does not consider the influence of time lag.Chinese invention patent is based on the approximate time lag power trains of Pad é
Characteristic value of uniting calculates and Convenient stable criterion .201210271783.8:[P] approaches time lag ring using Pade approximation polynomials
Section, 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 EIGD
Extensive time-lag power system feature value calculating method .201510055743.3. [P] propose it is a kind of based on display IGD
The extensive time-lag power system characteristic value of (Explicit IGD, EIGD) calculates.The calculated system rightmost side
Critical eigenvalue, it can be determined that stability of the system under fixed time lag.These time lag Convenient stable criterions, are required to pass through
Critical eigenvalue in Multiple-Scan [0.1,2.5] Hz low-frequency oscillations frequency range, close to the imaginary axis, could judge the time lag of system
Stability.
Invention content
In order to solve the deficiencies in the prior art, the present invention provides the time-lag power system stabilities based on IGD-IRK to sentence
Other method；
To achieve the goals above, the present invention adopts the following technical scheme that：
Time-lag power system stability method of discrimination based on IGD-IRK, including：
Step (1)：Establish time-lag power system model；Characteristic value and time lag power train according to time-lag power system model
Spectrum mapping relations between the infinitesimal generator characteristic value of model of uniting convert the characteristic value for calculating time-lag power system model
At the characteristic value for calculating infinitesimal generator；It is infinitely small that the problem of to will determine that time-lag power system stability, is converted into calculating
Generate the maximum eigenvalue problem of modulus value of member；
Step (2)：Discretization is carried out to infinitesimal generator using implicit Runge-Kutta methods, obtains infinitely small generation
The discretization matrix of member；Convert infinite dimensional eigenvalue problem to the eigenvalue problem of finite dimension；
Step (3)：The discretization matrix for the infinitesimal generator that step (2) obtains is calculated using implicit Arnoldi algorithm
The approximation of the maximum characteristic value of modulus value；Sparse realization is carried out using the property of Kronecker products in calculating process, using induction
Dimension-reduction algorithm iteratively solves the product problem of matrix inversion and vector；
Step (4)：According to spectrum mapping relations, by the approximation of the maximum characteristic value of infinitesimal generator discretization matrix modulus value
Value is converted into the approximate eigenvalue of time-lag power system model；
Step (5)：It is modified using Newton iteration method pairing approximation characteristic value, obtains the accurate profile of time-lag power system
Value；
Step (6)：The stability of time-lag power system is judged according to the size of accurate profile value.
The time-lag power system model of the step (1) is：
In formula：For system mode.τi>0, τiFor time lag constant.
Assuming thatWherein τmaxFor maximum time lag.
It is dense matrix for systematic observation matrix；It is height sparse matrix for system time lags state matrix.
IfTo be defined on the n dimensional linear vector spaces of complex field, if state spaceIt is
By delay interval [- τmax, 0] and arrive n dimension real number spacesThe spaces Banach Banach that the continuous function of mapping is constituted, and assign
There is supremum norm
Infinitesimal generatorIt is defined as：
In formula：It is a linear functional：
The step of step (2) is：
According to the definition of infinitesimal generator, its approximate matrix, i.e. the discretization matrix of infinitesimal generator are obtained：
Step (21)：Given positive integer N, utilizes section [- τmax, 0] on the different discrete points of N+1 form set omegasN,
ΩN={ θi, i=0,1..., N }, and then convert continuous state space X to separate manufacturing firms
Step (22)：Given continuous functionIf its discrete approximation isIfIts discrete approximation is
It calculatesAccurate derivative(i.e.) approximation ψ.In discrete point θiPlace, utilizes ψiTo approachAt this
The functional value of point：
In discrete point θ0At=0, obtained using formula (1.5)Accurate derivativeIt willIt is approximately
I.e.
Step (23)：It willApproximation ψiIt is expressed asLinear combination, obtain：
In formula：ajFor constant, dijFor constant.
It is write formula (1.6) as matrix equation, is obtained：
The coefficient of equationAs infinitesimal generatorDiscretization matrix.
For single time lag system, have：θ0=0, θN=-τmax.At this point, matrixFirst block row be simplified.
So far, the thought of the method for infinitesimal generator discretization is discussed, specific derivation formula is as follows：
Provide the discretization method of single time lag system infinitesimal generator based on Radau-IIA first, and then by the party
Method extends to Systems with Multiple Time-Delays.
The implicit Runge-Kutta method recurrence formula of s grades of p ranks is as follows：
It enablesbT=(b1,…,bs), cT=(c1,…,cs)。
Then, in formula (1.8) Runge-Kutta method coefficient (A, b, c) Butcher table Butcher ' s Tableau tables
Show.
There are following features for the coefficient (A, b, c) of Radau-IIA Runge-Kutta methods：
(1)0<c1<…<cs=1；
(2) if method is from being in harmony (consistent),I=1,2 ..., s；
(3) b=[as1,…,ass]T。
Single time lag situation
First, section [- τ, 0] is divided into the subinterval that N number of length is h, h=τ/N；
Then, further division is done to each subinterval with the abscissa of s grades of Runge-Kutta methods；
Finally, the set omega with Ns+1 discrete point is obtainedN。
Utilize set omegaN, convert continuous space X to and discrete turn to space
Given continuous functionIts discrete approximation vector is Φ ∈ XN。
It enablesIts discrete approximation vector is Ψ ∈ XN。
Wherein：
Due to cs=1, have：
In θ0At=0, functionThe exact value ψ of derivative0, obtained by formula (1.3)：
The discrete point θ on+1 subinterval of jth [- (j+1) h ,-jh] (j=1 ..., N-1)j-clH (j=0,1 ..., N-1；
L=1 ..., s) at, by the y in IRK iterative formulasnIt replaces withBy yn+1It replaces with：
It transplants, obtains for formula (1.13)：
Wherein,
By the k in formula (1.15)lIt replaces withBy kiIt replaces with：
It willIt is write as vector form, is obtained：
It enablesUsing variable-definition listed by formula (1.101) and formula (1.102), then formula
(1.17) it is abbreviated as：
Formula (1.18) both sides while premultiplication：
It enables：
By [ψ]j+1It is write as vector form, is obtained：
In view of formula (1.11), then formula (1.21) is written as following reduced form：
In formula：MatrixEnable ωiIndicate i-th of element of vector ω, wiAnd wijRespectively in representing matrix W
I-th row and (i, j) a element, thenExplicit representation is：
Simultaneous formula (1.12) and formula (1.22), obtain infinitesimal generatorDiscretization matrix
In formula：
Multiple time delay situation
First in section [- τmax, 0] on establish include Ns+1 discrete point set omegaN：
In formula： For section [- τi,-τi-1] on NiThe set that s discrete point is constituted.
Utilize set omegaN, continuous space X is converted into discrete space
Given continuous functionIts discrete approximation vector is Φ ∈ XN。
It enablesIts discrete approximation vector is Ψ ∈ XN。
Wherein：
For section [- τi,-τi-1] discrete point on (i=1 ..., m), due to cs=1, then there is following relationship：
Further, since two adjacent interval [- τi+1,-τi] and [- τi,-τi-1] endpoint overlap, then have following relationship：
In θ0At=0, functionThe exact value ψ of derivative0, obtained by formula (1.3).
In i-th of delay interval [- τi,-τi-1]+1 subinterval of jth [- (j+1) hi,-jhi] on discrete point θj,i-clh
Place, by the y in IRK iterative formulasnIt replaces withyn+1It replaces with：
It transplants, obtains for formula (1.29)：
In formula：
By the k in formula (1.31)lIt replaces withkrIt replaces with：
It willIt is write as vector form, is obtained：
It enablesUsing variable-definition listed by formula (1.251) and formula (1.252), then formula
(1.33) it is abbreviated as：
Formula (1.34) both sides while premultiplication：
It enables：
By [ψ]j+1,i(j=0,1 ..., Ni- 1) it is write as vector form, is obtained：
In view of formula (1.26), then formula (1.37) is written as following reduced form：
In formula：MatrixEnable ωiIndicate i-th of element of vector ω, wiI-th row member in representing matrix W
Element, thenIt can explicitly be expressed as：
Formula (1.39) is applied to all delay interval [- τi,-τi-1] (i=1 ..., m), and consider formula (1.27),
：
In formula：Matrix
For i=1 ..., m is enabledThenI.e.：
Simultaneous formula (1.28) and formula (1.40), are derived by the relational expression between Φ and Ψ in the case of multiple time delay：
In formula：The discretization matrix of infinitesimal generator is tieed up for (Ns+1) n × (Ns+1) n.It first
A block rowWrite as unit vectorAnd systematic observation matrix(i=0,1 ..., m)
The sum of Kronecker products.
The step of step (3) is：
First, the eigenvalue λ that time-lag power system is substituted with λ '+s, then can be obtained the characteristic equation after displacement, i.e.,：
In formula：
After displacement operation, infinitesimal generator discretization approximate matrix that IGD-IRK methods obtainIt is mapped asIn turn, inverse matrix is represented by：
In formula：
Then, it is sought using implicit Arnoldi algorithmThe maximum characteristic value of modulus value.
In IRA algorithms, the maximum operation of calculation amount is to utilizeKrylov subspace is formed with vector product.
If k-th of Krylov vector isThen+1 Krylov vectors q ' of kthk+1It calculates as follows：
Due to matrixWithout special logical construction,Without explicit expression form.
For extensive time-lag power system, calculated using direct inversion technique (such as LU is decomposed and Gaussian elimination method)
Inverse matrix when, it is on the one hand very high to request memory, and memory overflow problem may be caused；On the other hand, it cannot make full use of
The sparse characteristic of system augmented state matrix.
In order to avoid direct solutionQ ' is calculated using alternative mannerk+1.Then, formula (1.50) is converted to：
In formula：It is q ' after the l times iterationk+1Approximation.
The advantage of iterative solution is during solving system of linear equations, not increasing any element, maintain's
Sparse characteristic.
(q ' is calculated using induction dimension reduction methodk+1)(l), it is as follows.
First, by (q 'k+1)(l)In element rearranged according to the direction of row, obtain matrixThat is (q 'k+1)(l)=vec (Q ').In turn, the property accumulated using Kronecker, formula
(1.51) left end is calculated as：
In formula:
In formula (1.52), the maximum operation of calculation amount is matrix-vector product
It is calculated using the power method of sparse realization, reduces computation burden, improves computational efficiency.
It is calculated using the power method of sparse realization：
Given convergence precision ε1, then (q ' is solvedk+1)(l)The conditions of convergence of IDR (s) algorithms be：
The step of step (4) is：If IRA algorithms are calculatedCharacteristic value be λ ", then time lag electric power
The approximate eigenvalue of system modelFor：
The step of step (5) is：The vector that the preceding n element of formula (1.54) Krylov vectors corresponding with λ " is formedIt is the approximation of the corresponding feature vector v of accurate profile value λ.WithWithFor initial value, essence is obtained by iteration using Newton method
True eigenvalue λ and corresponding feature vector v.
The step of step (6) is：The stability of time-lag power system is judged according to the size of accurate profile value.
If the damping ratio of some eigenvalue λ is less than 0, time-lag power system is in small interference unsure state；
If the minimum damping ratio of all characteristic values is equal to 0, time-lag power system is in the state of neutrality；
If the damping ratio of all characteristic values is all higher than 0, time-lag power system is in the state of asymptotically stability.
Compared with prior art, the beneficial effects of the invention are as follows：
The first, IGD-IRK algorithms proposed by the present invention are corresponding crucial special for calculating real system electromechanic oscillation mode
When value indicative, the scale of real system and the influence of communication delay have been fully considered.
The second, IGD-IRK algorithms proposed by the present invention calculate the partial feature value near designated displacement point, so that it may to obtain
The corresponding characteristic value of time-lag power system electromechanic oscillation mode, computational accuracy are more accurate.
Third, IGD-IRK algorithms proposed by the present invention utilize implicit Runge-Kutta discretization scheme, have expanded extensive
The computational methods of time-lag power system critical eigenvalue, and its obtained critical eigenvalue contributes to wide area damping control
Design.
4th, the present invention is based on the time-lag power system key feature value calculating methods of IGD-IRK so that follow conventional lines electricity
The Eigenvalues analysis theoretical method and frame of Force system small signal stability, analyze time-lag power system small signal stability and
Design wide area damping control is possibly realized.This is theoretical for the analysis on Small Disturbance Stability of perfect and abundant feature based value,
Have great importance and is worth.
Description of the drawings
The accompanying drawings which form a part of this application are used for providing further understanding of the present application, and the application's shows
Meaning property embodiment and its explanation do not constitute the improper restriction to the application for explaining the application.
Fig. 1 is the flow chart of the present invention.
In the case of the mono- time lags of Fig. 2, the discrete point set omega in IGD-IRK methodsN；
In the case of Fig. 3 multiple time delay, the discrete point set omega in IGD-IRK methodsN。
Specific implementation mode
It is noted that following detailed description is all illustrative, it is intended to provide further instruction to the application.Unless another
It indicates, all technical and scientific terms used herein has usual with the application person of an ordinary skill in the technical field
The identical meanings of understanding.
It should be noted that term used herein above is merely to describe specific implementation mode, and be not intended to restricted root
According to the illustrative embodiments of the application.As used herein, unless the context clearly indicates otherwise, otherwise singulative
It is also intended to include plural form, additionally, it should be understood that, when in the present specification using term "comprising" and/or " packet
Include " when, indicate existing characteristics, step, operation, device, component and/or combination thereof.
As shown in Figure 1, the time-lag power system stability method of discrimination based on IGD-IRK, including：
Step (1)：Establish time-lag power system model；Characteristic value and time lag power train according to time-lag power system model
Spectrum mapping relations between the infinitesimal generator characteristic value of model of uniting convert the characteristic value for calculating time-lag power system model
At the characteristic value for calculating infinitesimal generator；It is infinitely small that the problem of to will determine that time-lag power system stability, is converted into calculating
Generate the maximum eigenvalue problem of modulus value of member；
Step (2)：Discretization is carried out to infinitesimal generator using implicit Runge-Kutta methods, obtains infinitely small generation
The discretization matrix of member；Convert infinite dimensional eigenvalue problem to the eigenvalue problem of finite dimension；
Step (3)：The discretization matrix for the infinitesimal generator that step (2) obtains is calculated using implicit Arnoldi algorithm
The approximation of the maximum characteristic value of modulus value；Sparse realization is carried out using the property of Kronecker products in calculating process, using induction
Dimension-reduction algorithm iteratively solves the product problem of matrix inversion and vector；
Step (4)：According to spectrum mapping relations, by the approximation of the maximum characteristic value of infinitesimal generator discretization matrix modulus value
Value is converted into the approximate eigenvalue of time-lag power system model；
Step (5)：It is modified using Newton iteration method pairing approximation characteristic value, obtains the accurate profile of time-lag power system
Value；
Step (6)：The stability of time-lag power system is judged according to the size of accurate profile value.
The general type of the implicit Runge-Kutta method recurrence formula of s grades of p ranks is as follows：
It enablesbT=(b1,…,bs), cT=(c1,…,cs).Then, Runge-Kutta method in formula (1.56)
Coefficient (A, b, c) can use Butcher tables (Butcher ' s Tableau) to indicate.
There are following features for the coefficient (A, b, c) of Radau-IIA Runge-Kutta methods：
(1)0<c1<…<cs=1；
(2) if method is from being in harmony (consistent),I=1,2 ..., s；
(3) b=[as1,…,ass]T。
Single time lag situation
1. discrete point set
First, section [- τ, 0] is divided into the subinterval that N number of length is h, h=τ/N；Then, with s grades of Runge-Kuttas
The abscissa of method does further division to each subinterval；Finally, the set omega with Ns+1 discrete point is obtainedN, such as Fig. 2 institutes
Show.
2. variable-definition
Utilize set omegaN, convert continuous space X to and discrete turn to spaceIt is given to connect
Continuous functionIts discrete approximation vector is Φ ∈ XN.It enablesIts discrete approximation vector is Ψ ∈ XN。
In formula：
Due to cs=1, have：
3. in θ0Function at=0Derivative ψ0Estimation
In θ0At=0, functionThe exact value ψ of derivative0, can be obtained by splicing conditional (1.3).
4. in θj-clH (j=0,1 ..., N-1；L=1 ..., s) at functionDerivative [ψ]j+1Estimation
The discrete point θ on+1 subinterval of jth [- (j+1) h ,-jh] (j=1 ..., N-1)j-clH (j=0,1 ..., N-1；
L=1 ..., s) at, by the y in IRK iterative formulasnAnd yn+1It replaces with respectivelyWith：
It transplants, obtains for formula (1.61)：
In formula：
By the k in (1.63) formulal(l=1 ..., s) and ki(i=1 ..., s) it replaces with respectivelyWith
：
It willIt is write as vector form, is obtained：
It enablesUsing variable-definition listed by formula (1.58), then formula (1.65) can be abbreviated as：
Formula (1.66) both sides while premultiplication：
It enablesIt can obtain：
By [ψ]j+1(j=0,1 ..., N-1) is write as vector form, obtains：
In view of formula (1.59), then it is following reduced form that formula (1.69) is rewritable：
In formula：MatrixEnable ωiIndicate i-th of element of vector ω, wiAnd wijRespectively in representing matrix W
I-th row and (i, j) a element, thenCan explicit representation be：
5. the Radau-IIA discretization matrixes of infinitesimal generator
Simultaneous formula (1.60) and formula (1.70), can obtain the discretization matrix of infinitesimal generator
In formula：
7.1.1 multiple time delay situation
The IRK discretization methods of infinitesimal generator in the case of single time lag are expanded to Systems with Multiple Time-Delays by this section.
1. discrete point set
First in section [- τmax, 0] on establish include Ns+1 discrete point set omegaN：
In formula： (i=1 ..., m) it is section [- τi,-τi-1] on NiThe set that s discrete point is constituted,
As shown in Figure 3.
In the case of 1 multiple time delay, the discrete point set omega in IGD-IRK methodsN
2. variable-definition
Utilize set omegaN, continuous space X is converted into discrete space
Given continuous functionIts discrete approximation vector is Φ ∈ XN。
It enablesIts discrete approximation vector is Ψ ∈ XN。
In formula：
For section [- τi,-τi-1] discrete point on (i=1 ..., m), due to cs=1, then there is following relationship：
Further, since two adjacent interval [- τi+1,-τi] and [- τi,-τi-1] endpoint overlap, then have following relationship：
3. in θ0Function at=0Derivative ψ0Estimation
In θ0At=0, functionThe exact value ψ of derivative0, can be obtained by splicing conditional (1.3).
4. i-th of delay interval [- τi,-τi-1] discrete point θj,i-clFunction at h (l=1 ..., s)Derivative [ψ]j,i(j=
0,1,…,Ni- 1) estimation
As shown in figure 3, in i-th of delay interval [- τi,-τi-1]+1 subinterval of jth [- (j+1) hi,-jhi] (j=
0,1,…,Ni- 1) discrete point θ onj,i-clAt h (l=1 ..., s), by the y in IRK iterative formulasnAnd yn+1It replaces with respectivelyWith：
It transplants, obtains for formula (1.77)：
In formula：
By the k in formula (1.79)l(l=1 ..., s) and kr(r=1 ..., s) it replaces with respectivelyWith：
It willIt is write as vector form, is obtained：
It enablesUsing variable-definition listed by formula (1.73), then formula (1.81) can be abbreviated as：
Formula (1.82) both sides while premultiplication：
It enablesIt can obtain：
By [ψ]j+1,i(j=0,1 ..., Ni- 1) it is write as vector form, is obtained：
In view of formula (1.74), then it is following reduced form that formula (1.85) is rewritable：
In formula：MatrixEnable ωiIndicate i-th of element of vector ω, wiI-th row member in representing matrix W
Element, thenIt can explicitly be expressed as：
5. all delay interval [- τi,-τi-1] (i=1 ..., m) discrete point θj,i-clFunction at h (l=1 ..., s)It leads
Number [ψ]j,i(j=0,1 ..., Ni- 1) estimation
Formula (1.87) is applied to all delay interval [- τi,-τi-1] (i=1 ..., m), and consider formula (1.75),
It can obtain：
In formula：Matrix
For i=1 ..., m is enabledThenI.e.：
5. the Radau-IIA discretization matrixes of infinitesimal generator
Simultaneous formula (1.76) and formula (1.88), can be derived by the relational expression between Φ and Ψ in the case of multiple time delay：
In formula：The discretization matrix of infinitesimal generator is tieed up for (Ns+1) n × (Ns+1) n.Its first
Block rowIt can be write as unit vectorAnd systematic observation matrix(i=0,1 ..., m)
The sum of Kronecker products.
7.2 sparse eigenvalues calculate
7.2.1 displacement-inverse transformation
First, the eigenvalue λ that DCPPS is substituted with λ '+s, then can be obtained the characteristic equation after displacement, i.e.,：
In formula：
After displacement operation, infinitesimal generator discretization approximate matrix that IGD-IRK methods obtainIt is mapped as
In turn, inverse matrix is represented by：
In formula：
7.2.2 sparse eigenvalue calculates
It utilizes implicit restarted Arnoldi algorithm (Implicitly Restarted Arnoldi algorithm, IRA)
Etc. Corresponding Sparse Algorithms askModulus value the best part characteristic value.
In IRA algorithms, the maximum operation of calculation amount is to utilizeKrylov subspace is formed with vector product.If
K-th of Krylov vector be respectivelyThen+1 Krylov vectors q ' of kthk+1It can calculate as follows：
Due to matrixWithout special logical construction,Without explicit expression form.
For extensive time-lag power system, calculated using direct inversion technique (such as LU is decomposed and Gaussian elimination method)
Inverse matrix when, it is on the one hand very high to request memory, and memory overflow problem may be caused；On the other hand, it cannot make full use of
The sparse characteristic of system augmented state matrix.
In order to avoid direct solutionUse alternative manner to calculate q ' herek+1.Then, formula (1.98) are converted
For：
In formula：It is q ' after the l times iterationk+1Approximation.
The advantage of iterative solution is during solving system of linear equations, not increasing any element, maintain's
Sparse characteristic.(q ' is calculated using induction dimensionality reduction (Induced Dimension Reduction, IDR (s)) method hereink+1
)(l), it is as follows：
First, by (q 'k+1)(l)In element rearranged according to the direction of row, obtain matrixThat is (q 'k+1)(l)=vec (Q ').In turn, the property accumulated using Kronecker, formula (1.99)
Left end can be calculated as：
In formula:
In formula (1.100), the maximum operation of calculation amount is matrix-vector productCan be used the power method of sparse realization into
Row calculates, and reduces computation burden, improves computational efficiency.It is as follows：
Given convergence precision ε1, then (q ' is solvedk+1)(l)The conditions of convergence of IDR (s) algorithms be：
7.2.3 newton verifies
If IRA algorithms are calculatedCharacteristic value be λ ", thenApproximate eigenvalue be：
" the vector that the preceding n element of corresponding Krylov vectors is formed with λIt is the corresponding feature vectors of accurate profile value λ
The good approximation of v.WithWithFor initial value, accurate profile value λ and corresponding feature can be obtained by iteration using Newton method
Vector v.
The preferred embodiment of upper described only the application, is not intended to limit this application, for the technology of this field
For personnel, the application can have various modifications and variations.Within the spirit and principles of this application, any made by repair
Change, equivalent replacement, improvement etc., should be included within the protection domain of the application.
Claims (8)
1. the time-lag power system stability method of discrimination based on IGD-IRK, characterized in that including：
Step (1)：Establish time-lag power system model；Characteristic value and time-lag power system mould according to time-lag power system model
The characteristic value for calculating time-lag power system model is converted to meter by the spectrum mapping relations between the infinitesimal generator characteristic value of type
Calculate the characteristic value of infinitesimal generator；The problem of to will determine that time-lag power system stability, is converted into the infinitely small generation of calculating
The maximum eigenvalue problem of modulus value of member；
Step (2)：Discretization is carried out to infinitesimal generator using implicit Runge-Kutta methods, obtains infinitesimal generator
Discretization matrix；Convert infinite dimensional eigenvalue problem to the eigenvalue problem of finite dimension；
Step (3)：The discretization matrix modulus value for the infinitesimal generator that step (2) obtains is calculated using implicit Arnoldi algorithm
The approximation of maximum characteristic value；Sparse realization is carried out using the property of Kronecker products in calculating process, using induction dimensionality reduction
Algorithm iteration solution matrix is inverse with vectorial product problem；
Step (4)：According to spectrum mapping relations, the approximation of the maximum characteristic value of infinitesimal generator discretization matrix modulus value is turned
Turn to the approximate eigenvalue of time-lag power system model；
Step (5)：It is modified using Newton iteration method pairing approximation characteristic value, obtains the accurate profile value of time-lag power system；
Step (6)：The stability of time-lag power system is judged according to the size of accurate profile value.
2. the time-lag power system stability method of discrimination based on IGD-IRK as described in claim 1, characterized in that described
The time-lag power system model of step (1) is：
In formula：For system mode；τi>0, τiFor time lag constant；
Assuming thatWherein τmaxFor maximum time lag；
It is dense matrix for systematic observation matrix；It is height sparse matrix for system time lags state matrix；
IfTo be defined on the n dimensional linear vector spaces of complex field, if state space X:=C is by delay interval [- τmax,0]
Real number space is tieed up to nThe spaces Banach Banach that the continuous function of mapping is constituted, and possess supremum norm
3. the time-lag power system stability method of discrimination based on IGD-IRK as claimed in claim 2, characterized in that infinite
Small generation memberIt is defined as：
In formula：It is a linear functional：
The step of step (2) is：
According to the definition of infinitesimal generator, the discretization matrix of infinitesimal generator is obtained：
Step (21)：Given positive integer N, utilizes section [- τmax, 0] on the different discrete points of N+1 form set omegasN, ΩN=
{θi, i=0,1..., N }, and then convert continuous state space X to separate manufacturing firms
Step (22)：Given continuous functionIf its discrete approximation is
IfIts discrete approximation is
It calculatesAccurate derivativeApproximation ψ；In discrete point θiPlace, utilizes ψiTo approachIn the functional value of the point
：
In discrete point θ0At=0, obtained using formula (1.5)Accurate derivativeIt willIt is approximatelyI.e.
Step (23)：It willApproximation ψiIt is expressed asLinear combination, obtain：
In formula：ajFor constant, dijFor constant；
It is write formula (1.6) as matrix equation, is obtained：
The coefficient of equationAs infinitesimal generatorDiscretization matrix；
For single time lag system, have：θ0=0, θN=-τmax；At this point, matrixFirst block row be simplified；
So far, the thought of the method for infinitesimal generator discretization is discussed.
4. the time-lag power system stability method of discrimination based on IGD-IRK as claimed in claim 3, characterized in that
The discretization method of single time lag system infinitesimal generator based on Radau-IIA is provided first, and then this method is pushed away
Extensively to Systems with Multiple Time-Delays；
The implicit Runge-Kutta method recurrence formula of s grades of p ranks is as follows：
It enablesbT=(b1,…,bs), cT=(c1,…,cs)；
Then, the coefficient (A, b, c) of Runge-Kutta method is indicated with Butcher table Butcher ' s Tableau in formula (1.8)；
There are following features for the coefficient (A, b, c) of Radau-IIA Runge-Kutta methods：
(1)0<c1<…<cs=1；
(2) if method is from being in harmony (consistent),
(3) b=[as1,…,ass]T；
Single time lag situation
First, section [- τ, 0] is divided into the subinterval that N number of length is h, h=τ/N；
Then, further division is done to each subinterval with the abscissa of s grades of Runge-Kutta methods；
Finally, the set omega with Ns+1 discrete point is obtainedN；
Utilize set omegaN, convert continuous space X to and discrete turn to space
Given continuous functionIts discrete approximation vector is Φ ∈ XN；
It enablesIts discrete approximation vector is Ψ ∈ XN；
Wherein：
Due to cs=1, have：
In θ0At=0, functionThe exact value ψ of derivative0, obtained by formula (1.3)：
The discrete point θ on+1 subinterval of jth [- (j+1) h ,-jh]j-clAt h, by the y in IRK iterative formulasnIt replaces with
By yn+1It replaces with：
It transplants, obtains for formula (1.13)：
Wherein,
By the k in formula (1.15)lIt replaces withBy kiIt replaces with：
It willIt is write as vector form, is obtained：
It enablesUtilize variable-definition listed by formula (1.101) and formula (1.102), then formula (1.17) letter
It is written as：
Formula (1.18) both sides while premultiplication：
It enables：
By [ψ]j+1It is write as vector form, is obtained：
In view of formula (1.11), then formula (1.21) is written as following reduced form：
In formula：MatrixEnable ωiIndicate i-th of element of vector ω, wiAnd wijRespectively i-th in representing matrix W
Row and (i, j) a element, thenExplicit representation is：
Simultaneous formula (1.12) and formula (1.22), obtain infinitesimal generatorDiscretization matrixXN→XN：
In formula：
Multiple time delay situation
First in section [- τmax, 0] on establish include Ns+1 discrete point set omegaN：
In formula： For section [- τi,-τi- 1] N oniThe set that s discrete point is constituted；
Utilize set omegaN, continuous space X is converted into discrete space
Given continuous functionIts discrete approximation vector is Φ ∈ XN；
It enablesIts discrete approximation vector is Ψ ∈ XN；
Wherein：
For section [- τi,-τi-1] on discrete point, due to cs=1, then there is following relationship：
Further, since two adjacent interval [- τi+1,-τi] and [- τi,-τi-1] endpoint overlap, then have following relationship：
In θ0At=0, functionThe exact value ψ of derivative0, obtained by formula (1.3)；
In i-th of delay interval [- τi,-τi-1]+1 subinterval of jth [- (j+1) hi,-jhi] on discrete point θj,i-clIt, will at h
Y in IRK iterative formulasnIt replaces withyn+1It replaces with：
It transplants, obtains for formula (1.29)：
In formula：
By the k in formula (1.31)lIt replaces withkrIt replaces with：
It willIt is write as vector form, is obtained：
It enablesUtilize variable-definition listed by formula (1.251) and formula (1.252), then formula (1.33) letter
It is written as：
Formula (1.34) both sides while premultiplication：
It enables：
By [ψ]j+1,iIt is write as vector form, is obtained：
In view of formula (1.26), then formula (1.37) is written as following reduced form：
In formula：MatrixEnable ωiIndicate i-th of element of vector ω, wiI-th column element in representing matrix W, thenIt can explicitly be expressed as：
Formula (1.39) is applied to all delay interval [- τi,-τi-1], and consider formula (1.27), it obtains：
In formula：Matrix
For i=1 ..., m is enabledThenI.e.：
Simultaneous formula (1.28) and formula (1.40), are derived by the relational expression between Φ and Ψ in the case of multiple time delay：
In formula：XN→XNThe discretization matrix of infinitesimal generator is tieed up for (Ns+1) n × (Ns+1) n；Its first block rowWrite as unit vectorAnd systematic observation matrixKronecker product the sum of；
5. the time-lag power system stability method of discrimination based on IGD-IRK as claimed in claim 3, characterized in that
The step of step (3) is：
First, the eigenvalue λ that time-lag power system is substituted with λ '+s, then can be obtained the characteristic equation after displacement, i.e.,：
In formula：
After displacement operation, infinitesimal generator discretization approximate matrix that IGD-IRK methods obtainIt is mapped asInto
And inverse matrix is expressed as：
In formula：
Then, it is sought using implicit Arnoldi algorithmThe maximum characteristic value of modulus value；
In IRA algorithms, the maximum operation of calculation amount is to utilizeKrylov subspace is formed with vector product；
If k-th of Krylov vector isThen+1 Krylov vectors q ' of kthk+1It calculates as follows：
In order to avoid direct solutionQ ' is calculated using alternative mannerk+1；Then, formula (1.50) is converted to：
In formula：It is q ' after the l times iterationk+1Approximation；
(q ' is calculated using induction dimension reduction methodk+1)(l), it is as follows：
First, by (q 'k+1)(l)In element rearranged according to the direction of row, obtain matrixThat is (q 'k+1)(l)=vec (Q ')；In turn, the property accumulated using Kronecker, formula
(1.51) left end is calculated as：
In formula:
In formula (1.52), the maximum operation of calculation amount is matrix-vector product
It is calculated using the power method of sparse realization：
Given convergence precision ε1, then (q ' is solvedk+1)(l)The conditions of convergence of IDR (s) algorithms be：
6. the time-lag power system stability method of discrimination based on IGD-IRK as claimed in claim 5, characterized in that
The step of step (4) is：If IRA algorithms are calculatedCharacteristic value be λ ", then time-lag power system
The approximate eigenvalue of modelFor：
7. the time-lag power system stability method of discrimination based on IGD-IRK as claimed in claim 6, characterized in that
The step of step (5) is：The vector that the preceding n element of formula (1.54) Krylov vectors corresponding with λ " is formedIt is
The approximation of the corresponding feature vector v of accurate profile value λ；WithWithFor initial value, essence is obtained by iteration using Newton method
True eigenvalue λ and corresponding feature vector v.
8. the time-lag power system stability method of discrimination based on IGD-IRK as claimed in claim 7, characterized in that
The step of step (6) is：The stability of time-lag power system is judged according to the size of accurate profile value：
If the damping ratio of some eigenvalue λ is less than 0, time-lag power system is in small interference unsure state；
If the minimum damping ratio of all characteristic values is equal to 0, time-lag power system is in the state of neutrality；
If the damping ratio of all characteristic values is all higher than 0, time-lag power system is in the state of asymptotically stability.
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810157511.2A CN108280593A (en) | 2018-02-24 | 2018-02-24 | Time-lag power system stability method of discrimination based on IGD-IRK |
| CN201910132829.XA CN109685400B (en) | 2018-02-24 | 2019-02-22 | Time-lag power system stability discrimination method based on time integral IGD |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810157511.2A CN108280593A (en) | 2018-02-24 | 2018-02-24 | Time-lag power system stability method of discrimination based on IGD-IRK |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108280593A true CN108280593A (en) | 2018-07-13 |
Family
ID=62808683
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810157511.2A Withdrawn CN108280593A (en) | 2018-02-24 | 2018-02-24 | Time-lag power system stability method of discrimination based on IGD-IRK |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108280593A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109726490A (en) * | 2019-01-02 | 2019-05-07 | 华南理工大学 | A kind of more sinusoidal signal design methods of low-frequency range for the identification of POWER SYSTEM STATE spatial model |
-
2018
- 2018-02-24 CN CN201810157511.2A patent/CN108280593A/en not_active Withdrawn
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109726490A (en) * | 2019-01-02 | 2019-05-07 | 华南理工大学 | A kind of more sinusoidal signal design methods of low-frequency range for the identification of POWER SYSTEM STATE spatial model |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108462165B (en) | Load characteristic evaluation method for new energy access power system | |
| CN108321821A (en) | Time-lag power system stability method of discrimination based on SOD-IRK | |
| Tang et al. | Convergence analysis for stochastic collocation methods to scalar hyperbolic equations with a random wave speed | |
| Li et al. | Nonlinear predictors and hybrid corrector for fast continuation power flow | |
| Bleher et al. | Topological expansion in the cubic random matrix model | |
| CN108647906A (en) | Time-lag power system stability analysis method based on low order EIGD | |
| Datta et al. | Higher spin quasinormal modes and one-loop determinants in the BTZ black hole | |
| CN108242808A (en) | Time-lag power system stability method of discrimination based on IGD-LMS | |
| Dong et al. | Effective method to determine time‐delay stability margin and its application to power systems | |
| CN105468909A (en) | Time delay power system electromechanical oscillation mode computing method based on SOD-PS-R R algorithm | |
| CN108054757A (en) | A kind of embedded idle and voltage N-1 Close loop security check methods | |
| Ferlaino et al. | Asymptotic symmetries and thermodynamics of higher spin black holes in AdS 3 | |
| Ye et al. | Iterative infinitesimal generator discretization-based method for eigen-analysis of large delayed cyber-physical power system | |
| CN105046588A (en) | Improved DC (Direct Current) dynamic optimal power flow calculating method based on network loss iteration | |
| CN109685400A (en) | Time-lag power system stability method of discrimination based on time integral IGD | |
| Zhang et al. | Solving ordinary differential equations with adaptive differential evolution | |
| Zacur et al. | Left-invariant riemannian geodesics on spatial transformation groups | |
| CN106372440B (en) | A kind of adaptive robust state estimation method of the power distribution network of parallel computation and device | |
| CN108280593A (en) | Time-lag power system stability method of discrimination based on IGD-IRK | |
| Lahooti et al. | A new fourth order central WENO method for 3D hyperbolic conservation laws | |
| CN108808706B (en) | Time-lag power system electromechanical oscillation mode calculation method based on SOD-PS-II-R algorithm | |
| Gomes et al. | Modal analysis applied to s-domain models of ac networks | |
| Gunes et al. | Proper orthogonal decomposition reconstruction of a transitional boundary layer with and without control | |
| CN108808702A (en) | Time-lag power system electromechanic oscillation mode computational methods based on low order IGD-LMS algorithms | |
| CN110020506B (en) | Differential format selection method based on operation optimization of electric heating type comprehensive energy system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WW01 | Invention patent application withdrawn after publication |
Application publication date: 20180713 |
|
| WW01 | Invention patent application withdrawn after publication |