1.删除680这个没负荷的节点。

2.修复一个50%量测权重置0的bug。
3.不考虑电压约束
4.恢复20%的边界设置。

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

FormG.m

@@ -3,6 +3,6 @@ function Mat_G=FormG(Volt,PD,QD,Loadi)
 Mat_G=[
     sparse(PD);
     sparse(QD);
-        Volt;
+%         Volt;
     ];
 end

FormLw.m

@@ -26,10 +26,10 @@ function Lw=FormLw(Mat_G,Init_U,Loadi,mPD3P,mQD3P,mVoltABCV)
 % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
 % PDU=1.8*ones(length(Loadi),1);
 % QDU=1.8*ones(length(Loadi),1);
-PDU=1.05*mPD3P;
-QDU=1.05*mQD3P;
+PDU=1.20*mPD3P;
+QDU=1.20*mQD3P;
 % t1=([PDU',QDU',1*ones(length(Volt),1)'])';
-t1=([PDU',QDU',1.1*mVoltABCV'])';
+t1=([PDU',QDU'])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
 

FormLz.m

@@ -23,10 +23,10 @@ function Lz=FormLz(Mat_G,Init_L,Loadi,mPD3P,mQD3P,mVoltABCV)
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
 % PDL=0*ones(length(Loadi),1);
 % QDL=0*ones(length(Loadi),1);
-PDL=0.95*mPD3P;
-QDL=0.95*mQD3P;
+PDL=0.8*mPD3P;
+QDL=0.8*mQD3P;
 % t1=([PDL',QDL',0.88*ones(length(Volt),1)'])';
-t1=([PDL',QDL',0.9*mVoltABCV'])';
+t1=([PDL',QDL'])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
 

OPF.m

@@ -225,15 +225,16 @@ wQD=1./(abs(mQD3P*sigma).^2);
 
 %% 没有量测量
 noLoadi=[1,5,6,10,11];
-noLoadi=Loadi;
+% noLoadi=Loadi;
 % noLoadi=[1,11];
 % noLoadi=[];
 noPQi3P=zeros(length(noLoadi)*3,1);
 noPQi3P(1:3:end)=(noLoadi-1)*3+1;
 noPQi3P(2:3:end)=(noLoadi-1)*3+2;
 noPQi3P(3:3:end)=(noLoadi-1)*3+3;
-% noLoadi=noPQi3P;
-wVolt( ismember( Loadi,noLoadi) )=0;
+noLoadi=noPQi3P;
+wVolt(setdiff(1:length(wVolt),Loadi))=0;%只有负荷处才有电压量测。
+wVolt( noLoadi)=0;
 wPD(ismember( Loadi,noLoadi))=0;
 wQD(ismember( Loadi,noLoadi))=0;
 
@@ -242,8 +243,8 @@ wQD(ismember( Loadi,noLoadi))=0;
 % wPD(5)=0;
 % wQD(5)=0;
 %%
-RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %约束条件数
-% RestraintCount=size(Loadi,1)*2; %约束条件数
+% RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %约束条件数
+RestraintCount=size(Loadi,1)*2; %约束条件数
 Init_Z=sparse(ones(1,RestraintCount));
 Init_W=sparse(-1*ones(1,RestraintCount));
 Init_L=1*sparse(ones(1,RestraintCount));
@@ -252,12 +253,13 @@ Init_Y=sparse(1,2*Busnum*3);%
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 ContrlCount=size(Loadi,1)*2+Busnum*6;
-kmax=2000;
+kmax=500;
 %%
 Precision=1e-5;
 CenterA=0.1;
 %% 加误差
 Volt=1*Vp3m;
+Volt=mVoltABCV;
 UAngel=Vp3a;
 while(abs(Gap)>Precision)
     if KK>kmax
@@ -290,8 +292,8 @@ while(abs(Gap)>Precision)
     Mat_G=FormG(Volt,PD3P,QD3P,Loadi);
     Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
     Ly=Mat_H;
-    Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P,rVoltABCV);
-    Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P,rVoltABCV);
+    Lz=FormLz(Mat_G,Init_L,Loadi,mPD3P,mQD3P,rVoltABCV);
+    Lw=FormLw(Mat_G,Init_U,Loadi,mPD3P,mQD3P,rVoltABCV);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     %     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
@@ -323,6 +325,13 @@ AME_mPD=sum(sum(abs(rPD3P-mPD3P)));
 AME_mQD=sum(sum(abs(rQD3P-mQD3P)));
 %返回收敛信息
 isConverged=1;
+if KK>=kmax
+    isConverged=0;
+end
+if abs(Gap)>Precision
+    isConverged=0;
+end
+
 fprintf('迭代次数%d\n',KK);
 toc
 end

feeder13/data1.txt

@@ -8,14 +8,15 @@
 602	684	652	800
 601	632	671	2000
 600	671	684	300
-602	671	680	1000
 600	671	692	300
 602	684	611	300
 601	692	675	500
 0
+652 20 20 0
+671 20 220 120
 684 200 200 200
-671 200 320 320
-652 150 300 200
+675 120 220 220
+634 50 130 120
 0
 634	120	90	120	90	160	110
 645	40	30	50	40	60	50

func_deltG.m

@@ -26,7 +26,7 @@ deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
        dg3_dQD,dg4_dQD,dg5_dQD;
        dg3_dx,dg4_dx,dg5_dx;
 ];
-% deltG=[dg3_dPD,dg4_dPD;
-%        dg3_dQD,dg4_dQD;
-%        dg3_dx,dg4_dx;];
+deltG=[dg3_dPD,dg4_dPD;
+       dg3_dQD,dg4_dQD;
+       dg3_dx,dg4_dx;];
 end

test.m

@@ -9,7 +9,7 @@ AME_QDSum=0;
 AME_mVoltSum=0;
 AME_mPDSum=0;
 AME_mQDSum=0;
-N=1000;
+N=500;
 for I=1:N
     [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged]=OPF();
     if isConverged==0