CN102868157A - Robust estimation state estimating method based on maximum index absolute value target function - Google Patents

Abstract

The invention provides a robust estimation state estimating method based on a maximum index absolute value target function. The robust estimation state estimating method comprises the steps of: providing a robust state estimation basic model based on the maximum index absolute value target function; introducing an auxiliary variable into the robust state estimation basic model based on the maximum index absolute value target function to obtain a robust state estimation equivalent model based on the maximum index absolute value target function; and by using a primal-dual interior point algorithm, solving the robust state estimation equivalent model based on the maximum index absolute value target function. Proved through example analysis, the robust estimation state estimating method has the advantages of strong robustness and high computing efficiency, and better engineering application prospect.

Description

A kind of anti-poor method for estimating state based on maximal index absolute value target function

Technical field

The invention belongs to the dispatching automation of electric power systems field, be specifically related to a kind of anti-poor method for estimating state based on maximal index absolute value target function.

Background technology

Power system state estimation is basis and the core of EMS.Almost state estimator has been installed by each large-scale control centre now, and state estimation has become the foundation stone of electric power netting safe running.Since 1970 foreign scholars proposed state estimation first, people had had the history in more than 40 year to research and the application of state estimation, have emerged various method for estimating state during this.

At present, the state estimation that at home and abroad is most widely used is weighted least-squares method (Weighted least squares, WLS).The WLS model simple find the solution easily, but its anti-poor property is very poor.In order to strengthen anti-poor property, two kinds of methods are arranged generally.The first is to add bad data identification link after WLS estimates, such as maximum regularization residual test method (LNR) or estimation discrimination method etc.; Another kind is to adopt anti-poor method for estimating state (Robust state estimation).At present, the anti-poor method for estimating state that Chinese scholars has proposed comprises weighting least absolute value estimation (Weighted least absolute value, WLAV), Non quadratic criteria method (QL, QC etc.), be state estimation (the Maximum normal measurement rate of target to the maximum with qualification rate, MNMR) and exponential type target function state estimation (Maximum exponential square, MES) etc.But the not high enough characteristics of these anti-poor method for estimating state ubiquity computational efficiencies, thereby their application in real system have been affected to a certain extent.

Summary of the invention

The present invention one of is intended to solve the problems of the technologies described above at least to a certain extent or provides at least a kind of useful commerce to select.For this reason, one object of the present invention is to propose the anti-poor method for estimating state based on maximal index absolute value target function that a kind of anti-poor property is good, computational efficiency is high (Maximum exponential absolute value state estimation, MEAV).

The anti-poor method for estimating state based on maximal index absolute value target function according to the embodiment of the invention comprises step: steps A. the anti-poor state estimation basic model based on maximal index absolute value target function is provided; Step B. introduces auxiliary variable to described anti-poor state estimation basic model based on maximal index absolute value target function, and conversion obtains the anti-poor state estimation equivalence model based on maximal index absolute value target function; And step C. utilizes former-dual interior point, and described anti-poor state estimation equivalence model based on maximal index absolute value target function is found the solution.

In one embodiment of the invention, described anti-poor state estimation basic model based on maximal index absolute value target function is:

S.t.g (x)=0, r=z-h (x), wherein: z ∈ R ^mFor measuring vector, comprise that the node injection is meritorious and idle, branch road meritorious and idle and node voltage amplitude measurement; X ∈ R ⁿBe state vector, comprise other each node phase angles except node voltage amplitude and the balance node; H:R ⁿ→ R ^mFor by state vector to the Nonlinear Mapping that measures vector; r _iI element for residual error vector r; G (x): R ⁿ→ R ^cIt is zero injecting power equality constraint; w _iBe the weight of i measurement amount, be both and be the window width parameter.

In one embodiment of the invention, described step B comprises: introduce non-negative slack variable u, v ∈ R ^m, the described anti-poor state estimation equivalence model based on maximal index absolute value target function that conversion obtains is:

S.t.g (x)=0, z-h (x)-u+v=0, u, v 〉=0.

