CN103279676B - A kind of power system WLAV Robust filter method based on substitution of variable - Google Patents

Abstract

The invention discloses a kind of power system WLAV Robust filter method based on substitution of variable, the method is by introducing intermediate variable, it is achieved the non-linear measurement equation substep linearisation of state estimation.Interior point method is introduced the measurement equation after substep linearisation and carries out WLAV Robust filter, it is not necessary to form extra large gloomy matrix, reduce the dimension of coefficient matrix in Kuhn column gram (Karush Kuhn Tucker, KKT) equation group simultaneously.Multiple Simulation Example result verification set forth herein the validity of method, carries out Comparative result with traditional interior point method WLAV Robust filter, shows that method in this paper, in the advantage calculated on speed and robustness, has future in engineering applications.

Description

A kind of power system WLAV Robust filter method based on substitution of variable

Technical field

The present invention relates to a kind of power system WLAV Robust filter method based on substitution of variable, belong to Operation of Electric Systems With control technical field.

Background technology

As the core of EMS (Energy Management System, EMS), Power system state estimation By the process to raw data, it is thus achieved that the best estimate of quantity of state.Traditional weighted least-squares method (Weighted Least Squares, WLS) state estimation algorithm is simple, it is fast to calculate speed, but robustness is poor, without robustness, easily by umber of defectives According to impact.

Robust filter model has higher robustness, in the case of bad data is inevitable, it is possible to reduce bad data Impact, draw the best estimate of quantity of state.Compared with WLS estimation theory, the model that Robust filter is set up more conforms to number Border distribution pattern factually.It is current that weighting least absolute value (Weighted Least Absolute Values, WLAV) is estimated Study more Robust filter method.Research to WLAV robust state estimation is concentrated mainly on two aspects: one is to reduce meter Evaluation time, meets application request；Two is how to improve estimated accuracy, it is ensured that obtain optimal solution.

The advantages such as interior point method has initial value insensitive, and robustness is good, thus obtain extensively in Optimization Problems In Power Systems Application.In WLAV Robust filter, eliminate the absolute value amount in object function, by WLAV Robust filter by adding inequality constraints It is converted into the nonlinear optimal problem containing inequality constraints, and interior point method is the effective side of one processing problems Method.In state estimation, measurement equation is nonlinear function formula, and interior point method, when processing WLAV Robust filter, needs to form sea Gloomy matrix, adds programming difficulty, have impact on calculating speed, limits the engineer applied of WLAV Robust filter.Herein by change Amount replacement, is converted to system of linear equations by non-linear measurement equation group, then utilizes interior point method to process linear optimization problem, the party Method avoids the calculating of extra large gloomy matrix, solves above-mentioned problem, has engineer applied and is worth.

Summary of the invention

Goal of the invention: the present invention proposes a kind of power system WLAV Robust filter method based on substitution of variable, simplifies Calculating process, improves calculating speed.

Technical scheme: the technical solution used in the present invention is a kind of power system WLAV Robust filter based on substitution of variable Method, comprises the following steps:

1) network parameter and the measurement of power system are obtained；

2) initialize；

3) intermediate variable y is set:

Y=[M L O]^T

In formula: M_ij=U_iU_jcosθ_ij,L_ij=U_iU_jsinθ_ij,

4) intermediate variable utilizing step 3) set up following distribution estimate model:

Z=Ay+ ε_z

Y=Bx+ ε_x；

5) WLAV robust state estimation model is set up:

$\left\{\begin{array}{cc}\min .& w^T(l+u)\\ s.t.& g\left(x\right)=0\\ & z-h\left(x\right)+l-u=0\end{array}\right.$

6) Lagrangian is set up according to interior point method:

L=w^T(l+u)-η^T[z-h(x)+l-u]-α^Tl-β^Tu-λ^Tg(x)；

7) first equation of substep state estimation model is combined WLAV and estimates model, form corresponding Lagrangian letter Number, lists Ku En-Tucker type, carries out interior point method WLAV Robust filter and calculates, and iteration solves intermediate variable y；

8) the intermediate variable y obtained by step 7) is brought into second side of substep state estimation model in step 4) Journey, is carried out second equation and operation that first equation is same, solves quantity of state, exports result.

As a further improvement on the present invention, in described step 7), Ku En-Tucker type is according to Taylor expansion, and retains one Rank item, it is thus achieved that update equation is as follows:

