把原来的OPF重新加入这里面，收敛了，但总觉得有地方不太对。

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

FormG.m

@@ -1,8 +1,14 @@
 function Mat_G=FormG(I1r,I1i)
-% Mat_G=[
+%t1=PG(PVi);
-%     I1r;
+%GP=t1;%发电机P
-%     I1i;
+%GP=[4.5 4.5]';
-%     ];
+%%线路
+%发电机Q
+% t1=Volt'*Volt;
+% t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
+% t3=t1.*t2;
+% t4=sum(t3,2);%发电机Q
+%GQ=t4;
 Mat_G=[
     I1r.^2+I1i.^2;
     ];

FormH.m

@@ -1,15 +1,6 @@
 function Mat_H=FormH(fsY11,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance)
-%% 
+
-%QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt');
+%%%%一下是学姐给的公式
-%QD(QD~=0)=PD(QD~=0)./tan(QDcos);
-%QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
-%% 
-%PD(Loadi)=QD(Loadi)./tan(acos(0.98));
-% AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
-% dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt;
-% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt;
-% 
-% Mat_H=[dP;dQ;];
 H=[
     real(fsY11),-imag(fsY11);
     imag(fsY11),real(fsY11);
@@ -21,4 +12,5 @@ Mat_H=H*[V1r;V1i];
 Mat_H=Mat_H+sparse(Loadi,1,-I1r,busNum*2,1)+sparse(Loadi+busNum,1,-I1i,busNum*2,1);
 Mat_H(Balance)=Mat_H(Balance)-BalI1r;
 Mat_H(Balance+busNum)=Mat_H(Balance+busNum)-BalI1i;
+
 end

FormLw.m

@@ -1,12 +1,5 @@
 function Lw=FormLw(Loadi,Mat_G,Init_U)
-% VoltU=sparse(1.2*ones(Busnum,1));
+
-% PDU=sparse(2*ones(length(Loadi),1));
-% QDU=sparse(2*ones(length(Loadi),1));
-% t1=[PDU;
-%     QDU;
-%     VoltU];
-% t2=Mat_G+Init_U-t1;
-% Lw=t2;
 upper=2*sparse(ones(length(Loadi)*1,1));
 Lw=Mat_G+Init_U-upper;
 end

FormLz.m

@@ -1,14 +1,5 @@
 function Lz=FormLz(Loadi,Mat_G,Init_L)
-% VoltL=sparse(0.8*ones(Busnum,1));
+
-% PDL=sparse(-2*ones(length(Loadi),1));
-% QDL=sparse(-2*ones(length(Loadi),1));
-% t1=[PDL;
-%     QDL;
-%     VoltL
-%     ];
-% t2=Mat_G-Init_L-t1;
-% Lz=t2;
 lower=0*sparse(ones(length(Loadi)*1,1));
 Lz=Mat_G-Init_L-lower;
-% Lz=-[]
 end

IPMLoop.m

@@ -2,8 +2,8 @@ function [ V1r,V1i,I1r,I1i ] = IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fs
 %把每个序的循环写在这个函数中。其实也就是内点法循环。
 V1r=1*ones(busNum,1);
 V1i=0*ones(busNum,1);
-I1r=0.01*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
+I1r=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
-I1i=0.01*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
+I1i=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
 KK=0;
 plotGap=zeros(1,60);
 %初始化
@@ -12,32 +12,36 @@ RestraintCount=length(Loadi)*1;
 ContrlCount=busNum*2+length(Loadi)*2;
 CenterA=0.1;
 Init_Z=sparse(ones(RestraintCount,1));
-Init_W=sparse(-1*ones(RestraintCount,1));
+Init_W=sparse(-1.0*ones(RestraintCount,1));
 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);
-% PG=sparse(busNum,1);
+kmax=60;
-% PG(Balance)=0.1105;
+while(abs(Gap)>0.000001)
-% QG=sparse(busNum,1);
+    if KK>kmax
-% QG(Balance)=0.0984;
+        break;
-% SEPD(2)=0.1105;
+    end
-% SEQD(2)=0.0984;
+    plotGap(KK+1)=Gap;
-while Gap>1e-5 && KK<20
-    KK=KK+1;
     Init_u=Gap/2/RestraintCount*CenterA;
