ANL-CEEESA
/
powersas.m
mirror of https://github.com/ANL-CEEESA/powersas.m
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

implement DT solvers

Browse Source
pull/7/head
Yang Liu 2 years ago
parent a2c73763da
commit 587bedaee4
19 changed files with 1088 additions and 135 deletions
Show Stats Download Patch File Download Diff File

7
calCPF.m
Unescape View File

@ -64,7 +64,12 @@ else
snapshot=[]; snapshot=[];
hotStart=0; hotStart=0;
end end
simSettings.eventList=[simSettings.eventList;[2 0.0000 options.simLen 50 1 0.0 0.0000]]; if options.method==4
simSettings.eventList(1,6) = 4;
simSettings.eventList=[simSettings.eventList;[2 0.0000 options.simLen 50 1 4.0 0.0000]];
else
simSettings.eventList=[simSettings.eventList;[2 0.0000 options.simLen 50 1 0.0 0.0000]];
end
simSettings.evtDynZipRamp=zipRamp; simSettings.evtDynZipRamp=zipRamp;
simSettings.evtDynPV=pvRamp; simSettings.evtDynPV=pvRamp;
simSettings.evtDyn=zeros(1,23); simSettings.evtDyn(1,1)=1; simSettings.evtDyn=zeros(1,23); simSettings.evtDyn(1,1)=1;

4
calTSA.m
Unescape View File

@ -78,10 +78,14 @@ else
nAggFault=size(aggFaultIdxAux,1)-1; nAggFault=size(aggFaultIdxAux,1)-1;
simSettings.evtFault=[(1:nAggFault)',aggFaultIdxAux(1:nAggFault),aggFaultIdxAux(2:nAggFault+1)-1]; simSettings.evtFault=[(1:nAggFault)',aggFaultIdxAux(1:nAggFault),aggFaultIdxAux(2:nAggFault+1)-1];
simSettings.evtFaultSpec=faultEventsSort(:,1:end-1); simSettings.evtFaultSpec=faultEventsSort(:,1:end-1);
simSettings.eventList(1,6) = options.method;
simSettings.eventList=[simSettings.eventList;[(2:nAggFault+1)',faultEventsSort(aggFaultIdxAux(1:nAggFault),end),... simSettings.eventList=[simSettings.eventList;[(2:nAggFault+1)',faultEventsSort(aggFaultIdxAux(1:nAggFault),end),...
zeros(nAggFault,1),6*ones(nAggFault,1),(1:nAggFault)',repmat(simSettings.eventList(1,[6,7]),nAggFault,1)]]; zeros(nAggFault,1),6*ones(nAggFault,1),(1:nAggFault)',repmat(simSettings.eventList(1,[6,7]),nAggFault,1)]];
simSettings.eventList=[simSettings.eventList;[nAggFault+2,options.simLen,0,99,0,simSettings.eventList(1,[6,7])]]; simSettings.eventList=[simSettings.eventList;[nAggFault+2,options.simLen,0,99,0,simSettings.eventList(1,[6,7])]];
res=runDynamicSimulationExec(caseName,SysData,simSettings,options,snapshot); res=runDynamicSimulationExec(caseName,SysData,simSettings,options,snapshot);
end end

11
internal/calculateInitialState.m
Unescape View File

@ -43,6 +43,7 @@ function [SysDataUpd,x0,finalAlpha,alphaList,diff,t,stateCurve]=calculateInitial
MatGV0,MatGV1,MatGRhs0,MatGRhs1,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara); MatGV0,MatGV1,MatGRhs0,MatGRhs1,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara);
ltc=[];sup=[];dem=[];busName=[]; ltc=[];sup=[];dem=[];busName=[];
% pv=[]; % with synchronous generators, the PV buses are not useful % pv=[]; % with synchronous generators, the PV buses are not useful
[bus,pv,pq,sw,line,shunt,ind,zip,syn,exc,tg,agc,cac,cluster,pm,oldToNew,newToOld]=regulateSystemData(bus,pv,pq,sw,line,shunt,ind,zip,syn,exc,tg,agc,cac,cluster);
[nIslands,islands]=searchIslands(bus(:,1),line(:,[1,2])); [nIslands,islands]=searchIslands(bus(:,1),line(:,[1,2]));
if nIslands>1 if nIslands>1
addLog('There are more than 1 islands in the system. Only the largest island will be analyzed.','INFO'); addLog('There are more than 1 islands in the system. Only the largest island will be analyzed.','INFO');
@ -50,7 +51,9 @@ if nIslands>1
[bus,~,pq,sw,line,ltc,sup,dem,shunt,ind,zip,syn,busName,newToOld,oldToNew,pvInd,pqInd,swInd,lineInd,ltcInd,supInd,demInd,shuntInd,indInd,zipInd,synInd] ... [bus,~,pq,sw,line,ltc,sup,dem,shunt,ind,zip,syn,busName,newToOld,oldToNew,pvInd,pqInd,swInd,lineInd,ltcInd,supInd,demInd,shuntInd,indInd,zipInd,synInd] ...
=forgeLargestIsland(maxIsland,nIslands,islands,bus,[],pq,sw,line,ltc,sup,dem,shunt,ind,zip,syn,busName); =forgeLargestIsland(maxIsland,nIslands,islands,bus,[],pq,sw,line,ltc,sup,dem,shunt,ind,zip,syn,busName);
end end
if ~isempty(Pm1)
pm=Pm1;
end
nbus=size(bus,1); nbus=size(bus,1);
nline=size(line,1); nline=size(line,1);
nIslands=1; nIslands=1;
@ -79,11 +82,7 @@ if isempty(syn)&&isempty(sw)
addLog('Steady-state analysis needs at Swing bus or synchronous machine(s) in the system.','ERROR'); addLog('Steady-state analysis needs at Swing bus or synchronous machine(s) in the system.','ERROR');
return; return;
end end
% [bus,pv,pq,sw,line,shunt,ind,zip,syn,exc,tg,agc,cac,cluster,pm,oldToNew,newToOld]=regulateSystemData(bus,pv,pq,sw,line,shunt,ind,zip,syn,exc,tg,agc,cac,cluster);
[bus,pv,pq,sw,line,shunt,ind,zip,syn,exc,tg,agc,cac,cluster,pm,oldToNew,newToOld]=regulateSystemData(bus,pv,pq,sw,line,shunt,ind,zip,syn,exc,tg,agc,cac,cluster);
if ~isempty(Pm1)
pm=Pm1;
end
nbus=size(bus,1); nbus=size(bus,1);
nline=size(line,1); nline=size(line,1);
% nlvl=10; % nlvl=10;

11
internal/formatEventList.m
Unescape View File

@ -91,8 +91,15 @@ if ~isempty(sstTime)
end end
end end
end end
eventAdd=[zeros(size(stTime,1),1),stTime,endTime,EvtType.DYN_SIM*ones(size(stTime,1),1),zeros(size(stTime,1),1),MethodType.FULL_HE*ones(size(stTime,1),1),SimSet.DEFAULT_DT*ones(size(stTime,1),1)]; %khuang 5 Jul
if event(1,6) ==0
eventAdd=[zeros(size(stTime,1),1),stTime,endTime,EvtType.DYN_SIM*ones(size(stTime,1),1),zeros(size(stTime,1),1),MethodType.FULL_HE*ones(size(stTime,1),1),SimSet.DEFAULT_DT*ones(size(stTime,1),1)];
elseif event(1,7) >0
eventAdd=[zeros(size(stTime,1),1),stTime,endTime,EvtType.DYN_SIM*ones(size(stTime,1),1),zeros(size(stTime,1),1),event(1,6)*ones(size(stTime,1),1),event(1,7)*ones(size(stTime,1),1)];
else
eventAdd=[zeros(size(stTime,1),1),stTime,endTime,EvtType.DYN_SIM*ones(size(stTime,1),1),zeros(size(stTime,1),1),event(1,6)*ones(size(stTime,1),1),SimSet.DEFAULT_DT*ones(size(stTime,1),1)];
end
%%%%%%%%%%%%%%%%
eventExt=[eventSort;eventAdd]; eventExt=[eventSort;eventAdd];
[~,iSort]=sort(eventExt(:,2)+(eventExt(:,4)==EvtType.DYN_SIM)*minDiffTime/2+(eventExt(:,4)==EvtType.DYN_SIM).*eventExt(:,3)/eSortTemp(end)*minDiffTime/4); [~,iSort]=sort(eventExt(:,2)+(eventExt(:,4)==EvtType.DYN_SIM)*minDiffTime/2+(eventExt(:,4)==EvtType.DYN_SIM).*eventExt(:,3)/eSortTemp(end)*minDiffTime/4);
eventExt=eventExt(iSort,:); eventExt=eventExt(iSort,:);

4
internal/generalDynSimulation.m
Unescape View File

@ -118,6 +118,10 @@ if DEmethod==0
[stateUpd,t,finalAlpha,alphaList,diffList,exitFlag]=simulationTimeDomainHem(SimData,SysDataTmp,SysPara,x0); [stateUpd,t,finalAlpha,alphaList,diffList,exitFlag]=simulationTimeDomainHem(SimData,SysDataTmp,SysPara,x0);
alphaList=cumsum(alphaList); alphaList=cumsum(alphaList);
SimDataUpd=SimData; SimDataUpd=SimData;
elseif DEmethod==4
[stateUpd,t,finalAlpha,alphaList,diffList,exitFlag]=simulationTimeDomainDT(SimData,SysDataTmp,SysPara,x0);
alphaList=cumsum(alphaList);
SimDataUpd=SimData;
else else
[stateUpd,t,diffList,nxtDt,exitFlag]=simulationTimeDomainNI(SimData,SysDataTmp,SysPara,x0); [stateUpd,t,diffList,nxtDt,exitFlag]=simulationTimeDomainNI(SimData,SysDataTmp,SysPara,x0);
finalAlpha=t(end); finalAlpha=t(end);

1
internal/getSimMethodType.m
Unescape View File

@ -30,6 +30,7 @@ function MethodType=getSimMethodType()
% %
MethodType.FULL_HE=0.0; MethodType.FULL_HE=0.0;
MethodType.FULL_DT=4.0;
MethodType.ME_HE=1.0; MethodType.ME_HE=1.0;
MethodType.ME_NR=1.1; MethodType.ME_NR=1.1;
MethodType.RK4_HE=2.0; MethodType.RK4_HE=2.0;

7
internal/getseq.m
Unescape View File

@ -28,6 +28,13 @@ function seq=getseq(n,d)
% *************************************************************************************************** % ***************************************************************************************************
% %
%NCHOOSEK(N,K) where N and K are non-negative integers returns N!/K!(N-K)!.
%NCHOOSEK(V,K) where V is a vector of length N, produces a matrix
%with N!/K!(N-K)! rows and K columns. Each row of the result has K of
%the elements in the vector V. This syntax is only practical for
%situations where N is less than about 15.
aux=nchoosek(1:(n+d-1),d-1); aux=nchoosek(1:(n+d-1),d-1);
aaux=[zeros(size(aux,1),1),aux,(n+d)*ones(size(aux,1),1)]; aaux=[zeros(size(aux,1),1),aux,(n+d)*ones(size(aux,1),1)];
seq=aaux(:,2:end)-aaux(:,1:end-1)-1; seq=aaux(:,2:end)-aaux(:,1:end-1)-1;

398
internal/hemMachinePFSalientcontinueDyn.m
Unescape View File

