加其它方法。

Signed-off-by: dmy@lab <dmy@lab.lab>

@Opti/equ.asv

@@ -38,27 +38,35 @@ this.cl=zeros(length(output_args),1);
 %% 开始增加不等式约束-电压\PD\QD
 %电压不等式约束
 output_args=[output_args;Volt];
-this.cu=[this.cu;1.07*ones(length(Volt),1)];%9½Úµã
+% this.cu=[this.cu;1.07*ones(length(Volt),1)];%9½Úµã
-this.cl=[this.cl;0.93*ones(length(Volt),1)];%9½Úµã
+% this.cl=[this.cl;0.93*ones(length(Volt),1)];%9½Úµã
-% this.cu=[this.cu;1.00*ones(length(Volt),1)];
+this.cu=[this.cu;1.00*ones(length(Volt),1)];
-% this.cl=[this.cl;0.90*ones(length(Volt),1)];
+this.cl=[this.cl;0.90*ones(length(Volt),1)];
 %% PD
 rPD=this.rPD;
 output_args=[output_args;PD];
 % PDU=[0.124;0.315;0.5;1;1;0.5;0.63;0.4];%原始数据
 % PDU=[0.63;0.4;0.5;1;0.8;1;0.63;0.4];%偏差比较大
-this.cu=[this.cu;1.2*rPD];
+%Generate values from the uniform distribution on the interval [a, b].
-this.cl=[this.cl;0.8*rPD];
+global loadFlag;
-this.cu=[this.cu;PDU];
+if loadFlag==0
-this.cl=[this.cl;0*PDU];
+    r =-0.201  + (-0.01-(-0.2)).*rand(length(rPD),1);
+else
+    
+end
+
+this.cu=[this.cu;(r+1.4).*rPD];
+this.cl=[this.cl;(r+1).*rPD];
+% this.cu=[this.cu;PDU];
+% this.cl=[this.cl;0*PDU];
 %% QD
 rQD=this.rQD;
 output_args=[output_args;QD];
-QDU=PDU;
+% QDU=PDU;
-% this.cu=[this.cu;1.2*rQD];
+this.cu=[this.cu;(r+1.4).*rQD];
-% this.cl=[this.cl;0.8*rQD];
+this.cl=[this.cl;(r+1).*rQD];
-this.cu=[this.cu;QDU];
+% this.cu=[this.cu;QDU];
-this.cl=[this.cl;0*QDU];
+% this.cl=[this.cl;0*QDU];
 %% 稠密化
 output_args=full(output_args);
 this.cu=full(this.cu);

@Opti/equ.m

@@ -47,16 +47,27 @@ rPD=this.rPD;
 output_args=[output_args;PD];
 % PDU=[0.124;0.315;0.5;1;1;0.5;0.63;0.4];%原始数据
 % PDU=[0.63;0.4;0.5;1;0.8;1;0.63;0.4];%偏差比较大
-this.cu=[this.cu;1.2*rPD];
+%Generate values from the uniform distribution on the interval [a, b].
-this.cl=[this.cl;0.8*rPD];
+global loadFlag;
+if loadFlag==1
+    r=load('rLD');
+    r=r.r;
+    
+else
+    %r=-0.2;
+    r =-0.201  + (-0.001-(-0.2)).*rand(length(rPD),1);
+end
+
+this.cu=[this.cu;(r+1.4).*rPD];
+this.cl=[this.cl;(r+1).*rPD];
 % this.cu=[this.cu;PDU];
 % this.cl=[this.cl;0*PDU];
 %% QD
 rQD=this.rQD;
 output_args=[output_args;QD];
 % QDU=PDU;
-this.cu=[this.cu;1.2*rQD];
+this.cu=[this.cu;(r+1.4).*rQD];
-this.cl=[this.cl;0.8*rQD];
+this.cl=[this.cl;(r+1).*rQD];
 % this.cu=[this.cu;QDU];
 % this.cl=[this.cl;0*QDU];
 %% 稠密化

LineP.m

@@ -0,0 +1,7 @@
+function [ output_args ] = LineP( Linei,Linej,Liner,Linex,SEVolt,SEVAngel )
+L=LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
+SEVolt=SEVolt.*exp(1j*SEVAngel);
+output_args=real(SEVolt(Linei).*conj(L));
+
+end
+

LineQ.m

@@ -0,0 +1,7 @@
+function [ output_args ] = LineQ( Linei,Linej,Liner,Linex,SEVolt,SEVAngel )
+L=LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
+SEVolt=SEVolt.*exp(1j*SEVAngel);
+output_args=imag(SEVolt(Linei).*conj(L));
+
+end
+

LoadCurve.m