In one embodiment of the invention, described step C comprises: step C1: make that x is flat starting state variable; Select λ ⁽⁰⁾=π ⁽⁰⁾=0 and u ⁽⁰⁾, v ⁽⁰⁾, α ⁽⁰⁾, β ⁽⁰⁾0; Make Center Parameter ρ ∈ (0,1) and convergence criterion ε=10 ^-3, put iteration count k=0; Step C2: calculate duality gap Gap=α ^TV+ β ^TU judges whether convergence, if Gap＜ε then turns step C7, otherwise enters step C3; Step C3: find the solution update equation, to finish the correction to former variable and dual variable, obtain [dx ^TD λ ^TD π ^T] ^T, dv, du, d α and d β; Step C4: the correction step-length θ that calculates former problem and dual problem _PAnd θ _D, wherein: ${\mathrm{\&theta;}}_P=0.9995\min \{\underset{i}{\min }(-\frac{v_i}{d{v}_i}:d{v}_i<0;-\frac{u_i}{d{u}_i}:d{u}_i<0),1\},$ ${\mathrm{\&theta;}}_D=0.9995\min \{\underset{i}{\min }(-\frac{{\mathrm{\&alpha;}}_i}{d{\mathrm{\&alpha;}}_i}:d{\mathrm{\&alpha;}}_i<0;-\frac{{\mathrm{\&beta;}}_i}{d{\mathrm{\&beta;}}_i}:d{\mathrm{\&beta;}}_i<0),1\};$ Step C5: the variable of revising respectively former problem and dual problem is: $\left[\begin{array}{c}{x}^{(k+1)}\\ v^{(k+1)}\\ u^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{x}^{\left(k\right)}\\ v^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_P\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{dv}\\ \mathrm{du}\end{array}\right],$ $\left[\begin{array}{c}{\mathrm{\&lambda;}}^{(k+1)}\\ {\mathrm{\&pi;}}^{(k+1)}\\ {\mathrm{\&alpha;}}^{(k+1)}\\ {\mathrm{\&beta;}}^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&lambda;}}^{\left(k\right)}\\ {\mathrm{\&pi;}}^{\left(k\right)}\\ {\mathrm{\&alpha;}}^{\left(k\right)}\\ {\mathrm{\&beta;}}^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_D\left[\begin{array}{c}\mathrm{d\&lambda;}\\ \mathrm{d\&pi;}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\end{array}\right];$ Step C6: make iteration count k=k+1, enter step C2; And step C7: the output optimal solution, finish.

In one embodiment of the invention, described step C3 comprises: step C31: calculation perturbation parameter μ=ρ Gap/2m; Step C32: form measurement equation and Jacobian matrix corresponding to zero injecting power constraint

And

Form measurement equation and extra large gloomy matrix corresponding to zero injecting power constraint And

