由于是直角坐标下写的，平衡节点不要置0置1，否则有问题。

现在好了，只把电压放目标函数中，电压误差就为0了。

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

AssignXX.m

@@ -1,5 +1,7 @@
-function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum)
+function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum,Balance)
 deltX=XX(1:ContrlCount);
+deltX(Balance)=0;
+deltX(Balance+Busnum)=0;
 k1=ContrlCount+2*Busnum;
 deltY=XX(ContrlCount+1:k1);
 k2=k1;

IPMLoop.m

@@ -29,6 +29,7 @@ while(abs(Gap)>0.000001)
     deltH=func_deltH(busNum,fsY1,Loadi,Balance);
     %% 形成不等式约束的雅克比
     deltG=func_deltG(busNum,Loadi,I1r,I1i);
+    deltG=zeros(size(deltG));
     %%
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
@@ -40,6 +41,7 @@ while(abs(Gap)>0.000001)
     ddh=0;
     %%  开始构建ddg
     ddg=func_ddg(busNum,Loadi,Init_Z,Init_W);
+    ddg=zeros(size(ddg));
     %%  开始构建deltF
     deltF=func_deltF(Vmeasurement,Imeasurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i);
 %     deltF=0;
@@ -47,6 +49,7 @@ while(abs(Gap)>0.000001)
     Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
     Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
     Mat_G=FormG(I1r,I1i);
+    Mat_G=zeros(size(Mat_G));
     Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
     Ly=Mat_H;
     Lz=FormLz(Loadi,Mat_G,Init_L);
@@ -56,13 +59,13 @@ while(abs(Gap)>0.000001)
     %% 开始解方程
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,busNum);
     %%取各分量
-    [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum);
+    [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);
+    Gap=(Init_L'*Init_Z-Init_U'*Init_W)+max(abs(deltX));
     fprintf('Gap %f\n',full(Gap));
     KK=KK+1;
 end
-f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V1r(Loadi)).^2)+sum((real(Vmeasurement)-V1i(Loadi)).^2);
+f=sum(([real(Imeasurement);imag(Imeasurement)]-[-I1r;-I1i]).^2)+sum((real(Vmeasurement)-V1r(Loadi)).^2)+sum((imag(Vmeasurement)-V1i(Loadi)).^2);
 
 end
 

Modification.m

@@ -12,9 +12,9 @@ Init_Y=Init_Y+AlphaD*deltY;
 %PG(PGi)=PG(PGi)+deltX(1:size(PGi,1));
 % PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
 V1r=V1r+AlphaP*deltX(1:Busnum);
-V1r(Balance)=Vref;
+% V1r(Balance)=Vref;
 V1i=V1i+AlphaP*deltX(Busnum+1:2*Busnum);
-V1i(Balance)=0;
+% V1i(Balance)=0;
 %QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 % QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 t=deltX(Busnum*2+1:ContrlCount);

SolveIt.m

@@ -11,15 +11,15 @@ yy=[LxComa;-Ly];
 
 % t=size(PVi,1)+size(PGi,1);%电压不变
 % yy(t+Balance)=0;
-aa(Balance,:)=0;
-aa(:,Balance)=0;
-aa=aa+sparse(Balance,Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
+% aa(Balance,:)=0;
+% aa(:,Balance)=0;
+% aa=aa+sparse(Balance,Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
 
 %暂时改一下
-t=Busnum;%相角不变
-aa(t+Balance,:)=0;
-aa(:,t+Balance)=0;
-aa(t+Balance,t+Balance)=1;
+% t=Busnum;%相角不变
+% aa(t+Balance,:)=0;
+% aa(:,t+Balance)=0;
+% aa(t+Balance,t+Balance)=1;
 dxdy=aa\yy;
 dX=dxdy(1:ContrlCount);
 dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);

func_deltF.m

@@ -11,7 +11,7 @@ deltF=[
     %zeros(busNum*2,1);
     sparse(Loadi,1,-2*(real(Vmeasurement)-V1r(Loadi)),busNum,1);
     sparse(Loadi,1,-2*(imag(Vmeasurement)-V1i(Loadi)),busNum,1);
-    -2*( [real(Imeasurement);imag(Imeasurement)]-[I1r;I1i]);
+    -0*( [real(Imeasurement);imag(Imeasurement)]-[I1r;I1i]);
 ];
 
 end

func_deltdeltF.m

@@ -4,7 +4,7 @@ deltdeltF=[
   %zeros(busNum*2,busNum*2+length(Loadi)*2);
   sparse(Loadi,Loadi,2,busNum,busNum*2+length(Loadi)*2);
   sparse(Loadi,busNum+Loadi,2,busNum,busNum*2+length(Loadi)*2);
-  zeros(length(Loadi)*2,busNum*2),2*eye(length(Loadi)*2);
+  zeros(length(Loadi)*2,busNum*2),0*eye(length(Loadi)*2);
   %sparse(Loadi,busNum*2+Loadi,2,busNum,busNum*2+length(Loadi)*2);
 ];
 end