@@ -1,10 +1,13 @@
+clear
+clc
+close all
 txt=dlmread('E:\书籍\论文\配电网状态估计\负荷.txt');
 txt=txt./10000;
 x=1:length(txt);
-plot(x,txt,'k');
+plot(x,txt,'k','linewidth',2);
 hold on
-plot(x,txt*0.8,'k--');
+plot(x,txt*0.8,'k--','linewidth',2);
-plot(x,txt*1.2,'k--');
+plot(x,txt*1.2,'k--','linewidth',2);
 legend('典型负荷曲线','考虑波动后的负荷范围')
 
 ylabel('负荷值/kW')

OPF_Init.asv

@@ -80,7 +80,7 @@ wVolt(notLoadi)=0;
 %% 50%
 % noLoad=[10 11 14 15 25 17 22  20];
 % noLoad=[2 5 6 9];
-%noLoad=[3,5,6,9,11,12,15,17,20,21,22,23,27,29,30,31];%33节点
+% noLoad=[3,5,6,9,11,12,15,17,20,21,22,23,27,29,30,31];%33节点
 % wPD(noLoad)=0;
 % wQD(noLoad)=0;
 % wVolt(noLoad)=0;

ReadMe.txt

@@ -1 +1,2 @@
+2015.3.18 依据电工学报专家的意见，添加其它状态估计方法做对比。
 2013.10.15 准备加入支路功率量测量试试，原来只有电压和负荷的效果太不好了，主要是负荷的估计效果不好。

Run_YALMIP2.asv

@@ -2,11 +2,14 @@ clc
 clear
 close all
 % yalmip('clear')
+global loadFlag;
+LineCurs=zeros(32,100);
 tic
