尝试了半天，这一组参数是最好的

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

OPF.m

@@ -28,7 +28,7 @@ mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
 % mVolt=load('mVolt');
 % mVolt=mVolt.mVolt;
 
-mVolt(2)=rVolt(2)*(1-sigma*4);
+mVolt(2)=rVolt(2)*(1+0.2);
 %% 画Case A的图
 % figure('Color',[1 1 1]);
 [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有

subOPF.m

@@ -90,64 +90,14 @@ lVolt(noMeasurei)=0.7*mVolt(noMeasurei);%
 uVolt(noMeasurei)=0.7*mVolt(noMeasurei);
 %错误数据
 %mVolt(2)=5;
-bigM=2;
+bigM=0.3;
-Vbi=sparse(0.9*ones(Busnum,1));
+Vbi=sparse(0.0*ones(Busnum,1));
-PDbi=sparse(0.9*ones(length(Loadi),1));
+% Vbi(2)=1;
-QDbi=sparse(0.9*ones(length(Loadi),1));
+PDbi=sparse(0.0*ones(length(Loadi),1));
-eps=10;
+QDbi=sparse(0.0*ones(length(Loadi),1));
+eps=1;
 % 第一遍，算连续的值
 fprintf('第1次迭代，算连续量。\n');
-while(abs(Gap)>Precision*1)
-    if KK>kmax
-        break;
-    end
-    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;
-    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);
-    %%取各分量
-    [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
-% 第二遍，算离散
-% [~,~,Init_Z,Init_W,Init_L,Init_U,~,~,~,RestraintCount,~,~,~,~,~,~,~,~,~]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
-% Gap=(Init_L*Init_Z'-Init_U*Init_W');
-Gap=1000;
-% KK=0;
-eps=1;
-fprintf('\n');
-fprintf('第2次迭代，算离散量。\n');
-while eps>0.001
 while(abs(Gap)>Precision*10)
     if KK>kmax
         break;
@@ -174,9 +124,67 @@ while eps>0.001
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
+    bigM=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;
+    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);
+    %%取各分量
+    [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
+% 第二遍，算离散
+% [~,~,Init_Z,Init_W,Init_L,Init_U,~,~,~,RestraintCount,~,~,~,~,~,~,~,~,~]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
+% Gap=(Init_L*Init_Z'-Init_U*Init_W');
+Gap=1;
+% KK=0;
+eps=1;
+fprintf('\n');
+fprintf('第2次迭代，算离散量。\n');
+% while eps>0.0001
+    while(abs(Gap)>Precision*10 || eps>0.0001)
+        if KK>kmax
+            break;
+        end
+        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);
+        bigM=0.2;
+        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 eps>0.01
+        eps=Gap*2
+        else
+            eps=eps*0.2
+        end
         %     if KK>20
         %         eps=eps*0.1;
         %         eps=Gap;
@@ -184,9 +192,9 @@ while eps>0.001
         %             eps=1e-5;
         %         end
         %         eps;
-        % %        if any(Vbi>0.002)
+%                if any(Vbi>0.5)
-        % %            Vbi(Vbi>0.002)=1;
+%                    Vbi(Vbi>0.5)=1;
-        % %        end
+%                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);
@@ -204,9 +212,10 @@ while eps>0.001
         Gap=(Init_L*Init_Z'-Init_U*Init_W');
         KK=KK+1;
     end
-    eps=eps*0.9
+%     eps=eps*0.95;
-    Gap=100;
+%     eps=Gap
-end
+%     Gap=100;
+% end
 %% 计算最大不平衡量
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';