修改为只只有不存在实时数据的时候才加到不等式中，收敛速度快很多。

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

FormGstate1.m

@@ -1,4 +1,4 @@
-function Mat_G=FormGstate1(I1r,I1i)
+function Mat_G=FormGstate1(I1r,I1i,Loadi)
 %t1=PG(PVi);
 %GP=t1;%发电机P
 %GP=[4.5 4.5]';
@@ -10,7 +10,7 @@ function Mat_G=FormGstate1(I1r,I1i)
 % t4=sum(t3,2);%发电机Q
 %GQ=t4;
 Mat_G=[
-    I1r;
-    I1i;
+    I1r(Loadi);
+    I1i(Loadi);
     ];
 end

IPMLoop.m

@@ -8,12 +8,15 @@ I1r=real(guessIf1);
 I1i=imag(guessIf1);
 KK=0;
 plotGap=zeros(1,60);
-%³õʼ»¯
-%״̬Á¿Îª SEPD SEQD SEVmf1 SEVaf1
-RestraintCount=length(Loadi)*2;
+
 % ContrlCount=busNum*2+length(Loadi)*2;
 noBoundedLoadi=setdiff(Loadi,noLoadi);
 % noBoundedLoadi=[];
+boundedLoadi=setdiff(Loadi,noBoundedLoadi);
+%³õʼ»¯
+%״̬Á¿Îª SEPD SEQD SEVmf1 SEVaf1
+% RestraintCount=length(Loadi)*2;
+RestraintCount=length(boundedLoadi)*2;
 ContrlCount=busNum*2+length(Loadi)*2;
 CenterA=0.1;
 Init_Z=1*sparse(ones(RestraintCount,1));
@@ -22,7 +25,7 @@ Init_L=1*sparse(ones(RestraintCount,1));
 Init_U=1*sparse(ones(RestraintCount,1));
 Init_Y=sparse(2*busNum,1);%等式约束乘子
 Gap=(Init_L'*Init_Z-Init_U'*Init_W);
-kmax=100;
+kmax=20;
 isSetBound=0;
 largerBound=0;
 realBound=1;
@@ -37,7 +40,8 @@ while(abs(Gap)>1e-4)
     %% 形成等式约束的雅克比
     deltH=func_deltH(busNum,fsY1,Loadi,Balance);
     %% 形成不等式约束的雅克比
-     deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
+%      deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
+deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
     %%
 %     L_1Z=diag(Init_Z./Init_L);
 %     U_1W=diag(Init_W./Init_U);
@@ -55,22 +59,24 @@ while(abs(Gap)>1e-4)
     %%
     Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
     Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
-    Mat_G=FormGstate1(I1r,I1i);
+%     Mat_G=FormGstate1(I1r,I1i);
+    Mat_G=FormGstate1(I1r,I1i,ismember(Loadi,boundedLoadi) );
     Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
     Ly=Mat_H;
     
 %     if isSetBound==0
     
-        Ir=real(guessIf1);
+        Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
         pIr=find(Ir>0);
         nIr=find(Ir<0);
-        Ii=imag(guessIf1);
+        Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
         pIi=find(Ii>0);
         nIi=find(Ii<0);
-        lower=ones(length(Loadi)*2,1);
+%         lower=ones(length(Loadi)*2,1);
+lower=ones(length(boundedLoadi)*2,1);
 %     if abs(Gap)<1 || isSetBound==1
 
-%     if abs(Gap)<0.01
+%     if abs(Gap)<0.001
 %         largerBound=0;
 %         realBound=1;
 %     end
@@ -85,8 +91,8 @@ while(abs(Gap)>1e-4)
         lower(nIr)=(1.300)*Ir(nIr);
         lower(pIi+length(Ir))=(0.7)*Ii(pIi);
         lower(nIi+length(Ir))=(1.30)*Ii(nIi);
-        lower( ismember(Loadi,noBoundedLoadi))=-0.9;
-        lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-0.9;
+%         lower( ismember(Loadi,noBoundedLoadi))=-1.9;
+%         lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-1.9;
      end
         
     if largerBound==1
@@ -112,13 +118,13 @@ while(abs(Gap)>1e-4)
 %         I1r(I1r./lowerR>0.9998)=lowerR(I1r./lowerR>0.9998);
 % I1i(I1i./lowerI>0.9998)=lowerI(I1i./lowerI>0.9998);
     