Wherein h (x) is the measurement estimated value for the mapping of state vector to the measurement vector, and g (x)=0 is zero injecting power constraint; Step C33: calculate L _x=G ^Tλ-H ^Tπ, L _λ=g (x), L _π=z-h (x)-u+v, $L_{v_i}=-{\mathrm{\&omega;}}_i-{\mathrm{\&pi;}}_i-{\mathrm{\&alpha;}}_i,$ $L_{u_i}=-{\mathrm{\&omega;}}_i-{\mathrm{\&pi;}}_i-{\mathrm{\&beta;}}_i,$ $L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}}={\mathrm{\&alpha;}}_i{v}_i-\mathrm{\&mu;}$ And $L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}={\mathrm{\&beta;}}_i{u}_i-\mathrm{\&mu;},$ ω wherein _i=exp ((u _i+ v _i)/w _i), w _iBe i the weight that the measurement amount is corresponding; Step C34: calculate γ=z-h (x)-u+v+AA-BB, AA wherein, BB ∈ R ^m, ${\mathrm{AA}}_i=-a_i(v_i{L}_{v_i}+L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}})-b_i(u_i{L}_{u_i}+L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}),$ ${\mathrm{BB}}_i=-c_i(v_i{L}_{v_i}+L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}})-d_i(u_i{L}_{u_i}+L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}),$ $\left[\begin{array}{cc}{a}_i& b_i\\ c_i& d_i\end{array}\right]={\left[\begin{array}{cc}{v}_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&alpha;}}_i& v_i{\mathrm{\&omega;}}_i{w}_i^{-1}\\ u_i{\mathrm{\&omega;}}_i{w}_i^{-1}& u_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&beta;}}_i\end{array}\right]}^{-1},$ Z ∈ R ^mFor measuring vector; Step C35: solving equation $\left[\begin{array}{ccc}{\&dtri;}^{2}g\left(x\right)\mathrm{\&lambda;}-{\&dtri;}^{2}h\left(x\right)\mathrm{\&pi;}& G^T& -H^T\\ H& 0& \mathrm{<figure-callout\; id="0"\; label="Q\; G"\; filenames="HDA00002126017800012.png,HDA00002126017800013.png"\; state="\{\{state\}\}">Q</figure-callout>}& \\ G& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{d\&lambda;}\\ \mathrm{d\&pi;}\end{array}\right]=\left[\begin{array}{c}-L_x\\ \mathrm{\&gamma;}\\ -L_{\mathrm{\&lambda;}}\end{array}\right]$ Obtain [dx ^TD λ ^TD π ^T] ^TStep C36: find the solution dv _i=k _1iD π _i+ AA _iAnd du _i=k _2iD π _i+ BB _i, k wherein _1i=a _iv _i-b _iu _i, k _2i=c _iv _i-d _iu _iAnd step C37: find the solution $d{\mathrm{\&alpha;}}_i={\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i-d{\mathrm{\&pi;}}_i+L_{v_i}$ And $d{\mathrm{\&beta;}}_i={\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i+d{\mathrm{\&pi;}}_i+L_{u_i}.$

The anti-poor method for estimating state based on maximal index absolute value target function of the embodiment of the invention (Maximum exponential absolute value state estimation, MEAV) but establishment comprises a plurality of bad datas of consistency bad data in estimation procedure, shown good anti-poor property, and have very high computational efficiency, be suitable for very much practical engineering application.

Additional aspect of the present invention and advantage in the following description part provide, and part will become obviously from the following description, or recognize by practice of the present invention.

Description of drawings

Above-mentioned and/or additional aspect of the present invention and advantage are from obviously and easily understanding becoming the description of embodiment in conjunction with following accompanying drawing, wherein:

Fig. 1 is the schematic flow sheet of the anti-poor method for estimating state based on maximal index absolute value target function of the present invention;

Fig. 2 is IEEE 300 system node voltage magnitude maximum estimated error curves;

Fig. 3 is IEEE 300 system node voltage phase angle maximum estimated error curves: and

Fig. 4 is that the qualification rate of four kinds of state estimators of IEEE 300 systems compares.

Embodiment

The below describes embodiments of the invention in detail, and the example of described embodiment is shown in the drawings, and wherein same or similar label represents same or similar element or the element with identical or similar functions from start to finish.Be exemplary below by the embodiment that is described with reference to the drawings, be intended to for explaining the present invention, and can not be interpreted as limitation of the present invention.

As shown in Figure 1, the anti-poor method for estimating state based on maximal index absolute value target function of the embodiment of the invention comprises the following steps:

Steps A: the basic model of the anti-poor state estimation (Maximum exponential absolute value state estimation, MEAV) based on maximal index absolute value target function is provided.

Particularly, the basic model of the MEAV of the present invention's proposition is as follows

$\underset{x}{\max }J\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{w}_i\exp (-\frac{\left|r_i\right|}{w_i})---\left(1\right)$

s.t.g(x)=0 (2)

r=z-h(x) (3)

In the formula: z ∈ R ^mFor measuring vector, often comprise that the node injection is meritorious and idle, branch road meritorious and idle and the measurement of node voltage amplitude etc.; X ∈ R ⁿIt is the state vector (except the balance node phase angle) that comprises node voltage amplitude and phase angle; H:R ⁿ→ R ^mFor by state vector to the Nonlinear Mapping that measures vector; r _iI the element of residual error vector r; G (x): R ⁿ→ R ^cIt is zero injecting power equality constraint; w _iBe the weight of i measurement amount, be both and be the window width parameter.

Step B: the anti-poor state estimation basic model based on maximal index absolute value target function is introduced auxiliary variable, and conversion obtains the anti-poor state estimation equivalence model based on maximal index absolute value target function.

Particularly, although the target function everywhere continuous of MEAV basic model can not lead at 0 place, thereby direct solution relatively the difficulty.Model (1) ~ (3) can be converted into the equivalence model that can lead everywhere.

Introduce variable ξ ∈ R ^m, it is satisfied

|r _i|≤ξ _i i＝1,…,m (4)

Can be got by formula (1) and (4),

Maximum is equivalent to

Maximum.

Introduce non-negative slack variable l, k ∈ R ^m, inequality (4) be converted into two equality constraints be

r _i-l _i=-ξ _i i=1,…,m (5)

r _i+k _i＝ξ _i i=1，…,m (6)

Introduce non-negative slack variable u, v ∈ R ^m, it is satisfied

u _i＝l _i/2 i=1,…,m (7)

v _i＝k _i/2 i=1,…,m (8)

By formula (5) ~ (8), can get

r _i＝u _i-v _i i=1,…,m (9)

ξ _i=u _i+v _i i=1,…,m (10)

Bring formula (9) into formula (3), can get measurement constraints of equal value and be

z-h(x)-u+v=0 (11)

Then the equivalence model of the MEAV basic model that provides of formula (1) ~ (3) is

$\underset{x}{\max }J\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{w}_i\exp (-\frac{u_i+v_i}{w_i})---\left(12\right)$

s.t.g(x)=0 (13)

z-h(x)-u+v=0 (14)

u,v ≥0 (15)

Model (12) ~ (15) are the MEAV equivalence model, and this model everywhere continuous can be led, and can find the solution with the method based on gradient.Step C: utilize former-dual interior point, described anti-poor state estimation equivalence model based on maximal index absolute value target function is found the solution.

(1) method for solving of MEAV equivalence model

Equivalence model (12) ~ (15) of noticing MEAV are optimization problems that contains equality constraint and inequality constraints, suit to find the solution with former-dual interior point.For making those skilled in the art understand better the present invention, the derivation that given first is detailed is as follows:

Introduce Lagrangian

$L\&equiv;\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\left[w_i\exp (-\frac{u_i+v_i}{w_i})\right]-{\mathrm{\&lambda;}}^{T}g\left(x\right)-{\mathrm{\&pi;}}^T(z-h\left(x\right)-u+v)-{\mathrm{\&alpha;}}^{T}v-{\mathrm{\&beta;}}^{T}u---\left(16\right)$

In the formula: λ ∈ R ^cAnd π, α, β ∈ R ^mBe the Lagrange multiplier vector.

For obtaining optimal value, according to the KKT condition, can get

$L_x\&equiv;\&PartialD;L/\&PartialD;x=G^T\mathrm{\&lambda;}-H^T\mathrm{\&pi;}=0---\left(17\right)$

$L_{\mathrm{\&lambda;}}\&equiv;\&PartialD;L/\&PartialD;\mathrm{\&lambda;}=g\left(x\right)=0---\left(18\right)$

$L_{\mathrm{\&pi;}}\&equiv;\&PartialD;L/\&PartialD;\mathrm{\&pi;}=z-h\left(x\right)-u+v=0---\left(19\right)$

$L_{v_i}\&equiv;\&PartialD;L/\&PartialD;v_i=-{\mathrm{\&omega;}}_i-{\mathrm{\&pi;}}_i-{\mathrm{\&alpha;}}_i=0---\left(20\right)$

$L_{u_i}\&equiv;\&PartialD;L/\&PartialD;u_i=-{\mathrm{\&omega;}}_i+{\mathrm{\&pi;}}_i-{\mathrm{\&beta;}}_i=0---\left(21\right)$

$L_{{\mathrm{\&alpha;}}_i}\&equiv;\&PartialD;L/\&PartialD;{\mathrm{\&alpha;}}_i={\mathrm{\&alpha;}}_i{v}_i=0---\left(22\right)$

$L_{{\mathrm{\&beta;}}_i}\&equiv;\&PartialD;L/\&PartialD;{\mathrm{\&beta;}}_i={\mathrm{\&beta;}}_i{u}_i=0---\left(23\right)$

In the formula: ω _i=exp ((u _i+ v _i)/w _i), $H=\&PartialD;h\left(x\right)/\&PartialD;x,$ $G=\&PartialD;g\left(x\right)/\&PartialD;x.$

For effectively overcoming the above problems, the modern interior-point method is introduced disturbance parameter μ〉0 pair of formula (22), (23) relax, thereby

$L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}}\&equiv;{\mathrm{\&alpha;}}_i{v}_i-\mathrm{\&mu;}=0---\left(24\right)$

