1.修复了三相有电容计算错误的bug。

2.找到一个电压还算合理的算例，不过迭代次数太多。

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

FormG.m

@@ -1,7 +1,8 @@
-function Mat_G=FormG(PD,QD,Loadi)
+function Mat_G=FormG(Volt,PD,QD,Loadi)
 
 Mat_G=[
     sparse(PD);
     sparse(QD);
+        Volt;
     ];
 end

FormLw.m

@@ -1,4 +1,4 @@
-function Lw=FormLw(Mat_G,Init_U,Loadi,mPD3P,mQD3P)
+function Lw=FormLw(Mat_G,Init_U,Loadi,mPD3P,mQD3P,mVoltABCV)
 % KK=999;
 % %PU=GenU(:,2);%发电机有功上界
 % 
@@ -26,9 +26,10 @@ function Lw=FormLw(Mat_G,Init_U,Loadi,mPD3P,mQD3P)
 % 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.2*mPD3P;
+PDU=1.05*mPD3P;
-QDU=1.2*mQD3P;
+QDU=1.05*mQD3P;
-t1=([PDU',QDU'])';
+% t1=([PDU',QDU',1*ones(length(Volt),1)'])';
+t1=([PDU',QDU',1.1*mVoltABCV'])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
 

FormLz.m

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

OPF.m

@@ -41,7 +41,7 @@ PD=Pabc/3;
 QD=Qabc/3;
 Loadi=find(PD~=0);
 maxD=100000;% 最大不平衡量
-EPS=1e-5;
+EPS=1e-6;
 k=0;
 kmax=20;
 fsY11=fsY1;
@@ -68,12 +68,12 @@ tic
 VoltpA=sparse(ones(busNum,1));
 VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
 VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
-while(k<=kmax && maxD> EPS)
+while(k<=kmax+10 && maxD> EPS)
     k=k+1;
     %把补偿电容看作负荷
     SA=VoltpA.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,1),busNum,1));
-    SB=VoltpB.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1));
+    SB=VoltpB.*conj(VoltpB.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1));
-    SC=VoltpC.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
+    SC=VoltpC.*conj(VoltpC.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
     %先不要电容
     %     SA=0;
     %     SB=0;
@@ -180,6 +180,7 @@ Balance3P(3:3:end)=(Balance-1)*3+3;
 % QGC=sum(QD3P(3:3:end))*1.03;
 % PG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[PGA,PGB,PGC],Busnum*3,1);
 % QG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[QGA,QGB,QGC],Busnum*3,1);
+
 PG3P=real(PGQG);
 QG3P=imag(PGQG);
 
@@ -210,7 +211,7 @@ rQD3P(3:3:end)=phaseCSpotLoadQ;
 rPD3P=rPD3P(Loadi);
 rQD3P=rQD3P(Loadi);
 %量测量
-sigma=0.01;
+sigma=0.03;
 VoltSigma=(1+normrnd(0,sigma,length(rVoltABCV),1));
 mVoltABCV=rVoltABCV.*VoltSigma;
 PD3PSigma=(1+normrnd(0,sigma,length(rPD3P),1));
@@ -221,11 +222,28 @@ mQD3P=rQD3P.*QD3PSigma;
 wVolt=1./(abs(mVoltABCV*sigma).^2);
 wPD=1./(abs(mPD3P*sigma).^2);
 wQD=1./(abs(mQD3P*sigma).^2);
+
+%% 没有量测量
+noLoadi=[1,5,6,10,11];
+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;
+wPD(ismember( Loadi,noLoadi))=0;
+wQD(ismember( Loadi,noLoadi))=0;
+
+
 % wVolt(5)=0;
 % wPD(5)=0;
 % wQD(5)=0;
 %%
-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));
@@ -234,12 +252,12 @@ 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=20;
+kmax=2000;
 %%
 Precision=1e-5;
 CenterA=0.1;
 %% 加误差
-Volt=Vp3m;
+Volt=1*Vp3m;
 UAngel=Vp3a;
 while(abs(Gap)>Precision)
     if KK>kmax
@@ -256,7 +274,7 @@ while(abs(Gap)>Precision)
 %     L_1Z=diag(Init_Z./Init_L);
 %     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount);
+    deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount,Loadi);
 %         deltdeltF=0;
     %% 形成ddHy
     ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Yangle,Loadi,ContrlCount);
@@ -269,11 +287,11 @@ while(abs(Gap)>Precision)
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    Mat_G=FormG(PD3P,QD3P,Loadi);
+    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,mPD3P,mQD3P);
+    Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P,rVoltABCV);
-    Lw=FormLw(Mat_G,Init_U,Loadi,mPD3P,mQD3P);
+    Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P,rVoltABCV);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     %     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程

checkSSatisfied.m

@@ -17,7 +17,7 @@ pLoadABC(2:3:end)=phaseBSpotLoadP+1j*phaseBSpotLoadQ;
 pLoadABC(3:3:end)=phaseCSpotLoadP+1j*phaseCSpotLoadQ;
 ck=Sp3+pLoadABC;
 ck(3*(Balance-1)+1:3*(Balance-1)+3)=0;
-if any(abs(ck)>1e-5)
+if any(real(ck)>1e-5) ||any(imag(ck)>1e-5)
     fprintf('反推回的功率不匹配。潮流方程约束不满足。\n');
 end
 output_args=ck;