@ -42,20 +42,32 @@ function [V,Q,s,d,w,eq1,eq2,ed1,ed2,psid,psiq,Pm,Ef,Vavrm,Vavrr,Vavrf,Vavrref,tg
global IS_OCTAVE; global IS_OCTAVE;
% import system data khuang, 8 Jul
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData); [bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
% nbus:Total number of buses khuang, 8 Jul
nbus=size(bus,1); nbus=size(bus,1);
nline=size(line,1); nline=size(line,1);
%determine islanding, this is for AGC part khuang 8 JUl
if isfield(SysPara,'nIslands')&&isfield(SysPara,'islands')&&isfield(SysPara,'refs') if isfield(SysPara,'nIslands')&&isfield(SysPara,'islands')&&isfield(SysPara,'refs')
nIslands=SysPara.nIslands;islands=SysPara.islands;refs=SysPara.refs; nIslands=SysPara.nIslands;islands=SysPara.islands;refs=SysPara.refs;
else else
[nIslands,islands,refs]=searchIslands(bus(:,1),line(:,1:2)); [nIslands,islands,refs]=searchIslands(bus(:,1),line(:,1:2));
end end
% improt initial condition of state variables khuang 8 JUl
[V0,Q0,s0,d0,w0,eq10,eq20,ed10,ed20,psid0,psiq0,Pm0,Ef0,Vavrm0,Vavrr0,Vavrf0,Vavrref0,tgovg0,tgovm0,tgovmech0,f0,dpg0,qplt0,vg0]=unfoldX(x0,SysData); [V0,Q0,s0,d0,w0,eq10,eq20,ed10,ed20,psid0,psiq0,Pm0,Ef0,Vavrm0,Vavrr0,Vavrf0,Vavrref0,tgovg0,tgovm0,tgovmech0,f0,dpg0,qplt0,vg0]=unfoldX(x0,SysData);
% import simualtion data khuang 8 JUl
[~,~,~,nlvl,taylorN,~,~,~,~]=unfoldSimData(SimData); [~,~,~,nlvl,taylorN,~,~,~,~]=unfoldSimData(SimData);
% import system parameters khuang 8 JUl
[pqIncr,pvIncr,Rind0,Rind1,Reind0,Reind1,Rzip0,Rzip1,Ytr0,Ytr1,Ysh0,Ysh1,VspSq2,~,~,~,~,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara); [pqIncr,pvIncr,Rind0,Rind1,Reind0,Reind1,Rzip0,Rzip1,Ytr0,Ytr1,Ysh0,Ysh1,VspSq2,~,~,~,~,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara);
% %
% PQ incremental is given for every pq bus.
% need to figure how to utilize them, now i suppose this part is for CPF
% problem khuang 8 JUL
Pls=zeros(nbus,2);Pls(pq(:,1),1)=pqIncr(:,1);if ~isempty(pv);Pls(pv(:,1),1)=Pls(pv(:,1),1)-pvIncr;end Pls=zeros(nbus,2);Pls(pq(:,1),1)=pqIncr(:,1);if ~isempty(pv);Pls(pv(:,1),1)=Pls(pv(:,1),1)-pvIncr;end
Qls=zeros(nbus,2);Qls(pq(:,1),1)=pqIncr(:,2); Qls=zeros(nbus,2);Qls(pq(:,1),1)=pqIncr(:,2);
if size(pqIncr,2)>=4 if size(pqIncr,2)>=4
@ -63,10 +75,14 @@ if size(pqIncr,2)>=4
Qls(pq(:,1),2)=pqIncr(:,4); Qls(pq(:,1),2)=pqIncr(:,4);
end end
% formatting Ymatrix, here the default value of fault is empty, khuang 8 JUL
if isempty(Ytr0) if isempty(Ytr0)
[Y,Ytr0,Ysh,ytrfr,ytrto,yshfr,yshto]=getYMatrix(nbus,line); [Y,Ytr0,Ysh,ytrfr,ytrto,yshfr,yshto]=getYMatrix(nbus,line);
end end
% reshape the pv, pq shunt and swing buses as zeros if they are empty
% khuang 8JUL
busType=zeros(nbus,1); busType=zeros(nbus,1);
if isempty(pv) if isempty(pv)
pv=zeros(0,6); pv=zeros(0,6);
@ -80,33 +96,54 @@ end
if isempty(sw) if isempty(sw)
sw=zeros(0,13); sw=zeros(0,13);
end end
% label pv and swing buses khuang 8 JUL
% 1: PV bus, 0: PQ bus
busType(pv(:,1))=1; busType(pv(:,1))=1;
busType(sw(:,1))=2; busType(sw(:,1))=2;
% zip(busType(zip(:,1))~=0,10)=0; % zip(busType(zip(:,1))~=0,10)=0;
% index of slack bus (isw), pv bus(ipv), and pq bus(ipq)
% is given: isw, ipv, and ipq khuang 8 JUL
% Additionally, number of pv and pq buses are given:
%npv, npq respectively khuang 8 JUL
isw=find(busType==2); isw=find(busType==2);
ipv=find(busType==1); ipv=find(busType==1);
ipq=find(busType==0); ipq=find(busType==0);
npq=size(ipq,1); npq=size(ipq,1);
npv=size(ipv,1); npv=size(ipv,1);
% shunt capacitator is initialized as yShunt which is a complex number.
% for every bus khuang 8 JUL
yShunt=zeros(nbus,1); yShunt=zeros(nbus,1);
yShunt(shunt(:,1))=shunt(:,5)+1j*shunt(:,6); yShunt(shunt(:,1))=shunt(:,5)+1j*shunt(:,6);
% check if zip load exists in the system
% and initialized zip load khuang 8 JUL
if ~isempty(zip)%zipMode=0 if ~isempty(zip)%zipMode=0
Ysh0=Ysh0+accumarray(zip(:,1),Rzip0.*(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]); 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]); Ysh1=Ysh1+accumarray(zip(:,1),Rzip1.*(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]);
end end
% now zip load + shunt khuang 8 JUL
Ysh0=Ysh0+yShunt; Ysh0=Ysh0+yShunt;
% Y=Y+sparse(1:nbus,1:nbus,yShunt,nbus,nbus); % Y=Y+sparse(1:nbus,1:nbus,yShunt,nbus,nbus);
% now zip load + shunt+ network Y matrix khuang 8 JUL
Y=Ytr0+sparse(1:nbus,1:nbus,Ysh0,nbus,nbus); Y=Ytr0+sparse(1:nbus,1:nbus,Ysh0,nbus,nbus);
%initialize P and Q for every bus khuang 8 JUL
pVec=zeros(nbus,1); pVec=zeros(nbus,1);
qVec=zeros(nbus,1); qVec=zeros(nbus,1);
% vSp=zeros(nbus,1); % vSp=zeros(nbus,1);
% need to figure out the meaning of index 1, 4,5
% based on Kaiyang's understanding, 1 is load side and 4&5 are generators'
% output khuang 8 JUL
pVec(pv(:,1))=pVec(pv(:,1))+pv(:,4); pVec(pv(:,1))=pVec(pv(:,1))+pv(:,4);
pVec(pq(:,1))=pVec(pq(:,1))-pq(:,4); pVec(pq(:,1))=pVec(pq(:,1))-pq(:,4);
qVec(pq(:,1))=qVec(pq(:,1))-pq(:,5); qVec(pq(:,1))=qVec(pq(:,1))-pq(:,5);
% account the zip load, i.e dynamic load khuang 8 JUL
if ~isempty(zip)%zipMode=0, account for the PQ components in ZIP loads 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]); pVec=pVec-accumarray(zip(:,1),Rzip0.*zip(:,7).*zip(:,12),[nbus,1]);
qVec=qVec-accumarray(zip(:,1),Rzip0.*zip(:,10).*zip(:,12),[nbus,1]); qVec=qVec-accumarray(zip(:,1),Rzip0.*zip(:,10).*zip(:,12),[nbus,1]);
@ -114,21 +151,39 @@ end
% qVec(ipv)=qVec(ipv)+Q0(ipv); % qVec(ipv)=qVec(ipv)+Q0(ipv);
% vSp(ipv)=pv(:,5); % vSp(ipv)=pv(:,5);
% so far, initialization of PQ for every bus and Y matrix is ready khuang 8 JUL
% initiliza voltage V and W = 1/V for every bus khuang 8 JUL
V=zeros(nbus,nlvl+1); V=zeros(nbus,nlvl+1);
V(:,1)=V0; V(:,1)=V0;
W=zeros(nbus,nlvl+1); W=zeros(nbus,nlvl+1);
W(:,1)=1./V0; W(:,1)=1./V0;
% initiliza magnitude of voltage V for every bus khuang 8 JUL
Vmag=zeros(nbus,nlvl+1); Vmag=zeros(nbus,nlvl+1);
Vmag(:,1)=abs(V0); Vmag(:,1)=abs(V0);
% Power is initilized as we already cooked pVec and qVec khuang 8 JUL
P=zeros(nbus,nlvl+1); P=zeros(nbus,nlvl+1);
P(:,1)=pVec; P(:,1)=pVec;
% P(isw,2:end)=0; % P(isw,2:end)=0;
% here we need to figure out what Q extra mean, and difference from Q
% notice that Q0 is initialized with sysmdata but not P0 khuang 8 JUL
Q=zeros(nbus,nlvl+1); Q=zeros(nbus,nlvl+1);
Qxtra=zeros(size(Q)); Qxtra=zeros(size(Q));
Q(:,1)=Q0; Q(:,1)=Q0;
Qxtra(:,1)=qVec; Qxtra(:,1)=qVec;
% Also, the meaning of Pls and Qls need to be verified, i assume they are
% for CPF khuang 8 JUL
P(:,2:(size(Pls,2)+1))=-Pls; P(:,2:(size(Pls,2)+1))=-Pls;
Qxtra(:,2:(size(Qls,2)+1))=-Qls; Qxtra(:,2:(size(Qls,2)+1))=-Qls;
% In the previous, pVec and qvec are considered zip load, here Pls and Qls
% are not, so we need to do it.khuang 8 JUL
if ~isempty(zip) if ~isempty(zip)
P(:,2)=P(:,2)-accumarray(zip(:,1),Rzip1.*zip(:,7).*zip(:,12),[nbus,1]); 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]); Qxtra(:,2)=Qxtra(:,2)-accumarray(zip(:,1),Rzip1.*zip(:,10).*zip(:,12),[nbus,1]);
@ -136,11 +191,17 @@ end
% Qxtra(busType~=0,2:end)=Q(busType~=0,2:end); % Qxtra(busType~=0,2:end)=Q(busType~=0,2:end);
% Q(busType~=0,2:end)=0; % Q(busType~=0,2:end)=0;
% seperate real and image part of voltages and their inverse khuang 8 JUL
% here V = C+1i*D khuang 8 JUL
% and W = 1./V = E + 1i*F khuang 8 JUL
C0=real(V(:,1)); C0=real(V(:,1));
D0=imag(V(:,1)); D0=imag(V(:,1));
E0=real(W(:,1)); E0=real(W(:,1));
F0=imag(W(:,1)); F0=imag(W(:,1));
% Construct sparse matrix individually for C,D,E,F,P,Q. khuang 8 JUL
% Notice that Q = Q(:,1)+Qxtra(:,1) which is different from P. khuang 8 JUL
C0M=sparse(1:nbus,1:nbus,C0,nbus,nbus); C0M=sparse(1:nbus,1:nbus,C0,nbus,nbus);
D0M=sparse(1:nbus,1:nbus,D0,nbus,nbus); D0M=sparse(1:nbus,1:nbus,D0,nbus,nbus);
E0M=sparse(1:nbus,1:nbus,E0,nbus,nbus); E0M=sparse(1:nbus,1:nbus,E0,nbus,nbus);
@ -148,10 +209,15 @@ F0M=sparse(1:nbus,1:nbus,F0,nbus,nbus);
P0M=sparse(1:nbus,1:nbus,P(:,1),nbus,nbus); P0M=sparse(1:nbus,1:nbus,P(:,1),nbus,nbus);
Q0M=sparse(1:nbus,1:nbus,Q(:,1)+Qxtra(:,1),nbus,nbus); Q0M=sparse(1:nbus,1:nbus,Q(:,1)+Qxtra(:,1),nbus,nbus);
% get real part and image part of Y matrix, not sure why do this. khuang 8 JUL
G=real(Y); G=real(Y);
B=imag(Y); B=imag(Y);
% so Y = G + 1i*B. khuang 8 JUL
% Determine the frequency model of each island % this part is for AGC khuang 8 JUL
%--------------------------------------------------------------
% Determine the frequency model of each island
% 0:sw,1:syn,2:steady-state f
freqTypeTag=zeros(nIslands,1);%0:sw,1:syn,2:steady-state f freqTypeTag=zeros(nIslands,1);%0:sw,1:syn,2:steady-state f
freqKeptTag=zeros(nbus,1); freqKeptTag=zeros(nbus,1);
frefs=refs; frefs=refs;
@ -177,16 +243,26 @@ end
freqKeptTagxRef=freqKeptTag; freqKeptTagxRef=freqKeptTag;
freqKeptTagxRef(frefs)=0; freqKeptTagxRef(frefs)=0;
nFreqKept=sum(freqKeptTag); nFreqKept=sum(freqKeptTag);
%-----------------------------------------------------------
if ~isempty(ind)
nInd=size(ind,1);
indIdx=ind(:,1);
s=zeros(nInd,nlvl+1); %------------------------------------------------------------------
s(:,1)=s0; % this part is for initialling inductor. khuang 8 JUL
IL=zeros(nInd,nlvl+1); if ~isempty(ind) % check if there is any inductor% khuang 8 JUL
IR=zeros(nInd,nlvl+1); nInd=size(ind,1); % determine the number of inductors % khuang 8 JUL
Vm=zeros(nInd,nlvl+1); indIdx=ind(:,1); % store the index of inductors among all buses% khuang 8 JUL
s=zeros(nInd,nlvl+1); % slip% khuang 8 JUL
s(:,1)=s0; % initialize slip% khuang 8 JUL
IL=zeros(nInd,nlvl+1); % |
IR=zeros(nInd,nlvl+1); % |
Vm=zeros(nInd,nlvl+1); % initialization finished 0 value% khuang 8 JUL
%-----------------parameters of inductors---------% khuang 8 JUL
%-----------------START----------------% khuang 8 JUL
R1=ind(:,7); R1=ind(:,7);
X1=ind(:,8); X1=ind(:,8);
Z1=ind(:,7)+1j*ind(:,8); Z1=ind(:,7)+1j*ind(:,8);
@ -197,11 +273,15 @@ if ~isempty(ind)
T1=-ind(:,16)-2*ind(:,17); T1=-ind(:,16)-2*ind(:,17);
T2=ind(:,17); T2=ind(:,17);
Hm=ind(:,14); Hm=ind(:,14);
%-----------------parameters of inductors---------% khuang 8 JUL
%-----------------END----------------% khuang 8 JUL
Rm=zeros(nInd,1); Rm=zeros(nInd,1);
Am=sparse(indIdx,(1:nInd)',ones(1,nInd),nbus,nInd); Am=sparse(indIdx,(1:nInd)',ones(1,nInd),nbus,nInd);
% first order value of induction motor IL,VM,IR % khuang 8 JUL
IL(:,1)=V0(indIdx)./(Z1+Ze.*(R2+1j*X2.*s0)./(R2.*Reind0+(1j*X2.*Reind0+Ze).*s0)); IL(:,1)=V0(indIdx)./(Z1+Ze.*(R2+1j*X2.*s0)./(R2.*Reind0+(1j*X2.*Reind0+Ze).*s0));
Vm(:,1)=V0(indIdx)-IL(:,1).*Z1; Vm(:,1)=V0(indIdx)-IL(:,1).*Z1;
IR(:,1)=Vm(:,1).*s0./(R2+1j*X2.*s0); IR(:,1)=Vm(:,1).*s0./(R2+1j*X2.*s0);
@ -211,6 +291,7 @@ if ~isempty(ind)
JL0=real(IL(:,1)); JL0=real(IL(:,1));
KL0=imag(IL(:,1)); KL0=imag(IL(:,1));
% prepare the algebric matrix % khuang 8 JUL
Yeind0=Reind0./Ze; Yeind0=Reind0./Ze;
Yeind1=Reind1./Ze; Yeind1=Reind1./Ze;
Ye1ind0=Reind0.*Z1./Ze; Ye1ind0=Reind0.*Z1./Ze;
@ -290,7 +371,7 @@ if ~isempty(ind)
RHS_C_Shr_sqz=[temp2_c12_inv_sqz,... RHS_C_Shr_sqz=[temp2_c12_inv_sqz,...
-[temp2_c12_inv_sqz(:,1).*temp1_sqz(:,1)+temp2_c12_inv_sqz(:,3).*temp1_sqz(:,2),temp2_c12_inv_sqz(:,2).*temp1_sqz(:,1)+temp2_c12_inv_sqz(:,4).*temp1_sqz(:,2),... -[temp2_c12_inv_sqz(:,1).*temp1_sqz(:,1)+temp2_c12_inv_sqz(:,3).*temp1_sqz(:,2),temp2_c12_inv_sqz(:,2).*temp1_sqz(:,1)+temp2_c12_inv_sqz(:,4).*temp1_sqz(:,2),...
temp2_c12_inv_sqz(:,1).*temp1_sqz(:,3)+temp2_c12_inv_sqz(:,3).*temp1_sqz(:,4),temp2_c12_inv_sqz(:,2).*temp1_sqz(:,3)+temp2_c12_inv_sqz(:,4).*temp1_sqz(:,4)]];% R\[I,-CA^-1] temp2_c12_inv_sqz(:,1).*temp1_sqz(:,3)+temp2_c12_inv_sqz(:,3).*temp1_sqz(:,4),temp2_c12_inv_sqz(:,2).*temp1_sqz(:,3)+temp2_c12_inv_sqz(:,4).*temp1_sqz(:,4)]];% R\[I,-CA^-1]
% will be used to calculate algebric variabls for motors% khuang 8 JUL
LHS_MatInd_Bus_sqz=zeros(nbus,4); % \sum{-R\S} by buses LHS_MatInd_Bus_sqz=zeros(nbus,4); % \sum{-R\S} by buses
LHS_MatInd_Bus_sqz(:,1)=accumarray(indIdx,LHS_MatInd_Shr_sqz(:,1),[nbus,1]); LHS_MatInd_Bus_sqz(:,1)=accumarray(indIdx,LHS_MatInd_Shr_sqz(:,1),[nbus,1]);
LHS_MatInd_Bus_sqz(:,2)=accumarray(indIdx,LHS_MatInd_Shr_sqz(:,2),[nbus,1]); LHS_MatInd_Bus_sqz(:,2)=accumarray(indIdx,LHS_MatInd_Shr_sqz(:,2),[nbus,1]);
@ -299,24 +380,29 @@ if ~isempty(ind)
else else
s=zeros(0,nlvl+1); s=zeros(0,nlvl+1);
end end
% Initialization of inductors is finished % khuang 8 JUL
%------------------------------------------------------------------
%------------------------------Initialization of ZIP load---------% khuang 8 JUL
if ~isempty(zip) if ~isempty(zip)
nZip=size(zip,1); nZip=size(zip,1);
zipIdx=zip(:,1); zipIdx=zip(:,1);
IiL=zeros(nZip,nlvl+1); IiL=zeros(nZip,nlvl+1);
BiL=zeros(nZip,nlvl+1); BiL=zeros(nZip,nlvl+1);
% prepare the necessary matrix by blocks% khuang 8 JUL
Bi0=abs(V0(zipIdx)); Bi0=abs(V0(zipIdx));
JI=zip(:,6); JI=zip(:,6);
KI=-zip(:,9); KI=-zip(:,9);
% current % khuang 8 JUL
Ii0L=Rzip0.*(JI+1j*KI).*V0(zipIdx)./Bi0; Ii0L=Rzip0.*(JI+1j*KI).*V0(zipIdx)./Bi0;
Ji0L=real(Ii0L); Ji0L=real(Ii0L);
Ki0L=imag(Ii0L); Ki0L=imag(Ii0L);
IiL(:,1)=Ii0L; IiL(:,1)=Ii0L;
BiL(:,1)=Bi0; BiL(:,1)=Bi0;
% voltage% khuang 8 JUL
Ci0=real(V0(zipIdx)); Ci0=real(V0(zipIdx));
Di0=imag(V0(zipIdx)); Di0=imag(V0(zipIdx));
@ -326,38 +412,43 @@ if ~isempty(zip)
else else
IiL=zeros(0,nlvl+1); IiL=zeros(0,nlvl+1);
end end
%------------------------------Initialization of ZIP load------------------% khuang 8 JUL
%------------------------------Initialization of ZIP load is finished---------------- % khuang 8 JUL
%------------------------------Initialization of GEN------------------% khuang 8 JUL
%------------------------------Start------------------------% khuang 8 JUL
nSyn=size(syn,1); nSyn=size(syn,1);
if ~isempty(syn) if ~isempty(syn)
synIdx=syn(:,1); synIdx =syn(:,1);% index number of Generators% khuang 8 JUL
wgb=syn(:,4); wgb =syn(:,4);% maybe the base value% khuang 8 JUL
modSyn=syn(:,5); modSyn =syn(:,5);% the order of generator models% khuang 8 JUL
Xgl=syn(:,6); Xgl =syn(:,6);
Rga=syn(:,7); Rga =syn(:,7);
Xgd=syn(:,8); Xgd =syn(:,8);
Xgd1=syn(:,9); Xgd1 =syn(:,9);
Xgd2=syn(:,10); Xgd2 =syn(:,10);
Tgd1=syn(:,11); Tgd1 =syn(:,11);
Tgd2=syn(:,12); Tgd2 =syn(:,12);
Xgq=syn(:,13); Xgq =syn(:,13);
Xgq1=syn(:,14); Xgq1 =syn(:,14);
Xgq2=syn(:,15); Xgq2 =syn(:,15);
Tgq1=syn(:,16); Tgq1 =syn(:,16);
Tgq2=syn(:,17); Tgq2 =syn(:,17);
Mg=syn(:,18); Mg =syn(:,18);
Dg=syn(:,19); Dg =syn(:,19);
TgAA=syn(:,24); TgAA =syn(:,24);
gammad=Tgd2./Tgd1.*Xgd2./Xgd1.*(Xgd-Xgd1); gammad =Tgd2./Tgd1.*Xgd2./Xgd1.*(Xgd-Xgd1);
gammaq=Tgq2./Tgq1.*Xgq2./Xgq1.*(Xgq-Xgq1); gammaq =Tgq2./Tgq1.*Xgq2./Xgq1.*(Xgq-Xgq1);
d=zeros(nSyn,nlvl+1); d=zeros(nSyn,nlvl+1); % delta% khuang 8 JUL
w=zeros(nSyn,nlvl+1); w=zeros(nSyn,nlvl+1); % omega% khuang 8 JUL
eq1=zeros(nSyn,nlvl+1); eq1=zeros(nSyn,nlvl+1); %eq'% khuang 8 JUL
eq2=zeros(nSyn,nlvl+1); eq2=zeros(nSyn,nlvl+1); %eq''% khuang 8 JUL
ed1=zeros(nSyn,nlvl+1); ed1=zeros(nSyn,nlvl+1); %ed'% khuang 8 JUL
ed2=zeros(nSyn,nlvl+1); ed2=zeros(nSyn,nlvl+1); %ed''% khuang 8 JUL
psiq=zeros(nSyn,nlvl+1); psiq=zeros(nSyn,nlvl+1); % not sure, only in 8th order model% khuang 8 JUL
psid=zeros(nSyn,nlvl+1); psid=zeros(nSyn,nlvl+1); % not sure, only in 8th order model% khuang 8 JUL
JG=zeros(nSyn,nlvl+1); JG=zeros(nSyn,nlvl+1);
KG=zeros(nSyn,nlvl+1); KG=zeros(nSyn,nlvl+1);
IGq=zeros(nSyn,nlvl+1); IGq=zeros(nSyn,nlvl+1);
@ -373,7 +464,7 @@ if ~isempty(syn)
sind=sin(d0); sind=sin(d0);
CG0=C0(synIdx); CG0=C0(synIdx);
DG0=D0(synIdx); DG0=D0(synIdx);
% the first value is given here, notice all are 8th order model% khuang 8 JUL
d(:,1)=d0; d(:,1)=d0;
w(:,1)=w0; w(:,1)=w0;
eq1(:,1)=eq10; eq1(:,1)=eq10;
@ -382,13 +473,19 @@ if ~isempty(syn)
ed2(:,1)=ed20; ed2(:,1)=ed20;
psiq(:,1)=psiq0; psiq(:,1)=psiq0;
psid(:,1)=psid0; psid(:,1)=psid0;
% transform between grid side and dq side% khuang 8 JUL
VGd(:,1)=sind.*CG0-cosd.*DG0; VGd(:,1)=sind.*CG0-cosd.*DG0;
VGq(:,1)=cosd.*CG0+sind.*DG0; VGq(:,1)=cosd.*CG0+sind.*DG0;
Cd(:,1)=cosd; % now they are under dq side% khuang 8 JUL
Sd(:,1)=sind;
Cd(:,1)=cosd; % first order of cos(delta)% khuang 8 JUL
Sd(:,1)=sind; % first order of sin(delta)% khuang 8 JUL
Ef(:,1)=Ef0; Ef(:,1)=Ef0;
Pm(:,1)=Pm0; Pm(:,1)=Pm0;
%check if controller exists% khuang 8 JUL
if ~isempty(Ef1) if ~isempty(Ef1)
Ef(:,2)=Ef1; Ef(:,2)=Ef1;
end end
@ -399,26 +496,47 @@ if ~isempty(syn)
Pm(:,2)=Pm1; Pm(:,2)=Pm1;
end end
% notice that here truncated taylor is applied % khuang 8 JUL
% and this is the key differnet from Dt rule% khuang 8 JUL
% Here only at most 5 th order taylor series are considered for sin% khuang 8 JUL
% and cos function % khuang 8 JUL
[cosp,sinp,taylorN]=getTaylorPolynomials(d0,taylorN); % taylorN may be truncated [cosp,sinp,taylorN]=getTaylorPolynomials(d0,taylorN); % taylorN may be truncated
Mats=zeros(nSyn,4); Mats=zeros(nSyn,4);
MatsR=zeros(nSyn,4); MatsR=zeros(nSyn,4);
MatsRs=zeros(nSyn,4); MatsRs=zeros(nSyn,4);
% count the number for different kinds models % khuang 8 JUL
% ex: modelTag = [ 0 0 0 0 0 10 0 0].' % khuang 8 JUL
% ex: there are 10 gens using 6th order model % khuang 8 JUL
modelTag=accumarray(modSyn,ones(nSyn,1),[8,1]); modelTag=accumarray(modSyn,ones(nSyn,1),[8,1]);
% determine the order of the model % khuang 8 JUL
% Do we really need for loop? % khuang 8 JUL
% the answer is yes since different gen may use different% khuang 8 JUL
% order model% khuang 8 JUL
for i=1:nSyn for i=1:nSyn
% 8th order, no need to change% khuang 8 JUL
if modSyn(i)==8 if modSyn(i)==8
IGd(i,1)=(eq20(i)-psid0(i))/Xgd2(i); IGd(i,1)=(eq20(i)-psid0(i))/Xgd2(i);
IGq(i,1)=(-ed20(i)-psiq0(i))/Xgq2(i); IGq(i,1)=(-ed20(i)-psiq0(i))/Xgq2(i);
Mats(i,:)=[sind(i),cosd(i),-cosd(i),sind(i)]; Mats(i,:)=[sind(i),cosd(i),-cosd(i),sind(i)];
% 6th order% khuang 8 JUL
elseif modSyn(i)==6 elseif modSyn(i)==6
% algebric equation to solve Id, Iq% khuang 8 JUL
IGd(i,1)=((ed20(i)-VGd(i,1))*Rga(i)+(eq20(i)-VGq(i,1))*Xgq2(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i)); IGd(i,1)=((ed20(i)-VGd(i,1))*Rga(i)+(eq20(i)-VGq(i,1))*Xgq2(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
IGq(i,1)=(-(ed20(i)-VGd(i,1))*Xgd2(i)+(eq20(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i)); IGq(i,1)=(-(ed20(i)-VGd(i,1))*Xgd2(i)+(eq20(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
% transform matrix (inverse version)% khuang 8 JUL
Mats(i,:)=[sind(i),cosd(i),-cosd(i),sind(i)]; Mats(i,:)=[sind(i),cosd(i),-cosd(i),sind(i)];
% Here matrix is the inverse matrix, to understand this% khuang 8 JUL
% We have A*Ixy+B*Vxy = f => MatsR = A^-1, MatsRs = A^-1*B = MatsRs*B% khuang 8 JUL
% so Ixy = MatsR*f-MatsRs*Vxy, which is used later to% khuang 8 JUL
% eliminate Ixy when disturbance happens% khuang 8 JUL
MatsR(i,:)=[sind(i)*Rga(i)-cosd(i)*Xgd2(i),sind(i)*Xgq2(i)+cosd(i)*Rga(i),-cosd(i)*Rga(i)-sind(i)*Xgd2(i),-cosd(i)*Xgq2(i)+sind(i)*Rga(i)]/... MatsR(i,:)=[sind(i)*Rga(i)-cosd(i)*Xgd2(i),sind(i)*Xgq2(i)+cosd(i)*Rga(i),-cosd(i)*Rga(i)-sind(i)*Xgd2(i),-cosd(i)*Xgq2(i)+sind(i)*Rga(i)]/...
(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i)); (Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),... MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),...
MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)]; MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)];
% 5th order% khuang 8 JUL
elseif modSyn(i)==5 elseif modSyn(i)==5
IGd(i,1)=((ed20(i)-VGd(i,1))*Rga(i)+(eq20(i)-VGq(i,1))*Xgq2(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i)); IGd(i,1)=((ed20(i)-VGd(i,1))*Rga(i)+(eq20(i)-VGq(i,1))*Xgq2(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
IGq(i,1)=(-(ed20(i)-VGd(i,1))*Xgd2(i)+(eq20(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i)); IGq(i,1)=(-(ed20(i)-VGd(i,1))*Xgd2(i)+(eq20(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
@ -427,6 +545,7 @@ if ~isempty(syn)
(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i)); (Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),... MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),...
MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)]; MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)];
% 4th order% khuang 8 JUL
elseif modSyn(i)==4 elseif modSyn(i)==4
IGd(i,1)=((ed10(i)-VGd(i,1))*Rga(i)+(eq10(i)-VGq(i,1))*Xgq1(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i)); IGd(i,1)=((ed10(i)-VGd(i,1))*Rga(i)+(eq10(i)-VGq(i,1))*Xgq1(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i));
IGq(i,1)=(-(ed10(i)-VGd(i,1))*Xgd1(i)+(eq10(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i)); IGq(i,1)=(-(ed10(i)-VGd(i,1))*Xgd1(i)+(eq10(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i));
@ -435,6 +554,7 @@ if ~isempty(syn)
(Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i)); (Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i));
MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),... MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),...
MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)]; MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)];
% 3rd order% khuang 8 JUL
elseif modSyn(i)==3 elseif modSyn(i)==3
IGd(i,1)=((-VGd(i,1))*Rga(i)+(eq10(i)-VGq(i,1))*Xgq(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq(i)); IGd(i,1)=((-VGd(i,1))*Rga(i)+(eq10(i)-VGq(i,1))*Xgq(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq(i));
IGq(i,1)=(-(-VGd(i,1))*Xgd1(i)+(eq10(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq(i)); IGq(i,1)=(-(-VGd(i,1))*Xgd1(i)+(eq10(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd1(i)*Xgq(i));
@ -443,6 +563,7 @@ if ~isempty(syn)
(Rga(i)*Rga(i)+Xgd1(i)*Xgq(i)); (Rga(i)*Rga(i)+Xgd1(i)*Xgq(i));
MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),... MatsRs(i,:)=[MatsR(i,1)*sind(i)+MatsR(i,2)*cosd(i),-MatsR(i,1)*cosd(i)+MatsR(i,2)*sind(i),...
MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)]; MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)];
% classical model% khuang 8 JUL
elseif modSyn(i)==2 elseif modSyn(i)==2
IGd(i,1)=((-VGd(i,1))*Rga(i)+(Ef0(i)-VGq(i,1))*Xgq(i))/(Rga(i)*Rga(i)+Xgd(i)*Xgq(i)); IGd(i,1)=((-VGd(i,1))*Rga(i)+(Ef0(i)-VGq(i,1))*Xgq(i))/(Rga(i)*Rga(i)+Xgd(i)*Xgq(i));
IGq(i,1)=(-(-VGd(i,1))*Xgd(i)+(Ef0(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd(i)*Xgq(i)); IGq(i,1)=(-(-VGd(i,1))*Xgd(i)+(Ef0(i)-VGq(i,1))*Rga(i))/(Rga(i)*Rga(i)+Xgd(i)*Xgq(i));
@ -453,9 +574,12 @@ if ~isempty(syn)
MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)]; MatsR(i,3)*sind(i)+MatsR(i,4)*cosd(i),-MatsR(i,3)*cosd(i)+MatsR(i,4)*sind(i)];
end end
end end
JG(:,1)=sind.*IGd(:,1)+cosd.*IGq(:,1); % not sure how to use them now, but they are zeroth order of Ix and Iy% khuang 8 JUL
JG(:,1)= sind.*IGd(:,1)+cosd.*IGq(:,1);
KG(:,1)=-cosd.*IGd(:,1)+sind.*IGq(:,1); KG(:,1)=-cosd.*IGd(:,1)+sind.*IGq(:,1);
% put previous matrix in a right place in all buses instead of only% khuang 8 JUL
% generator buses% khuang 8 JUL
MatGCD=-[sparse(synIdx,synIdx,MatsRs(:,1),nbus,nbus),sparse(synIdx,synIdx,MatsRs(:,2),nbus,nbus);... MatGCD=-[sparse(synIdx,synIdx,MatsRs(:,1),nbus,nbus),sparse(synIdx,synIdx,MatsRs(:,2),nbus,nbus);...
sparse(synIdx,synIdx,MatsRs(:,3),nbus,nbus),sparse(synIdx,synIdx,MatsRs(:,4),nbus,nbus)]; sparse(synIdx,synIdx,MatsRs(:,3),nbus,nbus),sparse(synIdx,synIdx,MatsRs(:,4),nbus,nbus)];
else else
@ -478,43 +602,57 @@ else
Ef=zeros(0,nlvl+1); Ef=zeros(0,nlvl+1);
Pm=zeros(0,nlvl+1); Pm=zeros(0,nlvl+1);
end end
%------------------------------Initialization of GEN------------------% khuang 8 JUL
%------------------------------EnD------------------------------------% khuang 8 JUL
%------------------------------Initialization of Exciter------------------% khuang 8 JUL
%------------------------------START------------------------------------% khuang 8 JUL
if ~isempty(exc) if ~isempty(exc)
nExc=size(exc,1); nExc =size(exc,1);
% All Type 3 AVR % All Type 3 AVR, a 3rd order controller
% for Type 3 AVR, avr0(:,1:3) are Vavrm, Vavrr, Vavrf, % for Type 3 AVR, avr0(:,1:3) are Vavrm, Vavrr, Vavrf,
% and avr0(:,4) is reference Vref (input for secondary voltage control). % and avr0(:,4) is reference Vref (input for secondary voltage control).
excIdx=exc(:,1); excIdx = exc(:,1);
VavrMax=exc(:,3); VavrMax = exc(:,3);
VavrMin=exc(:,4); VavrMin = exc(:,4);
muavr0=exc(:,5); muavr0 = exc(:,5);
Tavr1=exc(:,7); Tavr1 = exc(:,7);
Tavr2=exc(:,6); Tavr2 = exc(:,6);
vavrf0=exc(:,8); vavrf0 = exc(:,8);
Vavr0=exc(:,9); Vavr0 = exc(:,9);
Tavre=exc(:,10); Tavre = exc(:,10);
Tavrr=exc(:,11); Tavrr = exc(:,11);
Vavrm=zeros(nExc,nlvl+1); %here I need to check why Vavrref is time varing instead of constant% khuang 8 JUL
Vavrr=zeros(nExc,nlvl+1); % memory is given to state variables of EXC% khuang 8 JUL
Vavrf=zeros(nExc,nlvl+1); Vavrm = zeros(nExc,nlvl+1);
Vavrref=zeros(nExc,nlvl+1); Vavrr = zeros(nExc,nlvl+1);
Vavrf = zeros(nExc,nlvl+1);
Vavrref= zeros(nExc,nlvl+1);
% zeroth order value is given% khuang 8 JUL
Vavrm(:,1)=real(Vavrm0); Vavrm(:,1)=real(Vavrm0);
Vavrr(:,1)=real(Vavrr0); Vavrr(:,1)=real(Vavrr0);
Vavrf(:,1)=real(Vavrf0); Vavrf(:,1)=real(Vavrf0);
Vavrref(:,1)=real(Vavrref0); Vavrref(:,1)=real(Vavrref0);
% here Varref1 is given with syspara.% khuang 8 JUL
if ~isempty(Varref1) if ~isempty(Varref1)
Vavrref(:,2)=Varref1; Vavrref(:,2)=Varref1;
end end
% non-windup limiter, check the limit.% khuang 8 JUL
tavrMaxDiff=Vavrf(:,1)-VavrMax; tavrMaxDiff=Vavrf(:,1)-VavrMax;
tavrMinDiff=Vavrf(:,1)-VavrMin; tavrMinDiff=Vavrf(:,1)-VavrMin;
% label values in different interval.% khuang 8 JUL
avrSt=zeros(nExc,1); avrSt=zeros(nExc,1);
avrSt(tavrMaxDiff>0)=1; avrSt(tavrMaxDiff>0)=1;
avrSt(tavrMinDiff<0)=-1; avrSt(tavrMinDiff<0)=-1;
% output after the limiter.% khuang 8 JUL
Ef(excIdx(avrSt==-1),1)=VavrMin(avrSt==-1); Ef(excIdx(avrSt==-1),1)=VavrMin(avrSt==-1);
Ef(excIdx(avrSt== 1),1)=VavrMax(avrSt== 1); Ef(excIdx(avrSt== 1),1)=VavrMax(avrSt== 1);
Ef(excIdx(avrSt== 0),1)=Vavrf(avrSt==0,1); Ef(excIdx(avrSt== 0),1)=Vavrf(avrSt==0,1);
@ -525,37 +663,48 @@ else
Vavrf=zeros(0,nlvl+1); Vavrf=zeros(0,nlvl+1);
Vavrref=zeros(0,nlvl+1); Vavrref=zeros(0,nlvl+1);
end end
%------------------------------Initialization of Exciter------------------% khuang 8 JUL
%------------------------------END------------------------------------% khuang 8 JUL
%------------------------------Initialization of Turbing Governor------------------% khuang 8 JUL
%------------------------------START------------------------------------% khuang 8 JUL
if ~isempty(tg) if ~isempty(tg)
nTg=size(tg,1); nTg = size(tg,1);
% Type 2 Turbing governor. % Type 2 Turbing governor.
tgIdx=tg(:,1); % one DE, one AE and one limiter
tgIdx = tg(:,1);
wtgref=tg(:,3); wtgref = tg(:,3);
Rtg=tg(:,4); Rtg = tg(:,4);
Ttgmax=tg(:,5); Ttgmax = tg(:,5);
Ttgmin=tg(:,6); Ttgmin = tg(:,6);
Ttg2=tg(:,7); Ttg2 = tg(:,7);
Ttg1=tg(:,8); Ttg1 = tg(:,8);
tgovg=zeros(nTg,nlvl+1); tgovg = zeros(nTg,nlvl+1); % tg % khuang 8 JUL
tgovm=zeros(nTg,nlvl+1); tgovm = zeros(nTg,nlvl+1); % Tmi* % khuang 8 JUL
Tmech=zeros(nTg,nlvl+1); Tmech = zeros(nTg,nlvl+1); % Tmi0 % khuang 8 JUL
% zeroth value is given % khuang 8 JUL
tgovg(:,1)=real(tgovg0); tgovg(:,1)=real(tgovg0);
tgovm(:,1)=real(tgovm0); tgovm(:,1)=real(tgovm0);
Tmech(:,1)=real(tgovmech0); Tmech(:,1)=real(tgovmech0);
if ~isempty(Tmech1) if ~isempty(Tmech1)
Tmech(:,2)=Tmech1; Tmech(:,2)=Tmech1;
end end
% check if limit is approached % khuang 8 JUL
tgovMaxDiff=tgovm(:,1)-Ttgmax; tgovMaxDiff=tgovm(:,1)-Ttgmax;
tgovMinDiff=tgovm(:,1)-Ttgmin; tgovMinDiff=tgovm(:,1)-Ttgmin;
govSt=zeros(nTg,1); govSt=zeros(nTg,1);
govSt(tgovMaxDiff>0)=1; govSt(tgovMaxDiff>0)=1;
govSt(tgovMinDiff<0)=-1; govSt(tgovMinDiff<0)=-1;
% if limit is approached, set Pm to constant value. % khuang 8 JUL
Pm(tgIdx(govSt==0),1)=tgovm(govSt==0,1); Pm(tgIdx(govSt==0),1)=tgovm(govSt==0,1);
Pm(tgIdx(govSt==1),1)=Ttgmax(govSt==1,1); Pm(tgIdx(govSt==1),1)=Ttgmax(govSt==1,1);
Pm(tgIdx(govSt==-1),1)=Ttgmin(govSt==-1,1); Pm(tgIdx(govSt==-1),1)=Ttgmin(govSt==-1,1);
@ -564,7 +713,12 @@ else
tgovm=zeros(0,nlvl+1); tgovm=zeros(0,nlvl+1);
Tmech=zeros(0,nlvl+1); Tmech=zeros(0,nlvl+1);
end end
%------------------------------Initialization of Turbing Governor------------------% khuang 8 JUL
%------------------------------END------------------------------------% khuang 8 JUL
% this part i don't quite understand. It looks like f denotes frequency % khuang 1 JUL
% on every bus, is it relevant with frequency dependant load? % khuang 1 JUL
% now i find that this is for dynamics of agc. % khuang 8 JUL
f=zeros(nbus,nlvl+1); f=zeros(nbus,nlvl+1);
f(:,1)=f0; f(:,1)=f0;
synTag=zeros(nbus,1); synTag=zeros(nbus,1);
@ -581,6 +735,7 @@ for islIdx=1:nIslands
end end
end end
%AGC part % khuang 8 JUL
if ~isempty(agc) if ~isempty(agc)
agcExt=zeros(nbus,size(agc,2)); agcExt=zeros(nbus,size(agc,2));
agcExt(agc(:,1),:)=agc; agcExt(agc(:,1),:)=agc;
@ -592,6 +747,7 @@ else
fdk=zeros(nbus,1); fdk=zeros(nbus,1);
end end
% this is for long term dynamic, it seems that not considered here % khuang 8 JUL
if ~isempty(cac)&&~isempty(cluster) if ~isempty(cac)&&~isempty(cluster)
else else
@ -599,14 +755,16 @@ else
vg=zeros(0,nlvl+1); vg=zeros(0,nlvl+1);
end end
% freq % freq relevant part induced by AGC % khuang 8 JUL
FreqReal=sparse(1:nbus,1:nbus,-freqKeptTag.*fdk.*E0,nbus,nbus); FreqReal=sparse(1:nbus,1:nbus,-freqKeptTag.*fdk.*E0,nbus,nbus);
FreqImag=sparse(1:nbus,1:nbus,-freqKeptTag.*fdk.*F0,nbus,nbus); FreqImag=sparse(1:nbus,1:nbus,-freqKeptTag.*fdk.*F0,nbus,nbus);
Freq2freq=sparse([1:nbus,1:nbus],[1:nbus,frefs(islands)'],[ones(1,nbus),-ones(1,nbus)],nbus,nbus); Freq2freq=sparse([1:nbus,1:nbus],[1:nbus,frefs(islands)'],[ones(1,nbus),-ones(1,nbus)],nbus,nbus);
Y11=-G;Y12=B;Y21=-B;Y22=-G; Y11=-G;Y12=B;Y21=-B;Y22=-G;
% Influence to Origianl Power flow % khuang 8 JUL
YEF11=P0M+sparse(1:nbus,1:nbus,freqKeptTag.*(-fdk.*f0+dpg0),nbus,nbus);YEF12=-Q0M;YEF21=-Q0M;YEF22=-P0M-sparse(1:nbus,1:nbus,freqKeptTag.*(-fdk.*f0+dpg0),nbus,nbus); YEF11=P0M+sparse(1:nbus,1:nbus,freqKeptTag.*(-fdk.*f0+dpg0),nbus,nbus);YEF12=-Q0M;YEF21=-Q0M;YEF22=-P0M-sparse(1:nbus,1:nbus,freqKeptTag.*(-fdk.*f0+dpg0),nbus,nbus);
% Counting influence of ZIP load into Y matrix. % khuang 8 JUL
if ~isempty(zip) if ~isempty(zip)
Y11=Y11-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,1),[nbus,1]),nbus,nbus); 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); Y12=Y12-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,2),[nbus,1]),nbus,nbus);
@ -615,12 +773,14 @@ if ~isempty(zip)
end end
YLHS=[Y11,Y12;Y21,Y22]; YLHS=[Y11,Y12;Y21,Y22];
% Counting influence of Motors into small Y matrix % khuang 8 JUL
if ~isempty(ind) if ~isempty(ind)
YLHS=YLHS-... YLHS=YLHS-...
[sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,1),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,3),nbus,nbus);... [sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,1),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,3),nbus,nbus);...
sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,2),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,4),nbus,nbus)]; sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,2),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,4),nbus,nbus)];
end end
% Counting influence of generators into small Y matrix % khuang 8 JUL
if ~isempty(syn) if ~isempty(syn)
YLHS=YLHS+MatGCD; YLHS=YLHS+MatGCD;
end end
@ -629,6 +789,7 @@ idxNonSw=find(busType~=2);
idxNonSwxD=find(fswTagxD==0); idxNonSwxD=find(fswTagxD==0);
idxNonSwD=find(busType~=2&fswTagxD==1); idxNonSwD=find(busType~=2&fswTagxD==1);
% This is the left hand side matrix totally % khuang 8 JUL
LHS_mat=[YLHS([idxNonSw;idxNonSw+nbus],[idxNonSw;idxNonSw+nbus]),... LHS_mat=[YLHS([idxNonSw;idxNonSw+nbus],[idxNonSw;idxNonSw+nbus]),...
[YEF11(idxNonSw,idxNonSw),YEF12(idxNonSw,idxNonSw),-F0M(idxNonSw,ipv),FreqReal(idxNonSw,freqKeptTag==1);... [YEF11(idxNonSw,idxNonSw),YEF12(idxNonSw,idxNonSw),-F0M(idxNonSw,ipv),FreqReal(idxNonSw,freqKeptTag==1);...
YEF21(idxNonSw,idxNonSw),YEF22(idxNonSw,idxNonSw),-E0M(idxNonSw,ipv),-FreqImag(idxNonSw,freqKeptTag==1)];... YEF21(idxNonSw,idxNonSw),YEF22(idxNonSw,idxNonSw),-E0M(idxNonSw,ipv),-FreqImag(idxNonSw,freqKeptTag==1)];...
@ -642,6 +803,11 @@ LHS_mat=[YLHS([idxNonSw;idxNonSw+nbus],[idxNonSw;idxNonSw+nbus]),...
% [L_LHS_mat,U_LHS_mat,p_LHS_mat]=lu(LHS_mat,'vector'); % [L_LHS_mat,U_LHS_mat,p_LHS_mat]=lu(LHS_mat,'vector');
% end % end
% deterine if we use LU factoration
% for this part, i assume the system algebrac equation is under a good
% condition number and the dimension is not very high, otherwise LU will
% be time consuming % khuang 8 JUL
useLU=isfield(SysPara,'iden')&&isfield(SysPara,'p_amd'); useLU=isfield(SysPara,'iden')&&isfield(SysPara,'p_amd');
if useLU if useLU
@ -660,8 +826,22 @@ if useLU
end end
end end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%this is the recursive part for computing high order of time series%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% khuang 8 JUL
% strat interations nlvl: order of Taylor series. % khuang 8 JUL
for i=1:nlvl for i=1:nlvl
% khuang 8 JUL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% seq2 provides two columns from 0 to i, and i to 0
% seq2p provides two columns from 0 to i+1, and i+1 to 0
% seq3 provides 3 columns, the summary of each row is equal to i(binominal coefficients)
% khuang 8 JUL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
seq2=getseq(i,2); seq2=getseq(i,2);
seq2p=getseq(i+1,2); seq2p=getseq(i+1,2);
seq3=getseq(i,3); seq3=getseq(i,3);
@ -670,6 +850,15 @@ for i=1:nlvl
idxSeq2p=sum(seq2p>=i,2); idxSeq2p=sum(seq2p>=i,2);
idxSeq3=sum(seq3==i,2); idxSeq3=sum(seq3==i,2);
idxSeq3x=sum(seq3(:,[2,3])==i,2); idxSeq3x=sum(seq3(:,[2,3])==i,2);
% khuang 8 JUL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% seq2R is usually used in constructing algebric equations
% seq2R provides two columns from 1 to i-1, and i-1 to 1
% seq2x provides two columns from 1 to i, and i-1 to 0
% seq2m provides two columns from 0 to i-1, and i-1 to 0
% seq2mm provides two columns from 0 to i-2, and i-2 to 0
% khuang 8 JUL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
seq2R=seq2(idxSeq2==0,:); seq2R=seq2(idxSeq2==0,:);
seq2x=seq2(idxSeq2x==0,:); seq2x=seq2(idxSeq2x==0,:);
seq2m=getseq(i-1,2); seq2m=getseq(i-1,2);
@ -677,7 +866,10 @@ for i=1:nlvl
RHSILr=zeros(nbus,1); RHSILr=zeros(nbus,1);
RHSILi=zeros(nbus,1); RHSILi=zeros(nbus,1);
% This part is for induction motor % khuang 8 JUL
if ~isempty(ind) if ~isempty(ind)
% package right hand side vector at every iteration. % khuang 8 JUL
rhsM=sum(Vm(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2)-1j*X2.*sum(IR(:,seq2R(:,1)+1).*s(:,seq2R(:,2)+1),2); 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))+... % 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))./R2+... % (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))./R2+...
@ -692,6 +884,8 @@ for i=1:nlvl
% ./(2*Hm.*s(:,1)*i); % ./(2*Hm.*s(:,1)*i);
% end % end
% update the high order of slip, a special setting is required for
% low order when i<2. % khuang 8 JUL
s(:,i+1)=(Rind0.*(T1.*s(:,i)+T2.*sum(s(:,seq2m(:,1)+1).*s(:,seq2m(:,2)+1),2))-real(sum(Vm(:,seq2m(:,1)+1).*conj(IR(:,seq2m(:,2)+1)),2)))./(2*Hm*i); s(:,i+1)=(Rind0.*(T1.*s(:,i)+T2.*sum(s(:,seq2m(:,1)+1).*s(:,seq2m(:,2)+1),2))-real(sum(Vm(:,seq2m(:,1)+1).*conj(IR(:,seq2m(:,2)+1)),2)))./(2*Hm*i);
if i>=2 if i>=2
s(:,i+1)=s(:,i+1)+... s(:,i+1)=s(:,i+1)+...
@ -704,6 +898,7 @@ for i=1:nlvl
if i==2 if i==2
s(:,i+1)=s(:,i+1)+Rind1.*T0./(2*Hm*i); s(:,i+1)=s(:,i+1)+Rind1.*T0./(2*Hm*i);
end end
% for dynamic model, Right hand side vector is required a update. % khuang 8 JUL
addenRhs=Vm(:,1).*s(:,i+1)-1j*X2.*IR(:,1).*s(:,i+1); addenRhs=Vm(:,1).*s(:,i+1)-1j*X2.*IR(:,1).*s(:,i+1);
% rhsBus=zeros(2,nInd); % rhsBus=zeros(2,nInd);
@ -711,10 +906,12 @@ for i=1:nlvl
% rhsBus(:,j)=RHS_C_Shr{j}*[real(rhsM(j)+addenRhs(j));imag(rhsM(j)+addenRhs(j));0;0]; % rhsBus(:,j)=RHS_C_Shr{j}*[real(rhsM(j)+addenRhs(j));imag(rhsM(j)+addenRhs(j));0;0];
% end % end
% count the influence of dynamic of slip into rigt hand side
% vector.% khuang 8 JUL
tempRhsInd=rhsM+addenRhs; tempRhsInd=rhsM+addenRhs;
rhsBus=[RHS_C_Shr_sqz(:,1).*real(tempRhsInd)+RHS_C_Shr_sqz(:,3).*imag(tempRhsInd),RHS_C_Shr_sqz(:,2).*real(tempRhsInd)+RHS_C_Shr_sqz(:,4).*imag(tempRhsInd)]'; rhsBus=[RHS_C_Shr_sqz(:,1).*real(tempRhsInd)+RHS_C_Shr_sqz(:,3).*imag(tempRhsInd),RHS_C_Shr_sqz(:,2).*real(tempRhsInd)+RHS_C_Shr_sqz(:,4).*imag(tempRhsInd)]';
%accumulate currents %accumulate currents IL.% khuang 8 JUL
RHSILr=accumarray(indIdx,rhsBus(1,:)',[nbus,1]); RHSILr=accumarray(indIdx,rhsBus(1,:)',[nbus,1]);
RHSILi=accumarray(indIdx,rhsBus(2,:)',[nbus,1]); RHSILi=accumarray(indIdx,rhsBus(2,:)',[nbus,1]);
@ -734,6 +931,7 @@ for i=1:nlvl
% RHSILi=accumarray(indIdx,rhsBus(4,:)',[nbus,1]); % RHSILi=accumarray(indIdx,rhsBus(4,:)',[nbus,1]);
end end
% strat update ZIP load into currents.% khuang 8 JUL
RHSIiLr=zeros(nbus,1); RHSIiLr=zeros(nbus,1);
RHSIiLi=zeros(nbus,1); RHSIiLi=zeros(nbus,1);
if ~isempty(zip) if ~isempty(zip)
@ -745,6 +943,7 @@ for i=1:nlvl
RHSIiLi=accumarray(zipIdx,RHSILi_full,[nbus,1]); RHSIiLi=accumarray(zipIdx,RHSILi_full,[nbus,1]);
end end
% Start update generators.% khuang 8 JUL
RHSIGr=zeros(nbus,1); RHSIGr=zeros(nbus,1);
RHSIGi=zeros(nbus,1); RHSIGi=zeros(nbus,1);
if ~isempty(syn) if ~isempty(syn)
@ -752,7 +951,7 @@ for i=1:nlvl
RhsEq=zeros(nSyn,1); RhsEq=zeros(nSyn,1);
IGdAdd=zeros(nSyn,1); IGdAdd=zeros(nSyn,1);
IGqAdd=zeros(nSyn,1); IGqAdd=zeros(nSyn,1);
% select different models for generators.% khuang 8 JUL
if modelTag(8)>0 if modelTag(8)>0
d(modSyn==8,i+1)=(wgb(modSyn==8).*w(modSyn==8,i))/i; d(modSyn==8,i+1)=(wgb(modSyn==8).*w(modSyn==8,i))/i;
w(modSyn==8,i+1)=(Pm(modSyn==8,i)-... w(modSyn==8,i+1)=(Pm(modSyn==8,i)-...
@ -822,10 +1021,11 @@ for i=1:nlvl
RhsEd(modSyn==2)=0; RhsEd(modSyn==2)=0;
RhsEq(modSyn==2)=eq1(modSyn==2,i+1); RhsEq(modSyn==2)=eq1(modSyn==2,i+1);
end end
% this part may be different from DT.% khuang 8 JUL
AG0=cosp(:,2).*d(:,i+1); AG0=cosp(:,2).*d(:,i+1);
BG0=sinp(:,2).*d(:,i+1); BG0=sinp(:,2).*d(:,i+1);
% here multi-convolution is utilized as sine function is
% approxiamted as a taylor series of delta.% khuang 8 JUL
if taylorN>=2 if taylorN>=2
AG0=AG0+cosp(:,3).*sum(d(:,seq2(:,1)+1).*d(:,seq2(:,2)+1),2); AG0=AG0+cosp(:,3).*sum(d(:,seq2(:,1)+1).*d(:,seq2(:,2)+1),2);
BG0=BG0+sinp(:,3).*sum(d(:,seq2(:,1)+1).*d(:,seq2(:,2)+1),2); BG0=BG0+sinp(:,3).*sum(d(:,seq2(:,1)+1).*d(:,seq2(:,2)+1),2);
@ -840,14 +1040,19 @@ for i=1:nlvl
BG0=BG0+sinp(:,5).*sum(d(:,seq4(:,1)+1).*d(:,seq4(:,2)+1).*d(:,seq4(:,3)+1).*d(:,seq4(:,4)+1),2); BG0=BG0+sinp(:,5).*sum(d(:,seq4(:,1)+1).*d(:,seq4(:,2)+1).*d(:,seq4(:,3)+1).*d(:,seq4(:,4)+1),2);
end end
Cd(:,i+1)=AG0; % now Sd store high order terms of sin(dta) and Cd for cos(dta),
Sd(:,i+1)=BG0; % the following sentence is for DT, i commentde it for HE. % khuang 8 JUL
%Sd(:,i+1) = 1/(i)*sum(repmat([i:-1:1],nSyn,1).*Cd(:,1:i).*d(:,i+1:-1:2),2);
%Cd(:,i+1) =-1/(i)*sum(repmat([i:-1:1],nSyn,1).*Sd(:,1:i).*d(:,i+1:-1:2),2);
% high order coefficients of cos(delta) and sin(delta).% khuang 8 JUL
Cd(:,i+1)=AG0;
Sd(:,i+1)=BG0;
VGdCr=sum(Cd(:,seq2x(:,1)+1).*VGd(:,seq2x(:,2)+1),2); VGdCr=sum(Cd(:,seq2x(:,1)+1).*VGd(:,seq2x(:,2)+1),2);% Vd*cosdta% khuang 8 JUL
VGqCr=sum(Cd(:,seq2x(:,1)+1).*VGq(:,seq2x(:,2)+1),2); VGqCr=sum(Cd(:,seq2x(:,1)+1).*VGq(:,seq2x(:,2)+1),2);% Vq*cosdta% khuang 8 JUL
VGdSr=sum(Sd(:,seq2x(:,1)+1).*VGd(:,seq2x(:,2)+1),2); VGdSr=sum(Sd(:,seq2x(:,1)+1).*VGd(:,seq2x(:,2)+1),2);% Vd*sindta% khuang 8 JUL
VGqSr=sum(Sd(:,seq2x(:,1)+1).*VGq(:,seq2x(:,2)+1),2); VGqSr=sum(Sd(:,seq2x(:,1)+1).*VGq(:,seq2x(:,2)+1),2);% Vq*sindta% khuang 8 JUL
JCr=sum(Cd(:,seq2x(:,1)+1).*JG(:,seq2x(:,2)+1),2); JCr=sum(Cd(:,seq2x(:,1)+1).*JG(:,seq2x(:,2)+1),2);% similar, for currents% khuang 8 JUL
KCr=sum(Cd(:,seq2x(:,1)+1).*KG(:,seq2x(:,2)+1),2); KCr=sum(Cd(:,seq2x(:,1)+1).*KG(:,seq2x(:,2)+1),2);
JSr=sum(Sd(:,seq2x(:,1)+1).*JG(:,seq2x(:,2)+1),2); JSr=sum(Sd(:,seq2x(:,1)+1).*JG(:,seq2x(:,2)+1),2);
KSr=sum(Sd(:,seq2x(:,1)+1).*KG(:,seq2x(:,2)+1),2); KSr=sum(Sd(:,seq2x(:,1)+1).*KG(:,seq2x(:,2)+1),2);
@ -856,10 +1061,11 @@ for i=1:nlvl
(MatsR(:,1).*RhsEd+MatsR(:,2).*RhsEq)-(Mats(:,1).*(JSr-KCr)+Mats(:,2).*(JCr+KSr))+(Mats(:,1).*IGdAdd+Mats(:,2).*IGqAdd); (MatsR(:,1).*RhsEd+MatsR(:,2).*RhsEq)-(Mats(:,1).*(JSr-KCr)+Mats(:,2).*(JCr+KSr))+(Mats(:,1).*IGdAdd+Mats(:,2).*IGqAdd);
RHSIGxi=-(MatsRs(:,3).*(-VGdSr-VGqCr)+MatsRs(:,4).*(VGdCr-VGqSr))+... RHSIGxi=-(MatsRs(:,3).*(-VGdSr-VGqCr)+MatsRs(:,4).*(VGdCr-VGqSr))+...
(MatsR(:,3).*RhsEd+MatsR(:,4).*RhsEq)-(Mats(:,3).*(JSr-KCr)+Mats(:,4).*(JCr+KSr))+(Mats(:,3).*IGdAdd+Mats(:,4).*IGqAdd); (MatsR(:,3).*RhsEd+MatsR(:,4).*RhsEq)-(Mats(:,3).*(JSr-KCr)+Mats(:,4).*(JCr+KSr))+(Mats(:,3).*IGdAdd+Mats(:,4).*IGqAdd);
% current injections from generators IG.% khuang 8 JUL
RHSIGr=accumarray(synIdx,RHSIGxr,[nbus,1]); RHSIGr=accumarray(synIdx,RHSIGxr,[nbus,1]);
RHSIGi=accumarray(synIdx,RHSIGxi,[nbus,1]); RHSIGi=accumarray(synIdx,RHSIGxi,[nbus,1]);
end end
% update exciter, 3 state variables.% khuang 8 JUL
if ~isempty(exc) if ~isempty(exc)
Vavrm(:,i+1)=(Vmag(synIdx(excIdx),i)-Vavrm(:,i))./Tavrr/i; Vavrm(:,i+1)=(Vmag(synIdx(excIdx),i)-Vavrm(:,i))./Tavrr/i;
Vavrr(:,i+1)=(muavr0.*(1-Tavr1./Tavr2).*(Vavrref(:,i)-Vavrm(:,i))-Vavrr(:,i))./Tavr2/i; Vavrr(:,i+1)=(muavr0.*(1-Tavr1./Tavr2).*(Vavrref(:,i)-Vavrm(:,i))-Vavrr(:,i))./Tavr2/i;
@ -871,6 +1077,7 @@ for i=1:nlvl
Ef(excIdx(avrSt== 0),i+1)=Vavrf(avrSt==0,i+1); Ef(excIdx(avrSt== 0),i+1)=Vavrf(avrSt==0,i+1);
end end
% update agc, one state variables.% khuang 8 JUL
if ~isempty(agc) if ~isempty(agc)
dpg(:,i+1)=-f(:,i).*agcExt(:,4)/i; dpg(:,i+1)=-f(:,i).*agcExt(:,4)/i;
for islIdx=1:nIslands for islIdx=1:nIslands
@ -882,6 +1089,7 @@ for i=1:nlvl
end end
end % TODO: steady-state model end % TODO: steady-state model
% update generator participation part from agc.% khuang 8 JUL
if ~isempty(syn) %dynamic model (synchronous generators) if ~isempty(syn) %dynamic model (synchronous generators)
if ~isempty(tg) if ~isempty(tg)
Tmech(:,i+1)=Tmech(:,i+1)+dpg(syn(tg(:,1),1),i+1)./numSynOnBus(syn(tg(:,1),1)); Tmech(:,i+1)=Tmech(:,i+1)+dpg(syn(tg(:,1),1),i+1)./numSynOnBus(syn(tg(:,1),1));
@ -889,7 +1097,7 @@ for i=1:nlvl
Pm(:,i+1)=Pm(:,i+1)+dpg(syn(:,1),i+1)./numSynOnBus(syn(:,1)); Pm(:,i+1)=Pm(:,i+1)+dpg(syn(:,1),i+1)./numSynOnBus(syn(:,1));
end end
end end
% update Turbine, 2 state variables.% khuang 8 JUL
if ~isempty(tg) if ~isempty(tg)
tgovg(:,i+1)=(-(1-Ttg1./Ttg2).*w(tgIdx,i)./Rtg-tgovg(:,i))./Ttg2/i; tgovg(:,i+1)=(-(1-Ttg1./Ttg2).*w(tgIdx,i)./Rtg-tgovg(:,i))./Ttg2/i;
tgovm(:,i+1)=tgovg(:,i+1)-Ttg1./Ttg2.*w(tgIdx,i+1)./Rtg+Tmech(:,i+1); tgovm(:,i+1)=tgovg(:,i+1)-Ttg1./Ttg2.*w(tgIdx,i+1)./Rtg+Tmech(:,i+1);
@ -906,12 +1114,15 @@ for i=1:nlvl
RHS2=-0.5*real(sum(V(:,seq2R(:,1)+1).*conj(V(:,seq2R(:,2)+1)),2)); 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); RHS3=sum(-W(:,seq2R(:,1)+1).*V(:,seq2R(:,2)+1),2);
if i==1 if i==1
RHS2=RHS2+0.5*VspSq2(:,2); RHS2=RHS2+0.5*VspSq2(:,2);
end end
compactRHS1=RHS1(busType~=2); compactRHS1=RHS1(busType~=2);
compactRHS1=compactRHS1+Y(busType~=2,isw)*real(V(isw,i+1)); compactRHS1=compactRHS1+Y(busType~=2,isw)*real(V(isw,i+1));
% combine all current injection involing Motor, zip load, and
% Generators.% khuang 8 JUL
RHS=[real(compactRHS1)+RHSILr(busType~=2)+RHSIiLr(busType~=2)-RHSIGr(busType~=2);... RHS=[real(compactRHS1)+RHSILr(busType~=2)+RHSIiLr(busType~=2)-RHSIGr(busType~=2);...
imag(compactRHS1)+RHSILi(busType~=2)+RHSIiLi(busType~=2)-RHSIGi(busType~=2);... imag(compactRHS1)+RHSILi(busType~=2)+RHSIiLi(busType~=2)-RHSIGi(busType~=2);...
RHS2(ipv);... RHS2(ipv);...
@ -919,7 +1130,9 @@ for i=1:nlvl
imag(RHS3(busType~=2));... imag(RHS3(busType~=2));...
zeros(sum(freqKeptTagxRef),1);... zeros(sum(freqKeptTagxRef),1);...
zeros(size(idxNonSwD,1),1)]; zeros(size(idxNonSwD,1),1)];
% solve AE, notice that every time we need to solve Ax(k) =b(k), which
% means that A in invariant for every order. so we only need to rebulid
% b every iteration.% khuang 8 JUL
if useLU if useLU
if IS_OCTAVE if IS_OCTAVE
x = real(MxQ * MxQx* (MxU \ (MxL \ (MxP * RHS)))) ; x = real(MxQ * MxQx* (MxU \ (MxL \ (MxP * RHS)))) ;
@ -930,6 +1143,7 @@ for i=1:nlvl
x=real(LHS_mat\RHS); x=real(LHS_mat\RHS);
end end
% x= [V;W;Q_pv;f]
xC=real(V(:,i+1)); xC=real(V(:,i+1));
xD=imag(V(:,i+1)); xD=imag(V(:,i+1));
xC(idxNonSw)=x(1:(npq+npv)); xC(idxNonSw)=x(1:(npq+npv));
@ -942,6 +1156,7 @@ for i=1:nlvl
Vmag(:,i+1)=(sum(V(:,seq2(:,1)+1).*conj(V(:,seq2(:,2)+1)),2)-sum(Vmag(:,seq2R(:,1)+1).*Vmag(:,seq2R(:,2)+1),2))./Vmag(:,1)/2; % Calculate voltage magnitude Vmag(:,i+1)=(sum(V(:,seq2(:,1)+1).*conj(V(:,seq2(:,2)+1)),2)-sum(Vmag(:,seq2R(:,1)+1).*Vmag(:,seq2R(:,2)+1),2))./Vmag(:,1)/2; % Calculate voltage magnitude
% now update the Algebric variables for motors:IL,IR,VM.% khuang 8 JUL
if ~isempty(ind) if ~isempty(ind)
% for j=1:nInd % for j=1:nInd
% tempIL=squeeze(LHS_MatInd_Shr(j,:,:))*[real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+rhsBus(:,j); % tempIL=squeeze(LHS_MatInd_Shr(j,:,:))*[real(V(indIdx(j),i+1));imag(V(indIdx(j),i+1))]+rhsBus(:,j);
@ -959,11 +1174,13 @@ for i=1:nlvl
Vm(:,i+1)=V(indIdx,i+1)-IL(:,i+1).*Z1; Vm(:,i+1)=V(indIdx,i+1)-IL(:,i+1).*Z1;
end end
% now update the Algebric variables for ZIP loads.% khuang 8 JUL
if ~isempty(zip) 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); 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; BiL(:,i+1)=Mat_BZip(:,1).*real(V(zipIdx,i+1))+Mat_BZip(:,2).*imag(V(zipIdx,i+1))+RHS_BZip;
end end
% now update the Algebric variables for Generators: Vd,Vq,Id,Iq.% khuang 8 JUL
if ~isempty(syn) if ~isempty(syn)
JG(:,i+1)=-MatsRs(:,1).*real(V(synIdx,i+1))-MatsRs(:,2).*imag(V(synIdx,i+1))+RHSIGxr; JG(:,i+1)=-MatsRs(:,1).*real(V(synIdx,i+1))-MatsRs(:,2).*imag(V(synIdx,i+1))+RHSIGxr;
KG(:,i+1)=-MatsRs(:,3).*real(V(synIdx,i+1))-MatsRs(:,4).*imag(V(synIdx,i+1))+RHSIGxi; KG(:,i+1)=-MatsRs(:,3).*real(V(synIdx,i+1))-MatsRs(:,4).*imag(V(synIdx,i+1))+RHSIGxi;
@ -976,6 +1193,7 @@ for i=1:nlvl
end end
end end
% Output value: coefficients for every order.% khuang 8 JUL
Q=real(Q); Q=real(Q);
s=real(s); s=real(s);
d=real(d); d=real(d);

4
internal/hemMachinePFmultiStage.m
Unescape View File

@ -92,8 +92,8 @@ while alphaConfirm<1-alphaTol/1000
x=zeros(nState,1); x=zeros(nState,1);
x(idxs.vIdx)=V0x;x(idxs.sIdx)=s0x;x(idxs.deltaIdx)=d0x;x(idxs.qIdx)=Q0x;x(idxs.efIdx)=Ef0x;x(idxs.pgIdx)=Pm0x; x(idxs.vIdx)=V0x;x(idxs.sIdx)=s0x;x(idxs.deltaIdx)=d0x;x(idxs.qIdx)=Q0x;x(idxs.efIdx)=Ef0x;x(idxs.pgIdx)=Pm0x;
SysDatax=foldSysData(bus,sw,pvx,pqx,shuntx,line,indx,zipx,syn,exc,tg,agc,cac,cluster); SysDatax =foldSysData(bus,sw,pvx,pqx,shuntx,line,indx,zipx,syn,exc,tg,agc,cac,cluster);
SysParax=foldSysPara(pqIncr,pvIncr,Rind0x,Rind1,Reind0x,Reind1,Rzip0x,Rzip1,Ytr,[],Ysh0x,Ysh1,[Vsp0x,VspSq2(:,2)],[],[],[],[],[],[],Ef1,Pm1); SysParax =foldSysPara(pqIncr,pvIncr,Rind0x,Rind1,Reind0x,Reind1,Rzip0x,Rzip1,Ytr,[],Ysh0x,Ysh1,[Vsp0x,VspSq2(:,2)],[],[],[],[],[],[],Ef1,Pm1);
SysParax.iden=iden; SysParax.iden=iden;
if exist([iden,'.mat'],'file') if exist([iden,'.mat'],'file')
if ~exist('p_amd','var') if ~exist('p_amd','var')

108
internal/hemMachinePFmultiStageDyn.m
Unescape View File

@ -146,6 +146,7 @@ nZip=size(zip,1);
% psiqc=zeros(nSyn,1);dpsiqc=zeros(nSyn,1); % psiqc=zeros(nSyn,1);dpsiqc=zeros(nSyn,1);
alphaList=0; alphaList=0;
ratiob4 = NaN;
diffList=0; diffList=0;
diffMul=1.05; diffMul=1.05;
alphaPreMult=0.91; alphaPreMult=0.91;
@ -198,7 +199,7 @@ if ~isempty(varOpt) && isfield(varOpt,'useNomoveCtrl')
else else
useNomoveCtrl=1; useNomoveCtrl=1;
end end
alphab4 = NaN;
while alphaConfirm<maxAlpha-alphaTol/1000 while alphaConfirm<maxAlpha-alphaTol/1000
alpha=min([segAlpha,maxAlpha-alphaConfirm]); alpha=min([segAlpha,maxAlpha-alphaConfirm]);
alphax=alpha/alphaPreMult; alphax=alpha/alphaPreMult;
@ -263,8 +264,12 @@ while alphaConfirm<maxAlpha-alphaTol/1000
SysParax.refs=refs; SysParax.refs=refs;
[V,Q,s,d,w,eq1,eq2,ed1,ed2,psid,psiq,Pm,Ef,Vavrm,Vavrr,Vavrf,Vavrref,tgovg,tgovm,Tmech,f,dpg,qplt,vg]=... [V,Q,s,d,w,eq1,eq2,ed1,ed2,psid,psiq,Pm,Ef,Vavrm,Vavrr,Vavrf,Vavrref,tgovg,tgovm,Tmech,f,dpg,qplt,vg]=...
hemMachinePFSalientcontinueDyn(SimData,SysDatax,SysParax,xx); hemMachinePFSalientcontinueDyn(SimData,SysDatax,SysParax,xx);
%check the dirivative of every state variables w.r.t. time
ds=s(:,2:end).*repmat(1:(size(s,2)-1),size(s,1),1); ds=s(:,2:end).*repmat(1:(size(s,2)-1),size(s,1),1);
dd=d(:,2:end).*repmat(1:(size(d,2)-1),size(d,1),1); dd=d(:,2:end).*repmat(1:(size(d,2)-1),size(d,1),1);
dw=w(:,2:end).*repmat(1:(size(w,2)-1),size(w,1),1); dw=w(:,2:end).*repmat(1:(size(w,2)-1),size(w,1),1);
@ -285,7 +290,8 @@ while alphaConfirm<maxAlpha-alphaTol/1000
df=f(:,2:end).*repmat(1:(size(f,2)-1),size(f,1),1); %agc df=f(:,2:end).*repmat(1:(size(f,2)-1),size(f,1),1); %agc
ddpg=dpg(:,2:end).*repmat(1:(size(dpg,2)-1),size(dpg,1),1); %agc ddpg=dpg(:,2:end).*repmat(1:(size(dpg,2)-1),size(dpg,1),1); %agc
% TODO: avc % TODO: avc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main part of pade
[Va,Vb]=getPadeApproxHelm(V,0); [Va,Vb]=getPadeApproxHelm(V,0);
idxNanV=find(isnan(Va(:,1))|isnan(Vb(:,1)));Va(idxNanV,:)=0;Vb(idxNanV,:)=0;Va(idxNanV,1:3)=V(idxNanV,1:3); idxNanV=find(isnan(Va(:,1))|isnan(Vb(:,1)));Va(idxNanV,:)=0;Vb(idxNanV,:)=0;Va(idxNanV,1:3)=V(idxNanV,1:3);
[Qa,Qb]=getPadeApproxHelm(Q,find(busType==0)); [Qa,Qb]=getPadeApproxHelm(Q,find(busType==0));
@ -333,7 +339,8 @@ while alphaConfirm<maxAlpha-alphaTol/1000
[qplta,qpltb]=getPadeApproxHelm(qplt,0); [qplta,qpltb]=getPadeApproxHelm(qplt,0);
[vga,vgb]=getPadeApproxHelm(vg,0); [vga,vgb]=getPadeApproxHelm(vg,0);
% main part of pade
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(ind) if ~isempty(ind)
idxNanS=find(isnan(sa(:,1))|isnan(sb(:,1)));sa(idxNanS,:)=0;sb(idxNanS,:)=0;sa(idxNanS,1:2)=s(idxNanS,1:2); idxNanS=find(isnan(sa(:,1))|isnan(sb(:,1)));sa(idxNanS,:)=0;sb(idxNanS,:)=0;sa(idxNanS,1:2)=s(idxNanS,1:2);
idxNandS=find(isnan(dsa(:,1))|isnan(dsb(:,1)));dsa(idxNandS,:)=0;dsb(idxNandS,:)=0;dsa(idxNandS,1:2)=ds(idxNandS,1:2); idxNandS=find(isnan(dsa(:,1))|isnan(dsb(:,1)));dsa(idxNandS,:)=0;dsb(idxNandS,:)=0;dsa(idxNandS,1:2)=ds(idxNandS,1:2);
@ -385,11 +392,53 @@ while alphaConfirm<maxAlpha-alphaTol/1000
end end
alphaLeft=0; alphaLeft=0;
alphaRight=alpha/alphaPreMult; alphaRight=alpha/alphaPreMult;
%khuangJune14
% here we utilize PI control strategy instead of bisection search.
% local truncation error is required from all state variable
%-------------------------------------------------------------------------------------------
local_x=foldX(SysDatax,V,Q,s,d,w,eq1,eq2,ed1,ed2,psid,psiq,Pm,Ef,Vavrm,Vavrr,Vavrr,Vavrref,tgovg,tgovm,Tmech,f,dpg,qplt,vg);
% if alphax>0
% if size(alphaList,2)>1&&~isnan(alphab4)
%
% [K,h,~] =VSOO(LTE,LTE2,LTE3,SimData,alphab4,ratiob4);
% % if isnan(h)
% % h
% % end
% % SimData.nlvl = K;
% % alphax = h;
% % alphaRight = h ;
% % alphaLeft = 0.8*h;
%
% % K
% % h
% end
% else
% alphaRight = 0;
% end
% % % LTE=0;
% %------------------------------------------------------------------------------------------
%
%
% alphax=0.2; % alphax=0.2;
while 1 while 1
%pade
vc=polyvalVec(Va(:,end:-1:1),alphax)./polyvalVec([Vb(:,end:-1:1),ones(size(Vb,1),1)],alphax); vc=polyvalVec(Va(:,end:-1:1),alphax)./polyvalVec([Vb(:,end:-1:1),ones(size(Vb,1),1)],alphax);
qc=polyvalVec(Qa(:,end:-1:1),alphax)./polyvalVec([Qb(:,end:-1:1),ones(size(Qb,1),1)],alphax); qc=polyvalVec(Qa(:,end:-1:1),alphax)./polyvalVec([Qb(:,end:-1:1),ones(size(Qb,1),1)],alphax);
@ -439,6 +488,7 @@ while alphaConfirm<maxAlpha-alphaTol/1000
qpltc=polyvalVec(qplta(:,end:-1:1),alphax)./polyvalVec([qpltb(:,end:-1:1),ones(size(qplt,1),1)],alphax); qpltc=polyvalVec(qplta(:,end:-1:1),alphax)./polyvalVec([qpltb(:,end:-1:1),ones(size(qplt,1),1)],alphax);
vgc=polyvalVec(vga(:,end:-1:1),alphax)./polyvalVec([vgb(:,end:-1:1),ones(size(vg,1),1)],alphax); vgc=polyvalVec(vga(:,end:-1:1),alphax)./polyvalVec([vgb(:,end:-1:1),ones(size(vg,1),1)],alphax);
% Y=Ytr0x+alphax*Ytr1+sparse(1:nbus,1:nbus,Ysh0x+Ysh1*alphax,nbus,nbus); % Y=Ytr0x+alphax*Ytr1+sparse(1:nbus,1:nbus,Ysh0x+Ysh1*alphax,nbus,nbus);
%pade b4
Vsp2=Vsp0x+alphax*VspSq2(:,2); Vsp2=Vsp0x+alphax*VspSq2(:,2);
pqxx=pqx; pqxx=pqx;
@ -471,8 +521,9 @@ while alphaConfirm<maxAlpha-alphaTol/1000
diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xc,dxc); diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xc,dxc);
absDiff=max(abs(diff)); absDiff=max(abs(diff));
%bisection search
if absDiff<diffTol&&alphax==alpha/alphaPreMult % if absDiff<diffTol&&alphax==alpha/alphaPreMult
if absDiff<diffTol&&alphax==alpha/alphaPreMult
alphaLeft=alphaRight; alphaLeft=alphaRight;
else else
if absDiff<diffTol if absDiff<diffTol
@ -489,14 +540,35 @@ while alphaConfirm<maxAlpha-alphaTol/1000
end end
end end
diffParadigm=absDiff+1e-12; diffParadigm=absDiff+1e-12;
countCheckAlpha=1; countCheckAlpha=1;
alphax=alphaLeft*alphaPreMult; alphax=alphaLeft*alphaPreMult;
if alpha-alphax<alphaTol/1000 if alpha-alphax<alphaTol/1000
alphax=alpha; alphax=alpha;
end end
maxCount=10; maxCount=10;
while countCheckAlpha<maxCount
while countCheckAlpha<maxCount
%pade
vc=polyvalVec(Va(:,end:-1:1),alphax)./polyvalVec([Vb(:,end:-1:1),ones(size(Vb,1),1)],alphax); vc=polyvalVec(Va(:,end:-1:1),alphax)./polyvalVec([Vb(:,end:-1:1),ones(size(Vb,1),1)],alphax);
qc=polyvalVec(Qa(:,end:-1:1),alphax)./polyvalVec([Qb(:,end:-1:1),ones(size(Qb,1),1)],alphax); qc=polyvalVec(Qa(:,end:-1:1),alphax)./polyvalVec([Qb(:,end:-1:1),ones(size(Qb,1),1)],alphax);
@ -545,6 +617,7 @@ while alphaConfirm<maxAlpha-alphaTol/1000
qpltc=polyvalVec(qplta(:,end:-1:1),alphax)./polyvalVec([qpltb(:,end:-1:1),ones(size(qplt,1),1)],alphax); qpltc=polyvalVec(qplta(:,end:-1:1),alphax)./polyvalVec([qpltb(:,end:-1:1),ones(size(qplt,1),1)],alphax);
vgc=polyvalVec(vga(:,end:-1:1),alphax)./polyvalVec([vgb(:,end:-1:1),ones(size(vg,1),1)],alphax); vgc=polyvalVec(vga(:,end:-1:1),alphax)./polyvalVec([vgb(:,end:-1:1),ones(size(vg,1),1)],alphax);
% pade b4
% Y=Ytr0x+alphax*Ytr1+sparse(1:nbus,1:nbus,Ysh0x+Ysh1*alphax,nbus,nbus); % Y=Ytr0x+alphax*Ytr1+sparse(1:nbus,1:nbus,Ysh0x+Ysh1*alphax,nbus,nbus);
Vsp2=Vsp0x+alphax*VspSq2(:,2); Vsp2=Vsp0x+alphax*VspSq2(:,2);
@ -578,13 +651,13 @@ while alphaConfirm<maxAlpha-alphaTol/1000
SysParaxx.refs=refs; SysParaxx.refs=refs;
diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xc,dxc); diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xc,dxc);
if max(abs(diff))<1.05*diffParadigm if max(abs(diff))<1.05*diffParadigm
break; break;
end end
alphax=alphax*alphaPreMult; alphax=alphax*alphaPreMult;
countCheckAlpha=countCheckAlpha+1; countCheckAlpha=countCheckAlpha+1;
end end
if countCheckAlpha>=maxCount;alphax=alphax/alphaPreMult;end if countCheckAlpha>=maxCount;alphax=alphax/alphaPreMult;end
alpha=alphax; alpha=alphax;
@ -606,7 +679,7 @@ while alphaConfirm<maxAlpha-alphaTol/1000
else else
strPre=['Step=',num2str(alphaConfirm,'%7.5f'),', ']; strPre=['Step=',num2str(alphaConfirm,'%7.5f'),', '];
end end
addLog([strPre,'dt=',num2str(alpha,'%7.5f'),', (maxdiff<',num2str(diffTol),').'],'INFO'); addLog([strPre,'dt=',num2str(alpha,'%7.5f'),', K=',num2str(SimData.nlvl,'%7.1f'),', (maxdiff<',num2str(diffTol),').'],'INFO');
% addLog(['Fluctuation=',num2str(max(fluctuation)),' rate=',num2str(max(fluctuation)/maxTcoeff),'.'],'INFO'); % addLog(['Fluctuation=',num2str(max(fluctuation)),' rate=',num2str(max(fluctuation)/maxTcoeff),'.'],'INFO');
if alpha==0 if alpha==0
@ -828,8 +901,13 @@ while alphaConfirm<maxAlpha-alphaTol/1000
Vavrref0x=Vavrrefcx(:,end);tgovg0x=tgovgcx(:,end);tgovm0x=tgovmcx(:,end);Tmech0x=Tmechcx(:,end); Vavrref0x=Vavrrefcx(:,end);tgovg0x=tgovgcx(:,end);tgovm0x=tgovmcx(:,end);Tmech0x=Tmechcx(:,end);
f0x=fcx(:,end);dpg0x=dpgcx(:,end);qplt0x=qpltcx(:,end);vg0x=vgcx(:,end); f0x=fcx(:,end);dpg0x=dpgcx(:,end);qplt0x=qpltcx(:,end);vg0x=vgcx(:,end);
alphaList=[alphaList,alpha];
diffList=[diffList,max(abs(diff))];
% [LTE,LTE2,LTE3] = Cal_locl_er(local_x,SimData.nlvl,alpha);
% ratiob4 = LTE./alpha;
% alphab4 = alpha;
% alphaList=[alphaList,alpha];
% diffList=[diffList,max(abs(diff))];
if allowModelSwitch if allowModelSwitch
if exitFlag==DynSimFlag.STEADY||exitFlag==DynSimFlag.QSS if exitFlag==DynSimFlag.STEADY||exitFlag==DynSimFlag.QSS

