From 9c5239bd1e9620949e04c73a5f226f5fccdfa415 Mon Sep 17 00:00:00 2001 From: facat Date: Thu, 24 May 2012 10:34:04 +0800 Subject: [PATCH] =?UTF-8?q?=E5=B7=B2=E5=85=A8=E9=83=A8=E7=A8=80=E7=96=8F?= =?UTF-8?q?=E5=8C=96?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- FormG.m | 2 +- FormH.m | 21 +-- Modification.m | 4 +- OPF.m | 9 +- OPF_Init.m | 24 ++- func_ddg.m | 2 +- func_ddh3.asv | 375 +------------------------------------- func_ddh3.m | 445 +++++---------------------------------------- func_deltF.m | 8 +- func_deltG.m | 64 ++++--- func_deltH.m | 33 ++-- func_deltdeltF.asv | 21 +-- func_deltdeltF.m | 12 +- jacobian_M3.m | 82 +-------- 14 files changed, 166 insertions(+), 936 deletions(-) diff --git a/FormG.m b/FormG.m index 003f7c4..d8e17e8 100644 --- a/FormG.m +++ b/FormG.m @@ -13,7 +13,7 @@ Mat_G=[ %GP; PG(PGi); QG(PVi); - PD; + sparse(PD); %GQ(PVi); %[0 1.45]'; Volt'; diff --git a/FormH.m b/FormH.m index b43f9b5..f38d2e8 100644 --- a/FormH.m +++ b/FormH.m @@ -1,29 +1,14 @@ -function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND) -t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t2=Volt'*Volt; -t3=t1.*t2; -t4=sum(-t3,2);%P -t5=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t6=t2.*t5; -t7=sum(-t6,2);%Q -t8=PG-PD; -t9=QG-QD; -%Mat_H=([(PG-PD)',(QG-QD)'])'+([t4',t7'])'; -Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum); -Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum); -Mat_H=Mat_H'; +function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND) %% QDcos=textread('1047glys.txt'); QD(QD~=0)=PD(QD~=0)./tan(QDcos); QD(QD_NON_ZERO_IND)=QD_NON_ZERO; %% -%%%%一下是学姐给的公式 -AngleIJ=AngleIJMat-angle(GB); +AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); %dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt'; dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt'; dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; -%Mat_H(1:2:2*Busnum)=dP(1:Busnum);暂时改一下 20111227 -%Mat_H(2:2:2*Busnum)=dQ(1:Busnum);暂时改一下 20111227 + Mat_H=[dP;dQ;]; end \ No newline at end of file diff --git a/Modification.m b/Modification.m index 2cfb824..9dfbfa6 100644 --- a/Modification.m +++ b/Modification.m @@ -1,8 +1,8 @@ function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD) AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); -fprintf('AlphaP %f\n',AlphaP); +fprintf('AlphaP %f\n',full(AlphaP)); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); -fprintf('AlphaD %f\n',AlphaD); +fprintf('AlphaD %f\n',full(AlphaD)); Init_Z=Init_Z+AlphaD*deltZ'; Init_L=Init_L+AlphaP*deltL'; diff --git a/OPF.m b/OPF.m index e9f60e6..d77ff27 100644 --- a/OPF.m +++ b/OPF.m @@ -28,10 +28,11 @@ while(abs(Gap)>Precision) end plotGap(KK+1)=Gap; Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=repmat(UAngel',1,Busnum)-repmat(UAngel,Busnum,1); + %AngleIJMat=repmat(UAngel',1,Busnum)-repmat(UAngel,Busnum,1); + AngleIJMat=0; %% 开始计算OPF %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi); + deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle); %% 形成不等式约束的雅克比 deltG=func_deltG(Busnum,PVi,PGi); %% @@ -40,7 +41,7 @@ while(abs(Gap)>Precision) %% 形成海森阵 deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0); %% 形成ddHy - ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y); + ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle); %% 开始构建ddg ddg=func_ddg(PGi,PVi,Busnum,RestraintCount); %% 开始构建deltF @@ -50,7 +51,7 @@ 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(Volt,PVi,PGi,PG,QG,PD); - Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND); + Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND); Ly=Mat_H; Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0); diff --git a/OPF_Init.m b/OPF_Init.m index f9eb03f..28eaa75 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -1,23 +1,21 @@ function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0) RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数 t_Bal_volt=Volt(Balance); -Volt=ones(1,Busnum); +Volt=sparse(ones(1,Busnum)); Volt(Balance)=t_Bal_volt; %Volt(Balance)=1; -UAngel=zeros(1,Busnum); -Init_Z=ones(1,RestraintCount); -Init_W=-1*ones(1,RestraintCount); -Init_L=ones(1,RestraintCount); -Init_U=ones(1,RestraintCount); -%Init_Y=zeros(1,2*Busnum); -%Init_Y=ones(1,2*Busnum); -Init_Y=zeros(1,2*Busnum);%与学姐一致 +UAngel=sparse(1,Busnum); +Init_Z=sparse(ones(1,RestraintCount)); +Init_W=sparse(-1*ones(1,RestraintCount)); +Init_L=sparse(ones(1,RestraintCount)); +Init_U=sparse(ones(1,RestraintCount)); +Init_Y=sparse(1,2*Busnum);%与学姐一致 %Init_Y(1:2:2*Busnum)=1e-10; %Init_Y(2:2:2*Busnum)=-1e-10; -tPU=GenU(:,2);% 发电机有功上限 -tQU=PVQU(:,1);% 无功上限 -tPL=GenL(:,2);% 发电机有功下限 -tQL=PVQL(:,1);% 无功下限 +tPU=sparse(GenU(:,2));% 发电机有功上限 +tQU=sparse(PVQU(:,1));% 无功上限 +tPL=sparse(GenL(:,2));% 发电机有功下限 +tQL=sparse(PVQL(:,1));% 无功下限 %PG(4:5)=[4.5 4.5]; PG(PGi)=(tPU+tPL)/2; %QG(4:5)=[0 1.45]; diff --git a/func_ddg.m b/func_ddg.m index fe38d77..9e778fa 100644 --- a/func_ddg.m +++ b/func_ddg.m @@ -1,6 +1,6 @@ function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount) -t=zeros(size(PVi,1)+size(PGi,1)+3*Busnum,RestraintCount); +t=sparse(size(PVi,1)+size(PGi,1)+2*Busnum,RestraintCount); ddg=t; end \ No newline at end of file diff --git a/func_ddh3.asv b/func_ddh3.asv index 5704364..a1632f6 100644 --- a/func_ddh3.asv +++ b/func_ddh3.asv @@ -1,368 +1,23 @@ function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y) %决定用循环重写 -ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2; -%deltaPi/deltaThytai_deltaThytaj 非对角元素 -dPidTidTj=zeros(Busnum); -t1=-Volt'*Volt; -t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t3=t1.*t2; -for I=1:Busnum - for J=1:Busnum - dPidTidTj(I,J)=dPidTidTj(I,J)+Init_Y(2*I-1)*t3(I,J); - end -end -t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -t3=t1.*t2; -dPjdTidTj=zeros(Busnum); +ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; -for I=1:Busnum - for J=1:Busnum - dPjdTidTj(I,J)=dPjdTidTj(I,J)+Init_Y(2*J-1)*t3(I,J); - end -end -dPdTidTj=dPidTidTj+dPjdTidTj;%最终非对角元素 - - - - %% deltaP/deltaThyta_deltaThyta 对角元素 -t1=Volt'*Volt; -t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t3=t1.*t2; -t4=t3-diag(diag(t3)); -t5=sum(t4,2); -t6=t5'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -dPidTidTi=diag(t6); -t1=Volt'*Volt; -t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t2=t1.*t2; -dPidTjdTj=zeros(Busnum); -for J=1:Busnum - for I=1:Busnum - - if I==J - continue; - end - dPidTjdTj(J,J)=dPidTjdTj(J,J)+Init_Y(2*I-1)*t2(I,J); - - end -end - -dPdTidTi=dPidTidTi+dPidTjdTj; -%%%%%%%%% -dPdTidTj=dPdTidTj-diag(diag(dPdTidTj)); -hh=dPdTidTj+dPdTidTi; %#ok - %% deltaP/deltaThytai_dVi 对角元素 - t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -t2=diag(t1); -t3=t1-diag(t2);%去掉了对角元素的 -t4=sum(t3,2); -t4=t4'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -dPidTidVi=diag(t4); - -t1=-Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -dPjdTidVi=zeros(Busnum); -for J=1:Busnum - for I=1:Busnum - - if I==J - continue; - end - dPjdTidVi(J,J)=dPjdTidVi(J,J)+Init_Y(2*I-1)*t1(I,J); - - end -end -dPdTidVi=dPidTidVi+dPjdTidVi; -%% deltaP/deltaThytai_dVj 非对角元素 - -t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));%20111224 存疑与学姐给的公式不一致 -%t1=ones(1,Busnum)'*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -dPidTidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dPidTidVj(I,J)=dPidTidVj(I,J)+Init_Y(2*I-1)*t1(I,J); - end -end - -dPjdTidVj=zeros(Busnum); -t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -for I=1:Busnum - for J=1:Busnum - dPjdTidVj(I,J)=dPjdTidVj(I,J)+Init_Y(2*J-1)*t1(I,J); - end -end - -dPdTidVj=dPidTidVj+dPjdTidVj; -%%%%%%%%%%%%%%%%%%%%%% -dPdTidVj=dPdTidVj-diag(diag(dPdTidVj)); -hh=dPdTidVj+dPdTidVi; %#ok -%% deltaP/dVi_deltaThytaj 非对角元素 -t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -dPidVidTj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dPidVidTj(I,J)=dPidVidTj(I,J)+Init_Y(2*I-1)*t1(I,J); - end -end -t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -dPjdVidTj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dPjdVidTj(I,J)=dPjdVidTj(I,J)+Init_Y(2*J-1)*t1(I,J); - end -end -dPdVidTj=dPidVidTj+dPjdVidTj; -%% deltaPi/dVi_deltaThyta 对角元素 - -dPdVidTi=dPdTidVi; -%%%%%%%%%%%%%%%% -dPdVidTj=dPdVidTj-diag(diag(dPdVidTj)); -hh=dPdVidTj+dPdVidTi; -%% deltaP/dVi_dVj 非对角元素 -t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat)); -dPidVidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dPidVidVj(I,J)=dPidVidVj(I,J)+Init_Y(2*I-1)*t1(I,J); - end -end - - -t1=-(real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat)); -dPjdVidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dPjdVidVj(I,J)=dPjdVidVj(I,J)+Init_Y(2*J-1)*t1(I,J); - end -end -dPdVidVj=dPidVidVj+dPjdVidVj; - %% deltaP/dVi_dVi 对角元素 -t0=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t1=diag(t0); -t2=t1'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -dPidVidVi=-2*diag(t2); -dPidVjdVj=0; -dPdVidVi=dPidVidVi+dPidVjdVj; -%%%%%%%%%%% -dPdVidVj=dPdVidVj-diag(diag(dPdVidVj)); -hh=dPdVidVj+dPdVidVi; -%% 生成APi -APi=zeros(2*Busnum,2*Busnum); -dPdTidTj=dPdTidTj-diag(diag(dPdTidTj)); -dPdTidVj=dPdTidVj-diag(diag(dPdTidVj)); -dPdVidTj=dPdVidTj-diag(diag(dPdVidTj)); -dPdVidVj=dPdVidVj-diag(diag(dPdVidVj)); - -APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTj;%%非对角 TT -APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVj;%%非对角 TV -APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTj;%%非对角 VT -APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVj;%%非对角 VV - -APi(1:2:2*Busnum,1:2:2*Busnum)=APi(1:2:2*Busnum,1:2:2*Busnum)+dPdTidTi;%%对角 -APi(1:2:2*Busnum,2:2:2*Busnum)=APi(1:2:2*Busnum,2:2:2*Busnum)+dPdTidVi;%%对角 -APi(2:2:2*Busnum,1:2:2*Busnum)=APi(2:2:2*Busnum,1:2:2*Busnum)+dPdVidTi;%%对角 -APi(2:2:2*Busnum,2:2:2*Busnum)=APi(2:2:2*Busnum,2:2:2*Busnum)+dPdVidVi;%%对角 -%% deltaQ/deltaThyta_deltaThyta 非对角元素 -t1=-Volt'*Volt; -% %t1=Volt'*Volt; -t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t3=t1.*t2; -dQidTidTj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQidTidTj(I,J)=dQidTidTj(I,J)+Init_Y(2*I)*t3(I,J); - end -end - -t2=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat); -t3=-t1.*t2; - -dQjdTidTj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQjdTidTj(I,J)=dQjdTidTj(I,J)+Init_Y(2*J)*t3(I,J); - end -end - -dQdTidTj=dQidTidTj+dQjdTidTj; - %% deltaQ/deltaThyta_deltaThyta 对角元素 -t1=Volt'*Volt; -t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t3=t1.*t2; -t4=t3-diag(diag(t3)); -t5=sum(t4,2); -t6=t5'.*Init_Y(2:2:size(Init_Y,2)); -dQidTidTi=diag(t6); -t1=Volt'*Volt; - t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); - t3=t1.*t2; - dQjdTidTi=zeros(Busnum); - for J=1:Busnum - for I=1:Busnum - - if I==J - continue; - end - dQjdTidTi(J,J)=dQjdTidTi(J,J)+Init_Y(2*I)*t3(I,J); - - end - end - -dQdTidTi=dQjdTidTi+dQidTidTi; -%%%%%%%%%%%% -dQdTidTj=dQdTidTj-diag(diag(dQdTidTj)); -hh=dQdTidTj+dQdTidTi; -%% deltaQ/deltaThyta_deltaV 非对角元素 -t1=-Volt; -t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t3=t1'*ones(1,Busnum).*t2; -dQidTidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQidTidVj(I,J)=dQidTidVj(I,J)+Init_Y(2*I)*t3(I,J); - end -end -t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -t3=Volt'*ones(1,Busnum).*t2; -dQjdTidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - %dQjdTidVj(I,J)=dQidTidVj(I,J)+Init_Y(2*J)*t3(I,J); 20111225 - dQjdTidVj(I,J)=dQjdTidVj(I,J)+Init_Y(2*J)*t3(I,J); - end -end - -dQdTidVj=dQidTidVj+dQjdTidVj; - %% deltaQ/deltaThyta_deltaV 对角元素 - -t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t2=-ones(Busnum,1)*Volt.*t1; -t2=t2-diag(diag(t2)); -t3=sum(t2,2); -t3=t3'.*Init_Y(2:2:size(Init_Y,2)); -dQidTidVi=diag(t3); - - - -t1=Volt'*ones(1,Busnum).*(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat)); -t2=t1-diag(diag(t1)); -dQjdTidVi=zeros(Busnum); -for J=1:Busnum - for I=1:Busnum - - if I==J - continue; - end - dQjdTidVi(J,J)=dQjdTidVi(J,J)+Init_Y(2*I)*t2(I,J); - - end -end - -dQdTidVi=dQidTidVi+dQjdTidVi; -dQdTidVj=dQdTidVj-diag(diag(dQdTidVj)); -hh=dQdTidVj+dQdTidVi; -%% deltaQ/deltaV_deltaV 非对角元素 -t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -t2=-t1; -dQidVidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQidVidVj(I,J)=dQidVidVj(I,J)+Init_Y(2*I)*t2(I,J); - end -end - - t1=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat); -dQjdVidVj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQjdVidVj(I,J)=dQjdVidVj(I,J)+Init_Y(2*J)*t1(I,J); - end -end - -dQdVidVj=dQidVidVj+dQjdVidVj; -%% deltaQ/deltaV_deltaV 对角元素 -t1=-2*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -t2=diag(t1); -t3=t2'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -dQidVidVi=diag(t3); - -dQjdVidVi=0; -dQdVidVi=dQidVidVi+dQjdVidVi; -%%%%%%%%%%%%%% -dQdVidVj=dQdVidVj-diag(diag(dQdVidVj)); -hh=dQdVidVi+dQdVidVj; - %% deltaQ/deltaV_deltaThyta 非对角元素 - -t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -t2=ones(Busnum,1)*Volt.*(t1); -dQidVidTj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQidVidTj(I,J)=dQidVidTj(I,J)+Init_Y(2*I)*t2(I,J); - end -end - - -t1=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -t2=-ones(Busnum,1)*Volt.*(t1); -dQjdVidTj=zeros(Busnum); -for I=1:Busnum - for J=1:Busnum - dQjdVidTj(I,J)=dQjdVidTj(I,J)+Init_Y(2*J)*t2(I,J); - end -end - - dQdVidTj=dQidVidTj+dQjdVidTj; -%% deltaQ/deltaV_deltaThyta 对角元素 -dQdVidTi=dQdTidVi;% @ -dQdVidTj=dQdVidTj-diag(diag(dQdVidTj)); -hhd=dQdVidTi+dQdVidTj; -%% 生成AQi -AQi=zeros(2*Busnum,2*Busnum); -dQdTidTj=dQdTidTj-diag(diag(dQdTidTj)); -dQdTidVj=dQdTidVj-diag(diag(dQdTidVj)); -dQdVidTj=dQdVidTj-diag(diag(dQdVidTj)); -dQdVidVj=dQdVidVj-diag(diag(dQdVidVj)); - -AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTj;%%非对角 TT -AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVj;%%非对角 TV -AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTj;%%非对角 VT -AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVj;%%非对角 VV - -AQi(1:2:2*Busnum,1:2:2*Busnum)=AQi(1:2:2*Busnum,1:2:2*Busnum)+dQdTidTi;%%对角 -AQi(1:2:2*Busnum,2:2:2*Busnum)=AQi(1:2:2*Busnum,2:2:2*Busnum)+dQdTidVi;%%对角 -AQi(2:2:2*Busnum,1:2:2*Busnum)=AQi(2:2:2*Busnum,1:2:2*Busnum)+dQdVidTi;%%对角 -AQi(2:2:2*Busnum,2:2:2*Busnum)=AQi(2:2:2*Busnum,2:2:2*Busnum)+dQdVidVi;%%对角 -%% 生成ddh -t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount); - zeros(2*Busnum,size(PVi,1)+size(PGi,1)),AQi+APi; -]; -ddh=t; -%% 以下是学姐给的公式 AngleIJ=AngleIJMat-angle(GB); -%yP=Init_Y(1:2:size(Init_Y,2)); -yP=Init_Y(1:2:size(Init_Y,2));%暂时改这里 20111227 -%yQ=Init_Y(2:2:size(Init_Y,2)); +yP=Init_Y(1:size(Init_Y,2)/2);%暂时改这里 20111227 +yQ=Init_Y(size(Init_Y,2)/2+1:size(Init_Y,2));%暂时改这里 20111227 t1=-diag(Y.*cos(AngleIJ')*diag(Volt)*yP'); t2=diag(diag(Volt)*yP')*Y.*cos(AngleIJ); t3=(t1+t2)*diag(Volt); t4=-(diag(Y.*cos(AngleIJ)*Volt') -diag(Volt)*Y.*cos(AngleIJ') )*diag(diag(Volt)*yP'); ddPdTdT=t3+t4;%ok1 -tttt=t2*diag(Volt);%ok1 -ttttt=diag(Volt)*Y.*cos(AngleIJ')*diag(diag(Volt)*yP');%ok1 -tttttt=tttt+ttttt;%ok1 -ttttttt=-diag(Y.*cos(AngleIJ)*Volt')*diag(diag(Volt)*yP')+tttt;%ok1 -tttttttt=diag(Volt)*Y.*cos(AngleIJ')*diag(diag(Volt)*yP')+t1*diag(Volt);%ok1 t1=(-diag(Y.*sin(AngleIJ)*Volt')+diag(Volt)*Y.*sin(AngleIJ') )*diag(yP); t2= -diag( diag(Volt)*yP' )*Y.*sin(AngleIJ)+diag(Y.*sin(AngleIJ')*diag(Volt)*yP'); ddPdVdT=t1+t2;%ok1 -tttt=-diag( diag(Volt)*yP' )*Y.*sin(AngleIJ); t1=diag( Y.*sin(AngleIJ')*diag(Volt)*yP'); t2=diag(yP)*Y.*sin(AngleIJ)*diag(Volt); t3=-diag(yP)*diag(Y.*sin(AngleIJ)*Volt'); t4=-Y.*sin(AngleIJ')*diag( diag(Volt)*yP' ); ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样 -tttt=t2; -ttttt=t4; t1=Y.*cos(AngleIJ')*diag(yP); t2=diag(yP)*Y.*cos(AngleIJ); ddPdVdV=t1+t2; @@ -374,14 +29,8 @@ t4=-diag( diag(Volt)*yQ' )*Y.*sin(AngleIJ); t5=diag(Y.*sin(AngleIJ')*diag(Volt)*yQ'); t6=-(t4+t5)*diag(Volt); ddQdTdT=t3+t6;%ok1 -tttt=-(t4)*diag(Volt); -ttttt=t2*diag( diag(Volt)*yQ' ); -tttttt=t1*diag( diag(Volt)*yQ' )+tttt; -ttttttt=-t5*diag(Volt)+t2*diag( diag(Volt)*yQ' ); t1=(diag(Y.*cos(AngleIJ)*Volt')-diag(Volt)*Y.*cos(AngleIJ') )*diag(yQ); t2=+diag( diag(Volt)*yQ' )*Y.*cos(AngleIJ)-diag(Y.*cos(AngleIJ')*diag(Volt)*yQ'); -tttt=diag( diag(Volt)*yQ' )*Y.*cos(AngleIJ); -ttttt=-diag(Volt)*Y.*cos(AngleIJ') *diag(yQ); ddQdVdT=t1+t2; t1=Y.*cos(AngleIJ')*diag(diag(Volt)*yQ'); t2=diag(yQ)*diag(Y.*cos(AngleIJ)*Volt'); @@ -391,19 +40,13 @@ ddQdTdV=t1+t2+t3+t4; t1=Y.*sin(AngleIJ')*diag(yQ); t2=diag(yQ)*Y.*sin(AngleIJ); ddQdVdV=t1+t2; -%%%% -t=zeros(2*Busnum); -% t(1:2:2*Busnum,1:2:2*Busnum)=ddPdTdT+ddQdTdT; -% %t(1:2:2*Busnum,2:2:2*Busnum)=ddPdTdV+ddQdTdV; -% %t(2:2:2*Busnum,1:2:2*Busnum)=ddPdVdT+ddQdVdT; -% t(1:2:2*Busnum,2:2:2*Busnum)=ddPdVdT+ddQdVdT; -% t(2:2:2*Busnum,1:2:2*Busnum)=ddPdTdV+ddQdTdV; -% t(2:2:2*Busnum,2:2:2*Busnum)=ddPdVdV+ddQdVdV;暂时改一下 20111227 -t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; - ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT; +%% %% + +t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ; + ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT; ]; -t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount); - zeros(2*Busnum,size(PVi,1)+size(PGi,1)),-t; +t=[zeros(size(PGi,1)+size(PVi,1)+Busnum,ContrlCount); + zeros(2*Busnum,size(PVi,1)+size(PGi,1)+Busnum),-t; ]; ddh=t; end \ No newline at end of file diff --git a/func_ddh3.m b/func_ddh3.m index 3a71958..c710f1a 100644 --- a/func_ddh3.m +++ b/func_ddh3.m @@ -1,413 +1,58 @@ -function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y) +function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle) %决定用循环重写 ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; -%deltaPi/deltaThytai_deltaThytaj 非对角元素 -% dPidTidTj=zeros(Busnum); -% t1=-Volt'*Volt; -% t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t3=t1.*t2; -% for I=1:Busnum -% for J=1:Busnum -% dPidTidTj(I,J)=dPidTidTj(I,J)+Init_Y(2*I-1)*t3(I,J); -% end -% end -% t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -% t3=t1.*t2; -% dPjdTidTj=zeros(Busnum); -% -% for I=1:Busnum -% for J=1:Busnum -% dPjdTidTj(I,J)=dPjdTidTj(I,J)+Init_Y(2*J-1)*t3(I,J); -% end -% end -% dPdTidTj=dPidTidTj+dPjdTidTj;%最终非对角元素 -% -% -% -% %% deltaP/deltaThyta_deltaThyta 对角元素 -% t1=Volt'*Volt; -% t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t3=t1.*t2; -% t4=t3-diag(diag(t3)); -% t5=sum(t4,2); -% t6=t5'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -% dPidTidTi=diag(t6); -% t1=Volt'*Volt; -% t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t2=t1.*t2; -% dPidTjdTj=zeros(Busnum); -% for J=1:Busnum -% for I=1:Busnum -% -% if I==J -% continue; -% end -% dPidTjdTj(J,J)=dPidTjdTj(J,J)+Init_Y(2*I-1)*t2(I,J); -% -% end -% end -% -% dPdTidTi=dPidTidTi+dPidTjdTj; -% %%%%%%%%% -% dPdTidTj=dPdTidTj-diag(diag(dPdTidTj)); -% hh=dPdTidTj+dPdTidTi; %#ok -% %% deltaP/deltaThytai_dVi 对角元素 -% t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -% t2=diag(t1); -% t3=t1-diag(t2);%去掉了对角元素的 -% t4=sum(t3,2); -% t4=t4'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -% dPidTidVi=diag(t4); -% -% t1=-Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -% dPjdTidVi=zeros(Busnum); -% for J=1:Busnum -% for I=1:Busnum -% -% if I==J -% continue; -% end -% dPjdTidVi(J,J)=dPjdTidVi(J,J)+Init_Y(2*I-1)*t1(I,J); -% -% end -% end -% dPdTidVi=dPidTidVi+dPjdTidVi; -% %% deltaP/deltaThytai_dVj 非对角元素 -% -% t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));%20111224 存疑与学姐给的公式不一致 -% %t1=ones(1,Busnum)'*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -% dPidTidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dPidTidVj(I,J)=dPidTidVj(I,J)+Init_Y(2*I-1)*t1(I,J); -% end -% end -% -% dPjdTidVj=zeros(Busnum); -% t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -% for I=1:Busnum -% for J=1:Busnum -% dPjdTidVj(I,J)=dPjdTidVj(I,J)+Init_Y(2*J-1)*t1(I,J); -% end -% end -% -% dPdTidVj=dPidTidVj+dPjdTidVj; -% %%%%%%%%%%%%%%%%%%%%%% -% dPdTidVj=dPdTidVj-diag(diag(dPdTidVj)); -% hh=dPdTidVj+dPdTidVi; %#ok -% %% deltaP/dVi_deltaThytaj 非对角元素 -% t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -% dPidVidTj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dPidVidTj(I,J)=dPidVidTj(I,J)+Init_Y(2*I-1)*t1(I,J); -% end -% end -% t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat)); -% dPjdVidTj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dPjdVidTj(I,J)=dPjdVidTj(I,J)+Init_Y(2*J-1)*t1(I,J); -% end -% end -% dPdVidTj=dPidVidTj+dPjdVidTj; -% %% deltaPi/dVi_deltaThyta 对角元素 -% -% dPdVidTi=dPdTidVi; -% %%%%%%%%%%%%%%%% -% dPdVidTj=dPdVidTj-diag(diag(dPdVidTj)); -% hh=dPdVidTj+dPdVidTi; -% %% deltaP/dVi_dVj 非对角元素 -% t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat)); -% dPidVidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dPidVidVj(I,J)=dPidVidVj(I,J)+Init_Y(2*I-1)*t1(I,J); -% end -% end -% -% -% t1=-(real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat)); -% dPjdVidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dPjdVidVj(I,J)=dPjdVidVj(I,J)+Init_Y(2*J-1)*t1(I,J); -% end -% end -% dPdVidVj=dPidVidVj+dPjdVidVj; -% %% deltaP/dVi_dVi 对角元素 -% t0=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t1=diag(t0); -% t2=t1'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数 -% dPidVidVi=-2*diag(t2); -% dPidVjdVj=0; -% dPdVidVi=dPidVidVi+dPidVjdVj; -% %%%%%%%%%%% -% dPdVidVj=dPdVidVj-diag(diag(dPdVidVj)); -% hh=dPdVidVj+dPdVidVi; -% %% 生成APi -% APi=zeros(2*Busnum,2*Busnum); -% dPdTidTj=dPdTidTj-diag(diag(dPdTidTj)); -% dPdTidVj=dPdTidVj-diag(diag(dPdTidVj)); -% dPdVidTj=dPdVidTj-diag(diag(dPdVidTj)); -% dPdVidVj=dPdVidVj-diag(diag(dPdVidVj)); -% -% APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTj;%%非对角 TT -% APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVj;%%非对角 TV -% APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTj;%%非对角 VT -% APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVj;%%非对角 VV -% -% APi(1:2:2*Busnum,1:2:2*Busnum)=APi(1:2:2*Busnum,1:2:2*Busnum)+dPdTidTi;%%对角 -% APi(1:2:2*Busnum,2:2:2*Busnum)=APi(1:2:2*Busnum,2:2:2*Busnum)+dPdTidVi;%%对角 -% APi(2:2:2*Busnum,1:2:2*Busnum)=APi(2:2:2*Busnum,1:2:2*Busnum)+dPdVidTi;%%对角 -% APi(2:2:2*Busnum,2:2:2*Busnum)=APi(2:2:2*Busnum,2:2:2*Busnum)+dPdVidVi;%%对角 -% %% deltaQ/deltaThyta_deltaThyta 非对角元素 -% t1=-Volt'*Volt; -% % %t1=Volt'*Volt; -% t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -% t3=t1.*t2; -% dQidTidTj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQidTidTj(I,J)=dQidTidTj(I,J)+Init_Y(2*I)*t3(I,J); -% end -% end -% -% t2=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat); -% t3=-t1.*t2; -% -% dQjdTidTj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQjdTidTj(I,J)=dQjdTidTj(I,J)+Init_Y(2*J)*t3(I,J); -% end -% end -% -% dQdTidTj=dQidTidTj+dQjdTidTj; -% %% deltaQ/deltaThyta_deltaThyta 对角元素 -% t1=Volt'*Volt; -% t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -% t3=t1.*t2; -% t4=t3-diag(diag(t3)); -% t5=sum(t4,2); -% t6=t5'.*Init_Y(2:2:size(Init_Y,2)); -% dQidTidTi=diag(t6); -% t1=Volt'*Volt; -% t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -% t3=t1.*t2; -% dQjdTidTi=zeros(Busnum); -% for J=1:Busnum -% for I=1:Busnum -% -% if I==J -% continue; -% end -% dQjdTidTi(J,J)=dQjdTidTi(J,J)+Init_Y(2*I)*t3(I,J); -% -% end -% end -% -% dQdTidTi=dQjdTidTi+dQidTidTi; -% %%%%%%%%%%%% -% dQdTidTj=dQdTidTj-diag(diag(dQdTidTj)); -% hh=dQdTidTj+dQdTidTi; -% %% deltaQ/deltaThyta_deltaV 非对角元素 -% t1=-Volt; -% t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t3=t1'*ones(1,Busnum).*t2; -% dQidTidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQidTidVj(I,J)=dQidTidVj(I,J)+Init_Y(2*I)*t3(I,J); -% end -% end -% t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -% t3=Volt'*ones(1,Busnum).*t2; -% dQjdTidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% %dQjdTidVj(I,J)=dQidTidVj(I,J)+Init_Y(2*J)*t3(I,J); 20111225 -% dQjdTidVj(I,J)=dQjdTidVj(I,J)+Init_Y(2*J)*t3(I,J); -% end -% end -% -% dQdTidVj=dQidTidVj+dQjdTidVj; -% %% deltaQ/deltaThyta_deltaV 对角元素 -% -% t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t2=-ones(Busnum,1)*Volt.*t1; -% t2=t2-diag(diag(t2)); -% t3=sum(t2,2); -% t3=t3'.*Init_Y(2:2:size(Init_Y,2)); -% dQidTidVi=diag(t3); -% -% -% -% t1=Volt'*ones(1,Busnum).*(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat)); -% t2=t1-diag(diag(t1)); -% dQjdTidVi=zeros(Busnum); -% for J=1:Busnum -% for I=1:Busnum -% -% if I==J -% continue; -% end -% dQjdTidVi(J,J)=dQjdTidVi(J,J)+Init_Y(2*I)*t2(I,J); -% -% end -% end -% -% dQdTidVi=dQidTidVi+dQjdTidVi; -% dQdTidVj=dQdTidVj-diag(diag(dQdTidVj)); -% hh=dQdTidVj+dQdTidVi; -% %% deltaQ/deltaV_deltaV 非对角元素 -% t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat); -% t2=-t1; -% dQidVidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQidVidVj(I,J)=dQidVidVj(I,J)+Init_Y(2*I)*t2(I,J); -% end -% end -% -% t1=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat); -% dQjdVidVj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQjdVidVj(I,J)=dQjdVidVj(I,J)+Init_Y(2*J)*t1(I,J); -% end -% end -% -% dQdVidVj=dQidVidVj+dQjdVidVj; -% %% deltaQ/deltaV_deltaV 对角元素 -% t1=-2*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat)); -% t2=diag(t1); -% t3=t2'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数 -% dQidVidVi=diag(t3); -% -% dQjdVidVi=0; -% dQdVidVi=dQidVidVi+dQjdVidVi; -% %%%%%%%%%%%%%% -% dQdVidVj=dQdVidVj-diag(diag(dQdVidVj)); -% hh=dQdVidVi+dQdVidVj; -% %% deltaQ/deltaV_deltaThyta 非对角元素 -% -% t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); -% t2=ones(Busnum,1)*Volt.*(t1); -% dQidVidTj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQidVidTj(I,J)=dQidVidTj(I,J)+Init_Y(2*I)*t2(I,J); -% end -% end -% -% -% t1=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat); -% t2=-ones(Busnum,1)*Volt.*(t1); -% dQjdVidTj=zeros(Busnum); -% for I=1:Busnum -% for J=1:Busnum -% dQjdVidTj(I,J)=dQjdVidTj(I,J)+Init_Y(2*J)*t2(I,J); -% end -% end -% -% dQdVidTj=dQidVidTj+dQjdVidTj; -% %% deltaQ/deltaV_deltaThyta 对角元素 -% dQdVidTi=dQdTidVi;% @ -% dQdVidTj=dQdVidTj-diag(diag(dQdVidTj)); -% hhd=dQdVidTi+dQdVidTj; -% %% 生成AQi -% AQi=zeros(2*Busnum,2*Busnum); -% dQdTidTj=dQdTidTj-diag(diag(dQdTidTj)); -% dQdTidVj=dQdTidVj-diag(diag(dQdTidVj)); -% dQdVidTj=dQdVidTj-diag(diag(dQdVidTj)); -% dQdVidVj=dQdVidVj-diag(diag(dQdVidVj)); -% -% AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTj;%%非对角 TT -% AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVj;%%非对角 TV -% AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTj;%%非对角 VT -% AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVj;%%非对角 VV -% -% AQi(1:2:2*Busnum,1:2:2*Busnum)=AQi(1:2:2*Busnum,1:2:2*Busnum)+dQdTidTi;%%对角 -% AQi(1:2:2*Busnum,2:2:2*Busnum)=AQi(1:2:2*Busnum,2:2:2*Busnum)+dQdTidVi;%%对角 -% AQi(2:2:2*Busnum,1:2:2*Busnum)=AQi(2:2:2*Busnum,1:2:2*Busnum)+dQdVidTi;%%对角 -% AQi(2:2:2*Busnum,2:2:2*Busnum)=AQi(2:2:2*Busnum,2:2:2*Busnum)+dQdVidVi;%%对角 -% %% 生成ddh -% t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount); -% zeros(2*Busnum,size(PVi,1)+size(PGi,1)),AQi+APi; -% ]; -% ddh=t; -%% 以下是学姐给的公式 -AngleIJ=AngleIJMat-angle(GB); -%yP=Init_Y(1:2:size(Init_Y,2)); -yP=Init_Y(1:size(Init_Y,2)/2);%暂时改这里 20111227 -%yQ=Init_Y(2:2:size(Init_Y,2)); -yQ=Init_Y(size(Init_Y,2)/2+1:size(Init_Y,2));%暂时改这里 20111227 -t1=-diag(Y.*cos(AngleIJ')*diag(Volt)*yP'); -t2=diag(diag(Volt)*yP')*Y.*cos(AngleIJ); -t3=(t1+t2)*diag(Volt); -t4=-(diag(Y.*cos(AngleIJ)*Volt') -diag(Volt)*Y.*cos(AngleIJ') )*diag(diag(Volt)*yP'); -ddPdTdT=t3+t4;%ok1 -% tttt=t2*diag(Volt);%ok1 -% ttttt=diag(Volt)*Y.*cos(AngleIJ')*diag(diag(Volt)*yP');%ok1 -% tttttt=tttt+ttttt;%ok1 -% ttttttt=-diag(Y.*cos(AngleIJ)*Volt')*diag(diag(Volt)*yP')+tttt;%ok1 -% tttttttt=diag(Volt)*Y.*cos(AngleIJ')*diag(diag(Volt)*yP')+t1*diag(Volt);%ok1 -t1=(-diag(Y.*sin(AngleIJ)*Volt')+diag(Volt)*Y.*sin(AngleIJ') )*diag(yP); -t2= -diag( diag(Volt)*yP' )*Y.*sin(AngleIJ)+diag(Y.*sin(AngleIJ')*diag(Volt)*yP'); -ddPdVdT=t1+t2;%ok1 -%tttt=-diag( diag(Volt)*yP' )*Y.*sin(AngleIJ); -t1=diag( Y.*sin(AngleIJ')*diag(Volt)*yP'); -t2=diag(yP)*Y.*sin(AngleIJ)*diag(Volt); -t3=-diag(yP)*diag(Y.*sin(AngleIJ)*Volt'); -t4=-Y.*sin(AngleIJ')*diag( diag(Volt)*yP' ); -ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样 -% tttt=t2; -% ttttt=t4; -t1=Y.*cos(AngleIJ')*diag(yP); -t2=diag(yP)*Y.*cos(AngleIJ); -ddPdVdV=t1+t2; -t1=-diag(Y.*sin(AngleIJ)*Volt'); -t2=diag(Volt)*Y.*sin(AngleIJ'); -t3=(t1+t2)*diag( diag(Volt)*yQ' ); -t4=-diag( diag(Volt)*yQ' )*Y.*sin(AngleIJ); -t5=diag(Y.*sin(AngleIJ')*diag(Volt)*yQ'); +%AngleIJ=AngleIJMat-angle(GB); +mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); +mat_INV_AngleIJ=mat_AngleIJ'; +yP=Init_Y(1:size(Init_Y,2)/2);%暂时改这里 20111227 +yQ=Init_Y(size(Init_Y,2)/2+1:size(Init_Y,2));%暂时改这里 20111227 +t1=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yP'); +t2=diag(diag(Volt)*yP')*Y.*cos(mat_AngleIJ); +t3=(t1+t2)*diag(Volt); +t4=-(diag(Y.*cos(mat_AngleIJ)*Volt') -diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(diag(Volt)*yP'); +ddPdTdT=t3+t4;%ok1 +t1=(-diag(Y.*sin(mat_AngleIJ)*Volt')+diag(Volt)*Y.*sin(mat_INV_AngleIJ) )*diag(yP); +t2= -diag( diag(Volt)*yP' )*Y.*sin(mat_AngleIJ)+diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP'); +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 -% tttt=-(t4)*diag(Volt); -% ttttt=t2*diag( diag(Volt)*yQ' ); -% tttttt=t1*diag( diag(Volt)*yQ' )+tttt; -% ttttttt=-t5*diag(Volt)+t2*diag( diag(Volt)*yQ' ); -t1=(diag(Y.*cos(AngleIJ)*Volt')-diag(Volt)*Y.*cos(AngleIJ') )*diag(yQ); -t2=+diag( diag(Volt)*yQ' )*Y.*cos(AngleIJ)-diag(Y.*cos(AngleIJ')*diag(Volt)*yQ'); -% tttt=diag( diag(Volt)*yQ' )*Y.*cos(AngleIJ); -% ttttt=-diag(Volt)*Y.*cos(AngleIJ') *diag(yQ); +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(AngleIJ')*diag(diag(Volt)*yQ'); -t2=diag(yQ)*diag(Y.*cos(AngleIJ)*Volt'); -t3=-diag(Y.*cos(AngleIJ')*diag(Volt)*yQ'); -t4=-diag(yQ)*Y.*cos(AngleIJ)*diag(Volt); +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(AngleIJ')*diag(yQ); -t2=diag(yQ)*Y.*sin(AngleIJ); +t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ); +t2=diag(yQ)*Y.*sin(mat_AngleIJ); ddQdVdV=t1+t2; -%%%% -%t=zeros(2*Busnum); -% t(1:2:2*Busnum,1:2:2*Busnum)=ddPdTdT+ddQdTdT; -% %t(1:2:2*Busnum,2:2:2*Busnum)=ddPdTdV+ddQdTdV; -% %t(2:2:2*Busnum,1:2:2*Busnum)=ddPdVdT+ddQdVdT; -% t(1:2:2*Busnum,2:2:2*Busnum)=ddPdVdT+ddQdVdT; -% t(2:2:2*Busnum,1:2:2*Busnum)=ddPdTdV+ddQdTdV; -% t(2:2:2*Busnum,2:2:2*Busnum)=ddPdVdV+ddQdVdV;暂时改一下 20111227 -% t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT; -% ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT; -% ];再改20111227 t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ; ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT; ]; -t=[zeros(size(PGi,1)+size(PVi,1)+Busnum,ContrlCount); - zeros(2*Busnum,size(PVi,1)+size(PGi,1)+Busnum),-t; -]; +sizePGi=size(PGi,1); +sizePVi=size(PVi,1); +% t=[zeros(size(PGi,1)+size(PVi,1),ContrlCount); +% zeros(2*Busnum,size(PVi,1)+size(PGi,1)),-t; +% ]; +t=[ + sparse(sizePGi+sizePVi+Busnum,ContrlCount); + sparse(2*Busnum,sizePVi+sizePGi+Busnum),-t; + ]; ddh=t; end \ No newline at end of file diff --git a/func_deltF.m b/func_deltF.m index ca22864..1de0334 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -10,10 +10,10 @@ function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum) t1=2*wG.*(PG(PGi) - PG0(PGi) ); t2=2*wD.*(PD-PD0); deltF=[ - t1; - zeros(size(PVi)); - t2; - zeros(2*Busnum,1); + sparse(t1); + sparse(size(PVi,1),size(PVi,2)); + sparse(t2); + sparse(2*Busnum,1); ]; end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m index 23c3719..56ce9e4 100644 --- a/func_deltG.m +++ b/func_deltG.m @@ -1,27 +1,49 @@ function deltG=func_deltG(Busnum,PVi,PGi) -dg1_dPg=eye(size(PGi,1)); -dg2_dPg=zeros(size(PGi,1),size(PVi,1)); -dg3_dPg=zeros(size(PGi,1),Busnum); -dg4_dPg=zeros(size(PGi,1),Busnum); -%% -dg1_dQr=zeros(size(PVi,1),size(PGi,1)); -dg2_dQr=eye(size(PVi,1)); -dg3_dQr=zeros(size(PVi,1),Busnum); -dg4_dQr=zeros(size(PVi,1),Busnum); %% -dg1_dPD=zeros(Busnum,size(PGi,1)); -dg2_dPD=zeros(Busnum,size(PVi,1)); -dg3_dPD=eye(Busnum,Busnum); -dg4_dPD=zeros(Busnum,Busnum); +sizePGi=size(PGi,1); +sizePVi=size(PVi,1); +%% +% dg1_dPg=eye(size(PGi,1)); +% dg2_dPg=zeros(size(PGi,1),size(PVi,1)); +% dg3_dPg=zeros(size(PGi,1),Busnum); +% dg4_dPg=zeros(size(PGi,1),Busnum); +dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi); +dg2_dPg=sparse(sizePGi,sizePVi); +dg3_dPg=sparse(sizePGi,Busnum); +dg4_dPg=sparse(sizePGi,Busnum); %% -dg1_dx=zeros(2*Busnum,size(PGi,1)); -dg2_dx=zeros(2*Busnum,size(PVi,1)); -dg3_dx=zeros(2*Busnum,Busnum); -dg4_dx=zeros(2*Busnum,Busnum); -for I=1:Busnum - %dg3_dx(2*I,I)=1;暂时改一下 - dg4_dx(I,I)=1; -end +% dg1_dQr=zeros(size(PVi,1),size(PGi,1)); +% dg2_dQr=eye(size(PVi,1)); +% dg3_dQr=zeros(size(PVi,1),Busnum); +% dg4_dQr=zeros(size(PVi,1),Busnum); +dg1_dQr=sparse(sizePVi,sizePGi); +dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi); +dg3_dQr=sparse(sizePVi,Busnum); +dg4_dQr=sparse(sizePVi,Busnum); +%% +% dg1_dPD=zeros(Busnum,size(PGi,1)); +% dg2_dPD=zeros(Busnum,size(PVi,1)); +% dg3_dPD=eye(Busnum,Busnum); +% dg4_dPD=zeros(Busnum,Busnum); +dg1_dPD=sparse(Busnum,size(PGi,1)); +dg2_dPD=sparse(Busnum,size(PVi,1)); +dg3_dPD=sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); +dg4_dPD=sparse(Busnum,Busnum); +%% +% dg1_dx=zeros(2*Busnum,size(PGi,1)); +% dg2_dx=zeros(2*Busnum,size(PVi,1)); +% dg3_dx=zeros(2*Busnum,Busnum); +% dg4_dx=zeros(2*Busnum,Busnum); +% for I=1:Busnum +% %dg3_dx(2*I,I)=1;暂时改一下 +% dg4_dx(I,I)=1; +% end +dg1_dx=sparse(2*Busnum,sizePGi); +dg2_dx=sparse(2*Busnum,sizePVi); +dg3_dx=sparse(2*Busnum,Busnum); +dg4_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum); + sparse(Busnum,Busnum); + ]; %% deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg; dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr; diff --git a/func_deltH.m b/func_deltH.m index ced053c..34f7ab8 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -1,21 +1,24 @@ -function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi) -dH_dPg=zeros(size(PGi,1),2*Busnum); +function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle) +% dH_dPg=zeros(size(PGi,1),2*Busnum); +% +% for I=1:size(PGi,1) +% %dH_dPg(I,2*PVi(I)-1)=-1;王锡凡书上的公式 +% %dH_dPg(I,2*PGi(I)-1)=1;暂时改一下20111227 +% dH_dPg(I,PGi(I))=1; +% end -for I=1:size(PGi,1) - %dH_dPg(I,2*PVi(I)-1)=-1;王锡凡书上的公式 - %dH_dPg(I,2*PGi(I)-1)=1;暂时改一下20111227 - dH_dPg(I,PGi(I))=1; -end +dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum); -dH_dQr=zeros(size(PVi,1),2*Busnum); -for I=1:size(PVi,1) - %dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式 - %dH_dQr(I,2*PVi(I))=1;暂时改一下20111227 - dH_dQr(I,PVi(I)+Busnum)=1; -end +% dH_dQr=zeros(size(PVi,1),2*Busnum); +% for I=1:size(PVi,1) +% %dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式 +% %dH_dQr(I,2*PVi(I))=1;暂时改一下20111227 +% dH_dQr(I,PVi(I)+Busnum)=1; +% end +dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum); dH_dPD=[-eye(Busnum) zeros(Busnum)]; -Angle=angle(GB); -dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat); %形成雅克比矩阵 +%Angle=angle(GB); +dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c); %形成雅克比矩阵 %deltH=[dH_dPg;dH_dQr;dH_dx'];%dH_dx 需要使用一下转置 暂时改一下 deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx']; end \ No newline at end of file diff --git a/func_deltdeltF.asv b/func_deltdeltF.asv index 10ebf6b..3a69884 100644 --- a/func_deltdeltF.asv +++ b/func_deltdeltF.asv @@ -1,16 +1,9 @@ -function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi) -%t1=deltG*[L_1Z-U_1W]*deltG'; -%TotalDim=2*Busnum+2*size(PVi,1); -% deltdeltF=[diag(GenC(:,2))*2,zeros(size(GenC,1),TotalDim-size(GenC,1)); -% zeros(TotalDim-size(GenC,1),TotalDim) +function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0) + +ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数 +C=[wG' zeros(size(PVi))' wD']; +% deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2)); +% zeros(ContrlCount-size(C,2),ContrlCount); % ]; -% deltdeltF=[diag(GenC(:,2))*2,zeros(2,12); -% zeros(12,14); -% ]; -ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2; %P,Q,Volt theta这些控制变量数 -t=zeros(size(PGi,1)); -t(PGi,PGi)=GenC(:,2); -deltdeltF=[t*2,zeros()); - zeros(ContrlCount-size(GenC,1),ContrlCount); - ]; + end \ No newline at end of file diff --git a/func_deltdeltF.m b/func_deltdeltF.m index 4661b57..163883d 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -2,8 +2,14 @@ function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0) ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数 C=[wG' zeros(size(PVi))' wD']; -deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2)); - zeros(ContrlCount-size(C,2),ContrlCount); +sizeC=size(C,2); +diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC); +% deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2)); +% zeros(ContrlCount-size(C,2),ContrlCount); +% ]; +%sizeGenC=size(GenC(:,2),1); +deltdeltF=[ + diagC*2,sparse(sizeC,ContrlCount-sizeC); + sparse(ContrlCount-sizeC,ContrlCount); ]; - end \ No newline at end of file diff --git a/jacobian_M3.m b/jacobian_M3.m index 1f3a9a9..013275a 100644 --- a/jacobian_M3.m +++ b/jacobian_M3.m @@ -1,4 +1,4 @@ -function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat) +function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c) %************************************************************************** % 程序功能 : 子函数——形成雅可比矩阵Jacobian % 编 者： @@ -6,79 +6,13 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat) %************************************************************************** %%参照图书馆6楼的书编写 %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q -AngleIJ=AngleIJMat-Angle; -% t1=Volt'*Volt; -% H=t1.*Y.*sin(AngleIJ); -% N=-t1.*Y.*cos(AngleIJ); -% %J=Volt'*ones(1,Busnum).*cos(AngleIJ);这里错了 -% J=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y; -% -% %L=Volt'*ones(1,Busnum).*sin(AngleIJ);这里错了 -% L=Volt'*ones(1,Busnum).*sin(AngleIJ).*Y; -% -% %%对角 -% t1=Volt'*Volt; -% t2=t1.*Y.*sin(AngleIJ); -% t3=diag(t2); -% t4=t2-diag(t3); -% t5=sum(t4,2); -% HH=-diag(t5); -% t2=t1.*Y.*cos(AngleIJ); -% t3=diag(t2); -% t4=t2-diag(t3); -% t5=sum(t4,2); -% NN=diag(t5); -% %t1=ones(Busnum,1)*Volt; -% t1=ones(Busnum,1)*Volt.*Y; -% t2=t1.*cos(AngleIJ); -% t3=sum(t2,2); -% JJ=diag(t3); -% t1=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y; -% %t1=Volt'*ones(1,Busnum).*cos(AngleIJ); -% t2=diag(t1);% -% JJ=JJ+diag(t2); -% t1=ones(Busnum,1)*Volt.*Y; -% %t1=ones(Busnum,1)*Volt; -% t2=t1.*sin(AngleIJ); -% t3=sum(t2,2); -% LL=diag(t3); -% %t1=Volt'*ones(1,Busnum).*sin(AngleIJ); -% t1=Volt'*ones(1,Busnum).*sin(AngleIJ).*Y; -% t2=diag(t1);% -% %LL=LL-diag(t2);这里错了 -% LL=LL+diag(t2); -% -% H=H-diag(diag(H)); -% N=N-diag(diag(N)); -% J=J-diag(diag(J)); -% L=L-diag(diag(L)); -% -% H=H+HH; -% J=J+JJ; -% N=N+NN; -% L=L+LL; -% -% t1=zeros(2*Busnum); -% t1(1:2:2*Busnum,1:2:2*Busnum)=H; -% t1(1:2:2*Busnum,2:2:2*Busnum)=N; -% t1(2:2:2*Busnum,1:2:2*Busnum)=J; -% t1(2:2:2*Busnum,2:2:2*Busnum)=L; -% Jacob=-t1; -%%以下是学姐给的公式 -H=diag(Volt)*Y.*sin(AngleIJ')*diag(Volt)-diag(Y.*sin(AngleIJ)*Volt')*diag(Volt); -N=-diag(Volt)*Y.*cos(AngleIJ')*diag(Volt)+diag(Y.*cos(AngleIJ)*Volt')*diag(Volt); -J=diag(Y.*cos(AngleIJ)*Volt')+Y.*cos(AngleIJ')*diag(Volt); -L=diag(Y.*sin(AngleIJ)*Volt')+Y.*sin(AngleIJ')*diag(Volt); -H=H; -N=N; -J=J; -L=L; -t1=zeros(2*Busnum); -% t1(1:2:2*Busnum,1:2:2*Busnum)=-H;%10111227 -% t1(1:2:2*Busnum,2:2:2*Busnum)=N; -% t1(2:2:2*Busnum,1:2:2*Busnum)=-J;%10111227 -% t1(2:2:2*Busnum,2:2:2*Busnum)=L; -%暂时改一下 +AngleIJ=UAngel(r)-UAngel(c)-Angle'; +mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum); +mat_IvAngleIJ=mat_AngleIJ'; +H=diag(Volt)*Y.*sin(mat_IvAngleIJ)*diag(Volt)-diag(Y.*sin(mat_AngleIJ)*Volt')*diag(Volt); +N=-diag(Volt)*Y.*cos(mat_IvAngleIJ)*diag(Volt)+diag(Y.*cos(mat_AngleIJ)*Volt')*diag(Volt); +J=diag(Y.*cos(mat_AngleIJ)*Volt')+Y.*cos(mat_IvAngleIJ)*diag(Volt); +L=diag(Y.*sin(mat_AngleIJ)*Volt')+Y.*sin(mat_IvAngleIJ)*diag(Volt); t1=[J,L; H,N; ]';