CN105403834B - A kind of generator dynamic state estimator method - Google Patents
A kind of generator dynamic state estimator method Download PDFInfo
- Publication number
- CN105403834B CN105403834B CN201510973903.2A CN201510973903A CN105403834B CN 105403834 B CN105403834 B CN 105403834B CN 201510973903 A CN201510973903 A CN 201510973903A CN 105403834 B CN105403834 B CN 105403834B
- Authority
- CN
- China
- Prior art keywords
- generator
- moment
- state
- measurement
- pmu
- 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.)
- Active
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 kind of generator dynamic state estimator method, time-varying multidimensional observation noise scale factor is introduced into CKF dynamic state estimator by this method, is carried out on-line tuning to error in measurement according to new breath is measured, is made its more approaching to reality noise.Error calculation filtering gain adjusted is recycled, quantity of state predicted value can accurately be corrected in the case where measurement contains bad data, accurate quantity of state estimated value is obtained；Meanwhile giving the specific configuration method of time-varying multidimensional observation noise scale factor；And invert for filtering gain and unusual problem occurs, propose solution；It, being capable of influence of the effective inhibitory amount survey bad data to generator dynamic state estimator under the premise of guaranteeing real-time demand by using method disclosed by the invention.
Description
Technical field
The present invention relates to Electrical Power System Dynamic State Estimation field more particularly to a kind of generator dynamic state estimators
Method.
Background technique
The appearance of synchronous phasor measurement unit (phasor measurement unit, PMU) is electric power system transient stability
Analysis provides new technological means with control.However, when interference, measuring device failure, synchronization signal are lost, it is past
Toward the appearance for leading to bad data.Bad data, which may make PMU in application process, leads to the analysis result and control strategy of mistake.
State estimation can reject bad data present in measurement, therefore, study in the electric system electromechanics transient process based on PMU
Dynamic state estimator is most important.
In recent years, the electromechanical transient process generator dynamic state estimator problem based on PMU is Electrical Power System Dynamic state
The hot spot in estimation field.For with nonlinear generator dynamical equation, dynamic state estimator problem is one typical
Nonlinear problem, using be basic scheme with Kalman filtering algorithm is a kind of universal resolving ideas, is such as based on expansion card
The dynamic state estimator of Kalman Filtering (EKF), Unscented kalman filtering (UKF).The linearization procedure of EKF leads to truncated error mistake
Greatly, UKF is it needs to be determined that parameter value, flexibility is poor, using inconvenience.For this problem, Canadian scholar Simon Haykin in
It is proposed volume Kalman filtering (CKF) algorithm in 2009.However, no matter UKF or CKF, equivalent measurement estimates there are when bad data
Meter precision can all be affected to a certain extent, even result in estimated result substantial deviation actual value, and dynamic state estimator loses
It loses.
Summary of the invention
The object of the present invention is to provide a kind of generator dynamic state estimator method, this method can be measured in PMU
Amount obtains accurate quantity of state estimated value containing accurately being corrected in the case where bad data to Generator Status amount predicted value.
The purpose of the present invention is what is be achieved through the following technical solutions:
A kind of generator dynamic state estimator method, comprising:
It is calculated according to generator dynamic state estimator result of the system equation to the k-1 moment, obtains the power generation of kth moment
Motor-driven state state estimation result；
Kth moment generator dynamic state estimator result is forecast and is filtered using CKF filtering algorithm；Its
In, the generator dynamical state forecast result for obtaining the k+1 moment is handled by forecast；In filtering processing, time-varying multidimensional is introduced
Observation noise scale factor simultaneously combines the real-time amount measured value of PMU to filter the generator dynamical state forecast result at k+1 moment
Wave, to realize the accurate estimation of generator dynamical state result.
Further, the system equation are as follows:
In formula, subscript k+1 and k corresponding expression k+1 moment and k moment；F and H is respectively state equation function and amount
Equation functions are surveyed, x, u and z are respectively quantity of state, control amount and measurement；V and w is respectively system noise and measures noise, accidentally
Poor variance matrix is respectively the normal distribution of Q and R；
Wherein, quantity of state x=[δ ω E 'qE′d]T, control amountThe concrete form of state equation are as follows:
In formula, δ is the absolute generator rotor angle per unit value of generator amature, and ω is angular rate per unit value；WithIt is to time diffusion
Simplification literary style, letter top expression d/dt differential operator；TJFor generator inertia time constant；TmFor machine torque；Ut
WithThe respectively amplitude and phase angle of generator outlet voltage phasor；X′qWith X 'dRespectively q axis and d axis transient reactance；E′qWith
E′dRespectively q axis and d axis transient state electromotive force；D is damped coefficient；T′q0With T 'd0Respectively q axis and d axis open circuit transition time is normal
Number；EfFor stator excitation electromotive force；XqAnd XdRespectively q axis and d axis synchronous reactance；
Measurement z=[δz ωz Pe]T, the concrete form of measurement equation are as follows:
In formula, δz、ωzAnd PeThe respectively absolute generator rotor angle of rotor, angular rate and the PMU measuring value for exporting electromagnetic power；
According to the concrete form and k moment system noise variance matrix Q of above-mentioned state equation and measurement equationkAnd k+1
Moment measuring noise square difference battle array Rk+1, then can be realized generator dynamic state estimator；
K+1 moment measuring noise square difference battle array R thereink+1Value, k moment system noise are carried out according to the practical error in measurement of PMU
Sound variance matrix QkIt indicates are as follows:
In formula,WithRespectively k moment δ, ω, E 'qWith E 'dSystem noise variance；
It is calculated by formula of error transmission:
In formula: σ is noise variance；For the amplitude PMU error in measurement standard deviation of generator outlet voltage phasor,For power generation
The phase angle PMU error in measurement standard deviation of machine exit potential phasor；
Then have:
In formula, Δ t is step-length.
Further, described the step of handling the generator dynamical state forecast result for obtaining the k+1 moment by forecast, wraps
It includes:
The kth moment generator dynamic state estimator result includes: the estimation of kth moment generator dynamic state quantity
ValueWith estimation error variance battle array Pk|k；
To Pk|kCholesky decomposition is carried out, the On Square-Rooting Matrices S of k moment estimation error variance battle array is obtainedk|k:
Using following formula to the estimated value of kth moment generator dynamic state quantityGenerate the volume point of the weights such as one group
Xi,k|k:
In formula, parameterN is quantity of state dimension；
Each volume point is converted using following formula, obtains the predicted value of all volume points
Summation is weighted to the predicted value of all Generator Status amount volume points, obtains quantity of state predicted value
And prediction error conariance battle array P is obtained by following formulak+1|k:
Further, described in filtering processing, introduce time-varying multidimensional observation noise scale factor and in conjunction with the real-time of PMU
Measuring value is filtered the generator dynamical state forecast result at k+1 moment, to realize generator dynamical state result
Accurate estimation includes:
The generator dynamical state forecast result at the k+1 moment includes: that the generator dynamic state quantity at k+1 moment is pre-
Report valueWith Generator Status amount prediction error conariance battle array Pk+1|k；
To prediction error conariance battle array Pk+1|kCholesky decomposition is carried out, the forecast of k+1 moment Generator Status amount is obtained and misses
The On Square-Rooting Matrices S of poor covariance matrixk+1|k:
Using following formula to the generator dynamic state quantity predicted value at k+1 momentThe quantity of state for generating the weights such as one group is pre-
Report value volume point Xi,k+1|k:
Each quantity of state predicted value volume point is converted using following formula, obtains the volume of PMU measurement predicted value
Point Zi,k+1|k:
Zi,k+1|k=H (Xi,k+1|k,uk)；
Weighted sum is carried out to the volume point of all PMU measurement predicted values, and then obtains PMU measurement predicted value
Calculate PMU measurement prediction error variance matrix Pvv,k+1:
Using obtaining the real-time amount measured value z of PMUk+1=[δzk+1 ωzk+1 Pek+1]TCalculate new breath vector ek+1, and then obtain
Time-varying multidimensional observation noise scale factor γk+1:
In formula, M is that the window of the windowing estimation technique is long；
Following formula is recycled to calculate diagonal matrix γ 'k+1:
γ′k+1=diag (γ '1,γ′2,…,γ′m)；
In formula, diagonal element γ '1Value are as follows: γ 'i=max { 1, γk+1,ii, i=1,2 ..., n；γk+1,iiFor γk+1
I-th of diagonal element；
The Cross-covariance between Generator Status amount predicted value and PMU measurement predicted value is calculated according to the following formula
Pxz,k+1|k:
Kalman filtering gain W is calculated againk+1:
Wk+1=Pxz,k+1|k(Pvv,k+1+γ′k+1Rk+1)-1；
Vector e is newly ceased using the k+1 momentk+1, and pass through Kalman filtering gain Wk+1To the generator dynamic shape at k+1 moment
State amount predicted valueIt is filtered, obtains the estimated value of k+1 moment generator dynamic state quantity
And it is calculate by the following formula and obtains generator dynamic state quantity estimation error variance battle array Pk+1|k+1:
As seen from the above technical solution provided by the invention, which leads for bad data in PMU measurement
The problem of causing generator dynamic state estimator result to deviate true value, introduces CKF for time-varying multidimensional observation noise scale factor
In filtering algorithm, the robustness of this method is not only increased, accurate quantity of state estimated value also can be obtained.
Detailed description of the invention
In order to illustrate the technical solution of the embodiments of the present invention more clearly, required use in being described below to embodiment
Attached drawing be briefly described, it should be apparent that, drawings in the following description are only some embodiments of the invention, for this
For the those of ordinary skill in field, without creative efforts, it can also be obtained according to these attached drawings other
Attached drawing.
Fig. 1 is a kind of outline flowchart of generator dynamic state estimator method provided in an embodiment of the present invention；
Fig. 2 is IEEE9 bus test system schematic diagram provided in an embodiment of the present invention；
Fig. 3 is generator G1 dynamic state estimator result figure in IEEE9 node emulation testing provided in an embodiment of the present invention；
Fig. 4 is generator G2 dynamic state estimator result figure in IEEE9 node emulation testing provided in an embodiment of the present invention；
Fig. 5 is generator G3 dynamic state estimator result figure in IEEE9 node emulation testing provided in an embodiment of the present invention；
Fig. 6 is that multidimensional time-varying observation noise scale factor becomes in IEEE9 node emulation testing provided in an embodiment of the present invention
Change situation map.
Specific embodiment
With reference to the attached drawing in the embodiment of the present invention, technical solution in the embodiment of the present invention carries out clear, complete
Ground description, it is clear that described embodiments are only a part of the embodiments of the present invention, instead of all the embodiments.Based on this
The embodiment of invention, every other implementation obtained by those of ordinary skill in the art without making creative efforts
Example, belongs to protection scope of the present invention.
The embodiment of the present invention provides a kind of generator dynamic state estimator method, mainly includes the following steps:
It is calculated according to generator dynamic state estimator result of the system equation to the k-1 moment, obtains the power generation of kth moment
Motor-driven state state estimation result；
Kth moment generator dynamic state estimator result is forecast and is filtered using CKF filtering algorithm；Its
In, the generator dynamical state forecast result for obtaining the k+1 moment is handled by forecast；In filtering processing, time-varying multidimensional is introduced
Observation noise scale factor simultaneously combines the real-time amount measured value of PMU to filter the generator dynamical state forecast result at k+1 moment
Wave, to realize the accurate estimation of generator dynamical state result.
The outline flowchart of generator dynamic state estimator is as shown in Figure 1.Generator port voltage can be directly direct by PMU
It measures, is decoupled with external network, it is not necessary to network topology structure is considered, to carry out distributed generator dynamic state estimator.Hair
The huge inertia in inherence of rotor makes generator amature generator rotor angle and angular rate not to dash forward in electromechanical transient process
Become, and meets the constraint condition of the generator equation of motion.PMU can directly measure generator's power and angle and angular speed,
Utilize the synchronous measure value of the available electromagnetic power of GPS；Previous moment quantity of state combination generator port voltage phasor, passes through
The estimated result of system equation acquisition subsequent time Generator Status amount；Then, the predicted value of quantity of state is obtained by forecast, then
In conjunction with PMU to the measuring value of quantity of state, there is robustness CKF algorithm by improved, finally obtain subsequent time generator
Dynamic state estimator result.
In order to make it easy to understand, being described in detail below to the present invention.
Generator dynamic state estimator problem is Nonlinear Filtering Problem, and the system equation of nonlinear dynamic system can be with table
It is shown as:
In formula, subscript k+1 and k corresponding expression k+1 moment and k moment；F and H is respectively state equation function and amount
Equation functions are surveyed, x, u and z are respectively quantity of state, control amount and measurement；V and w is respectively system noise and measures noise, accidentally
Poor variance matrix is respectively the normal distribution of Q and R；
Be difficult to make timely reaction to secondary transient process since existing control is standby, ignore time transient state, dynamical equation by
Six ranks are reduced to quadravalence；
Quantity of state x=[δ ω E 'qE′d]T, control amountThe concrete form of state equation are as follows:
In formula, δ is the absolute generator rotor angle per unit value of generator amature, and ω is angular rate per unit value；WithIt is to time diffusion
Simplification literary style, letter top expression d/dt differential operator；TJFor generator inertia time constant；TmFor machine torque；Ut
WithThe respectively amplitude and phase angle of generator outlet voltage phasor；X′qWith X 'dRespectively q axis and d axis transient reactance；E′qWith
E′dRespectively q axis and d axis transient state electromotive force；D is damped coefficient；T′q0With T 'd0Respectively q axis and d axis open circuit transition time is normal
Number；EfFor stator excitation electromotive force；XqAnd XdRespectively q axis and d axis synchronous reactance；
Measurement z=[δz ωz Pe]T, the concrete form of measurement equation are as follows:
In formula, δz、ωzAnd PeThe respectively absolute generator rotor angle of rotor, angular rate and the PMU measuring value for exporting electromagnetic power；
According to the concrete form and k moment system noise variance matrix Q of above-mentioned state equation and measurement equationkAnd k+1
Moment measuring noise square difference battle array Rk+1, then can be realized generator dynamic state estimator.
K+1 moment measuring noise square difference battle array R thereink+1Value is carried out according to the practical error in measurement of PMU.
For system noise, in the system equation of generator dynamic state estimator, it is assumed that model parameter is accurate；Quantity of state
It is all made of estimated value, precision is higher；T in control amount um、EfIt assume that as steady state value；UtWithUsing PMU measuring value, exist
Error in measurement.Therefore, system noise is mostly derived from UtWithError in measurement.The error in measurement is transmitted by system equation,
Eventually lead to the generation of system noise.
K moment system noise variance matrix QkIt indicates are as follows:
In formula,WithRespectively k moment δ, ω, E 'qWith E 'dSystem noise variance；
It is calculated by formula of error transmission:
In formula: σ is noise variance；σUFor the amplitude PMU error in measurement standard deviation of generator outlet voltage phasor,For hair
The phase angle PMU error in measurement standard deviation of motor exit potential phasor；Its settable value is respectively 0.2% and 0.2 °.
According to formula (5), partial derivative is asked respectively about generator outlet voltage magnitude and phase angle to system equation, and then derive
K moment quantity of state ω, E 'qWith E 'dThe calculation formula of process-noise variance.Since the predicted value of generator's power and angle δ needs to use electricity
The estimated value of angular velocity omega, and ω hair estimate error is relatively small, therefore, the system noise error variance of generator rotor angle is set to one
A lesser positive number, is selected as 10 herein-6。
Then have:
Q can be calculated by formula (7)-(9)kMiddle each element value.In formula, Δ t is step-length, U in formula (7)-(9)tWithValue
For the measuring value of k moment PMU；E′q、E′dValue with δ is the estimated value at k moment.
CKF filtering algorithm is recycled to forecast and be filtered kth moment generator dynamic state estimator result.
1, forecasting process.
The generator dynamic shape that forecast processing obtains the k+1 moment is carried out to kth moment generator dynamic state estimator result
State forecast result；
The kth moment generator dynamic state estimator result includes: the estimation of kth moment generator dynamic state quantity
ValueWith estimation error variance battle array Pk|k；
To Pk|kCholesky decomposition is carried out, the On Square-Rooting Matrices S of k moment estimation error variance battle array is obtainedk|k:
Using following formula to the estimated value of kth moment generator dynamic state quantityGenerate the volume point of the weights such as one group
Xi,k|k:
In formula, parameterN is quantity of state dimension；Illustratively, if generator dynamical state is estimated
Quantity of state dimension n=4 are counted, then volume point number is 8.
Each volume point is converted using following formula, obtains the predicted value of all volume points
Summation is weighted to the predicted value of all Generator Status amount volume points, obtains quantity of state predicted value
Illustratively, it if generator dynamic state estimator quantity of state dimension n=4, is generated using spherical surface-radial direction rule
Volume point weight is 1/8.
And prediction error conariance battle array P is obtained by following formulak+1|k:
2, filtering.
In filtering processing, when introducing time-varying multidimensional observation noise scale factor and combining the real-time amount measured value of PMU to k+1
The generator dynamical state forecast result at quarter is filtered, to realize the accurate estimation of generator dynamical state result.
The generator dynamical state forecast result at the k+1 moment includes: that the generator dynamic state quantity at k+1 moment is pre-
Report valueGenerator Status amount prediction error conariance battle array Pk+1|k；
To prediction error conariance battle array Pk+1|kCholesky decomposition is carried out, obtains the forecast of k+1 moment Generator Status amount and misses
The On Square-Rooting Matrices S of poor covariance matrixk+1|k:
Using following formula to the generator dynamic state quantity predicted value at k+1 momentThe quantity of state for generating the weights such as one group is pre-
Report value volume point Xi,k+1|k:
Each quantity of state predicted value volume point is converted using following formula, obtains the volume of PMU measurement predicted value
Point Zi,k+1|k:
Zi,k+1|k=H (Xi,k+1|k,uk)； (17)
Weighted sum is carried out to the volume point of all PMU measurement predicted values, and then obtains PMU measurement predicted value
Calculate PMU measurement prediction error variance matrix Pvv,k+1:
Using obtaining the real-time amount measured value z of PMUk+1=[δzk+1 ωzk+1 Pek+1]TCalculate new breath vector ek+1, and then obtain
Time-varying multidimensional observation noise scale factor γk+1:
In formula, M is that the window of the windowing estimation technique is long；
Following formula is recycled to calculate diagonal matrix γ 'k+1:
γ′k+1=diag (γ '1,γ′2,…,γ′m)； (22)
In formula, diagonal element γ '1Value are as follows: γ 'i=max { 1, γk+1,ii, i=1,2 ..., n；γk+1,iiFor γk+1
I-th of diagonal element；
The Cross-covariance between Generator Status amount predicted value and PMU measurement predicted value is calculated according to the following formula
Pxz,k+1|k:
Kalman filtering gain W is calculated againk+1:
Wk+1=Pxz,k+1|k(Pvv,k+1+γ′k+1Rk+1)-1； (24)
Vector e is newly ceased using the k+1 momentk+1, and pass through Kalman filtering gain Wk+1To the generator dynamic shape at k+1 moment
State amount predicted valueIt is filtered, obtains the estimated value of k+1 moment generator dynamic state quantity
And it is calculate by the following formula and obtains generator dynamic state quantity estimation error variance battle array Pk+1|k+1:
In above-mentioned filtering, time-varying multidimensional observation noise scale factor is introduced into filtering processing, it is new according to measuring
Breath carries out on-line tuning to error in measurement, makes its more approaching to reality noise；Error calculation filtering gain adjusted is recycled,
Quantity of state predicted value can accurately be corrected in the case where measurement contains bad data, obtain accurate quantity of state
Estimated value.
The principle that the embodiment of the present invention improves CKF algorithm is as follows:
Volume Kalman filtering is to calculate posterior probability density function according to determining volume point, without calculating complexity
Jacobian matrix is more easier than Extended Kalman filter realization, while also avoiding truncated error.With Unscented transform Kalman
Filtering is compared, and the volume point of volume Kalman filtering symmetrically occurs in lower one-dimensional subspace, and each volume point weight is big
Small equal, simpler on selection mode it is not necessary that parameter is arranged in advance as Unscented transform Kalman filtering, adaptability is more
By force.
As a kind of method for estimating nonlinear state based on Kalman filtering framework, volume Kalman filtering is needed with standard
Based on true system model and noise statistics.However, on the one hand obtaining accurate noise priori statistics in practical application
Characteristic is often relatively difficult；On the other hand, even if obtaining accurate priori noise statistics, system in the process of running by
To the influence of internal or external environment uncertain factor, the variation of noise statistics may be caused.For example, motor-driven generating electricity
In state state estimation procedure, PMU measurement is possible to bad data occur.Once bad data occurs in measurement, will to measure
Error covariance matrix R is not inconsistent with actual error, by conventional method come the meeting of calculating so that error in measurement population variance can not correctly reflect
The deviation of predicted value, error in measurement population variance P in conventional methodzz,k+1|k=Pvv,k+1+Rk+1；So that karr in conventional method
Graceful filtering gain can not correctly correct quantity of state predicted value, so that estimated result is inaccurate, Kalman in conventional method
Filtering gainTraditional volume Kalman filtering algorithm is more sensitive to initial error, for amount
Survey the adaptive ability that noise statistics lack on-line tuning.The factors such as suddenly change, the reduction of observation confidence level of noise
Algorithm estimated accuracy may be all influenced, filter value diverging is even resulted in.
Therefore, it is necessary to by changing filtering gain to adapt to the variation of noise statistics.Noise priori statistical property ratio
When more accurate, the new credibility that ceases is high, and corresponding weight is bigger, can effectively inhibit noise by a biggish filtering gain
Caused error；And when bad data occurs in PMU, noise statistics deviation is larger, and the new credibility that ceases reduces, and needs to lead to
It crosses reduction filtering gain and changes new breath weight, adjusting is made to state update, improves filtering accuracy.Therefore, change can be passed through
Filtering gain is to adapt to the variations of noise statistics.When noise priori statistical property is more accurate, the new credibility that ceases is high, right
The weight answered is bigger, can effectively inhibit error caused by noise by a biggish filtering gain；And when PMU appearance is bad
When data, noise statistics deviation is larger, and the new credibility that ceases reduces, and needs to change new breath power by reducing filtering gain
Weight makes adjusting to state update, improves filtering accuracy.
In the embodiment of the present invention, aiming at the problem that bad data, which occurs, in measurement causes error in measurement variance and actual value not to be inconsistent,
One time-varying multidimensional observation noise scale factor is introduced into CKF.From formula (25) as it can be seen that CKF filtering algorithm need according to newly cease to
AmountTo the predicted value of quantity of stateIt is modified, and then obtains quantity of state estimated valueCorrect journey
Degree is by newly ceasing vector ek+1With Kalman filtering gain Wk+1It codetermines.Specifically:
When using conventional method, when actually measurement noise is consistent with given measurement variance matrix R, ek+1And Wk+1Energy
Enough correctly to be corrected to predicted value, CKF can obtain accurate estimated result.However, when bad data occurs in equivalent measurement, newly
Cease vector ek+1The corresponding element of middle bad data increases suddenly, and Wk+1It is not adjusted therewith, quantity of state predicted value is repaired
Positive inaccuracy, causes estimated result accuracy to decline.
By the formula in conventional method described previously:
As it can be seen that Wk+1It is the function of measuring noise square difference battle array R, in order to make Wk+1Quantity of state predicted value can correctly be corrected, this hair
Time-varying multidimensional observation noise scale factor γ is constructed in bright embodimentk+1, on-line tuning is carried out to R.There is bad data in equivalent measurement
When, γk+1Value change therewith, and then adjust Wk+1Quantity of state predicted value can accurately be corrected.
That is, by traditional Pzz,k+1|k=Pvv,k+1+Rk+1Expression formula modification are as follows:
Pzz,k+1|k=Pvv,k+1+γk+1Rk+1 (26)
In formula, γk+1Time-varying multidimensional observation noise scale factor is tieed up for m, m is measurement dimension.Traditional expression as a result,
FormulaBecome:
In formula, right side of the equal sign inverts part for observation noise variance, to make Wk+1Quantity of state predicted value can be accurately corrected,
It needs to meet
When equivalent measurement does not have bad data, γk+1Unit matrix is tieed up for m.When bad data occurs in equivalent measurement, e is newly ceasedk+1Meeting
Increase, γk+1It is no longer unit matrix, specific value condition needs separately to calculate.
Real-time covariance is newly ceased using windowing estimation technique calculating:
In formula, Pe,k+1Newly to cease covariance；M is that the window of the windowing estimation technique is long.When measuring existing bad data, necessarily have
Pe,k+1> Pvv,k+1+Rk+1, then need to calculate γk+1, and then error in measurement variance is adjusted.γk+1Value principle be to make
It is equal with observation noise to obtain newly breath variance, i.e.,
It is solved by following formula (31) and obtains γk+1:
It is observed using the time-varying multidimensional that the new breath covariance that the windowing estimation technique is calculated may result in formula (31) calculating
The Noise Criterion factor is not diagonal matrix, so that the generation of formula (27) right side of the equal sign matrix inversion is unusual.For this problem, it defines
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 of diagonal element.Due to γ 'k+1For diagonal matrix and diagonal element is all larger than 0,
Therefore, by γ 'k+1Filtering gain as time-varying multidimensional observation noise scale factor calculating formula (27), which ensures that, inverts not
Occur unusual.Then filtering gain calculation formula is
Wk+1=Pxz,k+1|k(Pvv,k+1+γ′k+1Rk+1)-1 (34)
By formula (34) as it can be seen that introducing time-varying multidimensional observation noise scale factor in filtering gain, can newly be ceased according to measuring
Error in measurement variance matrix is adjusted.After there is bad data in equivalent measurement, under the action of observation noise scale factor, filter
Wave gain reduces, and then reduces influence of the bad data to estimated result.
In order to further illustrate the present invention, emulation testing is carried out to above scheme with specific example again below, it is specific next
It says:
IEEE9 node emulation testing:
Electric system distributed dynamic method for estimating state provided by the invention based on PMU is applied to as shown in Figure 2
IEEE9 bus test system (reference capacity of system be 100MVA), the relevant parameter of IEEE9 bus test system and power generation
Machine dynamic parameter and branch parameters are as shown in table 1- table 4.
1 IEEE of table, 9 bus test system node data table
Note: balance nodes are node 1, all to be not known that provide the node of voltage magnitude be PQ node in addition to balance nodes,
The point for clearly providing voltage magnitude is PV node.
2 IEEE of table, 9 bus test system circuit branch road tables of data
3 IEEE of table, 9 bus test system transformer branch tables of data
4 IEEE of table, 9 bus test system generator dynamic data tables
In this example, it is assumed that three-phase metallic short circuit failure occurs for route of the t=0 moment between node 5 and 6, then
Breaker actuation opens faulty line from system break.Simulation calculation is carried out to failure process using BPA, simulation step length is 1 week
Wave, i.e. 0.02s, simulation time 18s.Using simulation result as true value, white Gaussian noise conduct is superimposed on the basis of true value
Measuring value.Generator amature generator rotor angle and the PMU error in measurement standard deviation of angular rate are respectively 2 ° and 0.1%, generator outlet
The amplitude of voltage phasor and the PMU error in measurement standard deviation of phase angle are respectively 0.1% and 0.1 °.
The dynamic state estimator knot of generator G1, G2 and G3 in IEEE9 bus test system is set forth in Fig. 3-Fig. 5
Fruit, " --- robust CKF " therein be scheme of the present invention as a result, " ... CKF " be tradition CKF algorithm as a result, "-
True value " is actual value.
For generator G1, artificially it is superimposed in the PMU measuring value of the absolute generator rotor angle of generator in t=4s and t=8s respectively
10 continuous umber of defectives strong points.From figure 3, it can be seen that due to the presence of bad data, the power generation based on volume Kalman filtering
There is biggish fluctuation in machine dynamic state estimator.This is because bad in the PMU of error in measurement variance and absolute generator rotor angle measurement
The actual error difference of data is larger, so that the filtering gain of volume Kalman filtering can not carry out correctly quantity of state predicted value
Amendment eventually leads to the estimated value inaccuracy of quantity of state.And the solution of the present invention uses time-varying multidimensional observation noise scale factor
Error in measurement variance is adjusted in real time, can be changed with noise is measured.Thus the filtering gain energy being calculated
The predicted value of enough correct amendment Generator Status amounts.Therefore, when bad data occurs in the PMU measuring value of generator's power and angle, according to
Accurate quantity of state estimated value can so be obtained.
For generator G2, respectively in t=6s and t=12s into the PMU measuring value of generator amature angular rate people
To be superimposed 10 continuous bad datas, single-point bad data is superimposed in t=14s.From fig. 4, it can be seen that due to rotor angular rate
Measurement in there are bad data, biggish wave has occurred in the generator dynamic state estimator based on volume Kalman filtering algorithm
It is dynamic, and the solution of the present invention can still obtain accurate dynamic state estimator value.
For generator G3, in t=12s simultaneously in the PMU amount of generator amature angular rate and the absolute generator rotor angle of rotor
Continuous bad data is artificially superimposed in measured value.From figure 5 it can be seen that even if the absolutely PMU of generator rotor angle and angular rate measurement exists simultaneously
Continuous multiple spot bad data, the solution of the present invention can still obtain accurate state estimation.
Fig. 6 gives PMU measurement there are when bad data, the change of the time-varying multidimensional observation noise scale factor of measurement
Change situation.As can be seen that when bad data occurs in PMU measurement, the time-varying multidimensional observation noise scale factor meeting of measurement
Increase suddenly.This can effectively adjust measurement error variance, and thus calculated Kalman filtering gain can be repaired accurately
The predicted value of positive status amount, to improve precision of state estimation.
Through the above description of the embodiments, those skilled in the art can be understood that above-described embodiment can
The mode of necessary general hardware platform can also be added to realize by software by software realization.Based on this understanding,
The technical solution of above-described embodiment can be embodied in the form of software products, which can store non-easy at one
In the property lost storage medium (can be CD-ROM, USB flash disk, mobile hard disk etc.), including some instructions are with so that a computer is set
Standby (can be personal computer, server or the network equipment etc.) executes method described in each embodiment of the present invention.
The foregoing is only a preferred embodiment of the present invention, but scope of protection of the present invention is not limited thereto,
Within the technical scope of the present disclosure, any changes or substitutions that can be easily thought of by anyone skilled in the art,
It should be covered by the protection scope of the present invention.Therefore, protection scope of the present invention should be with the protection model of claims
Subject to enclosing.
Claims (3)
1. a kind of generator dynamic state estimator method characterized by comprising
It is calculated according to generator dynamic state estimator result of the system equation to the k-1 moment, it is motor-driven to obtain the power generation of kth moment
State state estimation result；
Kth moment generator dynamic state estimator result is forecast and is filtered using CKF filtering algorithm；Wherein, lead to
Cross the generator dynamical state forecast result that forecast processing obtains the k+1 moment；In filtering processing, using obtaining the real-time of PMU
Measuring value zk+1=[δzk+1 ωzk+1 Pek+1]TCalculate new breath vector ek+1, and then obtain time-varying multidimensional observation noise scale factor；
Subscript k+1 indicates k+1 moment, δz、ωzAnd PeThe respectively absolute generator rotor angle of rotor, angular rate and the PMU amount for exporting electromagnetic power
Measured value；It introduces time-varying multidimensional observation noise scale factor and combines the real-time amount measured value of PMU to the generator dynamic shape at k+1 moment
State forecast result is filtered, to realize the accurate estimation of generator dynamical state result；
Wherein, the system equation are as follows:
In formula, subscript k+1 and k corresponding expression k+1 moment and k moment；F and H is respectively state equation function and measurement side
Eikonal number, x, u and z are respectively quantity of state, control amount and measurement；V and w is respectively system noise and measurement noise, error side
Poor battle array is respectively the normal distribution of Q and R；
Wherein, quantity of state x=[δ ω E 'qE′d]T, control amountThe concrete form of state equation are as follows:
In formula, δ is the absolute generator rotor angle per unit value of generator amature, and ω is angular rate per unit value；WithIt is the letter to time diffusion
Change literary style, the expression d/dt differential operator of letter top；TJFor generator inertia time constant；TmFor machine torque；UtWith
The respectively amplitude and phase angle of generator outlet voltage phasor；X′qWith X 'dRespectively q axis and d axis transient reactance；E′qWith E 'dPoint
It Wei not q axis and d axis transient state electromotive force；D is damped coefficient；T′q0With T 'd0Respectively q axis and d axis open circuit transient time-constant；Ef
For stator excitation electromotive force；XqAnd XdRespectively q axis and d axis synchronous reactance；
Measurement z=[δz ωz Pe]T, the concrete form of measurement equation are as follows:
In formula, δz、ωzAnd PeThe respectively absolute generator rotor angle of rotor, angular rate and the PMU measuring value for exporting electromagnetic power；
According to the concrete form and k moment system noise variance matrix Q of above-mentioned state equation and measurement equationkWith the k+1 moment
Measuring noise square difference battle array Rk+1, then can be realized generator dynamic state estimator；
K+1 moment measuring noise square difference battle array R thereink+1Value, the moment system noise side k are carried out according to the practical error in measurement of PMU
Poor matrix QkIt indicates are as follows:
In formula,WithRespectively k moment δ, ω, E 'qWith E 'dSystem noise variance；
It is calculated by formula of error transmission:
In formula: σ is noise variance；σUFor the amplitude PMU error in measurement standard deviation of generator outlet voltage phasor,For generator
The phase angle PMU error in measurement standard deviation of exit potential phasor；
Then have:
In formula, Δ t is step-length.
2. the method according to claim 1, wherein described handle the generator for obtaining the k+1 moment by forecast
The step of dynamical state forecast result includes:
The kth moment generator dynamic state estimator result includes: the estimated value of kth moment generator dynamic state quantityWith estimation error variance battle array Pk|k；
To Pk|kCholesky decomposition is carried out, the On Square-Rooting Matrices S of k moment estimation error variance battle array is obtainedk|k:
Using following formula to the estimated value of kth moment generator dynamic state quantityGenerate the volume point X of the weights such as one groupi,k|k:
In formula, parameterN is quantity of state dimension；
Each volume point is converted using following formula, obtains the predicted value of all volume points
Summation is weighted to the predicted value of all Generator Status amount volume points, obtains quantity of state predicted value
And prediction error conariance battle array P is obtained by following formulak+1|k:
3. introducing time-varying multidimensional observation is made an uproar the method according to claim 1, wherein described in filtering processing
Sound scale factor simultaneously combines the real-time amount measured value of PMU to be filtered the generator dynamical state forecast result at k+1 moment, thus
Realize generator dynamical state result it is accurate estimate include:
The generator dynamical state forecast result at the k+1 moment includes: the generator dynamic state quantity predicted value at k+1 momentWith Generator Status amount prediction error conariance battle array Pk+1|k；
To prediction error conariance battle array Pk+1|kCholesky decomposition is carried out, k+1 moment Generator Status amount prediction error association is obtained
The On Square-Rooting Matrices S of variance matrixk+1|k:
Using following formula to the generator dynamic state quantity predicted value at k+1 momentGenerate the quantity of state predicted value of the weights such as one group
Volume point Xi,k+1|k:
In formula, parameterN is quantity of state dimension；
Each quantity of state predicted value volume point is converted using following formula, obtains the volume point of PMU measurement predicted value
Zi,k+1|k:
Zi,k+1|k=H (Xi,k+1|k,uk)；
Weighted sum is carried out to the volume point of all PMU measurement predicted values, and then obtains PMU measurement predicted value
Calculate PMU measurement prediction error variance matrix Pvv,k+1:
Using obtaining the real-time amount measured value z of PMUk+1=[δzk+1 ωzk+1 Pek+1]TCalculate new breath vector ek+1, and then obtain time-varying
Multidimensional observation noise scale factor γk+1:
In formula, M is that the window of the windowing estimation technique is long；
Following formula is recycled to calculate diagonal matrix γ 'k+1:
γ′k+1=diag (γ '1,γ′2,…,γ′m)；
In formula, diagonal element γ '1Value are as follows: γ 'i=max { 1, γk+1,ii, i=1,2 ..., n；γk+1,iiFor γk+1?
I diagonal element；
The Cross-covariance P between Generator Status amount predicted value and PMU measurement predicted value is calculated according to the following formulaxz,k+1|k:
Kalman filtering gain W is calculated againk+1:
Wk+1=Pxz,k+1|k(Pvv,k+1+γ′k+1Rk+1)-1；
Vector e is newly ceased using the k+1 momentk+1, and pass through Kalman filtering gain Wk+1To the generator dynamic state quantity at k+1 moment
Predicted valueIt is filtered, obtains the estimated value of k+1 moment generator dynamic state quantity
And it is calculate by the following formula and obtains generator dynamic state quantity estimation error variance battle array Pk+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 CN105403834A (en) | 2016-03-16 |
| CN105403834B true 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) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR102104117B1 (en) * | 2016-07-25 | 2020-04-23 | 미쓰비시덴키 가부시키가이샤 | Diagnostic device of electric 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 |
| CN107425548B (en) * | 2017-09-11 | 2020-06-16 | 河海大学 | Interpolation H∞Dynamic state estimation method for extended Kalman filter generator |
| CN107478990B (en) * | 2017-09-11 | 2019-11-12 | 河海大学 | A kind of generator electromechanical transient process method for dynamic estimation |
| CN109100649B (en) * | 2018-06-25 | 2020-10-16 | 南京南瑞继保电气有限公司 | Parameter estimation method for generator excitation system and speed regulation system based on phasor measurement |
| 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 |
| CN111948534B (en) * | 2020-07-31 | 2023-05-05 | 华北电力科学研究院有限责任公司 | 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 |
|---|
| 基于容积卡尔曼滤波的发电机动态状态估计;陈亮等;《中国电机工程学报》;20140605;第34卷(第16期);第2708页第2节 |
| 基于扩展集员滤波的发电机动态状态估计;郑维荣等;《电力科学与工程》;20141130;第30卷(第11期);第11-15页 |
| 电力系统分布式动态状态估计研究;陈亮;《中国博士学位论文全文数据库工程科技II辑》;20141215(第12期);第11页第2.2节、第17页第2.3.3节-2.3.4节、第34页第3.2.2、第37页第3.3节 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105403834A (en) | 2016-03-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105403834B (en) | A kind of generator dynamic state estimator method | |
| Mai et al. | A dynamic synchrophasor estimation algorithm for online application | |
| US9178386B2 (en) | State-matrix-independent dynamic process estimation method in real-time for weakly observable measurement nodes without PMU | |
| Zografos et al. | Estimation of power system inertia | |
| US8462004B2 (en) | Method and arrangement for generating an error signal | |
| CN106844952A (en) | Based on the generator dynamic state estimator method without mark Particle filtering theory | |
| CN102510263B (en) | Method for identifying practical parameters of synchronous generator on basis of load rejection test and numerical difference | |
| US10884060B1 (en) | Dynamic parameter estimation of generators | |
| CN107590317A (en) | A kind of generator method for dynamic estimation of meter and model parameter uncertainty | |
| Mei et al. | Clustering-based dynamic event location using wide-area phasor measurements | |
| CN103401238B (en) | A kind of power load modelling approach based on Measurement-based approach | |
| CN108574291A (en) | One kind being based on Ensemble Kalman Filter generator dynamic state estimator method | |
| CN107851995A (en) | For detecting the method and apparatus for supplying the voltage in power network | |
| CN107658881A (en) | Voltage stability critical point determination methods based on Thevenin's equivalence method | |
| US20220399721A1 (en) | Method for grid impedance and dynamics estimation | |
| CN103983847A (en) | Self-adaptive frequency tracking measurement method based on RLS (Recursive Least Squares) in synchronized phasor measurement | |
| Susuki et al. | Estimation of power system inertia using nonlinear koopman modes | |
| Zacharia et al. | Measurement errors and delays on wide-area control based on IEEE Std C37. 118.1-2011: impact and compensation | |
| Bila | Power system dynamic state estimation and load modeling | |
| CN112511056A (en) | Robust generator dynamic state estimation method based on phasor measurement | |
| CN111585285A (en) | Load modeling method and device based on voltage monitoring system | |
| Cokkinides et al. | Visualization and characterization of stability swings via GPS-synchronized data | |
| Xue et al. | Power system frequency measurement for frequency relaying | |
| Lee et al. | A two-stage data-driven algorithm to estimate the system inertia utilizing event-driven disturbed PMU measurements | |
| Al-Saadi et al. | Voltage stability evaluation by using extended modal analysis with dynamic internal impedance of synchronous generators |
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 |