删除多余的函数文件。

Signed-off-by: facat <dmy@dmy-PC.(none)>
This commit is contained in:
facat 2012-05-25 17:23:23 +08:00
parent e515e5fd4a
commit 97822756be
26 changed files with 100 additions and 1891 deletions

View File

@ -1,9 +1,23 @@
function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX)
deltZ=XX(1:14);
deltL=XX(15:28);
deltW=XX(29:42);
deltU=XX(43:56);
deltX=XX(57:70);
deltY=XX(71:80);
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

View File

@ -1,23 +0,0 @@
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);
k1=0;
k2=0;
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(53:66);
deltU=XX(67:80);
end

View File

@ -1,23 +0,0 @@
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

View File

@ -1,12 +0,0 @@
function AA=FormAA(L_1Z,deltG,U_1W,Hcoma,deltH)
tOnes=eye(14);
tZeros=zeros(14);
AA=[
tOnes,L_1Z,tZeros,tZeros,tZeros,zeros(14,10);
tZeros,tOnes,tZeros,tZeros,-deltG',zeros(14,10);
tZeros,tZeros,tOnes,U_1W,tZeros,zeros(14,10);
tZeros,tZeros,tZeros,tOnes,deltG',zeros(14,10);
tZeros,tZeros,tZeros,tZeros,Hcoma,deltH;
zeros(10,14),zeros(10,14),zeros(10,14),zeros(10,14),deltH',zeros(10,10);
];
end

View File

@ -1,17 +0,0 @@
function AA=FormAA1(deltG,deltdeltF,ddh,ddg,deltH,Init_L,Init_U,Init_W,Init_Z,Busnum,PVi,PGi,RestraintCount,Balance)
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2;
H=-deltdeltF+ddh;%+ddg;
AA=[
H,deltH,deltG,deltG,zeros(ContrlCount,RestraintCount),zeros(ContrlCount,RestraintCount);
deltH',zeros(2*Busnum,2*Busnum),zeros(2*Busnum,RestraintCount),zeros(2*Busnum,RestraintCount),zeros(2*Busnum,RestraintCount),zeros(2*Busnum,RestraintCount);
deltG',zeros(RestraintCount,2*Busnum),zeros(RestraintCount,RestraintCount),zeros(RestraintCount,RestraintCount),-eye(RestraintCount,RestraintCount),zeros(RestraintCount,RestraintCount);
deltG',zeros(RestraintCount,2*Busnum),zeros(RestraintCount),zeros(RestraintCount),zeros(RestraintCount),eye(RestraintCount);
zeros(RestraintCount,ContrlCount),zeros(RestraintCount,2*Busnum),diag(Init_L),zeros(RestraintCount),diag(Init_Z),zeros(RestraintCount);
zeros(RestraintCount,ContrlCount),zeros(RestraintCount,2*Busnum),zeros(RestraintCount),diag(Init_U),zeros(RestraintCount),diag(Init_W);
];
%´¦ÀíÆ½ºâ½Úµã
t=size(PVi,1)+size(PGi,1);
AA(t+2*Balance-1,:)=0;
AA(:,t+2*Balance-1)=0;
AA(t+2*Balance-1,t+2*Balance-1)=1;
end

View File

@ -1,11 +1,10 @@
function YY=FormYY(Init_L,Lul,Lz,Ly,Init_U,Luu,Lw,LxComa)
t=[
-inv(diag(Init_L))*Lul;
Lz;
-inv(diag(Init_U))*Luu;
-Lw;
LxComa;
function YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx)
YY=[
Lx;
-Ly;
-Lz;
-Lw;
-Lul;
-Luu;
];
YY=t;
end

View File

@ -1,10 +0,0 @@
function YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx)
YY=[
Lx;
-Ly;
-Lz;
-Lw;
-Lul;
-Luu;
];
end

6
OPF.m
View File

@ -36,7 +36,7 @@ while(abs(Gap)>Precision)
%%
deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
%% ddHy
ddh=func_ddh3(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
%% ddg
ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi);
%% deltF
@ -51,11 +51,11 @@ while(abs(Gap)>Precision)
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%%
XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
%%
[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum);
[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]=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,Loadi);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1;

View File

