CN109255541A

CN109255541A - A kind of power distribution network robust state estimation method the sum of multiplied based on least square and one

Info

Publication number: CN109255541A
Application number: CN201811097154.1A
Authority: CN
Inventors: 陈艳波; 沈玉兰; 张璞; 李翔宇; 方恒福; 周猛; 郎燕生
Original assignee: China Electric Power Research Institute Co Ltd CEPRI; North China Electric Power University; State Grid Economic and Technological Research Institute
Current assignee: China Electric Power Research Institute Co Ltd CEPRI; North China Electric Power University; State Grid Economic and Technological Research Institute
Priority date: 2018-09-20
Filing date: 2018-09-20
Publication date: 2019-01-22
Anticipated expiration: 2038-09-20
Also published as: CN109255541B

Abstract

The invention belongs to dispatching automation of electric power systems technical fields, a kind of more particularly to power distribution network robust state estimation method the sum of multiplied based on least square and one, comprising: step A: the power distribution network robust state estimation basic model that the building of vector sum state vector the sum of is multiplied based on least square and one is measured with power distribution network node；Step B: auxiliary variable is introduced by the basic model and converts to obtain guidable equivalence model everywhere；Step C: primal-dual interior point algorithm parity price model solution is utilized.Multiple bad datas including consistency bad data can be effectively suppressed in the present invention in estimation procedure, it is shown that good Robustness least squares, and there is very high computational efficiency, it is extremely suitable for practical engineering application.

Description

A kind of power distribution network robust state estimation method the sum of multiplied based on least square and one

Technical field

The invention belongs to dispatching automation of electric power systems technical fields, more particularly to a kind of least square and one that is based on to multiply it The power distribution network robust state estimation method of sum.

Background technique

Power system state estimation is basis and the core of Energy Management System.Present almost each large-scale control centre It is assembled with state estimator, state estimation has become the foundation stone of electric power netting safe running.Shape is put forward for the first time from 1970 foreign scholars Since state is estimated, people have had more than 40 years history to the research of state estimation and application, have emerged during this various The method for estimating state of various kinds.

Currently, the state estimation being at home and abroad most widely used is weighted least-squares method (Weighted least squares,WLS).WLS model simple is solved and is easy, but its Robustness least squares is very poor.In order to enhance Robustness least squares, it is general there are two types of Method.The first is that bad data recognition link, such as maximum regularization residual test method (LNR) are added after WLS estimation Or estimation discrimination method etc.；Another kind is using robust state estimation method.Currently, the robust shape that domestic and foreign scholars have proposed State estimation method (Robust state estimation) includes that weighting least absolute value estimates (Weighted least Absolute value, WLAV), Non quadratic criteria method (QL, QC etc.), the state estimation for being up to qualification rate target (Maximum normal measurement rate, MNMR) and exponential type objective function state estimation (Maximum Exponential square, MES) etc..But the estimation performance of these robust state estimation methods is still to be improved.

Summary of the invention

In view of the above-mentioned problems, being estimated the invention proposes a kind of based on the power distribution network robust state that least square and one the sum of multiply (Weighter Least Square and Absolute Value, the WLSAV) method of counting, comprising:

Step A: the distribution that the building of vector sum state vector the sum of is multiplied based on least square and one is measured with power distribution network node Net robust state estimation basic model；

Step B: auxiliary variable is introduced by the basic model and converts to obtain guidable equivalence model everywhere；

Step C: primal-dual interior point algorithm parity price model solution is utilized.

The basic model are as follows:S.t.g (x)=0, r=z-h (x)

Wherein: J (x) is objective function, z ∈ R^mTo measure vector, including the active and idle, branch of node injection it is active and Idle and node voltage amplitude measures；x∈RⁿFor state vector, including other except node voltage amplitude and balance nodes Each node phase angle；h:Rⁿ→R^mFor by state vector to the Nonlinear Mapping for measuring vector；R^mM dimensional vector is represented, m is amount The number of measurement, RⁿM dimensional vector is represented, n is the number of state variable, r_iFor i-th of element of residual error vector r；w_iIt is i-th The weight of a measurement；y_iFor the weight of the corresponding first order of i-th of measurement；g(x):Rⁿ→R^cBe zero injecting power equation about Beam, R^cRepresent c dimensional vector；R represents residual vector, and h (x), which is represented, measures function.

