pes2014-wronglowervoltagebound/func_ddh.m

187 lines
7.1 KiB
Mathematica
Raw Normal View History

2012-05-22 11:33:21 +08:00
function ddh=func_ddh(AngleIJMat,GB,Volt,Init_Y,Busnum)
%% deltaPi/deltaThytai_deltaThytaj <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%(<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>˶Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>)
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%ȥ<EFBFBD><EFBFBD><EFBFBD>˶Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>
t3=repmat(Init_Y',1,size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPidTidTj_2=dPidTidTj_2.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));
t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat);
t3=t1.*t2;
dPjdTidTj=t3-diag(diag(t3));%ȥ<EFBFBD><EFBFBD><EFBFBD>˶Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPdTidTj=dPidTidTj_2+dPjdTidTj;%<EFBFBD><EFBFBD><EFBFBD>շǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaP/deltaThyta_deltaThyta <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=sum(-dPidTidTj_2,2);
t2=diag(t1'.*Init_Y(1:2:size(Init_Y,2)));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t3=sum(-dPidTidTj_2,1);
t4=diag(t3.*Init_Y(1:2:size(Init_Y,2)));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPdTidTi=t2+t4;%%<EFBFBD><EFBFBD><EFBFBD>նԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaP/deltaThytai_dVi <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
t2=diag(t1);
t3=t1-diag(t2);%ȥ<EFBFBD><EFBFBD><EFBFBD>˶Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>
t4=sum(t3,2);
t4=t4'.*Init_Y(1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPidTidVi=diag(t4);
dPidTjdVj=-t3;
t6=sum(dPidTjdVj,1);
t6=t6.*Init_Y(1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPdTidVi=dPidTidVi+diag(t6);%%<EFBFBD><EFBFBD><EFBFBD>նԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaP/deltaThytai_dVj <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%ȥ<EFBFBD><EFBFBD><EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%ȥ<EFBFBD><EFBFBD><EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPdTidVj=dPidTidVj+dPjdTidVj;%<EFBFBD><EFBFBD><EFBFBD>շǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaP/dVi_deltaThytaj <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidVidTj=t1-diag(diag(t1));
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidTj=dPidVidTj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t1=-ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdVidTj=t1-diag(diag(t1));
t2=t2';
dPjdVidTj=dPjdVidTj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPdVidTj=dPidVidTj+dPjdVidTj;%<EFBFBD><EFBFBD><EFBFBD>շǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaPi/dVi_deltaThyta <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
dPdVidTi=dPdTidVi;%<EFBFBD><EFBFBD><EFBFBD>նԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@
%% deltaP/dVi_dVj <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%ȥ<EFBFBD><EFBFBD><EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD>ص<EFBFBD>
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t1=-(real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat));
dPjdVidVj=t1-diag(diag(t1));
t2=t2';
dPjdVidVj=dPjdVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPdVidVj=dPidVidVj+dPjdVidVj;%<EFBFBD><EFBFBD><EFBFBD>շǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@@
%% deltaP/dVi_dVi <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t0=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t1=diag(t0);
t2=t1'.*Init_Y(1:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dPidVidVi=-2*diag(t2);
dPidVjdVj=0;
dPdVidVi=dPidVidVi+dPidVjdVj;%<EFBFBD><EFBFBD><EFBFBD>նԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>APi
APi=zeros(2*Busnum,2*Busnum);
APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTj;%%<EFBFBD>ǶԽ<EFBFBD> TT
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVj;%%<EFBFBD>ǶԽ<EFBFBD> TV
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTj;%%<EFBFBD>ǶԽ<EFBFBD> VT
APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVj;%%<EFBFBD>ǶԽ<EFBFBD> VV
APi(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTi;%%<EFBFBD>Խ<EFBFBD>
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%<EFBFBD>Խ<EFBFBD>
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%<EFBFBD>Խ<EFBFBD>
APi(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVi;%%<EFBFBD>Խ<EFBFBD>
%% deltaQ/deltaThyta_deltaThyta <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
%t1=-Volt'*Volt;
t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
dQidTidTj=t1.*t2;%<EFBFBD><EFBFBD>ȥ<EFBFBD><EFBFBD><EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD>ˣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t2=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat);
dQjdTidTj=-t1.*t2;
t3=t3';
dQjdTidTj=dQjdTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQdTidTj=dQidTidTj+dQjdTidTj;%<EFBFBD><EFBFBD><EFBFBD>շǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaQ/deltaThyta_deltaThyta <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=dQidTidTj-diag(diag(dQidTidTj));%ȥ<EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t2=sum(t1,2);
t3=t2'.*Init_Y(2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
%dQidTidTi=diag(t3);
dQidTidTi=-diag(t3);
t1=-Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
t5=repmat(Init_Y',1,size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));
t7=sum(t6,1);
dQjdTidTi=diag(t7);
dQdTidTi=dQjdTidTi+dQidTidTi;%<EFBFBD><EFBFBD><EFBFBD>նԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaQ/deltaThyta_deltaV <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
%t1=-Volt;
t1=Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t3=t1'*ones(1,Busnum).*t2;
t4=repmat(Init_Y',1,size(Init_Y,2));
dQidTidVj=t3.*t4(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t2=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat);
t3=-Volt'*ones(1,Busnum).*t2;
t4=t4';
dQjdTidVj=t3.*t4(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQdTidVj=dQidTidVj+dQjdTidVj;%<EFBFBD><EFBFBD><EFBFBD>շǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD> @@@
%% deltaQ/deltaThyta_deltaV <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%ȥ<EFBFBD><EFBFBD><EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=ones(Busnum,1)*Volt.*t1;
t3=sum(t2,2);
t4=t3'.*Init_Y(2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQidTidVi=diag(t4);
t1=-Volt'*ones(1,Busnum).*(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
t2=t1-diag(diag(t1));
t3=repmat(Init_Y',1,size(Init_Y,2));
t4=t2.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
t5=sum(t4,1);
dQjdTidVi=diag(t5);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@
%% deltaQ/deltaV_deltaV <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t2=repmat(Init_Y',1,size(Init_Y,2));
t3=t1.*t2(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQidVidVj=t3;
t1=-real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t2=t2';
t3=t1.*t2(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQjdVidVj=t3;
dQdVidVj=dQidVidVj+dQjdVidVj;% @@@
%% deltaQ/deltaV_deltaV <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
t1=2*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
t2=diag(t1);
t3=t2'.*Init_Y(2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQidVidVi=diag(t3);
% t1=-real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
% t2=diag(t1);
% t3=t1-diag(t2);
% t4=sum(t3,1);
% t5=t4.*Init_Y(2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@
%% deltaQ/deltaV_deltaThyta <EFBFBD>ǶԽ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
%t1=real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat);
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*(t1);
t3=repmat(Init_Y',1,size(Init_Y,2));
t4=t2'.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQidVidTj=t4;
t1=real(GB).*cos(AngleIJMat)-imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*(t1);
t3=t3';
t4=t2.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%<EFBFBD><EFBFBD>y<EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @
%% deltaQ/deltaV_deltaThyta <EFBFBD>Խ<EFBFBD>Ԫ<EFBFBD><EFBFBD>
dQdVidTi=dQdTidVi;% @
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>AQi
AQi=zeros(2*Busnum,2*Busnum);
AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTj;%%<EFBFBD>ǶԽ<EFBFBD> TT
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVj;%%<EFBFBD>ǶԽ<EFBFBD> TV
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTj;%%<EFBFBD>ǶԽ<EFBFBD> VT
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVj;%%<EFBFBD>ǶԽ<EFBFBD> VV
AQi(1:2:2*Busnum,1:2:2*Busnum)=dQdTidTi;%%<EFBFBD>Խ<EFBFBD>
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%<EFBFBD>Խ<EFBFBD>
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%<EFBFBD>Խ<EFBFBD>
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVi;%%<EFBFBD>Խ<EFBFBD>
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ddh
t=[zeros(4,14);
zeros(2*5,4),AQi+APi;
];
ddh=t;
end