@ -1,185 +0,0 @@
function ddh=func_ddh(AngleIJMat,GB,Volt,Init_Y,Busnum)
%% deltaPi/deltaThytai_deltaThytaj 非对角元素
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%(保留了对角元素的)
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%去掉了对角元素的
t3=repmat(Init_Y',1,size(Init_Y,2));%乘y的系数
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));%去掉了对角元素的
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
dPdTidTj=dPidTidTj_2+dPjdTidTj;%最终非对角元素 @@@
%% deltaP/deltaThyta_deltaThyta 对角元素
t1=sum(-dPidTidTj_2,2);
t2=diag(t1'.*Init_Y(1:2:size(Init_Y,2)));%乘y的系数
t3=sum(-dPidTidTj_2,1);
t4=diag(t3.*Init_Y(1:2:size(Init_Y,2)));%乘y的系数
dPdTidTi=t2+t4;%%最终对角元素 @@@
%% 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);
dPidTjdVj=-t3;
t6=sum(dPidTjdVj,1);
t6=t6.*Init_Y(1:2:size(Init_Y,2));%乘y的系数
dPdTidVi=dPidTidVi+diag(t6);%%最终对角元素 @@@
%% deltaP/deltaThytai_dVj 非对角元素
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%去掉对角元素的
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%去掉对角元素的
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
dPdTidVj=dPidTidVj+dPjdTidVj;%最终非对角元素 @@@
%% deltaP/dVi_deltaThytaj 非对角元素
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));%乘y的系数
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));%乘y的系数
dPdVidTj=dPidVidTj+dPjdVidTj;%最终非对角元素 @@@
%% deltaPi/dVi_deltaThyta 对角元素
dPdVidTi=dPdTidVi;%最终对角元素 @@
%% deltaP/dVi_dVj 非对角元素
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%去掉对角元素的
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
dPdVidVj=dPidVidVj+dPjdVidVj;%最终非对角元素 @@@@
%% deltaP/dVi_dVi 对角元素
t0=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t1=diag(diag(t0));
t2=t1'.*Init_Y(1:2:size(Init_Y,2));%乘y的系数
dPidVidVi=-2*diag(t2);
dPidVjdVj=0;
dPdVidVi=dPidVidVi+dPidVjdVj;%最终对角元素 @@@
%% 生成APi
APi=zeros(2*Busnum,2*Busnum);
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)=dPdTidTi;%%对角
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%对角
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%对角
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);
dQidTidTj=t1.*t2;%不去掉对角元素了,反正最后是要修正的
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
dQdTidTj=dQidTidTj+dQjdTidTj;%最终非对角元素 @@@
%% deltaQ/deltaThyta_deltaThyta 对角元素
t1=dQidTidTj-diag(diag(dQidTidTj));%去对角元素
t2=sum(t1,2);
t3=t2'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数
%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));%乘y的系数
t6=t4.*t5;
t7=sum(t6,1);
dQjdTidTi=diag(t7);
dQdTidTi=dQjdTidTi+dQidTidTi;%最终对角元素 @@@
%% deltaQ/deltaThyta_deltaV 非对角元素
%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));%乘y的系数
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));%乘y的系数
dQdTidVj=dQidTidVj+dQjdTidVj;%最终非对角元素 @@@
%% deltaQ/deltaThyta_deltaV 对角元素
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%去掉对角元素
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=ones(Busnum,1)*Volt.*t
t2=t1'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数
dQidTidVi=diag(t2);
t1=-Volt'*ones(1,Busnum).*(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
t2=t1-diag(diag(t1));
t3=sum(t2,1);
t4=t3.*Init_Y(2:2:size(Init_Y,2));%乘y的系数
dQjdTidVi=diag(t4);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@
%% deltaQ/deltaV_deltaV 非对角元素
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));%乘y的系数
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));%乘y的系数
dQjdVidVj=t3;
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);
% 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));%乘y的系数
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@
%% deltaQ/deltaV_deltaThyta 非对角元素
%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));%乘y的系数
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));%乘y的系数
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @
%% deltaQ/deltaV_deltaThyta 对角元素
dQdVidTi=dQdTidVi;% @
%% 生成AQi
AQi=zeros(2*Busnum,2*Busnum);
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)=dQdTidTi;%%对角
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%对角
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%对角
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVi;%%对角
%% 生成ddh
t=[zeros(4,14);
zeros(2*5,4),AQi+APi;
];
ddh=t;
end

View File

@ -1,187 +1,55 @@
function ddh=func_ddh(AngleIJMat,GB,Volt,Init_Y,Busnum)
%% deltaPi/deltaThytai_deltaThytaj
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%()
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%
t3=repmat(Init_Y',1,size(Init_Y,2));%y
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));%
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
dPdTidTj=dPidTidTj_2+dPjdTidTj;% @@@
%% deltaP/deltaThyta_deltaThyta
t1=sum(-dPidTidTj_2,2);
t2=diag(t1'.*Init_Y(1:2:size(Init_Y,2)));%y
t3=sum(-dPidTidTj_2,1);
t4=diag(t3.*Init_Y(1:2:size(Init_Y,2)));%y
dPdTidTi=t2+t4;%% @@@
%% 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);
dPidTjdVj=-t3;
t6=sum(dPidTjdVj,1);
t6=t6.*Init_Y(1:2:size(Init_Y,2));%y
dPdTidVi=dPidTidVi+diag(t6);%% @@@
%% deltaP/deltaThytai_dVj
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
t1=ones(Busnum,1)*Volt.*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
dPdTidVj=dPidTidVj+dPjdTidVj;% @@@
%% deltaP/dVi_deltaThytaj
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));%y
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));%y
dPdVidTj=dPidVidTj+dPjdVidTj;% @@@
%% deltaPi/dVi_deltaThyta
dPdVidTi=dPdTidVi;% @@
%% deltaP/dVi_dVj
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
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));%y
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;% @@@
%% APi
APi=zeros(2*Busnum,2*Busnum);
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)=dPdTidTi;%%
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%
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);
dQidTidTj=t1.*t2;%
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
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));%y
dQdTidTj=dQidTidTj+dQjdTidTj;% @@@
%% deltaQ/deltaThyta_deltaThyta
t1=dQidTidTj-diag(diag(dQidTidTj));%
t2=sum(t1,2);
t3=t2'.*Init_Y(2:2:size(Init_Y,2));%y
%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));%y
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;% @@@
%% deltaQ/deltaThyta_deltaV
%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));%y
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));%y
dQdTidVj=dQidTidVj+dQjdTidVj;% @@@
%% deltaQ/deltaThyta_deltaV
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%
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));%y
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));%y
t5=sum(t4,1);
dQjdTidVi=diag(t5);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@
%% deltaQ/deltaV_deltaV
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));%y
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));%y
dQjdVidVj=t3;
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);
% 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));%y
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@
%% deltaQ/deltaV_deltaThyta
%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));%y
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));%y
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @
%% deltaQ/deltaV_deltaThyta
dQdVidTi=dQdTidVi;% @
%% AQi
AQi=zeros(2*Busnum,2*Busnum);
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)=dQdTidTi;%%
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVi;%%
%% ddh
t=[zeros(4,14);
zeros(2*5,4),AQi+APi;
function ddh=func_ddh3(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
%
%ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
%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=-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,ddPdTdV+ddQdTdV ;
ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT;
];
ddh=t;
sizePGi=size(PGi,1);
sizePVi=size(PVi,1);
sizeLoadi=size(Loadi,1);
ddh=[
sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount);
sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t;
];
end