The step 2 includes: to introduce non-negative slack variable u, v ∈ R^m, the equivalence model that converts are as follows:S.t.g (x)=0, z-h (x)-u+v=0, u, v >=0

Wherein, u_iFor non-zero auxiliary variable, v_iFor non-zero auxiliary variable.

The step C includes:

Step C1: enabling x is flat starting state variable；Select λ⁽⁰⁾=π⁽⁰⁾=0 and u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾,β⁽⁰⁾> 0, wherein λ ∈R^cAnd π, α, β ∈ R^mFor Lagrange multiplier vector, parameter upper right footmark indicates the number of iteration, and current is the 0th time；In order Heart parameter ρ ∈ (0,1) and convergence criterion ε=10^-3, set iteration count k=0；

Step C2: duality gap Gap=α is calculated^Tv+β^TU judges whether to restrain, if duality gap Gap < convergence criterion ε, C7 is then gone to step, C3 is otherwise entered step；

Step C3: solving update equation, to complete the amendment to former variable and dual variable, obtains [dx^T dλ^T dπ^T]^T, Dv, du, d α and d β；

Step C4: the amendment step-length θ of former problem and dual problem is calculated_PAnd θ_D, in which:

α_i Indicate the i-th dimension of vector α；β_iIndicate the i-th dimension of vector β.

Step C5: the variable of former problem and dual problem is corrected respectively are as follows: Parameter upper right footmark represents the number of iterations k；

Step C6: iteration count k=k+1 is enabled, C2 is entered step；

Step C7: output optimal solution.

The step C3 includes:

Step C31: calculation perturbation parameter μ=ρ Gap/2m, m is the number of measurement；

Step C32: forming measurement equation and zero injecting power constrains corresponding Jacobian matrixAndIt forms measurement equation and zero injecting power constrains corresponding Hessian matrixAndWherein h (x) For state vector to the mapping for measuring vector, as measurement estimated value, g (x)=0 is the constraint of zero injecting power；

Step C33: L is calculated_x=G^Tλ-H^Tπ, L_λ=g (x), L_π=z-h (x)-u+v, AndL_x、L_λ、L_π、For centre auxiliary Variable；

Step C34: calculating γ=z-h (x)-u+v+AA-BB, and wherein γ, AA, BB are intermediate auxiliary column vector, AA, BB ∈ R^m,AA_iFor intermediate auxiliary variable,BB_i For intermediate auxiliary variable,a_i、b_i、c_i、d_iFor intermediate auxiliary variable, p_i=2w_i, p_iFor in Between auxiliary variable；

Step C35: equation is solvedObtain [dx^T dλ^T dπ^T]^T, Q For the diagonal matrix of m row m column；

Step C36: dv is solved_i=k_1idπ_i+AA_iAnd du_i=k_2idπ_i+BB_i, wherein intermediate variable k_1i=a_iv_i-b_iu_i, Intermediate variable k_2i=c_iv_i-d_iu_i；

Step C37: it solvesAnd

Beneficial effects of the present invention: the present invention is can be effectively suppressed in estimation procedure including consistency bad data Multiple bad datas, it is shown that good Robustness least squares, and there is very high computational efficiency, it is extremely suitable for practical engineering application.

Detailed description of the invention

Fig. 1 is step figure of the invention.

Specific embodiment

With reference to the accompanying drawing, it elaborates to embodiment.

As shown in Figure 1, the power distribution network robust state estimation of the embodiment of the present invention the sum of multiplied based on least square and one (Weighter Least Square and Absolute Value, WLSAV) method includes the following steps:

Step A: power distribution network robust state estimation (the Weighter Least the sum of multiplied based on least square and one is provided Square and Absolute Value, WLSAV) method basic model.

Specifically, the basic model of WLSAV proposed by the present invention is as follows

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

R=z-h (x) (3)

