ANL-CEEESA
/
powersas.m
mirror of https://github.com/ANL-CEEESA/powersas.m
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'release'
${ noResults }
powersas.m/runDynamicSimulationExec.m

1305 lines
56 KiB
Raw Permalink Blame History

function res=runDynamicSimulationExec(caseName,SysDataBase,simSettings,options,snapshot)
% Main function for running dynamic simulation
% Author: Rui Yao <ruiyao@ieee.org>
%
% Copyright (C) 2021, UChicago Argonne, LLC. All rights reserved.
%
% OPEN SOURCE LICENSE
%
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
%
%
% ******************************************************************************************************
% DISCLAIMER
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
% OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% ***************************************************************************************************
%
% FUNCTION runDynamicSimulationExec
%
% INPUT
% caseName - Name of the studied case
% SysDataBase - Studied system data struct
% simSettings - Settings
% options - Other options
% snapshot - snap shot containing initial state
%
% OUTPUT
% res struct:
% flag -
% msg - Message returned
% t - time sequence
% stateCurve - list of system states corresponding to t
% timestamp - the timestamp of the analysis
% SysDataBase - Studied system data struct
% eventList - list of system events
% snapshot - the snapshot at the end of the analysis
% caseName - Name of the studied case
% simSettings - Settings
%
flag=0;
res=[];
if nargin<4
options=regulateOptions();
else
options=regulateOptions(options);
end
if options.dbstop
dbstop if error
end
if isfield(options,'timestamp')
timestamp=options.timestamp;
else
timestamp=datestr(clock,'yyyymmddTHHMMSSFFF');
end
addLog('================= START ==================','INFO');
% addLog('--------- WARNING ---------','WARN');
% addLog(disclaimerStr,'WARN');
% restSettingFile='restoration_plan_test_039';
comment=options.comment;
if isfield(simSettings,'nlvl');nlvl=simSettings.nlvl;else nlvl=options.nlvl;end
if isfield(simSettings,'segAlpha');segAlpha=simSettings.segAlpha;else segAlpha=options.segAlpha;end
if isfield(simSettings,'dAlpha');dAlpha=simSettings.dAlpha;else dAlpha=options.dAlpha;end
if isfield(simSettings,'taylorN');taylorN=simSettings.taylorN;else taylorN=options.taylorN;end
if isfield(simSettings,'alphaTol');alphaTol=simSettings.alphaTol;else alphaTol=options.alphaTol;end
if isfield(simSettings,'diffTol');diffTol=simSettings.diffTol;else diffTol=options.diffTol;end
if isfield(simSettings,'diffTolMax');diffTolMax=simSettings.diffTolMax;else diffTolMax=options.diffTolMax;end
if isfield(simSettings,'method');method=simSettings.method;else method=options.method;end
if isfield(simSettings,'diffTolCtrl');diffTolCtrl=simSettings.diffTolCtrl;else diffTolCtrl=options.diffTolCtrl;end
eventList=simSettings.eventList;
bsSyn=simSettings.bsSyn;
bsBus=simSettings.bsBus;
bsInd=simSettings.bsInd;
Efstd=simSettings.Efstd;
evtLine=simSettings.evtLine;
evtLineSpec=simSettings.evtLineSpec;
evtZip=simSettings.evtZip;
evtZipSpec=simSettings.evtZipSpec;
evtZipSpec2=simSettings.evtZipSpec2;
evtInd=simSettings.evtInd;
evtIndSpec=simSettings.evtIndSpec;
evtSyn=simSettings.evtSyn;
evtSynSpec=simSettings.evtSynSpec;
evtFault=simSettings.evtFault;
evtFaultSpec=simSettings.evtFaultSpec;
evtDyn=simSettings.evtDyn;
evtDynPQ=simSettings.evtDynPQ;
evtDynPV=simSettings.evtDynPV;
evtDynInd=simSettings.evtDynInd;
evtDynZip=simSettings.evtDynZip;
evtDynSh=simSettings.evtDynSh;
evtDynZipRamp=simSettings.evtDynZipRamp;
evtDynTmech=simSettings.evtDynTmech;
evtDynPm=simSettings.evtDynPm;
evtDynEf=simSettings.evtDynEf;
evtDynVref=simSettings.evtDynVref;
evtDynEq1=simSettings.evtDynEq1;
if exist('envDependentScript','file')
envDependentScript;
end
[eventList,evtDyn,evtDynInd,evtDynZip]=formatEventList(eventList,evtDyn,evtDynInd,evtDynZip);
EvtType=getSimEventType();
MethodType=getSimMethodType();
SimSet=getSimSettingDefault();
DynSimFlag=getDynSimFlags();
if ~options.inheritParaSettings
nlvl = options.nlvl;
taylorN = options.taylorN;
alphaTol =options.alphaTol;
diffTol =options.diffTol;
diffTolMax =options.diffTolMax;
method = options.method;
diffTolCtrl =options.diffTolCtrl;
end
tTol=options.tTol;
diffTolOrig=diffTol;
allowReverse=options.allowReverse;
allowSteadyDynSwitch=options.allowSteadyDynSwitch;
useDiffCtrl=options.useDiffCtrl;
useNoMoveCtrl=options.useNoMoveCtrl;
nPool=options.nPool;
switchMaskCondition='event(2)>10 && ~(eventList(item+1,4)==EvtType.DYN_SIM&&eventList(item+1,5)~=0&&eventList(item+1,3)-eventList(item+1,2)<10)';
if useDiffCtrl&&(~exist('diffTolCtrl','var')||isempty(diffTolCtrl));useDiffCtrl=0;end
usePar=license('test','Distrib_Computing_Toolbox');
if usePar;delete(gcp('nocreate'));end
if usePar&&exist('nPool','var')&&nPool>1
parpool(nPool);
pctRunOnAll warning('off','MATLAB:singularMatrix');
pctRunOnAll warning('off','MATLAB:nearlySingularMatrix');
pctRunOnAll warning('off','MATLAB:illConditionedMatrix');
else
warning('off','MATLAB:singularMatrix');
warning('off','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:illConditionedMatrix');
end
hotStart=options.hotStart;
if hotStart && nargin>=5 && isstruct(snapshot) && isfield(snapshot,'SysData') && isfield(snapshot,'state')
else
hotStart=0;
end
%% Temporary code here
% eventList(:,6)=4.00;
% eventList(:,7)=0.1;
% Efstd=1.2;
% evtDynZip(:,2)=0.2;
%%
startTag=tic;
initSuccess=1;
if ~hotStart
[busBase,swBase,pvBase,pqBase,shuntBase,lineBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase]=unfoldSysData(SysDataBase);
if eventList(1,5)==0 % Follow black start variables
% TODO: AGC, CAC and Cluster
nbusBase=size(busBase,1);
synTag=zeros(nbusBase,1);
if ~isempty(synBase);synTag(synBase(:,1))=1;end
swBase=swBase(synTag(swBase(:,1))==0,:);
if ~isempty(pvBase);pvBase=pvBase(synTag(pvBase(:,1))==0,:);end
[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);
SysDataBase=foldSysData(busBase,swBase,pvBase,pqBase,shuntBase,lineBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase);
[SysData,maps,Efstd,P0,pShare,Tm]=...
generateBlackStartNet(SysDataBase,bsBus,bsSyn,bsInd);
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
SimDataStatic=foldSimData([],segAlpha,dAlpha,nlvl,taylorN,alphaTol,diffTol,diffTolMax,[]);
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],(Efstd-1).*ones(size(syn,1),1),[],[]);
[SysData,x0,finalAlpha,alphaList,diff]=calculateInitialState(SysData,SimDataStatic,SysPara);
else % Full system steady-state
% [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);
% pvBase(:,4)=pm;
bus=busBase;sw=swBase;pv=pvBase;pq=pqBase;shunt=shuntBase;
line=lineBase;ind=indBase;zip=zipBase;syn=synBase;exc=excBase;tg=tgBase;
agc=agcBase;cac=cacBase;cluster=clusterBase;
maps.busMap=(1:size(busBase,1))';maps.invBusMap=maps.busMap;
maps.indMap=(1:size(indBase,1))';maps.invIndMap=maps.indMap;
maps.synMap=(1:size(synBase,1))';maps.invSynMap=maps.synMap;
maps.lineMap=(1:size(lineBase,1))';maps.invLineMap=maps.lineMap;
SimDataStatic=foldSimData([],segAlpha,[],nlvl,taylorN,alphaTol,diffTol,diffTolMax,[]);
SysData=foldSysData(bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster);
SysDataBase=foldSysData(busBase,swBase,pvBase,pqBase,shuntBase,lineBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase);
if size(Efstd,1)==1;Efstd=Efstd*ones(size(syn,1),1);end
if isempty(Efstd);Ef1=[];else Ef1=(Efstd-1).*ones(size(syn,1),1);end
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],Ef1,[],[]);
[SysData,x0,finalAlpha,alphaList,diff]=calculateInitialState(SysData,SimDataStatic,SysPara);
end
if finalAlpha<1
initSuccess=0;
msg='Fail to get initial steady-state point. Abort.';
flag=-1;
end
fault=zeros(0,4);
snapshot.SysData=SysData;
snapshot.state=x0;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst={snapshot};
SysDataList={SysData};
stateList={x0};
tList={0};
mapsList={maps};
faultList={fault};
diffList={diff};
stepsList={0};
currState=x0;
startIdx=1;
else
[busBase,swBase,pvBase,pqBase,shuntBase,lineBase,indBase,zipBase,synBase,excBase,tgBase,agcBase,cacBase,clusterBase]=unfoldSysData(SysDataBase);
if ~isfield(snapshot,'maps')
maps.busMap=(1:size(busBase,1))';maps.invBusMap=maps.busMap;
maps.indMap=(1:size(indBase,1))';maps.invIndMap=maps.indMap;
maps.synMap=(1:size(synBase,1))';maps.invSynMap=maps.synMap;
maps.lineMap=(1:size(lineBase,1))';maps.invLineMap=maps.lineMap;
snapshot.maps=maps;
end
if ~isfield(snapshot,'fault')
fault=zeros(0,4);
snapshot.fault=fault;
end
if ~isfield(snapshot,'diff')
diff=0;
snapshot.diff=diff;
end
if ~isfield(snapshot,'t')
snapshot.t=0;
end
if ~isfield(snapshot,'step')
snapshot.step=0;
end
currState=snapshot.state(:,end);
SysData=snapshot.SysData;
t=0;
maps=snapshot.maps;
fault=snapshot.fault;
diff=snapshot.diff;
startIdx=1;
SnapshotLilst={snapshot};
SysDataList={SysData};
stateList={currState};
tList={0};
mapsList={maps};
faultList={fault};
diffList={diff};
stepsList={0};
end
item=startIdx;
SysDataOrig=[];
lastSwitchT=0.0;
switchDelay=options.switchDelay;
notifyCtrlDiff=options.notifyCtrlDiff;
while initSuccess&&item<=size(eventList,1)
event=eventList(item,:);
eventType=event(4);
if eventType==EvtType.START %start / black start
elseif eventType==EvtType.ADD_LINE %add line
pos=event(5);
method=event(6);
evtLineSt=evtLine(pos,2);evtLineEnd=evtLine(pos,3);
eventSpec=evtLineSpec(evtLineSt:evtLineEnd,:);
eventSpec=eventSpec(eventSpec(:,2)==0,:);
lineAdd=lineBase(eventSpec(:,1),:);
lineBusNew=[maps.busMap(lineAdd(:,1)),maps.busMap(lineAdd(:,2))];
lineAddConnIdx=eventSpec(lineBusNew(:,1)~=0&lineBusNew(:,2)~=0,1);
lineAddNewIdx=eventSpec((lineBusNew(:,1)==0&lineBusNew(:,2)~=0)|(lineBusNew(:,1)~=0&lineBusNew(:,2)==0),1);
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
[stateNew,finalAlpha,alphaList,diff,SysData,maps]=...
simulationAddConnectionLine(SysData,SysDataBase,SysPara,SimData,maps,currState,lineAddConnIdx);
if finalAlpha<1-tTol
msg=['[',num2str(event(2)),']',' Add connection line causes singularity. Exit'];
flag=-1;
addLog(['Add connection line causes singularity. Exit'],'INFO');
break;
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
[stateUpd,finalAlpha,alphaList,diff,SysData,maps]=...
simulationAddNewLine(SysData,SysDataBase,SysPara,SimData,maps,stateNew,lineAddNewIdx);
if finalAlpha<1-tTol
msg=['[',num2str(event(2)),']',' Add new line causes singularity. Exit'];
flag=-1;
addLog(['Add new line causes singularity. Exit'],'INFO');
break;
end
snapshot.SysData=SysData;
snapshot.state=stateUpd;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateUpd;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateUpd;
addLog(['[E',num2str(event(1)),'] Added Line(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.ADD_LOAD %add static load
pos=event(5);
method=event(6);
evtZipSt=evtZip(pos,3);evtZipEnd=evtZip(pos,4);
if evtZip(pos,2)==0
eventSpec=evtZipSpec(evtZipSt:evtZipEnd,:);
eventSpec=eventSpec(eventSpec(:,2)==0,:);
zipAdd=zipBase(eventSpec(:,1),:);
else
zipAdd=evtZipSpec2(evtZipSt:evtZipEnd,:);
eventSpec=zipAdd;
end
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationAddZipLoadInstant(busTag,SysData,SysDataBase,SysPara,SimData,maps,currState,zipAdd);
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Add static load',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Added Zip(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.ADD_IND %add dynamic load
pos=event(5);
method=event(6);
evtIndSt=evtInd(pos,2);evtIndEnd=evtInd(pos,3);
eventSpec=evtIndSpec(evtIndSt:evtIndEnd,:);
eventSpecMod=eventSpec(eventSpec(:,2)==1,:);
eventSpec=eventSpec(eventSpec(:,2)==0,:);
TmAdd=eventSpec(:,3);
s0Add=eventSpec(:,4);
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationAddMotorInstant(busTag,SysData,SysDataBase,SysPara,SimData,maps,currState,eventSpec(:,1),TmAdd,s0Add);
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Add induction motor',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
if ~isempty(eventSpecMod)
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
indModIdx=maps.indMap(eventSpecMod(:,1));
eventSpecMod=eventSpecMod(indModIdx~=0,:);
indModIdx=indModIdx(indModIdx~=0);
ind(indModIdx,15)=eventSpecMod(:,3);
ind(indModIdx,16:17)=0;
SysData=foldSysData(bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster);
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Added Ind(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.ADD_SYN %add syn gen
pos=event(5);
method=event(6);
evtSynSt=evtSyn(pos,2);evtSynEnd=evtSyn(pos,3);
eventSpec=evtSynSpec(evtSynSt:evtSynEnd,:);
eventSpec=eventSpec(eventSpec(:,2)==0,:);
d0Add=eventSpec(:,3);
Pm0=eventSpec(:,4);
Ef0=eventSpec(:,5);
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationAddSynInstant(busTag,SysData,SysDataBase,SysPara,SimData,maps,...
currState,eventSpec(:,1),d0Add,Pm0,Ef0);
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Add syn gen',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Added Syn(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.DYN_SIM %dyn simulation
pos=event(5);
method=event(6);
AEmethod=round(10*mod(method,1));
maxTime=event(3)-event(2);
dt=event(7);
if ~isempty(SysDataOrig); dt=min([event(7)*50,1.0]);end
if item>1&&eventList(item-1,4)==EvtType.DYN_SIM; dt=event(7)*10;end
if dt<SimSet.MIN_DT;dt=SimSet.MIN_DT;end
if dt>SimSet.MAX_DT;dt=SimSet.MAX_DT;end
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
nline=size(line,1);
nInd=size(ind,1);
nZip=size(zip,1);
nSyn=size(syn,1);
nTg=size(tg,1);
nExc=size(exc,1);
nAgc=size(agc,1);
nCac=size(cac,1);
nCluster=size(cluster,1);
synToTg=zeros(size(syn,1),1);
synToTg(tg(:,1))=1:size(tg,1);
synToExc=zeros(size(syn,1),1);
synToExc(exc(:,1))=1:size(exc,1);
switchMask=eval(switchMaskCondition);
if pos==0
% switchMask=switchMask;
lastDiffCtrl=0;
if item>1&&eventList(item-1,4)==EvtType.DYN_SIM&&eventList(item-1,6)==MethodType.RK4_NR
lastDiffCtrl=1;
end
if useDiffCtrl&&AEmethod==0&&~lastDiffCtrl&&((diff(end)>diffTolCtrl&&maxTime>1)||notifyCtrlDiff)
eventListTmp=zeros(size(eventList,1)+1,size(eventList,2));
eventListTmp(1:item,:)=eventList(1:item,:);
if eventList(item-1,4)==EvtType.DYN_SIM
eventListTmp(item,6)=MethodType.RK4_NR;
eventListTmp(item,7)=SimSet.DEFAULT_DT/100;
duration=eventListTmp(item,7)*20;
if duration>1;duration=1;end
else
duration=maxTime/5;
if duration>2;duration=2;end
end
eventListTmp((item+1):end,:)=eventList(item:end,:);
eventListTmp(item,3)=eventListTmp(item,2)+duration;
eventListTmp(item+1,2)=eventListTmp(item,2)+duration;
eventList=eventListTmp;
event=eventList(item,:);
method=event(6);
AEmethod=round(10*mod(method,1));
dt=event(7);
maxTime=event(3)-event(2);
diffTolSim=diffTolOrig;
end
varOpt.allowModelSwitch=switchMask & allowSteadyDynSwitch;
varOpt.absT=event(2);
varOpt.lastSwitchT=lastSwitchT;
varOpt.switchDelay=switchDelay;
varOpt.useDiffCtrl=useDiffCtrl;
varOpt.diffTolCtrl=diffTolCtrl;
varOpt.useNomoveCtrl=useNoMoveCtrl|(~isempty(SysDataOrig)&switchMask & allowSteadyDynSwitch);
if ~isempty(SysDataOrig)
SimData=foldSimData(maxTime,5,dt,nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method,varOpt);
else
SimData=foldSimData(maxTime,1,dt,nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method,varOpt);
end
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
zipRampSpec=[];
else
% switchMask=switchMask;
ePara=evtDyn(pos,:);
pqSt=ePara(2);pqEnd=ePara(3);pvSt=ePara(4);pvEnd=ePara(5);
indSt=ePara(6);indEnd=ePara(7);zipSt=ePara(8);zipEnd=ePara(9);
shSt=ePara(10);shEnd=ePara(11);zipRampSt=ePara(12);zipRampEnd=ePara(13);
TmechSt=ePara(14);TmechEnd=ePara(15);PmSt=ePara(16);PmEnd=ePara(17);
EfSt=ePara(18);EfEnd=ePara(19);VrefSt=ePara(20);VrefEnd=ePara(21);Eq1St=ePara(22);Eq1End=ePara(23);
if pqSt>0&&pqEnd>0;pqSpec=evtDynPQ(pqSt:pqEnd,:);pqSpec(:,1)=maps.busMap(pqSpec(:,1));pqSpec=pqSpec(pqSpec(:,1)~=0,:);else pqSpec=zeros(0,3);end
if indSt>0&&indEnd>0;indSpec=evtDynInd(indSt:indEnd,:);indSpec(:,1)=maps.indMap(indSpec(:,1));indSpec=indSpec(indSpec(:,1)~=0,:);else indSpec=zeros(0,2);end
if zipSt>0&&zipEnd>0;zipSpec=evtDynZip(zipSt:zipEnd,:);zipSpec(:,1)=maps.busMap(zipSpec(:,1));zipSpec=zipSpec(zipSpec(:,1)~=0,:);else zipSpec=zeros(0,2);end
if shSt>0&&shEnd>0;shSpec=evtDynSh(shSt:shEnd,:);shSpec(:,1)=maps.busMap(shSpec(:,1));shSpec=shSpec(shSpec(:,1)~=0,:);else shSpec=zeros(0,2);end
if zipRampSt>0&&zipRampEnd>0;zipRampSpec=evtDynZipRamp(zipRampSt:zipRampEnd,:);zipRampSpec(:,1)=maps.busMap(zipRampSpec(:,1));zipRampSpec=zipRampSpec(zipRampSpec(:,1)~=0,:);else zipRampSpec=zeros(0,12);end
TmechSpecToPm=[];
TmechSpecToPv=[];
if TmechSt>0&&TmechEnd>0
if isempty(SysDataOrig) % full dynamic model
TmechSpec=evtDynTmech(TmechSt:TmechEnd,:);
TmechSpecIdx=synToTg(maps.synMap(TmechSpec(:,1)));
TmechSpecToPm=TmechSpec(TmechSpecIdx==0,:);
TmechSpecToPm(:,1)=maps.synMap(TmechSpecToPm(:,1));
TmechSpecToPm=TmechSpecToPm(TmechSpecToPm(:,1)~=0,:);
TmechSpec=TmechSpec(TmechSpecIdx~=0,:);
TmechSpec(:,1)=TmechSpecIdx;
else % QSS model
TmechSpecToPv=evtDynTmech(TmechSt:TmechEnd,:);
TmechSpecToPv(:,1)=synMapOrig(TmechSpecToPv(:,1));
TmechSpecToPv=TmechSpecToPv(TmechSpecToPv(:,1)~=0,:);
TmechSpecToPv(:,1)=SysDataOrig.syn(TmechSpecToPv(:,1),1);
TmechSpec=zeros(0,2);
end
else
TmechSpec=zeros(0,2);
end
PmSpecToPv=[];
if PmSt>0&&PmEnd>0
if isempty(SysDataOrig) % full dynamic model
PmSpec=evtDynPm(PmSt:PmEnd,:);
PmSpec(:,1)=maps.synMap(PmSpec(:,1));
PmSpec=PmSpec(PmSpec(:,1)~=0,:);
else % QSS model
PmSpecToPv=evtDynPm(PmSt:PmEnd,:);
PmSpecToPv(:,1)=synMapOrig(PmSpecToPv(:,1));
PmSpecToPv=PmSpecToPv(PmSpecToPv(:,1)~=0,:);
PmSpecToPv(:,1)=SysDataOrig.syn(PmSpecToPv(:,1),1);
PmSpec=zeros(0,2);
end
else
PmSpec=zeros(0,2);
end
if ~isempty(TmechSpecToPm)&&isempty(SysDataOrig)
PmSpecFull=zeros(size(syn,1),2);
PmSpecFull(:,1)=1:size(syn,1);
PmSpecFull(PmSpec(:,1),2)=PmSpec(:,2);
PmSpecFull(TmechSpecToPm(:,1),2)=PmSpecFull(TmechSpecToPm(:,1),2)+TmechSpecToPm(:,2);
PmSpec=PmSpecFull(PmSpecFull(:,2)~=0,:);
end
if pvSt>0&&pvEnd>0
pvSpec=evtDynPV(pvSt:pvEnd,:);
pvSpec(:,1)=maps.busMap(pvSpec(:,1));
pvSpec=pvSpec(pvSpec(:,1)~=0,:);
else
pvSpec=zeros(0,2);
end
if ~isempty(PmSpecToPv)
pvSpecFull=zeros(size(bus,1),2);
pvSpecFull(:,1)=1:size(bus,1);
pvSpecFull(pvSpec(:,1),2)=pvSpec(:,2);
pvSpecFull(PmSpecToPv(:,1),2)=PmSpecToPv(:,2);
pvSpec=pvSpecFull(pvSpecFull(:,2)~=0,:);
end
if ~isempty(TmechSpecToPv)
pvSpecFull=zeros(size(bus,1),2);
pvSpecFull(:,1)=1:size(bus,1);
pvSpecFull(pvSpec(:,1),2)=pvSpec(:,2);
pvSpecFull(TmechSpecToPv(:,1),2)=TmechSpecToPv(:,2);
pvSpec=pvSpecFull(pvSpecFull(:,2)~=0,:);
end
if EfSt>0&&EfEnd>0;EfSpec=evtDynEf(EfSt:EfEnd,:);EfSpec(:,1)=maps.synMap(EfSpec(:,1));EfSpec=EfSpec(EfSpec(:,1)~=0,:);else EfSpec=zeros(0,2);end
if VrefSt>0&&VrefEnd>0;VrefSpec=evtDynVref(VrefSt:VrefEnd,:);VrefSpec(:,1)=synToExc(maps.synMap(VrefSpec(:,1)));VrefSpec=VrefSpec(VrefSpec(:,1)~=0,:);else VrefSpec=zeros(0,2);end
if Eq1St>0&&Eq1End>0;Eq1Spec=evtDynEq1(Eq1St:Eq1End,:);Eq1Spec(:,1)=maps.synMap(Eq1Spec(:,1));Eq1Spec=Eq1Spec(Eq1Spec(:,1)~=0,:);else Eq1Spec=zeros(0,2);end
[~,Ytr0,Ysh0,~,~,~,~]=getYMatrix(nbus,line,fault);
pqBus=zeros(nbus,size(pqSpec,2)-1);pvBus=zeros(nbus,1);
pqBus(pqSpec(:,1),:)=pqSpec(:,2:end);pvBus(pvSpec(:,1),:)=pvSpec(:,2);
pqIncr=pqBus(pq(:,1),:);pvIncr=pvBus(pv(:,1),:);
zipBus=zeros(nbus,1);zipBus(zipSpec(:,1))=zipSpec(:,2);
Rzip0=ones(nZip,1);Rzip1=zipBus(zip(:,1),:);
Rind0=ones(nInd,1);Rind1=zeros(nInd,1);Rind1(indSpec(:,1),:)=indSpec(:,2);
Ysh1=zeros(nbus,1);Ysh1(shSpec(:,1),:)=shSpec(:,2);
Tmech1=zeros(nTg,1);Tmech1(TmechSpec(:,1),:)=TmechSpec(:,2);
Varref1=zeros(nExc,1);Varref1(VrefSpec(:,1),:)=VrefSpec(:,2);
Ef1=zeros(nSyn,1);Ef1(EfSpec(:,1),:)=EfSpec(:,2);
Pm1=zeros(nSyn,1);Pm1(PmSpec(:,1),:)=PmSpec(:,2);
Eq11=zeros(nSyn,1);Eq11(Eq1Spec(:,1),:)=Eq1Spec(:,2);
zipRampSpec(:,1)=maps.invBusMap(zipRampSpec(:,1));
lastDiffCtrl=0;
if item>1&&eventList(item-1,4)==EvtType.DYN_SIM&&eventList(item-1,6)==MethodType.RK4_NR
lastDiffCtrl=1;
end
if useDiffCtrl&&AEmethod==0&&~lastDiffCtrl&&((diff(end)>diffTolCtrl&&maxTime>1)||notifyCtrlDiff)
eventListTmp=zeros(size(eventList,1)+1,size(eventList,2));
eventListTmp(1:item,:)=eventList(1:item,:);
if eventList(item-1,4)==EvtType.DYN_SIM
eventListTmp(item,6)=MethodType.RK4_NR;
eventListTmp(item,7)=SimSet.DEFAULT_DT/100;
duration=eventListTmp(item,7)*20;
if duration>1;duration=1;end
else
duration=maxTime/5;
if duration>2;duration=2;end
end
eventListTmp(item,3)=event(2)+duration;
eventListTmp((item+1):end,:)=eventList(item:end,:);
eventListTmp(item+1,2)=event(2)+duration;
eventList=eventListTmp;
event=eventList(item,:);
if eventList(item+1,5)~=0
evtDynAdd=evtDyn(event(5),:);
if any(evtDynAdd(:,[6,7,8,9])>0)
if all(evtDynAdd(:,[6,7])>0) % Has induction motor item
dynIndAdd=evtDynInd(evtDynAdd(6):evtDynAdd(7),:);
dynIndAdd(:,2)=dynIndAdd(:,2)./(1+duration*dynIndAdd(:,2));
evtDynAdd(6)=size(evtDynInd,1)+1;
evtDynAdd(7)=size(evtDynInd,1)+size(dynIndAdd,1);
evtDynInd=[evtDynInd;dynIndAdd];
end
if all(evtDynAdd(:,[8,9])>0)
dynZipAdd=evtDynZip(evtDynAdd(8):evtDynAdd(9),:);
dynZipAdd(:,2)=dynZipAdd(:,2)./(1+duration*dynZipAdd(:,2));
evtDynAdd(8)=size(evtDynZip,1)+1;
evtDynAdd(9)=size(evtDynZip,1)+size(dynZipAdd,1);
evtDynZip=[evtDynZip;dynZipAdd];
end
evtDynAdd(:,1)=(size(evtDyn,1)+1):(size(evtDyn,1)+size(evtDynAdd,1));
eventList((item+1):(item+size(evtDynAdd,1)),5)=evtDynAdd(:,1);
evtDyn=[evtDyn;evtDynAdd];
end
end
method=event(6);
AEmethod=round(10*mod(method,1));
dt=event(7);
maxTime=event(3)-event(2);
diffTolSim=diffTolOrig;
end
varOpt.allowModelSwitch=switchMask & allowSteadyDynSwitch;
varOpt.absT=event(2);
varOpt.lastSwitchT=lastSwitchT;
varOpt.switchDelay=switchDelay;
varOpt.allowExit1=0;
varOpt.useDiffCtrl=useDiffCtrl;
varOpt.diffTolCtrl=diffTolCtrl;
varOpt.useNomoveCtrl=useNoMoveCtrl|(~isempty(SysDataOrig)&switchMask & allowSteadyDynSwitch);
if ~isempty(SysDataOrig)
SimData=foldSimData(maxTime,5,dt,nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method,varOpt);
else
SimData=foldSimData(maxTime,1,dt,nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method,varOpt);
end
SysPara=foldSysPara(pqIncr,pvIncr,Rind0,Rind1,[],[],Rzip0,Rzip1,Ytr0,0*Ytr0,Ysh0,Ysh1,[],[],[],[],[],Tmech1,Varref1,Ef1,Pm1,Eq11,fault);
end
if switchMask && allowSteadyDynSwitch && pos==0 && exist('skipTag','var') && skipTag==1
stateCurve=currState;
t=0.0;
finalAlpha=0.0;
alphaList=maxTime;
diff=0.0;
SysDataTmp=SysData;
exitFlag=DynSimFlag.STEADY;
else
skipTag=0;
[stateCurve,t,finalAlpha,alphaList,diff,SysDataTmp,exitFlag]=generalDynSimulation(SimData,SysData,SysPara,currState,zipRampSpec,maps.busMap);
end
terminateFlag=0;
reverseFlag=0;
switchReady=0;
if finalAlpha<maxTime-tTol
if (switchMask && allowSteadyDynSwitch&& (exitFlag==DynSimFlag.QSS || (exitFlag==DynSimFlag.STEADY && pos~=0) || exitFlag==DynSimFlag.NOMOVE_CTRL))||(useDiffCtrl&&exitFlag==DynSimFlag.DIFF_CTRL)
evtAdd=event;
curEndTime=event(2)+finalAlpha;
oriEndTime=event(3);
event(3)=curEndTime;
eventList(item,:)=event;
if item>=size(eventList,1)
evtAdd(1)=event(1)+1;
else
evtAdd(1)=(event(1)+eventList(item+1,1))/2;
end
evtAdd(2)=curEndTime;
evtAdd(3)=oriEndTime;
if evtAdd(5)~=0
evtDynAdd=evtDyn(evtAdd(5),:);
if any(evtDynAdd(:,[6,7,8,9])>0)
if all(evtDynAdd(:,[6,7])>0) % Has induction motor item
dynIndAdd=evtDynInd(evtDynAdd(6):evtDynAdd(7),:);
dynIndAdd(:,2)=dynIndAdd(:,2)./(1+finalAlpha*dynIndAdd(:,2));
evtDynAdd(6)=size(evtDynInd,1)+1;
evtDynAdd(7)=size(evtDynInd,1)+size(dynIndAdd,1);
evtDynInd=[evtDynInd;dynIndAdd];
end
if all(evtDynAdd(:,[8,9])>0)
dynZipAdd=evtDynZip(evtDynAdd(8):evtDynAdd(9),:);
dynZipAdd(:,2)=dynZipAdd(:,2)./(1+finalAlpha*dynZipAdd(:,2));
evtDynAdd(8)=size(evtDynZip,1)+1;
evtDynAdd(9)=size(evtDynZip,1)+size(dynZipAdd,1);
evtDynZip=[evtDynZip;dynZipAdd];
end
evtDynAdd(:,1)=(size(evtDyn,1)+1):(size(evtDyn,1)+size(evtDynAdd,1));
evtAdd(:,5)=evtDynAdd(:,1);
evtDyn=[evtDyn;evtDynAdd];
end
end
eventList=[eventList(1:item,:);evtAdd;eventList((item+1):end,:)];
if exitFlag==DynSimFlag.QSS || (exitFlag==DynSimFlag.STEADY && pos~=0)
switchReady=1;
end
if exitFlag==DynSimFlag.DIFF_CTRL
notifyCtrlDiff=1;
end
elseif allowSteadyDynSwitch&& (exitFlag==DynSimFlag.STEADY && pos==0)
alphaList=[alphaList,maxTime-t(end)];
diff=[diff,0.0];
stateCurve=[stateCurve,stateCurve(:,end)];
t=[t,maxTime];
skipTag=1;
addLog(['Simulation directly jumped to maxTime=',num2str(maxTime)],'INFO');
else
if method==0.0||~allowReverse
terminateFlag=1;
else
eventList(item,6)=0.0;
item=item-1;
reverseFlag=1;
end
end
end
if ~reverseFlag
SysData=SysDataTmp;
snapshot.SysData=SysData;
snapshot.state=stateCurve;
snapshot.t=t;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=alphaList;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysDataTmp;
stateList{end+1}=stateCurve;tList{end+1}=t;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=alphaList;
currState=stateCurve(:,end);
if allowSteadyDynSwitch
deltax=(stateCurve(:,end)-stateCurve(:,end-1))/(t(end)-t(end-1));
if switchMask && isempty(SysDataOrig) && exitFlag==DynSimFlag.QSS % Dynamic model -> QSS model
if switchReady && eventList(item+1,4)==EvtType.DYN_SIM
SysDataOrig=SysData;
stateOrig=stateCurve(:,end);
mapsOrig=maps;
[SysData,currState]=convertSteadyDynModels(SysData,[],currState,deltax,[],SysPara,1);
maps.synMap=zeros(size(synBase,1),1);maps.invSynMap=zeros(0,1);
else
lastSwitchT=event(3);
end
elseif (~isempty(SysDataOrig) &&...
(eventList(item+1,4)~=EvtType.DYN_SIM||(item<size(eventList,1)-1&&eventList(item+2,4)==EvtType.DYN_SIM&&eventList(item+2,5)~=0&&eventList(item+2,3)-eventList(item+2,2)<10)))||...
exitFlag==DynSimFlag.NOMOVE_CTRL % Steady-state model -> dynamic model
[SysData,currState]=convertSteadyDynModels(SysData,SysDataOrig,currState,deltax,stateOrig,SysPara,0);
maps=mapsOrig;
SysDataOrig=[];
lastSwitchT=event(3);
end
end
end
if eventList(item+1,4)~=EvtType.DYN_SIM
skipTag=0;
lastSwitchT=event(3);
notifyCtrlDiff=0;
end
if terminateFlag
msg=['[',num2str(event(2)+finalAlpha),']','Simulation terminated earlier. Exit'];
flag=-1;
addLog(['Simulation terminated earlier. Exit'],'INFO');
break;
end
elseif eventType==EvtType.FAULT % fault
pos=event(5);
method=event(6);
evtFaultSt=evtFault(pos,2);evtFaultEnd=evtFault(pos,3);
eventSpec=evtFaultSpec(evtFaultSt:evtFaultEnd,:);
addFault=eventSpec(eventSpec(:,5)==0,1:4);
clearFault=eventSpec(eventSpec(:,5)~=0,1:4);
addFault(:,1)=maps.lineMap(addFault(:,1));
clearFault(:,1)=maps.lineMap(clearFault(:,1));
faultPrev=fault;
faultTag=zeros(size(maps.lineMap));
faultTag(maps.invLineMap(fault(:,1)))=1;
faultTag(maps.invLineMap(clearFault(:,1)))=0;
fault=fault(faultTag(maps.invLineMap(fault(:,1)))==1,:);
addFault=addFault(faultTag(maps.invLineMap(addFault(:,1)))~=1,:);
fault=[fault;addFault];
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
while 1
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,faultPrev,fault,mapsTemp]=...
simulationChangeFault(busTag,SysData,SysDataBase,SysPara,SimData,maps,currState,faultPrev,fault);
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Changing fault status',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Applied fault(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.CUT_LINE %
pos=event(5);
method=event(6);
evtLineSt=evtLine(pos,2);evtLineEnd=evtLine(pos,3);
eventSpec=evtLineSpec(evtLineSt:evtLineEnd,:);
eventSpec=eventSpec(eventSpec(:,2)==1,:);
lineCutIdx=eventSpec(:,1);
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationCutLine(busTag,SysData,SysDataBase,SysPara,SimData,maps,currState,lineCutIdx);
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Cut line',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Cut line(s) ',num2str(lineCutIdx')],'INFO');
elseif eventType==EvtType.CUT_LOAD % Not necessary, cut load is adding negative load.
elseif eventType==EvtType.CUT_IND %
pos=event(5);
method=event(6);
evtIndSt=evtInd(pos,2);evtIndEnd=evtInd(pos,3);
eventSpec=evtIndSpec(evtIndSt:evtIndEnd,:);
eventSpecMod=eventSpec(eventSpec(:,2)==1,:);
eventSpec=eventSpec(eventSpec(:,2)==2,:);
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationCutMotorInstant(busTag,SysData,SysDataBase,SysPara,SimData,maps,currState,eventSpec(:,1));
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Cut induction motor',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
if ~isempty(eventSpecMod)
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
indModIdx=maps.indMap(eventSpecMod(:,1));
eventSpecMod=eventSpecMod(indModIdx~=0,:);
indModIdx=indModIdx(indModIdx~=0);
ind(indModIdx,15)=eventSpecMod(:,3);
ind(indModIdx,16:17)=0;
SysData=foldSysData(bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster);
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Cut ind(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.CUT_SYN %
pos=event(5);
method=event(6);
evtSynSt=evtSyn(pos,2);evtSynEnd=evtSyn(pos,3);
eventSpec=evtSynSpec(evtSynSt:evtSynEnd,:);
eventSpec=eventSpec(eventSpec(:,2)==1,:);
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationCutSynInstant(busTag,SysData,SysDataBase,SysPara,SimData,maps,...
currState,eventSpec(:,1));
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Cut synchronous generator',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Cut syn(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.ALT_SYN %
pos=event(5);
method=event(6);
evtSynSt=evtSyn(pos,2);evtSynEnd=evtSyn(pos,3);
eventSpec=evtSynSpec(evtSynSt:evtSynEnd,:);
eventSpec=eventSpec(eventSpec(:,2)==2,:);
if ~isempty(eventSpec)
dSynCompact=eventSpec(:,3:9);
newSynIdx=maps.synMap(eventSpec(:,1));
activeSynIdx=newSynIdx(newSynIdx~=0);
deltaSyn=zeros(size(syn,1),size(dSynCompact,2));
deltaSyn(activeSynIdx,:)=dSynCompact;
else
deltaSyn=zeros(size(syn,1),7);
end
if diff(end)>diffTol;diffTolSim=diff(end)+1e-9;else diffTolSim=diffTol;end
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
nbus=size(bus,1);
busTag=ones(nbus,1);
SimData=foldSimData([],1,[],nlvl,taylorN,alphaTol,diffTolSim,diffTolMax,method);
doContinue=1;
SysPara=foldSysPara([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],fault);
synNew=syn;
synNew(:,[7:10,13:15])=synNew(:,[7:10,13:15])+deltaSyn;
while 1
[stateNew,finalAlpha,alphaList,diff,SysDataTemp,SysPara,mapsTemp]=...
simulationAlterSynParameter(busTag,SysData,SysDataBase,SysPara,SimData,maps,currState,synNew);
if finalAlpha<1-tTol
[flag,doContinue,willBreak,msg]=processCollaps(SysDataTemp,stateNew(:,end),event,'Alter syn parameters',flag,doContinue);
if willBreak
break;
end
else
SysData=SysDataTemp;
maps=mapsTemp;
break;
end
end
if ~doContinue
break;
end
snapshot.SysData=SysData;
snapshot.state=stateNew;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
SnapshotLilst{end+1}=snapshot;
SysDataList{end+1}=SysData;
stateList{end+1}=stateNew;tList{end+1}=0;
mapsList{end+1}=maps;
faultList{end+1}=fault;
diffList{end+1}=diff;stepsList{end+1}=0;
currState=stateNew;
addLog(['[E',num2str(event(1)),'] Alter syn(s) ',num2str(eventSpec(:,1)')],'INFO');
elseif eventType==EvtType.END % end
msg=['[',num2str(event(2)),'] ','Simulation successfully finished.'];
flag=1;
snapshot.SysData=SysData;
snapshot.state=currState;
snapshot.t=0;
snapshot.maps=maps;
snapshot.fault=fault;
snapshot.diff=diff;
snapshot.step=0;
break;
else
% Undefined event types, do nothing
end
item=item+1;
end
timeSpent=toc(startTag);
addLog(['Total computation time:',num2str(timeSpent),'s.'],'INFO');
% disp(tp)
[t,stateCurve]=reformatCurves(SysDataBase,tList,stateList,SysDataList,mapsList);
clearAllTempFiles();
if options.output
save(['restoration/',caseName,'_',timestamp,'_simp.mat'],...
'caseName','simSettings','eventList','SysDataBase','t','stateCurve','timeSpent','comment',...
'snapshot','-v7.3');
addLog(['Curves saved at: ','restoration/',caseName,'_',timestamp,'_simp.mat'],'INFO');
end
if options.output&&~options.simpleOutput
save(['restoration/',caseName,'_',timestamp,'.mat'],...
'caseName','simSettings','eventList','SysDataBase',...
'SnapshotLilst','tList','stateList',...
'timeSpent','comment','snapshot','-v7.3');
addLog(['More detailed results saved at: ','restoration/',caseName,'_',timestamp,'.mat'],'INFO');
end
addLog('================== END ===================','INFO');
flushLogs();
res.flag=flag;
res.msg=msg;
res.t=t;
res.stateCurve=stateCurve;
res.simSettings=simSettings;
res.SysDataBase=SysDataBase;
res.caseName=caseName;
res.eventList=eventList;
res.timestamp=timestamp;
res.snapshot=snapshot;
end
function [flag,doContinue,willBreak,msg]=processCollaps(SysData,state,event,actionStr,flag,doContinue)
willBreak=0;
msg='';
[bus,sw,pv,pq,shunt,line,ind,zip,syn,exc,tg,agc,cac,cluster]=unfoldSysData(SysData);
[V,~,s,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~]=unfoldX(state(:,end),SysData);
[jac,jPIdx,jQIdx,jVIdx,jSIdx,busType]=getJacobianMod(size(bus,1),size(line,1),bus,sw,pv,pq,shunt,line,zip,ind,s,V,ones(size(ind,1)));
[Vx,Dx]=eig(full(jac));
dDx=diag(Dx);
[~,iSmall]=min(abs(dDx));
part=Vx(:,iSmall).*conj(Vx(:,iSmall));
partBus=zeros(size(bus,1),1);
partBus(busType~=2)=partBus(busType~=2)+part(jPIdx);
partBus(busType==0)=partBus(busType==0)+part(jQIdx);
[nIslands,islands,refs]=searchIslands(bus(:,1),line(:,1:2));
if nIslands==1
msg=['[',num2str(event(2)),'] ',actionStr,' causes singularity. Exit'];
flag=-1;
addLog([actionStr,' causes singularity. Exit'],'INFO');
doContinue=0;
willBreak=1;
else
addLog(['Part of the system collapses, purge collapsed parts.'],'INFO');
partIsland=accumarray(islands,partBus,[nIslands,1]);
[maxPartIsland,iMaxPart]=max(partIsland);
if maxPartIsland>0.8
busFilter=find(busTag==1);
busTag(busFilter(islands==iMaxPart(1)))=0;
else
msg=['[',num2str(event(2)),'] ','No dominant collapse island can be purged. Exit'];
flag=-1;
addLog(['No dominant collapse island can be purged. Exit'],'INFO');
doContinue=0;
willBreak=1;
end
end
end
function options=regulateOptions(options)
if nargin<1
options=[];
end
if ~isfield(options,'dbstop');options.dbstop=0;end
if ~isfield(options,'output');options.output=0;end
if ~isfield(options,'comment');options.comment='';end
if ~isfield(options,'tTol');options.tTol=1e-9;end
if ~isfield(options,'allowReverse');options.allowReverse=0;end
if ~isfield(options,'allowSteadyDynSwitch');options.allowSteadyDynSwitch=0;end
if ~isfield(options,'useDiffCtrl');options.useDiffCtrl=0;end
if ~isfield(options,'useNoMoveCtrl');options.useNoMoveCtrl=0;end
if ~isfield(options,'nPool');options.nPool=1;end
if ~isfield(options,'hotStart');options.hotStart=0;end
if ~isfield(options,'switchDelay');options.switchDelay=3.0;end
if ~isfield(options,'notifyCtrlDiff');options.notifyCtrlDiff=0;end
if ~isfield(options,'inheritParaSettings');options.inheritParaSettings=1;end
if ~isfield(options,'nlvl');options.nlvl=15;end
if ~isfield(options,'taylorN');options.taylorN=4;end
if ~isfield(options,'segAlpha');options.segAlpha=1;end
if ~isfield(options,'dAlpha');options.dAlpha=1;end
if ~isfield(options,'alphaTol');options.alphaTol=1.0000e-03;end
if ~isfield(options,'diffTol');options.diffTol=1.0000e-06;end
if ~isfield(options,'diffTolMax');options.diffTolMax=1.0000e-02;end
if ~isfield(options,'method');options.method=0;end
if ~isfield(options,'diffTolCtrl');options.diffTolCtrl=1.0000e-05;end
if ~isfield(options,'simpleOutput');options.simpleOutput=1;end
end
Reference in new issue View Git Blame Copy Permalink