From 99aab5404c3b947e1a6236786dbe7e8061f1e16b Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Wed, 21 May 2014 22:26:48 +0800 Subject: [PATCH] =?UTF-8?q?=E6=8A=8A=E4=B8=80=E6=AC=A1=E8=BF=AD=E4=BB=A3?= =?UTF-8?q?=E8=BF=87=E7=A8=8B=E7=9A=84=E4=BB=A3=E7=A0=81=E5=A4=8D=E5=88=B6?= =?UTF-8?q?=E8=BF=87=E6=9D=A5=E4=BA=86=EF=BC=8C=E8=BF=98=E6=B2=A1=E6=A3=80?= =?UTF-8?q?=E6=9F=A5=E6=98=AF=E4=B8=8D=E6=98=AF=E6=9C=89=E5=A4=8D=E5=88=B6?= =?UTF-8?q?=E9=94=99=E7=9A=84=E5=9C=B0=E6=96=B9=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- .gitignore | 3 +- AssignXX.m | 23 ++++ DistributionNetwork-OnlyFortescue | 1 + FormAlphaD.m | 9 ++ FormAlphaP.m | 9 ++ FormG.m | 7 ++ FormH.m | 14 +++ FormLw.m | 10 ++ FormLx.m | 4 + FormLxComa.m | 11 ++ FormLz.m | 12 +++ FormYY.m | 10 ++ Modification.m | 18 ++++ SolveIt.m | 46 ++++++++ func_ddg.m | 4 + func_ddh.m | 48 +++++++++ func_deltF.m | 8 ++ func_deltG.m | 22 ++++ func_deltH.m | 7 ++ func_deltdeltF.m | 5 + jacobian.m | 20 ++++ run.m | 168 +++++++++++++++++++----------- 22 files changed, 400 insertions(+), 59 deletions(-) create mode 100644 AssignXX.m create mode 160000 DistributionNetwork-OnlyFortescue create mode 100644 FormAlphaD.m create mode 100644 FormAlphaP.m create mode 100644 FormG.m create mode 100644 FormH.m create mode 100644 FormLw.m create mode 100644 FormLx.m create mode 100644 FormLxComa.m create mode 100644 FormLz.m create mode 100644 FormYY.m create mode 100644 Modification.m create mode 100644 SolveIt.m create mode 100644 func_ddg.m create mode 100644 func_ddh.m create mode 100644 func_deltF.m create mode 100644 func_deltG.m create mode 100644 func_deltH.m create mode 100644 func_deltdeltF.m create mode 100644 jacobian.m diff --git a/.gitignore b/.gitignore index 8f1013f..d1b71e6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.asv /feeder13 *.mexw32 -/DistributionNetwork-OnlyFortescue \ No newline at end of file +#/DistributionNetwork-OnlyFortescue +*.7z \ No newline at end of file diff --git a/AssignXX.m b/AssignXX.m new file mode 100644 index 0000000..2ff0f36 --- /dev/null +++ b/AssignXX.m @@ -0,0 +1,23 @@ +function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum) +% deltX=XX(1:14); +% deltY=XX(15:24); +% deltZ=XX(25:38); +% deltW=XX(39:52); +% deltL=XX(53:66); +% deltU=XX(67:80); +deltX=XX(1:ContrlCount); +k1=ContrlCount+2*Busnum; +deltY=XX(ContrlCount+1:k1); +k2=k1; +k1=k2+RestraintCount; +deltZ=XX(k2+1:k1); +k2=k1; +k1=k2+RestraintCount; +deltW=XX(k2+1:k1); +k2=k1; +k1=k2+RestraintCount; +deltL=XX(k2+1:k1); +k2=k1; +k1=k2+RestraintCount; +deltU=XX(k2+1:k1); +end \ No newline at end of file diff --git a/DistributionNetwork-OnlyFortescue b/DistributionNetwork-OnlyFortescue new file mode 160000 index 0000000..ecaed83 --- /dev/null +++ b/DistributionNetwork-OnlyFortescue @@ -0,0 +1 @@ +Subproject commit ecaed830b547d025369bd6137707f6a96856556a diff --git a/FormAlphaD.m b/FormAlphaD.m new file mode 100644 index 0000000..bd264d2 --- /dev/null +++ b/FormAlphaD.m @@ -0,0 +1,9 @@ +function AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW) +tdeltZinx=find(deltZ<0); +tdeltWinx=find(deltW>0); +t1=-Init_Z(tdeltZinx)./deltZ(tdeltZinx); +t2=-Init_W(tdeltWinx)./deltW(tdeltWinx); +t3=[t1,t2]; +t4=min(t3); +AlphaD=0.9995*min([t4 1]); +end \ No newline at end of file diff --git a/FormAlphaP.m b/FormAlphaP.m new file mode 100644 index 0000000..a3cf2be --- /dev/null +++ b/FormAlphaP.m @@ -0,0 +1,9 @@ +function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU) +tdeltLinx=find(deltL<0); +tdeltUinx=find(deltU<0); +t1=-Init_L(tdeltLinx)./deltL(tdeltLinx); +t2=-Init_U(tdeltUinx)./deltU(tdeltUinx); +t3=[t1,t2]; +t4=min(t3); +AlphaP=0.9995*min([t4 1]); +end \ No newline at end of file diff --git a/FormG.m b/FormG.m new file mode 100644 index 0000000..fb8bbbf --- /dev/null +++ b/FormG.m @@ -0,0 +1,7 @@ +function Mat_G=FormG(Volt,PD,QD) +Mat_G=[ + sparse(PD); + sparse(QD); + Volt; + ]; +end \ No newline at end of file diff --git a/FormH.m b/FormH.m new file mode 100644 index 0000000..7ce3fba --- /dev/null +++ b/FormH.m @@ -0,0 +1,14 @@ +function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle) +%% +%QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt'); +%QD(QD~=0)=PD(QD~=0)./tan(QDcos); +%QD(QD_NON_ZERO_IND)=QD_NON_ZERO; +%% +%PD(Loadi)=QD(Loadi)./tan(acos(0.98)); +AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); +dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt; +dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt; + +Mat_H=[dP;dQ;]; + +end \ No newline at end of file diff --git a/FormLw.m b/FormLw.m new file mode 100644 index 0000000..4685222 --- /dev/null +++ b/FormLw.m @@ -0,0 +1,10 @@ +function Lw=FormLw(Mat_G,Init_U,Busnum,Loadi) +VoltU=sparse(1.1*ones(Busnum,1)); +PDU=sparse(1.5*ones(length(Loadi),1)); +QDU=sparse(1.5*ones(length(Loadi),1)); +t1=[PDU; + QDU; + VoltU]; +t2=Mat_G+Init_U-t1; +Lw=t2; +end \ No newline at end of file diff --git a/FormLx.m b/FormLx.m new file mode 100644 index 0000000..826652b --- /dev/null +++ b/FormLx.m @@ -0,0 +1,4 @@ +function Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W) +t1=deltF-deltH*Init_Y-deltG*(Init_Z+Init_W); +Lx=t1; +end \ No newline at end of file diff --git a/FormLxComa.m b/FormLxComa.m new file mode 100644 index 0000000..9d25a12 --- /dev/null +++ b/FormLxComa.m @@ -0,0 +1,11 @@ +function LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx) +%t1=deltG*(Init_Z'+Init_W');%% +t2=Lul+diag(Init_Z)*Lz; +t3=inv(diag(Init_L)); +t4=t3*t2;% +t5=Luu-diag(Init_W)*Lw; +t6=inv(diag(Init_U)); +t7=t6*t5;% +t8=deltG*(t4+t7);%% +LxComa=Lx+t8; +end \ No newline at end of file diff --git a/FormLz.m b/FormLz.m new file mode 100644 index 0000000..5c373fb --- /dev/null +++ b/FormLz.m @@ -0,0 +1,12 @@ +function Lz=FormLz(Mat_G,Init_L,Busnum,Loadi) +VoltL=sparse(0.9*ones(Busnum,1)); +PDL=sparse(0*ones(length(Loadi),1)); +QDL=sparse(0*ones(length(Loadi),1)); +t1=[PDL; + QDL; + VoltL + ]; +t2=Mat_G-Init_L-t1; +Lz=t2; + +end \ No newline at end of file diff --git a/FormYY.m b/FormYY.m new file mode 100644 index 0000000..cbad106 --- /dev/null +++ b/FormYY.m @@ -0,0 +1,10 @@ +function YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx) +YY=[ + Lx; + -Ly; + -Lz; + -Lw; + -Lul; + -Luu; +]; +end \ No newline at end of file diff --git a/Modification.m b/Modification.m new file mode 100644 index 0000000..d55962c --- /dev/null +++ b/Modification.m @@ -0,0 +1,18 @@ +function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD,QD,Loadi) +AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); +AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); +Init_Z=Init_Z+AlphaD*deltZ; +Init_L=Init_L+AlphaP*deltL; +Init_W=Init_W+AlphaD*deltW; +Init_U=Init_U+AlphaP*deltU; +Init_Y=Init_Y+AlphaD*deltY; +t=deltX(1:size(Loadi,1)*2); +PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi)); +QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:length(Loadi)*2); +t=deltX(size(Loadi,1)*2+1:ContrlCount); +t(Busnum+Balance)=0; +balVolt=Volt(Balance); +Volt=Volt+AlphaP*t(1:Busnum); +Volt(Balance)=balVolt; +UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); +end \ No newline at end of file diff --git a/SolveIt.m b/SolveIt.m new file mode 100644 index 0000000..039d223 --- /dev/null +++ b/SolveIt.m @@ -0,0 +1,46 @@ +function XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,Busnum,Loadi) +LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx); +H=-deltdeltF+ddh; +t1=diag(Init_Z./Init_L-Init_W./Init_U); +t2=-deltG*( t1 )*deltG'; +aa=[ + (H+t2),deltH; + deltH',zeros(length(Init_Y)); + ]; +yy=[LxComa;-Ly]; +%% 平衡节点电压不变 +t=size(Loadi,1)*2; +aa(t+Balance,:)=0; +aa(:,t+Balance)=0; +%aa(t+Balance,t+Balance)=1; +aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); +deltG(t+Balance,:)=0; +%% +t=size(Loadi,1)*2+Busnum*1; +aa(t+Balance,:)=0; +aa(:,t+Balance)=0; +%aa(t+Balance,t+Balance)=1; +aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); +deltG(t+Balance,:)=0; +%% +dxdy=aa\yy; +%% KLU +%spy(aa) +%dxdy = klu(aa,'\',full(yy)); +%% +dX=dxdy(1:ContrlCount); +dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum); +dL=Lz+deltG'*dX; +dU=-Lw-deltG'*dX; +dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL; +dW=-diag(Init_U)\Luu-diag(Init_U)\diag(Init_W)*dU; +XX=[ + dX; + dY; + dZ; + dW; + dL; + dU; + +]; +end \ No newline at end of file diff --git a/func_ddg.m b/func_ddg.m new file mode 100644 index 0000000..7406e6d --- /dev/null +++ b/func_ddg.m @@ -0,0 +1,4 @@ +function ddg=func_ddg() + +ddg=0; +end \ No newline at end of file diff --git a/func_ddh.m b/func_ddh.m new file mode 100644 index 0000000..92f4f86 --- /dev/null +++ b/func_ddh.m @@ -0,0 +1,48 @@ +function ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Angle,Loadi,ContrlCount) +mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); +mat_INV_AngleIJ=mat_AngleIJ'; +yP=Init_Y(1:length(Init_Y)/2);%暂时改这里 20111227 +yQ=Init_Y(length(Init_Y)/2+1:length(Init_Y));%暂时改这里 20111227 +t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum); +t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum)*Y.*cos(mat_AngleIJ); +t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); +t4=-(sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum) -sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum); +ddPdTdT=t3+t4;%ok1 +t1=(-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,yP,Busnum,Busnum); +t2= -sparse(1:Busnum,1:Busnum, sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP ,Busnum,Busnum)*Y.*sin(mat_AngleIJ)+sparse(1:Busnum,1:Busnum,Y.*sin(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum); +ddPdVdT=t1+t2;%ok1 +t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP); +t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt); +t3=-diag(yP)*diag(Y.*sin(mat_AngleIJ)*Volt); +t4=-Y.*sin(mat_INV_AngleIJ)*diag( diag(Volt)*yP ); +ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样 +t1=Y.*cos(mat_INV_AngleIJ)*diag(yP); +t2=diag(yP)*Y.*cos(mat_AngleIJ); +ddPdVdV=t1+t2; +t1=-diag(Y.*sin(mat_AngleIJ)*Volt); +t2=diag(Volt)*Y.*sin(mat_INV_AngleIJ); +t3=(t1+t2)*diag( diag(Volt)*yQ ); +t4=-diag( diag(Volt)*yQ )*Y.*sin(mat_AngleIJ); +t5=diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yQ); +t6=-(t4+t5)*diag(Volt); +ddQdTdT=t3+t6;%ok1 +t1=(diag(Y.*cos(mat_AngleIJ)*Volt)-diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(yQ); +t2=+diag( diag(Volt)*yQ )*Y.*cos(mat_AngleIJ)-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ); +ddQdVdT=t1+t2; +t1=Y.*cos(mat_INV_AngleIJ)*diag(diag(Volt)*yQ); +t2=diag(yQ)*diag(Y.*cos(mat_AngleIJ)*Volt); +t3=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ); +t4=-diag(yQ)*Y.*cos(mat_AngleIJ)*diag(Volt); +ddQdTdV=t1+t2+t3+t4; +t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ); +t2=diag(yQ)*Y.*sin(mat_AngleIJ); +ddQdVdV=t1+t2; +t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; + ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT; +]; +sizeLoadi=size(Loadi,1)*2; +ddh=[ + sparse(sizeLoadi,ContrlCount); + sparse(2*Busnum,sizeLoadi),-t; + ]; +end \ No newline at end of file diff --git a/func_deltF.m b/func_deltF.m new file mode 100644 index 0000000..7bfb3c6 --- /dev/null +++ b/func_deltF.m @@ -0,0 +1,8 @@ +function deltF=func_deltF(SEPD,ContrlCount) +t=sparse(ones(length(SEPD),1)); +deltF=[ + t; + sparse(ContrlCount-length(t),1); +]; + +end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m new file mode 100644 index 0000000..0522de4 --- /dev/null +++ b/func_deltG.m @@ -0,0 +1,22 @@ +function deltG=func_deltG(busNum,Loadi) +dgV_dPD=sparse(length(Loadi),busNum); +dgV_dQD=sparse(length(Loadi),busNum); +dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum); + sparse(busNum,busNum); + ]; +% +dgPD_dPD=sparse(eye(length(Loadi))); +dgPD_dQD=sparse(length(Loadi),length(Loadi)); +dgPD_dx=sparse(busNum*2,length(Loadi)); +% +dgQD_dPD=sparse(length(Loadi),length(Loadi)); +dgQD_dQD=sparse(eye(length(Loadi))); +dgQD_dx=sparse(busNum*2,length(Loadi)); +%% + +deltG=[ + dgPD_dPD,dgQD_dPD,dgV_dPD; + dgPD_dQD,dgQD_dQD,dgV_dQD; + dgPD_dx,dgQD_dx,dgV_dx; +]; +end \ No newline at end of file diff --git a/func_deltH.m b/func_deltH.m new file mode 100644 index 0000000..fa65c0a --- /dev/null +++ b/func_deltH.m @@ -0,0 +1,7 @@ +function deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Angle,Loadi) +dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)]; +dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)]; +dH_dx = jacobian(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 +deltH=[dH_dPD;dH_dQD;dH_dx']; + +end \ No newline at end of file diff --git a/func_deltdeltF.m b/func_deltdeltF.m new file mode 100644 index 0000000..17fdbb2 --- /dev/null +++ b/func_deltdeltF.m @@ -0,0 +1,5 @@ +function deltdeltF=func_deltdeltF(SEPD,ContrlCount) +deltdeltF=[ + 0*zeros(ContrlCount); + ]; +end \ No newline at end of file diff --git a/jacobian.m b/jacobian.m new file mode 100644 index 0000000..45a7307 --- /dev/null +++ b/jacobian.m @@ -0,0 +1,20 @@ +function [Jacob]=jacobian(Busnum,Volt,Y,Angle,UAngel,r,c) +%************************************************************************** +% 程序功能 : 子函数——形成雅可比矩阵Jacobian +% 编 者： +% 编制时间：2010.12 +%************************************************************************** +%%参照图书馆6楼的书编写 +%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q +AngleIJ=UAngel(r)-UAngel(c)-Angle; +mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum); +mat_IvAngleIJ=mat_AngleIJ'; +H=sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); +N=-sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); +J=sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); +L=sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); +t1=[J,L; + H,N; +]'; +Jacob=-t1; +end \ No newline at end of file diff --git a/run.m b/run.m index 40a0f0a..4c671ea 100644 --- a/run.m +++ b/run.m @@ -1,9 +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 ... - cap]=dataRead(lineZ,'D:\share\feeder123\data.txt'); + cap]=dataRead(lineZ,'feeder13\data1.txt'); a=exp(1j*2*pi/3); Tp2f=1/3*[1 1 1; 1 a a^2; @@ -26,6 +26,7 @@ QG=sparse(busNum,1); QGi=[Balance]; PD=Pabc/3; QD=Qabc/3; +Loadi=find(PD~=0); maxD=100000;% 最大不平衡量 EPS=1e-5; k=0; @@ -35,7 +36,6 @@ fsY00=fsY0; fsY22=fsY2; Vf2=sparse(busNum,1); If2=sparse(busNum,1); -% Vf0=sparse(busNum,1); Vf0=sparse(busNum,1); If0=sparse(busNum,1); %准备序矩阵 @@ -115,58 +115,110 @@ ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ... phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ ); 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; -tic -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)); - fprintf('迭代时间%f\n',toc); -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)); - +%% 用牛顿法求解begin +% [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; +% tic +% 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)); +% fprintf('迭代时间%f\n',toc); +% 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)); +%% 用牛顿法求解end +%% 开始进入状态估计 +%状态量 +SEVoltpA=sparse(ones(busNum,1)); +SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi); +SEVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi); +SEphaseASpotLoadP=zeros(length(phaseASpotLoadP),1); +SEphaseBSpotLoadP=zeros(length(phaseBSpotLoadP),1); +SEphaseCSpotLoadP=zeros(length(phaseCSpotLoadP),1); +SEphaseASpotLoadQ=zeros(length(phaseASpotLoadQ),1); +SEphaseBSpotLoadQ=zeros(length(phaseBSpotLoadQ),1); +SEphaseCSpotLoadQ=zeros(length(phaseCSpotLoadQ),1); +SEVmf1=sparse(ones(busNum,1)); +SEVaf1=sparse(zeros(busNum,1)); +SEPD=sparse(length(Loadi),1); +SEQD=sparse(length(Loadi),1); +KK=0; +plotGap=zeros(1,60); +%初始化 +%状态量为 SEPD SEQD SEVmf1 SEVaf1 +RestraintCount=length(SEVmf1)+length(Loadi)*2; +ContrlCount=length(SEVmf1)*2+length(Loadi)*2; +CenterA=0.1; +Init_Z=sparse(ones(RestraintCount,1)); +Init_W=sparse(-1*ones(RestraintCount,1)); +Init_L=1*sparse(ones(RestraintCount,1)); +Init_U=1*sparse(ones(RestraintCount,1)); +Init_Y=sparse(2*busNum,1);%等式约束乘子 +Gap=(Init_L'*Init_Z-Init_U'*Init_W); +Init_u=Gap/2/RestraintCount*CenterA; +deltH=func_deltH(busNum,SEVmf1,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi); +deltG=func_deltG(busNum,Loadi); +L_1Z=diag(Init_Z./Init_L); +U_1W=diag(Init_W./Init_U); +deltdeltF=func_deltdeltF(SEPD,ContrlCount); +ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount); +ddg=func_ddg(); +deltF=func_deltF(SEPD,ContrlCount); +Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); +Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1); +Mat_G=FormG(SEVaf1,SEPD,SEQD); +Mat_H=FormH(busNum,SEVmf1,PG,SEPD,QG,SEQD,fsY1amp,SEVaf1,r,c,fsY1ang); +Ly=Mat_H; +Lz=FormLz(Mat_G,Init_L,busNum,Loadi); +Lw=FormLw(Mat_G,Init_U,busNum,Loadi); +Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); +YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); +%% 开始解方程 +fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); +XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi); +[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum); +[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,SEVmf1,SEVaf1,ContrlCount,Balance,busNum,PD,QD,Loadi); +Gap=(Init_L'*Init_Z-Init_U'*Init_W); \ No newline at end of file