In formula: z ∈ R^mIt often include that node injection is active and idle, branch is active and idle and node to measure vector Voltage magnitude measurement etc.；x∈RⁿIt is the state vector for including node voltage amplitude and phase angle (except balance nodes phase angle)；h:Rⁿ →R^mFor by state vector to the Nonlinear Mapping for measuring vector；r_iIt is i-th of element of residual error vector r；g(x):Rⁿ→R^cFor Zero injecting power equality constraint；w_iFor the weight of i-th of measurement.

In order to facilitate solution, minimum target function is changed to maximum target function:

Step B: auxiliary variable is introduced to WLSAV basic model, transformation obtains the equivalence model of WLSAV basic model.

Specifically, it although the objective function everywhere continuous of WLSAV basic model, can not lead at 0, thus directly asks It solves relatively difficult.Model (1)~(3) can be converted to a guidable equivalence model everywhere.

Introduce variable ξ ∈ R^m, make its satisfaction

|r_i|≤ξ_iI=1 ..., m (4)

It can be obtained by formula (1) and (4),Maximum is equivalent to It is maximum.

Introduce non-negative slack variable l, k ∈ R^m, converting two equality constraints for inequality (4) is

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, make its satisfaction

u_i=l_i/ 2i=1 ..., m (7)

v_i=k_i/ 2i=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 constraint condition of equal value is

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

The equivalence model for the WLSAV basic model that then formula (1)~(3) provide is

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

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

u,v≥0 (15)

Model (12)~(15) are the equivalence model of WLSAV basic model, which can lead, can be with being based on The method of gradient solves.

Step C: utilizing primal-dual interior point algorithm, to the state estimation the sum of multiplied based on least square and one The equivalence model of (Weighter Least Square and Absolute Value, WLSAV) model solves.

(1) method for solving of WLSAV equivalence model

Equivalence model (12)~(15) for noticing WLSAV are the optimizations containing equality constraint and inequality constraints Problem is suitable for being solved with primal-dual interior point algorithm.In order to enable those skilled in the art to better understand the present invention, it gives first Detailed derivation process is as follows out:

Introduce Lagrangian

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

It can be obtained to obtain optimal value according to KKT condition

In formula:

Effectively to solve problem above, modern interior point method introduces disturbance parameter μ > 0 pair formula (22), (23) relax, from And it obtains

Above equation can be obtained by Newton Algorithm

Gdx=-L_λ (27)

- Hdx-du+dv=-L_π (28)

Wherein, p_i=2w_i；

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

Bring formula (33), (34) into (31), (32) can obtain

It enablesIt can be obtained by formula (35) and (36)

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,

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

Hdx+Qd π=γ (39)

In formula: Q R^m×mDiagonal matrix, diagonal element Q_ii=-k_1i+k_2i；γ=z-h (x)-u+v+AA-BB, AA, BB∈R^m, AA_i,BB_iWith it is same in formula (37), (38).

According to formula (39), (26) and (27), can obtain update equation is

Solution formula (40) can obtain dx, d λ and d π；Dv and du can be obtained by formula (37), (38)；By acquired results bring into formula (33), (34) d α and d β can be obtained, then iteration, that is, sustainable progress.

(2) solution procedure of WLSAV equivalence model

After introducing the solution derivation process of WLSAV equivalence model, solution procedure is summarized as follows by inventor:

Step C1: being initialized, and enabling x is flat starting state variable；Select λ⁽⁰⁾=π⁽⁰⁾=0 and u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾, β⁽⁰⁾> 0；It enables Center Parameter ρ ∈ (0,1) and determines convergence criterion value, and setting iteration count is zero.

Specifically, x is enabled⁽⁰⁾∈RⁿRepresent the flat starting state variable (ginseng being made of all node voltage amplitudes and phase angle Except examining node phase angle)；Select λ⁽⁰⁾=π⁽⁰⁾=0 and u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾,β⁽⁰⁾> 0, wherein λ ∈ R^cAnd π, α, β ∈ R^mFor glug Bright day multiplier vector, m is the number of measurement, and c is the number of zero injecting power constraint；It enables Center Parameter ρ ∈ (0,1) and receives Hold back criterion ε=10^-3, set iteration count k=0.

