CN103886193A - Fuzzy self-adaptation robust estimation method of electric power system - Google Patents

Abstract

The invention discloses a fuzzy self-adaptation robust estimation method of an electric power system, and provides a fuzzy self-adaptation robust estimation algorithm with uncertainty metering and measuring weight of the electric power system. According to the algorithm, a measuring point inferior fuzzy membership conception is put forward, the standard deviation of a measuring point is corrected in an on-line mode according to the membership of the measuring point, and gross error self adaptation is achieved. A modern interior point method is used for solving, convergence is good, results cannot be affected by an initial value easily, the solving large-scale optimization problem is solved, bad data are effectively recognized, residual error contamination and residual error flooding are avoided, excellent robust estimation performance is achieved to the quantity of state and quantity measuring, and a good engineering application background is achieved.

Description

A kind of electric system fuzzy self-adaption robust method of estimation

Technical field

The present invention relates to a kind of electric system fuzzy self-adaption robust method of estimation, belong to power system monitoring, analysis and control technical field.

Background technology

State estimation is the basis of Power System Analysis and control, it is energy management system key foundation module, its main task is the real-time information providing according to data acquisition and monitoring system (SCADA), estimate the each busbar voltage of electrical network (amplitude, phase angle) and power, and comprise bad data detection and identification function.

Traditional state estimation algorithm is founded in 1970 by people such as Schweppe.According to the difference that solves objective function, state estimator comprises: weighted least-squares (WLS), Non quadratic criteria (QC, QL etc.), weighting least absolute value (WLAV), minimum median square (LMS), minimum truncated side (LTS) etc.At present most widely used is WLS estimation criterion, advantage is that model is simple, calculated amount is little, to desirable normal distribution measurement amount, can reach optimal estimation result, but easily being measured bad data, estimated result affects, need special bad data detection and identification program (as based on maximum weighted discrepancy principle) to reduce the impact of bad data on state estimation result, in the time there is extensive bad data in measurement, iterations is many, and may not necessarily effectively pick out all bad datas.The estimators such as QC, QL, WLAV, LMS, LTS all belong to robust and estimate category, its advantage is that state estimation algorithm is without extra bad data detection and identification program, but also exist need subjective to set that weighting factor, calculated amount are large, the defect such as robustness weakness under some particular case.

Domestic and international most state estimation algorithms are all that hypothesis measures variances sigma at present ²known, and in Practical Project, in large-scale electrical power system, there are numerous measuring instruments, measuring instrument is checked again and must pay very large economic cost, thereby be difficult to systematically again check measuring instrument; In addition along with the moment aging, running environment of equipment changes, the measurement accuracy of instrument is difficult to continue to keep stable.That is to say, it is time dependent measuring variance, is difficult to accurately estimate, measures weight and exists uncertain.

Summary of the invention

Goal of the invention: the present invention proposes a kind of electric system fuzzy self-adaption robust method of estimation, effectively identification bad data, avoided residual contamination and residual error to flood, have robust estimated performance capitally.

Technical scheme: the technical solution used in the present invention is a kind of electric system fuzzy self-adaption robust method of estimation, comprises the following steps:

1) fuzzy membership function of i measuring point poor quality of definition is:

$v_i\left(\right|r_i|,{\mathrm{\σ}}_i)=1-1/(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\σ}}_i}\right)}^b)\∀i=\mathrm{1,2},\mathrm{K\; m}$

In formula: r _ifor the residual error of measuring point i, σ _ifor the measuring standard of measuring point i poor, v _i(| r _i|, σ _i) be the fuzzy membership of measuring point i, a, b is greater than 0 constant;

2) taking the Weighted Fuzzy degree of membership sum that minimizes measuring point poor quality as optimization aim, following Optimized model is proposed:

$\left\{\begin{array}{c}\min J=\underset{i=1}{\overset{i=\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.r=z-h\left(x\right)\end{array}\right.$

In formula,

for the weight of measuring point i;

3) consider to measure weight uncertainty, i.e. the poor uncertainty of measuring standard, poor based on measuring point fuzzy membership correction measuring standard, to realize the self-adaptation to measuring rough error, for the k+1 time iteration, order:

${{\mathrm{\&sigma;}}_i}^{(k+1)}={{\mathrm{\&sigma;}}_i}^{\left(k\right)}\&CenterDot;f\left(v_i\right)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

4) introduce non-negative relaxation factor l, u, step 2) in Optimized model can be equivalent to:

$\left\{\begin{array}{c}\min J(l,u)=\underset{i=1}{\overset{i=\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.z-h\left(x\right)+l-u=0(l,u)>0\end{array}\right.$

5) equality constraint in step 3) is made as to barrier function, can be able to lower Lagrangian function:

$L(l,u,x,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;})=\underset{i=1}{\overset{i=\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">i</figure-callout>}}}^2}/(1-(1+{\left(\frac{|l_i+u_i|}{a{\mathrm{\&sigma;}}_i}\right)}^b))+\mathrm{\&lambda;}(z-h\left(x\right)+l-u)+{\mathrm{\&alpha;}}^{T}l+{\mathrm{\&beta;}}^{T}u$

In formula: λ, α, β are m dimension Lagrange multiplier, i.e. dual variable; L, u is former variable;

6) solution procedure 4) in the KKT condition of Lagrangian function obtain:

$\left\{\begin{array}{c}{L}_l=\frac{\&PartialD;J(l,u)}{\&PartialD;l}+\mathrm{\&lambda;}+\mathrm{\&alpha;}=0\\ L_u=\frac{\&PartialD;J(l,u)}{\&PartialD;l}-\mathrm{\&lambda;}+\mathrm{\&beta;}=0\\ L_x=\&dtri;h\left(x\right)\mathrm{\&lambda;}=0\\ L_{\mathrm{\&alpha;}}=l\&DoubleRightArrow;{L_{\mathrm{\&alpha;}}}^u=\mathrm{ALe}-\mathrm{\&mu;e}=0\\ L_{\mathrm{\&beta;}}=u\&DoubleRightArrow;{L_{\mathrm{\&beta;}}}^u=\mathrm{BUe}-\mathrm{\&mu;e}=0\\ L_{\mathrm{\&lambda;}}=z-h\left(x\right)+l-u=0\end{array}\right.$

In above formula, ▽ h (x) is the Jacobian matrix of h (x), L=diag (l ₁, l ₂..., l _m), U=diag (u ₁, u ₂, u _m); A=diag (α ₁, α ₂..., α _m), B=diag (β ₁, β ₂..., β _m), e=[1 ... 1] ^t, μ is the disturbance factor, and meets

8) will in step 5), for KKT condition, after newton-La Fusenfa linearization, obtain following update equation:

$\left[\begin{array}{cccccc}\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">l</figure-callout>}}^2}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}& I& \mathrm{<figure-callout\; id="0"\; label="I"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">I</figure-callout>}& 0& 0\\ \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u\&PartialD;l}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u^2}& -I& 0& I& 0\\ I& -\mathrm{<figure-callout\; id="0"\; label="I"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">I</figure-callout>}& 0& 0& 0& -{\&dtri;}^{T}h\left(x\right)\\ \mathrm{<figure-callout\; id="0"\; label="A"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">A</figure-callout>}& 0& 0& \mathrm{<figure-callout\; id="0"\; label="L"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">L</figure-callout>}& 0& 0\\ 0& \mathrm{<figure-callout\; id="0"\; label="B"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">B</figure-callout>}& 0& 0& \mathrm{<figure-callout\; id="0"\; label="U"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">U</figure-callout>}& 0\\ 0& 0& \&dtri;h\left(x\right)& 0& 0& {\&dtri;}^{2}h\left(x\right)\mathrm{\&lambda;}\end{array}\right]\&CenterDot;\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\\ \mathrm{d\&lambda;}\end{array}\right]=-\left[\begin{array}{c}{L}_l\\ L_u\\ L_x\\ L_{\mathrm{\&alpha;}}\\ L_{\mathrm{\&beta;}}\\ L_{\mathrm{\&lambda;}}\end{array}\right]$

In formula: ▽ ²h (x) is the gloomy matrix in sea of h (x);

8), according to the iterative step of former-dual interior point, upgrade quantity of state, until convergence.