(▽²h(x)η-▽²G (x) λ) dx+ h (x) d η-g (x) d λ=-L_x

-d η-d α=-L_l

D η-d α=-L_u

H (x) dx-dl+du=-L_η

G (x)=-L_λ

$\mathrm{Adl}+\mathrm{Ld\α}=-L_{\mathrm{\α}}^{\mathrm{\μ}}$

$\mathrm{Bdu}+\mathrm{Ud\β}=-L_{\mathrm{\β}}^{\mathrm{\μ}}$

Making η=0, α=β=w, then Ll=Lu=0, substituting into Ku En-Plutarch equation can obtain:

$\left\{\begin{array}{c}\mathrm{d\α}=-\mathrm{d\η}\\ \mathrm{d\β}=\mathrm{d\η}\\ \mathrm{dl}=-A^{-1}(-\mathrm{Ld\η}+L_{\mathrm{\α}}^{\mathrm{\μ}})\\ \mathrm{du}=-B^{-1}(-\mathrm{Ud\η}+L_{\mathrm{\β}}^{\mathrm{\μ}})\end{array}\right.$

And then can obtain iterative formula as follows:

$\left[\begin{array}{ccc}\hat{H}& -\▿h\left(x\right)& \▿g\left(x\right)\\ -\▿h{\left(x\right)}^T& S& 0\\ \▿g{\left(x\right)}^T& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{d\η}\\ \mathrm{d\λ}\end{array}\right]=\left[\begin{array}{c}-L_x\\ \mathrm{\γ}\\ -g\left(x\right)\end{array}\right]$

In formula: $\left\{\begin{array}{c}\hat{H}={\▿}^{2}g{\left(x\right)\mathrm{\λ}-\▿}^{2}h\left(x\right)\mathrm{\η}\\ \mathrm{\γ}=-L_{\mathrm{\η}}+A^{-1}{L}_{\mathrm{\α}}^{\mathrm{\μ}}-B^{-1}{L}_{\mathrm{\β}}^{\mathrm{\μ}}\end{array}\right.,$ D η back substitution by acquisition obtains the correction of its dependent variable.

Existing power system branch power equation expression formula is:

$\left\{\begin{array}{c}{P}_{\mathrm{ij}}=U_i^{2}g-U_i{U}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}g+U_i{U}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}b\\ Q_{\mathrm{ij}}=-U_i^2(b+y_c)+U_i{U}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}b-U_i{U}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}g\end{array}\right.$

Using after intermediate variable y, branch power equation can turn to following form:

P_ij=O_ig-M_ijg+L_ijb

Q_ij=-O_i(b+y_c)+M_ijb-L_ijg

Can be seen that in above formula, branch power and variable y are linear relationship, if form voltage measurement being set to square, The most all measurement and y are linear relationship.

If electrical network has b bar branch road, it is known that the dimension of y is N+2b, makes m > n, then the dimension of intermediate variable is less than measurement Dimension, sets up following measurement-intermediate variable functional equation:

Z=Ay+ ε_z

Calculated by a step and i.e. can get the weighted least-square solution of y:

Y=(A^TWA)^-1A^TWz

Beneficial effect: the present invention is by introducing intermediate variable, it is achieved the non-linear measurement equation substep of state estimation is linear Change.Interior point method is introduced the measurement equation after substep linearisation and carries out WLAV Robust filter, it is not necessary to form extra large gloomy matrix, subtract simultaneously The dimension of coefficient matrix in little Ku En-Plutarch (Karush-Kuhn-Tucker, KKT) equation group.Multiple Simulation Example results Demonstrate the validity of the method for set forth herein, carry out Comparative result with traditional interior point method WLAV Robust filter, show to carry herein The method gone out, in the advantage calculated on speed and robustness, has future in engineering applications.

Accompanying drawing explanation

Fig. 1 is the circulation figure of the present invention；

Fig. 2 is IEEE-14 node structure schematic diagram；

Fig. 3 is IEEE-57 node structure schematic diagram；

Fig. 4 is IEEE-118 node structure schematic diagram；

Fig. 5 is the voltage magnitude estimated accuracy comparison diagram of the present invention and IEEE-57 standard example；

Fig. 6 is the voltage phase angle estimated accuracy comparison diagram of the present invention and IEEE-57 standard example.