dataRead.m

@@ -46,7 +46,7 @@ phaseCSpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,7);
 %%  补偿电容
 cap=data(zeroEntries(2)+1:zeroEntries(3)-1,:);
 if ~isempty(cap)
-capNode=nodeNum(cap(:,1)==setIJ);
+capNode=nodeNum(ismember(setIJ,cap(:,1)));
 %数据文件中的补偿电容给的是补偿容量，要将其转换为导纳。
 capB=cap(:,2:4)/1000;
 cap=struct();
@@ -55,9 +55,9 @@ cap.capB=capB;
 %把补偿电容加到矩阵中。
 offSet=(capNode-1)*3+1;
 phaseABCY=phaseABCY+sparse(offSet,offSet,1j*capB(:,1),busNum*3,busNum*3);
-offSet=(capNode-1)*3+1;
+offSet=(capNode-1)*3+2;
 phaseABCY=phaseABCY+sparse(offSet,offSet,1j*capB(:,2),busNum*3,busNum*3);
-offSet=(capNode-1)*3+1;
+offSet=(capNode-1)*3+3;
 phaseABCY=phaseABCY+sparse(offSet,offSet,1j*capB(:,3),busNum*3,busNum*3);
 else
     cap=struct();

feeder13/data1.txt

@@ -13,14 +13,18 @@
 602	684	611	300
 601	692	675	500
 0
-632 0 0 0
+684 200 200 200
+671 200 320 320
+652 150 300 200
 0
-634 160 110 120 90 120 90
+634	120	90	120	90	160	110
-645 123 102 170 125 133 110
+645	40	30	50	40	60	50
-646 172 150 230 132 211 160
+646	100	50	70	40	60	40
-652 128 86 30 10 80 70
+652	20	10	60	30	40	40
-675 85 10 68 60 20 12
+671	385	220	385	220	385	220
-611 52 42 85 60 70 40
+675	68	60	290	212	285	190
+692	60	50	40	40	70	40
+611	80	20	60	40	30	20
 0
 600	632	645	500
 602	632	633	500
@@ -46,5 +50,11 @@
 692 142 140 153 121 170 151
 611 156 142 185 60 170 80
 0
-
+634 160 110 120 90 120 90
+645 123 102 170 125 133 110
+646 17 15 23 13 21 16
+652 128 86 30 10 80 70
+675 85 10 68 60 20 12
+611 52 42 85 60 70 40
+680 12 4 10 6 7 4
 

func_deltG.m

@@ -8,26 +8,25 @@ sizeLoadi=size(Loadi,1);
 %%
 dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
 dg4_dPD=sparse(size(Loadi,1),length(Loadi));
-% dg5_dPD=sparse(size(Loadi,1),Busnum*3);
+dg5_dPD=sparse(size(Loadi,1),Busnum*3);
 %% 
 dg3_dQD=sparse(length(Loadi),length(Loadi));
 dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
-% dg5_dQD=sparse(size(Loadi,1),Busnum*3);
+dg5_dQD=sparse(size(Loadi,1),Busnum*3);
 
 %% 
 dg3_dx=sparse(2*Busnum*3,sizeLoadi);
 dg4_dx=sparse(2*Busnum*3,length(Loadi));
-% dg5_dx=[sparse(1:Busnum*3,1:Busnum*3,zeros(Busnum*3,1),Busnum*3,Busnum*3);
+dg5_dx=[sparse(1:Busnum*3,1:Busnum*3,ones(Busnum*3,1),Busnum*3,Busnum*3);
-%     sparse(Busnum*3,Busnum*3);
+    sparse(Busnum*3,Busnum*3);
-%     ];
+    ];
 %% 
 
-% deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
+deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
-%        dg3_dQD,dg4_dQD,dg5_dQD;
+       dg3_dQD,dg4_dQD,dg5_dQD;
-%        dg3_dx,dg4_dx,dg5_dx;
+       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

func_deltdeltF.m

@@ -1,4 +1,4 @@
-function deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount)
+function deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount,Loadi)
 
 %ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
 C=[wPD' wQD'];

test.m

@@ -6,6 +6,9 @@ AME_VAngleSum=0;
 JMeasurementSum=0;
 AME_PDSum=0;
 AME_QDSum=0;
+AME_mVoltSum=0;
+AME_mPDSum=0;
+AME_mQDSum=0;
 N=1000;
 for I=1:N
     [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged]=OPF();
@@ -18,9 +21,15 @@ for I=1:N
     AME_VAngleSum=AME_VAngleSum+AME_VAngle;
     AME_PDSum=AME_PDSum+AME_PD;
     AME_QDSum=AME_QDSum+AME_QD;
+    AME_mVoltSum=AME_mVoltSum+AME_mVolt;
+    AME_mPDSum=AME_mPDSum+AME_mPD;
+    AME_mQDSum=AME_mQDSum+AME_mQD;
 end
 JMeasurementSum=JMeasurementSum/N;
 AME_VoltSum=AME_VoltSum/N;
 AME_VAngleSum=AME_VAngleSum/N;
 AME_PDSum=AME_PDSum/N;
 AME_QDSum=AME_QDSum/N;
+AME_mVoltSum=AME_mVoltSum/N;
+AME_mPDSum=AME_mPDSum/N;
+AME_mQDSum=AME_mQDSum/N;