ANL-CEEESA
/
powersas.m
mirror of https://github.com/ANL-CEEESA/powersas.m
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
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.
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'release'
${ noResults }
powersas.m/internal/hemMachinePFSalientcontinue...

782 lines
35 KiB
Raw Permalink Blame History

function [V,W,Q,s,d,JG,KG]=hemMachinePFSalientcontinueSep(SimData,SysData,SysPara,x,pShare,Vsp)
% HE computation of power flow (extended including synchronous generators)
%
% FUNCTION hemMachinePFSalientcontinueSep
%
% 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.
% ***************************************************************************************************
%
% INPUT (will be consolidated in a future version)
%
% OUTPUT (will be consolidated in a future version)
%
global IS_OCTAVE;
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
[V0,Q0,s0,d0,w,eq1,eq2,ed1,ed2,psid,psiq,Pm0,Ef0,Vavrm,Vavrr,Vavrf,Vavrref,tgovg,tgovm,tgovmech,f,dpg,qplt,vg]=unfoldX(x,SysData);
[pqIncr,pvIncr,Rind0,Rind1,Reind0,Reind1,Rzip0,Rzip1,Ytr,~,Ysh0,Ysh1,VspSq2,~,~,~,~,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara);
[maxAlpha,segAlpha,dAlpha,nlvl,taylorN,alphaTol,diffTol,diffTolMax,method]=unfoldSimData(SimData);
nbus=size(bus,1);
nline=size(line,1);
if isfield(SysPara,'nIslands')&&isfield(SysPara,'islands')&&isfield(SysPara,'refs')
nIslands=SysPara.nIslands;islands=SysPara.islands;refs=SysPara.refs;
else
[nIslands,islands,refs]=searchIslands(bus(:,1),line(:,1:2));
end
Pls=zeros(nbus,1);Pls(pq(:,1))=pqIncr(:,1);if ~isempty(pv);Pls(pv(:,1))=Pls(pv(:,1))-pvIncr;end
Qls=zeros(nbus,1);Qls(pq(:,1))=pqIncr(:,2);
[Y,Ytr,Ysh,ytrfr,ytrto,yshfr,yshto]=getYMatrix(nbus,line);
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
if isempty(sw)
sw=zeros(0,13);
end
busType(pv(:,1))=1;
busType(sw(:,1))=2;
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);
% yShunt=zeros(nbus,1);
% yShunt(shunt(:,1))=shunt(:,5)+1j*shunt(:,6);
% if ~isempty(zip)%zipMode=0
% yShunt=yShunt+accumarray(zip(:,1),(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]);
% end
% Ysh=Ysh+yShunt;
% Y=Y+sparse(1:nbus,1:nbus,yShunt,nbus,nbus);
if ~isempty(zip)%zipMode=0
Ysh0=Ysh0+accumarray(zip(:,1),Rzip0.*(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]);
Ysh1=Ysh1+accumarray(zip(:,1),Rzip1.*(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]);
end
Y=Ytr+sparse(1:nbus,1:nbus,Ysh0,nbus,nbus);
pVec=zeros(nbus,1);
qVec=zeros(nbus,1);
% vSp=zeros(nbus,1);
pVec(pv(:,1))=pVec(pv(:,1))+pv(:,4);
pVec(pq(:,1))=pVec(pq(:,1))-pq(:,4);
qVec(pq(:,1))=qVec(pq(:,1))-pq(:,5);
if ~isempty(zip)%zipMode=0, account for the PQ components in ZIP loads
pVec=pVec-accumarray(zip(:,1),Rzip0.*zip(:,7).*zip(:,12),[nbus,1]);
qVec=qVec-accumarray(zip(:,1),Rzip0.*zip(:,10).*zip(:,12),[nbus,1]);
end
% qVec(ipv)=Q0(ipv);
% vSp(ipv)=pv(:,5);
V=zeros(nbus,nlvl+1);
V(:,1)=V0;
V(isw,2)=Vsp(isw);
W=zeros(nbus,nlvl+1);
W(:,1)=1./V0;
P=zeros(nbus,nlvl+1);
P(:,1)=pVec;
% P(isw,2:end)=0;
Q=zeros(nbus,nlvl+1);
Qxtra=zeros(size(Q));
Q(:,1)=Q0;
Qxtra(:,1)=qVec;
P(:,2:(size(Pls,2)+1))=-Pls;
Qxtra(:,2:(size(Qls,2)+1))=-Qls;
if ~isempty(zip)
P(:,2)=P(:,2)-accumarray(zip(:,1),Rzip1.*zip(:,7).*zip(:,12),[nbus,1]);
Qxtra(:,2)=Qxtra(:,2)-accumarray(zip(:,1),Rzip1.*zip(:,10).*zip(:,12),[nbus,1]);
end
C0=real(V(:,1));
D0=imag(V(:,1));
E0=real(W(:,1));
F0=imag(W(:,1));
Pspec=P(:,1);
Qspec=Q(:,1)+Qxtra(:,1);
V0sq=C0.*C0+D0.*D0;
C0M=sparse(1:nbus,1:nbus,C0,nbus,nbus);
D0M=sparse(1:nbus,1:nbus,D0,nbus,nbus);
E0M=sparse(1:nbus,1:nbus,E0,nbus,nbus);
F0M=sparse(1:nbus,1:nbus,F0,nbus,nbus);
P0M=sparse(1:nbus,1:nbus,Pspec,nbus,nbus);
Q0M=sparse(1:nbus,1:nbus,Qspec,nbus,nbus);
G=real(Y);
B=imag(Y);
if ~isempty(ind)
nInd=size(ind,1);
indIdx=ind(:,1);
s=zeros(nInd,nlvl+1);
s(:,1)=s0;
IL=zeros(nInd,nlvl+1);
IR=zeros(nInd,nlvl+1);
Vm=zeros(nInd,nlvl+1);
R1=ind(:,7);
X1=ind(:,8);
Z1=ind(:,7)+1j*ind(:,8);
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);
Am=sparse(indIdx,(1:nInd)',ones(1,nInd),nbus,nInd);
IL(:,1)=V0(indIdx)./(Z1+Ze.*(R2+1j*X2.*s0)./(R2.*Reind0+(1j*X2.*Reind0+Ze).*s0));
Vm(:,1)=V0(indIdx)-IL(:,1).*Z1;
IR(:,1)=Vm(:,1).*s0./(R2+1j*X2.*s0);
J0=real(IR(:,1));
K0=imag(IR(:,1));
JL0=real(IL(:,1));
KL0=imag(IL(:,1));
Yeind0=Reind0./Ze;
Yeind1=Reind1./Ze;
Ye1ind0=Reind0.*Z1./Ze;
Ye1ind1=Reind1.*Z1./Ze;
Ge=real(Yeind0);
Be=imag(Yeind0);
kg1e=real(Ye1ind0);
kb1e=imag(Ye1ind0);
Ge1=real(Yeind1);
Be1=imag(Yeind1);
kg1e1=real(Ye1ind1);
kb1e1=imag(Ye1ind1);
% The structure:
% |C|D||x| |E|
% |A|B||y|=|F|
% Ly=(D-CA^-1B)y=E-CA^-1F
% x=-A^-1By+A^-1F
% x=[J;K;s];y=[JL;KL;C;D];
% L=[R,S][y1;y2];
% y1=[JL;KL];y2=[C;D];
% y1=-R\S*y2+R\(E-CA^-1F)
%
% idxE=[1,2,5];
% idxM=[3,4,6,7];
%
% LHS_MatInd_Shr=zeros(nInd,2,2);
% RHS_C_Shr=zeros(nInd,2,5);
% LHS_MatInd_Shr2=zeros(nInd,3,4); % A^-1B
% LHS_MatInd_Shr3=zeros(nInd,3,3); % A^-1
%
% LHS_MatInd_Shr_1=zeros(nInd,2,2);
% RHS_C_Shr_1=zeros(nInd,2,5);
% LHS_MatInd_Shr2_1=zeros(nInd,3,4); % A^-1B
% LHS_MatInd_Shr3_1=zeros(nInd,3,3); % A^-1
% for i=1:nInd
% if Rind0(i)~=0
% LHS_MatInd=[R2(i),-X2(i)*s0(i),R1(i)*s0(i),-X1(i)*s0(i),-K0(i)*X2(i)-C0(indIdx(i))+JL0(i)*R1(i)-KL0(i)*X1(i),-s0(i),0;...
% X2(i)*s0(i), R2(i),X1(i)*s0(i), R1(i)*s0(i), J0(i)*X2(i)-D0(indIdx(i))+JL0(i)*X1(i)+KL0(i)*R1(i),0,-s0(i);...
% 2*J0(i),2*K0(i),0,0,-(Rind0(i)*T0(i)+2*Rind0(i)*T1(i)*s0(i)+3*Rind0(i)*T2(i)*s0(i)*s0(i))/R2(i),0,0;...
% -1,0,1+kg1e(i),-kb1e(i),0,-Ge(i), Be(i);...
% 0,-1,kb1e(i), 1+kg1e(i),0,-Be(i),-Ge(i);];
% temp0=LHS_MatInd([3,4,5],idxE)\eye(3); % A^-1
% LHS_MatInd_Shr2(i,:,:)=temp0*LHS_MatInd([3,4,5],idxM); % A^-1B
% LHS_MatInd_Shr3(i,:,:)=temp0; % A^-1
% temp1=LHS_MatInd([1,2],idxE)/LHS_MatInd([3,4,5],idxE); % CA^-1
% temp2=LHS_MatInd([1,2],idxM)-temp1*LHS_MatInd([3,4,5],idxM); % L=D-CA^-1B
% LHS_MatInd_Shr(i,:,:)=-temp2(:,[1,2])\temp2(:,[3,4]); % -R\S
% RHS_C_Shr(i,:,:)=temp2(:,[1,2])\[eye(2),-temp1]; % R\[I,-CA^-1]
%
% LHS_MatInd_Shr_1(i,:,:)=LHS_MatInd_Shr(i,:,:);
% RHS_C_Shr_1(i,:,:)=RHS_C_Shr(i,:,:);
% LHS_MatInd_Shr2_1(i,:,:)=LHS_MatInd_Shr2(i,:,:);
% LHS_MatInd_Shr3_1(i,:,:)=LHS_MatInd_Shr3(i,:,:);
% else % singular case
% LHS_MatInd=[R2(i),-X2(i)*s0(i),R1(i)*s0(i),-X1(i)*s0(i),-K0(i)*X2(i)-C0(indIdx(i))+JL0(i)*R1(i)-KL0(i)*X1(i),-s0(i),0;...
% X2(i)*s0(i), R2(i),X1(i)*s0(i), R1(i)*s0(i), J0(i)*X2(i)-D0(indIdx(i))+JL0(i)*X1(i)+KL0(i)*R1(i),0,-s0(i);...
% 2*real(IR(i,2)),2*imag(IR(i,2)),0,0,-(Rind1(i)*T0(i)+2*Rind1(i)*T1(i)*s0(i)+3*Rind1(i)*T2(i)*s0(i)*s0(i))/R2(i),0,0;...
% -1,0,1+kg1e(i),-kb1e(i),0,-Ge(i), Be(i);...
% 0,-1,kb1e(i), 1+kg1e(i),0,-Be(i),-Ge(i);];
% temp0=LHS_MatInd([3,4,5],idxE)\eye(3); % A^-1
% LHS_MatInd_Shr2(i,:,:)=temp0*LHS_MatInd([3,4,5],idxM); % A^-1B
% LHS_MatInd_Shr3(i,:,:)=temp0; % A^-1
% temp1=LHS_MatInd([1,2],idxE)/LHS_MatInd([3,4,5],idxE); % CA^-1
% temp2=LHS_MatInd([1,2],idxM)-temp1*LHS_MatInd([3,4,5],idxM); % L=D-CA^-1B
% LHS_MatInd_Shr(i,:,:)=-temp2(:,[1,2])\temp2(:,[3,4]); % -R\S
% RHS_C_Shr(i,:,:)=temp2(:,[1,2])\[eye(2),-temp1]; % R\[I,-CA^-1]
%
% temp3=[1+kg1e(i),-kb1e(i);kb1e(i), 1+kg1e(i)];
% temp4=[Ge(i),-Be(i);Be(i),Ge(i)];
% LHS_MatInd_Shr_1(i,:,:)=temp3\temp4; % -R\S (i=1)
% RHS_C_Shr_1(i,:,:)=temp3\[eye(2),zeros(2,3)]; % R\[I,-CA^-1] (i=1)
% end
% end
% LHS_MatInd_Bus=zeros(nbus,2,2); % \sum{-R\S} by buses
% LHS_MatInd_Bus(:,1,1)=accumarray(indIdx,LHS_MatInd_Shr(:,1,1),[nbus,1]);
% LHS_MatInd_Bus(:,1,2)=accumarray(indIdx,LHS_MatInd_Shr(:,1,2),[nbus,1]);
% LHS_MatInd_Bus(:,2,1)=accumarray(indIdx,LHS_MatInd_Shr(:,2,1),[nbus,1]);
% LHS_MatInd_Bus(:,2,2)=accumarray(indIdx,LHS_MatInd_Shr(:,2,2),[nbus,1]);
%
% LHS_MatInd_Bus_1=zeros(nbus,2,2); % \sum{-R\S} by buses (i=1)
% LHS_MatInd_Bus_1(:,1,1)=accumarray(indIdx,LHS_MatInd_Shr_1(:,1,1),[nbus,1]);
% LHS_MatInd_Bus_1(:,1,2)=accumarray(indIdx,LHS_MatInd_Shr_1(:,1,2),[nbus,1]);
% LHS_MatInd_Bus_1(:,2,1)=accumarray(indIdx,LHS_MatInd_Shr_1(:,2,1),[nbus,1]);
% LHS_MatInd_Bus_1(:,2,2)=accumarray(indIdx,LHS_MatInd_Shr_1(:,2,2),[nbus,1]);
LHS_MatInd_Full=zeros(nInd,5,2);
LHS_MatInd_Shr=zeros(nInd,2,2);
% LHS_MatInd_Shr_1=zeros(nInd,2,2);
RHS_C_Shr=zeros(nInd,5,5);
for i=1:nInd
LHS_MatInd=[R2(i),-X2(i)*s0(i),R1(i)*s0(i),-X1(i)*s0(i),-K0(i)*X2(i)-C0(indIdx(i))+JL0(i)*R1(i)-KL0(i)*X1(i),-s0(i),0;...
X2(i)*s0(i), R2(i),X1(i)*s0(i), R1(i)*s0(i), J0(i)*X2(i)-D0(indIdx(i))+JL0(i)*X1(i)+KL0(i)*R1(i),0,-s0(i);...
C0(indIdx(i))-JL0(i)*R1(i)+KL0(i)*X1(i),D0(indIdx(i))-KL0(i)*R1(i)-JL0(i)*X1(i),-J0(i)*R1(i)-K0(i)*X1(i),-K0(i)*R1(i)+J0(i)*X1(i),-Rind0(i)*(T1(i)+2*T2(i)*s0(i)),J0(i),K0(i);...
-1,0,1+kg1e(i),-kb1e(i),0,-Ge(i), Be(i);...
0,-1,kb1e(i), 1+kg1e(i),0,-Be(i),-Ge(i);];
MatIndA=LHS_MatInd(:,1:5);
MatIndB=LHS_MatInd(:,[6,7]);
MatInvA=MatIndA\eye(5);
MatCDtoRest=-MatInvA*MatIndB;
RHS_C_Shr(i,:,:)=MatInvA;
LHS_MatInd_Full(i,:,:)=MatCDtoRest;
LHS_MatInd_Shr(i,:,:)=MatCDtoRest([3,4],:);
% LHS_MatInd_Shr_1(i,:,:)=LHS_MatInd_Shr(i,:,:);
% temp0=LHS_MatInd([3,4,5],idxE)\eye(3); % A^-1
% LHS_MatInd_Shr2(i,:,:)=temp0*LHS_MatInd([3,4,5],idxM); % A^-1B
% LHS_MatInd_Shr3(i,:,:)=temp0; % A^-1
% temp1=LHS_MatInd([1,2],idxE)/LHS_MatInd([3,4,5],idxE); % CA^-1
% temp2=LHS_MatInd([1,2],idxM)-temp1*LHS_MatInd([3,4,5],idxM); % L=D-CA^-1B
% LHS_MatInd_Shr(i,:,:)=-temp2(:,[1,2])\temp2(:,[3,4]); % -R\S
% RHS_C_Shr(i,:,:)=temp2(:,[1,2])\[eye(2),-temp1]; % R\[I,-CA^-1]
%
% LHS_MatInd_Shr_1(i,:,:)=LHS_MatInd_Shr(i,:,:);
% RHS_C_Shr_1(i,:,:)=RHS_C_Shr(i,:,:);
% LHS_MatInd_Shr2_1(i,:,:)=LHS_MatInd_Shr2(i,:,:);
% LHS_MatInd_Shr3_1(i,:,:)=LHS_MatInd_Shr3(i,:,:);
end
LHS_MatInd_Bus=zeros(nbus,2,2); % \sum{-R\S} by buses
LHS_MatInd_Bus(:,1,1)=accumarray(indIdx,LHS_MatInd_Shr(:,1,1),[nbus,1]);
LHS_MatInd_Bus(:,1,2)=accumarray(indIdx,LHS_MatInd_Shr(:,1,2),[nbus,1]);
LHS_MatInd_Bus(:,2,1)=accumarray(indIdx,LHS_MatInd_Shr(:,2,1),[nbus,1]);
LHS_MatInd_Bus(:,2,2)=accumarray(indIdx,LHS_MatInd_Shr(:,2,2),[nbus,1]);
% LHS_MatInd_Bus_1=zeros(nbus,2,2); % \sum{-R\S} by buses (i=1)
% LHS_MatInd_Bus_1(:,1,1)=accumarray(indIdx,LHS_MatInd_Shr_1(:,1,1),[nbus,1]);
% LHS_MatInd_Bus_1(:,1,2)=accumarray(indIdx,LHS_MatInd_Shr_1(:,1,2),[nbus,1]);
% LHS_MatInd_Bus_1(:,2,1)=accumarray(indIdx,LHS_MatInd_Shr_1(:,2,1),[nbus,1]);
% LHS_MatInd_Bus_1(:,2,2)=accumarray(indIdx,LHS_MatInd_Shr_1(:,2,2),[nbus,1]);
else
s=[];
end
if ~isempty(zip)
nZip=size(zip,1);
zipIdx=zip(:,1);
IiL=zeros(nZip,nlvl+1);
BiL=zeros(nZip,nlvl+1);
Bi0=abs(V0(zipIdx));
JI=zip(:,6);
KI=-zip(:,9);
Ii0L=Rzip0.*(JI+1j*KI).*V0(zipIdx)./Bi0;
Ji0L=real(Ii0L);
Ki0L=imag(Ii0L);
IiL(:,1)=Ii0L;
BiL(:,1)=Bi0;
Ci0=real(V0(zipIdx));
Di0=imag(V0(zipIdx));
LHS_MatZip=[Rzip0.*JI./Bi0-Ci0.*Ji0L./Bi0./Bi0,-Rzip0.*KI./Bi0-Di0.*Ji0L./Bi0./Bi0,...
Rzip0.*KI./Bi0-Ci0.*Ki0L./Bi0./Bi0,Rzip0.*JI./Bi0-Di0.*Ki0L./Bi0./Bi0];
Mat_BZip=[Ci0./Bi0,Di0./Bi0];
else
IiL=[];
end
if ~isempty(syn)
nSyn=size(syn,1);
synIdx=syn(:,1);
Rs=syn(:,7);
Xd=syn(:,8);
Xq=syn(:,13);
Ms=syn(:,18);
Ds=syn(:,19);
d=zeros(nSyn,nlvl+1);
JG=zeros(nSyn,nlvl+1);
KG=zeros(nSyn,nlvl+1);
Cd=zeros(nSyn,nlvl+1);
Sd=zeros(nSyn,nlvl+1);
Ef=zeros(nSyn,nlvl+1);
d(:,1)=d0;
Efd=zeros(nSyn,1);
Efq=Ef0;
cosd=cos(d0);
sind=sin(d0);
Cg=C0(synIdx);
Dg=D0(synIdx);
Vd=sind.*Cg-cosd.*Dg;
Vq=cosd.*Cg+sind.*Dg;
Id=(Rs.*(Efd-Vd)+Xq.*(Efq-Vq))./(Rs.*Rs+Xq.*Xd);
Iq=(-Xd.*(Efd-Vd)+Rs.*(Efq-Vq))./(Rs.*Rs+Xq.*Xd);
IG0=(sind.*Id+cosd.*Iq)+1j*(-cosd.*Id+sind.*Iq);
% IG0=(Ef0.*(cos(d0)+1j*sin(d0))-V0(synIdx))./(Rs+1j*Xs);
JG(:,1)=real(IG0);
KG(:,1)=imag(IG0);
Cd(:,1)=cos(d0);
Sd(:,1)=sin(d0);
Ef(:,[1,2])=[Ef0,Ef1];
CG0=C0(synIdx);
DG0=D0(synIdx);
JG0=JG(:,1);
KG0=KG(:,1);
Cd0=Cd(:,1);
Sd0=Sd(:,1);
[cosp,sinp,taylorN]=getTaylorPolynomials(d0,taylorN); % taylorN may be truncated
Pm=zeros(nSyn,nlvl+1);
Pm(:,1)=Pm0;
if ~isempty(Pm1)
Pm(:,2)=Pm1;
end
A1n=zeros(nSyn,1);
B1n=zeros(nSyn,1);
for i=1:taylorN
A1n=A1n.*d0+(taylorN+1-i)*cosp(:,taylorN+2-i);
B1n=B1n.*d0+(taylorN+1-i)*sinp(:,taylorN+2-i);
end
synByIsland=islands(synIdx);
sumShare=accumarray(synByIsland,pShare);
nIsland=length(sumShare);
pShare=pShare./sumShare(synByIsland);
mAuxIsland=sparse(1:nSyn,synByIsland,pShare,nSyn,nIsland);
[~,idxBal]=max(mAuxIsland);
idxBalSyn=idxBal(synByIsland);
MatG1C=sparse(1:nSyn,synIdx,Cd0,nSyn,nbus);
MatG1D=sparse(1:nSyn,synIdx,Sd0,nSyn,nbus);
MatG2C=sparse(1:nSyn,synIdx,Sd0,nSyn,nbus);
MatG2D=sparse(1:nSyn,synIdx,-Cd0,nSyn,nbus);
MatG5AC=sparse([1:nSyn,1:nSyn],[synIdx;synIdx(idxBalSyn)],[-pShare(idxBalSyn).*JG0;JG0(idxBalSyn).*pShare],nSyn,nbus);
MatG5AD=sparse([1:nSyn,1:nSyn],[synIdx;synIdx(idxBalSyn)],[-pShare(idxBalSyn).*KG0;KG0(idxBalSyn).*pShare],nSyn,nbus);
MatG5BC=sparse(nIsland,nbus);
MatG5BD=sparse(nIsland,nbus);
MatG5AC(idxBal,:)=MatG5BC;
MatG5AD(idxBal,:)=MatG5BD;
MatG1J=sparse(1:nSyn,1:nSyn,Rs.*Cd0+Xd.*Sd0,nSyn,nSyn);
MatG1K=sparse(1:nSyn,1:nSyn,Rs.*Sd0-Xd.*Cd0,nSyn,nSyn);
MatG1d=sparse(1:nSyn,1:nSyn,(CG0+Rs.*JG0-Xd.*KG0).*A1n+(DG0+Rs.*KG0+Xd.*JG0).*B1n,nSyn,nSyn);
MatG2J=sparse(1:nSyn,1:nSyn,Rs.*Sd0-Xq.*Cd0,nSyn,nSyn);
MatG2K=sparse(1:nSyn,1:nSyn,-Rs.*Cd0-Xq.*Sd0,nSyn,nSyn);
MatG2d=sparse(1:nSyn,1:nSyn,(-DG0-Rs.*KG0-Xq.*JG0).*A1n+(CG0+Rs.*JG0-Xq.*KG0).*B1n,nSyn,nSyn);
MatG5AJ=sparse([1:nSyn,1:nSyn],[1:nSyn,idxBalSyn'],[-pShare(idxBalSyn).*(CG0+2*JG0.*Rs);(CG0(idxBalSyn)+2*JG0(idxBalSyn).*Rs(idxBalSyn)).*pShare],nSyn,nSyn);
MatG5AK=sparse([1:nSyn,1:nSyn],[1:nSyn,idxBalSyn'],[-pShare(idxBalSyn).*(DG0+2*KG0.*Rs);(DG0(idxBalSyn)+2*KG0(idxBalSyn).*Rs(idxBalSyn)).*pShare],nSyn,nSyn);
MatG5Ad=sparse(nSyn,nSyn);
MatG5BJ=sparse(nIsland,nSyn);
MatG5BK=sparse(nIsland,nSyn);
MatG5Bd=sparse(synByIsland,1:nSyn,Ms,nIsland,nSyn);
MatG5AJ(idxBal,:)=MatG5BJ;
MatG5AK(idxBal,:)=MatG5BK;
MatG5Ad(idxBal,:)=MatG5Bd;
% [MatGA,MatGB]*[CD;JKd]=RHS
% JKd=-MatGA\MatGB*CD+MatGA\RHS=MatGBiA*CD+MatGA\RHS
MatGA=[MatG1C,MatG1D;MatG2C,MatG2D;MatG5AC,MatG5AD];
MatGB=[MatG1J,MatG1K,MatG1d;MatG2J,MatG2K,MatG2d;MatG5AJ,MatG5AK,MatG5Ad];
% MatGB2=[MatG1J,MatG1K,sparse(1:nSyn,1:nSyn,(CG0+Rs.*JG0-Xd.*KG0),nSyn,nSyn),sparse(1:nSyn,1:nSyn,(DG0+Rs.*KG0+Xd.*JG0),nSyn,nSyn),sparse(nSyn,nSyn);...
% MatG2J,MatG2K,sparse(1:nSyn,1:nSyn,(-DG0-Rs.*KG0-Xq.*JG0),nSyn,nSyn),sparse(1:nSyn,1:nSyn,(CG0+Rs.*JG0-Xq.*KG0),nSyn,nSyn),sparse(nSyn,nSyn);...
% sparse(nSyn,nSyn),sparse(nSyn,nSyn),sparse(1:nSyn,1:nSyn,ones(1,nSyn),nSyn,nSyn),sparse(nSyn,nSyn),sparse(1:nSyn,1:nSyn,-A1n,nSyn,nSyn);...
% sparse(nSyn,nSyn),sparse(nSyn,nSyn),sparse(nSyn,nSyn),sparse(1:nSyn,1:nSyn,ones(1,nSyn),nSyn,nSyn),sparse(1:nSyn,1:nSyn,-B1n,nSyn,nSyn);...
% MatG5AJ,MatG5AK,sparse(nSyn,nSyn),sparse(nSyn,nSyn),MatG5Ad];
MatGBiA=-MatGB\MatGA;
GTrMat=sparse(synIdx,1:nSyn,ones(nSyn,1),nbus,nSyn);
MatGTrans=[GTrMat,sparse(nbus,nSyn),sparse(nbus,nSyn);sparse(nbus,nSyn),GTrMat,sparse(nbus,nSyn)];
MatGCD=MatGTrans*MatGBiA;
else
d=[];
JG=[];
KG=[];
end
Y11=-G;Y12=B;Y21=-B;Y22=-G;
YEF11=P0M;YEF12=-Q0M;YEF21=-Q0M;YEF22=-P0M;
if ~isempty(zip)
Y11=Y11-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,1),[nbus,1]),nbus,nbus);
Y12=Y12-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,2),[nbus,1]),nbus,nbus);
Y21=Y21-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,3),[nbus,1]),nbus,nbus);
Y22=Y22-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,4),[nbus,1]),nbus,nbus);
end
YLHS=[Y11,Y12;Y21,Y22];
% YLHS_1=[Y11,Y12;Y21,Y22];
if ~isempty(ind)
YLHS=YLHS-...
[sparse(1:nbus,1:nbus,LHS_MatInd_Bus(:,1,1),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus(:,1,2),nbus,nbus);...
sparse(1:nbus,1:nbus,LHS_MatInd_Bus(:,2,1),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus(:,2,2),nbus,nbus)];
% YLHS_1=YLHS_1-...
% [sparse(1:nbus,1:nbus,LHS_MatInd_Bus_1(:,1,1),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_1(:,1,2),nbus,nbus);...
% sparse(1:nbus,1:nbus,LHS_MatInd_Bus_1(:,2,1),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_1(:,2,2),nbus,nbus)];
end
if ~isempty(syn)
YLHS=YLHS+MatGCD;
% YLHS_1=YLHS_1+MatGCD;
end
% Branch 0: smaller matrix
C0i=C0./V0sq;
D0i=D0./V0sq;
% C0Mi=sparse(1:nbus,1:nbus,C0./V0sq,nbus,nbus);
% D0Mi=sparse(1:nbus,1:nbus,D0./V0sq,nbus,nbus);
PCQD=Pspec.*C0i+Qspec.*D0i;
PDQC=Pspec.*D0i-Qspec.*C0i;
YLHSx=YLHS-...
[sparse(1:nbus,1:nbus,PCQD.*E0+PDQC.*F0,nbus,nbus),sparse(1:nbus,1:nbus,-PCQD.*F0+PDQC.*E0,nbus,nbus);...
sparse(1:nbus,1:nbus, PDQC.*E0-PCQD.*F0,nbus,nbus),sparse(1:nbus,1:nbus,-PDQC.*F0-PCQD.*E0,nbus,nbus)];
idxNonSw=find(busType~=2);
idxStackMat=[idxNonSw;idxNonSw+nbus];
C0Mx=sparse(1:npv,ipvInPvPq,C0(ipv),npv,npv+npq);
D0Mx=sparse(1:npv,ipvInPvPq,D0(ipv),npv,npv+npq);
LHS_matx=[YLHSx(idxStackMat,idxStackMat),[-F0M(busType~=2,ipv);-E0M(busType~=2,ipv)];...
C0Mx,D0Mx,sparse(npv,npv)];
useLU=isfield(SysPara,'iden')&&isfield(SysPara,'p_amd');
if useLU
if isempty(SysPara.p_amd)
p_amd = colamd (LHS_matx) ;
save([SysPara.iden,'.mat'],'p_amd');
else
p_amd=SysPara.p_amd;
end
MxI = speye (size(LHS_matx)) ;
MxQ = MxI (:, p_amd) ;
if IS_OCTAVE
[MxL,MxU,MxP,MxQx] = lu (LHS_matx*MxQ) ;
else
[MxL,MxU,MxP] = lu (LHS_matx*MxQ) ;
end
end
% END OF Branch 0: smaller matrix
% Branch 1: larger matrix
% idxNonSw=find(busType~=2);
% idxStackMat=[idxNonSw;idxNonSw+nbus];
%
% LHS_mat=[YLHS(idxStackMat,idxStackMat),...
% [YEF11(busType~=2,busType~=2),YEF12(busType~=2,busType~=2),-F0M(busType~=2,ipv);...
% YEF21(busType~=2,busType~=2),YEF22(busType~=2,busType~=2),-E0M(busType~=2,ipv)];...
% C0M(ipv,busType~=2),D0M(ipv,busType~=2),sparse(npv,2*npq+3*npv);...
% E0M(busType~=2,busType~=2),-F0M(busType~=2,busType~=2),C0M(busType~=2,busType~=2),-D0M(busType~=2,busType~=2),sparse(npq+npv,npv);...
% F0M(busType~=2,busType~=2),E0M(busType~=2,busType~=2),D0M(busType~=2,busType~=2),C0M(busType~=2,busType~=2),sparse(npq+npv,npv);];
%
% useLU=isfield(SysPara,'iden')&&isfield(SysPara,'p_amd');
%
% if useLU
% if isempty(SysPara.p_amd)
% p_amd = colamd (LHS_mat) ;
% save([SysPara.iden,'.mat'],'p_amd');
% else
% p_amd=SysPara.p_amd;
% end
% MxI = speye (size(LHS_mat)) ;
% MxQ = MxI (:, p_amd) ;
% [MxL,MxU,MxP] = lu (LHS_mat*MxQ) ;
% end
% END OF Branch 1: larger matrix
% if nbus<=LU_SIZE_TH
% [L_LHS_mat,U_LHS_mat,p_LHS_mat]=lu(LHS_mat,'vector');
% end
% LHS_mat_1=[YLHS_1,[YEF11,YEF12,-F0M(:,ipv);YEF21,YEF22,-E0M(:,ipv)];...
% C0M(ipv,:),D0M(ipv,:),sparse(npv,2*nbus+npv);...
% E0M,-F0M,C0M,-D0M,sparse(nbus,npv);...
% F0M,E0M,D0M,C0M,sparse(nbus,npv);];
for i=1:nlvl
seq2=getseq(i,2);
seq2p=getseq(i+1,2);
seq3=getseq(i,3);
idxSeq2=sum(seq2==i,2);
idxSeq2p=sum(seq2p>=i,2);
idxSeq3=sum(seq3==i,2);
idxSeq3x=sum(seq3(:,[2,3])==i,2);
seq2R=seq2(idxSeq2==0,:);
seq2px=seq2(idxSeq2p==0,:);
seq3R=seq3(idxSeq3==0,:);
seq3Rx=seq3(idxSeq3x==0,:);
seq2m=getseq(i-1,2);
seq3m=getseq(i-1,3);
seq2mm=getseq(i-2,2);
RHSILr=zeros(nbus,1);
RHSILi=zeros(nbus,1);
if ~isempty(ind)
% rhsBus=zeros(2,nInd);
% rhsM=sum(Vm(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)-1j*X2.*sum(IR(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2);
% rhsI=-real(sum(IR(:,seq2R(:,1)+1).*conj(IR(:,seq2R(:,2)+1)),2))+...
% (T1.*sum(s(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)+T2.*sum(s(:,seq3R(:,1)+1).*s(:,seq3R(:,2)+1).*s(:,seq3R(:,3)+1),2)).*Rind0./R2+...
% (T0.*s(:,i)+T1.*sum(s(:,seq2m(:,1)+1).*s(:,seq2m(:,2)+1),2)+T2.*sum(s(:,seq3m(:,1)+1).*s(:,seq3m(:,2)+1).*s(:,seq3m(:,3)+1),2)).*Rind1./R2;
% rhsIx=-real(sum(IR(:,seq2px(:,1)+1).*conj(IR(:,seq2px(:,2)+1)),2))+...
% (T1.*sum(s(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)+T2.*sum(s(:,seq3R(:,1)+1).*s(:,seq3R(:,2)+1).*s(:,seq3R(:,3)+1),2)).*Rind1./R2;
% rhsImod=Rind1.*(T1.*s(:,i)+T2.*sum(s(:,seq2m(:,1)+1).*s(:,seq2m(:,2)+1),2))+Rind0.*T2.*sum(s(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)-...
% real(sum(V(indIdx,seq2R(:,1)+1).*conj(IR(:,seq2R(:,2)+1)),2))+...
% real(sum(IL(:,seq2R(:,1)+1).*conj(IR(:,seq2R(:,2)+1)),2).*Z1);
% if i==1
% rhsImod=rhsImod+Rind1.*T0;
% end
% rhsIL=V(indIdx,i).*Yeind1-IL(:,i).*Ye1ind1;
% for j=1:nInd
% % if Rind0(j)~=0
% % rhsBus(:,j)=squeeze(RHS_C_Shr(j,:,:))*[real(rhsM(j));imag(rhsM(j));rhsI(j);real(rhsIL(j));imag(rhsIL(j))];
% % else
% % if i==1
% % s(j,i+1)=R2(j)*Rind1(j)*T0(j)/...
% % ((C0(indIdx(j))-R1(j)*JL0(j)+X1(j)*KL0(j))*(C0(indIdx(j))-R1(j)*JL0(j)+X1(j)*KL0(j))+...
% % (D0(indIdx(j))-X1(j)*JL0(j)-R1(j)*KL0(j))*(D0(indIdx(j))-X1(j)*JL0(j)-R1(j)*KL0(j)));
% % IR(j,i+1)=s(j,i+1)/R2(j)*((C0(indIdx(j))-R1(j)*JL0(j)+X1(j)*KL0(j))+1j*(D0(indIdx(j))-X1(j)*JL0(j)-R1(j)*KL0(j)));
% % rhsBus(:,j)=squeeze(RHS_C_Shr_1(j,:,:))*[real(IR(j,i+1)+rhsIL(j));imag(IR(j,i+1)+rhsIL(j));0;0;0];
% % else
% % rhsBus(:,j)=squeeze(RHS_C_Shr(j,:,:))*[real(rhsM(j));imag(rhsM(j));rhsIx(j);real(rhsIL(j));imag(rhsIL(j))];
% % end
% % end
% rhsBus(:,j)=squeeze(RHS_C_Shr(j,:,:))*[real(rhsM(j));imag(rhsM(j));rhsImod(j);real(rhsIL(j));imag(rhsIL(j))];
% end
% %accumulate currents
% RHSILr=accumarray(indIdx,rhsBus(1,:)',[nbus,1]);
% RHSILi=accumarray(indIdx,rhsBus(2,:)',[nbus,1]);
rhsBus=zeros(5,nInd);
rhsM=sum(Vm(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)-1j*X2.*sum(IR(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2);
rhsImod=Rind1.*(T1.*s(:,i)+T2.*sum(s(:,seq2m(:,1)+1).*s(:,seq2m(:,2)+1),2))+Rind0.*T2.*sum(s(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)-...
real(sum(V(indIdx,seq2R(:,1)+1).*conj(IR(:,seq2R(:,2)+1)),2))+...
real(sum(IL(:,seq2R(:,1)+1).*conj(IR(:,seq2R(:,2)+1)),2).*Z1);
if i==1
rhsImod=rhsImod+Rind1.*T0;
end
rhsIL=V(indIdx,i).*Yeind1-IL(:,i).*Ye1ind1;
for j=1:nInd
rhsBus(:,j)=squeeze(RHS_C_Shr(j,:,:))*[real(rhsM(j));imag(rhsM(j));rhsImod(j);real(rhsIL(j));imag(rhsIL(j))];
end
RHSILr=accumarray(indIdx,rhsBus(3,:)',[nbus,1]);
RHSILi=accumarray(indIdx,rhsBus(4,:)',[nbus,1]);
end
RHSIiLr=zeros(nbus,1);
RHSIiLi=zeros(nbus,1);
if ~isempty(zip)
RHS_BZip=(real(sum(V(zipIdx,seq2R(:,1)+1).*conj(V(zipIdx,seq2R(:,2)+1)),2))-sum(BiL(:,seq2R(:,1)+1).*BiL(:,seq2R(:,2)+1),2))./Bi0/2;
RHZ_BIConv=sum(IiL(:,seq2R(:,1)+1).*BiL(:,seq2R(:,2)+1),2);
RHSILr_full=Rzip1.*(JI.*real(V(zipIdx,i))-KI.*imag(V(zipIdx,i)))./Bi0-real(RHZ_BIConv)./Bi0-Ji0L.*RHS_BZip./Bi0;
RHSILi_full=Rzip1.*(KI.*real(V(zipIdx,i))+JI.*imag(V(zipIdx,i)))./Bi0-imag(RHZ_BIConv)./Bi0-Ki0L.*RHS_BZip./Bi0;
RHSIiLr=accumarray(zipIdx,RHSILr_full,[nbus,1]);
RHSIiLi=accumarray(zipIdx,RHSILi_full,[nbus,1]);
end
RHSIGr=zeros(nbus,1);
RHSIGi=zeros(nbus,1);
if ~isempty(syn)
AG0=zeros(nSyn,1);
BG0=zeros(nSyn,1);
if taylorN>=2
AG0=AG0+cosp(:,3).*sum(d(:,seq2R(:,1)+1).*d(:,seq2R(:,2)+1),2);
BG0=BG0+sinp(:,3).*sum(d(:,seq2R(:,1)+1).*d(:,seq2R(:,2)+1),2);
end
if taylorN>=3
AG0=AG0+cosp(:,4).*sum(d(:,seq3R(:,1)+1).*d(:,seq3R(:,2)+1).*d(:,seq3R(:,3)+1),2);
BG0=BG0+sinp(:,4).*sum(d(:,seq3R(:,1)+1).*d(:,seq3R(:,2)+1).*d(:,seq3R(:,3)+1),2);
end
if taylorN>=4
seq4=getseq(i,4);
idxSeq4=sum(seq4==i,2);
seq4R=seq4(idxSeq4==0,:);
AG0=AG0+cosp(:,5).*sum(d(:,seq4R(:,1)+1).*d(:,seq4R(:,2)+1).*d(:,seq4R(:,3)+1).*d(:,seq4R(:,4)+1),2);
BG0=BG0+sinp(:,5).*sum(d(:,seq4R(:,1)+1).*d(:,seq4R(:,2)+1).*d(:,seq4R(:,3)+1).*d(:,seq4R(:,4)+1),2);
end
CCr=sum(real(V(synIdx,seq2R(:,1)+1)).*Cd(:,seq2R(:,2)+1),2);
DCr=sum(imag(V(synIdx,seq2R(:,1)+1)).*Cd(:,seq2R(:,2)+1),2);
CSr=sum(real(V(synIdx,seq2R(:,1)+1)).*Sd(:,seq2R(:,2)+1),2);
DSr=sum(imag(V(synIdx,seq2R(:,1)+1)).*Sd(:,seq2R(:,2)+1),2);
JCr=sum(JG(:,seq2R(:,1)+1).*Cd(:,seq2R(:,2)+1),2);
KCr=sum(KG(:,seq2R(:,1)+1).*Cd(:,seq2R(:,2)+1),2);
JSr=sum(JG(:,seq2R(:,1)+1).*Sd(:,seq2R(:,2)+1),2);
KSr=sum(KG(:,seq2R(:,1)+1).*Sd(:,seq2R(:,2)+1),2);
RHSIG1=Ef(:,i+1)-(CCr+DSr+Rs.*(JCr+KSr)+Xd.*(JSr-KCr))-(CG0+Rs.*JG0-Xd.*KG0).*AG0-(DG0+Rs.*KG0+Xd.*JG0).*BG0;
RHSIG2=-(CSr-DCr+Rs.*(JSr-KCr)-Xq.*(JCr+KSr))-(-DG0-Rs.*KG0-Xq.*JG0).*AG0-(CG0+Rs.*JG0-Xq.*KG0).*BG0;
RHSIG3temp=-Pm(:,i+1)+...
sum(real(V(synIdx,seq2R(:,1)+1).*(JG(:,seq2R(:,2)+1)-1j*KG(:,seq2R(:,2)+1))),2)+...
(sum(JG(:,seq2R(:,1)+1).*JG(:,seq2R(:,2)+1),2)+sum(KG(:,seq2R(:,1)+1).*KG(:,seq2R(:,2)+1),2)).*Rs;
RHSIG3=pShare(idxBalSyn).*RHSIG3temp-RHSIG3temp(idxBalSyn).*pShare;
RHSIG3(idxBal)=0;
RHSIG=MatGB\[RHSIG1;RHSIG2;RHSIG3];
RHSIGJK=MatGTrans*RHSIG;
RHSIGr=RHSIGJK(1:nbus);
RHSIGi=RHSIGJK((nbus+1):(2*nbus));
end
% HEM Body
RHS1=sum((-P(:,seq2(:,1)+1)+1j*(Q(:,seq2(:,1)+1)+Qxtra(:,seq2(:,1)+1))).*conj(W(:,seq2(:,2)+1)),2)+Ysh1.*V(:,i);
RHS2=-0.5*real(sum(V(:,seq2R(:,1)+1).*conj(V(:,seq2R(:,2)+1)),2));
RHS3=sum(-W(:,seq2R(:,1)+1).*V(:,seq2R(:,2)+1),2);
if i==1
RHS2=RHS2+0.5*VspSq2(:,2);
end
compactRHS1=RHS1(busType~=2);
if ~isempty(isw);compactRHS1=compactRHS1+Y(busType~=2,isw)*V(isw,i+1);end
RHS3r=real(RHS3);RHS3i=imag(RHS3);
% RHS=[real(compactRHS1)+RHSILr(busType~=2)+RHSIiLr(busType~=2)-RHSIGr(busType~=2);...
% imag(compactRHS1)+RHSILi(busType~=2)+RHSIiLi(busType~=2)-RHSIGi(busType~=2);...
% RHS2(ipv);...
% RHS3r(busType~=2);...
% RHS3i(busType~=2);];
%
% if useLU
% x = MxQ * (MxU \ (MxL \ (MxP * RHS))) ;
% else
% x=LHS_mat\RHS;
% end
%
% V(busType~=2,i+1)=x(1:(npq+npv))+1j*x(((npq+npv)+1):(2*(npq+npv)));
% W(busType~=2,i+1)=x((2*(npq+npv)+1):(3*(npq+npv)))+1j*x((3*(npq+npv)+1):(4*(npq+npv)));
% Q(ipv,i+1)=x((4*(npq+npv)+1):end);
% Solve a smaller system
RHS3xr=PCQD.*RHS3r+PDQC.*RHS3i;
RHS3xi=PDQC.*RHS3r-PCQD.*RHS3i;
RHSx=[real(compactRHS1)+RHSILr(busType~=2)+RHSIiLr(busType~=2)-RHSIGr(busType~=2)-RHS3xr(busType~=2);...
imag(compactRHS1)+RHSILi(busType~=2)+RHSIiLi(busType~=2)-RHSIGi(busType~=2)-RHS3xi(busType~=2);...
RHS2(ipv);];
if useLU
if IS_OCTAVE
x = MxQ * MxQx* (MxU \ (MxL \ (MxP * RHSx))) ;
else
x = MxQ * (MxU \ (MxL \ (MxP * RHSx))) ;
end
else
x=LHS_matx\RHSx;
end
V(busType~=2,i+1)=x(1:(npq+npv))+1j*x(((npq+npv)+1):(2*(npq+npv)));
Q(ipv,i+1)=x((2*(npq+npv)+1):end);
Cx=real(V(:,i+1));Dx=imag(V(:,i+1));
RHS3xxr=RHS3r-E0.*Cx+F0.*Dx;RHS3xxi=RHS3i-F0.*Cx-E0.*Dx;
Ex=C0i.*RHS3xxr+D0i.*RHS3xxi;Fx=-D0i.*RHS3xxr+C0i.*RHS3xxi;
W(busType~=2,i+1)=Ex(busType~=2)+1j*Fx(busType~=2);
if ~isempty(ind)
for j=1:nInd
% if Rind0(j)~=0
% tempIL=squeeze(LHS_MatInd_Shr(j,:,:))*[real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+rhsBus(:,j);
% tempIRs=-squeeze(LHS_MatInd_Shr2(j,:,:))*[tempIL;real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+...
% squeeze(LHS_MatInd_Shr3(j,:,:))*[rhsI(j);real(rhsIL(j));imag(rhsIL(j))];
% IL(j,i+1)=tempIL(1)+1j*tempIL(2);
% IR(j,i+1)=tempIRs(1)+1j*tempIRs(2);
% s(j,i+1)=tempIRs(3);
% Vm(j,i+1)=V(indIdx(j),i+1)-IL(j,i+1)*Z1(j);
% else % singular case
% if i==1
% % s and IR already generated.
% % s(j,i+1)=R2(j)*Rind1(j)*T0(j)/...
% % ((C0(indIdx(j))-R1(j)*JL0(j)+X1(j)*KL0(j))*(C0(indIdx(j))-R1(j)*JL0(j)+X1(j)*KL0(j))+...
% % (D0(indIdx(j))-X1(j)*JL0(j)-R1(j)*KL0(j))*(D0(indIdx(j))-X1(j)*JL0(j)-R1(j)*KL0(j)));
% % IR(j,i+1)=s(j,i+1)/R2(j)*((C0(indIdx(j))-R1(j)*JL0(j)+X1(j)*KL0(j))+1j*(D0(indIdx(j))-X1(j)*JL0(j)-R1(j)*KL0(j)));
% IL(j,i+1)=(V(indIdx(j),i+1)*Yeind0(j)+IR(j,i+1)+V(indIdx(j),i)*Yeind1(j)-IL(j,i)*Ye1ind1(j))/(1+Ye1ind0(j));
% Vm(j,i+1)=V(indIdx(j),i+1)-IL(j,i+1)*Z1(j);
% else
% tempIL=squeeze(LHS_MatInd_Shr(j,:,:))*[real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+rhsBus(:,j);
% tempIRs=-squeeze(LHS_MatInd_Shr2(j,:,:))*[tempIL;real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+...
% squeeze(LHS_MatInd_Shr3(j,:,:))*[rhsIx(j);real(rhsIL(j));imag(rhsIL(j))];
% IL(j,i+1)=tempIL(1)+1j*tempIL(2);
% IR(j,i+1)=tempIRs(1)+1j*tempIRs(2);
% s(j,i+1)=tempIRs(3);
% Vm(j,i+1)=V(indIdx(j),i+1)-IL(j,i+1)*Z1(j);
% end
% end
tempx=squeeze(LHS_MatInd_Full(j,:,:))*[real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+rhsBus(:,j);
IL(j,i+1)=tempx(3)+1j*tempx(4);
IR(j,i+1)=tempx(1)+1j*tempx(2);
s(j,i+1)=tempx(5);
Vm(j,i+1)=V(indIdx(j),i+1)-IL(j,i+1)*Z1(j);
end
end
if ~isempty(zip)
IiL(:,i+1)=(LHS_MatZip(:,1)+1j*LHS_MatZip(:,3)).*real(V(zipIdx,i+1))+(LHS_MatZip(:,2)+1j*LHS_MatZip(:,4)).*imag(V(zipIdx,i+1))+(RHSILr_full+1j*RHSILi_full);
BiL(:,i+1)=Mat_BZip(:,1).*real(V(zipIdx,i+1))+Mat_BZip(:,2).*imag(V(zipIdx,i+1))+RHS_BZip;
end
if ~isempty(syn)
IGJKd=MatGBiA*[real(V(:,i+1));imag(V(:,i+1))]+RHSIG;
JG(:,i+1)=IGJKd(1:nSyn);
KG(:,i+1)=IGJKd((nSyn+1):(2*nSyn));
d(:,i+1)=IGJKd((2*nSyn+1):(3*nSyn));
Cd(:,i+1)=A1n.*d(:,i+1)+AG0;
Sd(:,i+1)=B1n.*d(:,i+1)+BG0;
end
end
end
Reference in new issue View Git Blame Copy Permalink