Beneficial effect: the present invention proposes one and takes into account measurement weight probabilistic electric system fuzzy self-adaption robust algorithm for estimating (FARE), this algorithm has proposed the fuzzy membership concept of measuring point poor quality, and revise online the standard deviation of measuring point according to the degree of membership of measuring point, realize the self-adaptation to rough error.Solve with modern interior-point method, convergence is good, and result is not subject to initial value affecting, is applicable to solve Large-scale Optimization Problems.Effectively identification bad data, avoided residual contamination and residual error to flood, quantity of state and measurement amount are had to robust estimated performance capitally, there is engineering application background well.

Brief description of the drawings

Fig. 1 is workflow diagram of the present invention;

Fig. 2 is taking IEEE14 as standard testing node, the test design sketch of the inventive method FARE and WLAV, QC, WLS algorithm measurement amount robustness, and the weighted residual average of measurement | r/ σ | value distribution plan;

Fig. 3 is taking IEEE118 as standard testing node, the test design sketch of the inventive method FARE and WLAV, QC, WLS algorithm measurement amount robustness, and the weighted residual average of measurement | r/ σ | value distribution plan.

Embodiment

Below in conjunction with the drawings and specific embodiments, further illustrate the present invention, should understand these embodiment is only not used in and limits the scope of the invention for the present invention is described, after having read the present invention, those skilled in the art all fall within the application's claims limited range to the amendment of various equivalents of the present invention.

1. set up electric system fuzzy self-adaption robust estimation model

Power system state estimation measurement equation is:

z＝h(x)+e

In formula: x is state vector, the system nodes of setting up departments is n, amplitude and phase angle that x comprises node voltage, and its dimension is 2n-1; Z is that m dimension measures vector; H (x) is measurement function vectors; E is error in measurement vector.

The residual equation of state estimation is:

$r=z-h\left(\hat{x}\right)$

In formula: r is m dimension residual vector,

for the estimated value of quantity of state x.

For i measuring point, if it estimates residual error r _ivery little, can think that this measuring point is high-quality measuring point; Otherwise, if it estimates residual error r _ivery large, can think that this measuring point is measuring point inferior.Because the measurement accuracy of each measuring point is not quite similar, definition σ _ifor the measuring standard of measuring point i poor, thereby with the weighted residual r of measuring point _i/ σ _ithe good and bad measuring point of size discrimination more reasonable.

The quality that definition event y is measuring point, if set Y has comprised all event y, element y and the available fundamental function of relation of gathering Y---membership function v (y) represents, for classical data acquisition theory, has so:

But the quality of actual measuring point is relative concept, there are not absolute superiority and inferiority, than classical data acquisition theory, fuzzy set allows degree of membership to get [0,1] any value on, choose continuously differentiable bell membership function herein, for i measuring point, the fuzzy membership function of its poor quality is:

$v_i\left(\right|r_i|,{\mathrm{\&sigma;}}_i)=1-1/(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

In formula: a, b are greater than 0 fuzzy membership characteristic parameter.

Taking weighted least-squares (WLS) as main Length Factor Method in Power System State is taking minimum weight residual sum of squares (RSS) as optimization aim, similarly, fuzzy membership based on measuring point poor quality herein, taking the Weighted Fuzzy degree of membership sum that minimizes measuring point poor quality as optimization aim, following Optimized model is proposed:

$\left\{\begin{array}{c}\min J=\underset{i=1}{\overset{i=\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.r=z-h\left(x\right)\end{array}\right.$

In formula:

for the weight of measuring point i.

Domestic and international most state estimation algorithms are all that hypothesis measures variances sigma at present ²known, and in Practical Project, in large-scale electrical power system, there are numerous measuring instruments, measuring instrument is checked again and must pay very large economic cost, thereby be difficult to systematically again check measuring instrument; In addition due to the increase of service time, the external environment condition that the moment changes, the measurement accuracy of instrument is difficult to continue to keep stable.That is to say, measure variance and be difficult to accurately estimate, measure weight and exist uncertain.The general subjectivity of state estimation of domestic real system arranges the weight factor of each measurement, debugging and safeguard complicatedly, and the subjective weight factor of setting may not measure variance with reality and coincide.

For this reason, poor based on measuring point fuzzy membership correction measuring standard herein, for the k+1 time iteration in solving-optimizing object procedure, order:

${{\mathrm{\&sigma;}}_i}^{(k+1)}={{\mathrm{\&sigma;}}_i}^{\left(k\right)}\&CenterDot;f\left(v_i\right)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

The poor correction function f of measuring standard (v) should follow 2 principles: 1. as measuring point fuzzy membership v _iapproach at 0 o'clock, now this measuring point is high-quality measuring point, should keep σ _iapproximate constant; 2. as measuring point fuzzy membership v _iapproach at 1 o'clock, now this measuring point is measuring point inferior, should increase σ _i, while occurring measuring rough error, should reduce the weight of this measuring point in iteration, reduce its impact on state estimation result.

Based on above-mentioned criterion, can choose correction function f (v)=[1-(1/v)] ^1/b.

2. solving of the FARE model based on former dual interior point

With reference to Fig. 1, first the present invention obtains network parameter, topological structure, the measurement parameters of electric system, network parameter comprises: bus numbering, the reactance of node shunt capacitance, branch road head end numbering and end numbering, the reactance of branch road resistance, reactance, charging capacitor over the ground, transformer resistance, reactance, standard no-load voltage ratio; Topological structure mainly comprises the on off state connecting between two electrical equipments; Measurement parameters comprises: node injects meritorious, idle measurement, and branch road head end and end are meritorious, idle measurement, and busbar voltage amplitude measures.

Obtaining after the data of state estimation program needs, carry out program initialization, comprising: quantity of state voltage magnitude and phase angle, Lagrange multiplier, former dual relaxation variable, former antithesis Center Parameter, maximum iteration time and convergence precision.And form bus admittance matrix.

Introduce the non-negative relaxation factor l of m dimension, u, optimization aim can be equivalent to:

$\left\{\begin{array}{c}\min J(l,u)=\underset{i=1}{\overset{i=\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.z-h\left(x\right)+l-u=0(l,u)>0\end{array}\right.$

Equality constraint in above formula is made as to barrier function, can be able to lower Lagrangian function:

$L(l,u,x,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;})=\underset{i=1}{\overset{i=\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">m</figure-callout>}}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{|l_i+u_i|}{a{\mathrm{\&sigma;}}_i}\right)}^b))+\mathrm{\&lambda;}(z-h\left(x\right)+l-u)+{\mathrm{\&alpha;}}^{T}l+{\mathrm{\&beta;}}^{T}u$

In formula: λ, α, β are m dimension Lagrange multiplier, i.e. dual variable; L, u is former variable.

Solving its KKT condition can obtain:

$L_l=\frac{\&PartialD;J(l,u)}{\&PartialD;l}+\mathrm{\&lambda;}+\mathrm{\&alpha;}=0$

$L_u=\frac{\&PartialD;J(l,u)}{\&PartialD;l}-\mathrm{\&lambda;}+\mathrm{\&beta;}=0$

L _x＝▽h(x)λ＝0

$L_{\mathrm{\&alpha;}}=l\&DoubleRightArrow;{L_{\mathrm{\&alpha;}}}^u=\mathrm{ALe}-\mathrm{\&mu;e}=0$

$L_{\mathrm{\&beta;}}=u\&DoubleRightArrow;{L_{\mathrm{\&beta;}}}^u=\mathrm{BUe}-\mathrm{\&mu;e}=0$

L _λ＝z-h(x)+l-u＝0

In above formula, ▽ h (x) is the Jacobian matrix of h (x), L=diag (l ₁, l ₂..., l _m), U=diag (u ₁, u ₂, u _m); A=diag (α ₁, α ₂..., α _m), B=diag (β ₁, β ₂..., β _m), e=[1 ... 1] ^t, μ is the disturbance factor, and meets

KKT condition is Nonlinear System of Equations, and available newton-La Fusenfa solves, and will after its linearization, obtain update equation group:

$\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">l</figure-callout>}}^2}\mathrm{dl}+\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}\mathrm{du}+\mathrm{d\&lambda;}+\mathrm{d\&alpha;}=-L_l$

$\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}\mathrm{dl}+\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">u</figure-callout>}}^2}\mathrm{du}-\mathrm{d\&lambda;}+\mathrm{d\&beta;}=-L_u$

▽ ²h(x)λdx+▽h(x)dλ＝-L _x

-▽ ^Th(x)dx+dl-du＝-L _λ

$\mathrm{Ld\&alpha;}+\mathrm{Adl}=-L_{\mathrm{\&alpha;}}^u$

$\mathrm{Ud\&beta;}+\mathrm{Bdu}=-L_{\mathrm{\&beta;}}^u$

In formula: ▽ ²h (x) is the gloomy matrix in sea of h (x).Being write as matrix form can obtain:

$\left[\begin{array}{cccccc}\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">l</figure-callout>}}^2}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}& I& I& 0& 0\\ \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u\&PartialD;l}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u^2}& -I& 0& I& 0\\ I& -\mathrm{<figure-callout\; id="0"\; label="I"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">I</figure-callout>}& 0& 0& 0& -{\&dtri;}^{T}h\left(x\right)\\ A& 0& 0& \mathrm{<figure-callout\; id="0"\; label="L"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">L</figure-callout>}& 0& 0\\ 0& \mathrm{<figure-callout\; id="0"\; label="B"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">B</figure-callout>}& 0& 0& \mathrm{<figure-callout\; id="0"\; label="U"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">U</figure-callout>}& 0\\ 0& 0& \&dtri;h\left(x\right)& 0& 0& {\&dtri;}^{2}h\left(x\right)\mathrm{\&lambda;}\end{array}\right]\&CenterDot;\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\\ \mathrm{d\&lambda;}\end{array}\right]=-\left[\begin{array}{c}{L}_l\\ L_u\\ L_x\\ L_{\mathrm{\&alpha;}}\\ L_{\mathrm{\&beta;}}\\ L_{\mathrm{\&lambda;}}\end{array}\right]$

In formula: $\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}=\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u\&PartialD;l},\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">u</figure-callout>}}^2}=\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">l</figure-callout>}}^2},$ And can calculate [dl du dx] by above formula ^t, [d α d β d λ] ^t.

Calculate duality gap Gap=α ^tl+ β ^tu and disturbance factor mu=δ Gap/2m.In addition, permanent in non-negative in order to ensure relaxation factor, the iteration step length (θ of, dual variable former according to following formula _r, θ _d):

${\mathrm{\&theta;}}_r=0.9995\min \{\underset{i}{\min }(-\frac{l_i}{{\mathrm{dl}}_i}:{\mathrm{dl}}_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\}$

Revise former, dual variable:

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

Correction measuring standard is poor: σ ^(k+1)=σ ^(k)(1/v) ^1/b.

Repeat above-mentioned makeover process until duality gap Gap reaches the maximum iteration time that convergence precision or iterations exceed setting.

3. sample calculation analysis

Test example of the present invention comprises IEEE14, IEEE30, IEEE57, IEEE118 node, Poland 2383,2746 nodes are (with WP-2383, WP-2746 represents), metric data adds random error in measurement by strict trend result and obtains, and considers the poor uncertainty of measuring standard, the poor obedience [0.002 of voltage measuring standard, being uniformly distributed 0.005], power measurement standard deviation is obeyed being uniformly distributed of [0.004,0.01].Measure redundance

between 3.5～4.5.It is poor that bad data is added and subtracted [5,20] maximum measuring standard doubly on the basis of metric data, and add at random 10% bad data.

The superiority of FARE algorithm proposing for checking the present invention, by FARE with proposed in the past state estimation algorithm, comprise WLS, WLAV, QC, carry out comparative study from aspects such as bad data identification, robustness, constringency performances.The voltage of WLS, WLAV, QC algorithm measures weight and gets 4, and meritorious, idle measurement weight gets 1.For the poor initial value σ of measuring standard of FARE algorithm ⁽⁰⁾, voltage measures 0.001, and meritorious, idle measurement gets 0.002; Fuzzy membership parameter is got a=2.5, b=3.

3.1 measurement amount robustness

For the bad data identification of more various algorithms and resist the performance that residual contamination floods, measuring true value

the poor true value of measuring standard

in known situation, definition:

${\mathrm{\&tau;}}_i=|\widehat{z_i}-{\stackrel{\&OverBar;}{z}}_i|/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_i$

In formula:

be the estimated value of i measuring point, index τ has weighed measurement estimated value and has departed from the degree that measures true value.

Taking IEEE14, IEEE118 as test example, under 4 kinds of state estimation algorithms, the percent profile of τ is shown in Figure of description 2, Fig. 3.

From Fig. 2, Fig. 3, the number percent of τ > 3 under the estimation of WLS algorithm is more than 20%, be only 10% and measure bad data ratio, this explanation has occurred that very serious residual contamination floods phenomenon, thereby QC algorithm based on WLS is also difficult to effective identification bad data and resists residual contamination flood.Than WLAV, QC, WLS algorithm, FARE algorithm estimates that the distribution of lower τ mainly concentrates between 0～1, and the ratio of τ > 3 is obviously less than all the other 3 kinds of algorithms, thereby measurement estimated value under FARE algorithm is closer to measuring true value, has bad data identification well and resists residual contamination and flood performance.

3.2 quantity of state robustness

In order to compare the quantity of state robustness of 4 kinds of algorithms, adopt following 2 indexs to evaluate the estimated performance of different conditions algorithm for estimating to quantity of state:

1) square error of voltage magnitude:

$\mathrm{<figure-callout\; id="1"\; label="MSE"\; filenames="HDA0000476422930000011.png,HDA0000476422930000021.png"\; state="\{\{state\}\}">MSE</figure-callout>}1=\underset{i=1}{\overset{i=n}{\mathrm{\&Sigma;}}}{(V_i-\stackrel{}{\widehat{V_i}})}^2/n$

In formula: V _ifor the true value of node i voltage magnitude, for the estimated value of node i voltage magnitude.

2) square error of voltage phase angle:

$\mathrm{<figure-callout\; id="2"\; label="MSE"\; filenames="HDA0000476422930000021.png"\; state="\{\{state\}\}">MSE</figure-callout>}2=\underset{i=2}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&theta;}}_i-{\widehat{\mathrm{\&theta;}}}_i)}^2/(n-1)$

