编写了基于牛顿法的三相配网潮流

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

.gitignore

@@ -1,3 +1,4 @@
 *.asv
 /feeder13
-*.mexw32
+*.mexw32
+/DistributionNetwork-OnlyFortescue

Unbalance.m

@@ -1,11 +1,5 @@
 function [dP, dQ, YdotSinVmf1, YdotCosVmf1, diag_Vmf1_YdotSin, diag_Vmf1_YdotCos]=Unbalance(Balance,busNum, ... 
     PQi,PG,QG,QGi,PD,QD,Vmf1,Vaf1,Y,Yangle,r,c,Vf2,If2,Vf0,If0)
-%real(diag(Vmf1.*exp(1j*Vaf1))*(fsY1*(Vmf1.*exp(1j*Vaf1))));
-% Y=abs(fsY1);
-% [r,c,Yangle]=find(fsY1);
-% Yangle=angle(Yangle);
-% Vmf1=abs(Vf1);
-% Vaf1=angle(Vf1);
 t1=sparse(r,c,Vaf1(r)-Vaf1(c) -Yangle,busNum,busNum) ;
 YdotSin=Y.* ( spfun(@sin,t1) );
 YdotCos=Y.* ( spfun (@cos, t1 ) );
@@ -19,9 +13,6 @@ diag_Vmf1_YdotCosVmf1=diag_Vmf1_YdotCos*Vmf1;
 diag_Vmf1_YdotSinVmf1=diag_Vmf1_YdotSin*Vmf1;
 dP=diag_Vmf1_YdotCosVmf1+orderedPDPG+real(Vf2.*conj(If2)+Vf0.*conj(If0));
 dQ=diag_Vmf1_YdotSinVmf1+orderedQDQG+imag(Vf2.*conj(If2)+Vf0.*conj(If0));
-%orderedPDPG+real(Vf2.*conj(If2)+Vf0.*conj(If0))+real(diag(Vf1)*conj(fsY1*Vf1));
-%dP=diag_Vmf1_YdotCosVmf1+orderedPDPG;
-%dQ=diag_Vmf1_YdotSinVmf1+orderedQDQG;
 dP(Balance)=0; 
 dQ(QGi)=0;  
 end

run.m

@@ -92,6 +92,53 @@ disp([setIJ,nodeNum ])
 ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
     phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
-fprintf('最大不平衡量为%f\n',full(max(abs(ub))))
+fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
+fprintf('开始牛顿法迭代\n');
+%%用牛顿法求解
+[r,c,GB]=find(phaseABCY);
+Y=abs(phaseABCY);
+Yangle=angle(GB);
+Vp3=sparse(ones(busNum*3,1));%给电压赋初值
+Vp3(2:3:end)=Vp3(2:3:end)*exp(1j*-120/180*pi);
+Vp3(3:3:end)=Vp3(3:3:end)*exp(1j*+120/180*pi);
+PQi3P=zeros(length(PQi)*3,1);
+PQi3P(1:3:end)=(PQi-1)*3+1;
+PQi3P(2:3:end)=(PQi-1)*3+2;
+PQi3P(3:3:end)=(PQi-1)*3+3;
+PG=0;
+QG=0;
+PD3P=sparse(busNum*3,1);
+QD3P=sparse(busNum*3,1);
+PD3P(1:3:end)=phaseASpotLoadP;
+PD3P(2:3:end)=phaseBSpotLoadP;
+PD3P(3:3:end)=phaseCSpotLoadP;
+QD3P(1:3:end)=phaseASpotLoadQ;
+QD3P(2:3:end)=phaseBSpotLoadQ;
+QD3P(3:3:end)=phaseCSpotLoadQ;
+QGi3P=zeros(length(QGi)*3,1);
+QGi3P(1:3:end)=(QGi-1)*3+1;
+QGi3P(2:3:end)=(QGi-1)*3+2;
+QGi3P(3:3:end)=(QGi-1)*3+3;
+Vp3m=abs(Vp3);
+Vp3a=angle(Vp3);
+Balance3P=zeros(length(Balance)*3,1);
+Balance3P(1:3:end)=(Balance-1)*3+1;
+Balance3P(2:3:end)=(Balance-1)*3+2;
+Balance3P(3:3:end)=(Balance-1)*3+3;
+Vp3a((Balance-1)*3+1)=0;
+Vp3a((Balance-1)*3+2)=-120/180*pi;
+Vp3a((Balance-1)*3+3)=+120/180*pi;
+k=0;
+maxD=10000;
+while(k<=kmax && maxD> EPS)
+    k=k+1;
+[dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance3P,busNum*3, ... 
+    PQi3P,PG,QG,QGi3P,PD3P,QD3P,Vp3m,Vp3a,Y,Yangle,r,c,0,0,0,0);
+maxD=max(abs([dP;dQ;]));
+jaco=Jacobi(Balance3P,busNum*3,QGi3P,Vp3m,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
+[dV, dVangle]=Solv(busNum*3,jaco,dP,dQ);%解出修正量
+[Vp3m, Vp3a]=Modify(Vp3m,Vp3a,dV,dVangle,1);
+fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
+end