$L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}\&equiv;{\mathrm{\&beta;}}_i{u}_i-\mathrm{\&mu;}=0---\left(25\right)$

Above equation can be got by Newton Algorithm

$[{\&dtri;}^{2}g\left(x\right)\mathrm{\&lambda;}-{\&dtri;}^{2}h\left(x\right)\mathrm{\&pi;}]\mathrm{dx}+G^T\mathrm{d\&lambda;}-H^T\mathrm{d\&pi;}=-L_x---\left(26\right)$

Gdx=-L _λ (27)

-Hdx-du+dv=-L _π (28)

${\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i-d{\mathrm{\&pi;}}_i-d{\mathrm{\&alpha;}}_i=-L_{v_i}---\left(29\right)$

${\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i+d{\mathrm{\&pi;}}_i-d{\mathrm{\&beta;}}_i=-L_{u_i}---\left(30\right)$

${\mathrm{\&alpha;}}_{i}d{v}_i+v_{i}d{\mathrm{\&alpha;}}_i=-L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}}---\left(31\right)$

${\mathrm{\&beta;}}_{i}d{u}_i+u_{i}d{\mathrm{\&beta;}}_i=-L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}---\left(32\right)$

By formula (29) and (30), can get

$d{\mathrm{\&alpha;}}_i={\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i-d{\mathrm{\&pi;}}_i+L_{v_i}---\left(33\right)$

$d{\mathrm{\&beta;}}_i={\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i+d{\mathrm{\&pi;}}_i+L_{u_i}---\left(34\right)$

With formula (33), (34) bring (31) into, (32) can get

$[v_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&alpha;}}_i]d{v}_i+v_i{\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i=v_{i}d{\mathrm{\&pi;}}_i-v_i{L}_{v_i}-L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}}---\left(35\right)$

$u_i{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i+[u_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&beta;}}_i]d{u}_i=-u_{i}d{\mathrm{\&pi;}}_i-u_i{L}_{u_i}-L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}---\left(36\right)$