3
internal/readSimulationSettings.m
Unescape View File

@ -49,6 +49,9 @@ try
if exist('diffTolMax','var');simSettings.diffTolMax=diffTolMax;end if exist('diffTolMax','var');simSettings.diffTolMax=diffTolMax;end
if exist('method','var');simSettings.method=method;end if exist('method','var');simSettings.method=method;end
if exist('diffTolCtrl','var');simSettings.diffTolCtrl=diffTolCtrl;end if exist('diffTolCtrl','var');simSettings.diffTolCtrl=diffTolCtrl;end
%%%%%%khuang Sep 5
if exist('dAlpha','var');simSettings.dAlpha=dAlpha;end
%%%%%%%%%%%%%%%%
simSettings.eventList=eventList; simSettings.eventList=eventList;
simSettings.bsSyn=bsSyn; simSettings.bsSyn=bsSyn;
simSettings.bsBus=bsBus; simSettings.bsBus=bsBus;

6
internal/regulateSystemData.m
Unescape View File

@ -127,9 +127,9 @@ if ~isempty(syn)
if ~isempty(pv) if ~isempty(pv)
pg(pv(:,1))=pv(:,4); pg(pv(:,1))=pv(:,4);
end end
% if ~isempty(sw) % if ~isempty(sw)
% pg(sw(:,1))=sw(:,10); % pg(sw(:,1))=sw(:,10);
% end % end
pm=pg(syn(:,1)); pm=pg(syn(:,1));
nm=accumarray(syn(:,1),ones(size(syn,1),1),[nbus,1]); nm=accumarray(syn(:,1),ones(size(syn,1),1),[nbus,1]);
pm=pm./nm(syn(:,1)); pm=pm./nm(syn(:,1));