View File

@ -1,203 +0,0 @@
function ddh=func_ddh1(AngleIJMat,GB,Volt,Init_Y,Busnum)
%% deltaPi/deltaThytai_deltaThytaj 非对角元素
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%(保留了对角元素的)
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%去掉了对角元素的
t3=repmat(Init_Y',1,size(Init_Y,2));%乘y的系数
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));%去掉了对角元素的
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
dPdTidTj=dPidTidTj_2+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=t2-diag(diag(t2));
t3=repmat(Init_Y',1,size(Init_Y,2));%乘y的系数
t4=t1.*t2.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));
t5=sum(t4,1);
dPidTjdTj=diag(t5);
dPdTidTi=dPidTidTi+dPidTjdTj;%%最终对角元素 @@@@@@
%% 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));
t1=t1-diag(diag(t1));%去掉对角元素
t2=repmat(Init_Y',1,size(Init_Y,2));
t3=t1.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
t6=sum(t3,1);
dPdTidVi=dPidTidVi+diag(t6);%%最终对角元素 @@@@@@
%% deltaP/deltaThytai_dVj 非对角元素
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%去掉对角元素的
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%去掉对角元素的
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
dPdTidVj=dPidTidVj+dPjdTidVj;%最终非对角元素 @@@@@
%% deltaP/dVi_deltaThytaj 非对角元素
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));%乘y的系数
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));%乘y的系数
dPdVidTj=dPidVidTj+dPjdVidTj;%最终非对角元素 @@@@
%% deltaPi/dVi_deltaThyta 对角元素
dPdVidTi=dPdTidVi;%最终对角元素 @@
%% deltaP/dVi_dVj 非对角元素
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%去掉对角元素的
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
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;%最终对角元素 @@@@
%% 生成APi
APi=zeros(2*Busnum,2*Busnum);
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)=dPdTidTi;%%对角
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%对角
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%对角
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);
dQidTidTj=t1.*t2;%不去掉对角元素了,反正最后是要修正的
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
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=repmat(Init_Y',1,size(Init_Y,2));
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
t7=sum(t6,2);
dQidTidTi=diag(t7);
t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
%t5=t5';
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
t7=sum(t6,1);
dQjdTidTi=diag(t7);
dQdTidTi=dQjdTidTi+dQidTidTi;%最终对角元素 @@@@
%% deltaQ/deltaThyta_deltaV 非对角元素
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));%乘y的系数
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));%乘y的系数
dQdTidVj=dQidTidVj+dQjdTidVj;%最终非对角元素 @@@@@
%% deltaQ/deltaThyta_deltaV 对角元素
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%去掉对角元素
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*t1;
t2=t2-diag(diag(t2));
t3=sum(t2,2);
t4=t3'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
t5=sum(t4,1);
dQjdTidVi=diag(t5);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@@@
%% deltaQ/deltaV_deltaV 非对角元素
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));%乘y的系数
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));%乘y的系数
dQjdVidVj=t3;
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);
% 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));%乘y的系数
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@@@
%% deltaQ/deltaV_deltaThyta 非对角元素
%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));%乘y的系数
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));%乘y的系数
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @@@
%% deltaQ/deltaV_deltaThyta 对角元素
dQdVidTi=dQdTidVi;% @
%% 生成AQi
AQi=zeros(2*Busnum,2*Busnum);
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)=dQdTidTi;%%对角
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%对角
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%对角
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVi;%%对角
%% 生成ddh
t=[zeros(4,14);
zeros(2*Busnum,2*size(PVi,1)),AQi+APi;
];
ddh=t;
end

View File