Order $\left[\begin{array}{cc}{a}_i& b_i\\ c_i& d_i\end{array}\right]={\left[\begin{array}{cc}{v}_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&alpha;}}_i& v_i{\mathrm{\&omega;}}_i{w}_i^{-1}\\ u_i{\mathrm{\&omega;}}_i{w}_i^{-1}& u_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&beta;}}_i\end{array}\right]}^{-1},$ Can be got by formula (35) and (36)

dv _i＝k _1idπ _i+AA _i (37)

du _i=k _2idπ _i+BB _i (38)

In the formula: k _1i=a _iv _i-b _iu _i, k _2i=c _iv _i-d _iu _i, ${\mathrm{AA}}_i=-a_i(v_i{L}_{v_i}+L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}})-b_i(u_i{L}_{u_i}+L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}),{\mathrm{BB}}_i=-c_i(v_i{L}_{v_i}+L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}})-d_i(u_i{L}_{u_i}+L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}).$

Formula (37), (38) are brought into (28) and can be got

Hdx+Qdπ＝γ (39)

In the formula: Q is R ^{M * m}Diagonal matrix, its diagonal element is Q _Ii=-k _1i+ k _2iγ=z-h (x)-u+v+AA-BB, AA, BB ∈ R ^m, AA _i, BB _iWith same in formula (37), (38).

According to formula (39), (26) and (27), can get update equation and be

$\left[\begin{array}{ccc}{\&dtri;}^{2}g\left(x\right)\mathrm{\&lambda;}-{\&dtri;}^{2}h\left(x\right)\mathrm{\&pi;}& G^T& -H^T\\ H& 0& \mathrm{<figure-callout\; id="0"\; label="Q\; G"\; filenames="HDA00002126017800012.png,HDA00002126017800013.png"\; state="\{\{state\}\}">Q</figure-callout>}& \\ G& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{d\&lambda;}\\ \mathrm{d\&pi;}\end{array}\right]=\left[\begin{array}{c}-L_x\\ \mathrm{\&gamma;}\\ -L_{\mathrm{\&lambda;}}\end{array}\right]---\left(40\right)$

Find the solution formula (40) and can get dx, d λ and d π; Can get dv and du by formula (37), (38); With acquired results bring formula (33) into, (34) can get d α and d β, then iteration is sustainable carrying out.

(2) solution procedure of MEAV equivalence model

Introducing the finding the solution after the derivation of MEAV equivalence model, the inventor is summarized as follows solution procedure:

Step C1: carry out initialization, make that x is flat starting state variable; Select λ ⁽⁰⁾=π ⁽⁰⁾=0 and u ⁽⁰⁾, v ⁽⁰⁾, α ⁽⁰⁾, β ⁽⁰⁾0; Make Center Parameter ρ ∈ (0,1) and determine the convergence criterion value, and to put iteration count be zero.

Particularly, make x ⁽⁰⁾∈ R ⁿRepresentative by all node voltage amplitudes and phase angle form flat starting state variable (except the reference node phase angle); Select λ ⁽⁰⁾=π ⁽⁰⁾=0 and u ⁽⁰⁾, v ⁽⁰⁾, α ⁽⁰⁾, β ⁽⁰⁾0, λ ∈ R wherein ^cAnd π, α, β ∈ R ^mBe the Lagrange multiplier vector, m is the number of measurement amount, and c is the number of zero injecting power constraint; Make Center Parameter ρ ∈ (0,1) and convergence criterion ε=10 ^-3, put iteration count k=0.

Step C2: calculate duality gap Gap=α ^TV+ β ^TU judges whether convergence.Particularly, if Gap＜ε then thinks convergence can directly enter step C7; Otherwise for not restraining, enter step C3.

Step C3: find the solution update equation, to finish the correction to former variable and dual variable, obtain [dx ^TD λ ^TD π ^T] ^T, dv, du, d α and d β.

Particularly, at first then calculation perturbation parameter μ=ρ Gap/2m finds the solution formula (40) and gets [dx ^TD λ ^TD π ^T] ^T; Find the solution formula (37), (38) dv, du; Find the solution (33), (34) d α, d β.

Step C4: the correction step-length θ that calculates former problem and dual problem _PAnd θ _D, wherein: ${\mathrm{\&theta;}}_P=0.9995\min \{\underset{i}{\min }(-\frac{v_i}{d{v}_i}:d{v}_i<0;-\frac{u_i}{d{u}_i}:d{u}_i<0),1\},$ ${\mathrm{\&theta;}}_D=0.9995\min \{\underset{i}{\min }(-\frac{{\mathrm{\&alpha;}}_i}{d{\mathrm{\&alpha;}}_i}:d{\mathrm{\&alpha;}}_i<0;-\frac{{\mathrm{\&beta;}}_i}{d{\mathrm{\&beta;}}_i}:d{\mathrm{\&beta;}}_i<0),1\}$