Detailed description of the invention

Below in conjunction with the accompanying drawings and specific embodiment, it is further elucidated with the present invention, it should be understood that these embodiments are merely to illustrate The present invention rather than limit the scope of the present invention, after having read the present invention, each to the present invention of those skilled in the art The amendment planting equivalents all falls within the application claims limited range.

As it is shown in figure 1, first the first step of the present invention obtains network parameter and the measurement of power system.Network parameter bag Include: bus numbering, title, compensation electric capacity, the branch road number of transmission line of electricity, headend node and endpoint node numbering, series resistance, string Connection reactance, shunt conductance, shunt susceptance, transformer voltage ratio and impedance.Measurement includes that node voltage amplitude, branch road head end are meritorious Power and reactive power, branch road end active power and reactive power.

After obtaining above-mentioned parameter, start program initializes, including iteration precision, maximum iteration time, former antithesis The slack variable of interior point method, Lagrange multiplier and penalty factor initial value, form bus admittance matrix.

General, the measurement equation of Power system state estimation is:

Z=h (x)+ε

In formula: x is quantity of state (dimension n=2N-1, N are nodes)；Z be measurement (dimension m, m > n)；H is that m ties up non-thread Property measure function；ε is that m ties up error in measurement.

Present invention introduces intermediate variable y:y=[M L O]^T, in formula: M_ij=U_iU_jcosθ_ij,L_ij=U_iU_jsinθ_ij,Utilize intermediate variable y to set up following distribution and estimate model:

Z=Ay+ ε_z

Y=Bx+ ε_x。

The object function of WLAV state estimation is:

$\min .\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{-1}\left|{\mathrm{\ϵ}}_i\right|$

In formula: ω_iAnd ε_iIt is respectively weight vectors and the i-th component of residual vector.

Above-mentioned object function is turned to the optimization problem containing constraints:

$\left\{\begin{array}{cc}\min .& w^T\mathrm{\ϵ}\\ s.t.& g\left(x\right)=0\\ & z-h\left(x\right)-\mathrm{\ϵ}\≤0\\ & -z+h\left(x\right)-\mathrm{\ϵ}\≤0\end{array}\right.$

Slack variable is added in the inequality constraints of above formula, is translated into equality constraint:

$\left\{\begin{array}{c}z-h\left(x\right)-\mathrm{\ϵ}+2l=0,l\≥0\\ -z+h\left(x\right)-\mathrm{\ϵ}+2u=0,u\≥0\end{array}\right.$

In above formula two equality constraints are carried out plus and minus calculation turn to:

$\left\{\begin{array}{c}\mathrm{\ϵ}=l+u\\ z-h\left(x\right)+l-u=0\end{array}\right.$

New WLAV estimates to be modeled as:

$\left\{\begin{array}{cc}\min .& w^T(l+u)\\ s.t.& g\left(x\right)=0\\ & z-h\left(x\right)+l-u=0\end{array}\right.$

The Lagrangian of structure above formula is:

L=w^T(l+u)-η^T[z-h(x)+l-u]-

α^Tl-β^Tu-λ^Tg(x)

In formula: η and λ is Lagrange multiplier, l and u is slack variable, α and β is barrier function.

On Ku En-Plutarch equation be:

$\left\{\begin{array}{c}{L}_x={\mathrm{\η}}^T\▿h\left(x\right)-{\mathrm{\λ}}^T\▿g\left(x\right)=0\\ L_l=w-\mathrm{\η}-\mathrm{\α}=0\\ L_u=w+\mathrm{\η}-\mathrm{\β}=0\\ L_{\mathrm{\η}}=-z+h\left(x\right)-l+u=0\\ L_{\mathrm{\λ}}=g\left(x\right)=0\\ L_{\mathrm{\α}}^{\mathrm{\μ}}=\mathrm{ALe}-\mathrm{\μe}=0\\ L_{\mathrm{\β}}^{\mathrm{\μ}}=\mathrm{BUe}-\mathrm{\μe}=0\end{array}\right.$

In formula: (A, L, B, U) is so that (α, l, β u) are the diagonal matrix of diagonal element respectively；μ is Discontinuous Factors.

By Ku En-Plutarch equation Taylor expansion, retain single order item, it is thus achieved that update equation is as follows:

(▽²h(x)η-▽²G (x) λ) dx+ h (x) d η-g (x) d λ=-L_x

-d η-d α=-L_l

D η-d α=-L_u

H (x) dx-dl+du=-L_η

G (x)=-L_λ

$\mathrm{Adl}+\mathrm{Ld\α}=-L_{\mathrm{\α}}^{\mathrm{\μ}}$

$\mathrm{Bdu}+\mathrm{Ud\β}=-L_{\mathrm{\β}}^{\mathrm{\μ}}$

Making η=0, α=β=w, then Ll=Lu=0, substituting into Ku En-Plutarch equation can obtain:

$\left\{\begin{array}{c}\mathrm{d\α}=-\mathrm{d\η}\\ \mathrm{d\β}=\mathrm{d\η}\\ \mathrm{dl}=-A^{-1}(-\mathrm{Ld\η}+L_{\mathrm{\α}}^{\mathrm{\μ}})\\ \mathrm{du}=-B^{-1}(-\mathrm{Ud\η}+L_{\mathrm{\β}}^{\mathrm{\μ}})\end{array}\right.$

And then can obtain iterative formula as follows:

$\left[\begin{array}{ccc}\hat{H}& -\▿h\left(x\right)& \▿g\left(x\right)\\ -\▿h{\left(x\right)}^T& S& 0\\ \▿g{\left(x\right)}^T& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{dx}\\ \mathrm{d\η}\\ \mathrm{d\λ}\end{array}\right]=\left[\begin{array}{c}-L_x\\ \mathrm{\γ}\\ -g\left(x\right)\end{array}\right]$

In formula: $\left\{\begin{array}{c}\hat{H}={\▿}^{2}g\left(x\right)\mathrm{\λ}-{\▿}^{2}h\left(x\right)\mathrm{\η}\\ \mathrm{\γ}=-L_{\mathrm{\η}}+A^{-1}{L}_{\mathrm{\α}}^{\mathrm{\μ}}-B^{-1}{L}_{\mathrm{\β}}^{\mathrm{\μ}}\end{array}\right.,$ D η back substitution by acquisition obtains the correction of its dependent variable.

Existing power system branch power equation expression formula is:

$\left\{\begin{array}{c}{P}_{\mathrm{ij}}=U_i^{2}g-U_i{U}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}g+U_i{U}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}b\\ Q_{\mathrm{ij}}=-U_i^2(b+y_c)+U_i{U}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}b-U_i{U}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}g\end{array}\right.$

Using after intermediate variable y, branch power equation can turn to following form:

P_ij=O_ig-M_ijg+L_ijb

Q_ij=-O_i(b+y_c)+M_ijb-L_ijg

Can be seen that in above formula, branch power and variable y are linear relationship, if form voltage measurement being set to square, The most all measurement and y are linear relationship.

If electrical network has b bar branch road, it is known that the dimension of y is N+2b, makes m > n, then the dimension of intermediate variable is less than measurement Dimension, sets up following measurement-intermediate variable functional equation:

Z=Ay+ ε_z

Calculated by a step and i.e. can get the weighted least-square solution of y:

Y=(A^TWA)^-1A^TWz

Y is done suitable conversion as follows:

$y^{\′}=\left[\begin{array}{c}\ln (M_{\mathrm{ij}}^2+N_{\mathrm{ij}}^2)\\ \arctan (M_{\mathrm{ij}}/N_{\mathrm{ij}})\\ \ln \left(O_i\right)\end{array}\right]=\left[\begin{array}{c}2\ln \left(U_i\right)+2\ln \left(U_i\right)\\ {\mathrm{\θ}}_i-{\mathrm{\θ}}_j\\ 2\ln \left(U_i\right)\end{array}\right]$

Voltage magnitude in quantity of state is taken its logarithm simultaneously, is turned to:

$x^{\′}=\left[\begin{array}{c}\ln \left(U\right)\\ \mathrm{\θ}\end{array}\right]$

Y ' can be obtained, and x ', function expression be:

$y^{\′}=\left[\begin{array}{ccc}I& 0& 0\\ C^T& 0& 0\\ 0& C_R^T& 0\end{array}\right]x^{\′}$

$={\mathrm{Bx}}^{\′}$

In formula: I is unit battle array；C is branch node incidence matrix, C_RFor deleting the node branch road association square of balance node Battle array.Then y ' and x ' is linear relationship, and the least square solution of x ' is:

X '=(B^TW_y′B)B^TW_y′y′

After obtaining x ', its logarithm voltage segment is taken the logarithm and can obtain the estimate of voltage magnitude.

To z=Ay+ ε_zCarry out interior point method WLAV Robust filter with y '=Bx ' to calculate, owing to two formulas are all systems of linear equations, Belong to linear programming problem, it is to avoid the calculating of extra large gloomy matrix, reduce KKT equation group coefficient matrix dimension simultaneously, simplify Matrix operation, improves calculating speed.

The present invention uses the IEEE-14 node shown in Fig. 2 to Fig. 4, IEEE-57 node, the standard of IEEE-118 node to calculate Example, in order to contrast the efficiency of algorithm of two kinds of methods, simulation result is as shown in the table:

1 two kinds of method computational efficiencies of table compare

As can be seen from Table 1, substitution of variable interior point method WLAV robust state estimation is on iterations and traditional interior point Method WLAV Robust filter is close, and a moment is for 1-2 time；For calculating the time, the time of substitution of variable method is shorter, has by contrast Preferably practical value.

Meanwhile, present invention employs the standard example of IEEE-57 node, in order to verify the validity of the method for set forth herein, Estimated accuracy being carried out emulation compare, result is as shown in Figure 5 and Figure 6.

Claims (1)

1. a power system WLAV Robust filter method based on substitution of variable, it is characterised in that comprise the following steps:

1) network parameter and the measurement of power system are obtained；

2) measurement equation of Power system state estimation is:

Z=h (x)+ε

In formula: x is quantity of state, dimension n=2N-1, N are nodes, and z is measurement, dimension m, m > n, h be m tie up non-linear measurement Function, ε is that m ties up error in measurement；

3) object function of WLAV state estimation is:

In formula: ω_iAnd ε_iIt is respectively weight vectors and the i-th component of residual vector；

4) by step 3) in object function be converted into new WLAV robust state estimation model:

5) Lagrangian is set up according to interior point method:

L=w^T(l+u)-η^T[z-h(x)+l-u]-α^Tl-β^Tu-λ^Tg(x)；

In formula: η and λ is Lagrange multiplier, l and u is slack variable, α and β is barrier function；

6) by step 5) in Lagrangian be converted into Ku En-Plutarch equation:

In formula: (A, L, B, U) is so that (α, l, β u) are the diagonal matrix of diagonal element respectively；μ is Discontinuous Factors；

7) by Ku En-Plutarch equation Taylor expansion, single order item is retained, it is thus achieved that update equation:

-d η-d α=-L_l

D η-d α=-L_u

Making η=0, α=β=w, then Ll=Lu=0, substituting into Ku En-Plutarch equation can obtain:

And then can obtain iterative formula as follows:

In formula:D η back substitution by acquisition obtains the correction of its dependent variable；

8) take intermediate variable and be set to following form:

Y=[M L O]^T

In formula: M_ij=U_iU_j cosθ_ij,L_ij=U_iU_j sinθ_ij,

Utilize intermediate variable y set up following distribution estimate model:

Z=Ay+ ε_z

Y=Bx+ ε_x,

Existing power system branch power equation expression formula is:

Using after intermediate variable, branch power equation can turn to following form:

P_ij=O_ig-M_ijg+L_ijb

Q_ij=-O_i(b+y_c)+M_ijb-L_ij；

9) set electrical network and have b bar branch road, it is known that the dimension of y is N+2b, makes m n, then the dimension of intermediate variable is tieed up less than measurement Number, sets up following measurement-intermediate variable functional equation:

Z=Ay+ ε_z

Calculated by a step and i.e. can get the weighted least-square solution of y:

Y=(A^TWA)^-1A^TWz

Y is done suitable conversion as follows:

Voltage magnitude in quantity of state is taken its logarithm simultaneously, is turned to:

The function expression that can obtain y ' ' and x ' ' is:

In formula: I is unit battle array；C is branch node incidence matrix, C_RFor deleting the node branch road incidence matrix of balance node, then y ' Being linear relationship with x ', the least square solution of x ' is:

X '=(B^TW_y′B)B^TW_y′y′

After obtaining x ', its logarithm voltage segment is taken the logarithm and can obtain the estimate of voltage magnitude.