In formula: θ _ifor the true value of node i voltage phase angle,

for the estimated value of node i phase angle, acquiescence node 1 is balance node.

At 6 service systems, all containing 10% bad data in the situation that, table 1 has provided the MSE1 of 4 kinds of state estimation algorithms, the index of MSE2.

The size of MSE1 and MSE2 in comparison sheet 1, the estimated performance of known each algorithm to quantity of state: FARE>WLAV>QC>W LS, and FARE is obviously better than all the other 3 kinds of algorithms.In addition, FARE algorithm is better than the estimated performance to voltage magnitude to the estimated performance of voltage phase angle.

The robustness of comprehensive measurement amount and quantity of state is known, and in 4 kinds of state estimation algorithms choosing, the FARE that the present invention proposes is top quality Length Factor Method in Power System State.

Table 1: containing MSE1 and the MSE2 index performance of different algorithm for estimating under 10% bad data

Claims (1)

1. an electric system fuzzy self-adaption robust method of estimation, is characterized in that, comprises the following steps:

1) fuzzy membership function of i measuring point poor quality of definition is:

$v_i\left(\right|r_i|,{\mathrm{\&sigma;}}_i)=1-1/(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

In formula: r _ifor the residual error of measuring point i, σ _ifor the measuring standard of measuring point i poor, v _i(| r _i|, σ _i) be the fuzzy membership of measuring point i, a, b is greater than 0 constant;

2) taking the Weighted Fuzzy degree of membership sum that minimizes measuring point poor quality as optimization aim, following Optimized model is proposed:

$\left\{\begin{array}{c}\min J=\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.r=z-h\left(x\right)\end{array}\right.$

In formula,

for the weight of measuring point i;

3) consider to measure weight uncertainty, i.e. the poor uncertainty of measuring standard, poor based on measuring point fuzzy membership correction measuring standard, to realize the self-adaptation to measuring rough error, for the k+1 time iteration, order:

${{\mathrm{\&sigma;}}_i}^{(k+1)}={{\mathrm{\&sigma;}}_i}^{\left(k\right)}\&CenterDot;f\left(v_i\right)\&ForAll;i=\mathrm{1,2},\mathrm{K\; m}$

4) introduce non-negative relaxation factor l, u, step 2) in Optimized model can be equivalent to:

$\left\{\begin{array}{c}\min J(l,u)=\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{\left|r_i\right|}{a{\mathrm{\&sigma;}}_i}\right)}^b))\\ s.t.z-h\left(x\right)+l-u=0(l,u)>0\end{array}\right.$

5) equality constraint in step 3) is made as to barrier function, can be able to lower Lagrangian function:

$L(l,u,x,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;})=\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\frac{1}{{{\mathrm{\&sigma;}}_i}^2}/(1-(1+{\left(\frac{|l_i+u_i|}{a{\mathrm{\&sigma;}}_i}\right)}^b))+\mathrm{\&lambda;}(z-h\left(x\right)+l-u)+{\mathrm{\&alpha;}}^{T}l+{\mathrm{\&beta;}}^{T}u$

