CN102868157B - 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 robust method for estimating state based on maximal index absolute value target function

Technical field

The invention belongs to dispatching automation of electric power systems field, be specifically related to a kind of robust 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 propose state estimation first, people have had the history of more than 40 year to the research of state estimation and application, have emerged various method for estimating state during this.

At present, the state estimation being at home and abroad most widely used is weighted least-squares method (Weighted least squares, WLS).WLS model simple, solve easily, but its robust is very poor.In order to strengthen robust, generally there are two kinds of methods.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 robust method for estimating state (Robust state estimation).At present, the robust 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.), with qualification rate, be state estimation (the Maximum normalmeasurement rate of target to the maximum, MNMR) and exponential type target function state estimation (Maximum exponential square, MES) etc.But the not high enough feature of these robust 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 at least provides a kind of useful business to select.For this reason, one object of the present invention is to propose the robust method for estimating state based on maximal index absolute value target function that a kind of robust is good, computational efficiency is high (Maximum exponential absolute value state estimation, MEAV).

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

In one embodiment of the invention, the described robust 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 node injection is meritorious and idle, branch road meritorious and idle and node voltage amplitude measurement; X ∈ R ⁿfor state vector, comprise other each node phase angles except node voltage amplitude and balance node; h:R ⁿ→ R ^mfor the Nonlinear Mapping to measurement vector by state vector; r _ii the 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 as window width parameter.

In one embodiment of the invention, described step B comprises: introduce non-negative slack variable u, v ∈ R ^m, the described robust 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: making 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 < is ε, goes to step C7, otherwise enters step C3; Step C3: solve update equation, to complete 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: output optimal solution, finishes.

In one embodiment of the invention, described step C3 comprises: step C31: calculation perturbation parameter μ=ρ Gap/2m; Step C32: formation measurement equation and zero injecting power retrain corresponding Jacobian matrix and formation measurement equation and zero injecting power retrain the gloomy matrix in corresponding sea and wherein h (x), for state vector is to the mapping that measures vector, is measurement estimated value, 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 that corresponding weight is measured in i measurement; 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: solve 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 _i; And step C37: solve $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 robust method for estimating state based on maximal index absolute value target function of the embodiment of the present invention (Maximumexponential absolute value state estimation, MEAV) in estimation procedure, can effectively suppress to comprise a plurality of bad datas of consistency bad data, shown good robust, and there is 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.

Accompanying drawing explanation

Above-mentioned and/or additional aspect of the present invention and advantage accompanying drawing below combination obviously and is easily understood becoming the description of embodiment, wherein:

Fig. 1 is the schematic flow sheet of the robust 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 the qualification rate comparison of four kinds of state estimators of IEEE 300 systems.

Embodiment

Describe embodiments of the invention below in detail, the example of described embodiment is shown in the drawings, and wherein same or similar label represents same or similar element or has the element of identical or similar functions from start to finish.Below by the embodiment being described with reference to the drawings, be exemplary, 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 robust method for estimating state based on maximal index absolute value target function of the embodiment of the present invention comprises the following steps:

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

Particularly, the basic model of the MEAV that the present invention proposes 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 formula: z ∈ R ^mfor measuring vector, often comprise that 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 balance node phase angle) that comprises node voltage amplitude and phase angle; h:R ⁿ→ R ^mfor the Nonlinear Mapping to measurement vector by state 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 as window width parameter.

Step B: the robust state estimation basic model based on maximal index absolute value target function is introduced to auxiliary variable, and conversion obtains the robust 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 is more difficult.Model (1) ~ (3) can be converted into the equivalence model that can lead everywhere.

Introduce variable ξ ∈ R ^m, it is met

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

By formula (1) and (4), can be obtained, 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=-ξ _ii=1,…,m (5)

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

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

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

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

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

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

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

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

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