2
internal/setting_func_simulationTimeDomainNI.m
Unescape View File

@ -29,5 +29,5 @@
maxStep=0.02; maxStep=0.02;
minStep=5e-5; minStep=5e-5;
stepCntThreshold=20; stepCntThreshold=5;
stepSizeEps=1e-9; stepSizeEps=1e-9;

2
internal/simulationTimeDomainHem.m
Unescape View File

@ -52,7 +52,7 @@ function [stateCurve,t,finalAlpha,alphaList,diffList,exitFlag]=simulationTimeDom
% 1 - Modified Euler % 1 - Modified Euler
% 2 - RK4 % 2 - RK4
% 3 - Trapezoidal % 3 - Trapezoidal
% % 4 - DT
% Y - Method for solving algebraic equations (AE) % Y - Method for solving algebraic equations (AE)
% 0 - HE % 0 - HE
% 1 - NR % 1 - NR

11
internal/simulationTimeDomainNI.m
Unescape View File

@ -75,9 +75,9 @@ else
end end
[~,~,~,~,~,~,~,~,method]=unfoldSimData(SimData); [~,~,~,~,~,~,~,~,method]=unfoldSimData(SimData);
DEmethod=floor(method); DEmethod =floor(method);
AEmethod=round(10*mod(method,1)); AEmethod =round(10*mod(method,1));
varStep=round(100*mod(method,0.1)); varStep =round(100*mod(method,0.1));
eval('setting_func_simulationTimeDomainNI'); eval('setting_func_simulationTimeDomainNI');
if varStep~=0 if varStep~=0
@ -85,7 +85,7 @@ if varStep~=0
dt=min([maxStep,max([minStep,dt])]); dt=min([maxStep,max([minStep,dt])]);
% %
end end
% dt = 0.01;
tt=0; tt=0;
diffList=0; diffList=0;
@ -230,7 +230,6 @@ while tt(end)+dt<=maxTime+stepSizeEps
SysDataOrig=foldSysData(bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster); SysDataOrig=foldSysData(bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster);
[xTemp,diff,exitflag]=singleStepRK4(SimDatax,SysDataOrig,SysParaOrig,xt,tt(end),dt,iden,AEmethod); [xTemp,diff,exitflag]=singleStepRK4(SimDatax,SysDataOrig,SysParaOrig,xt,tt(end),dt,iden,AEmethod);
if exitflag~=0 if exitflag~=0
if varStep if varStep
dt=max([minStep,dt/10]); dt=max([minStep,dt/10]);
@ -257,7 +256,7 @@ while tt(end)+dt<=maxTime+stepSizeEps
maxErrX=max(abs(errX)); maxErrX=max(abs(errX));
errPtoX=[errPtoX;maxErrX/diffList(end)]; errPtoX=[errPtoX;maxErrX/diffList(end)];
stepRatio=min([max([(diffTol*mean(errPtoX)/2/maxErrX)^0.2,0.3]),2]); stepRatio=min([max([(diffTol*mean(errPtoX)/2/maxErrX)^0.2,0.3]),2]);
nxtDt=correctRatio^1.5*0.95*dt*stepRatio; nxtDt=correctRatio^1.9*0.99*dt*stepRatio;
correctRatio=sqrt(correctRatio); correctRatio=sqrt(correctRatio);
disp(''); disp('');
end end