@ -1,204 +0,0 @@
function ddh=func_ddh1(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi)
%% deltaPi/deltaThytai_deltaThytaj
ContrlCount=size(PVi,1)*2+Busnum*2;
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%()
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%
t3=repmat(Init_Y',1,size(Init_Y,2));%y
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));%
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
dPdTidTj=dPidTidTj_2+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=t2-diag(diag(t2));
t3=repmat(Init_Y',1,size(Init_Y,2));%y
t4=t1.*t2.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));
t5=sum(t4,1);
dPidTjdTj=diag(t5);
dPdTidTi=dPidTidTi+dPidTjdTj;%% @@@@@@
%% 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));
t1=t1-diag(diag(t1));%
t2=repmat(Init_Y',1,size(Init_Y,2));
t3=t1.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
t6=sum(t3,1);
dPdTidVi=dPidTidVi+diag(t6);%% @@@@@@
%% deltaP/deltaThytai_dVj
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
dPdTidVj=dPidTidVj+dPjdTidVj;% @@@@@
%% deltaP/dVi_deltaThytaj
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));%y
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));%y
dPdVidTj=dPidVidTj+dPjdVidTj;% @@@@
%% deltaPi/dVi_deltaThyta
dPdVidTi=dPdTidVi;% @@
%% deltaP/dVi_dVj
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
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));%y
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;% @@@@
%% APi
APi=zeros(2*Busnum,2*Busnum);
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)=dPdTidTi;%%
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%
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);
dQidTidTj=t1.*t2;%
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
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));%y
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=repmat(Init_Y',1,size(Init_Y,2));
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
t7=sum(t6,2);
dQidTidTi=diag(t7);
t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
%t5=t5';
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
t7=sum(t6,1);
dQjdTidTi=diag(t7);
dQdTidTi=dQjdTidTi+dQidTidTi;% @@@@
%% deltaQ/deltaThyta_deltaV
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));%y
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));%y
dQdTidVj=dQidTidVj+dQjdTidVj;% @@@@@
%% deltaQ/deltaThyta_deltaV
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*t1;
t2=t2-diag(diag(t2));
t3=sum(t2,2);
t4=t3'.*Init_Y(2:2:size(Init_Y,2));%y
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));%y
t5=sum(t4,1);
dQjdTidVi=diag(t5);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@@@
%% deltaQ/deltaV_deltaV
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));%y
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));%y
dQjdVidVj=t3;
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);
% 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));%y
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@@@
%% deltaQ/deltaV_deltaThyta
%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));%y
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));%y
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @@@
%% deltaQ/deltaV_deltaThyta
dQdVidTi=dQdTidVi;% @
%% AQi
AQi=zeros(2*Busnum,2*Busnum);
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)=dQdTidTi;%%
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%
AQi(2:2:2*Busnum,2:2:2*Busnum)=dQdVidVi;%%
%% ddh
t=[zeros(2*size(PVi,1),ContrlCount);
zeros(2*Busnum,2*size(PVi,1)),AQi+APi;
];
ddh=t;
end

View File

@ -1,205 +0,0 @@
function ddh=func_ddh2(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi)
%% deltaPi/deltaThytai_deltaThytaj 非对角元素
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2;
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%(保留了对角元素的)
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%去掉了对角元素的
t3=repmat(Init_Y',1,size(Init_Y,2));%乘y的系数
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));%去掉了对角元素的
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
dPdTidTj=dPidTidTj_2+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;
t2=t2-diag(diag(t2));
t3=repmat(Init_Y',1,size(Init_Y,2));%乘y的系数
t4=t2.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));
t5=sum(t4,1);
dPidTjdTj=diag(t5);
dPdTidTi=dPidTidTi+dPidTjdTj;%%最终对角元素 @@@@@@@
%% 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));
t1=t1-diag(diag(t1));%去掉对角元素
t2=repmat(Init_Y',1,size(Init_Y,2));
t3=t1.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
t6=sum(t3,1);
dPdTidVi=dPidTidVi+diag(t6);%%最终对角元素 @@@@@@
%% deltaP/deltaThytai_dVj 非对角元素
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%去掉对角元素的
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%去掉对角元素的
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
dPdTidVj=dPidTidVj+dPjdTidVj;%最终非对角元素 @@@@@@
%% deltaP/dVi_deltaThytaj 非对角元素
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));%乘y的系数
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));%乘y的系数
dPdVidTj=dPidVidTj+dPjdVidTj;%最终非对角元素 @@@@@
%% deltaPi/dVi_deltaThyta 对角元素
dPdVidTi=dPdTidVi;%最终对角元素 @@
%% deltaP/dVi_dVj 非对角元素
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%去掉对角元素的
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
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;%最终对角元素 @@@@@
%% 生成APi
APi=zeros(2*Busnum,2*Busnum);
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)=dPdTidTi;%%对角
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%对角
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%对角
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);
dQidTidTj=t1.*t2;%不去掉对角元素了,反正最后是要修正的
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
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=repmat(Init_Y',1,size(Init_Y,2));
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
t7=sum(t6,2);
dQidTidTi=diag(t7);
t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
%t5=t5';
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%乘y的系数
t7=sum(t6,1);
dQjdTidTi=diag(t7);
dQdTidTi=dQjdTidTi+dQidTidTi;%最终对角元素 @@@@@
%% deltaQ/deltaThyta_deltaV 非对角元素
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));%乘y的系数
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));%乘y的系数
dQdTidVj=dQidTidVj+dQjdTidVj;%最终非对角元素 @@@@@@
%% deltaQ/deltaThyta_deltaV 对角元素
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%去掉对角元素
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*t1;
t2=t2-diag(diag(t2));
t3=sum(t2,2);
t4=t3'.*Init_Y(2:2:size(Init_Y,2));%乘y的系数
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));%乘y的系数
t5=sum(t4,1);
dQjdTidVi=diag(t5);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@@@
%% deltaQ/deltaV_deltaV 非对角元素
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));%乘y的系数
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));%乘y的系数
dQjdVidVj=t3;
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);
% 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));%乘y的系数
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@@@
%% deltaQ/deltaV_deltaThyta 非对角元素
%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));%乘y的系数
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));%乘y的系数
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @@@
%% deltaQ/deltaV_deltaThyta 对角元素
dQdVidTi=dQdTidVi;% @
%% 生成AQi
AQi=zeros(2*Busnum,2*Busnum);
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)=dQdTidTi;%%对角
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%对角
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%对角
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;
end