The equivalence model of the MEAV basic model that formula (1) ~ (3) provide 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 MEAV equivalence model, and this model everywhere continuous can be led, and can solve by the method based on gradient.Step C: utilize former-dual interior point, the described robust state estimation equivalence model based on maximal index absolute value target function is solved.

(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, and suitable use is former-and dual interior point solves.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 formula: λ ∈ R ^cand π, α, β ∈ R ^mfor Lagrange multiplier vector.

For obtaining optimal value, according to KKT condition, can obtain

$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 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, modern interior-point method is introduced disturbance parameter μ >0 formula (22), (23) is relaxed, 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 obtained 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 obtain

$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)$

By formula (33), (34) bring (31) into, (32) can obtain

$[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},$ By formula (35) and (36), can be obtained

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

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

In 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 obtained

Hdx+Qdπ＝γ (39)

In 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 obtain 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& 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]---\left(40\right)$

Solve formula (40) and can obtain dx, d λ and d π; By formula (37), (38), can obtain dv and du; By acquired results bring formula (33) into, (34) can obtain d α and d β, iteration is sustainable carrying out.

(2) solution procedure of MEAV equivalence model

Introduce MEAV equivalence model solve derivation after, inventor is summarized as follows solution procedure:

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

Particularly, make x ⁽⁰⁾∈ R ⁿrepresentative by all node voltage amplitudes and phase angle, formed flat starting state variable (except reference node phase angle); Select λ ⁽⁰⁾=π ⁽⁰⁾=0 and u ⁽⁰⁾, v ⁽⁰⁾, α ⁽⁰⁾, β ⁽⁰⁾>0, wherein λ ∈ R ^cand π, α, β ∈ R ^mfor 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< is ε, thinks and convergence can directly enter step C7; Otherwise for not restraining, enter step C3.

Step C3: solve update equation, to complete the correction to former variable and dual variable, obtain [dx ^td λ ^td π ^t] ^t, dv, du, d α and d β.

Particularly, first calculation perturbation parameter μ=ρ Gap/2m, then solves formula (40) and obtains [dx ^td λ ^td π ^t] ^t; Solve formula (37), (38) dv, du; Solve (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: output optimal solution, finishes.

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

Setting utilizes the performance of the MEAV of ieee standard system test based on 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 τ) in the result of calculating in trend.For voltage, measure, 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

Inventor compares MEAV of the present invention and other state estimators, tests the robust of MEAV.

(1) IEEE-14 system

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

The identification of table 1MEAV to IEEE14 system conformance bad data

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

As a comparison, first with the WLS of widely application, estimate, and with LNR, carry out the identification (brief note is WLS+LNR) of bad data.The result of identification is first: 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; Now find P ₂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.Visible, WLS+LNR can not identification consistency bad data.

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

