From e00c8da8d514ad3903ed5a37a66b2f4264ee1682 Mon Sep 17 00:00:00 2001
From: facat <facat@ipso.laptop>
Date: Tue, 22 Apr 2014 15:50:28 +0800
Subject: [PATCH] =?UTF-8?q?=E7=BC=96=E5=86=99=E4=BA=86=E5=9F=BA=E4=BA=8E?=
 =?UTF-8?q?=E7=89=9B=E9=A1=BF=E6=B3=95=E7=9A=84=E4=B8=89=E7=9B=B8=E9=85=8D?=
 =?UTF-8?q?=E7=BD=91=E6=BD=AE=E6=B5=81?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: facat <facat@ipso.laptop>
---
 .gitignore  |  3 ++-
 Unbalance.m |  9 ---------
 run.m       | 49 ++++++++++++++++++++++++++++++++++++++++++++++++-
 3 files changed, 50 insertions(+), 11 deletions(-)

diff --git a/.gitignore b/.gitignore
index b6392e9..8f1013f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 *.asv
 /feeder13
-*.mexw32
\ No newline at end of file
+*.mexw32
+/DistributionNetwork-OnlyFortescue
\ No newline at end of file
diff --git a/Unbalance.m b/Unbalance.m
index 014670c..e8bfcf6 100644
--- a/Unbalance.m
+++ b/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
\ No newline at end of file
diff --git a/run.m b/run.m
index ec990da..a6e6a93 100644
--- a/run.m
+++ b/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