Step C5: the variable of revising respectively former problem and dual problem is:

$\left[\begin{array}{c}{x}^{(k+1)}\\ v^{(k+1)}\\ u^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{x}^{\left(k\right)}\\ v^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_P\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{dv}\\ \mathrm{du}\end{array}\right],\left[\begin{array}{c}{\mathrm{\&lambda;}}^{(k+1)}\\ {\mathrm{\&pi;}}^{(k+1)}\\ {\mathrm{\&alpha;}}^{(k+1)}\\ {\mathrm{\&beta;}}^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&lambda;}}^{\left(k\right)}\\ {\mathrm{\&pi;}}^{\left(k\right)}\\ {\mathrm{\&alpha;}}^{\left(k\right)}\\ {\mathrm{\&beta;}}^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_D\left[\begin{array}{c}\mathrm{d\&lambda;}\\ \mathrm{d\&pi;}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\end{array}\right];$

Step C6: make iteration count k=k+1, enter step C2; And

Step C7: the output optimal solution, finish.

Be the advantage that makes those skilled in the art understand better the present invention and understand the relative prior art of the present invention, the applicant further explains in conjunction with specific embodiments.

Setting utilizes ieee standard system test based on the performance of the MEAV of former-dual interior point.Test adopts full dose to survey, and measuring value obtains by Additive White Noise (average is 0, and standard deviation is τ) on the result who calculates in trend.Measure for voltage, get τ _V=0.005p.u; For power measurement, get τ _PQ=1MW/MVar.Test environment is PC, and CPU is that Intel (R) Core (TM) i3M370, dominant frequency are 2.40GHz, internal memory 2.00GB.

1. the comparison of robustness

The inventor compares MEAV of the present invention and other state estimators, tests the anti-poor property of MEAV.

(1) IEEE-14 system

In the IEEE-14 system 4 consistency bad data (P are set _1-2, Q _1-2, P ₁, Q ₁).Set bad measuring value and the right value of measurement amount are as shown in table 1.

Table 1MEAV is to the identification of IEEE14 system conformance bad data

Table1 Estimation of conforming bad data for the IEEE 14-bus system by MEAV

As a comparison, at first estimate with the WLS that widely uses, and carry out the identification (brief note is WLS+LNR) of bad data with LNR.The result of first identification is: the standardized residual of 10 measurement amounts is greater than threshold value (3.0), and these 10 measurement amounts are considered to suspicious data; Wherein the amount of standardized residual maximum is measured as P _2-1, rerun WLS after leaving out this measurement; Find P this moment ₂Standardized residual maximum.Above process circulation 4 times, 4 good measurement amounts are thought by mistake be suspicious data and being left out by LNR, but real bad data still exists.As seen, WLS+LNR can not identification consistency bad data.

The estimated result of using the MEAV method is as shown in table 1.Can find, even there is the consistency bad data in the measurement amount, the estimated value of MEAV and true value also can be coincide well.Test of many times at the IEEE other system shows that also MEAV can suppress bad data automatically in the process of estimation, have good anti-poor property.

(2) IEEE-300 system

In order to test the estimated performance of MEAV in larger system, this section is based on IEEE 300 node systems, respectively to the bad data ratio be 0% ~ 10% totally 11 kinds of situations test, include l-G simulation test under each ratio 50 times.Bad data is the error of iteration 50% on correct measuring value and producing.Fig. 2,3 has provided under the bad data of different proportion, the peaked mean value of the node voltage amplitude that WLS+LNR and MEAV obtain and phase angle evaluated error absolute value (50 tests), wherein, | Δ U _Max| and | Δ θ _Max| represent respectively the maximum of node voltage amplitude and phase angle evaluated error absolute value, η represents the ratio of bad data.

By Fig. 2,3 as seen, when not having bad data (η=0), the evaluated error of MEAV and WLS+LNR is all less.Yet, along with the increase of bad data ratio, the evaluated error of WLS+LNR increases rapidly, and the amplitude error that MEAV estimates is then substantially constant, and phase angle error only slightly increases (even also like this when the bad data ratio is increased to 10%), has embodied good robustness.

Further, we test the qualification rate of four kinds of state estimators (WLS+LNR, WLAV, MNMR and MEAV) when the different bad data ratio.50 times the test mean value as shown in Figure 4, wherein η and Ω represent respectively the ratio of bad data and the qualification rate of state estimation.As seen from Figure 4, along with the rising of bad data ratio, the qualification rate of WLS+LNR descends rapidly; And the decline of the qualification rate of the anti-poor estimator of its excess-three kind is relatively less.When same bad data ratio, the qualification rate of WLS+LNR is minimum, and the qualification rate of MEAV is the highest, and the qualification rate of MNMR is taken second place.This shows, in above four kinds of state estimators, the anti-poor property of MEAV is best.

