115 lines
3.2 KiB
Plaintext
115 lines
3.2 KiB
Plaintext
function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat)
|
||
%**************************************************************************
|
||
% 程序功能 : 子函数——形成雅可比矩阵Jacobian
|
||
% 编 者：
|
||
% 编制时间：2010.12
|
||
%**************************************************************************
|
||
%%参照图书馆6楼的书编写
|
||
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
|
||
AngleIJ=AngleIJMat-Angle;
|
||
t1=Volt'*Volt;
|
||
H=t1.*Y.*sin(AngleIJ);
|
||
N=-t1.*Y.*cos(AngleIJ);
|
||
%J=Volt'*ones(1,Busnum).*cos(AngleIJ);这里错了
|
||
J=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y;
|
||
|
||
%L=Volt'*ones(1,Busnum).*sin(AngleIJ);这里错了
|
||
L=Volt'*ones(1,Busnum).*sin(AngleIJ).*Y;
|
||
|
||
%%对角
|
||
t1=Volt'*Volt;
|
||
t2=t1.*Y.*sin(AngleIJ);
|
||
t3=diag(t2);
|
||
t4=t2-diag(t3);
|
||
t5=sum(t4,2);
|
||
HH=-diag(t5);
|
||
t2=t1.*Y.*cos(AngleIJ);
|
||
t3=diag(t2);
|
||
t4=t2-diag(t3);
|
||
t5=sum(t4,2);
|
||
NN=diag(t5);
|
||
%t1=ones(Busnum,1)*Volt;
|
||
t1=ones(Busnum,1)*Volt.*Y;
|
||
t2=t1.*cos(AngleIJ);
|
||
t3=sum(t2,2);
|
||
JJ=diag(t3);
|
||
t1=Volt'*ones(1,Busnum).*cos(AngleIJ).*Y;
|
||
%t1=Volt'*ones(1,Busnum).*cos(AngleIJ);
|
||
t2=diag(t1);%
|
||
JJ=JJ+diag(t2);
|
||
t1=ones(Busnum,1)*Volt.*Y;
|
||
%t1=ones(Busnum,1)*Volt;
|
||
t2=t1.*sin(AngleIJ);
|
||
t3=sum(t2,2);
|
||
LL=diag(t3);
|
||
%t1=Volt'*ones(1,Busnum).*sin(AngleIJ);
|
||
t1=Volt'*ones(1,Busnum).*sin(AngleIJ).*Y;
|
||
t2=diag(t1);%
|
||
%LL=LL-diag(t2);这里错了
|
||
LL=LL+diag(t2);
|
||
|
||
H=H-diag(diag(H));
|
||
N=N-diag(diag(N));
|
||
J=J-diag(diag(J));
|
||
L=L-diag(diag(L));
|
||
|
||
H=H+HH;
|
||
J=J+JJ;
|
||
N=N+NN;
|
||
L=L+LL;
|
||
|
||
t1=zeros(2*Busnum);
|
||
t1(1:2:2*Busnum,1:2:2*Busnum)=H;
|
||
t1(1:2:2*Busnum,2:2:2*Busnum)=N;
|
||
t1(2:2:2*Busnum,1:2:2*Busnum)=J;
|
||
t1(2:2:2*Busnum,2:2:2*Busnum)=L;
|
||
Jacob=-t1;
|
||
%%以下是学姐给的公式
|
||
H=diag(Volt)*Y.*sin(AngleIJ')*diag(Volt)-diag(Y.*sin(AngleIJ)*Volt')*diag(Volt);
|
||
N=-diag(Volt)*Y.*cos(AngleIJ')*diag(Volt)+diag(Y.*cos(AngleIJ)*Volt')*diag(Volt);
|
||
J=diag(Y.*cos(AngleIJ)*Volt')+Y.*cos(AngleIJ')*diag(Volt);
|
||
L=diag(Y.*sin(AngleIJ)*Volt')+Y.*sin(AngleIJ')*diag(Volt);
|
||
H=H;
|
||
N=N;
|
||
J=J;
|
||
L=L;
|
||
t1=zeros(2*Busnum);
|
||
% t1(1:2:2*Busnum,1:2:2*Busnum)=-H;%10111227
|
||
% t1(1:2:2*Busnum,2:2:2*Busnum)=N;
|
||
% t1(2:2:2*Busnum,1:2:2*Busnum)=-J;%10111227
|
||
% t1(2:2:2*Busnum,2:2:2*Busnum)=L;
|
||
%暂时改一下
|
||
t1=[J,L;
|
||
H,N;
|
||
]';
|
||
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 |