View File

@ -1,206 +0,0 @@
function ddh=func_ddh2(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi)
%% deltaPi/deltaThytai_deltaThytaj
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*2;
t1=-Volt'*Volt;
t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
dPdTidTj=t1.*t2; %%()
dPidTidTj_2=dPdTidTj-diag(diag(dPdTidTj));%
t3=repmat(Init_Y',1,size(Init_Y,2));%y
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));%
t3=repmat(Init_Y,size(Init_Y,2),1);
dPjdTidTj=dPjdTidTj.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
dPdTidTj=dPidTidTj_2+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;
t2=t2-diag(diag(t2));
t3=repmat(Init_Y',1,size(Init_Y,2));%y
t4=t2.*t3(1:2:size(Init_Y,2),1:2:size(Init_Y,2));
t5=sum(t4,1);
dPidTjdTj=diag(t5);
dPdTidTi=dPidTidTi+dPidTjdTj;%% @@@@@@@
%% 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));
t1=t1-diag(diag(t1));%
t2=repmat(Init_Y',1,size(Init_Y,2));
t3=t1.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
t6=sum(t3,1);
dPdTidVi=dPidTidVi+diag(t6);%% @@@@@@
%% deltaP/deltaThytai_dVj
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat));
dPidTidVj=t1-diag(diag(t1));%%
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidTidVj=dPidTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
t1=Volt'*ones(1,Busnum).*(real(GB).*sin(AngleIJMat)+imag(GB).*cos(AngleIJMat));
dPjdTidVj=t1-diag(diag(t1));%
t2=t2';
dPjdTidVj=dPjdTidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
dPdTidVj=dPidTidVj+dPjdTidVj;% @@@@@@
%% deltaP/dVi_deltaThytaj
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));%y
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));%y
dPdVidTj=dPidVidTj+dPjdVidTj;% @@@@@
%% deltaPi/dVi_deltaThyta
dPdVidTi=dPdTidVi;% @@
%% deltaP/dVi_dVj
t1=-(real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat));
dPidVidVj=t1-diag(diag(t1));%
t2=repmat(Init_Y',1,size(Init_Y,2));
dPidVidVj=dPidVidVj.*t2(1:2:size(Init_Y,2),1:2:size(Init_Y,2));%y
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));%y
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;% @@@@@
%% APi
APi=zeros(2*Busnum,2*Busnum);
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)=dPdTidTi;%%
APi(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%
APi(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%
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);
dQidTidTj=t1.*t2;%
t3=repmat(Init_Y',1,size(Init_Y,2));
dQidTidTj=dQidTidTj.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
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));%y
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=repmat(Init_Y',1,size(Init_Y,2));
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
t7=sum(t6,2);
dQidTidTi=diag(t7);
t1=Volt'*Volt;
t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t3=t1.*t2;
t4=t3-diag(diag(t3));
%t5=t5';
t6=t4.*t5(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
t7=sum(t6,1);
dQjdTidTi=diag(t7);
dQdTidTi=dQjdTidTi+dQidTidTi;% @@@@@
%% deltaQ/deltaThyta_deltaV
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));%y
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));%y
dQdTidVj=dQidTidVj+dQjdTidVj;% @@@@@@
%% deltaQ/deltaThyta_deltaV
%t1=sum(dQidTidVj,2)-diag(dQidTidVj);%
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=-ones(Busnum,1)*Volt.*t1;
t2=t2-diag(diag(t2));
t3=sum(t2,2);
t4=t3'.*Init_Y(2:2:size(Init_Y,2));%y
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));%y
t5=sum(t4,1);
dQjdTidVi=diag(t5);
dQdTidVi=dQidTidVi+dQjdTidVi;% @@@@
%% deltaQ/deltaV_deltaV
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));%y
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));%y
dQjdVidVj=t3;
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);
% 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));%y
% dQjdVidVi=diag(t5);
dQjdVidVi=0;
dQdVidVi=dQidVidVi+dQjdVidVi; % @@@@@
%% deltaQ/deltaV_deltaThyta
%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));%y
t4=t2.*t3(2:2:size(Init_Y,2),2:2:size(Init_Y,2));%y
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));%y
dQjdVidTj=t4;
dQdVidTj=dQidVidTj+dQjdVidTj;% @@@@
%% deltaQ/deltaV_deltaThyta
dQdVidTi=dQdTidVi;% @
%% AQi
AQi=zeros(2*Busnum,2*Busnum);
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)=dQdTidTi;%%
AQi(1:2:2*Busnum,2:2:2*Busnum)=dQdTidVi;%%
AQi(2:2:2*Busnum,1:2:2*Busnum)=dQdVidTi;%%
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;
end

