CN105403834A - Dynamic state evaluation method of generator - Google Patents
Dynamic state evaluation method of generator Download PDFInfo
- Publication number
- CN105403834A CN105403834A CN201510973903.2A CN201510973903A CN105403834A CN 105403834 A CN105403834 A CN 105403834A CN 201510973903 A CN201510973903 A CN 201510973903A CN 105403834 A CN105403834 A CN 105403834A
- Authority
- CN
- China
- Prior art keywords
- generator
- moment
- state
- measurement
- value
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R31/00—Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
- G01R31/34—Testing dynamo-electric machines
- G01R31/343—Testing dynamo-electric machines in operation
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The invention discloses a dynamic state evaluation method of a generator. The method comprises steps of introducing time-varying multi-dimensional observation noise scale factors into CKF dynamic state evaluation; carrying out on-line adjustment on error measurement according to measurement information so as to allow the noise to be closer to actual noise; and using the adjusted errors to calculate filtering gain so as to precisely correct a state quantity predicated value when bad data is measured, thereby obtaining a precise state quantity evaluation value. The invention also provides a specific construction method for the time-varying multi-dimensional observation noise scale factors and provides solutions for abnormity of inversion of filtering gain. In this way, on the premise that requirements on real-time performance are satisfied, effects on dynamic state evaluation of the generator imposed by measurement bad data can be effectively restrained.
Description
Technical field
The present invention relates to Electrical Power System Dynamic State Estimation field, particularly relate to a kind of generator dynamic state estimator method.
Background technology
The transient stability analysis of power system that appears as of synchronous phasor measurement unit (phasormeasurementunit, PMU) provides new technological means with control.But, in the situation such as interference, measuring equipment fault, synchronizing signal loss, often cause the appearance of bad data.The analysis result that bad data may make PMU lead to errors in application process and control strategy.State estimation can reject the bad data existed in measurement amount, therefore, studies dynamic state estimator in the dynamo-electric transient state process of electric system based on PMU most important.
In recent years, based on the electromechanical transient process generator dynamic state estimator problem of PMU be the focus in Electrical Power System Dynamic state estimation field.For having nonlinear generator dynamic equation, its dynamic state estimator problem is a typical nonlinear problem, adopt based on Kalman filtering algorithm scheme be a kind of general resolving ideas, as the dynamic state estimator based on EKF (EKF), Unscented kalman filtering (UKF).The linearization procedure of EKF causes truncation error excessive, and UKF needs to determine parameter value, very flexible, application inconvenience.For this problem, Canadian scholar SimonHaykin proposed volume Kalman filtering (CKF) algorithm in 2009.But no matter UKF or CKF, when equivalent measurement exists bad data, estimated accuracy all can be affected to a certain extent, even causes estimated result substantial deviation actual value, dynamic state estimator failure.
Summary of the invention
The object of this invention is to provide a kind of generator dynamic state estimator method, the method accurately can be revised Generator Status amount predicted value when PMU measurement amount contains bad data, obtains quantity of state estimated value accurately.
The object of the invention is to be achieved through the following technical solutions:
A kind of generator dynamic state estimator method, comprising:
Calculate according to the generator dynamic state estimator result of system equation to the k-1 moment, obtain kth moment generator dynamic state estimator result;
CKF filtering algorithm is utilized to forecast and filtering process kth moment generator dynamic state estimator result; Wherein, the generator dynamical state forecast result in k+1 moment is obtained by forecast process; When filtering process, become multidimensional observation noise scale factor during introducing and the generator dynamical state forecast result of real-time measurement value to the k+1 moment in conjunction with PMU carries out filtering, thus realize the accurate estimation of generator dynamical state result.
Further, described system equation is:
In formula, the expression k+1 moment that subscript k+1 and k is all corresponding and k moment; F and H is respectively state equation function and measurement equation function, and x, u and z are respectively quantity of state, controlled quentity controlled variable and measurement amount; V and w is respectively system noise and measurement noise, and error covariance matrix is respectively the normal distribution of Q and R;
Wherein, quantity of state x=[δ ω E '
qe '
d]
t, controlled quentity controlled variable
the concrete form of state equation is:
In formula, δ and ω is respectively the absolute merit angle of generator amature and angular rate perunit value;
with
the simplification literary style to time diffusion, the expression d/dt differentiating operator of letter top; T
jfor generator inertia time constant; T
mfor machine torque; U
twith
be respectively amplitude and the phase angle of generator outlet voltage phasor; X '
qwith X '
dbe respectively q axle and d axle transient reactance; E '
qwith E '
dbe respectively q axle and d axle transient state electromotive force; D is ratio of damping; T '
q0with T '
d0be respectively q axle and d axle open circuit transient time-constant; E
ffor stator excitation electromotive force; X
qand X
dbe respectively q axle and d axle synchronous reactance;
Measurement amount z=[δ
zω
zp
e]
t, the concrete form of measurement equation is:
In formula, δ
z, ω
zand P
ebe respectively the PMU measuring value of the absolute merit angle of rotor, angular rate and output electromagnetic power;
According to the concrete form of above-mentioned state equation and measurement equation, and k moment system noise variance matrix Q
kwith k+1 moment measuring noise square difference battle array R
k+1, then generator dynamic state estimator can be realized;
K+1 moment measuring noise square difference battle array R wherein
k+1value is carried out, k moment system noise variance matrix Q according to the actual error in measurement of PMU
kbe expressed as:
In formula,
with
be respectively k moment δ, ω, E '
qwith E '
dsystem noise variance;
Calculated by formula of error transmission:
In formula: σ is noise variance; σ
uwith
be respectively amplitude and the phase angle PMU error in measurement standard deviation of generator outlet voltage phasor;
Then have:
In formula, Δ t is step-length.
Further, the step of the described generator dynamical state forecast result by the forecast process acquisition k+1 moment comprises:
Described kth moment generator dynamic state estimator result comprises: the estimated value of kth moment generator dynamic state quantity
with estimation error variance battle array P
k|k;
To P
k|kcarry out Cholesky decomposition, obtain the On Square-Rooting Matrices S of k moment estimation error variance battle array
k|k:
Utilize following formula to the estimated value of kth moment generator dynamic state quantity
generate the volume point X of the weights such as a group
i, k|k:
In formula, parameter
n is quantity of state dimension;
Utilize following formula to convert each volume point, obtain the predicted value of all volumes point
Summation is weighted to the predicted value of all Generator Status amount volume points, obtains quantity of state predicted value
And obtain prediction error conariance battle array P by following formula
k+1|k:
Further, described when filtering process, become multidimensional observation noise scale factor during introducing and the generator dynamical state forecast result of real-time measurement value to the k+1 moment in conjunction with PMU carries out filtering, thus the accurate estimation realizing generator dynamical state result comprises:
The generator dynamical state forecast result in described k+1 moment comprises: the generator dynamic state quantity predicted value in k+1 moment
with Generator Status amount prediction error conariance battle array P
k+1|k;
To prediction error conariance battle array P
k+1|kcarry out Cholesky decomposition, obtain the On Square-Rooting Matrices S of k+1 moment Generator Status amount prediction error conariance battle array
k+1|k:
Utilize following formula to the generator dynamic state quantity predicted value in k+1 moment
generate the quantity of state predicted value volume point X of the weights such as a group
i, k+1|k:
Utilize following formula to convert each quantity of state predicted value volume point, obtain the volume point Z of PMU measurement amount predicted value
i, k+1|k:
Z
i,k+1|k=H(X
i,k+1|k,u
k);
Summation is weighted to the volume point of all PMU measurement amount predicted values, and then obtains PMU measurement amount predicted value
Calculate PMU measurement amount prediction error variance matrix P
vv, k+1:
Utilize the real-time measurement value z obtaining PMU
k+1=[δ
zk+1ω
zk+1p
ek+1]
tcalculate the vectorial e of new breath
k+1, and then when obtaining, become multidimensional observation noise scale factor γ
k+1:
In formula, M is that the window window of the estimation technique is long;
Recycling following formula calculates diagonal matrix γ '
k+1:
γ′
k+1=diag(γ′
1,γ′
2,…,γ′
m);
In formula, diagonal element γ '
1value be: γ '
i=max{1, γ
k+1, ii, i=1,2 ..., n; γ
k+1, iifor γ
k+1i-th diagonal element;
According to the Cross-covariance P between following formula calculating generator quantity of state predicted value and PMU measurement amount predicted value
xz, k+1|k:
Computer card Kalman Filtering gain W again
k+1:
W
k+1=P
xz,k+1|k(P
vv,k+1+γ′
k+1R
k+1)
-1;
The k+1 moment is utilized newly to cease vectorial e
k+1, and by Kalman filtering gain W
k+1to the generator dynamic state quantity predicted value in k+1 moment
carry out filtering, obtain the estimated value of k+1 moment generator dynamic state quantity
And calculate acquisition generator dynamic state quantity estimation error variance battle array P by following formula
k+1|k+1:
As seen from the above technical solution provided by the invention, this method of estimation causes generator dynamic state estimator result to depart from the problem of actual value for bad data in PMU measurement amount, will time become multidimensional observation noise scale factor and be incorporated in CKF filtering algorithm, not only increase the robustness of the method, also can obtain quantity of state estimated value accurately.
Accompanying drawing explanation
In order to be illustrated more clearly in the technical scheme of the embodiment of the present invention, below the accompanying drawing used required in describing embodiment is briefly described, apparently, accompanying drawing in the following describes is only some embodiments of the present invention, for those of ordinary skill in the art, under the prerequisite not paying creative work, other accompanying drawings can also be obtained according to these accompanying drawings.
The outline flowchart of a kind of generator dynamic state estimator method that Fig. 1 provides for the embodiment of the present invention;
The IEEE9 bus test system schematic diagram that Fig. 2 provides for the embodiment of the present invention;
Generator G1 dynamic state estimator result figure in the IEEE9 node emulation testing that Fig. 3 provides for the embodiment of the present invention;
Generator G2 dynamic state estimator result figure in the IEEE9 node emulation testing that Fig. 4 provides for the embodiment of the present invention;
Generator G3 dynamic state estimator result figure in the IEEE9 node emulation testing that Fig. 5 provides for the embodiment of the present invention;
Observation noise scale factor situation of change figure is become during multidimensional in the IEEE9 node emulation testing that Fig. 6 provides for the embodiment of the present invention.
Embodiment
Below in conjunction with the accompanying drawing in the embodiment of the present invention, be clearly and completely described the technical scheme in the embodiment of the present invention, obviously, described embodiment is only the present invention's part embodiment, instead of whole embodiments.Based on embodiments of the invention, those of ordinary skill in the art, not making the every other embodiment obtained under creative work prerequisite, belong to protection scope of the present invention.
The embodiment of the present invention provides a kind of generator dynamic state estimator method, and it mainly comprises the steps:
Calculate according to the generator dynamic state estimator result of system equation to the k-1 moment, obtain kth moment generator dynamic state estimator result;
CKF filtering algorithm is utilized to forecast and filtering process kth moment generator dynamic state estimator result; Wherein, the generator dynamical state forecast result in k+1 moment is obtained by forecast process; When filtering process, become multidimensional observation noise scale factor during introducing and the generator dynamical state forecast result of real-time measurement value to the k+1 moment in conjunction with PMU carries out filtering, thus realize the accurate estimation of generator dynamical state result.
The outline flowchart of generator dynamic state estimator as shown in Figure 1.Generator port voltage can directly directly be measured by PMU, with external network decoupling zero, need not consider network topology structure, thus carry out distributed generator dynamic state estimator.The huge inertia in inherence of generator amature makes generator amature merit angle and angular rate can not undergo mutation in electromechanical transient process, and meets the constraint condition of the generator equation of motion.PMU directly can measure generator's power and angle and angular velocity, utilizes GPS can obtain the synchronous measure value of electromagnetic power; Previous moment quantity of state, in conjunction with generator port voltage phasor, obtains the estimated result of subsequent time Generator Status amount by system equation; Then, obtained the predicted value of quantity of state by forecast, then in conjunction with the measuring value of PMU to quantity of state, by what improve, there is robustness CKF algorithm, finally obtain the dynamic state estimator result of subsequent time generator.
For the ease of understanding, below the present invention is described in detail.
Generator dynamic state estimator problem is Nonlinear Filtering Problem, and the system equation of nonlinear dynamic system can be expressed as:
In formula, the expression k+1 moment that subscript k+1 and k is all corresponding and k moment; F and H is respectively state equation function and measurement equation function, and x, u and z are respectively quantity of state, controlled quentity controlled variable and measurement amount; V and w is respectively system noise and measurement noise, and error covariance matrix is respectively the normal distribution of Q and R;
Be difficult to make timely reaction to secondary transient state process because existing control is standby, therefore, ignore time transient state, dynamic equation reduces to quadravalence by six rank;
Quantity of state x=[δ ω E '
qe '
d]
t, controlled quentity controlled variable
the concrete form of state equation is:
In formula, δ and ω is respectively the absolute merit angle of generator amature and angular rate perunit value;
with
the simplification literary style to time diffusion, the expression d/dt differentiating operator of letter top; T
jfor generator inertia time constant; T
mfor machine torque; U
twith
be respectively amplitude and the phase angle of generator outlet voltage phasor; X '
qwith X '
dbe respectively q axle and d axle transient reactance; E '
qwith E '
dbe respectively q axle and d axle transient state electromotive force; D is ratio of damping; T '
q0with T '
d0be respectively q axle and d axle open circuit transient time-constant; E
ffor stator excitation electromotive force; X
qand X
dbe respectively q axle and d axle synchronous reactance;
Measurement amount z=[δ
zω
zp
e]
t, the concrete form of measurement equation is:
In formula, δ
z, ω
zand P
ebe respectively the PMU measuring value of the absolute merit angle of rotor, angular rate and output electromagnetic power;
According to the concrete form of above-mentioned state equation and measurement equation, and k moment system noise variance matrix Q
kwith k+1 moment measuring noise square difference battle array R
k+1, then generator dynamic state estimator can be realized.
K+1 moment measuring noise square difference battle array R wherein
k+1value is carried out according to the actual error in measurement of PMU.
For system noise, in the system equation of generator dynamic state estimator, hypothesized model parameter is accurate; Quantity of state all adopts estimated value, and precision is higher; T in controlled quentity controlled variable u
m, E
fsteady state value can be assumed to be; U
twith
adopt PMU measuring value, there is error in measurement.Therefore, system noise mainly comes from U
twith
error in measurement.This error in measurement is transmitted by system equation, finally causes the generation of system noise.
K moment system noise variance matrix Q
kbe expressed as:
In formula,
with
be respectively k moment δ, ω, E '
qwith E '
dsystem noise variance;
Calculated by formula of error transmission:
In formula: σ is noise variance;
with
be respectively amplitude and the phase angle PMU error in measurement standard deviation of generator outlet voltage phasor; Its value can be set and be respectively 0.2% and 0.2 °.
According to formula (5), about generator outlet voltage magnitude and phase angle, partial derivative is asked respectively to system equation, and then derivation k moment quantity of state ω, E '
qwith E '
dthe computing formula of process-noise variance.Predicted value due to generator's power and angle δ needs the estimated value using angular rate ω, and ω to send out estimated value error relatively little, therefore, the system noise error variance at merit angle is decided to be a less positive number, elects 10 as herein
-6.
Then have:
Just Q can be calculated by formula (7)-(9)
kin each element value.In formula, Δ t is step-length, U in formula (7)-(9)
twith
value is the measuring value of k moment PMU; E '
q, E '
dwith the estimated value that the value of δ is the k moment.
Recycling CKF filtering algorithm forecasts and filtering process kth moment generator dynamic state estimator result.
1, forecasting process.
Carry out forecasting that process obtains the generator dynamical state forecast result in k+1 moment to kth moment generator dynamic state estimator result;
Described kth moment generator dynamic state estimator result comprises: the estimated value of kth moment generator dynamic state quantity
with estimation error variance battle array P
k|k;
To P
k|kcarry out Cholesky decomposition, obtain the On Square-Rooting Matrices S of k moment estimation error variance battle array
k|k:
Utilize following formula to the estimated value of kth moment generator dynamic state quantity
generate the volume point X of the weights such as a group
i, k|k:
In formula, parameter
n is quantity of state dimension; Exemplary, if generator dynamic state estimator quantity of state dimension n=4, then volume point number is 8.
Utilize following formula to convert each volume point, obtain the predicted value of all volumes point
Summation is weighted to the predicted value of all Generator Status amount volume points, obtains quantity of state predicted value
Exemplary, if generator dynamic state estimator quantity of state dimension n=4, then the volume point weight utilizing sphere-radial generate rule is 1/8.
And obtain prediction error conariance battle array P by following formula
k+1|k:
2, filtering.
When filtering process, become multidimensional observation noise scale factor during introducing and the generator dynamical state forecast result of real-time measurement value to the k+1 moment in conjunction with PMU carries out filtering, thus realize the accurate estimation of generator dynamical state result.
The generator dynamical state forecast result in described k+1 moment comprises: the generator dynamic state quantity predicted value in k+1 moment
with Generator Status amount prediction error conariance battle array P
k+1|k;
To prediction error conariance battle array P
k+1|kcarry out Cholesky decomposition, obtain the On Square-Rooting Matrices S of k+1 moment Generator Status amount prediction error conariance battle array
k+1|k:
Utilize following formula to the generator dynamic state quantity predicted value in k+1 moment
generate the quantity of state predicted value volume point X of the weights such as a group
i, k+1|k:
Utilize following formula to convert each quantity of state predicted value volume point, obtain the volume point Z of PMU measurement amount predicted value
i, k+1|k:
Z
i,k+1|k=H(X
i,k+1|k,u
k);(17)
Summation is weighted to the volume point of all PMU measurement amount predicted values, and then obtains PMU measurement amount predicted value
Calculate PMU measurement amount prediction error variance matrix P
vv, k+1:
Utilize the real-time measurement value z obtaining PMU
k+1=[δ
zk+1ω
zk+1p
ek+1]
tcalculate the vectorial e of new breath
k+1, and then when obtaining, become multidimensional observation noise scale factor γ
k+1:
In formula, M is that the window window of the estimation technique is long;
Recycling following formula calculates diagonal matrix γ '
k+1:
γ′
k+1=diag(γ′
1,γ′
2,…,γ′
m);(22)
In formula, diagonal element γ '
1value be: γ '
i=max{1, γ
k+1, ii, i=1,2 ..., n; γ
k+1, iifor γ
k+1i-th diagonal element;
According to the Cross-covariance P between following formula calculating generator quantity of state predicted value and PMU measurement amount predicted value
xz, k+1|k:
Computer card Kalman Filtering gain W again
k+1:
W
k+1=P
xz,k+1|k(P
vv,k+1+γ′
k+1R
k+1)
-1;(24)
The k+1 moment is utilized newly to cease vectorial e
k+1, and by Kalman filtering gain W
k+1to the generator dynamic state quantity predicted value in k+1 moment
carry out filtering, obtain the estimated value of k+1 moment generator dynamic state quantity
And calculate acquisition generator dynamic state quantity estimation error variance battle array P by following formula
k+1|k+1:
In above-mentioned filtering, will time become multidimensional observation noise scale factor and be incorporated in filtering process, according to measurement new breath, on-line tuning is carried out to error in measurement, makes its approaching to reality noise more; Error calculation filter gain after recycling adjustment, accurately can revise quantity of state predicted value when measurement amount contains bad data, obtain accurate quantity of state estimated value.
The principle that the embodiment of the present invention is improved CKF algorithm is as follows:
Volume Kalman filtering calculates posterior probability density function according to the volume point determined, without the need to the Jacobian matrix of calculation of complex, realizes being more prone to, it also avoid truncation error simultaneously than EKF.Compared with Unscented transform Kalman filtering, the volume point of volume Kalman filtering is symmetrical in the subspace of lower one dimension to be occurred, and each volume point weights equal and opposite in direction, without the need to shifting to an earlier date parameters as Unscented transform Kalman filtering, more simple on selection mode, adaptability is stronger.
As a kind of method for estimating nonlinear state based on Kalman filtering framework, volume Kalman filtering needs based on system model accurately and noise statistics.But, in practical application, obtain noise priori statistical property accurately on the one hand often more difficult; On the other hand, even if obtain accurate priori noise statistics, system is subject to the impact of inside or external environment condition uncertain factor in operational process, may cause the change of noise statistics.Such as, in generator dynamic state estimator process, PMU measures and likely occurs bad data.There is bad data once measure, error in measurement variance matrix R and actual error will be made not to be inconsistent, calculated the deviation that error in measurement population variance can be made correctly cannot to reflect predicted value by classic method, error in measurement population variance P in classic method
zz, k+1|k=P
vv, k+1+ R
k+1; Thus Kalman filtering gain in classic method correctly cannot be revised quantity of state predicted value, make estimated result inaccurate, Kalman filtering gain in classic method
traditional volume Kalman filtering algorithm is more responsive to initial error, measurement noise statistical property is lacked to the adaptive ability of on-line tuning.The factor such as unexpected change, the reduction of observed reading confidence level of noise all may affect algorithm estimated accuracy, even causes filter value to be dispersed.
Therefore, need by changing filter gain to adapt to the change of noise statistics.When noise priori statistical property is more accurate, new breath credibility is high, and corresponding weight ratio is comparatively large, can the effective error that causes of restraint speckle by a larger filter gain; And when bad data appears in PMU, noise statistics deviation is comparatively large, new breath credibility reduces, and needs to be changed by reduction filter gain newly to cease weight, makes adjustment, improve filtering accuracy to state updating.Therefore, can by changing filter gain to adapt to the change of noise statistics.When noise priori statistical property is more accurate, new breath credibility is high, and corresponding weight ratio is comparatively large, can the effective error that causes of restraint speckle by a larger filter gain; And when bad data appears in PMU, noise statistics deviation is comparatively large, new breath credibility reduces, and needs to be changed by reduction filter gain newly to cease weight, makes adjustment, improve filtering accuracy to state updating.
In the embodiment of the present invention, for measuring the problem occurring that bad data causes error in measurement variance and actual value and is not inconsistent, multidimensional observation noise scale factor introducing CKF when one, will be become.From formula (25), CKF filtering algorithm needs according to newly ceasing vector
to the predicted value of quantity of state
revise, and then obtain quantity of state estimated value
correction degree is by newly ceasing vectorial e
k+1with Kalman filtering gain W
k+1common decision.Specifically:
When adopting classic method, when actual measurement noise conforms to given error in measurement variance matrix R, e
k+1and W
k+1correctly can revise predicted value, CKF can obtain estimated result accurately.But, when there is bad data in equivalent measurement, the vectorial e of new breath
k+1the element that middle bad data is corresponding increases suddenly, and W
k+1do not adjust thereupon, inaccurate to the correction of quantity of state predicted value, cause estimated result accuracy to decline.
Formula by described classic method above: P
zz, k+1|k=P
vv, k+1+ R
k+1,
visible, W
k+1the function of measuring noise square difference battle array R, in order to make W
k+1correctly can revise quantity of state predicted value, when constructing in the embodiment of the present invention, become multidimensional observation noise scale factor γ
k+1, on-line tuning is carried out to R.When there is bad data in equivalent measurement, γ
k+1value change thereupon, and then adjustment W
k+1accurately can revise quantity of state predicted value.
That is, by traditional P
zz, k+1|k=P
vv, k+1+ R
k+1expression formula is revised as:
P
zz,k+1|k=P
vv,k+1+γ
k+1R
k+1(26)
In formula, γ
k+1for becoming multidimensional observation noise scale factor when m ties up, m is measurement amount dimension.Thus, traditional expression formula
Become:
W
k+1=P
xz,k+1|k(P
vv,k+1+γ
k+1R
k+1)
-1(27)
In formula, part of inverting on the right of equal sign is observation noise variance, for making W
k+1can accurately revise quantity of state predicted value, demand fulfillment
When equivalent measurement does not have a bad data, γ
k+1for m ties up unit matrix.When bad data appears in equivalent measurement, newly cease e
k+1can increase, γ
k+1be no longer unit matrix, concrete value condition needs to calculate separately.
Utilization is windowed, and estimation technique calculating is new ceases real-time covariance:
In formula, P
e, k+1for newly ceasing covariance; M is that the window window of the estimation technique is long.When measuring existing bad data, P must be had
e, k+1> P
vv, k+1+ R
k+1, then need to calculate γ
k+1, and then error in measurement variance is adjusted.γ
k+1value principle be that to make newly to cease variance equal with observation noise, namely
Solved by following formula (31) and obtain γ
k+1:
It is not diagonal matrix that the new breath covariance using the estimation technique of windowing to calculate becomes multidimensional observation noise scale factor when formula (31) may be caused to calculate, thus makes matrix inversion on the right of formula (27) equal sign occur unusual.For this problem, definition diagonal matrix
γ′
k+1=diag(γ′
1,γ′
2,…,γ′
m)(32)
Diagonal element γ ' in formula (32)
ivalue be
γ′
i=max{1,γ
k+1,ii},i=1,2,…,n(33)
In formula, γ
k+1, iifor γ
k+1i-th diagonal element.Due to γ '
k+1for diagonal matrix and diagonal element is all greater than 0, therefore, by γ '
k+1as time become multidimensional observation noise scale factor calculating formula (27) filter gain just can ensure to invert do not occur unusual.Then filter gain computing formula is
W
k+1=P
xz,k+1|k(P
vv,k+1+γ′
k+1R
k+1)
-1(34)
From formula (34), when introducing in filter gain, become multidimensional observation noise scale factor, can adjust error in measurement variance matrix according to the new breath of measurement.After bad data appears in equivalent measurement, under the effect of observation noise scale factor, filter gain reduces, and then reduces bad data to the impact of estimated result.
In order to further illustrate the present invention, with concrete example, emulation testing is carried out to such scheme more below, specifically:
The emulation testing of IEEE9 node:
Electric system distributed dynamic method for estimating state based on PMU provided by the invention is applied to IEEE9 bus test system (reference capacity of system is 100MVA) as shown in Figure 2, the correlation parameter of IEEE9 bus test system and generator dynamic parameter and branch parameters are as shown in table 1-table 4.
Table 1IEEE9 bus test system node data table
Note: balance node is node 1, except balance node, all nodes clearly not providing voltage magnitude are PQ node, and the point clearly providing voltage magnitude is PV node.
Table 2IEEE9 bus test system circuit branch data table
Table 3IEEE9 bus test system transformer branch tables of data
Table 4IEEE9 bus test system generator dynamic data tables
In this example, be engraved in the circuit generation three-phase metallic short circuit fault between node 5 and 6 when supposing t=0, faulty line is opened from system break by breaker actuation subsequently.Utilize BPA to carry out simulation calculation to failure process, simulation step length is 1 cycle, i.e. 0.02s, and simulation time is 18s.Using simulation result as true value, the basis of true value superposes white Gaussian noise as measuring value.The PMU error in measurement standard deviation of generator amature merit angle and angular rate is respectively 2 ° and 0.1%, and the amplitude of generator outlet voltage phasor and the PMU error in measurement standard deviation of phase angle are respectively 0.1% and 0.1 °.
Fig. 3-Fig. 5 sets forth the dynamic state estimator result of generator G1, G2 and G3 in IEEE9 bus test system, the result that "---robust CKF " is wherein scheme of the present invention, the result that " ... CKF " is traditional C KF algorithm, "-true value " is actual value.
For generator G1, respectively when t=4s and t=8s generator absolute merit angle PMU measuring value in artificial superposition 10 continuous umber of defectives strong points.As can be seen from Figure 3, due to the existence of bad data, there is larger fluctuation in the generator dynamic state estimator based on volume Kalman filtering.This is because error in measurement variance measure with the PMU at absolute merit angle in the actual error of bad data differ comparatively large, the filter gain of volume Kalman filtering correctly cannot be revised quantity of state predicted value, finally causes the estimated value of quantity of state inaccurate.And change multidimensional observation noise scale factor adjusts in real time to error in measurement variance during the solution of the present invention employing, can change with measurement noise.The filter gain calculated thus correctly can revise the predicted value of Generator Status amount.Therefore, when bad data appears in the PMU measuring value of generator's power and angle, quantity of state estimated value accurately still can be obtained.
For generator G2, in the PMU measuring value of generator amature angular rate, artificially superpose 10 continuous bad datas when t=6s and t=12s respectively, superpose single-point bad data when t=14s.As can be seen from Figure 4, bad data is there is in measurement due to rotor angular rate, generator dynamic state estimator based on volume Kalman filtering algorithm there occurs larger fluctuation, and the solution of the present invention still can obtain dynamic state estimator value comparatively accurately.
For generator G3, the simultaneously continuous bad data of artificial superposition in the PMU measuring value at generator amature angular rate and rotor absolute merit angle when t=12s.As seen from Figure 5, even if definitely the PMU measurement of merit angle and angular rate exists continuous multiple spot bad data simultaneously, the solution of the present invention still can obtain state estimation accurately.
Fig. 6 gives PMU measurement amount when there is bad data, measurement amount time become the situation of change of multidimensional observation noise scale factor.Can find out, when there is bad data in PMU measurement amount, measurement amount time become multidimensional observation noise scale factor and can increase suddenly.This just can effectively regulate measurement error variance, and the Kalman filtering gain calculated thus just accurately can revise the predicted value of quantity of state, thus improves precision of state estimation.
Through the above description of the embodiments, those skilled in the art can be well understood to above-described embodiment can by software simulating, and the mode that also can add necessary general hardware platform by software realizes.Based on such understanding, the technical scheme of above-described embodiment can embody with the form of software product, it (can be CD-ROM that this software product can be stored in a non-volatile memory medium, USB flash disk, portable hard drive etc.) in, comprise some instructions and perform method described in each embodiment of the present invention in order to make a computer equipment (can be personal computer, server, or the network equipment etc.).
The above; be only the present invention's preferably embodiment, but protection scope of the present invention is not limited thereto, is anyly familiar with those skilled in the art in the technical scope that the present invention discloses; the change that can expect easily or replacement, all should be encompassed within protection scope of the present invention.Therefore, protection scope of the present invention should be as the criterion with the protection domain of claims.
Claims (4)
1. a generator dynamic state estimator method, is characterized in that, comprising:
Calculate according to the generator dynamic state estimator result of system equation to the k-1 moment, obtain kth moment generator dynamic state estimator result;
CKF filtering algorithm is utilized to forecast and filtering process kth moment generator dynamic state estimator result; Wherein, the generator dynamical state forecast result in k+1 moment is obtained by forecast process; When filtering process, become multidimensional observation noise scale factor during introducing and the generator dynamical state forecast result of real-time measurement value to the k+1 moment in conjunction with PMU carries out filtering, thus realize the accurate estimation of generator dynamical state result.
2. method according to claim 1, is characterized in that, described system equation is:
In formula, the expression k+1 moment that subscript k+1 and k is all corresponding and k moment; F and H is respectively state equation function and measurement equation function, and x, u and z are respectively quantity of state, controlled quentity controlled variable and measurement amount; V and w is respectively system noise and measurement noise, and error covariance matrix is respectively the normal distribution of Q and R;
Wherein, quantity of state x=[δ ω E '
qe '
d]
t, controlled quentity controlled variable
the concrete form of state equation is:
In formula, δ and ω is respectively the absolute merit angle of generator amature and angular rate perunit value;
with
the simplification literary style to time diffusion, the expression d/dt differentiating operator of letter top; T
jfor generator inertia time constant; T
mfor machine torque; U
twith
be respectively amplitude and the phase angle of generator outlet voltage phasor; X '
qwith X '
dbe respectively q axle and d axle transient reactance; E '
qwith E '
dbe respectively q axle and d axle transient state electromotive force; D is ratio of damping; T '
q0with T '
d0be respectively q axle and d axle open circuit transient time-constant; E
ffor stator excitation electromotive force; X
qand X
dbe respectively q axle and d axle synchronous reactance;
Measurement amount z=[δ
zω
zp
e]
t, the concrete form of measurement equation is:
In formula, δ
z, ω
zand P
ebe respectively the PMU measuring value of the absolute merit angle of rotor, angular rate and output electromagnetic power;
According to the concrete form of above-mentioned state equation and measurement equation, and k moment system noise variance matrix Q
kwith k+1 moment measuring noise square difference battle array R
k+1, then generator dynamic state estimator can be realized;
K+1 moment measuring noise square difference battle array R wherein
k+1value is carried out, k moment system noise variance matrix Q according to the actual error in measurement of PMU
kbe expressed as:
In formula,
with
be respectively k moment δ, ω, E '
qwith E '
dsystem noise variance;
Calculated by formula of error transmission:
In formula: σ is noise variance; σ
uwith
be respectively amplitude and the phase angle PMU error in measurement standard deviation of generator outlet voltage phasor;
Then have:
In formula, Δ t is step-length.
3. method according to claim 2, is characterized in that, the step of the described generator dynamical state forecast result by the forecast process acquisition k+1 moment comprises:
Described kth moment generator dynamic state estimator result comprises: the estimated value of kth moment generator dynamic state quantity
with estimation error variance battle array P
k|k;
To P
k|kcarry out Cholesky decomposition, obtain the On Square-Rooting Matrices S of k moment estimation error variance battle array
k|k:
Utilize following formula to the estimated value of kth moment generator dynamic state quantity
generate the volume point X of the weights such as a group
i, k|k:
In formula, parameter
i=1,2 ..., 2n, n are quantity of state dimension;
Utilize following formula to convert each volume point, obtain the predicted value of all volumes point
Summation is weighted to the predicted value of all Generator Status amount volume points, obtains quantity of state predicted value
And obtain prediction error conariance battle array P by following formula
k+1|k:
4. method according to claim 2, it is characterized in that, described when filtering process, become multidimensional observation noise scale factor during introducing and the generator dynamical state forecast result of real-time measurement value to the k+1 moment in conjunction with PMU carries out filtering, thus the accurate estimation realizing generator dynamical state result comprises:
The generator dynamical state forecast result in described k+1 moment comprises: the generator dynamic state quantity predicted value in k+1 moment
with Generator Status amount prediction error conariance battle array P
k+1|k;
To prediction error conariance battle array P
k+1|kcarry out Cholesky decomposition, obtain the On Square-Rooting Matrices S of k+1 moment Generator Status amount prediction error conariance battle array
k+1|k:
Utilize following formula to the generator dynamic state quantity predicted value in k+1 moment
generate the quantity of state predicted value volume point X of the weights such as a group
i, k+1|k:
Utilize following formula to convert each quantity of state predicted value volume point, obtain the volume point Z of PMU measurement amount predicted value
i, k+1|k:
Z
i,k+1|k=H(X
i,k+1|k,u
k);
Summation is weighted to the volume point of all PMU measurement amount predicted values, and then obtains PMU measurement amount predicted value
Calculate PMU measurement amount prediction error variance matrix P
vv, k+1:
Utilize the real-time measurement value z obtaining PMU
k+1=[δ
zk+1ω
zk+1p
ek+1]
tcalculate the vectorial e of new breath
k+1, and then when obtaining, become multidimensional observation noise scale factor γ
k+1:
In formula, M is that the window window of the estimation technique is long;
Recycling following formula calculates diagonal matrix γ '
k+1:
γ′
k+1=diag(γ′
1,γ′
2,…,γ′
m);
In formula, diagonal element γ '
1value be: γ '
i=max{1, γ
k+1, ii, i=1,2 ..., n; γ
k+1, iifor γ
k+1i-th diagonal element;
According to the Cross-covariance P between following formula calculating generator quantity of state predicted value and PMU measurement amount predicted value
xz, k+1|k:
Computer card Kalman Filtering gain W again
k+1:
W
k+1=P
xz,k+1|k(P
vv,k+1+γ′
k+1R
k+1)
-1;
The k+1 moment is utilized newly to cease vectorial e
k+1, and by Kalman filtering gain W
k+1to the generator dynamic state quantity predicted value in k+1 moment
carry out filtering, obtain the estimated value of k+1 moment generator dynamic state quantity
And calculate acquisition generator dynamic state quantity estimation error variance battle array P by following formula
k+1|k+1:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510973903.2A CN105403834B (en) | 2015-12-22 | 2015-12-22 | A kind of generator dynamic state estimator method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510973903.2A CN105403834B (en) | 2015-12-22 | 2015-12-22 | A kind of generator dynamic state estimator method |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105403834A true CN105403834A (en) | 2016-03-16 |
CN105403834B CN105403834B (en) | 2019-04-12 |
Family
ID=55469437
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510973903.2A Active CN105403834B (en) | 2015-12-22 | 2015-12-22 | A kind of generator dynamic state estimator method |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105403834B (en) |
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106844952A (en) * | 2017-01-20 | 2017-06-13 | 河海大学 | Based on the generator dynamic state estimator method without mark Particle filtering theory |
CN107425548A (en) * | 2017-09-11 | 2017-12-01 | 河海大学 | A kind of interpolation H ∞ EKFs generator dynamic state estimator method |
CN107478990A (en) * | 2017-09-11 | 2017-12-15 | 河海大学 | A kind of generator electromechanical transient process method for dynamic estimation |
CN107765179A (en) * | 2017-06-26 | 2018-03-06 | 河海大学 | It is a kind of to be applied to measure the generator dynamic state estimator method lost |
CN109100649A (en) * | 2018-06-25 | 2018-12-28 | 南京南瑞继保电气有限公司 | Parameter estimation method for generator excitation system and speed regulation system based on phasor measurement |
CN109270455A (en) * | 2018-10-24 | 2019-01-25 | 郑州轻工业学院 | Induction machine state monitoring method based on hyposensitiveness Ensemble Kalman Filter |
CN109477869A (en) * | 2016-07-25 | 2019-03-15 | 三菱电机株式会社 | The diagnostic device of motor |
CN109711012A (en) * | 2018-12-14 | 2019-05-03 | 华北电力大学 | A kind of PMU single channel based on singular spectrum analysis loses the restoration methods of data |
CN109918862A (en) * | 2019-04-28 | 2019-06-21 | 河海大学 | A kind of generator method for dynamic estimation filtered based on robust without mark H infinity |
CN111948534A (en) * | 2020-07-31 | 2020-11-17 | 华北电力科学研究院有限责任公司 | Generator state early warning method and system |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104777426A (en) * | 2015-04-17 | 2015-07-15 | 河海大学 | Power generator dynamic state estimation method based on unscented transformation strong tracking filtering |
-
2015
- 2015-12-22 CN CN201510973903.2A patent/CN105403834B/en active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104777426A (en) * | 2015-04-17 | 2015-07-15 | 河海大学 | Power generator dynamic state estimation method based on unscented transformation strong tracking filtering |
Non-Patent Citations (3)
Title |
---|
郑维荣等: "基于扩展集员滤波的发电机动态状态估计", 《电力科学与工程》 * |
陈亮: "电力***分布式动态状态估计研究", 《中国博士学位论文全文数据库工程科技II辑》 * |
陈亮等: "基于容积卡尔曼滤波的发电机动态状态估计", 《中国电机工程学报》 * |
Cited By (16)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109477869A (en) * | 2016-07-25 | 2019-03-15 | 三菱电机株式会社 | The diagnostic device of motor |
CN109477869B (en) * | 2016-07-25 | 2021-02-26 | 三菱电机株式会社 | Diagnostic device for motor |
CN106844952A (en) * | 2017-01-20 | 2017-06-13 | 河海大学 | Based on the generator dynamic state estimator method without mark Particle filtering theory |
CN107765179A (en) * | 2017-06-26 | 2018-03-06 | 河海大学 | It is a kind of to be applied to measure the generator dynamic state estimator method lost |
CN107478990A (en) * | 2017-09-11 | 2017-12-15 | 河海大学 | A kind of generator electromechanical transient process method for dynamic estimation |
CN107478990B (en) * | 2017-09-11 | 2019-11-12 | 河海大学 | A kind of generator electromechanical transient process method for dynamic estimation |
CN107425548B (en) * | 2017-09-11 | 2020-06-16 | 河海大学 | Interpolation H∞Dynamic state estimation method for extended Kalman filter generator |
CN107425548A (en) * | 2017-09-11 | 2017-12-01 | 河海大学 | A kind of interpolation H ∞ EKFs generator dynamic state estimator method |
CN109100649A (en) * | 2018-06-25 | 2018-12-28 | 南京南瑞继保电气有限公司 | Parameter estimation method for generator excitation system and speed regulation system based on phasor measurement |
CN109100649B (en) * | 2018-06-25 | 2020-10-16 | 南京南瑞继保电气有限公司 | Parameter estimation method for generator excitation system and speed regulation system based on phasor measurement |
CN109270455A (en) * | 2018-10-24 | 2019-01-25 | 郑州轻工业学院 | Induction machine state monitoring method based on hyposensitiveness Ensemble Kalman Filter |
CN109270455B (en) * | 2018-10-24 | 2021-02-02 | 郑州轻工业学院 | Induction motor state monitoring method based on weak-sensitivity ensemble Kalman filtering |
CN109711012A (en) * | 2018-12-14 | 2019-05-03 | 华北电力大学 | A kind of PMU single channel based on singular spectrum analysis loses the restoration methods of data |
CN109918862A (en) * | 2019-04-28 | 2019-06-21 | 河海大学 | A kind of generator method for dynamic estimation filtered based on robust without mark H infinity |
CN111948534A (en) * | 2020-07-31 | 2020-11-17 | 华北电力科学研究院有限责任公司 | Generator state early warning method and system |
CN111948534B (en) * | 2020-07-31 | 2023-05-05 | 华北电力科学研究院有限责任公司 | Generator state early warning method and system |
Also Published As
Publication number | Publication date |
---|---|
CN105403834B (en) | 2019-04-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105403834A (en) | Dynamic state evaluation method of generator | |
CN107577870B (en) | Power distribution network voltage power sensitivity robust estimation method based on synchronous phasor measurement | |
CN101404412B (en) | Method for static electric voltage stability analysis | |
CN106844952A (en) | Based on the generator dynamic state estimator method without mark Particle filtering theory | |
CN107590317A (en) | A kind of generator method for dynamic estimation of meter and model parameter uncertainty | |
CN110543720B (en) | State estimation method based on SDAE-ELM pseudo-measurement model | |
CN103474992B (en) | Real-time on-line identification criterion of electric system node voltage steady state | |
US10884060B1 (en) | Dynamic parameter estimation of generators | |
CN103972884A (en) | Electric system state estimation method | |
CN105184027B (en) | A kind of power load modelling approach based on interacting multiple model algorithm | |
CN107132772A (en) | Real-time simulation system and method for AC/DC power grid | |
CN107658881A (en) | Voltage stability critical point determination methods based on Thevenin's equivalence method | |
CN105375484A (en) | PMU-based electric power system distributed dynamic-state estimation method | |
CN104156542A (en) | Implicit-projection-based method for simulating stability of active power distribution system | |
Zacharia et al. | Measurement errors and delays on wide-area control based on IEEE Std C37. 118.1-2011: impact and compensation | |
CN102280884B (en) | Power grid equivalence method | |
Chouhan et al. | A literature review on optimal placement of PMU and voltage stability | |
CN112511056A (en) | Robust generator dynamic state estimation method based on phasor measurement | |
Bila | Power system dynamic state estimation and load modeling | |
CN108054768A (en) | Transient stability evaluation in power system method based on principal component analysis | |
CN109638811B (en) | Power distribution network voltage power sensitivity robust estimation method based on model equivalence | |
CN107765179A (en) | It is a kind of to be applied to measure the generator dynamic state estimator method lost | |
CN111585285A (en) | Load modeling method and device based on voltage monitoring system | |
CN104779613B (en) | Test-based equivalent modeling method for electric element comprising converter | |
Xie et al. | Distributed quasi-dynamic state estimation with both GPS-synchronized and non-synchronized data |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |