You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
202 lines
8.0 KiB
202 lines
8.0 KiB
function [jac,pIdx,qIdx,vIdx,sIdx,busType]=getJacobianMod(nbus,nline,bus,sw,pv,pq,shunt,line,zip,ind,s,V,dyn,Ytr,Ysh)
|
|
% Generate jacobian matrix
|
|
%
|
|
% FUNCTION getJacobianMod
|
|
%
|
|
% Author: Rui Yao <ruiyao@ieee.org>
|
|
%
|
|
% Copyright (C) 2021, UChicago Argonne, LLC. All rights reserved.
|
|
%
|
|
% OPEN SOURCE LICENSE
|
|
%
|
|
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
|
|
%
|
|
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
|
|
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
|
|
% 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
|
|
%
|
|
%
|
|
% ******************************************************************************************************
|
|
% DISCLAIMER
|
|
%
|
|
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
|
|
% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
|
|
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
|
|
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
|
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
|
|
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
|
|
% OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
% ***************************************************************************************************
|
|
%
|
|
|
|
[~,Ytrx,Yshx,ytrfr,ytrto,yshfr,yshto]=getYMatrix(nbus,line);
|
|
|
|
if nargin<14
|
|
Ytr=Ytrx;
|
|
end
|
|
if nargin<15
|
|
Ysh=Yshx;
|
|
end
|
|
|
|
yBr=[-ytrfr;-ytrto];
|
|
|
|
busType=zeros(nbus,1);
|
|
if isempty(pv)
|
|
pv=zeros(0,6);
|
|
end
|
|
if isempty(pq)
|
|
pq=zeros(0,6);
|
|
end
|
|
if isempty(shunt)
|
|
shunt=zeros(0,7);
|
|
end
|
|
busType(pv(:,1))=1;
|
|
busType(sw(:,1))=2;
|
|
|
|
yShunt=zeros(nbus,1);
|
|
yShunt(shunt(:,1))=shunt(:,5)+1j*shunt(:,6);
|
|
Ysh=Ysh+yShunt;
|
|
if ~isempty(zip)%zipMode=0
|
|
Ysh=Ysh+accumarray(zip(:,1),(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]);
|
|
end
|
|
Y=Ytr+sparse(1:nbus,1:nbus,Ysh,nbus,nbus);
|
|
|
|
if ~isempty(ind)
|
|
nInd=size(ind,1);
|
|
indIdx=ind(:,1);
|
|
R1=ind(:,7);
|
|
X1=ind(:,8);
|
|
Z1=ind(:,7)+1j*ind(:,8);
|
|
Xm=ind(:,13);
|
|
Ze=1j*ind(:,13);
|
|
R2=ind(:,9);
|
|
X2=ind(:,10);
|
|
|
|
Zind=Z1+Ze.*(R2+1j*X2.*s)./(R2+(1j*X2+Ze).*s);
|
|
|
|
Y=Y+sparse(indIdx,indIdx,1./Zind,nbus,nbus);
|
|
end
|
|
|
|
isw=find(busType==2);
|
|
ipv=find(busType==1);
|
|
ipq=find(busType==0);
|
|
busTypePvPq=busType(busType~=2);
|
|
ipvInPvPq=find(busTypePvPq==1);
|
|
npq=size(ipq,1);
|
|
npv=size(ipv,1);
|
|
|
|
i_fr=line(:,1);
|
|
i_to=line(:,2);
|
|
yDiag=diag(Y);
|
|
|
|
C=real(V);
|
|
D=imag(V);
|
|
vMag=abs(V);
|
|
|
|
G=real(Y);
|
|
B=imag(Y);
|
|
|
|
gcbdx=G*C-B*D;
|
|
gcbd=sparse(1:nbus,1:nbus,gcbdx,nbus,nbus);
|
|
bcgdx=B*C+G*D;
|
|
bcgd=sparse(1:nbus,1:nbus,bcgdx,nbus,nbus);
|
|
diagc=sparse(1:nbus,1:nbus,C,nbus,nbus);
|
|
diagd=sparse(1:nbus,1:nbus,D,nbus,nbus);
|
|
cg=diagc*G;
|
|
cb=diagc*B;
|
|
dg=diagd*G;
|
|
db=diagd*B;
|
|
|
|
Hmat=gcbd+cg+db; % P/C
|
|
Nmat=bcgd-cb+dg; % P/D
|
|
Jmat=-bcgd-cb+dg; % Q/C
|
|
Lmat=gcbd-cg-db; % Q/D
|
|
|
|
if ~isempty(zip)
|
|
Czip=C(zip(:,1));
|
|
Dzip=D(zip(:,1));
|
|
vMagzip=vMag(zip(:,1));
|
|
Hmat=Hmat+sparse(zip(:,1),zip(:,1),zip(:,6).*Czip./vMagzip,nbus,nbus);
|
|
Nmat=Nmat+sparse(zip(:,1),zip(:,1),zip(:,6).*Dzip./vMagzip,nbus,nbus);
|
|
Jmat=Jmat+sparse(zip(:,1),zip(:,1),zip(:,9).*Czip./vMagzip,nbus,nbus);
|
|
Lmat=Lmat+sparse(zip(:,1),zip(:,1),zip(:,9).*Dzip./vMagzip,nbus,nbus);
|
|
end
|
|
|
|
% Hdiag=sparse(ones(2*nline+nbus,1),[i_to;i_fr;(1:nbus)'],[-(real(yBr).*C([i_to;i_fr])-imag(yBr).*D([i_to;i_fr]));-2*real(yDiag).*C],1,nbus);
|
|
% Ndiag=sparse(ones(2*nline+nbus,1),[i_to;i_fr;(1:nbus)'],[-(real(yBr).*D([i_to;i_fr])+imag(yBr).*C([i_to;i_fr]));-2*real(yDiag).*D],1,nbus);
|
|
% Jdiag=sparse(ones(2*nline+nbus,1),[i_to;i_fr;(1:nbus)'],[(real(yBr).*D([i_to;i_fr])+imag(yBr).*C([i_to;i_fr]));2*imag(yDiag).*C],1,nbus);
|
|
% Ldiag=sparse(ones(2*nline+nbus,1),[i_to;i_fr;(1:nbus)'],[-(real(yBr).*C([i_to;i_fr])-imag(yBr).*D([i_to;i_fr]));2*imag(yDiag).*D],1,nbus);
|
|
%
|
|
% Hmat=sparse([i_fr;i_to;(1:nbus)'],[i_to;i_fr;(1:nbus)'],[-(real(yBr).*C([i_fr;i_to])+imag(yBr).*D([i_fr;i_to]));Hdiag.'],nbus,nbus);
|
|
% Nmat=sparse([i_fr;i_to;(1:nbus)'],[i_to;i_fr;(1:nbus)'],[(imag(yBr).*C([i_fr;i_to])-real(yBr).*D([i_fr;i_to]));Ndiag.'],nbus,nbus);
|
|
% Jmat=sparse([i_fr;i_to;(1:nbus)'],[i_to;i_fr;(1:nbus)'],[(imag(yBr).*C([i_fr;i_to])-real(yBr).*D([i_fr;i_to]));Jdiag.'],nbus,nbus);
|
|
% Lmat=sparse([i_fr;i_to;(1:nbus)'],[i_to;i_fr;(1:nbus)'],[(real(yBr).*C([i_fr;i_to])+imag(yBr).*D([i_fr;i_to]));Ldiag.'],nbus,nbus);
|
|
|
|
% VdC=sparse(1:nbus,1:nbus,-2*C,nbus,nbus);
|
|
% VdD=sparse(1:nbus,1:nbus,-2*D,nbus,nbus);
|
|
C0Mx=sparse(1:npv,ipvInPvPq,-2*C(ipv),npv,npv+npq);
|
|
D0Mx=sparse(1:npv,ipvInPvPq,-2*D(ipv),npv,npv+npq);
|
|
|
|
jac=[Hmat(busType~=2,busType~=2),Nmat(busType~=2,busType~=2);...
|
|
Jmat(busType==0,busType~=2),Lmat(busType==0,busType~=2);...
|
|
C0Mx,D0Mx];
|
|
|
|
pIdx=1:(npv+npq);
|
|
qIdx=npv+npq+(1:npq);
|
|
vIdx=npv+npq+npq+(1:npv);
|
|
sIdx=[];
|
|
|
|
ind=ind(dyn==0,:);
|
|
nInd=size(ind,1);
|
|
if nInd>1
|
|
indIdx=ind(:,1);
|
|
s=s(dyn==0);
|
|
R1=ind(:,7);
|
|
X1=ind(:,8);
|
|
Z1=ind(:,7)+1j*ind(:,8);
|
|
Xm=ind(:,13);
|
|
Ze=1j*ind(:,13);
|
|
R2=ind(:,9);
|
|
X2=ind(:,10);
|
|
T0=ind(:,15)+ind(:,16)+ind(:,17);
|
|
T1=-ind(:,16)-2*ind(:,17);
|
|
T2=ind(:,17);
|
|
H=ind(:,14);
|
|
|
|
Vind=V(indIdx);
|
|
C1=real(Vind);
|
|
C2=imag(Vind);
|
|
|
|
y2=-(R2.*Xm.*Xm)./(s.*(R2./s - X2*1i + (Xm.*(R1 - X1*1i)*1i)./...
|
|
(- R1 + X1*1i + Xm*1i)).*(R1 + X1*1i + Xm*1i).*(X2*1i + R2./s + (Xm.*(R1 + X1*1i)*1i)./(R1 + X1*1i + Xm*1i)).*(- R1 + X1*1i + Xm*1i)); % Te=V^2.*y2
|
|
y2ds=(R2.*Xm.*Xm)./(s.*s.*(R2./s - X2*1i + (Xm.*(R1 - X1*1i)*1i)./(- R1 + X1*1i + Xm*1i)).*(R1 + X1*1i + Xm*1i).*(X2*1i + R2./s + (Xm.*(R1 + X1*1i)*1i)./(R1 + X1*1i + Xm*1i))...
|
|
.*(- R1 + X1*1i + Xm*1i)) - (R2.*R2.*Xm.*Xm)./(s.*s.*s.*(R2./s - X2*1i + (Xm.*(R1 - X1*1i)*1i)./(- R1 + X1*1i + Xm*1i)).*(R1 + X1*1i + Xm*1i)...
|
|
.*(X2*1i + R2./s + (Xm.*(R1 + X1*1i)*1i)./(R1 + X1*1i + Xm*1i)).^2.*(- R1 + X1*1i + Xm*1i)) - (R2.*Xm.*R2.*Xm)./(s.*s.*s.*(R2./s - X2*1i + (Xm.*(R1 - X1*1i)*1i)...
|
|
./(- R1 + X1*1i + Xm*1i)).^2.*(R1 + X1*1i + Xm*1i).*(X2*1i + R2./s + (Xm.*(R1 + X1*1i)*1i)./(R1 + X1*1i + Xm*1i)).*(- R1 + X1*1i + Xm*1i));
|
|
dTds=(T1+2*T2.*s-abs(Vind).*abs(Vind).*y2ds)/2./H;
|
|
dTdC1=-2*C1.*real(y2)/2./H;
|
|
dTdC2=-2*C2.*real(y2)/2./H;
|
|
|
|
yIL=(R2 + s.*(X2 + Xm)*1i)./(R1.*R2 - s.*(X1.*X2 + X1.*Xm + X2.*Xm) + R2.*(X1 + Xm)*1i + R1.*s.*(X2 + Xm)*1i);
|
|
yILds=(X2*1i + Xm*1i)./(R1.*R2 - s.*(X1.*X2 + X1.*Xm + X2.*Xm) + R2.*(X1 + Xm)*1i + R1.*s.*(X2 + Xm)*1i) ...
|
|
+ (R2 + s.*(X2 + Xm)*1i).*(X1.*X2 + X1.*Xm + X2.*Xm - R1.*(X2 + Xm)*1i)./(R1.*R2 - s.*(X1.*X2 + X1.*Xm + X2.*Xm) + R2.*(X1 + Xm)*1i + R1.*s.*(X2 + Xm)*1i).^2;
|
|
% VILdC1=-2*C1.*conj(yIL);
|
|
% VILdC2=-2*C2.*conj(yIL);
|
|
% Hmat=Hmat+sparse(indIdx,indIdx,real(VILdC1),nbus,nbus);
|
|
% Nmat=Nmat+sparse(indIdx,indIdx,real(VILdC2),nbus,nbus);
|
|
% Jmat=Jmat+sparse(indIdx,indIdx,imag(VILdC1),nbus,nbus);
|
|
% Lmat=Lmat+sparse(indIdx,indIdx,imag(VILdC2),nbus,nbus);
|
|
|
|
VILds=abs(Vind).*abs(Vind).*conj(yILds);
|
|
PILdsmat=sparse(indIdx,1:nInd,-real(VILds),nbus,nInd);
|
|
QILdsmat=sparse(indIdx,1:nInd,-imag(VILds),nbus,nInd);
|
|
dTdC1mat=sparse(1:nInd,indIdx,dTdC1,nInd,nbus);
|
|
dTdC2mat=sparse(1:nInd,indIdx,dTdC2,nInd,nbus);
|
|
ind2ind=sparse(1:nInd,1:nInd,dTds,nInd,nInd);
|
|
|
|
jac=[jac,[PILdsmat(busType~=2,:);QILdsmat(busType==0,:);sparse(npv,nInd)]...
|
|
dTdC1mat(:,busType~=2),dTdC2mat(:,busType~=2),ind2ind];
|
|
sIdx=vIdx(end)+(1:nInd);
|
|
end
|
|
|
|
end |