2. the comparison of computational efficiency

The inventor tests four kinds of state estimator WLS, WLAV, MNMR and MEAV respectively under normal measurement condition in order to carry out efficient relatively, and wherein rear three kinds belong to anti-poor state estimator.In test, WLS adopts Newton Algorithm, and other three kinds of state estimation adopt interior point method to find the solution; And MNMR adopts two-phase method, and namely the phase I is carried out the WLS estimation, and second stage is calculated the estimated value of WLS as the initial value that MNMR estimates.

Carry out altogether l-G simulation test 50 times, iterations and average computation during the state estimation convergence are consuming time as shown in table 2.By as seen from Table 2, in these four kinds of state estimators, the computational efficiency of WLS is the highest; And in rear three kinds of anti-poor state estimators, the computational efficiency of MEAV is the highest; And along with the increase of system scale, the iterations of MEAV and calculate the very slow of growth consuming time, thereby MEAV is applicable to the estimation of actual large scale system.

The iterations of four kinds of state estimators of table 2 and calculate consuming time

Table2 Iterations and CPU time of the four estimators

In sum, but MEAV establishment in estimation procedure that the present invention proposes comprises a plurality of bad datas of consistency bad data, has shown good anti-poor property, and has had very high computational efficiency, is suitable for very much practical engineering application.

In the description of this specification, the description of reference term " embodiment ", " some embodiment ", " example ", " concrete example " or " some examples " etc. means to be contained at least one embodiment of the present invention or the example in conjunction with specific features, structure, material or the characteristics of this embodiment or example description.In this manual, the schematic statement of above-mentioned term not necessarily referred to identical embodiment or example.And the specific features of description, structure, material or characteristics can be with suitable mode combinations in any one or more embodiment or example.

Although the above has illustrated and has described embodiments of the invention, be understandable that, above-described embodiment is exemplary, can not be interpreted as limitation of the present invention, those of ordinary skill in the art is not in the situation that break away from principle of the present invention and aim can change above-described embodiment within the scope of the invention, modification, replacement and modification.

Claims (5)

1. the anti-poor method for estimating state based on maximal index absolute value target function is characterized in that, comprises step:

Steps A: the anti-poor state estimation basic model based on maximal index absolute value target function is provided;

Step B: described anti-poor state estimation basic model based on maximal index absolute value target function is introduced auxiliary variable, and conversion obtains the anti-poor state estimation equivalence model based on maximal index absolute value target function; And

Step C: utilize former-dual interior point, described anti-poor state estimation equivalence model based on maximal index absolute value target function is found the solution.

2. the anti-poor method for estimating state based on maximal index absolute value target function as claimed in claim 1 is characterized in that, described anti-poor state estimation basic model based on maximal index absolute value target function is:

S.t.g (x)=0, r=z-h (x), wherein: z ∈ R ^mFor measuring vector, comprise that the node injection is meritorious and idle, branch road meritorious and idle and node voltage amplitude measurement; X ∈ R ⁿBe state vector, comprise other each node phase angles except node voltage amplitude and the balance node; H:R ⁿ→ R ^mFor by state vector to the Nonlinear Mapping that measures vector; r _iI element for residual error vector r; G (x): R ⁿ→ R ^cIt is zero injecting power equality constraint; w _iBe the weight of i measurement amount, be both and be the window width parameter.

3. the anti-poor method for estimating state based on maximal index absolute value target function as claimed in claim 2 is characterized in that, described step B comprises: introduce non-negative slack variable u, v ∈ R ^m, the described anti-poor state estimation equivalence model based on maximal index absolute value target function that conversion obtains is:

S.t.g (x)=0; Z-h (x)-u+v=0; U, v 〉=0.

4. the anti-poor method for estimating state based on maximal index absolute value target function as claimed in claim 3 is characterized in that, described step C comprises:

Step C1: make that x is flat starting state variable; Select λ ⁽⁰⁾=π ⁽⁰⁾=0 and u ⁽⁰⁾, v ⁽⁰⁾, α ⁽⁰⁾, β ⁽⁰⁾0; Make Center Parameter ρ ∈ (0,1) and convergence criterion ε=10 ^-3, put iteration count k=0;

Step C2: calculate duality gap Gap=α ^TV+ β ^TU judges whether convergence, if Gap＜ε then turns step C7, otherwise enters step C3;

Step C3: find the solution update equation, to finish the correction to former variable and dual variable, obtain [dx ^TD λ ^TD π ^T] ^T, dv, du, d α and d β;