-        Ir=real(guessIf1);
+        Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
         pIr=find(Ir>0);
         nIr=find(Ir<0);
-        Ii=imag(guessIf1);
+        Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
         pIi=find(Ii>0);
         nIi=find(Ii<0);
-        upper=ones(length(Loadi)*2,1);
+        upper=ones(length(boundedLoadi)*2,1);
 %     if abs(Gap)<1 || isSetBound==1
    if realBound==1
 %         upper(pIr)=(1.10+3/(KK+1))*Ir(pIr);
@@ -130,8 +136,8 @@ while(abs(Gap)>1e-4)
         upper(nIr)=(0.7)*Ir(nIr);
         upper(pIi+length(Ir))=(1.30)*Ii(pIi);
         upper(nIi+length(Ir))=(0.7)*Ii(nIi);
-        upper(ismember(Loadi,noBoundedLoadi))=0.9;
-        upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=0.9;
+%         upper(ismember(Loadi,noBoundedLoadi))=1.9;
+%         upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=1.9;
    end
 
     if largerBound==1
@@ -160,13 +166,14 @@ while(abs(Gap)>1e-4)
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum,Balance);
     [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref);
     Gap=(Init_L'*Init_Z-Init_U'*Init_W);+max(abs(deltX));
+%     Gap=max([max(abs(Mat_H)),max(Lz),max(Lw) ]);
     fprintf('Gap %f :%d\n',full(Gap),KK+1);
     KK=KK+1;
 end
 toc
 calcuTime=toc/KK;
 % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
-if abs(Gap)<0.00001 && KK<kmax
+if abs(Gap)<1e-4 && KK<kmax
     isConverged=1;
 else
     isConverged=0;

func_deltGstate1.m

@@ -1,7 +1,10 @@
-function deltG=func_deltGstate1(busNum,Loadi,I1r,I1i)
+function deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i)
+ind_=find(ismember(Loadi,boundedLoadi));
 deltG=[
-    zeros(busNum*2,length(Loadi)*2);
-    1*eye(length(Loadi)*1),zeros(length(Loadi)*1);%I1r
-    zeros(length(Loadi)*1),1*eye(length(Loadi)*1);%I1i
+    zeros(busNum*2,length(boundedLoadi)*2);
+    sparse(ind_,1:length(boundedLoadi),1,length(I1r),length(boundedLoadi)),zeros(length(I1i)*1,length(boundedLoadi));
+    zeros(length(I1r)*1,length(boundedLoadi)),sparse(ind_,1:length(boundedLoadi),1,length(I1i),length(boundedLoadi));
+%     1*eye(length(Loadi)*1),zeros(length(I1i)*1);%I1r
+%     zeros(length(Loadi)*1),1*eye(length(Loadi)*1);%I1i
 ];
 end

run.m

@@ -4,10 +4,10 @@ function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_m
 close all
 clc
 clear
-lineZ=readLineZ('feeder26\lineParameter.txt');
+lineZ=readLineZ('feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
-    cap]=dataRead(lineZ,'feeder26\data1.txt');
+    cap]=dataRead(lineZ,'feeder13\data1.txt');
 a=exp(1j*2*pi/3);
 Tp2f=1/3*[1 1 1;
     1 a a^2;
@@ -291,11 +291,11 @@ wV0i=wV0i(Loadi);
 
 %% 不同量测配置
 noLoadi=[1,8,10,12];
-noLoadi=[1:7:max(Loadi)];
-noLoadi=[1,19];
+% noLoadi=[1:3:max(Loadi)];
+% noLoadi=[1,19];
 % noLoadi=[4,3,2,5,6,7,11,14,22];
 % noLoadi=Loadi;
-% noLoadi=[];
+noLoadi=[1];
 wI0r(ismember(Loadi,noLoadi))=1e100000;
 wI0i(ismember(Loadi,noLoadi))=1e100000;
 wI1r(ismember(Loadi,noLoadi))=1e100000;

test.m

@@ -14,7 +14,7 @@ totalTime2Sum=10000;
 totalTime0Sum=10000;
 timeOfPrepare=0.00272333333;%准备数据的时间
 timeOfPrepare=0;
-N=500;
+N=200;
 loopN=1;
 while 1
     [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0]=run();