View File

@ -1,52 +0,0 @@
function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y)
%决定用循环重写
ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3;
AngleIJ=AngleIJMat-angle(GB);
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
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
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;%存疑与我的不一样
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');
t6=-(t4+t5)*diag(Volt);
ddQdTdT=t3+t6;%ok1
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');
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);
ddQdTdV=t1+t2+t3+t4;
t1=Y.*sin(AngleIJ')*diag(yQ);
t2=diag(yQ)*Y.*sin(AngleIJ);
ddQdVdV=t1+t2;
%% %%
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;
];
ddh=t;
end

View File

@ -1,55 +0,0 @@
function ddh=func_ddh3(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
%
%ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
%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=-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,ddPdTdV+ddQdTdV ;
ddPdVdT+ddQdVdT,ddPdTdT+ddQdTdT;
];
sizePGi=size(PGi,1);
sizePVi=size(PVi,1);
sizeLoadi=size(Loadi,1);
ddh=[
sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount);
sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t;
];
end

View File

@ -3,6 +3,6 @@ function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c); %
dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c); %
deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
end

View File

@ -1,60 +0,0 @@
function Jacob=jacobian_M(Busnum,Volt,Y,Angle,AngleIJMat)
%**************************************************************************
% 程序功能 : 子函数——形成雅可比矩阵Jacobian
% 编 者:
% 编制时间2010.12
%**************************************************************************
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
temp1=-Volt'*Volt.*Y;
AngleIJ=AngleIJMat-Angle;
temp11=Volt'*ones(1,Busnum).*Y;
temp2=sum(temp1.*sin(AngleIJ),2);
temp22 = sum(temp11.*sin(AngleIJ),2);
temp3 = sum(temp1.*cos(AngleIJ),2);
temp33 = sum(temp11.*cos(AngleIJ),2);
temp4=diag(temp2);
temp44=diag(temp22);
temp5=diag(temp3);
temp55=diag(temp33);
%计算Lii的累加项
t1=ones(Busnum,1)*Volt.*Y;
t2=sum(t1.*sin(AngleIJ),2);
t3=sum(t1.*cos(AngleIJ),2);
t4=diag(t2);
H = temp1.*sin(AngleIJ)-temp4;%
L = -temp11.*sin(AngleIJ);%
%L(1:Busnum,1:Busnum)=-temp44+;
L=L-t3.*
N=-temp11.*cos(AngleIJ);%
N(1:Busnum,1:Busnum)=-temp55-diag(diag(temp11.*cos(Angle) ) );
J = -temp1.*cos(AngleIJ)+temp5;%
%%
%Q = Q0+temp2'; %求有功分量P
%P = P0+temp3'; %求无功分量Q
%% 处理平衡节点和pv节点
% H(:,Balance) = 0;
% H(Balance,:) = 0;
% H(Balance,Balance) = 100; % 平衡节点对应的对角元素置一个有限数
% L(:,PVi) = 0;
% L(PVi,:) = 0;
% L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV节点对应的对角元素置为1
% J(:,Balance) = 0;
% J(PVi,:) = 0;
% N(:,PVi) = 0;
% N(Balance,:) = 0;
% Q(PVi) = 0; % 将pv节点的无功不平衡分量置零
% P(Balance) = 0; % 平衡节点的有功功率不平衡分量置零
%% 合成PQ和雅可比矩阵
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;
% t1(1:)
% PQ = cat(2,P,Q); % 形成功率不平衡分量列向量
%Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % 形成Jacobian矩阵
Jacob=t1;
end

View File

@ -1,61 +1,20 @@
function Jacob=jacobian_M(Busnum,Volt,Y,Angle,AngleIJMat)
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
%**************************************************************************
% : Jacobian
%
% 2010.12
%**************************************************************************
%%²ÎÕÕͼÊé¹İ6Â¥µÄÊé±àĞ´
%% H,L,N,JP,Q
temp1=-Volt'*Volt.*Y;
AngleIJ=AngleIJMat-Angle;
temp11=Volt'*ones(1,Busnum).*Y;
temp2=sum(temp1.*sin(AngleIJ),2);
temp22 = sum(temp11.*sin(AngleIJ),2);
temp3 = sum(temp1.*cos(AngleIJ),2);
temp33 = sum(temp11.*cos(AngleIJ),2);
temp4=diag(temp2);
temp44=diag(temp22);
temp5=diag(temp3);
temp55=diag(temp33);
%Lii
t1=ones(Busnum,1)*Volt.*Y;
t2=sum(t1.*sin(AngleIJ),2);
t3=sum(t1.*cos(AngleIJ),2);
t4=diag(t2);
t5=diag(t3);
H = temp1.*sin(AngleIJ)-temp4;%
L = -temp11.*sin(AngleIJ);%
%L(1:Busnum,1:Busnum)=-temp44+;
L=L-t4;
N=-temp11.*cos(AngleIJ);%
%N(1:Busnum,1:Busnum)=-temp55-diag(diag(temp11.*cos(Angle) ) );
N=N-t5;
J = -temp1.*cos(AngleIJ)+temp5;%
%%
%Q = Q0+temp2'; %P
%P = P0+temp3'; %Q
%% pv
% H(:,Balance) = 0;
% H(Balance,:) = 0;
% H(Balance,Balance) = 100; %
% L(:,PVi) = 0;
% L(PVi,:) = 0;
% L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV1
% J(:,Balance) = 0;
% J(PVi,:) = 0;
% N(:,PVi) = 0;
% N(Balance,:) = 0;
% Q(PVi) = 0; % pv
% P(Balance) = 0; %
%% PQ
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;
% t1(1:)
% PQ = cat(2,P,Q); %
%Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % Jacobian
Jacob=t1;
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