Step C2: duality gap Gap=α is calculated^Tv+β^TU judges whether to restrain.Specifically, if Gap < ε, then it is assumed that receive It holds back, step C7 can be directly entered；Otherwise not restrain, C3 is entered step.

Step C3: solving update equation, to complete the amendment to former variable and dual variable, obtains [dx^T dλ^T dπ^T]^T, Dv, du, d α and d β.

Specifically, then calculation perturbation parameter μ=ρ Gap/2m first solves formula (40) and obtains [dx^T dλ^T dπ^T]^T；It asks Solution formula (37), (38) obtain dv, du；Solve (33), (34) obtain d α, d β.

Step C5: the variable of former problem and dual problem is corrected respectively are as follows:

Step C6: iteration count k=k+1 is enabled, C2 is entered step；And

Step C7: output optimal solution terminates.

In order to enable those skilled in the art to better understand the present invention and understand the present invention compared with the advantages of the prior art, Shen It asks someone further to be illustrated in conjunction with specific embodiments.

Setting utilizes the performance of WLSAV of the ieee standard system test based on primal-dual interior point algorithm.Test uses full dose Survey, measuring value by the result of Load flow calculation Additive White Noise (mean value 0, standard deviation τ) obtain.For voltage It surveys, takes τ_V=0.005p.u.；For power measurement, τ is taken_PQ=1MW/MVar.Test environment is PC machine, and CPU is Intel (R) Core (TM) i3M370, dominant frequency 2.40GHz, memory 2.00GB.

1. the comparison of robustness

WLSAV of the invention is compared by inventor with other state estimators, to test the Robustness least squares of WLSAV.

4 consistency bad data (P are set in IEEE-39 system_1-2、Q_1-2、P₁、Q₁).Set bad measuring value And the right value of measurement is as shown in table 1.

Identification of 1 WLSAV of table to 39 system conformance bad data of IEEE

As a comparison, estimated first with widely used WLS, and (be abbreviated as with the identification that LNR carries out bad data WLS+LNR).The result recognized for the first time are as follows: the standardized residual of 10 measurements is greater than threshold value (3.0), this 10 measurements It is considered as suspicious data；Wherein the maximum measurement of standardized residual is P_2-1, WLS is reruned after leaving out the measurement；At this time It was found that P₂Standardized residual it is maximum.Above procedure recycles 4 times, 4 good measurements suspicious data is mistakenly considered by LNR and Left out, but really bad data still has.As it can be seen that WLS+LNR cannot recognize consistency bad data.

Estimated result using WLSAV method is as shown in table 1.It can be found that bad even if there are consistency in measurement Data, estimated value and the true value of WLSAV can also coincide well.Also indicate that WLSAV exists in the test of many times of IEEE other systems It can inhibit bad data during estimation automatically, there are good Robustness least squares.

2. the comparison of computational efficiency

Inventor in order to carry out efficiency comparison, under the conditions of normal measure respectively to four kinds of state estimator WLS, WLAV, MNMR and WLSAV are tested, wherein latter three kinds belong to robust state estimator.In test, WLS is asked using Newton method Solution, other three kinds of state estimations are solved using interior point method；And MNMR uses two-phase method, i.e. the first stage carries out WLS estimation, the Two-stage calculates the estimated value of WLS as the MNMR initial value estimated.

50 l-G simulation tests are carried out altogether, and the number of iterations and average computation when state estimation restrains are time-consuming as shown in table 2. As can be seen from Table 2, in these four state estimators, the computational efficiency highest of WLS；And in rear three kinds of robust state estimators, The computational efficiency highest of WLSAV；And increase with the increase of system scale, the number of iterations and calculating time-consuming of WLSAV It is very slow, thus WLSAV is suitable for the estimation of actual large scale system.

The number of iterations of 2 four kinds of state estimators of table and calculating are time-consuming

