加入了补偿电容的考虑。

Signed-off-by: facat <facat@ipso.laptop>

dataRead.m

@@ -1,11 +1,11 @@
 function [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
-    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY] = dataRead(lineZ,dataFile )
+    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
+    cap] = dataRead(lineZ,dataFile )
 data=dlmread(dataFile);
 baseHighU=data(1,2);
 fsY0=sparse(0);
 fsY1=sparse(0);
 fsY2=sparse(0);
-
 zeroEntries=find(data(:,1)==0);
 lines=data(zeroEntries(1)+1:zeroEntries(2)-1,:);
 [setIJ,nodeNum]=numberNode(lines);
@@ -43,5 +43,20 @@ phaseCSpotLoadP( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,6);
 phaseASpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,3);
 phaseBSpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,5);
 phaseCSpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,7);
+%%  补偿电容
+cap=data(zeroEntries(2)+1:zeroEntries(3)-1,:);
+capNode=nodeNum(cap(:,1)==setIJ);
+%数据文件中的补偿电容给的是补偿容量，要将其转换为导纳。
+capB=cap(:,2:4)/1000;
+cap=struct();
+cap.capNode=capNode;
+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;
+phaseABCY=phaseABCY+sparse(offSet,offSet,1j*capB(:,2),busNum*3,busNum*3);
+offSet=(capNode-1)*3+1;
+phaseABCY=phaseABCY+sparse(offSet,offSet,1j*capB(:,3),busNum*3,busNum*3);
 end
 

run.m

@@ -1,8 +1,9 @@
 clc
 clear
-lineZ=readLineZ('D:\share\feeder123\lineParameter.txt');
+lineZ=readLineZ('.\feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
-    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY]=dataRead(lineZ,'D:\share\feeder123\data.txt');
+    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
+    cap]=dataRead(lineZ,'.\feeder13\data1.txt');
 a=exp(1j*2*pi/3);
 Tp2f=1/3*[1 1 1;
     1 a a^2;
@@ -18,8 +19,6 @@ busNum=length(phaseASpotLoadP);
 %给序电压赋初值
 Vmf1=sparse(ones(busNum,1));
 Vaf1=sparse(zeros(busNum,1));
-% Vf1=Vmf1.*exp(1j*Vaf1);
-%Vaf1=sparse(zeros(busNum,1));
 %先求解正序的
 PQi=nodeNum;
 PG=sparse(busNum,1);
@@ -49,11 +48,27 @@ fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
 %%LU分解
 [fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
 [fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
+%算初始补偿功率
 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)
     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));
+SC=VoltpC.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
+iterPD=PD+real(SA+SB+SC)/3;
+iterQD=QD+imag(SA+SB+SC)/3;
+iterPhaseASpotLoadP=phaseASpotLoadP+real(SA);
+iterPhaseBSpotLoadP=phaseBSpotLoadP+real(SB);
+iterPhaseCSpotLoadP=phaseCSpotLoadP+real(SC);
+iterPhaseASpotLoadQ=phaseASpotLoadQ+imag(SA);
+iterPhaseBSpotLoadQ=phaseBSpotLoadQ+imag(SB);
+iterPhaseCSpotLoadQ=phaseCSpotLoadQ+imag(SC);
     [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance,busNum, ...
-        PQi,PG,QG,QGi,PD,QD,Vmf1,Vaf1,fsY1amp,fsY1ang,r,c,Vf2,If2,Vf0,If0);%˛ťĆ˝şâÁż
+        PQi,PG,QG,QGi,iterPD,iterQD,Vmf1,Vaf1,fsY1amp,fsY1ang,r,c,Vf2,If2,Vf0,If0);%不平衡量
     maxD=max(abs([dP;dQ;]));
     jaco=Jacobi(Balance,busNum,QGi,Vmf1,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
     [dV, dVangle]=Solv(busNum,jaco,dP,dQ);%解出修正量
@@ -62,11 +77,11 @@ while(k<=kmax && maxD> EPS)
     %转换为三相电压
     VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%用Tp2f\ 代替Tf2p*
     VoltpA=conj(VoltpABC(1,:)');
-    CurpA=-conj((phaseASpotLoadP+1j*phaseASpotLoadQ)./VoltpA);
+    CurpA=-conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
     VoltpB=conj(VoltpABC(2,:)');
-    CurpB=-conj((phaseBSpotLoadP+1j*phaseBSpotLoadQ)./VoltpB);
+    CurpB=-conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
     VoltpC=conj(VoltpABC(3,:)');
-    CurpC=-conj((phaseCSpotLoadP+1j*phaseCSpotLoadQ)./VoltpC);
+    CurpC=-conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
     f012=Tp2f*conj([CurpA';CurpB';CurpC']);
     If0=conj(f012(1,:)');
     If1=conj(f012(2,:)');
@@ -78,12 +93,14 @@ while(k<=kmax && maxD> EPS)
     %Vf2=fsY2\If2;
     Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
     fprintf('迭代时间%f\n',toc);
+    
+    %
 end
 FortiscueToc=toc;
 fprintf('Fortiscue法计算时间 %f\n',FortiscueToc);
 Vf1=Vmf1.*exp(1j*Vaf1);
 %%
-(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;
+(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%包含补偿电容的功率
 conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
 IpABC=Tf2p*conj([If0';If1';If2']);
 %转换回三相电压
@@ -149,5 +166,7 @@ end
 NewtonToc=toc;
 fprintf('牛顿法计算时间 %f\n',NewtonToc);
 fprintf('加速比为%f\n',NewtonToc/FortiscueToc);
-
+VoltpA=Vp3m(1:3:end).*exp(1j*Vp3a(1:3:end));
+VoltpB=Vp3m(2:3:end).*exp(1j*Vp3a(2:3:end));
+VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end));