+    %% 开始计算OPF
+    %% 形成等式约束的雅克比
     deltH=func_deltH(busNum,fsY1,Loadi,Balance);
+    %% 形成不等式约束的雅克比
     deltG=func_deltG(busNum,Loadi,I1r,I1i);
-%     deltG=0;
+    %%
-%     L_1Z=diag(Init_Z./Init_L);
+    L_1Z=diag(Init_Z./Init_L);
-%     U_1W=diag(Init_W./Init_U);
+    U_1W=diag(Init_W./Init_U);
+    %% 形成海森阵
     deltdeltF=func_deltdeltF(busNum,fsY1,Loadi);
-%     deltdeltF=0;
+    %% 形成ddHy
-%     ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount);
+%     ddh=func_ddh(busNum,Loadi,Init_Z,Init_W);
     ddh=0;
+    %%  开始构建ddg
     ddg=func_ddg(busNum,Loadi,Init_Z,Init_W);
+    %%  开始构建deltF
     deltF=func_deltF(measurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i);
-%     deltF=0;
+    %%
     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);
@@ -46,14 +50,15 @@ while Gap>1e-5 && KK<20
     Lz=FormLz(Loadi,Mat_G,Init_L);
     Lw=FormLw(Loadi,Mat_G,Init_U);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-%     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
+    YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
-    plotGap(KK)=Gap;
+    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);
-    fprintf('迭代次数 %d Gap %f\n',KK,plotGap(KK));
+    %%取各分量
-    XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi);
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum);
-    [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);
+    [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);
     Gap=(Init_L'*Init_Z-Init_U'*Init_W);
+    fprintf('Gap %f\n',full(Gap));
+    KK=KK+1;
 end
 f=sum(([real(measurement);imag(measurement)]-[-I1r;-I1i]).^2);
 

Modification.m

@@ -1,22 +1,27 @@
-function [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, ... 
+function [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)
-    V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
+% fprintf('AlphaP %f\n',full(AlphaP));
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
+% fprintf('AlphaD %f\n',full(AlphaD));
+
 Init_Z=Init_Z+AlphaD*deltZ;
 Init_L=Init_L+AlphaP*deltL;
 Init_W=Init_W+AlphaD*deltW;
 Init_U=Init_U+AlphaP*deltU;
 Init_Y=Init_Y+AlphaD*deltY;
-t=deltX(1:busNum*2);
+%PG(PGi)=PG(PGi)+deltX(1:size(PGi,1));
-V1r=V1r+AlphaP*t(1:busNum);
+% PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
-V1r(Balance)=Vref;
+V1r=V1r+AlphaP*deltX(1:Busnum);
-V1i=V1i+AlphaP*t(busNum+1:end);
+V1r(Balance)=1;
+V1i=V1i+AlphaP*deltX(Busnum+1:2*Busnum);
 V1i(Balance)=0;
-t=deltX(busNum*2+1:ContrlCount);
+%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);
 I1r=I1r+AlphaP*t(1:length(Loadi));
 I1i=I1i+AlphaP*t(length(Loadi)+1:end);
-%balVolt=1;%Volt(Balance);
+
-%Volt=Volt+AlphaP*t(1:Busnum);
+%Volt=Volt+AlphaP*t(2:2:2*Busnum);ÔÝʱ¸ÄÒ»ÏÂ20111227
-%Volt(Balance)=balVolt;
+%UAngel=UAngel+AlphaP*t(1:2:2*Busnum);ÔÝʱ¸ÄÒ»ÏÂ20111227
-%UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum);
+
 end

SolveIt.m

@@ -1,31 +1,26 @@
-function XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,Busnum,Loadi)
+function XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,ConstraintCount,Lx,Balance,Busnum)
 LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
-H=-deltdeltF+ddh+ddg;
+H=-deltdeltF+ddh+ddg;%+ddg*(Init_Z'+Init_W');
-t1=diag(Init_Z./Init_L-Init_W./Init_U);
+t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
 t2=-deltG*( t1 )*deltG';
 aa=[
     (H+t2),deltH;
-    deltH',zeros(length(Init_Y));
+    deltH',sparse(length(Init_Y),length(Init_Y));
     ];
 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);