-for I=1:1
+for I=1:100
+    close all;
     [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB, ...
         Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL, ...
-        Liner,Linex,Lineb,Transforr,Transforx,Transfork0]=pf('E:/ËãÀý/feeder33/feeder33.txt');
+        Liner,Linex,Lineb,Transforr,Transforx,Transfork0]=pf('E:\算例\feeder33\feeder33.txt');
     %pf('C:\bpa\长虹世纪线_556844_2013-10-28\长虹世纪线_556844_2013-10-28_iPso_newFil
     %e.txt');
      %   pf('E:/算例/东际911_2751267_2012-09-05/pyth-增加3倍.txt'); 20131125
@@ -57,6 +60,7 @@ for I=1:1
 %     load('QD0');
 %     load('mVolt');
 %     load('mLoadCurrent');
+    loadFlag=0;
     %% 17号负荷上加负数
     % PD0(17)=-10;
     %% 负荷用Pd=(1+k)Pd0添加负荷
@@ -65,7 +69,7 @@ for I=1:1
     % load('I');
     sigma=0.03;
     %
-%     [ mVolt,PD0,QD0,mLoadCurrent ] = GetMeasure(sigma,rVolt,rPD,rQD,rLoadCurrent );
+    [ mVolt,PD0,QD0,mLoadCurrent ] = GetMeasure(sigma,rVolt,rPD,rQD,rLoadCurrent );
     % sigma=sigma*(I-1)/50;
     % 排除大于3倍sigma的数据
     % save('mVolt','mVolt');
@@ -124,6 +128,11 @@ for I=1:1
     % sum([SEVolt;PD;QD]>cu(length(SEVolt)*2+2:end));
     % sum([SEVolt;PD;QD]<cl(length(SEVolt)*2+2:end));
 %     RealValuePlot(mPD(PDi),mQD(QDi),rPD(PDi),rQD(QDi),mVolt,rVolt,rVAngel);
+    %measurementN=sum(find(wVolt)>0)+sum(find(wPD)>0)+sum(find(wQD)>0);
+    stE=StErrorS(SEVolt,SEVAngel,PD,QD,rVolt,rVAngel,rPD(PDi),rQD(QDi));
+    fprintf('真实值统计误差为: %f',stE);
+    stE=StErrorZ(SEVolt,PD,QD,rVolt,mPD(PDi),mQD(QDi),noLoadi);
+    fprintf('测量值统计误差为: %f',stE);
     NormalizedResiduals(x,sigma,PDi,QDi,Volt0,mPD,mQD);
     SECurrent=LoadCurrent( SEVolt,SEVAngel,PD,QD,PDi,QDi );
     [flag,t1,t2]=MaxSigma( x,PDi,QDi,mPD,mQD,mVolt,mLoadCurrent,sigma,rPD,rQD,rVolt,rLoadCurrent );
@@ -132,11 +141,14 @@ for I=1:1
     [ok,msg] = checkSol(Opt);
     figure();
     plotError( SEVolt,SEVAngel,rVolt,rVAngel,mVolt,PD,mPD(PDi),rPD(PDi),QD,mQD(QDi),rQD(QDi),PDi );
-%     LineCurrent( Linei,Linej,Liner,Linex,Volt,VAngle )
+%     LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngel )
+LineCurs(:,I)=LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
+    LineP( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
+    LineQ( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
     if flag==1
         I
         t1
         t2
-        break;
+        %break;
     end
 end

Run_YALMIP2.m

@@ -2,8 +2,11 @@ clc
 clear
 close all
 % yalmip('clear')
+global loadFlag;
+LineCurs=zeros(32,100);
 tic
 for I=1:1
+    close all;
     [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB, ...
         Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL, ...
         Liner,Linex,Lineb,Transforr,Transforx,Transfork0]=pf('E:\算例\feeder33\feeder33.txt');
@@ -53,10 +56,11 @@ for I=1:1
     [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,wLoadCurrent,PD,PD0,QD,randPDind,Loadi,noLoadi]=OPF_Init(Busnum,Balance,PG,QG,rVolt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,rQD,rPD,rLoadCurrent);
     %% 开始
     %% 加噪音
-    load('PD0');
+%     load('PD0');
-    load('QD0');
+%     load('QD0');
-    load('mVolt');
+%     load('mVolt');
-    load('mLoadCurrent');
+%     load('mLoadCurrent');
+    loadFlag=0;
     %% 17号负荷上加负数
     % PD0(17)=-10;
     %% 负荷用Pd=(1+k)Pd0添加负荷
@@ -65,7 +69,7 @@ for I=1:1
     % load('I');
     sigma=0.03;
     %
-%     [ mVolt,PD0,QD0,mLoadCurrent ] = GetMeasure(sigma,rVolt,rPD,rQD,rLoadCurrent );
+    [ mVolt,PD0,QD0,mLoadCurrent ] = GetMeasure(sigma,rVolt,rPD,rQD,rLoadCurrent );
     % sigma=sigma*(I-1)/50;
     % 排除大于3倍sigma的数据
     % save('mVolt','mVolt');
@@ -89,7 +93,7 @@ for I=1:1
     opts.tolafun=1e-4;
     opts.warnings='all';
     opts.display='off';
-    x0=[0.1*rPD(PDi);0.1*rQD(QDi); ...
+    x0=[rand()*rPD(PDi);rand()*rQD(QDi); ...
         ones(length(Volt),1); ...
         zeros(length(Volt),1)];
     % x0=[PD(PDi);QD(QDi);xVolt';xUAngel'];
@@ -127,7 +131,7 @@ for I=1:1
     %measurementN=sum(find(wVolt)>0)+sum(find(wPD)>0)+sum(find(wQD)>0);
     stE=StErrorS(SEVolt,SEVAngel,PD,QD,rVolt,rVAngel,rPD(PDi),rQD(QDi));
     fprintf('真实值统计误差为: %f',stE);
-    stE=StErrorZ(SEVolt,PD,QD,rVolt,mPD(PDi),mQD(QDi),noLoadi);
+    stE=StErrorZ(SEVolt,PD,QD,mVolt,mPD(PDi),mQD(QDi),noLoadi);
     fprintf('测量值统计误差为: %f',stE);
     NormalizedResiduals(x,sigma,PDi,QDi,Volt0,mPD,mQD);
     SECurrent=LoadCurrent( SEVolt,SEVAngel,PD,QD,PDi,QDi );
@@ -137,11 +141,14 @@ for I=1:1
     [ok,msg] = checkSol(Opt);
     figure();
     plotError( SEVolt,SEVAngel,rVolt,rVAngel,mVolt,PD,mPD(PDi),rPD(PDi),QD,mQD(QDi),rQD(QDi),PDi );
-%     LineCurrent( Linei,Linej,Liner,Linex,Volt,VAngle )
+%     LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngel )
+% LineCurs(:,I)=LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
+    LineP( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
+    LineQ( Linei,Linej,Liner,Linex,SEVolt,SEVAngel );
     if flag==1
-        I
+        I;
-        t1
+        t1;
-        t2
+        t2;
-        break;
+        %break;
     end
 end

StErrorS.m

@@ -3,8 +3,10 @@ function [ output_args ] = StErrorS( SEVolt,SEVangle,SEPD,SEQD,rVolt,rVangle,rSE
 %   Detailed explanation goes here
 t1=[SEVolt;SEVangle;SEPD;SEQD];
 t2=[rVolt;rVangle;rSEPD;rSEQD];
+t1=t1(t2~=0);
+t2=t2(t2~=0);
 m=length(t1);
-t3=sum( ((t1-t2)./0.03).^2)/m;
+t3=sum( ((t1-t2)./0.999./t2).^2)/m;
 output_args=t3^.5;
 
 end

graph/9/201406163/ReadME.txt

@@ -0,0 +1,2 @@
+这个比较好
+2014.6.16 17：01