CN105403834A - Dynamic state evaluation method of generator - Google Patents

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

A kind of generator dynamic state estimator method

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:

$\left\{\begin{array}{c}{x}_{k+1}=F(x_k,u_k,v_k)\\ z_{k+1}=H(x_{k+1},u_{k+1},w_{k+1})\end{array}\right.$

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:

$Q_k=\left[\begin{array}{cccc}{\mathrm{\σ}}_{m\mathrm{\δ}k}^2& 0& 0& 0\\ 0& {\mathrm{\σ}}_{m\mathrm{\ω}k}^2& 0& 0\\ 0& 0& {\mathrm{\σ}}_{{\mathrm{mE}}_q^{\′}k}^2& 0\\ 0& 0& 0& {\mathrm{\σ}}_{{\mathrm{mE}}_d^{\′}k}^2\end{array}\right];$

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:

${\mathrm{\σ}}_{m\mathrm{\δ}k}^2={10}^{-6};$

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:

$S_{k|k}=\sqrt{P_{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

$X_{i,k+1|k}^{*}=F(X_{i,k|k},u_k);$

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:

$S_{k+1|k}=\sqrt{P_{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:

${\mathrm{\γ}}_{k+1}=(\frac{1}{M}\underset{j=0}{\overset{M-1}{\Σ}}{e}_{k+1-j}{e}_{k+1-j}^T-P_{vv,k+1})R_{k+1}^{-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(γ′ ₁,γ′ ₂,…,γ′ _m)；

In formula, diagonal element γ ' ₁value be: γ ' _i=max{1, γ _{k+1, ii}, i=1,2 ..., n; γ _{k+1, ii}for γ _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:

$P_{k+1|k+1}=P_{k+1|k}-W_{k+1}{P}_{zz,k+1|k}{W}_{k+1}^T.$

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:

$\left\{\begin{array}{c}{x}_{k+1}=F(x_k,u_k,v_k)\\ z_{k+1}=H(x_{k+1},u_{k+1},w_{k+1})\end{array}\right.---\left(1\right)$

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:

$Q_k=\left[\begin{array}{cccc}{\mathrm{\σ}}_{m\mathrm{\δ}k}^2& 0& 0& 0\\ 0& {\mathrm{\σ}}_{m\mathrm{\ω}k}^2& 0& 0\\ 0& 0& {\mathrm{\σ}}_{{\mathrm{mE}}_q^{\′}k}^2& 0\\ 0& 0& 0& {\mathrm{\σ}}_{{\mathrm{mE}}_d^{\′}k}^2\end{array}\right];---\left(4\right)$

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:

${\mathrm{\σ}}_{m\mathrm{\δ}k}^2={10}^{-6};---\left(6\right)$

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:

$S_{k|k}=\sqrt{P_{k|k}};---\left(10\right)$

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

$X_{i,k+1|k}^{*}=F(X_{i,k|k},u_k);---\left(12\right)$

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:

$S_{k+1|k}=\sqrt{P_{k+1|k}};---\left(15\right)$

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:

${\mathrm{\γ}}_{k+1}=(\frac{1}{M}\underset{j=0}{\overset{M-1}{\Σ}}{e}_{k+1-j}{e}_{k+1-j}^T-P_{vv,k+1})R_{k+1}^{-1};---\left(21\right)$

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(γ′ ₁,γ′ ₂,…,γ′ _m)；(22)

In formula, diagonal element γ ' ₁value be: γ ' _i=max{1, γ _{k+1, ii}, i=1,2 ..., n; γ _{k+1, ii}for γ _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:

$P_{k+1|k+1}=P_{k+1|k}-W_{k+1}{P}_{zz,k+1|k}{W}_{k+1}^T.---\left(26\right)$

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 $W_{k+1}=P_{xz,k+1|k}{P}_{zz,k+1|k}^{-1}$ 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

$e_{k+1}{e}_{k+1}^T=P_{vv,k+1}+{\mathrm{\γ}}_{k+1}{R}_{k+1}---\left(28\right)$

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:

$P_{e,k+1}=\frac{1}{M}\underset{j=0}{\overset{M-1}{\Σ}}{e}_{k+1-j}{e}_{k+1-j}^T---\left(29\right)$

In formula, P _{e, k+1}for 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

$\frac{1}{M}\underset{j=0}{\overset{M-1}{\Σ}}{e}_{k+1-j}{e}_{k+1-j}^T=P_{vv,k+1}+{\mathrm{\γ}}_{k+1}{R}_{k+1}---\left(30\right)$

Solved by following formula (31) and obtain γ _k+1:

${\mathrm{\γ}}_{k+1}=(\frac{1}{M}\underset{j=0}{\overset{M-1}{\Σ}}{e}_{k+1-j}{e}_{k+1-j}^T-P_{vv,k+1})R_{k+1}^{-1}---\left(31\right)$

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(γ′ ₁,γ′ ₂,…,γ′ _m)(32)

Diagonal element γ ' in formula (32) _ivalue be

γ′ _i＝max{1,γ _k+1,ii}，i＝1,2,…,n(33)

In formula, γ _{k+1, ii}for γ _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:

$\left\{\begin{array}{c}{x}_{k+1}=F(x_k,u_k,v_k)\\ z_{k+1}=H(x_{k+1},u_{k+1},w_{k+1})\end{array}\right.$

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:

$Q_k=\left[\begin{array}{cccc}{\mathrm{\σ}}_{m\mathrm{\δ}k}^2& 0& 0& 0\\ 0& {\mathrm{\σ}}_{m\mathrm{\ω}k}^2& 0& 0\\ 0& 0& {\mathrm{\σ}}_{{\mathrm{mE}}_q^{\′}k}^2& 0\\ 0& 0& 0& {\mathrm{\σ}}_{{\mathrm{mE}}_d^{\′}k}^2\end{array}\right];$

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:

${\mathrm{\σ}}_{m\mathrm{\δ}k}^2={10}^{-6};$

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:

$S_{k|k}=\sqrt{P_{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

$X_{i,k+1|k}^{*}=F(X_{i,k|k},u_k);$

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:

$S_{k+1|k}=\sqrt{P_{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:

$y_{k+1}=(\frac{1}{M}\underset{j=0}{\overset{M-1}{\Σ}}{e}_{k+1-j}{e}_{k+1-j}^T-P_{vv,k+1})R_{k+1}^{-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(γ′ ₁,γ′ ₂,…,γ′ _m)；

In formula, diagonal element γ ' ₁value be: γ ' _i=max{1, γ _{k+1, ii}, i=1,2 ..., n; γ _{k+1, ii}for γ _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:

$P_{k+1|k+1}=P_{k+1|k}-W_{k+1}{P}_{zz,k+1|k}{W}_{k+1}^T.$