-deltG(Balance,:)=0;
+
-%% 平衡节点电压虚部
+%暂时改一下
-t=Busnum*1;
+t=Busnum;%相角不变
 aa(t+Balance,:)=0;
 aa(:,t+Balance)=0;
-aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
+aa(t+Balance,t+Balance)=1;
-deltG(t+Balance,:)=0;
-%%
 dxdy=aa\yy;
-%% KLU
-%spy(aa)
-%dxdy = klu(aa,'\',full(yy));
-%%
 dX=dxdy(1:ContrlCount);
 dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);
 dL=Lz+deltG'*dX;

feeder13/data1.txt

@@ -5,8 +5,8 @@
 0
 632 30 0 0
 0
-632 120 120 120 110 90 90
 631 90 100 110 80 90 90
+632 120 120 120 110 90 90
 0
 
 

func_ddg.m

@@ -1,5 +1,5 @@
 function ddg=func_ddg(busNum,Loadi,Init_Z,Init_W)
-% ddg=0;
+
 Init_ZW=Init_Z+Init_W;
 t=[2*diag( ones(length(Loadi),1).*Init_ZW),zeros(length(Loadi));
     zeros(length(Loadi)),2*diag( ones(length(Loadi),1).*Init_ZW);

func_deltF.m

@@ -1,15 +1,15 @@
 function deltF=func_deltF(measurement,busNum,fsY11,Loadi,V1r,V1i,I1r,I1i)
-% H=[
+%t1=PG(setdiff(PVi,Balance));
-%     real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2);
+% t2=Volt'*Volt;
-%     imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2);
+% t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-%     zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi));
+% t4=t2.*t3;
-%     zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi));
+% t5=sum(t4,2);
-% ];
+% PBal=t5(Balance);
-% X=[V1r;V1i;I1r;I1i];
+% PPG=([PQ(1),PBal])';%暂时用土办法处理一下
-% deltF=-H'*(measurement-H*X);
+%% 
 deltF=[
     zeros(busNum*2,1);
     -2*( [real(measurement);imag(measurement)]-[I1r;I1i]);
 ];
-% deltF=0;
+
 end

func_deltG.m

@@ -1,19 +1,4 @@
 function deltG=func_deltG(busNum,Loadi,I1r,I1i)
-% dgV_dPD=sparse(length(Loadi),busNum);
-% dgV_dQD=sparse(length(Loadi),busNum);
-% dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum);
-%     sparse(busNum,busNum);
-%     ];
-% %
-% dgPD_dPD=sparse(eye(length(Loadi)));
-% dgPD_dQD=sparse(length(Loadi),length(Loadi));
-% dgPD_dx=sparse(busNum*2,length(Loadi));
-% %
-% dgQD_dPD=sparse(length(Loadi),length(Loadi));
-% dgQD_dQD=sparse(eye(length(Loadi)));
-% dgQD_dx=sparse(busNum*2,length(Loadi));
-% %% 
-
 deltG=[
     zeros(busNum*2,length(Loadi)*1);
     2*eye(length(Loadi)*1)*diag(I1r);

func_deltH.m

@@ -1,8 +1,4 @@
 function deltH=func_deltH(busNum,fsY11,Loadi,Balance)
-% dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
-% dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
-% dH_dx = jacobian(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
-% deltH=[dH_dPD;dH_dQD;dH_dx'];
 H=[
     real(fsY11),-imag(fsY11);
     imag(fsY11),real(fsY11);

func_deltdeltF.m

@@ -1,22 +1,7 @@
 function deltdeltF=func_deltdeltF(busNum,fsY11,Loadi)
-% H=[
-%     real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2);
-%     imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2);
-%     zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi));
-%     zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi));
-% ];
 
-% H=[
-%     real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2);
-%     imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2);
-%     zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi));
-%     zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi));
-% ];
-
-% deltdeltF=H'*H;
 deltdeltF=[
   zeros(busNum*2,busNum*2+length(Loadi)*2);
   zeros(length(Loadi)*2,busNum*2),2*eye(length(Loadi)*2);
 ];
-% deltdeltF=0;
 end