Step C4: the correction step-length θ that calculates former problem and dual problem _PAnd θ _D, wherein: ${\mathrm{\&theta;}}_P=0.9995\min \{\underset{i}{\min }(-\frac{v_i}{d{v}_i}:d{v}_i<0;-\frac{u_i}{d{u}_i}:d{u}_i<0),1\},$ ${\mathrm{\&theta;}}_D=0.9995\min \{\underset{i}{\min }(-\frac{{\mathrm{\&alpha;}}_i}{d{\mathrm{\&alpha;}}_i}:d{\mathrm{\&alpha;}}_i<0;-\frac{{\mathrm{\&beta;}}_i}{d{\mathrm{\&beta;}}_i}:d{\mathrm{\&beta;}}_i<0),1\};$

Step C5: the variable of revising respectively former problem and dual problem is:

$\left[\begin{array}{c}{x}^{(k+1)}\\ v^{(k+1)}\\ u^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{x}^{\left(k\right)}\\ v^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_P\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{dv}\\ \mathrm{du}\end{array}\right],$ $\left[\begin{array}{c}{\mathrm{\&lambda;}}^{(k+1)}\\ {\mathrm{\&pi;}}^{(k+1)}\\ {\mathrm{\&alpha;}}^{(k+1)}\\ {\mathrm{\&beta;}}^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&lambda;}}^{\left(k\right)}\\ {\mathrm{\&pi;}}^{\left(k\right)}\\ {\mathrm{\&alpha;}}^{\left(k\right)}\\ {\mathrm{\&beta;}}^{\left(k\right)}\end{array}\right]+{\mathrm{\&theta;}}_D\left[\begin{array}{c}\mathrm{d\&lambda;}\\ \mathrm{d\&pi;}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\end{array}\right];$

Step C6: make iteration count k=k+1, enter step C2; And

Step C7: the output optimal solution, finish.

5. the anti-poor method for estimating state based on maximal index absolute value target function as claimed in claim 4 is characterized in that, described step C3 comprises:

Step C31: calculation perturbation parameter μ=ρ Gap/2m;

Step C32: form measurement equation and Jacobian matrix corresponding to zero injecting power constraint

And

Form measurement equation and extra large gloomy matrix corresponding to zero injecting power constraint

And

Wherein h (x) is the measurement estimated value for the mapping of state vector to the measurement vector, and g (x)=0 is zero injecting power constraint;

Step C33: calculate L _x=G ^Tλ-H ^Tπ, L _λ=g (x), L _π=z-h (x)-u+v,

And

ω wherein _i=exp ((u _i+ v _i)/w _i), w _iBe i the weight that the measurement amount is corresponding;

Step C34: calculate γ=z-h (x)-u+v+AA-BB, AA wherein, BB ∈ R ^m, ${\mathrm{AA}}_i=-a_i(v_i{L}_{v_i}+L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}})-b_i(u_i{L}_{u_i}+L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}),$ ${\mathrm{BB}}_i=-c_i(v_i{L}_{v_i}+L_{{\mathrm{\&alpha;}}_i}^{\mathrm{\&mu;}})-d_i(u_i{L}_{u_i}+L_{{\mathrm{\&beta;}}_i}^{\mathrm{\&mu;}}),$ $\left[\begin{array}{cc}{a}_i& b_i\\ c_i& d_i\end{array}\right]={\left[\begin{array}{cc}{v}_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&alpha;}}_i& v_i{\mathrm{\&omega;}}_i{w}_i^{-1}\\ u_i{\mathrm{\&omega;}}_i{w}_i^{-1}& u_i{\mathrm{\&omega;}}_i{w}_i^{-1}+{\mathrm{\&beta;}}_i\end{array}\right]}^{-1},$ Z ∈ R ^mFor measuring vector;

Step C35: solving equation $\left[\begin{array}{ccc}{\&dtri;}^{2}g\left(x\right)\mathrm{\&lambda;}-{\&dtri;}^{2}h\left(x\right)\mathrm{\&pi;}& G^T& -H^T\\ H& 0& Q\\ G& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{d\&lambda;}\\ \mathrm{d\&pi;}\end{array}\right]=\left[\begin{array}{c}-L_x\\ \mathrm{\&gamma;}\\ -L_{\mathrm{\&lambda;}}\end{array}\right]$ Obtain [dx ^TD λ ^TD π ^T] ^T

Step C36: find the solution dv _i=k _1iD π _i+ AA _iAnd du _i=k _2iD π _i+ BB _i, k wherein _Ii=a _iv _i-b _iu _i, k _2i=c _iv _i-d _iu _iAnd

Step C37: find the solution $d{\mathrm{\&alpha;}}_i={\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i-d{\mathrm{\&pi;}}_i+L_{v_i},$ And $d{\mathrm{\&beta;}}_i={\mathrm{\&omega;}}_i{w}_i^{-1}d{u}_i+{\mathrm{\&omega;}}_i{w}_i^{-1}d{v}_i+d{\mathrm{\&pi;}}_i+L_{u_i}.$

Images