互补理论调得不好。

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

FormLxComa.m

@@ -4,8 +4,9 @@ t2=Lul+diag(Init_Z)*Lz;
 t3=inv(diag(Init_L));
 t4=t3*t2;%
 t5=Luu-diag(Init_W)*Lw;
-t6=inv(diag(Init_U));
+% t6=inv(diag(Init_U));
-t7=t6*t5;%
+% t7=t6*t5;%
+t7=diag(Init_U)\t5;
 t8=deltG*(t4+t7);%%
 LxComa=Lx+t8;
 end

OPF.m

@@ -5,7 +5,7 @@ arrayA=zeros(21,10);
 sumCaseA_SE=0;
 sumCaseB_SE=0;
 sumCaseC_SE=0;
-for I=1:20
+for I=1:1
 close
 [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]= ...
     pf('E:\算例\feeder33\feeder33.txt');
@@ -28,14 +28,18 @@ mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
 % mVolt=load('mVolt');
 % mVolt=mVolt.mVolt;
 
-% mVolt(19)=rVolt(19)*(1+sigma*4);
+mVolt(2)=rVolt(2)*(1-sigma*4);
 %% 画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);%全部有
+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;
+%真实值的
 plot(1:length(CaseAREV),abs(CaseAREV),'k.:','Marker','diamond');
+%测量值的
+% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond');
 box off;
 set(gca,'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])
 xlabel('节点号');
@@ -43,6 +47,7 @@ ylabel('
 subplot(4,1,2);
 CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
 % CaseAREA=CaseAREA(2:end)*100;
+%真实值的
 plot(1:length(CaseAREA),abs(CaseAREA),'k:','Marker','diamond');
 box off;
 set(gca,'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])
@@ -51,7 +56,10 @@ ylabel('
 subplot(4,1,3);
 CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
 % CaseAREP=CaseAREP(2:end)*100;
+%真实值的
 plot(1:length(CaseAREP),abs(CaseAREP),'k:','Marker','diamond');
+%测量值的
+% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond');
 box off;
 set(gca,'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])
 xlabel('节点号');
@@ -59,7 +67,10 @@ ylabel('
 subplot(4,1,4);
 CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
 % CaseAREQ=CaseAREQ(2:end)*100;
+%真实值的
 plot(1:length(CaseAREQ),abs(CaseAREQ),'k:','Marker','diamond');
+%测量值的
+% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond');
 box off;
 set(gca,'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])
 xlabel('节点号');
@@ -73,7 +84,6 @@ CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5;
 sumCaseA_SE=sumCaseA_SE+CaseA_SE;
 % arrayA(1:19,I)=Vbi;
 % arrayA(21,I)=CaseAE*1000;
-
 %% Case B
 % [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);%
 [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),PD0,QD0,mVolt,sigma);%
@@ -156,8 +166,6 @@ fprintf('Case A    Case B    Case C
 fprintf('%f    %f    %f    \n',CaseAE,CaseBE,CaseCE);
 fprintf('三个Case目标值\n')
 fprintf('%f\t%f\t%f \n',objA,objB,objC)
-
-end
 %% 画测量值
 % subplot(4,1,1);
 % plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram')
@@ -309,3 +317,10 @@ fz=find(abs(plotGapA)==0);
 fz=fz(1);
 figure('Name','互补曲线')
 plot(1:fz-1,plotGapA(1:fz-1));
+figure('Name','最大不平衡量');
+%% 最大不平衡量
+% maxDismatchPQ = [0.3123e-10  0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10];
+% bar(maxDismatchPQ);
+%% 计算时间
+calTime=[87 61; 211 87; 108 59;];
+bar(calTime,'stacked');

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=100000;
 Precision=Precision/1;
 %% 加误差
 %找DG
@@ -90,10 +90,10 @@ lVolt(noMeasurei)=0.7*mVolt(noMeasurei);%
 uVolt(noMeasurei)=0.7*mVolt(noMeasurei);
 %错误数据
 %mVolt(2)=5;
-bigM=10;
+bigM=2;
-Vbi=sparse(0.5*ones(Busnum,1));
+Vbi=sparse(0.9*ones(Busnum,1));
-PDbi=sparse(0.5*ones(length(Loadi),1));
+PDbi=sparse(0.9*ones(length(Loadi),1));
-QDbi=sparse(0.5*ones(length(Loadi),1));
+QDbi=sparse(0.9*ones(length(Loadi),1));
 eps=10;
 % 第一遍，算连续的值
 fprintf('第1次迭代，算连续量。\n');
@@ -144,62 +144,68 @@ end
 % Gap=(Init_L*Init_Z'-Init_U*Init_W');
 Gap=1000;
 % KK=0;
+eps=1;
 fprintf('\n');
 fprintf('第2次迭代，算离散量。\n');
-while(abs(Gap)>Precision*1)
+while eps>0.001
-    if KK>kmax
+    while(abs(Gap)>Precision*10)
-        break;
+        if KK>kmax
-    end
+            break;
-    plotGap(KK+1)=Gap;
-    Init_u=Gap/2/RestraintCount*CenterA;
-    AngleIJMat=0;
-    %% 开始计算OPF
-    %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
-    %% 形成不等式约束的雅克比
-    deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
-    %%
-    L_1Z=diag(Init_Z./Init_L);
-    U_1W=diag(Init_W./Init_U);
-    %% 形成海森阵
-    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
-    %% 形成ddHy
-    ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
-    %%  开始构建ddg
-    ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
-    %%  开始构建deltF
-    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
-    %% 形成方程矩阵
-    Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
-    Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,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);
-    Ly=Mat_H;
-    if KK>0
-        eps=eps*0.1;
-        if abs(eps)<1e-6
-            eps=1e-5;
         end
-        eps;
+        plotGap(KK+1)=Gap;
-%        if any(Vbi>0.002)
+        Init_u=Gap/2/RestraintCount*CenterA;
-%            Vbi(Vbi>0.002)=1;
+        AngleIJMat=0;
-%        end
+        %% 开始计算OPF
+        %% 形成等式约束的雅克比
+        deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
+        %% 形成不等式约束的雅克比
+        deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
+        %%
+        L_1Z=diag(Init_Z./Init_L);
+        U_1W=diag(Init_W./Init_U);
+        %% 形成海森阵
+        deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
+        %% 形成ddHy
+        ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
+        %%  开始构建ddg
+        ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
+        %%  开始构建deltF
+        deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
+        %% 形成方程矩阵
+        Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
+        Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,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);
+        Ly=Mat_H;
+        %     if KK>20
+        %         eps=eps*0.1;
+        %         eps=Gap;
+        %         if abs(eps)<1e-6
+        %             eps=1e-5;
+        %         end
+        %         eps;
+        % %        if any(Vbi>0.002)
+        % %            Vbi(Vbi>0.002)=1;
+        % %        end
+        %     end
+        Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
+        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
+        Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
+        YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
+        %% 开始解方程
+        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;
+        [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');
+        KK=KK+1;
     end
-    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
+    eps=eps*0.9
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
+    Gap=100;
-    Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-    YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
-    %% 开始解方程
-    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;
-    [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');
-    KK=KK+1;
 end
 %% 计算最大不平衡量
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);