facat
/
pes2014
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

已全部稀疏化 

Signed-off-by: facat <dmy@dmy-PC.(none)>
This commit is contained in:
facat 2012-05-24 10:34:04 +08:00
parent ca24c4f8d5
commit 9c5239bd1e
14 changed files with 166 additions and 936 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
FormG.m
View File

@ -13,7 +13,7 @@ Mat_G=[
%GP; %GP;
PG(PGi); PG(PGi);
QG(PVi); QG(PVi);
PD; sparse(PD);
%GQ(PVi); %GQ(PVi);
%[0 1.45]'; %[0 1.45]';
Volt'; Volt';

21
FormH.m
View File

@ -1,29 +1,14 @@
function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND) function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,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';
%% %%
QDcos=textread('1047glys.txt'); QDcos=textread('1047glys.txt');
QD(QD~=0)=PD(QD~=0)./tan(QDcos); QD(QD~=0)=PD(QD~=0)./tan(QDcos);
QD(QD_NON_ZERO_IND)=QD_NON_ZERO; QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
%% %%
%%%%一下是学姐给的公式 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
AngleIJ=AngleIJMat-angle(GB);
%dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt'; %dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
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'; 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;]; Mat_H=[dP;dQ;];
end end

4
Modification.m
View File

@ -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) 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); 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); 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_Z=Init_Z+AlphaD*deltZ';
Init_L=Init_L+AlphaP*deltL'; Init_L=Init_L+AlphaP*deltL';

9
OPF.m
View File

@ -28,10 +28,11 @@ while(abs(Gap)>Precision)
end end
plotGap(KK+1)=Gap; plotGap(KK+1)=Gap;
Init_u=Gap/2/RestraintCount*CenterA; 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 %% 开始计算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); deltG=func_deltG(Busnum,PVi,PGi);
%% %%
@ -40,7 +41,7 @@ while(abs(Gap)>Precision)
%% 形成海森阵 %% 形成海森阵
deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0); deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
%% 形成ddHy %% 形成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
ddg=func_ddg(PGi,PVi,Busnum,RestraintCount); ddg=func_ddg(PGi,PVi,Busnum,RestraintCount);
%% 开始构建deltF %% 开始构建deltF
@ -50,7 +51,7 @@ while(abs(Gap)>Precision)
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-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_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; Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0); Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);

24
OPF_Init.m
View File

@ -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) 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; %Ô¼ÊøÌõ¼þÊý RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %Ô¼ÊøÌõ¼þÊý
t_Bal_volt=Volt(Balance); t_Bal_volt=Volt(Balance);
Volt=ones(1,Busnum); Volt=sparse(ones(1,Busnum));
Volt(Balance)=t_Bal_volt; Volt(Balance)=t_Bal_volt;
%Volt(Balance)=1; %Volt(Balance)=1;
UAngel=zeros(1,Busnum); UAngel=sparse(1,Busnum);
Init_Z=ones(1,RestraintCount); Init_Z=sparse(ones(1,RestraintCount));
Init_W=-1*ones(1,RestraintCount); Init_W=sparse(-1*ones(1,RestraintCount));
Init_L=ones(1,RestraintCount); Init_L=sparse(ones(1,RestraintCount));
Init_U=ones(1,RestraintCount); Init_U=sparse(ones(1,RestraintCount));
%Init_Y=zeros(1,2*Busnum); Init_Y=sparse(1,2*Busnum);%与学姐一致
%Init_Y=ones(1,2*Busnum);
Init_Y=zeros(1,2*Busnum);%与学姐一致
%Init_Y(1:2:2*Busnum)=1e-10; %Init_Y(1:2:2*Busnum)=1e-10;
%Init_Y(2:2:2*Busnum)=-1e-10; %Init_Y(2:2:2*Busnum)=-1e-10;
tPU=GenU(:,2);% 发电机有功上限 tPU=sparse(GenU(:,2));% 发电机有功上限
tQU=PVQU(:,1);% 无功上限 tQU=sparse(PVQU(:,1));% 无功上限
tPL=GenL(:,2);% 发电机有功下限 tPL=sparse(GenL(:,2));% 发电机有功下限
tQL=PVQL(:,1);% 无功下限 tQL=sparse(PVQL(:,1));% 无功下限
%PG(4:5)=[4.5 4.5]; %PG(4:5)=[4.5 4.5];
PG(PGi)=(tPU+tPL)/2; PG(PGi)=(tPU+tPL)/2;
%QG(4:5)=[0 1.45]; %QG(4:5)=[0 1.45];

2
func_ddg.m
View File

@ -1,6 +1,6 @@
function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount) 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; ddg=t;
end end

375
func_ddh3.asv
View File

@ -1,368 +1,23 @@
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)
%决定用循环重写 %决定用循环重写
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2; 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<NASGU>
%% 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<NASGU>
%% 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); AngleIJ=AngleIJMat-angle(GB);
%yP=Init_Y(1:2:size(Init_Y,2)); yP=Init_Y(1:size(Init_Y,2)/2);%暂时改这里 20111227
yP=Init_Y(1:2:size(Init_Y,2));%暂时改这里 20111227 yQ=Init_Y(size(Init_Y,2)/2+1:size(Init_Y,2));%暂时改这里 20111227
%yQ=Init_Y(2:2:size(Init_Y,2));
t1=-diag(Y.*cos(AngleIJ')*diag(Volt)*yP'); t1=-diag(Y.*cos(AngleIJ')*diag(Volt)*yP');
t2=diag(diag(Volt)*yP')*Y.*cos(AngleIJ); t2=diag(diag(Volt)*yP')*Y.*cos(AngleIJ);
t3=(t1+t2)*diag(Volt); t3=(t1+t2)*diag(Volt);
t4=-(diag(Y.*cos(AngleIJ)*Volt') -diag(Volt)*Y.*cos(AngleIJ') )*diag(diag(Volt)*yP'); t4=-(diag(Y.*cos(AngleIJ)*Volt') -diag(Volt)*Y.*cos(AngleIJ') )*diag(diag(Volt)*yP');
ddPdTdT=t3+t4;%ok1 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); 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'); t2= -diag( diag(Volt)*yP' )*Y.*sin(AngleIJ)+diag(Y.*sin(AngleIJ')*diag(Volt)*yP');
ddPdVdT=t1+t2;%ok1 ddPdVdT=t1+t2;%ok1
tttt=-diag( diag(Volt)*yP' )*Y.*sin(AngleIJ);
t1=diag( Y.*sin(AngleIJ')*diag(Volt)*yP'); t1=diag( Y.*sin(AngleIJ')*diag(Volt)*yP');
t2=diag(yP)*Y.*sin(AngleIJ)*diag(Volt); t2=diag(yP)*Y.*sin(AngleIJ)*diag(Volt);
t3=-diag(yP)*diag(Y.*sin(AngleIJ)*Volt'); t3=-diag(yP)*diag(Y.*sin(AngleIJ)*Volt');
t4=-Y.*sin(AngleIJ')*diag( diag(Volt)*yP' ); t4=-Y.*sin(AngleIJ')*diag( diag(Volt)*yP' );
ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样 ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样
tttt=t2;
ttttt=t4;
t1=Y.*cos(AngleIJ')*diag(yP); t1=Y.*cos(AngleIJ')*diag(yP);
t2=diag(yP)*Y.*cos(AngleIJ); t2=diag(yP)*Y.*cos(AngleIJ);
ddPdVdV=t1+t2; ddPdVdV=t1+t2;
@ -374,14 +29,8 @@ t4=-diag( diag(Volt)*yQ' )*Y.*sin(AngleIJ);
t5=diag(Y.*sin(AngleIJ')*diag(Volt)*yQ'); t5=diag(Y.*sin(AngleIJ')*diag(Volt)*yQ');
t6=-(t4+t5)*diag(Volt); t6=-(t4+t5)*diag(Volt);
ddQdTdT=t3+t6;%ok1 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); 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'); 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; ddQdVdT=t1+t2;
t1=Y.*cos(AngleIJ')*diag(diag(Volt)*yQ'); t1=Y.*cos(AngleIJ')*diag(diag(Volt)*yQ');
t2=diag(yQ)*diag(Y.*cos(AngleIJ)*Volt'); t2=diag(yQ)*diag(Y.*cos(AngleIJ)*Volt');
@ -391,19 +40,13 @@ ddQdTdV=t1+t2+t3+t4;
t1=Y.*sin(AngleIJ')*diag(yQ); t1=Y.*sin(AngleIJ')*diag(yQ);
t2=diag(yQ)*Y.*sin(AngleIJ); t2=diag(yQ)*Y.*sin(AngleIJ);
ddQdVdV=t1+t2; ddQdVdV=t1+t2;
%%%% %% %%
t=zeros(2*Busnum);
% t(1:2:2*Busnum,1:2:2*Busnum)=ddPdTdT+ddQdTdT; t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ;
% %t(1:2:2*Busnum,2:2:2*Busnum)=ddPdTdV+ddQdTdV; ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT;
% %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=[zeros(size(PGi,1)+size(PVi,1),ContrlCount); t=[zeros(size(PGi,1)+size(PVi,1)+Busnum,ContrlCount);
zeros(2*Busnum,size(PVi,1)+size(PGi,1)),-t; zeros(2*Busnum,size(PVi,1)+size(PGi,1)+Busnum),-t;
]; ];
ddh=t; ddh=t;
end end

445
func_ddh3.m
View File

@ -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; 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<NASGU>
% %% 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<NASGU>
% %% 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); t6=-(t4+t5)*diag(Volt);
ddQdTdT=t3+t6;%ok1 ddQdTdT=t3+t6;%ok1
% tttt=-(t4)*diag(Volt); t1=(diag(Y.*cos(mat_AngleIJ)*Volt')-diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(yQ);
% ttttt=t2*diag( diag(Volt)*yQ' ); t2=+diag( diag(Volt)*yQ' )*Y.*cos(mat_AngleIJ)-diag(Y.*cos(mat_INV_AngleIJ)*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; ddQdVdT=t1+t2;
t1=Y.*cos(AngleIJ')*diag(diag(Volt)*yQ'); t1=Y.*cos(mat_INV_AngleIJ)*diag(diag(Volt)*yQ');
t2=diag(yQ)*diag(Y.*cos(AngleIJ)*Volt'); t2=diag(yQ)*diag(Y.*cos(mat_AngleIJ)*Volt');
t3=-diag(Y.*cos(AngleIJ')*diag(Volt)*yQ'); t3=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ');
t4=-diag(yQ)*Y.*cos(AngleIJ)*diag(Volt); t4=-diag(yQ)*Y.*cos(mat_AngleIJ)*diag(Volt);
ddQdTdV=t1+t2+t3+t4; ddQdTdV=t1+t2+t3+t4;
t1=Y.*sin(AngleIJ')*diag(yQ); t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ);
t2=diag(yQ)*Y.*sin(AngleIJ); t2=diag(yQ)*Y.*sin(mat_AngleIJ);
ddQdVdV=t1+t2; 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 ; t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ;
ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT; ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT;
]; ];
t=[zeros(size(PGi,1)+size(PVi,1)+Busnum,ContrlCount); sizePGi=size(PGi,1);
zeros(2*Busnum,size(PVi,1)+size(PGi,1)+Busnum),-t; 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; ddh=t;
end end

8
func_deltF.m
View File

@ -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) ); t1=2*wG.*(PG(PGi) - PG0(PGi) );
t2=2*wD.*(PD-PD0); t2=2*wD.*(PD-PD0);
deltF=[ deltF=[
t1; sparse(t1);
zeros(size(PVi)); sparse(size(PVi,1),size(PVi,2));
t2; sparse(t2);
zeros(2*Busnum,1); sparse(2*Busnum,1);
]; ];
end end

64
func_deltG.m
View File

@ -1,27 +1,49 @@
function deltG=func_deltG(Busnum,PVi,PGi) 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)); sizePGi=size(PGi,1);
dg2_dPD=zeros(Busnum,size(PVi,1)); sizePVi=size(PVi,1);
dg3_dPD=eye(Busnum,Busnum); %%
dg4_dPD=zeros(Busnum,Busnum); % 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)); % dg1_dQr=zeros(size(PVi,1),size(PGi,1));
dg2_dx=zeros(2*Busnum,size(PVi,1)); % dg2_dQr=eye(size(PVi,1));
dg3_dx=zeros(2*Busnum,Busnum); % dg3_dQr=zeros(size(PVi,1),Busnum);
dg4_dx=zeros(2*Busnum,Busnum); % dg4_dQr=zeros(size(PVi,1),Busnum);
for I=1:Busnum dg1_dQr=sparse(sizePVi,sizePGi);
%dg3_dx(2*I,I)=1;ÔÝĘą¸ÄҝϠdg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
dg4_dx(I,I)=1; dg3_dQr=sparse(sizePVi,Busnum);
end 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; deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg;
dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr; dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr;

33
func_deltH.m
View File

@ -1,21 +1,24 @@
function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi) function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle)
dH_dPg=zeros(size(PGi,1),2*Busnum); % 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=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
%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_dQr=zeros(size(PVi,1),2*Busnum); % dH_dQr=zeros(size(PVi,1),2*Busnum);
for I=1:size(PVi,1) % for I=1:size(PVi,1)
%dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式 % %dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式
%dH_dQr(I,2*PVi(I))=1;暂时改一下20111227 % %dH_dQr(I,2*PVi(I))=1;暂时改一下20111227
dH_dQr(I,PVi(I)+Busnum)=1; % dH_dQr(I,PVi(I)+Busnum)=1;
end % 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)]; dH_dPD=[-eye(Busnum) zeros(Busnum)];
Angle=angle(GB); %Angle=angle(GB);
dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat); %形成雅克比矩阵 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_dx'];%dH_dx 需要使用一下转置 暂时改一下
deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx']; deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
end end

21
func_deltdeltF.asv
View File

@ -1,16 +1,9 @@
function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi) function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0)
%t1=deltG*[L_1Z-U_1W]*deltG';
%TotalDim=2*Busnum+2*size(PVi,1); ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数
% deltdeltF=[diag(GenC(:,2))*2,zeros(size(GenC,1),TotalDim-size(GenC,1)); C=[wG' zeros(size(PVi))' wD'];
% zeros(TotalDim-size(GenC,1),TotalDim) % 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 end

12
func_deltdeltF.m
View File

@ -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这些控制变量数 ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数
C=[wG' zeros(size(PVi))' wD']; C=[wG' zeros(size(PVi))' wD'];
deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2)); sizeC=size(C,2);
zeros(ContrlCount-size(C,2),ContrlCount); 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 end

82
jacobian_M3.m
View File

@ -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 % 程序功能 : 子函数——形成雅可比矩阵Jacobian
% 编 者： % 编 者：
@ -6,79 +6,13 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat)
%************************************************************************** %**************************************************************************
%%参照图书馆6楼的书编写 %%参照图书馆6楼的书编写
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
AngleIJ=AngleIJMat-Angle; AngleIJ=UAngel(r)-UAngel(c)-Angle';
% t1=Volt'*Volt; mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum);
% H=t1.*Y.*sin(AngleIJ); mat_IvAngleIJ=mat_AngleIJ';
% N=-t1.*Y.*cos(AngleIJ); H=diag(Volt)*Y.*sin(mat_IvAngleIJ)*diag(Volt)-diag(Y.*sin(mat_AngleIJ)*Volt')*diag(Volt);
% %J=Volt'*ones(1,Busnum).*cos(AngleIJ);这里错了 N=-diag(Volt)*Y.*cos(mat_IvAngleIJ)*diag(Volt)+diag(Y.*cos(mat_AngleIJ)*Volt')*diag(Volt);
% J=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y; 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);
% %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;
%暂时改一下
t1=[J,L; t1=[J,L;
H,N; H,N;
]'; ]';