1.加了一个画错误数据的图

2.开始一点点做测试

Signed-off-by: dugg@lab-desk <dugg@lab-desk>

OPF.m

@@ -5,34 +5,54 @@ arrayA=zeros(21,10);
 sumCaseA_SE=0;
 sumCaseB_SE=0;
 sumCaseC_SE=0;
-for I=1:1
+badDataResult=zeros(10,33);
-close
+badDataLocation=zeros(34,10);
-[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,Branchi,Branchg,Branchb,Transfork0]= ...
+maxDVolt=0;
-    pf('E:\算例\feeder33\feeder33.txt');
+maxDVAngle=0;
-sigma=0.01;
+% for badDataNode=1:33
-RealPD=PD;
+    for I=1:5
-RealQD=QD;
+        close
-rVolt=Volt;
+        [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,Branchi,Branchg,Branchb,Transfork0]= ...
-Loadi=find(PD~=0);
+            pf('E:\算例\feeder33\feeder33.txt');
-PD0=sparse(Busnum,1);
+        sigma=0.01;
-QD0=sparse(Busnum,1);
+        RealPD=PD;
-PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+        RealQD=QD;
-QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+        rVolt=Volt;
-mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
+        Loadi=find(PD~=0);
+        PD0=sparse(Busnum,1);
+        QD0=sparse(Busnum,1);
+        PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+        QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+        mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
         
-%加载保存的变量
+        %加载保存的变量
-% PD0=load('PD0');
+        % PD0=load('PD0');
-% PD0=PD0.PD0;
+        % PD0=PD0.PD0;
-% QD0=load('QD0');
+        % QD0=load('QD0');
-% QD0=QD0.QD0;
+        % QD0=QD0.QD0;
-% mVolt=load('mVolt');
+        % mVolt=load('mVolt');
-% mVolt=mVolt.mVolt;
+        % mVolt=mVolt.mVolt;
         
-mVolt(2)=rVolt(2)*(1+sigma*8);
+        mVolt(18)=rVolt(18)*(1-sigma*6);
-%% 画Case A的图
+        %% 画Case A的图
-% figure('Color',[1 1 1]);
+        % figure('Color',[1 1 1]);
-[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
+        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
+%         badDataResult(I,badDataNode)=sum(Vbi);
+%         badDataLocation(1:33,I)=Vbi;
+%         badDataLocation(34,I)=sum(abs((rVolt-Volt)./rVolt./length(rVolt)))+sum(abs( (UAngel(2:33)-rUAngel(2:33))./rUAngel(2:33)./length(rUAngel(2:33))));
+
+% 加最大偏差
+maxDVolt_=max(abs((rVolt-Volt)));
+if maxDVolt_>maxDVolt
+    maxDVolt=maxDVolt_;
 end
+maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33))));
+if maxDVAngle_>maxDVAngle
+    maxDVAngle=maxDVAngle_;
+end
+    
+end
+% end
 subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]);
 CaseAREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case A
 % CaseAREV=CaseAREV(2:end)*100;

graph/badDataGraph.m

@@ -0,0 +1,16 @@
+x=1:33;
+y=1:33;
+x=[x,15,33];
+y=[y,17,23];
+scatter(x,y);
+xlabel('存在坏数据的节点');
+ylabel('辨识出坏数据的点');
+hold on;
+for I=1:33
+% plot([1,1,1,1]-1+I,[[1:3]-2+I,I]);
+plot([I,I],[0,I]);
+end
+% plot([1,1,1,1]-1+15,[1:3]-2+17);
+% plot([1,1,1,1]-1+33,[1,2,12]-2+23);
+plot([15,15],[0,17]);
+plot([33,33],[0,23]);

subOPF.m

@@ -45,7 +45,7 @@ Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,60);
 ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2;
-kmax=100;
+kmax=500;
 Precision=Precision/1;
 %% 加误差
 %找DG
@@ -124,7 +124,7 @@ while(abs(Gap)>Precision*10)
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    bigM=1;
+    bigM=0.6;
     Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
     Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
     Ly=Mat_H;
@@ -150,15 +150,9 @@ Gap=1;
 fprintf('\n');
 
 fprintf('第2次迭代，算离散量。\n');
-% while eps>0.000001
-% Init_Z=sparse(ones(1,RestraintCount));
-% Init_W=sparse(-1*ones(1,RestraintCount));
-% Init_L=1*sparse(ones(1,RestraintCount));
-% Init_U=1*sparse(ones(1,RestraintCount));
-% Init_Y=sparse(1,2*Busnum);%Óëѧ½ãÒ»ÖÂ
 while(abs(Gap)>Precision*10)
     if KK>kmax
-        %             break;
+                    break;
     end
     plotGap(KK+1)=Gap;
     Init_u=Gap/2/RestraintCount*CenterA;
@@ -182,7 +176,7 @@ while(abs(Gap)>Precision*10)
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    bigM=0.7;
+    bigM=0.4;
     eps=Gap*0.1;
     Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
     Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
@@ -192,24 +186,14 @@ while(abs(Gap)>Precision*10)
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
+%         Vbi(Vbi>0.003)=.1;
     fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
-    %%È¡¸÷·ÖÁ¿
-    %     Vbi_=Vbi;
-    %     PDbi_=PDbi;
-    %     QDbi_=QDbi;
-    %                 Vbi=1./(1+exp( (0.06-Vbi) /0.01))*0.1;
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
     [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi);
     Gap=(Init_L*Init_Z'-Init_U*Init_W');
-    
-    %         Vbi(Vbi>0.01)=1;
     KK=KK+1;
 end
-    
-%     eps=Gap
-%     Gap=100;
-% end
 %% 计算最大不平衡量
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';