Exist in conclusion WLSAV proposed by the present invention can be effectively suppressed in estimation procedure including consistency bad data Interior multiple bad datas, it is shown that good Robustness least squares, and there is very high computational efficiency, it is extremely suitable for Practical Project and answers With.

This embodiment is merely preferred embodiments of the present invention, but scope of protection of the present invention is not limited thereto, In the technical scope disclosed by the present invention, any changes or substitutions that can be easily thought of by anyone skilled in the art, It should be covered by the protection scope of the present invention.Therefore, protection scope of the present invention should be with scope of protection of the claims Subject to.

Claims

1. a kind of power distribution network robust state estimation method the sum of multiplied based on least square and one characterized by comprising

Step A: it is anti-that the power distribution network that the building of vector sum state vector the sum of is multiplied based on least square and one is measured with power distribution network node Poor state estimation basic model；

2. method according to claim 1, which is characterized in that the basic model are as follows: S.t.g (x)=0, r=z-h (x)；

Wherein: J (x) is objective function, z ∈ R^mFor measure vector, including node injection it is active and idle, branch is active and idle And node voltage amplitude measures；x∈RⁿFor state vector, including except node voltage amplitude and balance nodes other are each Node phase angle；h:Rⁿ→R^mFor by state vector to the Nonlinear Mapping for measuring vector；R^mM dimensional vector is represented, m is measurement Number, RⁿN dimensional vector is represented, n is the number of state variable, r_iFor i-th of element of residual error vector r；w_iFor i-th of amount Measure the weight of corresponding quadratic term；y_iFor the weight of the corresponding first order of i-th of measurement；g(x):Rⁿ→R^cIt is zero injection function Rate equality constraint, R^cRepresent c dimensional vector；R represents residual vector, and h (x), which is represented, measures function.

3. method according to claim 2, which is characterized in that the step 2 includes: to introduce non-negative slack variable u, v ∈ R^m, Convert the obtained equivalence model are as follows:S.t.g (x)=0, z-h (x)-u + v=0, u, v >=0

4. method according to claim 3, which is characterized in that the step C includes:

Step C1: enabling x is flat starting state variable；Select λ⁽⁰⁾=π⁽⁰⁾=0 and u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾,β⁽⁰⁾> 0, wherein λ ∈ R^cAnd π,α,β∈R^mFor Lagrange multiplier vector, parameter upper right footmark indicates the number of iteration, and current is the 0th time；Enable Center Parameter ρ ∈ (0,1) and convergence criterion ε=10^-3, set iteration count k=0；

Step C2: duality gap Gap=α is calculated^Tv+β^TU judges whether to restrain, if duality gap Gap < convergence criterion ε, turns Step C7, otherwise enters step C3；

α_iIndicate the i-th dimension of vector α；β_iIndicate the i-th dimension of vector β；

Step C5: the variable of former problem and dual problem is corrected respectively are as follows: Parameter upper right footmark indicates the number of iterations k；

Step C6: iteration count k=k+1 is enabled, C2 is entered step；

Step C7: output optimal solution.

5. method according to claim 4, which is characterized in that the step C3 includes:

Step C33: L is calculated_x=G^Tλ-H^Tπ, L_λ=g (x), L_π=z-h (x)-u+v, AndL_x、L_λ、L_π、It is intermediate auxiliary Help variable；

Step C34: calculating γ=z-h (x)-u+v+AA-BB, and wherein γ, AA, BB are intermediate auxiliary column vector, AA, BB ∈ R^m,AA_iFor intermediate auxiliary variable,BB_iFor in Between auxiliary variable,a_i、b_i、c_i、d_iFor intermediate auxiliary variable, p_i=2w_i, p_iIt is intermediate auxiliary Help variable；

Step C35: equation is solvedObtain [dx^T dλ^T dπ^T]^T, Q is m row m The diagonal matrix of column, diagonal element Q_ii=-k_1i+k_2i, intermediate variable k_1i=a_iv_i-b_iu_i, intermediate variable k_2i=c_iv_i- d_iu_i；

Step C36: dv is solved_i=k_1idπ_i+AA_iAnd du_i=k_2idπ_i+BB_i；

Step C37: it solvesAnd