View File

@ -1,84 +0,0 @@
function [Jacob]=jacobian_M1(Busnum,Volt,Y,Angle,AngleIJMat)
%**************************************************************************
% : Jacobian
%
% 2010.12
%**************************************************************************
%% H,L,N,JP,Q
%Volt(PVi) = PVu;
temp1=Volt'*Volt.*Y;
AngleIJ=AngleIJMat-Angle;
temp2=sum(temp1.*sin(AngleIJ),2);
temp3 = sum(temp1.*cos(AngleIJ),2);
temp4=diag(temp2);
temp5=diag(temp3);
%t1=Volt'*ones(1,Busnum).*Y;
t1=ones(Busnum,1)*Volt.*Y;
%t1=Volt'*Volt.*Y;
t2=sum(t1.*sin(AngleIJ),2);
t3=sum(t1.*cos(AngleIJ),2);
t4=diag(t2);
t5=diag(t3);
H = -temp1.*sin(AngleIJ)+temp4;%
L = -t1.*sin(AngleIJ);%
%L(1:Busnum,1:Busnum)=-temp44+;
L=L-t4;
N=-t1.*cos(AngleIJ);%
%N(1:Busnum,1:Busnum)=-temp55-diag(diag(temp11.*cos(Angle) ) );
N=N-t5;
J = temp1.*cos(AngleIJ)-temp5;%
%%%%
%t=diag(Volt);
%N=t*N;%*t;
%L=t*L;%*t;
%%%%
%%
%Q
%
%% PQ
% t1(1:)
%
%Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % Jacobian
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;
end
% function Jacob=jacobian_M1(Busnum,PVi,PVu,U,Uangle,Y,Angle,r,c)
% AngleIJ = Uangle(r) - Uangle(c)- Angle';
% U(PVi) = PVu;
% temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum); %
% temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);
% temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);
% temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum);
% temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum);
% H = temp1.*sparse(r,c,sin(AngleIJ))-temp4;
% L = temp1.*sparse(r,c,sin(AngleIJ))+temp4;
% N = temp1.*sparse(r,c,cos(AngleIJ))+temp5;
% J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5;
%
%
% 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;
% % t1(1:)
% % PQ = cat(2,P,Q); %
% %Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % Jacobian
% Jacob=t1;
%
% end

View File

@ -1,84 +0,0 @@
function [Jacob]=jacobian_M2(Busnum,Volt,Y,Angle,AngleIJMat)
%**************************************************************************
% : Jacobian
%
% 2010.12
%**************************************************************************
%% H,L,N,JP,Q
temp1=Volt'*Volt.*Y;
AngleIJ=AngleIJMat-Angle;
temp2=sum(temp1.*sin(AngleIJ),2);
temp3 = sum(temp1.*cos(AngleIJ),2);
temp4=diag(temp2);
temp5=diag(temp3);
%t1=Volt'*ones(1,Busnum).*Y;
t1=ones(Busnum,1)*Volt.*Y;
t11=Volt'*ones(1,Busnum).*Y;
%t1=Volt'*Volt.*Y;
t2=sum(t1.*sin(AngleIJ),2);
t3=sum(t1.*cos(AngleIJ),2);
t4=diag(t2);
t5=diag(t3);
H = -temp1.*sin(AngleIJ)+temp4;%
L = -t11.*sin(AngleIJ);%
%L(1:Busnum,1:Busnum)=-temp44+;
L=L-t4;
N=-t11.*cos(AngleIJ);%
%N(1:Busnum,1:Busnum)=-temp55-diag(diag(temp11.*cos(Angle) ) );
N=N-t5;
J = temp1.*cos(AngleIJ)-temp5;%
%%%%
%t=diag(Volt);
%N=t*N;%*t;
%L=t*L;%*t;
%%%%
%%
%Q
%
%% PQ
% t1(1:)
%
%Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % Jacobian
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;
end
% function Jacob=jacobian_M1(Busnum,PVi,PVu,U,Uangle,Y,Angle,r,c)
% AngleIJ = Uangle(r) - Uangle(c)- Angle';
% U(PVi) = PVu;
% temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum); %
% temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);
% temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);
% temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum);
% temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum);
% H = temp1.*sparse(r,c,sin(AngleIJ))-temp4;
% L = temp1.*sparse(r,c,sin(AngleIJ))+temp4;
% N = temp1.*sparse(r,c,cos(AngleIJ))+temp5;
% J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5;
%
%
% 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;
% % t1(1:)
% % PQ = cat(2,P,Q); %
% %Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % Jacobian
% Jacob=t1;
%
% end

View File

@ -1,20 +0,0 @@
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,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: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;
]';
Jacob=-t1;
end

View File