4
internal/singleStepME.m
Unescape View File

@ -89,7 +89,9 @@ if AEmethod==0
exitflag=-1; exitflag=-1;
end end
else else
[stateNew,flag,diffRec,loop]=restorationAlgebraicNR(SimData,SysData,SysParax,xt,xTemp); %khuang Sep 6
[stateNew,flag,diffRec,loop]=solveAlgebraicNR(SimData,SysData,SysParax,xt,xTemp);
% [stateNew,flag,diffRec,loop]=restorationAlgebraicNR(SimData,SysData,SysParax,xt,xTemp);
if flag~=0 if flag~=0
addLog(['NR does not converge!'],'WARN'); addLog(['NR does not converge!'],'WARN');
exitflag=-1; exitflag=-1;

628
internal/singleStepTRAP.m
Unescape View File

@ -524,16 +524,24 @@ for iter=1:iterMax
end end
if ~isempty(syn) if ~isempty(syn)
if ~isempty(tg)
pActiveInSyn=ones(nSyn,1); pActiveInSyn=ones(nSyn,1);
tgovmActive=ones(nTg,1); tgovmActive=ones(nTg,1);
tgovmActive((tgovm0-Ttgmax)>0)=0; tgovmActive((tgovm0-Ttgmax)>0)=0;
tgovmActive((tgovm0-Ttgmin)<0)=0; tgovmActive((tgovm0-Ttgmin)<0)=0;
pActiveInSyn(tgIdx)=tgovmActive; pActiveInSyn(tgIdx)=tgovmActive;
% else
% tgIdx=[];
% tgovmActive=zeros(nTg,1);
% pActiveInSyn=zeros(nTg,1);
end
if ~isempty(exc)
texcActive=ones(nExc,1); texcActive=ones(nExc,1);
texcActive((Vavrf0-VavrMax)>0)=0; texcActive((Vavrf0-VavrMax)>0)=0;
texcActive((Vavrf0-VavrMin)<0)=0; texcActive((Vavrf0-VavrMin)<0)=0;
texcActiveInSyn=zeros(nSyn,1); texcActiveInSyn=zeros(nSyn,1);
texcActiveInSyn(tgIdx)=texcActive; texcActiveInSyn(tgIdx)=texcActive;
end
numSynOnBus=accumarray(syn(:,1),1,[nbus,1]); numSynOnBus=accumarray(syn(:,1),1,[nbus,1]);
%delta %delta
jffi=[jffi;idxs.deltaIdx]; jffi=[jffi;idxs.deltaIdx];
@ -633,4 +641,620 @@ dxc=(xTemp-xt)/dt;
SysParaxx=foldSysPara([],[],[],[],[],[],[],[],Ytr0+Ytr1*(t+dt),[],Ysh0+Ysh1*(t+dt),[],[VspSq2(:,1)+VspSq2(:,2)*(t+dt),zeros(size(VspSq2,1),1)],[],[],[],[]); SysParaxx=foldSysPara([],[],[],[],[],[],[],[],Ytr0+Ytr1*(t+dt),[],Ysh0+Ysh1*(t+dt),[],[VspSq2(:,1)+VspSq2(:,2)*(t+dt),zeros(size(VspSq2,1),1)],[],[],[],[]);
diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xTemp,dxc); diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xTemp,dxc);
end end
% % function [xTemp,diff,exitflag]=singleStepTRAP(SimData,SysData,SysPara,xt,t,dt,iden,AEmethod)
% % % Single step of trapezoidal formulation
% % % Copyright (C) Rui Yao <yaorui.thu@gmail.com>
% % %
% % % INPUT
% % % SimData - Simulation parameters
% % % SysData - System data for simulation
% % % SysPara - Parameters representing the events happening in the system
% % % xt - Initial system state
% % % t - Relative time compared with SysPara starting time
% % % dt - time step length
% % % iden - Aux identifier
% % % AEmethod - Method for solving algebraic equations
% % % 0 - HE
% % % 1 - NR
% % %
% % % OUTPUT
% % % xTemp - State at the new step
% % % diff - Error vector
% % % exitflag -
% % % 0 - Normally exit
% % % -1 - Fail
% % %
% %
% % [bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
% % % [nState,vIdx,qIdx,sIdx,deltaIdx,omegaIdx,eq1Idx,eq2Idx,ed1Idx,ed2Idx,psidIdx,psiqIdx,pgIdx,efIdx,...
% % % vavrmIdx,vavrrIdx,vavrfIdx,vavrrefIdx,tgovgIdx,tgovmIdx,tmechIdx,fIdx,dpgIdx,qpltIdx,vgIdx]...
% % % =getIndexDyn(SysData);
% % [nState,idxs]...
% % =getIndexDyn(SysData);
% % [pqIncr,pvIncr,Rind0,Rind1,Reind0,Reind1,Rzip0,Rzip1,Ytr0,Ytr1,Ysh0,Ysh1,VspSq2,~,~,~,~,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara);
% % [maxTime,segTime,~,nlvl,taylorN,alphaTol,diffTol,diffTolMax,method]=unfoldSimData(SimData);
% %
% % 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
% %
% % nbus=size(bus,1);
% %
% % 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~=0);
% % ipq=find(busType==0);
% % npq=size(ipq,1);
% % npv=size(ipv,1);
% %
% % % Determine the frequency model of each island
% % freqTypeTag=zeros(nIslands,1);%0:sw,1:syn,2:steady-state f
% % freqKeptTag=zeros(nbus,1);% corresponds to freqType=2
% % frefs=refs;
% % fswTag=zeros(nbus,1);
% % fsynTag=zeros(nbus,1);
% % fswTag(isw)=1;
% % fswTagxD=fswTag;
% % fsynTag(syn(:,1))=1;
% % D0=imag(xt(vIdx));
% % for isl=1:nIslands
% % if isempty(find(fswTag(islands==isl)==1, 1))
% % if isempty(find(fsynTag(islands==isl)==1, 1))
% % freqTypeTag(isl)=2;
% % busesInIsland=find(islands==isl);
% % [~,imin]=min(abs(D0(busesInIsland)));
% % frefs(isl)=busesInIsland(imin(1));
% % fswTagxD(frefs(isl))=1;
% % freqKeptTag(busesInIsland)=1;
% % else
% % freqTypeTag(isl)=1;
% % end
% % end
% % end
% % freqKeptTagxRef=freqKeptTag;
% % freqKeptTagxRef(frefs)=0;
% % nFreqKept=sum(freqKeptTag);
% %
% % if ~isempty(agc)
% % agcExt=zeros(nbus,size(agc,2));
% % agcExt(agc(:,1),:)=agc;
% % fdk=agcExt(:,2)+agcExt(:,3); %1/R+D
% % else
% % fdk=zeros(nbus,1);
% % end
% % fIngIdx=((2*nbus+1):(3*nbus))';
% %
% % stateFilterTag=ones(nState,1);
% % stateFilterTag(vIdx)=0;
% % stateFilterTag(qIdx)=0;
% % stateFilterTag(pgIdx)=0;
% % stateFilterTag(efIdx)=0;
% % stateFilterTag(tgovmIdx)=0;
% % stateFilterTag(fIdx)=0;
% % % stateFilterTag(dpgIdx)=0;
% % % busFilterTag=ones(nbus,1);
% % % busFilterTag(sw(:,1))=0;
% % % jacBusFilterTag=[busFilterTag;busFilterTag];
% %
% % busTag=ones(nbus,1);
% % busTag(sw(:,1))=0;
% % jggRowTag=[busTag;busTag;freqKeptTagxRef];
% % jggColTag=[busTag;fswTagxD==0;freqKeptTag];
% %
% % idxNonSw=find(busType~=2);
% % idxNonSwD=find(busType~=2&fswTagxD==1);
% %
% % nInd=size(ind,1);
% % if ~isempty(ind)
% % indIdx=ind(:,1);
% % Rm1=ind(:,7);
% % Xm1=ind(:,8);
% % Zm1=ind(:,7)+1j*ind(:,8);
% % Rm2=ind(:,9);
% % Xm2=ind(:,10);
% % Zmm=1j*ind(:,13);
% % Tm0=ind(:,15)+ind(:,16)+ind(:,17);
% % Tm1=-ind(:,16)-2*ind(:,17);
% % Tm2=ind(:,17);
% % Hm=ind(:,14);
% % end
% %
% % nSyn=size(syn,1);
% % if ~isempty(syn)
% % synIdx=syn(:,1);
% % wgb=syn(:,4);
% % model=syn(:,5);
% % Xgl=syn(:,6);
% % Rga=syn(:,7);
% % Xgd=syn(:,8);
% % Xgd1=syn(:,9);
% % Xgd2=syn(:,10);
% % Tgd1=syn(:,11);
% % Tgd2=syn(:,12);
% % Xgq=syn(:,13);
% % Xgq1=syn(:,14);
% % Xgq2=syn(:,15);
% % Tgq1=syn(:,16);
% % Tgq2=syn(:,17);
% % Mg=syn(:,18);
% % Dg=syn(:,19);
% % TgAA=syn(:,24);
% % gammad=Tgd2./Tgd1.*Xgd2./Xgd1.*(Xgd-Xgd1);
% % gammaq=Tgq2./Tgq1.*Xgq2./Xgq1.*(Xgq-Xgq1);
% % end
% %
% % nExc=size(exc,1);
% % if ~isempty(exc)
% % excIdx=exc(:,1);
% % VavrMax=exc(:,3);
% % VavrMin=exc(:,4);
% % muavr0=exc(:,5);
% % Tavr1=exc(:,7);
% % Tavr2=exc(:,6);
% % vavrf0=exc(:,8);
% % Vavr0=exc(:,9);
% % Tavre=exc(:,10);
% % Tavrr=exc(:,11);
% % end
% %
% % nTg=size(tg,1);
% % if ~isempty(tg)
% % tgIdx=tg(:,1);
% % wtgref=tg(:,3);
% % Rtg=tg(:,4);
% % Ttgmax=tg(:,5);
% % Ttgmin=tg(:,6);
% % Ttg2=tg(:,7);
% % Ttg1=tg(:,8);
% % end
% %
% % [V0,Q0,s0,d0,w0,eq10,eq20,ed10,ed20,psid0,psiq0,Pm0,Ef0,Vavrm0,Vavrr0,Vavrf0,Vavrref0,tgovg0,tgovm0,tgovmech0,f0,dpg0,qplt0,vg0]=unfoldX(xt,SysData);
% % if ~isempty(ind)
% % [YshInd0,Yshind1]=getLinearInterpolatorInd(nbus,ind,s0,s0);
% % Ysh0x=Ysh0+YshInd0;
% % Ysh1x=Ysh1+Yshind1;
% % else
% % Ysh0x=zeros(nbus,1);
% % Ysh1x=zeros(nbus,1);
% % end
% % Ytr=Ytr0+Ytr1*t;
% % Ysh=Ysh0x+Ysh1x*t;
% % YtrNew=Ytr0+Ytr1*(t+dt);
% % YshNew=Ysh0x+Ysh1x*(t+dt);
% %
% % SNew=-accumarray(pq(:,1),pq(:,4)+1j*pq(:,5),[nbus,1])+accumarray(pv(:,1),pv(:,4),[nbus,1]);
% % SNew=SNew-accumarray(pq(:,1),pqIncr(:,1)+1j*pqIncr(:,2),[nbus,1])*(t+dt)+accumarray(pv(:,1),pvIncr,[nbus,1])*(t+dt);
% % if size(pqIncr,2)>=4;SNew=SNew-accumarray(pq(:,1),pqIncr(:,3)+1j*pqIncr(:,4),[nbus,1])*(t+dt)*(t+dt);end
% % SNew=SNew-accumarray(zip(:,1),(Rzip0+Rzip1*(t+dt)).*(zip(:,7)+1j*zip(:,10)),[nbus,1]);
% %
% % INew=-accumarray(zip(:,1),(Rzip0+Rzip1*(t+dt)).*(zip(:,6)-1j*zip(:,9)),[nbus,1]);
% % JNew=real(INew);
% % KNew=imag(INew);
% %
% % YNew=YtrNew+sparse(1:nbus,1:nbus,YshNew,nbus,nbus);
% % GNew=real(YNew);BNew=imag(YNew);
% %
% % % Calculate diff and jac, do N-R
% % Vtemp=V0;
% % pqx=pq;pqx(:,[4,5])=pq(:,[4,5])+pqIncr(:,1:2)*t;if size(pqIncr,2)>=4;pqx(:,[4,5])=pqx(:,[4,5])+pqIncr(:,3:4)*t*t;end
% % pvx=pv;pvx(:,4)=pv(:,4)+pvIncr*t;
% % indx=ind;if ~isempty(ind);indx(:,15:17)=ind(:,15:17).*repmat(Rind0+Rind1*t,1,3);indx(:,13)=ind(:,13)./(Reind0+Reind1*t);end
% % zipx=zip;if ~isempty(zip);zipx(:,5:10)=zip(:,5:10).*repmat(Rzip0+Rzip1*t,1,6);end
% % SysDatax=foldSysData(bus,sw,pvx,pqx,shunt,line,indx,zipx,syn,exc,tg,agc,cac,cluster);
% %
% % pqxx=pq;pqxx(:,[4,5])=pq(:,[4,5])+pqIncr(:,1:2)*(t+dt);if size(pqIncr,2)>=4;pqxx(:,[4,5])=pqxx(:,[4,5])+pqIncr(:,3:4)*(t+dt)*(t+dt);end
% % pvxx=pv;pvxx(:,4)=pv(:,4)+pvIncr*(t+dt);
% % indxx=ind;if ~isempty(ind);indxx(:,15:17)=ind(:,15:17).*repmat(Rind0+Rind1*(t+dt),1,3);indxx(:,13)=ind(:,13)./(Reind0+Reind1*(t+dt));end
% % zipxx=zip;if ~isempty(zip);zipxx(:,5:10)=zip(:,5:10).*repmat(Rzip0+Rzip1*(t+dt),1,6);end
% % SysDataxx=foldSysData(bus,sw,pvxx,pqxx,shunt,line,indxx,zipxx,syn,exc,tg,agc,cac,cluster);
% %
% % SysData=foldSysData(bus,sw,pvx,pqx,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster);
% % dx1=integ(SysDatax,SysPara,xt,dt);
% % xTemp=xt;
% % xTemp=xTemp+dx1;
% %
% % iterMax=10;
% % exitflag=-1;
% % for iter=1:iterMax
% % % calculate df and dg
% % dx2=integ(SysDataxx,SysPara,xTemp,dt);
% % diffF=xt-xTemp+(dx1+dx2)/2;
% % diffF(stateFilterTag==0)=0;
% %
% % Ctemp=real(Vtemp);
% % Dtemp=imag(Vtemp);
% % Atemp=abs(Vtemp);
% % Vgtemp=Vtemp(syn(:,1));
% % ftemp=xTemp(fIdx);
% %
% % eq0=Ef0;ed0=zeros(size(Ef0));
% % [~,~,sNew,dNew,wNew,eq1New,eq2New,ed1New,ed2New,psidNew,psiqNew,EdNew,EfNew,VavrmNew,VavrrNew,VavrfNew,VavrrefNew,tgovgNew,tgovmNew,tgovmechNew,fNew,dpgNew,qpltNew,vgNew]=unfoldX(xTemp,SysData);
% % eqNew=EfNew;edNew=zeros(size(EfNew));
% %
% % [MatGV0,MatGV1,MatGRhs0,MatGRhs1]=getLinearInterpolatorSyn(syn,[],d0,dNew,ed0,ed10,ed20,edNew,ed1New,ed2New,eq0,eq10,eq20,eqNew,eq1New,eq2New,psid0,psiq0,psidNew,psiqNew);
% %
% % MatGVNew=MatGV0+MatGV1;
% % MatGRhsNew=MatGRhs0+MatGRhs1;
% % J1=GNew+sparse(syn(:,1),syn(:,1),MatGVNew(:,1),nbus,nbus);
% % J2=-BNew+sparse(syn(:,1),syn(:,1),MatGVNew(:,2),nbus,nbus);
% % J3=BNew+sparse(syn(:,1),syn(:,1),MatGVNew(:,3),nbus,nbus);
% % J4=GNew+sparse(syn(:,1),syn(:,1),MatGVNew(:,4),nbus,nbus);
% % P2freq=sparse(1:nbus,1:nbus,-freqKeptTag.*fdk,nbus,nbus);
% % Freq2freq=sparse([1:nbus,1:nbus],[1:nbus,frefs(islands)'],[ones(1,nbus),-ones(1,nbus)],nbus,nbus);
% %
% % Ig=accumarray(syn(:,1),MatGRhsNew(:,1)+1j*MatGRhsNew(:,2),[nbus,1])-...
% % accumarray(syn(:,1),(MatGVNew(:,1).*real(Vgtemp)+MatGVNew(:,2).*imag(Vgtemp))+1j*(MatGVNew(:,3).*real(Vgtemp)+MatGVNew(:,4).*imag(Vgtemp)),[nbus,1]);
% % % diffG=conj(SNew)+INew.*abs(Vtemp)+Ig.*conj(Vtemp)-(YNew*Vtemp).*conj(Vtemp)+freqKeptTag.*(dpgNew-fdk.*fNew);
% % % diffG(busFilterTag==0)=0;
% % % diffGf=zeros(nbus,1);
% %
% % diff=conj(SNew)+INew.*abs(Vtemp)+Ig.*conj(Vtemp)-(YNew*Vtemp).*conj(Vtemp);
% % diffP=real(diff)+freqKeptTag.*(-fdk.*ftemp+dpg0);
% % diffVQ=imag(diff);
% % diffGf=ftemp-ftemp(frefs(islands));
% % if ~isempty(pv)
% % diffVQ(pv(:,1))=abs(Vtemp(pv(:,1))).*abs(Vtemp(pv(:,1)))-(VspSq2(pv(:,1),1)+VspSq2(pv(:,1),2)*(t+dt));
% % end
% %
% % diffGx=[diffP;diffVQ;diffGf];
% % diffGx(jggRowTag==0)=0;
% %
% % if iter>1&&max(abs([diffF;diffGx]))<diffTol
% % exitflag=0;
% % break;
% % end
% %
% % % construct Jgg
% % dCtemp=sparse(1:nbus,1:nbus,Ctemp,nbus,nbus);
% % dDtemp=sparse(1:nbus,1:nbus,Dtemp,nbus,nbus);
% % Jac1=[dCtemp*J1+dDtemp*J3,dCtemp*J2+dDtemp*J4;...
% % dCtemp*J3-dDtemp*J1,dCtemp*J4-dDtemp*J2];
% % Jac2=[sparse(1:nbus,1:nbus,J1*Ctemp+J2*Dtemp,nbus,nbus),sparse(1:nbus,1:nbus,J3*Ctemp+J4*Dtemp,nbus,nbus);...
% % sparse(1:nbus,1:nbus,J3*Ctemp+J4*Dtemp,nbus,nbus),-sparse(1:nbus,1:nbus,J1*Ctemp+J2*Dtemp,nbus,nbus)];
% % Jac3=[sparse(1:nbus,1:nbus,JNew.*Ctemp./Atemp,nbus,nbus),sparse(1:nbus,1:nbus,JNew.*Dtemp./Atemp,nbus,nbus);...
% % sparse(1:nbus,1:nbus,KNew.*Ctemp./Atemp,nbus,nbus),sparse(1:nbus,1:nbus,KNew.*Dtemp./Atemp,nbus,nbus)];
% % Jac4=[sparse(syn(:,1),syn(:,1),MatGRhsNew(:,1),nbus,nbus),sparse(syn(:,1),syn(:,1),MatGRhsNew(:,2),nbus,nbus);...
% % sparse(syn(:,1),syn(:,1),MatGRhsNew(:,2),nbus,nbus),-sparse(syn(:,1),syn(:,1),MatGRhsNew(:,1),nbus,nbus)];
% % Jac=Jac3+Jac4-Jac1-Jac2;
% %
% % JacVc=sparse(1:nbus,1:nbus,2*Ctemp,nbus,nbus);
% % JacVd=sparse(1:nbus,1:nbus,2*Dtemp,nbus,nbus);
% %
% % Jac(nbus+ipv,:)=[JacVc(ipv,:),JacVd(ipv,:)];
% %
% % Jgg=[Jac,[P2freq;sparse(nbus,nbus)];...
% % sparse(nbus,2*nbus),Freq2freq];
% %
% % % construct Jgf
% %
% % jgfi=zeros(0,1);
% % jgfj=zeros(0,1);
% % jgfv=zeros(0,1);
% % jfgi=zeros(0,1);
% % jfgj=zeros(0,1);
% % jfgv=zeros(0,1);
% % jffi=zeros(0,1);
% % jffj=zeros(0,1);
% % jffv=zeros(0,1);
% %
% % if ~isempty(ind)
% % % s
% % deltasx=1e-3;
% % Zind0=Zm1+Zmm.*(Rm2+1j*Xm2.*sNew)./(Rm2+(Zmm+1j*Xm2).*sNew);
% % Yind0=1./Zind0;
% %
% % ZindNew=Zm1+Zmm.*(Rm2+1j*Xm2.*(sNew+deltasx))./(Rm2+(Zmm+1j*Xm2).*(sNew+deltasx));
% % YindNew=1./ZindNew;
% %
% % Yshind1x=YindNew-Yind0;
% % jgfi=[jgfi;indIdx;indIdx+nbus];
% % jgfj=[jgfj;sIdx;sIdx];
% % jgfv=[jgfv;Vtemp(indIdx).*conj(Vtemp(indIdx)).*real(Yshind1x/deltasx);Vtemp(indIdx).*conj(Vtemp(indIdx)).*imag(Yshind1x/deltasx)];
% % end
% %
% % if ~isempty(syn)
% % Vgtemp=Vtemp(synIdx);
% % Cgtemp=real(Vgtemp);
% % Dgtemp=imag(Vgtemp);
% %
% % cosd0=cos(dNew);
% % sind0=sin(dNew);
% % deltadNew=1e-5*abs(dNew)+1e-6;
% % cosdx=cos(dNew+deltadNew);
% % sindx=sin(dNew+deltadNew);
% %
% % MatGV0x=zeros(nSyn,4);
% % MatGV1x=zeros(nSyn,4);
% %
% % MatGRhs0x=zeros(nSyn,2);
% % MatGRhs1x=zeros(nSyn,2);
% %
% % for i=1:nSyn
% % Mg0=[sind0(i),-cosd0(i);cosd0(i),sind0(i)];
% % Mgx=[sindx(i),-cosdx(i);cosdx(i),sindx(i)];
% % if model(i)==6||model(i)==5
% % Yg=[Rga(i),Xgq2(i);-Xgd2(i),Rga(i)]/(Rga(i)*Rga(i)+Xgd2(i)*Xgq2(i));
% % Mxye=Mg0'*Yg;
% % Mdqe=Yg;
% % dIdq=-Yg*(Mgx-Mg0)*[Cgtemp(i);Dgtemp(i)]/deltadNew(i);
% % dVdq=-Yg*Mg0;
% % jgfi=[jgfi;synIdx(i);synIdx(i);synIdx(i)+nbus;synIdx(i)+nbus];
% % jgfj=[jgfj;ed2Idx(i);eq2Idx(i);ed2Idx(i);eq2Idx(i)];
% % jgfv=[jgfv;Cgtemp(i)*Mxye(1,1)+Dgtemp(i)*Mxye(2,1);Cgtemp(i)*Mxye(1,2)+Dgtemp(i)*Mxye(2,2);Cgtemp(i)*Mxye(2,1)-Dgtemp(i)*Mxye(1,1);Cgtemp(i)*Mxye(2,2)-Dgtemp(i)*Mxye(1,2)];
% %
% % jffi=[jffi;eq2Idx(i);eq2Idx(i);eq2Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);deltaIdx(i)];
% % jffv=[jffv;-(Xgd1(i)-Xgd2(i)+gammad(i))*Mdqe(1,1)/Tgd2(i);-(Xgd1(i)-Xgd2(i)+gammad(i))*Mdqe(1,2)/Tgd2(i);-(Xgd1(i)-Xgd2(i)+gammad(i))*dIdq(1)/Tgd2(i)];
% % jffi=[jffi;ed2Idx(i);ed2Idx(i);ed2Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);deltaIdx(i)];
% % jffv=[jffv;(Xgq1(i)-Xgq2(i)+gammaq(i))*Mdqe(2,1)/Tgq2(i);(Xgq1(i)-Xgq2(i)+gammaq(i))*Mdqe(2,2)/Tgq2(i);(Xgq1(i)-Xgq2(i)+gammaq(i))*dIdq(2)/Tgq2(i)];
% % jffi=[jffi;eq1Idx(i);eq1Idx(i);eq1Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);deltaIdx(i)];
% % jffv=[jffv;-(Xgd(i)-Xgd1(i)-gammad(i))*Mdqe(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i)-gammad(i))*Mdqe(1,2)/Tgd1(i);-(Xgd(i)-Xgd1(i)-gammad(i))*dIdq(1)/Tgd1(i)];
% % if model(i)==6
% % jffi=[jffi;ed1Idx(i);ed1Idx(i);ed1Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);deltaIdx(i)];
% % jffv=[jffv;(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,1)/Tgq1(i);(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,2)/Tgq1(i);(Xgq(i)-Xgq1(i)-gammaq(i))*dIdq(2)/Tgq1(i)];
% % end
% %
% % jfgi=[jfgi;eq2Idx(i);eq2Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;-(Xgd1(i)-Xgd2(i)+gammad(i))*dVdq(1,1)/Tgd2(i);-(Xgd1(i)-Xgd2(i)+gammad(i))*dVdq(1,2)/Tgd2(i)];
% % jfgi=[jfgi;ed2Idx(i);ed2Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;(Xgq1(i)-Xgq2(i)+gammaq(i))*dVdq(2,1)/Tgq2(i);(Xgq1(i)-Xgq2(i)+gammaq(i))*dVdq(2,2)/Tgq2(i)];
% % jfgi=[jfgi;eq1Idx(i);eq1Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;-(Xgd(i)-Xgd1(i)-gammad(i))*dVdq(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i)-gammad(i))*dVdq(1,2)/Tgd1(i)];
% % if model(i)==6
% % jfgi=[jfgi;ed1Idx(i);ed1Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,1)/Tgq1(i);(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,2)/Tgq1(i)];
% % end
% %
% % MatGRhs0x(i,:)=(Mg0'*Yg*[ed2New(i);eq2New(i)]).';
% % MatGRhs1x(i,:)=(Mgx'*Yg*[ed2New(i);eq2New(i)]).'-MatGRhs0x(i,:);
% % elseif model(i)==4
% % Yg=[Rga(i),Xgq1(i);-Xgd1(i),Rga(i)]/(Rga(i)*Rga(i)+Xgd1(i)*Xgq1(i));
% % Mxye=Mg0'*Yg;
% % Mdqe=Yg;
% % dVdq=-Yg*Mg0;
% % dIdq=-Yg*(Mgx-Mg0)*[Cgtemp(i);Dgtemp(i)]/deltadNew(i);
% % jgfi=[jgfi;synIdx(i);synIdx(i);synIdx(i)+nbus;synIdx(i)+nbus];
% % jgfj=[jgfj;ed1Idx(i);eq1Idx(i);ed1Idx(i);eq1Idx(i)];
% % jgfv=[jgfv;Cgtemp(i)*Mxye(1,1)+Dgtemp(i)*Mxye(2,1);Cgtemp(i)*Mxye(1,2)+Dgtemp(i)*Mxye(2,2);Cgtemp(i)*Mxye(2,1)-Dgtemp(i)*Mxye(1,1);Cgtemp(i)*Mxye(2,2)-Dgtemp(i)*Mxye(1,2)];
% %
% % jffi=[jffi;eq1Idx(i);eq1Idx(i);eq1Idx(i)];
% % jffj=[jffj;ed1Idx(i);eq1Idx(i);deltaIdx(i)];
% % jffv=[jffv;-(Xgd(i)-Xgd1(i))*Mdqe(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i))*Mdqe(1,2)/Tgd1(i);-(Xgd(i)-Xgd1(i))*dIdq(1)/Tgd1(i)];
% % jffi=[jffi;ed1Idx(i);ed1Idx(i);ed1Idx(i)];
% % jffj=[jffj;ed1Idx(i);eq1Idx(i);deltaIdx(i)];
% % jffv=[jffv;(Xgq(i)-Xgq1(i))*Mdqe(2,1)/Tgq1(i);(Xgq(i)-Xgq1(i))*Mdqe(2,2)/Tgq1(i);(Xgq(i)-Xgq1(i))*dIdq(2)/Tgq1(i)];
% %
% % jfgi=[jfgi;eq1Idx(i);eq1Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;-(Xgd(i)-Xgd1(i))*dVdq(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i))*dVdq(1,2)/Tgd1(i)];
% % jfgi=[jfgi;ed1Idx(i);ed1Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;(Xgq(i)-Xgq1(i))*Mdqe(2,1)/Tgq1(i);(Xgq(i)-Xgq1(i))*Mdqe(2,2)/Tgq1(i)];
% %
% % MatGRhs0x(i,:)=(Mg0'*Yg*[ed1New(i);eq1New(i)]).';
% % MatGRhs1x(i,:)=(Mgx'*Yg*[ed1New(i);eq1New(i)]).'-MatGRhs0x(i,:);
% % elseif model(i)==3||model(i)==2
% % Yg=[Rga(i),Xgq(i);-Xgd1(i),Rga(i)]/(Rga(i)*Rga(i)+Xgd1(i)*Xgq(i));
% % Mxye=Mg0'*Yg;
% % Mdqe=Yg;
% % dVdq=-Yg*Mg0;
% % dIdq=-Yg*(Mgx-Mg0)*[Cgtemp(i);Dgtemp(i)]/deltadNew(i);
% % if model(i)==3
% % jgfi=[jgfi;synIdx(i);synIdx(i)+nbus];
% % jgfj=[jgfj;eq1Idx(i);eq1Idx(i)];
% % jgfv=[jgfv;Cgtemp(i)*Mxye(1,2)+Dgtemp(i)*Mxye(2,2);Cgtemp(i)*Mxye(2,2)-Dgtemp(i)*Mxye(1,2)];
% % jffi=[jffi;eq1Idx(i);eq1Idx(i)];
% % jffj=[jffj;eq1Idx(i);deltaIdx(i)];
% % jffv=[jffv;-(Xgd(i)-Xgd1(i))*Mdqe(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i))*dIdq(1)/Tgd1(i)];
% % jfgi=[jfgi;eq1Idx(i);eq1Idx(i)];
% % jfgj=[jfgj;synIdx(i);synIdx(i)+nbus];
% % jfgv=[jfgv;-(Xgd(i)-Xgd1(i))*dVdq(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i))*dVdq(1,2)/Tgd1(i)];
% % end
% %
% % MatGRhs0x(i,:)=(Mg0'*Yg*[edNew(i);eq1New(i)]).';
% % MatGRhs1x(i,:)=(Mgx'*Yg*[edNew(i);eq1New(i)]).'-MatGRhs0x(i,:);
% % elseif model(i)==8
% % Yg=zeros(2);
% % Mxye=Mg0'*diag([1/Xgd2(i);1/Xgq2(i)]);
% % Mdqe=diag([1/Xgd2(i);1/Xgq2(i)]);
% % dVdq=-Yg*Mg0;
% % dIdq=-Yg*(Mgx-Mg0)*[Cgtemp(i);Dgtemp(i)]/deltadNew(i);
% % jgfi=[jgfi;synIdx(i);synIdx(i);synIdx(i)+nbus;synIdx(i)+nbus];
% % jgfj=[jgfj;ed2Idx(i);eq2Idx(i);ed2Idx(i);eq2Idx(i)];
% % jgfv=[jgfv;-(Cgtemp(i)*Mxye(1,1)+Dgtemp(i)*Mxye(2,1));Cgtemp(i)*Mxye(1,2)+Dgtemp(i)*Mxye(2,2);-(Cgtemp(i)*Mxye(2,1)-Dgtemp(i)*Mxye(1,1));Cgtemp(i)*Mxye(2,2)-Dgtemp(i)*Mxye(1,2)];
% % jgfi=[jgfi;synIdx(i);synIdx(i);synIdx(i)+nbus;synIdx(i)+nbus];
% % jgfj=[jgfj;psiqIdx(i);psidIdx(i);psiqIdx(i);psidIdx(i)];
% % jgfv=[jgfv;-(Cgtemp(i)*Mxye(1,1)+Dgtemp(i)*Mxye(2,1));-(Cgtemp(i)*Mxye(1,2)+Dgtemp(i)*Mxye(2,2));-(Cgtemp(i)*Mxye(2,1)-Dgtemp(i)*Mxye(1,1));-(Cgtemp(i)*Mxye(2,2)-Dgtemp(i)*Mxye(1,2))];
% %
% % jffi=[jffi;eq2Idx(i);eq2Idx(i);eq2Idx(i);eq2Idx(i);eq2Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);psiqIdx(i);psidIdx(i);deltaIdx(i)];
% % jffv=[jffv;-(Xgd1(i)-Xgd2(i)+gammad(i))*Mdqe(1,1)/Tgd2(i);-(Xgd1(i)-Xgd2(i)+gammad(i))*Mdqe(1,2)/Tgd2(i);...
% % -(Xgd1(i)-Xgd2(i)+gammad(i))*Mdqe(1,1)/Tgd2(i);(Xgd1(i)-Xgd2(i)+gammad(i))*Mdqe(1,2)/Tgd2(i);-(Xgd1(i)-Xgd2(i)+gammad(i))*dIdq(1)/Tgd2(i)];
% % jffi=[jffi;ed2Idx(i);ed2Idx(i);ed2Idx(i);ed2Idx(i);ed2Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);psiqIdx(i);psidIdx(i);deltaIdx(i)];
% % jffv=[jffv;(Xgq1(i)-Xgq2(i)+gammaq(i))*Mdqe(2,1)/Tgq2(i);(Xgq1(i)-Xgq2(i)+gammaq(i))*Mdqe(2,2)/Tgq2(i);...
% % (Xgq1(i)-Xgq2(i)+gammaq(i))*Mdqe(2,1)/Tgq2(i);-(Xgq1(i)-Xgq2(i)+gammaq(i))*Mdqe(2,2)/Tgq2(i);(Xgq1(i)-Xgq2(i)+gammaq(i))*dIdq(2)/Tgq2(i)];
% % jffi=[jffi;eq1Idx(i);eq1Idx(i);eq1Idx(i);eq1Idx(i);eq1Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);psiqIdx(i);psidIdx(i);deltaIdx(i)];
% % jffv=[jffv;-(Xgd(i)-Xgd1(i)-gammad(i))*Mdqe(1,1)/Tgd1(i);-(Xgd(i)-Xgd1(i)-gammad(i))*Mdqe(1,2)/Tgd1(i);...
% % -(Xgd(i)-Xgd1(i)-gammad(i))*Mdqe(1,1)/Tgd1(i);(Xgd(i)-Xgd1(i)-gammad(i))*Mdqe(1,2)/Tgd1(i);-(Xgd(i)-Xgd1(i)-gammad(i))*dIdq(1)/Tgd1(i)];
% % jffi=[jffi;ed1Idx(i);ed1Idx(i);ed1Idx(i);ed1Idx(i);ed1Idx(i)];
% % jffj=[jffj;ed2Idx(i);eq2Idx(i);psiqIdx(i);psidIdx(i);deltaIdx(i)];
% % jffv=[jffv;(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,1)/Tgq1(i);(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,2)/Tgq1(i);...
% % (Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,1)/Tgq1(i);-(Xgq(i)-Xgq1(i)-gammaq(i))*Mdqe(2,2)/Tgq1(i);(Xgq(i)-Xgq1(i)-gammaq(i))*dIdq(2)/Tgq1(i)];
% %
% % MatGRhs0x(i,:)=(Mg0'*[(eq2New(i)-psidNew(i))/Xgd2(i);(-ed2New(i)-psiqNew(i))/Xgq2(i)]).';
% % MatGRhs1x(i,:)=(Mgx'*[(eq2New(i)-psidNew(i))/Xgd2(i);(-ed2New(i)-psiqNew(i))/Xgq2(i)]).'-MatGRhs0x(i,:);
% % end
% % MatVg0x=Mg0'*Yg*Mg0;
% % MatVgNewx=Mgx'*Yg*Mgx;
% % MatGV0x(i,:)=[MatVg0x(1,1),MatVg0x(1,2),MatVg0x(2,1),MatVg0x(2,2)];
% % MatGV1x(i,:)=[MatVgNewx(1,1),MatVgNewx(1,2),MatVgNewx(2,1),MatVgNewx(2,2)]-MatGV0x(i,:);
% %
% % jgfi=[jgfi;synIdx(i);synIdx(i)+nbus];
% % jgfj=[jgfj;deltaIdx(i);deltaIdx(i)];
% % jgfv=[jgfv;...
% % (Cgtemp(i)*(MatGRhs1x(i,1)-MatGV1x(i,1)*Cgtemp(i)-MatGV1x(i,2)*Dgtemp(i))+Dgtemp(i)*(MatGRhs1x(i,2)-MatGV1x(i,3)*Cgtemp(i)-MatGV1x(i,4)*Dgtemp(i)))/deltadNew(i);...
% % (-Dgtemp(i)*(MatGRhs1x(i,1)-MatGV1x(i,1)*Cgtemp(i)-MatGV1x(i,2)*Dgtemp(i))+Cgtemp(i)*(MatGRhs1x(i,2)-MatGV1x(i,3)*Cgtemp(i)-MatGV1x(i,4)*Dgtemp(i)))/deltadNew(i)];
% % end
% %
% % end
% % Jgf=sparse(jgfi,jgfj,jgfv,3*nbus,nState);
% %
% % % construct Jff & Jfg
% %
% % if ~isempty(ind)
% % % s
% % jffi=[jffi;sIdx];
% % jffj=[jffj;sIdx];
% % divs=((Rm2+Rm1.*sNew).*(Rm2+Rm1.*sNew)+(Xm2+Xm1).*(Xm2+Xm1).*sNew.*sNew);
% % jffv=[jffv;(Tm1+2*Tm2.*sNew+Vtemp(indIdx).*conj(Vtemp(indIdx)).*Rm2.*(1./divs-2*sNew.*(sNew.*(Xm2+Xm1).*(Xm2+Xm1)+Rm1.*(Rm2+Rm1.*sNew))./divs./divs))/2./Hm];
% % % il=V(indIdx).*yind;
% % % ir=il-(V(indIdx)-il.*Zm1)./Zmm;
% % % dx(sIdx)=dt*((Tm0+Tm1.*s+Tm2.*s.*s)-ir.*conj(ir).*Rm2./s)/2./Hm;
% %
% % yind=1./(Zm1+Zmm.*(Rm2+1j*Xm2.*sNew)./(Rm2+sNew.*(Zmm+1j*Xm2)));
% % yInds=(1+Zm1./Zmm).*yind-1./Zmm;
% % yIndsRs=yInds.*conj(yInds).*Rm2./sNew;
% % yIndsRs(isnan(yIndsRs))=0;
% % jfgi=[jfgi;sIdx;sIdx];
% % jfgj=[jfgj;indIdx;indIdx+nbus];
% % jfgv=[jfgv;2*real(Vtemp(indIdx)).*yIndsRs/2./Hm;2*imag(Vtemp(indIdx)).*yIndsRs/2./Hm];
% % end
% %
% % if ~isempty(syn)
% % pActiveInSyn=ones(nSyn,1);
% % tgovmActive=ones(nTg,1);
% % tgovmActive((tgovm0-Ttgmax)>0)=0;
% % tgovmActive((tgovm0-Ttgmin)<0)=0;
% % pActiveInSyn(tgIdx)=tgovmActive;
% % texcActive=ones(nExc,1);
% % texcActive((Vavrf0-VavrMax)>0)=0;
% % texcActive((Vavrf0-VavrMin)<0)=0;
% % texcActiveInSyn=zeros(nSyn,1);
% % texcActiveInSyn(tgIdx)=texcActive;
% % numSynOnBus=accumarray(syn(:,1),1,[nbus,1]);
% % %delta
% % jffi=[jffi;deltaIdx];
% % jffj=[jffj;omegaIdx];
% % jffv=[jffv;wgb];
% % %omega
% % jffi=[jffi;omegaIdx;omegaIdx(tgIdx);omegaIdx(tgIdx);omegaIdx(tgIdx);omegaIdx];
% % jffj=[jffj;omegaIdx;tgovgIdx;omegaIdx(tgIdx);tmechIdx;dpgIdx(synIdx)];
% % jffv=[jffv;-Dg/2./Mg;tgovmActive;-tgovmActive.*Ttg1./Ttg2./Rtg;tgovmActive;pActiveInSyn./numSynOnBus(synIdx)];
% % %psiq
% % jffi=[jffi;psidIdx(model>=8)];
% % jffj=[jffj;psiqIdx(model>=8)];
% % jffv=[jffv;wgb(model>=8)];
% % %psid
% % jffi=[jffi;psiqIdx(model>=8)];
% % jffj=[jffj;psidIdx(model>=8)];
% % jffv=[jffv;-wgb(model>=8)];
% % %ed''
% % jffi=[jffi;ed2Idx(model>=5);ed2Idx(model>=5)];
% % jffj=[jffj;ed2Idx(model>=5);ed1Idx(model>=5)];
% % jffv=[jffv;-1./Tgq2(model>=5);1./Tgq2(model>=5)];
% % %eq''
% % jffi=[jffi;eq2Idx(model>=5);eq2Idx(model>=5);eq2Idx(model>=5&texcActiveInSyn==1)];
% % jffj=[jffj;eq2Idx(model>=5);eq1Idx(model>=5);vavrfIdx(texcActive==1&model(excIdx)>=5)];
% % jffv=[jffv;-1./Tgd2(model>=5);1./Tgd2(model>=5);TgAA(model>=5&texcActiveInSyn==1)./Tgd1(model>=5&texcActiveInSyn==1)./Tgd2(model>=5&texcActiveInSyn==1)];
% % %ed'
% % jffi=[jffi;ed1Idx(model>=4)];
% % jffj=[jffj;ed1Idx(model>=4)];
% % jffv=[jffv;-1./Tgq1(model>=4)];
% % %eq'
% % jffi=[jffi;eq1Idx(model>=3)];
% % jffj=[jffj;eq1Idx(model>=3)];
% % jffv=[jffv;-1./Tgd1(model>=3)];
% % jffi=[jffi;eq1Idx(model>=5&texcActiveInSyn==1)];
% % jffj=[jffj;vavrfIdx(texcActive==1&model(excIdx)>=5)];
% % jffv=[jffv;(1-TgAA(model>=5&texcActiveInSyn==1)./Tgd1(model>=5&texcActiveInSyn==1))./Tgd2(model>=5&texcActiveInSyn==1)];
% % jffi=[jffi;eq1Idx((model==3|model==4)&texcActiveInSyn==1)];
% % jffj=[jffj;vavrfIdx(texcActive==1&(model(excIdx)==3|model(excIdx)==4))];
% % jffv=[jffv;(1-TgAA((model==3|model==4)&texcActiveInSyn==1)./Tgd1((model==3|model==4)&texcActiveInSyn==1))./Tgd2((model==3|model==4)&texcActiveInSyn==1)];
% % end
% %
% % if ~isempty(tg)
% % jffi=[jffi;tgovgIdx;tgovgIdx];
% % jffj=[jffj;omegaIdx;tgovgIdx];
% % jffv=[jffv;-(1-Ttg1./Ttg2)./Rtg./Ttg2;-1./Ttg2];
% % end
% %
% % if ~isempty(exc)
% % VmagAvr=abs(Vtemp(synIdx(excIdx)));
% %
% % jffi=[jffi;vavrmIdx];
% % jffj=[jffj;vavrmIdx];
% % jffv=[jffv;-1./Tavrr];
% %
% % jffi=[jffi;vavrrIdx;vavrrIdx;vavrrIdx];
% % jffj=[jffj;vavrmIdx;vavrrefIdx;vavrrIdx];
% % jffv=[jffv;-muavr0.*(1-Tavr1./Tavr2)./Tavr2;muavr0.*(1-Tavr1./Tavr2)./Tavr2;-1./Tavr2];
% %
% % jffi=[jffi;vavrfIdx;vavrfIdx;vavrfIdx;vavrfIdx];
% % jffj=[jffj;vavrmIdx;vavrrefIdx;vavrrIdx;vavrfIdx];
% % jffv=[jffv;-muavr0.*Tavr1./Tavr2./Tavre.*VmagAvr./Vavr0;muavr0.*Tavr1./Tavr2./Tavre.*VmagAvr./Vavr0;VmagAvr./Vavr0./Tavre;-1./Tavre];
% % end
% %
% % if ~isempty(agc)
% % agcExt=zeros(nbus,size(agc,2));
% % agcExt(agc(:,1),:)=agc;
% %
% % jfgi=[jfgi;dpgIdx(freqKeptTag==1)];
% % jfgj=[jfgj;fIngIdx(freqKeptTag==1)];%f
% % jfgv=[jfgv;-agcExt(freqKeptTag==1,4)];
% %
% % numSynOnBus=accumarray(syn(:,1),1,[nbus,1]);
% % jffi=[jffi;dpgIdx(syn(:,1))];
% % jffj=[jffj;omegaIdx];%f
% % jffv=[jffv;-agcExt(syn(:,1),4)./numSynOnBus(syn(:,1))];
% % end
% %
% % Jff=dt*sparse(jffi,jffj,jffv,nState,nState)/2-speye(nState,nState);
% %
% % % construct Jfg
% % Jfg=dt*sparse(jfgi,jfgj,jfgv,nState,3*nbus)/2; %P,Q,f
% %
% % jac=[Jff(stateFilterTag==1,stateFilterTag==1),Jfg(stateFilterTag==1,jggColTag==1);Jgf(jggRowTag==1,stateFilterTag==1),Jgg(jggRowTag==1,jggColTag==1)];
% %
% % correction=-jac\[diffF(stateFilterTag==1);diffGx(jggRowTag==1)];
% % nStateCorr=sum(stateFilterTag);
% % xTemp(stateFilterTag==1)=xTemp(stateFilterTag==1)+correction(1:nStateCorr);
% % corrFull=zeros(3*nbus,1);
% % corrFull(jggRowTag==1)=correction((nStateCorr+1):end);
% % Vtemp=Vtemp+corrFull(1:nbus)+1j*corrFull((nbus+1):(2*nbus));
% % xTemp(vIdx)=Vtemp;
% % xTemp(fIdx)=corrFull(fIngIdx);
% % xTemp=adjustAlgebraic(SysDataxx,xTemp);
% % end
% %
% % dxc=(xTemp-xt)/dt;
% % SysParaxx=foldSysPara([],[],[],[],[],[],[],[],Ytr0+Ytr1*(t+dt),[],Ysh0+Ysh1*(t+dt),[],[VspSq2(:,1)+VspSq2(:,2)*(t+dt),zeros(size(VspSq2,1),1)],[],[],[],[]);
% % diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xTemp,dxc);
% %
% % end

7
runDynamicSimulationExec.m
Unescape View File

@ -163,7 +163,7 @@ else
end end
%% Temporary code here %% Temporary code here
% eventList(:,6)=0.0; % eventList(:,6)=4.00;
% eventList(:,7)=0.1; % eventList(:,7)=0.1;
% Efstd=1.2; % Efstd=1.2;
% evtDynZip(:,2)=0.2; % evtDynZip(:,2)=0.2;
@ -194,8 +194,9 @@ if ~hotStart
[SysData,x0,finalAlpha,alphaList,diff]=calculateInitialState(SysData,SimDataStatic,SysPara); [SysData,x0,finalAlpha,alphaList,diff]=calculateInitialState(SysData,SimDataStatic,SysPara);
else % Full system steady-state else % Full system steady-state
[busBase,pvBase,pqBase,swBase,lineBase,shuntBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase,pm]=... % [busBase,pvBase,pqBase,swBase,lineBase,shuntBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase,pm]=...
regulateSystemData(busBase,pvBase,pqBase,swBase,lineBase,shuntBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase); % regulateSystemData(busBase,pvBase,pqBase,swBase,lineBase,shuntBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase);
% pvBase(:,4)=pm;
bus=busBase;sw=swBase;pv=pvBase;pq=pqBase;shunt=shuntBase; bus=busBase;sw=swBase;pv=pvBase;pq=pqBase;shunt=shuntBase;
line=lineBase;ind=indBase;zip=zipBase;syn=synBase;exc=excBase;tg=tgBase; line=lineBase;ind=indBase;zip=zipBase;syn=synBase;exc=excBase;tg=tgBase;
agc=agcBase;cac=cacBase;cluster=clusterBase; agc=agcBase;cac=cacBase;cluster=clusterBase;

5
runPowerSAS.m
Unescape View File

@ -60,6 +60,7 @@ else
addLog(['OUTPUT IS DISABLED. Check ',options.outputPath,'/~default_log.txt for logs.'],'DEBUG'); addLog(['OUTPUT IS DISABLED. Check ',options.outputPath,'/~default_log.txt for logs.'],'DEBUG');
end end
addLog(['The session ID is ',timestamp,'.'],'INFO'); addLog(['The session ID is ',timestamp,'.'],'INFO');
%khuang 919
options.timestamp=timestamp; options.timestamp=timestamp;
homePath=pwd; homePath=pwd;
@ -167,9 +168,9 @@ function options=regulateOptions(options)
if nargin<1 if nargin<1
options=[]; options=[];
end end
if ~isfield(options,'dataPath');options.dataPath=[pwd,'\data'];end if ~isfield(options,'dataPath');options.dataPath=[pwd,'/data'];end
if ~isfield(options,'settingPath');options.settingPath=options.dataPath;end if ~isfield(options,'settingPath');options.settingPath=options.dataPath;end
if ~isfield(options,'dbstop');options.dbstop=0;end if ~isfield(options,'dbstop');options.dbstop=0;end
if ~isfield(options,'output');options.output=0;end if ~isfield(options,'output');options.output=0;end
if ~isfield(options,'outputPath');options.outputPath=[pwd,'\output'];end if ~isfield(options,'outputPath');options.outputPath=[pwd,'/output'];end
end end

Write Preview
Loading…
Cancel
Save