diff --git a/calCPF.m b/calCPF.m index dda7889..4b0a2fc 100644 --- a/calCPF.m +++ b/calCPF.m @@ -64,7 +64,12 @@ else snapshot=[]; hotStart=0; 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.evtDynPV=pvRamp; simSettings.evtDyn=zeros(1,23); simSettings.evtDyn(1,1)=1; diff --git a/calTSA.m b/calTSA.m index 75e1194..02900db 100644 --- a/calTSA.m +++ b/calTSA.m @@ -78,10 +78,14 @@ else nAggFault=size(aggFaultIdxAux,1)-1; simSettings.evtFault=[(1:nAggFault)',aggFaultIdxAux(1:nAggFault),aggFaultIdxAux(2:nAggFault+1)-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),... 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])]]; + + + res=runDynamicSimulationExec(caseName,SysData,simSettings,options,snapshot); end diff --git a/internal/calculateInitialState.m b/internal/calculateInitialState.m index a9a2d22..c746503 100644 --- a/internal/calculateInitialState.m +++ b/internal/calculateInitialState.m @@ -43,6 +43,7 @@ function [SysDataUpd,x0,finalAlpha,alphaList,diff,t,stateCurve]=calculateInitial MatGV0,MatGV1,MatGRhs0,MatGRhs1,Tmech1,Varref1,Ef1,Pm1,Eq11]=unfoldSysPara(SysPara); ltc=[];sup=[];dem=[];busName=[]; % 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])); if nIslands>1 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] ... =forgeLargestIsland(maxIsland,nIslands,islands,bus,[],pq,sw,line,ltc,sup,dem,shunt,ind,zip,syn,busName); end - +if ~isempty(Pm1) + pm=Pm1; +end nbus=size(bus,1); nline=size(line,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'); return; 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); -if ~isempty(Pm1) - pm=Pm1; -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); nbus=size(bus,1); nline=size(line,1); % nlvl=10; diff --git a/internal/formatEventList.m b/internal/formatEventList.m index 23a009c..034f315 100644 --- a/internal/formatEventList.m +++ b/internal/formatEventList.m @@ -91,8 +91,15 @@ if ~isempty(sstTime) 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]; [~,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,:); diff --git a/internal/generalDynSimulation.m b/internal/generalDynSimulation.m index 122de16..9e26ccf 100644 --- a/internal/generalDynSimulation.m +++ b/internal/generalDynSimulation.m @@ -118,6 +118,10 @@ if DEmethod==0 [stateUpd,t,finalAlpha,alphaList,diffList,exitFlag]=simulationTimeDomainHem(SimData,SysDataTmp,SysPara,x0); alphaList=cumsum(alphaList); SimDataUpd=SimData; +elseif DEmethod==4 + [stateUpd,t,finalAlpha,alphaList,diffList,exitFlag]=simulationTimeDomainDT(SimData,SysDataTmp,SysPara,x0); + alphaList=cumsum(alphaList); + SimDataUpd=SimData; else [stateUpd,t,diffList,nxtDt,exitFlag]=simulationTimeDomainNI(SimData,SysDataTmp,SysPara,x0); finalAlpha=t(end); diff --git a/internal/getSimMethodType.m b/internal/getSimMethodType.m index 49406f7..d647c19 100644 --- a/internal/getSimMethodType.m +++ b/internal/getSimMethodType.m @@ -30,6 +30,7 @@ function MethodType=getSimMethodType() % MethodType.FULL_HE=0.0; +MethodType.FULL_DT=4.0; MethodType.ME_HE=1.0; MethodType.ME_NR=1.1; MethodType.RK4_HE=2.0; diff --git a/internal/getseq.m b/internal/getseq.m index 25ce410..92a835a 100644 --- a/internal/getseq.m +++ b/internal/getseq.m @@ -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); aaux=[zeros(size(aux,1),1),aux,(n+d)*ones(size(aux,1),1)]; seq=aaux(:,2:end)-aaux(:,1:end-1)-1; diff --git a/internal/hemMachinePFSalientcontinueDyn.m b/internal/hemMachinePFSalientcontinueDyn.m index e720804..ef12e73 100644 --- a/internal/hemMachinePFSalientcontinueDyn.m +++ b/internal/hemMachinePFSalientcontinueDyn.m @@ -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; + + +% import system data khuang, 8 Jul [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); nline=size(line,1); + +%determine islanding, this is for AGC part khuang 8 JUl 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 +% 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); +% import simualtion data khuang 8 JUl [~,~,~,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); % +% 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 Qls=zeros(nbus,2);Qls(pq(:,1),1)=pqIncr(:,2); if size(pqIncr,2)>=4 @@ -63,10 +75,14 @@ if size(pqIncr,2)>=4 Qls(pq(:,1),2)=pqIncr(:,4); end +% formatting Ymatrix, here the default value of fault is empty, khuang 8 JUL if isempty(Ytr0) [Y,Ytr0,Ysh,ytrfr,ytrto,yshfr,yshto]=getYMatrix(nbus,line); end + +% reshape the pv, pq shunt and swing buses as zeros if they are empty +% khuang 8JUL busType=zeros(nbus,1); if isempty(pv) pv=zeros(0,6); @@ -80,33 +96,54 @@ end if isempty(sw) sw=zeros(0,13); end + +% label pv and swing buses khuang 8 JUL +% 1: PV bus, 0: PQ bus busType(pv(:,1))=1; busType(sw(:,1))=2; % 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); ipv=find(busType==1); ipq=find(busType==0); npq=size(ipq,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(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 Ysh0=Ysh0+accumarray(zip(:,1),Rzip0.*(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]); Ysh1=Ysh1+accumarray(zip(:,1),Rzip1.*(zip(:,5)+1j*zip(:,8)).*zip(:,12),[nbus,1]); end + +% now zip load + shunt khuang 8 JUL Ysh0=Ysh0+yShunt; % 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); +%initialize P and Q for every bus khuang 8 JUL pVec=zeros(nbus,1); qVec=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(pq(:,1))=pVec(pq(:,1))-pq(:,4); 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 pVec=pVec-accumarray(zip(:,1),Rzip0.*zip(:,7).*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); % 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(:,1)=V0; W=zeros(nbus,nlvl+1); W(:,1)=1./V0; + + +% initiliza magnitude of voltage V for every bus khuang 8 JUL Vmag=zeros(nbus,nlvl+1); Vmag(:,1)=abs(V0); + +% Power is initilized as we already cooked pVec and qVec khuang 8 JUL P=zeros(nbus,nlvl+1); P(:,1)=pVec; % 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); Qxtra=zeros(size(Q)); Q(:,1)=Q0; 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; 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) 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]); @@ -136,11 +191,17 @@ end % Qxtra(busType~=0,2:end)=Q(busType~=0,2:end); % 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)); D0=imag(V(:,1)); E0=real(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); D0M=sparse(1:nbus,1:nbus,D0,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); 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); 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 freqKeptTag=zeros(nbus,1); frefs=refs; @@ -177,16 +243,26 @@ end freqKeptTagxRef=freqKeptTag; freqKeptTagxRef(frefs)=0; nFreqKept=sum(freqKeptTag); +%----------------------------------------------------------- -if ~isempty(ind) - nInd=size(ind,1); - indIdx=ind(:,1); - s=zeros(nInd,nlvl+1); - s(:,1)=s0; - IL=zeros(nInd,nlvl+1); - IR=zeros(nInd,nlvl+1); - Vm=zeros(nInd,nlvl+1); + + + +%------------------------------------------------------------------ +% this part is for initialling inductor. khuang 8 JUL +if ~isempty(ind) % check if there is any inductor% khuang 8 JUL + nInd=size(ind,1); % determine the number of inductors % khuang 8 JUL + 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); X1=ind(:,8); Z1=ind(:,7)+1j*ind(:,8); @@ -197,11 +273,15 @@ if ~isempty(ind) T1=-ind(:,16)-2*ind(:,17); T2=ind(:,17); Hm=ind(:,14); - + %-----------------parameters of inductors---------% khuang 8 JUL + %-----------------END----------------% khuang 8 JUL + + Rm=zeros(nInd,1); 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)); Vm(:,1)=V0(indIdx)-IL(:,1).*Z1; IR(:,1)=Vm(:,1).*s0./(R2+1j*X2.*s0); @@ -211,6 +291,7 @@ if ~isempty(ind) JL0=real(IL(:,1)); KL0=imag(IL(:,1)); + % prepare the algebric matrix % khuang 8 JUL Yeind0=Reind0./Ze; Yeind1=Reind1./Ze; Ye1ind0=Reind0.*Z1./Ze; @@ -290,7 +371,7 @@ if ~isempty(ind) 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(:,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(:,1)=accumarray(indIdx,LHS_MatInd_Shr_sqz(:,1),[nbus,1]); LHS_MatInd_Bus_sqz(:,2)=accumarray(indIdx,LHS_MatInd_Shr_sqz(:,2),[nbus,1]); @@ -299,24 +380,29 @@ if ~isempty(ind) else s=zeros(0,nlvl+1); end +% Initialization of inductors is finished % khuang 8 JUL +%------------------------------------------------------------------ + +%------------------------------Initialization of ZIP load---------% khuang 8 JUL if ~isempty(zip) nZip=size(zip,1); zipIdx=zip(:,1); IiL=zeros(nZip,nlvl+1); BiL=zeros(nZip,nlvl+1); + % prepare the necessary matrix by blocks% khuang 8 JUL Bi0=abs(V0(zipIdx)); JI=zip(:,6); KI=-zip(:,9); - + % current % khuang 8 JUL Ii0L=Rzip0.*(JI+1j*KI).*V0(zipIdx)./Bi0; Ji0L=real(Ii0L); Ki0L=imag(Ii0L); IiL(:,1)=Ii0L; BiL(:,1)=Bi0; - + % voltage% khuang 8 JUL Ci0=real(V0(zipIdx)); Di0=imag(V0(zipIdx)); @@ -326,38 +412,43 @@ if ~isempty(zip) else IiL=zeros(0,nlvl+1); 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); if ~isempty(syn) - synIdx=syn(:,1); - wgb=syn(:,4); - modSyn=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); - - d=zeros(nSyn,nlvl+1); - w=zeros(nSyn,nlvl+1); - eq1=zeros(nSyn,nlvl+1); - eq2=zeros(nSyn,nlvl+1); - ed1=zeros(nSyn,nlvl+1); - ed2=zeros(nSyn,nlvl+1); - psiq=zeros(nSyn,nlvl+1); - psid=zeros(nSyn,nlvl+1); + synIdx =syn(:,1);% index number of Generators% khuang 8 JUL + wgb =syn(:,4);% maybe the base value% khuang 8 JUL + modSyn =syn(:,5);% the order of generator models% khuang 8 JUL + 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); + + d=zeros(nSyn,nlvl+1); % delta% khuang 8 JUL + w=zeros(nSyn,nlvl+1); % omega% khuang 8 JUL + eq1=zeros(nSyn,nlvl+1); %eq'% khuang 8 JUL + eq2=zeros(nSyn,nlvl+1); %eq''% khuang 8 JUL + ed1=zeros(nSyn,nlvl+1); %ed'% khuang 8 JUL + ed2=zeros(nSyn,nlvl+1); %ed''% khuang 8 JUL + psiq=zeros(nSyn,nlvl+1); % not sure, only in 8th order model% khuang 8 JUL + psid=zeros(nSyn,nlvl+1); % not sure, only in 8th order model% khuang 8 JUL JG=zeros(nSyn,nlvl+1); KG=zeros(nSyn,nlvl+1); IGq=zeros(nSyn,nlvl+1); @@ -373,7 +464,7 @@ if ~isempty(syn) sind=sin(d0); CG0=C0(synIdx); DG0=D0(synIdx); - + % the first value is given here, notice all are 8th order model% khuang 8 JUL d(:,1)=d0; w(:,1)=w0; eq1(:,1)=eq10; @@ -382,13 +473,19 @@ if ~isempty(syn) ed2(:,1)=ed20; psiq(:,1)=psiq0; psid(:,1)=psid0; + + % transform between grid side and dq side% khuang 8 JUL + VGd(:,1)=sind.*CG0-cosd.*DG0; VGq(:,1)=cosd.*CG0+sind.*DG0; - Cd(:,1)=cosd; - Sd(:,1)=sind; + % now they are under dq side% khuang 8 JUL + + 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; Pm(:,1)=Pm0; + %check if controller exists% khuang 8 JUL if ~isempty(Ef1) Ef(:,2)=Ef1; end @@ -399,26 +496,47 @@ if ~isempty(syn) Pm(:,2)=Pm1; 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 Mats=zeros(nSyn,4); MatsR=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]); + % 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 + % 8th order, no need to change% khuang 8 JUL if modSyn(i)==8 IGd(i,1)=(eq20(i)-psid0(i))/Xgd2(i); IGq(i,1)=(-ed20(i)-psiq0(i))/Xgq2(i); Mats(i,:)=[sind(i),cosd(i),-cosd(i),sind(i)]; + % 6th order% khuang 8 JUL 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)); 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)]; + % 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)]/... (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),... 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 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)); @@ -427,6 +545,7 @@ if ~isempty(syn) (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),... 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 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)); @@ -435,6 +554,7 @@ if ~isempty(syn) (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),... 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 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)); @@ -443,6 +563,7 @@ if ~isempty(syn) (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),... 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 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)); @@ -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)]; 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); + % 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);... sparse(synIdx,synIdx,MatsRs(:,3),nbus,nbus),sparse(synIdx,synIdx,MatsRs(:,4),nbus,nbus)]; else @@ -478,43 +602,57 @@ else Ef=zeros(0,nlvl+1); Pm=zeros(0,nlvl+1); end +%------------------------------Initialization of GEN------------------% khuang 8 JUL +%------------------------------EnD------------------------------------% khuang 8 JUL + + +%------------------------------Initialization of Exciter------------------% khuang 8 JUL +%------------------------------START------------------------------------% khuang 8 JUL if ~isempty(exc) - nExc=size(exc,1); - % All Type 3 AVR + nExc =size(exc,1); + % All Type 3 AVR, a 3rd order controller % for Type 3 AVR, avr0(:,1:3) are Vavrm, Vavrr, Vavrf, % and avr0(:,4) is reference Vref (input for secondary voltage control). - 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); - - Vavrm=zeros(nExc,nlvl+1); - Vavrr=zeros(nExc,nlvl+1); - Vavrf=zeros(nExc,nlvl+1); - Vavrref=zeros(nExc,nlvl+1); + 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); + + %here I need to check why Vavrref is time varing instead of constant% khuang 8 JUL + % memory is given to state variables of EXC% khuang 8 JUL + Vavrm = 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); Vavrr(:,1)=real(Vavrr0); Vavrf(:,1)=real(Vavrf0); Vavrref(:,1)=real(Vavrref0); + + % here Varref1 is given with syspara.% khuang 8 JUL if ~isempty(Varref1) Vavrref(:,2)=Varref1; end + % non-windup limiter, check the limit.% khuang 8 JUL tavrMaxDiff=Vavrf(:,1)-VavrMax; tavrMinDiff=Vavrf(:,1)-VavrMin; + % label values in different interval.% khuang 8 JUL avrSt=zeros(nExc,1); avrSt(tavrMaxDiff>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)=VavrMax(avrSt== 1); Ef(excIdx(avrSt== 0),1)=Vavrf(avrSt==0,1); @@ -525,37 +663,48 @@ else Vavrf=zeros(0,nlvl+1); Vavrref=zeros(0,nlvl+1); 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) - nTg=size(tg,1); + nTg = size(tg,1); % Type 2 Turbing governor. - tgIdx=tg(:,1); + % one DE, one AE and one limiter + tgIdx = tg(:,1); - wtgref=tg(:,3); - Rtg=tg(:,4); - Ttgmax=tg(:,5); - Ttgmin=tg(:,6); - Ttg2=tg(:,7); - Ttg1=tg(:,8); + wtgref = tg(:,3); + Rtg = tg(:,4); + Ttgmax = tg(:,5); + Ttgmin = tg(:,6); + Ttg2 = tg(:,7); + Ttg1 = tg(:,8); - tgovg=zeros(nTg,nlvl+1); - tgovm=zeros(nTg,nlvl+1); - Tmech=zeros(nTg,nlvl+1); + tgovg = zeros(nTg,nlvl+1); % tg % khuang 8 JUL + tgovm = zeros(nTg,nlvl+1); % Tmi* % khuang 8 JUL + Tmech = zeros(nTg,nlvl+1); % Tmi0 % khuang 8 JUL + % zeroth value is given % khuang 8 JUL tgovg(:,1)=real(tgovg0); tgovm(:,1)=real(tgovm0); Tmech(:,1)=real(tgovmech0); + if ~isempty(Tmech1) Tmech(:,2)=Tmech1; end + % check if limit is approached % khuang 8 JUL tgovMaxDiff=tgovm(:,1)-Ttgmax; tgovMinDiff=tgovm(:,1)-Ttgmin; govSt=zeros(nTg,1); govSt(tgovMaxDiff>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==1),1)=Ttgmax(govSt==1,1); Pm(tgIdx(govSt==-1),1)=Ttgmin(govSt==-1,1); @@ -564,7 +713,12 @@ else tgovm=zeros(0,nlvl+1); Tmech=zeros(0,nlvl+1); 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(:,1)=f0; synTag=zeros(nbus,1); @@ -581,6 +735,7 @@ for islIdx=1:nIslands end end +%AGC part % khuang 8 JUL if ~isempty(agc) agcExt=zeros(nbus,size(agc,2)); agcExt(agc(:,1),:)=agc; @@ -592,6 +747,7 @@ else fdk=zeros(nbus,1); end +% this is for long term dynamic, it seems that not considered here % khuang 8 JUL if ~isempty(cac)&&~isempty(cluster) else @@ -599,14 +755,16 @@ else vg=zeros(0,nlvl+1); end -% freq +% freq relevant part induced by AGC % khuang 8 JUL FreqReal=sparse(1:nbus,1:nbus,-freqKeptTag.*fdk.*E0,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); 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); +% Counting influence of ZIP load into Y matrix. % khuang 8 JUL if ~isempty(zip) Y11=Y11-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,1),[nbus,1]),nbus,nbus); Y12=Y12-sparse(1:nbus,1:nbus,accumarray(zipIdx,LHS_MatZip(:,2),[nbus,1]),nbus,nbus); @@ -615,12 +773,14 @@ if ~isempty(zip) end YLHS=[Y11,Y12;Y21,Y22]; +% Counting influence of Motors into small Y matrix % khuang 8 JUL if ~isempty(ind) 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(:,2),nbus,nbus),sparse(1:nbus,1:nbus,LHS_MatInd_Bus_sqz(:,4),nbus,nbus)]; end +% Counting influence of generators into small Y matrix % khuang 8 JUL if ~isempty(syn) YLHS=YLHS+MatGCD; end @@ -629,6 +789,7 @@ idxNonSw=find(busType~=2); idxNonSwxD=find(fswTagxD==0); 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]),... [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)];... @@ -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'); % 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'); if useLU @@ -660,8 +826,22 @@ if useLU 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 + % 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); seq2p=getseq(i+1,2); seq3=getseq(i,3); @@ -670,6 +850,15 @@ for i=1:nlvl idxSeq2p=sum(seq2p>=i,2); idxSeq3=sum(seq3==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,:); seq2x=seq2(idxSeq2x==0,:); seq2m=getseq(i-1,2); @@ -677,7 +866,10 @@ for i=1:nlvl RHSILr=zeros(nbus,1); RHSILi=zeros(nbus,1); + + % This part is for induction motor % khuang 8 JUL 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); % 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+... @@ -692,6 +884,8 @@ for i=1:nlvl % ./(2*Hm.*s(:,1)*i); % 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); if i>=2 s(:,i+1)=s(:,i+1)+... @@ -704,6 +898,7 @@ for i=1:nlvl if i==2 s(:,i+1)=s(:,i+1)+Rind1.*T0./(2*Hm*i); 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); % 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]; % end + % count the influence of dynamic of slip into rigt hand side + % vector.% khuang 8 JUL 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)]'; - %accumulate currents + %accumulate currents IL.% khuang 8 JUL RHSILr=accumarray(indIdx,rhsBus(1,:)',[nbus,1]); RHSILi=accumarray(indIdx,rhsBus(2,:)',[nbus,1]); @@ -734,6 +931,7 @@ for i=1:nlvl % RHSILi=accumarray(indIdx,rhsBus(4,:)',[nbus,1]); end + % strat update ZIP load into currents.% khuang 8 JUL RHSIiLr=zeros(nbus,1); RHSIiLi=zeros(nbus,1); if ~isempty(zip) @@ -745,6 +943,7 @@ for i=1:nlvl RHSIiLi=accumarray(zipIdx,RHSILi_full,[nbus,1]); end + % Start update generators.% khuang 8 JUL RHSIGr=zeros(nbus,1); RHSIGi=zeros(nbus,1); if ~isempty(syn) @@ -752,7 +951,7 @@ for i=1:nlvl RhsEq=zeros(nSyn,1); IGdAdd=zeros(nSyn,1); IGqAdd=zeros(nSyn,1); - + % select different models for generators.% khuang 8 JUL if modelTag(8)>0 d(modSyn==8,i+1)=(wgb(modSyn==8).*w(modSyn==8,i))/i; w(modSyn==8,i+1)=(Pm(modSyn==8,i)-... @@ -822,10 +1021,11 @@ for i=1:nlvl RhsEd(modSyn==2)=0; RhsEq(modSyn==2)=eq1(modSyn==2,i+1); end - + % this part may be different from DT.% khuang 8 JUL AG0=cosp(:,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 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); @@ -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); end - Cd(:,i+1)=AG0; - Sd(:,i+1)=BG0; + % now Sd store high order terms of sin(dta) and Cd for cos(dta), + % 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); - VGqCr=sum(Cd(:,seq2x(:,1)+1).*VGq(:,seq2x(:,2)+1),2); - VGdSr=sum(Sd(:,seq2x(:,1)+1).*VGd(:,seq2x(:,2)+1),2); - VGqSr=sum(Sd(:,seq2x(:,1)+1).*VGq(:,seq2x(:,2)+1),2); - JCr=sum(Cd(:,seq2x(:,1)+1).*JG(:,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);% Vq*cosdta% khuang 8 JUL + 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);% Vq*sindta% khuang 8 JUL + 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); JSr=sum(Sd(:,seq2x(:,1)+1).*JG(:,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); 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); + % current injections from generators IG.% khuang 8 JUL RHSIGr=accumarray(synIdx,RHSIGxr,[nbus,1]); RHSIGi=accumarray(synIdx,RHSIGxi,[nbus,1]); end - + % update exciter, 3 state variables.% khuang 8 JUL if ~isempty(exc) 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; @@ -871,6 +1077,7 @@ for i=1:nlvl Ef(excIdx(avrSt== 0),i+1)=Vavrf(avrSt==0,i+1); end + % update agc, one state variables.% khuang 8 JUL if ~isempty(agc) dpg(:,i+1)=-f(:,i).*agcExt(:,4)/i; for islIdx=1:nIslands @@ -882,6 +1089,7 @@ for i=1:nlvl end end % TODO: steady-state model + % update generator participation part from agc.% khuang 8 JUL if ~isempty(syn) %dynamic model (synchronous generators) if ~isempty(tg) 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)); end end - + % update Turbine, 2 state variables.% khuang 8 JUL if ~isempty(tg) 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); @@ -906,12 +1114,15 @@ for i=1:nlvl RHS2=-0.5*real(sum(V(:,seq2R(:,1)+1).*conj(V(:,seq2R(:,2)+1)),2)); RHS3=sum(-W(:,seq2R(:,1)+1).*V(:,seq2R(:,2)+1),2); + if i==1 RHS2=RHS2+0.5*VspSq2(:,2); end compactRHS1=RHS1(busType~=2); 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);... imag(compactRHS1)+RHSILi(busType~=2)+RHSIiLi(busType~=2)-RHSIGi(busType~=2);... RHS2(ipv);... @@ -919,7 +1130,9 @@ for i=1:nlvl imag(RHS3(busType~=2));... zeros(sum(freqKeptTagxRef),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 IS_OCTAVE x = real(MxQ * MxQx* (MxU \ (MxL \ (MxP * RHS)))) ; @@ -930,6 +1143,7 @@ for i=1:nlvl x=real(LHS_mat\RHS); end + % x= [V;W;Q_pv;f] xC=real(V(:,i+1)); xD=imag(V(:,i+1)); 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 + % now update the Algebric variables for motors:IL,IR,VM.% khuang 8 JUL if ~isempty(ind) % for j=1:nInd % 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; end + % now update the Algebric variables for ZIP loads.% khuang 8 JUL if ~isempty(zip) IiL(:,i+1)=(LHS_MatZip(:,1)+1j*LHS_MatZip(:,3)).*real(V(zipIdx,i+1))+(LHS_MatZip(:,2)+1j*LHS_MatZip(:,4)).*imag(V(zipIdx,i+1))+(RHSILr_full+1j*RHSILi_full); BiL(:,i+1)=Mat_BZip(:,1).*real(V(zipIdx,i+1))+Mat_BZip(:,2).*imag(V(zipIdx,i+1))+RHS_BZip; end + % now update the Algebric variables for Generators: Vd,Vq,Id,Iq.% khuang 8 JUL if ~isempty(syn) 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; @@ -976,6 +1193,7 @@ for i=1:nlvl end end +% Output value: coefficients for every order.% khuang 8 JUL Q=real(Q); s=real(s); d=real(d); diff --git a/internal/hemMachinePFmultiStage.m b/internal/hemMachinePFmultiStage.m index 58385a3..d411acd 100644 --- a/internal/hemMachinePFmultiStage.m +++ b/internal/hemMachinePFmultiStage.m @@ -92,8 +92,8 @@ while alphaConfirm<1-alphaTol/1000 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; - 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); + 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.iden=iden; if exist([iden,'.mat'],'file') if ~exist('p_amd','var') diff --git a/internal/hemMachinePFmultiStageDyn.m b/internal/hemMachinePFmultiStageDyn.m index 8e0c1d0..c6e198d 100644 --- a/internal/hemMachinePFmultiStageDyn.m +++ b/internal/hemMachinePFmultiStageDyn.m @@ -146,6 +146,7 @@ nZip=size(zip,1); % psiqc=zeros(nSyn,1);dpsiqc=zeros(nSyn,1); alphaList=0; +ratiob4 = NaN; diffList=0; diffMul=1.05; alphaPreMult=0.91; @@ -198,7 +199,7 @@ if ~isempty(varOpt) && isfield(varOpt,'useNomoveCtrl') else useNomoveCtrl=1; end - +alphab4 = NaN; while alphaConfirm0 +% 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; while 1 - + %pade 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); @@ -439,6 +488,7 @@ while alphaConfirm=maxCount;alphax=alphax/alphaPreMult;end alpha=alphax; @@ -606,7 +679,7 @@ while alphaConfirm0)=0; tgovmActive((tgovm0-Ttgmin)<0)=0; pActiveInSyn(tgIdx)=tgovmActive; +% else +% tgIdx=[]; +% tgovmActive=zeros(nTg,1); +% pActiveInSyn=zeros(nTg,1); + end + if ~isempty(exc) texcActive=ones(nExc,1); texcActive((Vavrf0-VavrMax)>0)=0; texcActive((Vavrf0-VavrMin)<0)=0; texcActiveInSyn=zeros(nSyn,1); - texcActiveInSyn(tgIdx)=texcActive; + texcActiveInSyn(tgIdx)=texcActive; + end numSynOnBus=accumarray(syn(:,1),1,[nbus,1]); %delta 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)],[],[],[],[]); diff=checkEquationBalanceSynDyn(SysDataxx,SysParaxx,xTemp,dxc); -end \ No newline at end of file +end +% % function [xTemp,diff,exitflag]=singleStepTRAP(SimData,SysData,SysPara,xt,t,dt,iden,AEmethod) +% % % Single step of trapezoidal formulation +% % % Copyright (C) Rui Yao +% % % +% % % 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]))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 \ No newline at end of file diff --git a/runDynamicSimulationExec.m b/runDynamicSimulationExec.m index 051735e..054d8c4 100644 --- a/runDynamicSimulationExec.m +++ b/runDynamicSimulationExec.m @@ -163,7 +163,7 @@ else end %% Temporary code here -% eventList(:,6)=0.0; +% eventList(:,6)=4.00; % eventList(:,7)=0.1; % Efstd=1.2; % evtDynZip(:,2)=0.2; @@ -194,8 +194,9 @@ if ~hotStart [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); +% [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; diff --git a/runPowerSAS.m b/runPowerSAS.m index 80a88ef..f7fc147 100644 --- a/runPowerSAS.m +++ b/runPowerSAS.m @@ -60,6 +60,7 @@ else addLog(['OUTPUT IS DISABLED. Check ',options.outputPath,'/~default_log.txt for logs.'],'DEBUG'); end addLog(['The session ID is ',timestamp,'.'],'INFO'); +%khuang 919 options.timestamp=timestamp; homePath=pwd; @@ -167,9 +168,9 @@ function options=regulateOptions(options) if nargin<1 options=[]; 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,'dbstop');options.dbstop=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