In formula: λ, α, β are m dimension Lagrange multiplier, i.e. dual variable; L, u is former variable;

6) solution procedure 4) in the KKT condition of Lagrangian function obtain:

$\left\{\begin{array}{c}{L}_l=\frac{\&PartialD;J(l,u)}{\&PartialD;l}+\mathrm{\&lambda;}+\mathrm{\&alpha;}=0\\ L_u=\frac{\&PartialD;J(l,u)}{\&PartialD;l}-\mathrm{\&lambda;}+\mathrm{\&beta;}=0\\ L_x=\&dtri;h\left(x\right)\mathrm{\&lambda;}=0\\ L_{\mathrm{\&alpha;}}=l\&DoubleRightArrow;{L_{\mathrm{\&alpha;}}}^u=\mathrm{ALe}-\mathrm{\&mu;e}=0\\ L_{\mathrm{\&beta;}}=u\&DoubleRightArrow;{L_{\mathrm{\&beta;}}}^u=\mathrm{BUe}-\mathrm{\&mu;e}=0\\ L_{\mathrm{\&lambda;}}=z-h\left(x\right)+l-u=0\end{array}\right.$

In above formula, ▽ h (x) is the Jacobian matrix of h (x), L=diag (l ₁, l ₂..., l _m), U=diag (u ₁, u ₂, u _m); A=diag (α ₁, α ₂..., α _m), B=diag (β ₁, β ₂..., β _m), e=[1 ... 1] ^t, μ is the disturbance factor, and meets

7) will in step 5), for KKT condition, after newton-La Fusenfa linearization, obtain following update equation:

$\left[\begin{array}{cccccc}\frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l^2}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;l\&PartialD;u}& I& I& 0& 0\\ \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u\&PartialD;l}& \frac{{\&PartialD;}^{2}J(l,u)}{\&PartialD;u^2}& -I& 0& I& 0\\ I& -I& 0& 0& 0& -{\&dtri;}^{T}h\left(x\right)\\ A& 0& 0& L& 0& 0\\ 0& B& 0& 0& U& 0\\ 0& 0& \&dtri;h\left(x\right)& 0& 0& {\&dtri;}^{2}h\left(x\right)\mathrm{\&lambda;}\end{array}\right]\&CenterDot;\left[\begin{array}{c}\mathrm{dl}\\ \mathrm{du}\\ \mathrm{dx}\\ \mathrm{d\&alpha;}\\ \mathrm{d\&beta;}\\ \mathrm{d\&lambda;}\end{array}\right]=-\left[\begin{array}{c}{L}_l\\ L_u\\ L_x\\ L_{\mathrm{\&alpha;}}\\ L_{\mathrm{\&beta;}}\\ L_{\mathrm{\&lambda;}}\end{array}\right]$

In formula: ▽ ²h (x) is the gloomy matrix in sea of h (x);

8), according to the iterative step of former-dual interior point, upgrade quantity of state, until convergence.

Images