@ -1,20 +0,0 @@
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
%**************************************************************************
% : Jacobian
%
% 2010.12
%**************************************************************************
%%6
%% H,L,N,JP,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

View File

@ -1,77 +0,0 @@
function [Jacob]=jacobian_M4(Busnum,Volt,Y,Angle,AngleIJMat)
%**************************************************************************
% 程序功能 : 子函数——形成雅可比矩阵Jacobian
% 编 者:
% 编制时间2010.12
%**************************************************************************
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
temp1=Volt'*Volt.*Y;
AngleIJ=AngleIJMat-Angle;
tt1=temp1.*sin(AngleIJ);
tt2=temp1.*cos(AngleIJ);
tt3=diag(tt1);
tt4=diag(tt2);
tt5=tt1-diag(tt3);
tt6=tt2-diag(tt4);
temp2=sum(tt5,2);
temp3 = sum(tt6,2);
HH=temp2;
JJ=temp3;
t1=ones(Busnum,1)*Volt.*Y;
t11=Volt'*ones(1,Busnum).*Y;
t2=sum(t1.*sin(AngleIJ),2);
t3=sum(t1.*cos(AngleIJ),2);
t4=diag(t1.*sin(AngleIJ));
t5=diag(t1.*cos(AngleIJ));
NN=-diag(t3)-diag(t5);
LL=-diag(t2)+diag(t4);
H = -temp1.*sin(AngleIJ);
L = -t11.*sin(AngleIJ);%
N=-t11.*cos(AngleIJ);%
J = temp1.*cos(AngleIJ);%
H=H-diag(diag(H));
N=N-diag(N));
J=J-diag(J);
L=L-diag(L);
H=H+diag(HH);
N=N+diag(NN);
J=J+diag(JJ);
L=L+diag(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;
end
% function Jacob=jacobian_M1(Busnum,PVi,PVu,U,Uangle,Y,Angle,r,c)
% AngleIJ = Uangle(r) - Uangle(c)- Angle';
% U(PVi) = PVu;
% temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum); % 计算雅克比矩阵可利用的中间变量
% temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);
% temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);
% temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum);
% temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum);
% H = temp1.*sparse(r,c,sin(AngleIJ))-temp4;
% L = temp1.*sparse(r,c,sin(AngleIJ))+temp4;
% N = temp1.*sparse(r,c,cos(AngleIJ))+temp5;
% J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5;
%
%
% 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;
% % t1(1:)
% % PQ = cat(2,P,Q); % 形成功率不平衡分量列向量
% %Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % 形成Jacobian矩阵
% Jacob=t1;
%
% end

View File

@ -1,77 +0,0 @@
function [Jacob]=jacobian_M4(Busnum,Volt,Y,Angle,AngleIJMat)
%**************************************************************************
% : Jacobian
%
% 2010.12
%**************************************************************************
%% H,L,N,JP,Q
temp1=Volt'*Volt.*Y;
AngleIJ=AngleIJMat-Angle;
tt1=temp1.*sin(AngleIJ);
tt2=temp1.*cos(AngleIJ);
tt3=diag(tt1);
tt4=diag(tt2);
tt5=tt1-diag(tt3);
tt6=tt2-diag(tt4);
temp2=sum(tt5,2);
temp3 = sum(tt6,2);
HH=temp2;
JJ=-temp3;
t1=ones(Busnum,1)*Volt.*Y;
t11=Volt'*ones(1,Busnum).*Y;
t2=sum(t1.*sin(AngleIJ),2);
t3=sum(t1.*cos(AngleIJ),2);
t4=diag(t1.*sin(AngleIJ));
t5=diag(t1.*cos(AngleIJ));
NN=-diag(t3)-diag(t5);
LL=-diag(t2)+diag(t4);
H = -temp1.*sin(AngleIJ);
L = -t11.*sin(AngleIJ);%
N=-t11.*cos(AngleIJ);%
J = temp1.*cos(AngleIJ);%
H=H-diag(diag(H));
N=N-diag(diag(N));
J=J-diag(diag(J));
L=L-diag(diag(L));
H=H+diag(HH);
N=N+NN;
J=J+diag(JJ);
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;
end
% function Jacob=jacobian_M1(Busnum,PVi,PVu,U,Uangle,Y,Angle,r,c)
% AngleIJ = Uangle(r) - Uangle(c)- Angle';
% U(PVi) = PVu;
% temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum); %
% temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);
% temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);
% temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum);
% temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum);
% H = temp1.*sparse(r,c,sin(AngleIJ))-temp4;
% L = temp1.*sparse(r,c,sin(AngleIJ))+temp4;
% N = temp1.*sparse(r,c,cos(AngleIJ))+temp5;
% J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5;
%
%
% 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;
% % t1(1:)
% % PQ = cat(2,P,Q); %
% %Jacob = cat(1,cat(2,H,N),cat(2,J,L)); % Jacobian
% Jacob=t1;
%
% end

14
sy.m
View File

@ -1,14 +0,0 @@
clear
clc
syms T11 T12 T21 T22;
syms V1 V2;
syms Y11 Y12 Y21 Y22;
yP=ones(1,2);
AngleIJ=[T11,T12;T21,T22];
Volt=[V1,V2];
Y=[Y11,Y12;Y21,Y22];
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=t1+t2+t3+t4