(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 bad data ratio be 0% ~ 10% totally 11 kinds of situations test, under each ratio, include l-G simulation test 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.

From Fig. 2,3, when there is no 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 substantially constant, and phase angle error only slightly increases (even also like this when bad data ratio is increased to 10%), embodied good robustness.

Further, we to four kinds of state estimators (WLS+LNR, WLAV, MNMR and MEAV) qualification rate when the different bad data ratio test.As shown in Figure 4, wherein η and Ω represent respectively the ratio of bad data and the qualification rate of state estimation to the mean value of 50 tests.As seen from Figure 4, along with the rising of bad data ratio, the qualification rate of WLS+LNR declines rapidly; And the decline of the qualification rate of its excess-three kind robust estimator 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.As can be seen here, in above four kinds of state estimators, the robust of MEAV is best.

2. the comparison of computational efficiency

Inventor, in order to carry out efficiency comparison, tests four kinds of state estimator WLS, WLAV, MNMR and MEAV respectively under normal measurement condition, and wherein latter three kinds belong to robust state estimator.In test, WLS adopts Newton Algorithm, and other three kinds of state estimation adopt interior point method to solve; And MNMR adopts two-phase method, the first stage is carried out WLS estimation, and the initial value that second stage is estimated the estimated value of WLS as MNMR calculates.

Carry out altogether l-G simulation test 50 times, iterations and average computation during state estimation convergence are consuming time as shown in table 2.From table 2, in these four kinds of state estimators, the computational efficiency of WLS is the highest; And in rear three kinds of robust 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, the MEAV that the present invention proposes can effectively suppress to comprise a plurality of bad datas of consistency bad data in estimation procedure, has shown good robust, 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 example in conjunction with specific features, structure, material or the feature of this embodiment or example description.In this manual, the schematic statement of above-mentioned term is not necessarily referred to identical embodiment or example.And the specific features of description, structure, material or feature can be with suitable mode combinations in any one or more embodiment or example.

Although illustrated and described embodiments of the invention above, 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 can change above-described embodiment within the scope of the invention in the situation that not departing from principle of the present invention and aim, modification, replacement and modification.

Claims (5)

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

Steps A: the basic model of the robust state estimation based on maximal index absolute value target function is provided;

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

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

2. the robust method for estimating state based on maximal index absolute value target function as claimed in claim 1, is characterized in that, the described robust 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, wherein m represents to measure the dimension of vector, comprises that node injection is meritorious and idle, branch road meritorious and idle and node voltage amplitude measurement; X ∈ R ⁿfor state vector, wherein n represents the dimension of state vector, comprises other each node phase angles except node voltage amplitude and balance node; h:R ⁿ→ R ^mfor the Nonlinear Mapping to measurement vector by state vector; r _ii the element for residual error vector r; G (x): R ⁿ→ R ^cbe zero injecting power equality constraint, wherein c represents the number of zero injecting power; w _iit is the weight of i measurement amount.

3. the robust 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 robust 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 robust 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 x ⁽⁰⁾for flat starting state variable, introduce λ ∈ R ^cand π, α, β ∈ R ^mas Lagrange multiplier vector, and to make its initial value be λ ⁽⁰⁾=π ⁽⁰⁾=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< is ε, goes to step C7, otherwise enters step C3;

Step C3: solve update equation, to complete 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}{{\mathrm{dv}}_i}:{\mathrm{dv}}_i<0;-\frac{u_i}{{\mathrm{du}}_i}:{\mathrm{du}}_i<0),1\},$ ${\mathrm{\&theta;}}_D=0.9995\min \{\underset{i}{\min }(-\frac{{\mathrm{\&alpha;}}_i}{{\mathrm{d\&alpha;}}_i}:{\mathrm{d\&alpha;}}_i<0;-\frac{{\mathrm{\&beta;}}_i}{{\mathrm{d\&beta;}}_i}:{\mathrm{d\&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: output optimal solution, finishes.

5. the robust 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: formation measurement equation and zero injecting power retrain corresponding Jacobian matrix and formation measurement equation and zero injecting power retrain the gloomy matrix in corresponding sea wherein h (x), for state vector is to the mapping that measures vector, is measurement estimated value, 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 that corresponding weight is measured in i measurement;

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, wherein Q is the diagonal matrix of m dimension, its diagonal element Q _ii=-k _1i+ k _2i;

Step C36: solve 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 _i; And

Step C37: solve ${\mathrm{d\&alpha;}}_i={\mathrm{\&omega;}}_i{w_i}^{-1}{\mathrm{du}}_i+{\mathrm{\&omega;}}_i{w_i}^{-1}{\mathrm{dv}}_i-{\mathrm{d\&pi;}}_i+L_{v_i},$ And ${\mathrm{d\&beta;}}_i={\mathrm{\&omega;}}_i{w_i}^{-1}{\mathrm{du}}_i+{\mathrm{\&omega;}}_i{w_i}^{-1}{\mathrm{dv}}_i-{\mathrm{d\&pi;}}_i+L_{v_i}.$