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

Merge pull request #5 from yangliu0101/doc

Browse Source 
Add basic version of EMT functionality
release
yangliu0101 2 years ago committed by GitHub
parent d1f3c13664 68c835ca1e
commit 59d62c45d2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
26 changed files with 1320 additions and 79 deletions
Show Stats Download Patch File Download Diff File

368
EMT/DTsolver.m
Unescape View File

@ -0,0 +1,368 @@
function [t_dtm,sol_dtm]=DTsolver(TimeRange,Resolution,x,N_DTM,Sys)
Gen_num=length(Sys.GenIdx);
%% Asign space to store the series
DTA=zeros(Gen_num,N_DTM+1); % machine number, term number, dimention of each term
THETA=DTA;OMG=DTA; LAMBDA_F=DTA;LAMBDA_D=DTA;LAMBDA_Q1=DTA;LAMBDA_Q2=DTA;
V1 =DTA; EFD=DTA; P1= DTA; P2= DTA; PM=DTA; W_R=DTA;
SIN=DTA; COS=DTA; ID =DTA; IQ =DTA;
LAMBDA_AD=DTA; LAMBDA_AQ=DTA; PE=DTA; VD_PP=DTA; VQ_PP=DTA;
VPP_A=DTA; VPP_B=DTA; VPP_C=DTA;
VT=DTA; VTP=DTA;
VBUS5 = zeros(3,N_DTM+1);
VBUS6 = VBUS5; VBUS7 =VBUS5; VBUS8 =VBUS5; VBUS9 =VBUS5; VBUS10 =VBUS5; VBUS11 =VBUS5;
ILINE1 = VBUS5; ILINE2 = VBUS5; ILINE3 = VBUS5; ILINE4 = VBUS5; ILINE5 = VBUS5; ILINE6 = VBUS5;
ILINE7 = VBUS5; ILINE8 = VBUS5; ILINE9 = VBUS5; ILINE10= VBUS5; ILINE11= VBUS5; ILINE12= VBUS5;
ILOAD1 = VBUS5; ILOAD2 = VBUS5;
%% Coefficients of t^0, initial value
DTA(:,1)=x(1:4,1); THETA(:,1)=x(5:8,1); OMG(:,1)=x(9:12,1);
LAMBDA_F(:,1)=x(13:16,1); LAMBDA_D(:,1)=x(17:20,1); LAMBDA_Q1(:,1)=x(21:24,1); LAMBDA_Q2(:,1)=x(25:28,1); V1(:,1)=x(29:32,1);
efd = x(33:36,1); efd = min(max(efd,Sys.Emin),Sys.Emax); EFD(:,1)= efd;
p1 = x(37:40,1); p1 = min(max(p1,Sys.GVmin),Sys.GVmax); P1(:,1) = p1;
P2(:,1)= x(41:44,1);
PM(:,1)= P2(:,1)-Sys.Gdt.*OMG(:,1);
W_R(:,1) =OMG(:,1)+1;
VBUS5(:,1) =x(45:47,1); VBUS6(:,1) =x(48:50,1); VBUS7(:,1) =x(51:53,1); VBUS8(:,1) =x(54:56,1);
VBUS9(:,1) =x(57:59,1); VBUS10(:,1)=x(60:62,1); VBUS11(:,1)=x(63:65,1);
ILINE1(:,1)=x(66:68,1); ILINE2(:,1)=x(69:71,1); ILINE3(:,1)=x(72:74,1); ILINE4(:,1)=x(75:77,1);
ILINE5(:,1)=x(78:80,1); ILINE6(:,1)=x(81:83,1); ILINE7(:,1)=x(84:86,1); ILINE8(:,1)=x(87:89,1);
ILINE9(:,1)=x(90:92,1); ILINE10(:,1)=x(93:95,1); ILINE11(:,1)=x(96:98,1); ILINE12(:,1)=x(99:101,1);
ILOAD1(:,1)=x(102:104,1); ILOAD2(:,1)=x(105:107,1);
VTP(:,1) =x(108:111,1);
% some other terms
SIN(:,1) =sin(THETA(:,1)); % sin(theta) no need to add '
COS(:,1) =cos(THETA(:,1)); % cos(theta)
ID(:,1) = 2/3*( [ILINE1(1,1); ILINE2(1,1); ILINE3(1,1); ILINE4(1,1)].*COS(:,1) + [ILINE1(2,1); ILINE2(2,1); ILINE3(2,1); ILINE4(2,1)].*( COS(:,1)*cos(2*pi/3)+SIN(:,1)*sin(2*pi/3) ) + [ILINE1(3,1); ILINE2(3,1); ILINE3(3,1); ILINE4(3,1)].*( COS(:,1)*cos(2*pi/3)-SIN(:,1)*sin(2*pi/3) ) );
IQ(:,1) = 2/3*( -[ILINE1(1,1); ILINE2(1,1); ILINE3(1,1); ILINE4(1,1)].*SIN(:,1) - [ILINE1(2,1); ILINE2(2,1); ILINE3(2,1); ILINE4(2,1)].*( SIN(:,1)*cos(2*pi/3)-COS(:,1)*sin(2*pi/3) ) - [ILINE1(3,1); ILINE2(3,1); ILINE3(3,1); ILINE4(3,1)].*( SIN(:,1)*cos(2*pi/3)+COS(:,1)*sin(2*pi/3) ) );
LAMBDA_AD(:,1) = Sys.Lad_pp.*(-ID(:,1)+LAMBDA_F(:,1)./Sys.Lfd+LAMBDA_D(:,1)./Sys.L1d);
LAMBDA_AQ(:,1) = Sys.Laq_pp.*(-IQ(:,1)+LAMBDA_Q1(:,1)./Sys.L1q+LAMBDA_Q2(:,1)./Sys.L2q);
PE(:,1) = LAMBDA_AD(:,1).*IQ(:,1)-LAMBDA_AQ(:,1).*ID(:,1);
VD_PP(:,1)=Sys.Vd_pp_coeff1.*ID(:,1)+Sys.Vd_pp_coeff2.*LAMBDA_F(:,1)+Sys.Vd_pp_coeff3.*LAMBDA_D(:,1)...
-W_R(:,1).*Sys.Laq_pp.*( LAMBDA_Q1(:,1)./Sys.L1q + LAMBDA_Q2(:,1)./Sys.L2q ) + Sys.Lad_pp./Sys.Lfd.*EFD(:,1) ;
VQ_PP(:,1)=Sys.Vq_pp_coeff1.*IQ(:,1)+Sys.Vq_pp_coeff2.*LAMBDA_Q1(:,1)+Sys.Vq_pp_coeff3.*LAMBDA_Q2(:,1)...
+W_R(:,1).*Sys.Lad_pp.*( LAMBDA_F(:,1)./Sys.Lfd + LAMBDA_D(:,1)./Sys.L1d );
VPP_A(:,1) = COS(:,1).*VD_PP(:,1) -SIN(:,1).*VQ_PP(:,1);
VPP_B(:,1) = ( COS(:,1)*cos(2*pi/3)+SIN(:,1)*sin(2*pi/3) ).*VD_PP(:,1) -( SIN(:,1)*cos(2*pi/3)-COS(:,1)*sin(2*pi/3) ).*VQ_PP(:,1);
VPP_C(:,1) = ( COS(:,1)*cos(2*pi/3)-SIN(:,1)*sin(2*pi/3) ).*VD_PP(:,1) -( SIN(:,1)*cos(2*pi/3)+COS(:,1)*sin(2*pi/3) ).*VQ_PP(:,1);
VT2(:,1) = abs(VD_PP(:,1)+1i*VQ_PP(:,1)).^2;
VT(:,1) = abs(VD_PP(:,1)+1i*VQ_PP(:,1));
%% preparation
SimuTime = TimeRange(2)-TimeRange(1);
t_dtm=(0:Resolution:SimuTime)'+TimeRange(1);
% hh = TimeRange(2)-t_dtm(end);
% if abs(hh)>1e-5
% t_dtm = [t_dtm;TimeRange(2)];
% end
%% store the results at each time step
% assign space
res_dta = zeros(Gen_num,size(t_dtm,1));
res_theta = res_dta;
res_omg = res_dta;
res_lambda_f = res_dta;
res_lambda_d = res_dta;
res_lambda_q1 = res_dta;
res_lambda_q2 = res_dta;
res_v1 = res_dta;
res_efd = res_dta;
res_p1 = res_dta;
res_p2 = res_dta;
res_vbus5 = zeros(3,size(t_dtm,1)); res_vbus6 = res_vbus5; res_vbus7 = res_vbus5;
res_vbus8 = res_vbus5; res_vbus9 = res_vbus5; res_vbus10 = res_vbus5; res_vbus11 = res_vbus5;
res_iline1 = res_vbus5; res_iline2 = res_vbus5; res_iline3 = res_vbus5; res_iline4 = res_vbus5;
res_iline5 = res_vbus5; res_iline6 = res_vbus5; res_iline7 = res_vbus5; res_iline8 = res_vbus5;
res_iline9 = res_vbus5; res_iline10= res_vbus5; res_iline11= res_vbus5; res_iline12= res_vbus5;
res_iload1 = res_vbus5; res_iload2 = res_vbus5;
res_vtp = res_dta;
% store the first one
res_dta(:,1) =DTA(:,1);
res_theta(:,1)=THETA(:,1);
res_omg(:,1) =OMG(:,1,1);
res_lambda_f(:,1) =LAMBDA_F(:,1);
res_lambda_d(:,1) =LAMBDA_D(:,1);
res_lambda_q1(:,1) =LAMBDA_Q1(:,1);
res_lambda_q2(:,1) =LAMBDA_Q2(:,1);
res_v1(:,1) =V1(:,1);
res_efd(:,1) =EFD(:,1);
res_p1(:,1) =P1(:,1);
res_p2(:,1) =P2(:,1);
res_vbus5(:,1) = VBUS5(:,1); res_vbus6(:,1)= VBUS6(:,1); res_vbus7(:,1)= VBUS7(:,1); res_vbus8(:,1)= VBUS8(:,1);
res_vbus9(:,1) = VBUS9(:,1); res_vbus10(:,1)= VBUS10(:,1); res_vbus11(:,1) = VBUS11(:,1);
res_iline1(:,1) = ILINE1(:,1); res_iline2(:,1) = ILINE2(:,1); res_iline3(:,1) = ILINE3(:,1); res_iline4(:,1) = ILINE4(:,1);
res_iline5(:,1) = ILINE5(:,1); res_iline6(:,1) = ILINE6(:,1); res_iline7(:,1) = ILINE7(:,1); res_iline8(:,1) = ILINE8(:,1);
res_iline9(:,1) = ILINE9(:,1); res_iline10(:,1)= ILINE10(:,1); res_iline11(:,1)= ILINE11(:,1); res_iline12(:,1)= ILINE12(:,1);
res_iload1(:,1) = ILOAD1(:,1); res_iload2(:,1) = ILOAD2(:,1);
res_vtp(:,1) = VTP(:,1);
%% begin the for loop to calculate
Tvec=zeros(1,N_DTM+1);
for Numstep=1:size(t_dtm,1)-1 % loop of time step
hhh = t_dtm(Numstep+1)-t_dtm(Numstep);
Tvec(1:N_DTM+1)=hhh.^(0:N_DTM);
idx = zeros(N_DTM+1,1); idx(1) =1; % for constant terms
% below are used to check limits
def_dt = ( V1(:,1).*Sys.Ek-EFD(:,1).*Sys.Lad./Sys.Rfd )./Sys.Ete.*Sys.Rfd./Sys.Lad;
dP1_dt = ( (Sys.Pref-OMG(:,1))./Sys.Gr-P1(:,1) )./Sys.Gt1;
for k = 1:N_DTM % loop of order
% synchronous generator differential equations
DTA(:,k+1) = 1/k*2*pi*60*OMG(:,k);
THETA(:,k+1) = 1/k*2*pi*60*(OMG(:,k)+1*idx(k));
OMG(:,k+1) = 1/k*Sys.w0*(PM(:,k) - PE(:,k) - Sys.D.*OMG(:,k)./Sys.w0)./Sys.H/2;
LAMBDA_F(:,k+1) = 1/k*120*pi*(EFD(:,k)-Sys.Rfd./Sys.Lfd.*(LAMBDA_F(:,k)-LAMBDA_AD(:,k)));
LAMBDA_D(:,k+1) = 1/k*(-120)*pi*Sys.R1d./Sys.L1d.*(LAMBDA_D(:,k) -LAMBDA_AD(:,k));
LAMBDA_Q1(:,k+1) = 1/k*(-120)*pi*Sys.R1q./Sys.L1q.*(LAMBDA_Q1(:,k)-LAMBDA_AQ(:,k));
LAMBDA_Q2(:,k+1) = 1/k*(-120)*pi*Sys.R2q./Sys.L2q.*(LAMBDA_Q2(:,k)-LAMBDA_AQ(:,k));
V1(:,k+1) = 1/k*( -VTP(:,k)+Sys.Vref*idx(k)-V1(:,k)-Sys.Eta.*(VT(:,k)-VTP(:,k))./Sys.T_vt )./Sys.Etb;
EFD_high = 1/k*( V1(:,k).*Sys.Ek-EFD(:,k).*Sys.Lad./Sys.Rfd )./Sys.Ete.*Sys.Rfd./Sys.Lad;
EFD_high(def_dt(EFD(:,1)>=Sys.Emax)>0)=0;
EFD_high(def_dt(EFD(:,1)<=Sys.Emin)<0)=0;
EFD(:,k+1) = EFD_high;
P1_high = 1/k*( (Sys.Pref*idx(k)-OMG(:,k))./Sys.Gr-P1(:,k) )./Sys.Gt1;
P1_high(dP1_dt(P1(:,1)>=Sys.GVmax)>0)=0;
P1_high(dP1_dt(P1(:,1)<=Sys.GVmin)<0)=0;
P1(:,k+1) = P1_high;
% P2(:,k+1)= 1/k*(P1(:,k)-P2(:,k)+P1(:,k+1).*Sys.Gt2)./Sys.Gt3;
P2(:,k+1)= 1/k*(P1(:,k)-P2(:,k)+k*P1(:,k+1).*Sys.Gt2)./Sys.Gt3;
% network differential equations
VBUS5(:,k+1) = 1/k*Sys.Cnode_3Phase(:,:,5)*(ILINE1(:,k)-ILINE5(:,k));
VBUS6(:,k+1) = 1/k*Sys.Cnode_3Phase(:,:,6)*(ILINE2(:,k)+ILINE5(:,k)-ILINE6(:,k));
VBUS7(:,k+1) = 1/k*Sys.Cnode_3Phase(:,:,7)*(ILINE6(:,k)-ILINE7(:,k)-ILINE8(:,k)-ILOAD1(:,k));
VBUS8(:,k+1) = 1/k*Sys.Cnode_3Phase(:,:,8)*(ILINE7(:,k)+ILINE8(:,k)-ILINE9(:,k)-ILINE10(:,k));
VBUS9(:,k+1) = 1/k*Sys.Cnode_3Phase(:,:,9)*(ILINE9(:,k)+ILINE10(:,k)-ILINE11(:,k)-ILOAD2(:,k));
VBUS10(:,k+1)= 1/k*Sys.Cnode_3Phase(:,:,10)*(ILINE4(:,k)+ILINE11(:,k)-ILINE12(:,k));
VBUS11(:,k+1)= 1/k*Sys.Cnode_3Phase(:,:,11)*(ILINE3(:,k)+ILINE12(:,k));
ILINE1(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,1)*( [VPP_A(1,k); VPP_B(1,k); VPP_C(1,k)]- VBUS5(:,k) -Sys.Gen_R(:,:,1)*ILINE1(:,k) ); % IBR voltage
ILINE2(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,2)*( [VPP_A(2,k); VPP_B(2,k); VPP_C(2,k)]- VBUS6(:,k) -Sys.Gen_R(:,:,2)*ILINE2(:,k) );
ILINE3(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,3)*( [VPP_A(3,k); VPP_B(3,k); VPP_C(3,k)]- VBUS11(:,k) -Sys.Gen_R(:,:,3)*ILINE3(:,k) );
ILINE4(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,4)*( [VPP_A(4,k); VPP_B(4,k); VPP_C(4,k)]- VBUS10(:,k) -Sys.Gen_R(:,:,4)*ILINE4(:,k) );
ILINE5(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,5)*( VBUS5(:,k)-VBUS6(:,k)-Sys.Rline_3Phase(:,:,5) * ILINE5(:,k) );
ILINE6(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,6)*( VBUS6(:,k)-VBUS7(:,k)-Sys.Rline_3Phase(:,:,6) * ILINE6(:,k) );
ILINE7(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,7)*( VBUS7(:,k)-VBUS8(:,k)-Sys.Rline_3Phase(:,:,7) * ILINE7(:,k) );
ILINE8(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,8)*( VBUS7(:,k)-VBUS8(:,k)-Sys.Rline_3Phase(:,:,8) * ILINE8(:,k) );
ILINE9(:,k+1)= 1/k*Sys.Lline_3Phase(:,:,9)*( VBUS8(:,k)-VBUS9(:,k)-Sys.Rline_3Phase(:,:,9) * ILINE9(:,k) );
ILINE10(:,k+1)=1/k*Sys.Lline_3Phase(:,:,10)*( VBUS8(:,k)-VBUS9(:,k)-Sys.Rline_3Phase(:,:,10) * ILINE10(:,k) );
ILINE11(:,k+1)=1/k*Sys.Lline_3Phase(:,:,11)*( VBUS9(:,k)-VBUS10(:,k)-Sys.Rline_3Phase(:,:,11) * ILINE11(:,k) );
ILINE12(:,k+1)=1/k*Sys.Lline_3Phase(:,:,12)*( VBUS10(:,k)-VBUS11(:,k)-Sys.Rline_3Phase(:,:,12) * ILINE12(:,k) );
ILOAD1(:,k+1)=1/k*Sys.Lload_3Phase(:,:,1)*( VBUS7(:,k)-Sys.Rload_3Phase(:,:,1) * ILOAD1(:,k) );
ILOAD2(:,k+1)=1/k*Sys.Lload_3Phase(:,:,2)*( VBUS9(:,k)-Sys.Rload_3Phase(:,:,2) * ILOAD2(:,k) );
VTP(:,k+1)=1/k*(VT(:,k)-VTP(:,k))./Sys.T_vt;
% below is algebric part
PM(:,k+1)= P2(:,k+1)-Sys.Gdt.*OMG(:,k+1);
W_R(:,k+1) = OMG(:,k+1)+1*idx(k+1);
COS(:,k+1)= -1/(k)*sum(repmat([k:-1:1],Gen_num,1).*SIN(:,1:k).*THETA(:,k+1:-1:2),2);
SIN(:,k+1)= 1/(k)*sum(repmat([k:-1:1],Gen_num,1).*COS(:,1:k).*THETA(:,k+1:-1:2),2);
ID(:,k+1) = 2/3*sum(( [ILINE1(1,1:k+1); ILINE2(1,1:k+1); ILINE3(1,1:k+1); ILINE4(1,1:k+1)].*COS(:,k+1:-1:1) + [ILINE1(2,1:k+1); ILINE2(2,1:k+1); ILINE3(2,1:k+1); ILINE4(2,1:k+1)].*( COS(:,k+1:-1:1)*cos(2*pi/3)+SIN(:,k+1:-1:1)*sin(2*pi/3) ) + [ILINE1(3,1:k+1); ILINE2(3,1:k+1); ILINE3(3,1:k+1); ILINE4(3,1:k+1)].*( COS(:,k+1:-1:1)*cos(2*pi/3)-SIN(:,k+1:-1:1)*sin(2*pi/3) ) ),2);
IQ(:,k+1) = 2/3*sum((-[ILINE1(1,1:k+1); ILINE2(1,1:k+1); ILINE3(1,1:k+1); ILINE4(1,1:k+1)].*SIN(:,k+1:-1:1) - [ILINE1(2,1:k+1); ILINE2(2,1:k+1); ILINE3(2,1:k+1); ILINE4(2,1:k+1)].*( SIN(:,k+1:-1:1)*cos(2*pi/3)-COS(:,k+1:-1:1)*sin(2*pi/3) ) - [ILINE1(3,1:k+1); ILINE2(3,1:k+1); ILINE3(3,1:k+1); ILINE4(3,1:k+1)].*( SIN(:,k+1:-1:1)*cos(2*pi/3)+COS(:,k+1:-1:1)*sin(2*pi/3) ) ),2);
LAMBDA_AD(:,k+1) = Sys.Lad_pp.*(-ID(:,k+1)+LAMBDA_F(:,k+1)./Sys.Lfd+LAMBDA_D(:,k+1)./Sys.L1d);
LAMBDA_AQ(:,k+1) = Sys.Laq_pp.*(-IQ(:,k+1)+LAMBDA_Q1(:,k+1)./Sys.L1q+LAMBDA_Q2(:,k+1)./Sys.L2q);
PE(:,k+1) = sum(LAMBDA_AD(:,1:k+1).*IQ(:,k+1:-1:1)-LAMBDA_AQ(:,1:k+1).*ID(:,k+1:-1:1),2);
VD_PP(:,k+1)=Sys.Vd_pp_coeff1.*ID(:,k+1)+Sys.Vd_pp_coeff2.*LAMBDA_F(:,k+1)+Sys.Vd_pp_coeff3.*LAMBDA_D(:,k+1)...
-sum(W_R(:,1:k+1).*Sys.Laq_pp.*( LAMBDA_Q1(:,k+1:-1:1)./Sys.L1q + LAMBDA_Q2(:,k+1:-1:1)./Sys.L2q ),2) + Sys.Lad_pp./Sys.Lfd.*EFD(:,k+1) ;
VQ_PP(:,k+1)=Sys.Vq_pp_coeff1.*IQ(:,k+1)+Sys.Vq_pp_coeff2.*LAMBDA_Q1(:,k+1)+Sys.Vq_pp_coeff3.*LAMBDA_Q2(:,k+1)...
+sum(W_R(:,1:k+1).*Sys.Lad_pp.*( LAMBDA_F(:,k+1:-1:1)./Sys.Lfd + LAMBDA_D(:,k+1:-1:1)./Sys.L1d ),2);
VPP_A(:,k+1) = sum(COS(:,1:k+1).*VD_PP(:,k+1:-1:1) -SIN(:,1:k+1).*VQ_PP(:,k+1:-1:1),2);
VPP_B(:,k+1) = sum(( COS(:,1:k+1)*cos(2*pi/3)+SIN(:,1:k+1)*sin(2*pi/3) ).*VD_PP(:,k+1:-1:1) -( SIN(:,1:k+1)*cos(2*pi/3)-COS(:,1:k+1)*sin(2*pi/3) ).*VQ_PP(:,k+1:-1:1),2);
VPP_C(:,k+1) = sum(( COS(:,1:k+1)*cos(2*pi/3)-SIN(:,1:k+1)*sin(2*pi/3) ).*VD_PP(:,k+1:-1:1) -( SIN(:,1:k+1)*cos(2*pi/3)+COS(:,1:k+1)*sin(2*pi/3) ).*VQ_PP(:,k+1:-1:1),2);
VT2(:,k+1) = sum(VD_PP(:,1:k+1).*VD_PP(:,k+1:-1:1)+ VQ_PP(:,1:k+1).*VQ_PP(:,k+1:-1:1),2);
VT(:,k+1) = (VT2(:,k+1) - sum(VT(:,2:k).*VT(:,k:-1:2),2))./(2*VT(:,1));
end
%% Evaluation, contain only differential equations
dta_2=sum(DTA.*repmat(Tvec,[Gen_num,1]),2);
theta_2=sum(THETA.*repmat(Tvec,[Gen_num,1]),2);
omg_2=sum(OMG.*repmat(Tvec,[Gen_num,1]),2);
lambda_f_2=sum(LAMBDA_F.*repmat(Tvec,[Gen_num,1]),2);
lambda_d_2=sum(LAMBDA_D.*repmat(Tvec,[Gen_num,1]),2);
lambda_q1_2=sum(LAMBDA_Q1.*repmat(Tvec,[Gen_num,1]),2);
lambda_q2_2=sum(LAMBDA_Q2.*repmat(Tvec,[Gen_num,1]),2);
vtp_2=sum(VTP.*repmat(Tvec,[Gen_num,1]),2);
v1_2=sum(V1.*repmat(Tvec,[Gen_num,1]),2);
efd_2=sum(EFD.*repmat(Tvec,[Gen_num,1]),2);
p1_2=sum(P1.*repmat(Tvec,[Gen_num,1]),2);
p2_2=sum(P2.*repmat(Tvec,[Gen_num,1]),2);
vbus5_2=sum(VBUS5.*repmat(Tvec,[3,1]),2);
vbus6_2=sum(VBUS6.*repmat(Tvec,[3,1]),2);
vbus7_2=sum(VBUS7.*repmat(Tvec,[3,1]),2);
vbus8_2=sum(VBUS8.*repmat(Tvec,[3,1]),2);
vbus9_2=sum(VBUS9.*repmat(Tvec,[3,1]),2);
vbus10_2=sum(VBUS10.*repmat(Tvec,[3,1]),2);
vbus11_2=sum(VBUS11.*repmat(Tvec,[3,1]),2);
iline1_2=sum(ILINE1.*repmat(Tvec,[3,1]),2);
iline2_2=sum(ILINE2.*repmat(Tvec,[3,1]),2);
iline3_2=sum(ILINE3.*repmat(Tvec,[3,1]),2);
iline4_2=sum(ILINE4.*repmat(Tvec,[3,1]),2);
iline5_2=sum(ILINE5.*repmat(Tvec,[3,1]),2);
iline6_2=sum(ILINE6.*repmat(Tvec,[3,1]),2);
iline7_2=sum(ILINE7.*repmat(Tvec,[3,1]),2);
iline8_2=sum(ILINE8.*repmat(Tvec,[3,1]),2);
iline9_2=sum(ILINE9.*repmat(Tvec,[3,1]),2);
iline10_2=sum(ILINE10.*repmat(Tvec,[3,1]),2);
iline11_2=sum(ILINE11.*repmat(Tvec,[3,1]),2);
iline12_2=sum(ILINE12.*repmat(Tvec,[3,1]),2);
iload1_2=sum(ILOAD1.*repmat(Tvec,[3,1]),2);
iload2_2=sum(ILOAD2.*repmat(Tvec,[3,1]),2);
% Save the results
res_dta(:,Numstep+1)=dta_2;
res_theta(:,Numstep+1)=theta_2;
res_omg(:,Numstep+1)=omg_2;
res_lambda_f(:,Numstep+1)=lambda_f_2;
res_lambda_d(:,Numstep+1)=lambda_d_2;
res_lambda_q1(:,Numstep+1)=lambda_q1_2;
res_lambda_q2(:,Numstep+1)=lambda_q2_2;
res_v1(:,Numstep+1)=v1_2;
res_efd(:,Numstep+1)=efd_2;
res_p1(:,Numstep+1)=p1_2;
res_p2(:,Numstep+1)=p2_2;
res_vtp(:,Numstep+1)=vtp_2;
res_vbus5(:,Numstep+1)=vbus5_2;
res_vbus6(:,Numstep+1)=vbus6_2;
res_vbus7(:,Numstep+1)=vbus7_2;
res_vbus8(:,Numstep+1)=vbus8_2;
res_vbus9(:,Numstep+1)=vbus9_2;
res_vbus10(:,Numstep+1)=vbus10_2;
res_vbus11(:,Numstep+1)=vbus11_2;
res_iline1(:,Numstep+1)=iline1_2;
res_iline2(:,Numstep+1)=iline2_2;
res_iline3(:,Numstep+1)=iline3_2;
res_iline4(:,Numstep+1)=iline4_2;
res_iline5(:,Numstep+1)=iline5_2;
res_iline6(:,Numstep+1)=iline6_2;
res_iline7(:,Numstep+1)=iline7_2;
res_iline8(:,Numstep+1)=iline8_2;
res_iline9(:,Numstep+1)=iline9_2;
res_iline10(:,Numstep+1)=iline10_2;
res_iline11(:,Numstep+1)=iline11_2;
res_iline12(:,Numstep+1)=iline12_2;
res_iload1(:,Numstep+1)=iload1_2;
res_iload2(:,Numstep+1)=iload2_2;
%% Update Initial value
% differential part
DTA(:,1)=dta_2;
THETA(:,1)=theta_2;
OMG(:,1)=omg_2;
LAMBDA_F(:,1)=lambda_f_2;
LAMBDA_D(:,1)=lambda_d_2;
LAMBDA_Q1(:,1)=lambda_q1_2;
LAMBDA_Q2(:,1)=lambda_q2_2;
V1(:,1)=v1_2;
EFD(:,1)=efd_2;
P1(:,1)=p1_2;
P2(:,1)=p2_2;
VTP(:,1)=vtp_2;
VBUS5(:,1)=vbus5_2;
VBUS6(:,1)=vbus6_2;
VBUS7(:,1)=vbus7_2;
VBUS8(:,1)=vbus8_2;
VBUS9(:,1)=vbus9_2;
VBUS10(:,1)=vbus10_2;
VBUS11(:,1)=vbus11_2;
ILINE1(:,1)=iline1_2;
ILINE2(:,1)=iline2_2;
ILINE3(:,1)=iline3_2;
ILINE4(:,1)=iline4_2;
ILINE5(:,1)=iline5_2;
ILINE6(:,1)=iline6_2;
ILINE7(:,1)=iline7_2;
ILINE8(:,1)=iline8_2;
ILINE9(:,1)=iline9_2;
ILINE10(:,1)=iline10_2;
ILINE11(:,1)=iline11_2;
ILINE12(:,1)=iline12_2;
ILOAD1(:,1)=iload1_2;
ILOAD2(:,1)=iload2_2;
% algebric part
PM(:,1)= P2(:,1)-Sys.Gdt.*OMG(:,1);
W_R(:,1) = OMG(:,1)+1;
SIN(:,1) =sin(THETA(:,1)); % sin(theta) no need to add '
COS(:,1) =cos(THETA(:,1)); % cos(theta)
ID(:,1) = 2/3*( [ILINE1(1,1); ILINE2(1,1); ILINE3(1,1); ILINE4(1,1)].*COS(:,1) + [ILINE1(2,1); ILINE2(2,1); ILINE3(2,1); ILINE4(2,1)].*( COS(:,1)*cos(2*pi/3)+SIN(:,1)*sin(2*pi/3) ) + [ILINE1(3,1); ILINE2(3,1); ILINE3(3,1); ILINE4(3,1)].*( COS(:,1)*cos(2*pi/3)-SIN(:,1)*sin(2*pi/3) ) );
IQ(:,1) = 2/3*( -[ILINE1(1,1); ILINE2(1,1); ILINE3(1,1); ILINE4(1,1)].*SIN(:,1) - [ILINE1(2,1); ILINE2(2,1); ILINE3(2,1); ILINE4(2,1)].*( SIN(:,1)*cos(2*pi/3)-COS(:,1)*sin(2*pi/3) ) - [ILINE1(3,1); ILINE2(3,1); ILINE3(3,1); ILINE4(3,1)].*( SIN(:,1)*cos(2*pi/3)+COS(:,1)*sin(2*pi/3) ) );
LAMBDA_AD(:,1) = Sys.Lad_pp.*(-ID(:,1)+LAMBDA_F(:,1)./Sys.Lfd+LAMBDA_D(:,1)./Sys.L1d);
LAMBDA_AQ(:,1) = Sys.Laq_pp.*(-IQ(:,1)+LAMBDA_Q1(:,1)./Sys.L1q+LAMBDA_Q2(:,1)./Sys.L2q);
PE(:,1) = LAMBDA_AD(:,1).*IQ(:,1)-LAMBDA_AQ(:,1).*ID(:,1);
VD_PP(:,1)=Sys.Vd_pp_coeff1.*ID(:,1)+Sys.Vd_pp_coeff2.*LAMBDA_F(:,1)+Sys.Vd_pp_coeff3.*LAMBDA_D(:,1)...
-W_R(:,1).*Sys.Laq_pp.*( LAMBDA_Q1(:,1)./Sys.L1q + LAMBDA_Q2(:,1)./Sys.L2q ) + Sys.Lad_pp./Sys.Lfd.*EFD(:,1) ;
VQ_PP(:,1)=Sys.Vq_pp_coeff1.*IQ(:,1)+Sys.Vq_pp_coeff2.*LAMBDA_Q1(:,1)+Sys.Vq_pp_coeff3.*LAMBDA_Q2(:,1)...
+W_R(:,1).*Sys.Lad_pp.*( LAMBDA_F(:,1)./Sys.Lfd + LAMBDA_D(:,1)./Sys.L1d );
VPP_A(:,1) = COS(:,1).*VD_PP(:,1) -SIN(:,1).*VQ_PP(:,1);
VPP_B(:,1) = ( COS(:,1)*cos(2*pi/3)+SIN(:,1)*sin(2*pi/3) ).*VD_PP(:,1) -( SIN(:,1)*cos(2*pi/3)-COS(:,1)*sin(2*pi/3) ).*VQ_PP(:,1);
VPP_C(:,1) = ( COS(:,1)*cos(2*pi/3)-SIN(:,1)*sin(2*pi/3) ).*VD_PP(:,1) -( SIN(:,1)*cos(2*pi/3)+COS(:,1)*sin(2*pi/3) ).*VQ_PP(:,1);
VT2(:,1) = abs(VD_PP(:,1)+1i*VQ_PP(:,1)).^2;
VT(:,1) = abs(VD_PP(:,1)+1i*VQ_PP(:,1));
end
sol_dtm = zeros(110,size(t_dtm,1)); % generator state variables
sol_dtm(1:4,:) = res_dta;
sol_dtm(5:8,:) = res_theta;
sol_dtm(9:12,:) = res_omg;
sol_dtm(13:16,:) = res_lambda_f;
sol_dtm(17:20,:) = res_lambda_d;
sol_dtm(21:24,:) = res_lambda_q1;
sol_dtm(25:28,:) = res_lambda_q2;
sol_dtm(29:32,:) = res_v1;
sol_dtm(33:36,:) = res_efd;
sol_dtm(37:40,:) = res_p1;
sol_dtm(41:44,:) = res_p2;
sol_dtm(45:47,:) = res_vbus5;
sol_dtm(48:50,:) = res_vbus6;
sol_dtm(51:53,:) = res_vbus7;
sol_dtm(54:56,:) = res_vbus8;
sol_dtm(57:59,:) = res_vbus9;
sol_dtm(60:62,:) = res_vbus10;
sol_dtm(63:65,:) = res_vbus11;
sol_dtm(66:68,:) = res_iline1;
sol_dtm(69:71,:) = res_iline2;
sol_dtm(72:74,:) = res_iline3;
sol_dtm(75:77,:) = res_iline4;
sol_dtm(78:80,:) = res_iline5;
sol_dtm(81:83,:) = res_iline6;
sol_dtm(84:86,:) = res_iline7;
sol_dtm(87:89,:) = res_iline8;
sol_dtm(90:92,:) = res_iline9;
sol_dtm(93:95,:) = res_iline10;
sol_dtm(96:98,:) = res_iline11;
sol_dtm(99:101,:) = res_iline12;
sol_dtm(102:104,:) = res_iload1;
sol_dtm(105:107,:) = res_iload2;
sol_dtm(108:111,:) = res_vtp;

54
EMT/calEMT.m
Unescape View File

@ -0,0 +1,54 @@
function [res] = calEMT(testSystem)
%%
[orderOfDT, timeStepOfDT] = setSASParameters();
[Base, w, FtBus, time_step, time_step_on, TimeRange1, FtTime, TotalTime, TimeRange2, TimeRange3] = set_simulation_parameters();
%%
[R_line_positive, wL_line_positive, B_line_positive, int_t, LC_fre, t_base] = setDefaultParameters();
[mpc, SUCESS] = runpf(testSystem);
Sys = myFunction(mpc);
[Sys,Laa0_pp,Lab0_pp,Laa2_pp] = transformParameters(Sys,t_base);
[Vt,dta,Id, Iq, Edp, Eqp, Edpp, Eqpp, Pm, Efd, Sys, omg, V1, efd, P1, P2, Epp] = generator_states(mpc, Sys, Base);
[Lambda_F, Lambda_D, Lambda_Q1, Lambda_Q2] = calculateLambda(Sys, Id, Iq, efd);
[R_line, L_line, C_line] = tline_to_matrix(R_line_positive, wL_line_positive, B_line_positive,LC_fre);
[L_tf_positive L_tf LoadPQ Zload R_load L_load Sys LoadShunt YShunt C_Shunt CShunt] = fun1(Sys,LC_fre,mpc,Vt,Base);
[Vbus_3Phase branch_number c Sys Iline_3Phase ind_line1 ind_line2 Zbranch Ibranch ...
Amp Ang Iload_3Phase ind_load theta] = fun2(Sys, mpc, dta, int_t, w, Zload);
L_temp = calculate_L_temp(L_tf, L_line, Laa0_pp, Lab0_pp);
[Sys] = calculateSystemMatrices(mpc, Sys, CShunt, R_line, wL_line_positive, branch_number, L_temp, C_line);
[Vbus, Iline, Iload] = convert_to_single_vector(Vbus_3Phase, Iline_3Phase, Iload_3Phase, Sys, branch_number);
Sys = update_system_parameters(Sys, Laa0_pp, Lab0_pp);
x_ini=[dta; theta; omg; Lambda_F; Lambda_D; Lambda_Q1; Lambda_Q2; V1; efd; P1; P2;...
Vbus(3*length(Sys.GenIdx)+1:end); Iline; Iload; abs(Epp)];
%%
options = odeset('RelTol',1e-4,'AbsTol',1e-4); % limit the error
[t_pre,x_pre]= ode23t(@funEMT_pre,TimeRange1(1):time_step:TimeRange1(2),x_ini,options,Sys);
x_ini_on = x_pre(end,:)';
[t_on,x_on] = ode23t(@funEMT_on,TimeRange2(1):time_step_on:TimeRange2(2),x_ini_on,options,Sys);
x_ini_post = x_on(end,:)';
[t_dt,x_dt]=DTsolver(TimeRange3,timeStepOfDT,x_ini_post,orderOfDT,Sys); %x_ini_post
t_post = t_dt;
x_post = x_dt';
res.t_final = [t_pre; t_on(2:end); t_post(2:end)];
res.sol_final = [x_pre; x_on(2:end,:); x_post(2:end,:)];
return

7
EMT/calculateLambda.m
Unescape View File

@ -0,0 +1,7 @@
function [Lambda_F, Lambda_D, Lambda_Q1, Lambda_Q2] = calculateLambda(Sys, Id, Iq, efd)
Ifd=efd./Sys.Rfd;
Lambda_F = -Sys.Lad.*Id+(Sys.Lfd+Sys.Lad).*Ifd; % steady state Id1=0
Lambda_D = -Sys.Lad.*Id+Sys.Lad.*Ifd;
Lambda_Q1 = -Sys.Laq.*Iq;
Lambda_Q2 = -Sys.Laq.*Iq;

68
EMT/calculateSystemMatrices.m
Unescape View File

@ -0,0 +1,68 @@
function [Sys] = calculateSystemMatrices(mpc, Sys, CShunt, R_line, wL_line_positive, branch_number, L_temp, C_line);
% Calculate System Matrices
% Inputs:
% - mpc: power system data
% - Sys: structure containing system matrices
% - CShunt: shunt capacitance matrix
% - L_tf: transformer leakage inductance matrix
% - R_line: transmission line resistance
% - L_line: transmission line inductance
% - wL_line_positive: positive sequence angular frequency
% - Laa0_pp: stator inductance matrix element Laa0, prime-prime frame
% - Lab0_pp: stator inductance matrix element Lab0, prime-prime frame
% Output:
% - Sys: structure containing calculated system matrices
L_tf = L_temp.L_tf;
L_line = L_temp.L_line;
Laa0_pp = L_temp.Laa0_pp;
Lab0_pp = L_temp.Lab0_pp;
Sys.Cnode_3Phase=zeros(3,3,Sys.bus_number);
% T line
for k=1:1:branch_number
ind_line1=mpc.branch(k,1);
ind_line2=mpc.branch(k,2);
% if it is transmission line branch
if ismember(ind_line1,Sys.NonGenIdx)
term=mpc.branch(k,4)/wL_line_positive*C_line;
Sys.Cnode_3Phase(:,:,ind_line1) = Sys.Cnode_3Phase(:,:,ind_line1) + term;
Sys.Cnode_3Phase(:,:,ind_line2) = Sys.Cnode_3Phase(:,:,ind_line2) + term;
end
end
% Fix shunt
for k=1:1:length(Sys.ShuntIdx)
Sys.Cnode_3Phase(:,:,Sys.ShuntIdx(k))= Sys.Cnode_3Phase(:,:,Sys.ShuntIdx(k))+ CShunt(:,:,k);
end
% transformed into inverse matrices
for k=1:1:Sys.bus_number
Sys.Cnode_3Phase(:,:,k)= inv(Sys.Cnode_3Phase(:,:,k));
end
% all the Branch R L, consider transformer, T line, and load
Sys.Rline_3Phase=zeros(3,3,branch_number);
Sys.Lline_3Phase=zeros(3,3,branch_number);
L_stator=zeros(3,3,length(Sys.GenIdx));
% T line
for k=1:1:branch_number
ind_line1=mpc.branch(k,1);
ind_line2=mpc.branch(k,2);
if ismember(ind_line1,Sys.GenIdx)
L_stator(:,:,k)= [Laa0_pp(k) Lab0_pp(k) Lab0_pp(k); Lab0_pp(k) Laa0_pp(k) Lab0_pp(k); Lab0_pp(k) Lab0_pp(k) Laa0_pp(k)]/120/pi;
Sys.Rline_3Phase(:,:,k) = zeros(3,3);
Sys.Lline_3Phase(:,:,k) = inv(L_tf+L_stator(:,:,k));
% if it is transmission line branch
elseif ismember(ind_line1,Sys.NonGenIdx)
Sys.Rline_3Phase(:,:,k) = mpc.branch(k,4)/wL_line_positive*R_line;
Sys.Lline_3Phase(:,:,k) = inv(mpc.branch(k,4)/wL_line_positive*L_line);
end
end
for k=1:1:length(Sys.GenIdx)
Sys.Gen_R(:,:,k)= [Sys.Ra(k), 0, 0; 0, Sys.Ra(k), 0; 0, 0, Sys.Ra(k)];
L_stator(:,:,k)= [Laa0_pp(k) Lab0_pp(k) Lab0_pp(k); Lab0_pp(k) Laa0_pp(k) Lab0_pp(k); Lab0_pp(k) Lab0_pp(k) Laa0_pp(k)]/120/pi;
% because Sys.Laa0_pp is already updated
Sys.Lline_3Phase(:,:,k) = Sys.Lline_3Phase(:,:,k)+ L_stator(:,:,k);
end

6
EMT/calculate_L_temp.m
Unescape View File

@ -0,0 +1,6 @@
function L_temp = calculate_L_temp(L_tf, L_line, Laa0_pp, Lab0_pp)
L_temp.L_tf = L_tf;
L_temp.L_line = L_line;
L_temp.Laa0_pp = Laa0_pp;
L_temp.Lab0_pp = Lab0_pp;
end

26
EMT/convert_to_single_vector.m
Unescape View File

@ -0,0 +1,26 @@
function [Vbus, Iline, Iload] = convert_to_single_vector(Vbus_3Phase, Iline_3Phase, Iload_3Phase, Sys, branch_number)
% This function converts the three-phase voltage, branch current, and load current
% vectors into a single vector format for each of them.
Vbus=zeros(3*Sys.bus_number,1);
for k=1:1:Sys.bus_number
for k1= 1:1:3
Vbus((k-1)*3+k1)=Vbus_3Phase(k1,1,k);
end
end
% change three phase branch current to a single vector
Iline=zeros(3*branch_number,1);
for k=1:1:branch_number
for k1= 1:1:3
Iline((k-1)*3+k1)=Iline_3Phase(k1,1,k);
end
end
% change three phase load current to a single vector
Iload=zeros(3*length(Sys.LoadIdx),1);
for k=1:1:length(Sys.LoadIdx)
for k1= 1:1:3
Iload((k-1)*3+k1)=Iload_3Phase(k1,1,k);
end
end

27
EMT/fun1.m
Unescape View File

@ -0,0 +1,27 @@
function [L_tf_positive L_tf LoadPQ Zload R_load L_load Sys LoadShunt YShunt C_Shunt CShunt] = fun1(Sys,LC_fre,mpc,Vt,Base)
%% transformer parameters, per unit
% Assume all the transformers have the same impedance
L_tf_positive=mpc.branch(Sys.GenIdx(1),4)/(LC_fre); %positive sequence leakage reactance
L_tf=[L_tf_positive 0 0; 0 L_tf_positive 0; 0 0 L_tf_positive];
%% Constant PQ load parameter, per unit
% PQ load
LoadPQ = (mpc.bus(Sys.LoadIdx,3)+1i*mpc.bus(Sys.LoadIdx,4))./Base;
Zload = conj((abs(Vt(Sys.LoadIdx,:)).^2)./LoadPQ); % constant Z load
R_load = real(Zload);
L_load = imag(Zload)/(LC_fre);
Sys.Rload_3Phase = zeros(3,3,length(R_load));
Sys.Lload_3Phase = Sys.Rload_3Phase;
for k=1:1:size(R_load)
Sys.Rload_3Phase(:,:,k)=[R_load(k) 0 0; 0 R_load(k) 0; 0 0 R_load(k)];
Sys.Lload_3Phase(:,:,k)=inv([L_load(k) 0 0; 0 L_load(k) 0; 0 0 L_load(k)]);
end
% fixed shunt
LoadShunt = (mpc.bus(Sys.ShuntIdx,5)-1i*mpc.bus(Sys.ShuntIdx,6))./Base; % fixed shunt load, use - because it generate reactive power
YShunt = conj(LoadShunt./1.^2 ); % constant Z load, do not use abs(Vt(LoadIdx,:)).^2, because voltage should be 1
C_Shunt = imag(YShunt)/(LC_fre);
CShunt = zeros(3,3,length(C_Shunt));
for k=1:1:length(C_Shunt)
CShunt(:,:,k)=[C_Shunt(k) 0 0; 0 C_Shunt(k) 0; 0 0 C_Shunt(k)];
end

42
EMT/fun2.m
Unescape View File

@ -0,0 +1,42 @@
function [Vbus_3Phase branch_number c Sys Iline_3Phase ind_line1 ind_line2 Zbranch Ibranch ...
Amp Ang Iload_3Phase ind_load theta] = fun2(Sys, mpc, dta, int_t, w, Zload)
%% current and voltage, three phase
% all the bus voltage, three phase, assume t=ini_t, and Va=sin(wt+angle)
Vbus_3Phase=zeros(3,1,Sys.bus_number);
for k=1:1:Sys.bus_number
Vbus_3Phase(:,:,k) = [mpc.bus(k,8).*sin(w*int_t+ mpc.bus(k,9)*pi/180);...
mpc.bus(k,8).*sin(w*int_t+(mpc.bus(k,9)-120)*pi/180);...
mpc.bus(k,8).*sin(w*int_t+(mpc.bus(k,9)+120)*pi/180)];
end
% all the branch circuit current, three phase, assume t=0, and Ia=sin(wt+angle)
[branch_number,c] = size(mpc.branch);
Sys.branch_number=branch_number;
Iline_3Phase=zeros(3,1,branch_number);
for k=1:1:branch_number % consider all branch and load
ind_line1=mpc.branch(k,1);
ind_line2=mpc.branch(k,2);
Zbranch = mpc.branch(k,3)+1i*mpc.branch(k,4);
Ibranch= ( mpc.bus(ind_line1,8)*exp(1i*mpc.bus(ind_line1,9)*pi/180)- ...
mpc.bus(ind_line2,8)*exp(1i*mpc.bus(ind_line2,9)*pi/180) )/Zbranch;
Amp=abs(Ibranch);
Ang=angle(Ibranch);
Iline_3Phase(:,:,k)=[Amp.*sin(w*int_t+Ang);...
Amp.*sin(w*int_t+Ang-120*pi/180);...
Amp.*sin(w*int_t+Ang+120*pi/180)];
end
% all load current
Iload_3Phase=zeros(3,1,length(Sys.LoadIdx));
for k=1:1:length(Sys.LoadIdx) % consider all branch and load
ind_load=Sys.LoadIdx(k);
Iload= mpc.bus(ind_load,8)*exp(1i*mpc.bus(ind_load,9)*pi/180)/Zload(k);
Amp=abs(Iload);
Ang=angle(Iload);
Iload_3Phase(:,:,k)=[Amp.*sin(w*int_t+Ang);...
Amp.*sin(w*int_t+Ang-120*pi/180);...
Amp.*sin(w*int_t+Ang+120*pi/180)];
end
%% Initialization of Theta angle
theta=dta-pi+120*pi*int_t; % I found that if int_t=0, theta=dta-pi

33
EMT/fun3.m
Unescape View File

@ -0,0 +1,33 @@
function [Sys, Laa0_pp, Lab0_pp, Laa2_pp, t_base] = fun3(Sys)
%% old, good version
t_base=1/(120*pi);
Sys.Lad=Sys.Xd-Sys.Xl;
Sys.Lfd=1./(1./(Sys.Xdp-Sys.Xl)-1./Sys.Lad);
Sys.Rfd=(Sys.Lad+Sys.Lfd)./(Sys.Td0p./t_base);
Sys.L1d=1./(1./(Sys.Xdpp-Sys.Xl)-1./Sys.Lfd-1./Sys.Lad);
Sys.R1d=(1./(1./Sys.Lad+1./Sys.Lfd)+Sys.L1d)./(Sys.Td0pp./t_base);
Sys.Laq=Sys.Xq-Sys.Xl;
Sys.L1q=1./(1./(Sys.Xqp-Sys.Xl)-1./Sys.Laq);
Sys.R1q=(Sys.Laq+Sys.L1q)./(Sys.Tq0p./t_base);
Sys.L2q=1./(1./(Sys.Xqpp-Sys.Xl)-1./Sys.L1q-1./Sys.Laq);
Sys.R2q=(1./(1./Sys.Laq+1./Sys.L1q)+Sys.L2q)./(Sys.Tq0pp./t_base);
Sys.Lad_pp=Sys.Xdpp-Sys.Xl;
Sys.Laq_pp=Sys.Xdpp-Sys.Xl;
Laa0_pp= Sys.Xl+( Sys.Lad_pp+ Sys.Laq_pp+Sys.L0-Sys.Xl )/3;
Lab0_pp= ( 2*Sys.L0-2*Sys.Xl-Sys.Lad_pp-Sys.Laq_pp )/6;
Laa2_pp= ( Sys.Lad_pp-Sys.Laq_pp )/3;
Sys.Vd_pp_coeff1=-(Sys.Rfd./(Sys.Lfd).^2+Sys.R1d./(Sys.L1d).^2).*(Sys.Lad_pp).^2;
Sys.Vd_pp_coeff2=-(Sys.Lad_pp.*Sys.Rfd./(Sys.Lfd).^2.*(1-Sys.Lad_pp./Sys.Lfd) - Sys.Lad_pp.^2.*Sys.R1d./(Sys.L1d).^2./Sys.Lfd);
Sys.Vd_pp_coeff3=(Sys.Lad_pp.^2.*Sys.Rfd./(Sys.Lfd).^2./Sys.L1d - Sys.Lad_pp.*Sys.R1d./(Sys.L1d).^2.*(1-Sys.Lad_pp./Sys.L1d));
Sys.Vq_pp_coeff1=-(Sys.R1q./(Sys.L1q).^2+Sys.R2q./(Sys.L2q).^2).*(Sys.Laq_pp).^2;
Sys.Vq_pp_coeff2=-(Sys.Laq_pp.*Sys.R1q./(Sys.L1q).^2.*(1-Sys.Laq_pp./Sys.L1q) - Sys.Laq_pp.^2.*Sys.R2q./(Sys.L2q).^2./Sys.L1q);
Sys.Vq_pp_coeff3=(Sys.Laq_pp.^2.*Sys.R1q./(Sys.L1q).^2./Sys.L2q - Sys.Laq_pp.*Sys.R2q./(Sys.L2q).^2.*(1-Sys.Laq_pp./Sys.L2q));

99
EMT/funEMT_on.m
Unescape View File

@ -0,0 +1,99 @@
function [xp] = funEMT_on(t,x,Sys)
xp = zeros(size(x));
R_ground = eye(3)*100000;
% global GLO
dta = x(1:4,1);
theta = x(4+1:2*4,1);
omg = x(2*4+1:3*4,1);
Lambda_F = x(3*4+1:4*4,1);
Lambda_D = x(4*4+1:5*4,1);
Lambda_Q1 = x(5*4+1:6*4,1);
Lambda_Q2 = x(6*4+1:7*4,1);
V1 = x(7*4+1:8*4,1);
efd = x(8*4+1:9*4,1);
% efd = min(max(efd,Sys.Emin),Sys.Emax);
P1 = x(9*4+1:10*4,1);
% P1 = min(max(P1,Sys.GVmin),Sys.GVmax);
P2 = x(10*4+1:11*4,1);
Pm = P2-Sys.Gdt.*omg;
Vbus = x(11*4+1:11*4+3*(11-4),1);
Iline = x(11*4+3*(11-4)+1:11*4+3*(11-4)+3*12,1);
Iload = x(11*4+3*(11-4)+3*12+1:11*4+3*(11-4)+3*12+3*2,1);
Vtp = x(107+1:107+4,1);
w_r = omg+1;
COS = cos(theta);
SIN = sin(theta);
Id = 2/3*( Iline(1:3:3*4).*COS + Iline(2:3:3*4).*( COS*cos(2*pi/3)+SIN*sin(2*pi/3) ) + Iline(3:3:3*4).*( COS*cos(2*pi/3)-SIN*sin(2*pi/3) ) );
Iq = 2/3*( -Iline(1:3:3*4).*SIN - Iline(2:3:3*4).*( SIN*cos(2*pi/3)-COS*sin(2*pi/3) ) - Iline(3:3:3*4).*( SIN*cos(2*pi/3)+COS*sin(2*pi/3) ) );
Lambda_ad = Sys.Lad_pp.*(-Id+Lambda_F./Sys.Lfd+Lambda_D./Sys.L1d);
Lambda_aq = Sys.Laq_pp.*(-Iq+Lambda_Q1./Sys.L1q+Lambda_Q2./Sys.L2q);
Pe = (Lambda_ad.*Iq-Lambda_aq.*Id);
Vd_pp=Sys.Vd_pp_coeff1.*Id+Sys.Vd_pp_coeff2.*Lambda_F+Sys.Vd_pp_coeff3.*Lambda_D...
-w_r.*Sys.Laq_pp.*( Lambda_Q1./Sys.L1q + Lambda_Q2./Sys.L2q ) + Sys.Lad_pp./Sys.Lfd.*efd ;
Vq_pp=Sys.Vq_pp_coeff1.*Iq+Sys.Vq_pp_coeff2.*Lambda_Q1+Sys.Vq_pp_coeff3.*Lambda_Q2...
+w_r.*Sys.Lad_pp.*( Lambda_F./Sys.Lfd + Lambda_D./Sys.L1d );
% Inv_Park(:,:,k)*[0; Vd_pp(k); Vq_pp(k)]
Vpp_a = COS.*Vd_pp -SIN.*Vq_pp;
Vpp_b = ( COS*cos(2*pi/3)+SIN*sin(2*pi/3) ).*Vd_pp -( SIN*cos(2*pi/3)-COS*sin(2*pi/3) ).*Vq_pp;
Vpp_c = ( COS*cos(2*pi/3)-SIN*sin(2*pi/3) ).*Vd_pp -( SIN*cos(2*pi/3)+COS*sin(2*pi/3) ).*Vq_pp;
Vt = sqrt(Vd_pp.^2+Vq_pp.^2); % in fact it is internal voltage
% ODEs
xp(1:4,1) = 2*pi*60*omg;
xp(4+1 : 2*4,1) = 2*pi*60*(omg+1);
xp(2*4+1: 3*4,1) = Sys.w0*(Pm - Pe - Sys.D.*omg./Sys.w0)./Sys.H/2;
xp(3*4+1: 4*4,1) = 120*pi*(efd-Sys.Rfd./Sys.Lfd.*(Lambda_F-Lambda_ad));
xp(4*4+1: 5*4,1) = -120*pi*Sys.R1d./Sys.L1d.*(Lambda_D-Lambda_ad);
xp(5*4+1: 6*4,1) = -120*pi*Sys.R1q./Sys.L1q.*(Lambda_Q1-Lambda_aq);
xp(6*4+1: 7*4,1) = -120*pi*Sys.R2q./Sys.L2q.*(Lambda_Q2-Lambda_aq);
xp(7*4+1: 8*4,1) = ( -Vtp+Sys.Vref-V1-Sys.Eta.*(Vt-Vtp)./Sys.T_vt )./Sys.Etb; % substitute d(Vt')=(Vt-Vtp)./Sys.T_vt
xp(8*4+1: 9*4,1) = ( V1.*Sys.Ek-efd.*Sys.Lad./Sys.Rfd )./Sys.Ete.*Sys.Rfd./Sys.Lad;
xp(9*4+1: 10*4,1) = ( (Sys.Pref-omg)./Sys.Gr-P1 )./Sys.Gt1;
xp(10*4+1: 11*4,1) = (P1-P2+( (Sys.Pref-omg)./Sys.Gr-P1 )./Sys.Gt1.*Sys.Gt2)./Sys.Gt3;
xp(44+1:44+3,1) = Sys.Cnode_3Phase(:,:,5)*([sum(Iline(([1 ]-1)*3+1,1));sum(Iline(([1 ]-1)*3+2,1));sum(Iline(([1 ]-1)*3+3,1))]-[sum(Iline(([5 ]-1)*3+1,1));sum(Iline(([5 ]-1)*3+2,1));sum(Iline(([5 ]-1)*3+3,1))]);
xp(47+1:47+3,1) = Sys.Cnode_3Phase(:,:,6)*([sum(Iline(([2 5 ]-1)*3+1,1));sum(Iline(([2 5 ]-1)*3+2,1));sum(Iline(([2 5 ]-1)*3+3,1))]-[sum(Iline(([6 ]-1)*3+1,1));sum(Iline(([6 ]-1)*3+2,1));sum(Iline(([6 ]-1)*3+3,1))]);
xp(50+1:50+3,1) = Sys.Cnode_3Phase(:,:,7)*([sum(Iline(([6 ]-1)*3+1,1));sum(Iline(([6 ]-1)*3+2,1));sum(Iline(([6 ]-1)*3+3,1))]-[sum(Iline(([7 8 ]-1)*3+1,1));sum(Iline(([7 8 ]-1)*3+2,1));sum(Iline(([7 8 ]-1)*3+3,1))]-[sum(Iload(([1 ]-1)*3+1,1));sum(Iload(([1 ]-1)*3+2,1));sum(Iload(([1 ]-1)*3+3,1))]-R_ground*Vbus((7-5)*3+1:(7-5)*3+3));
xp(53+1:53+3,1) = Sys.Cnode_3Phase(:,:,8)*([sum(Iline(([7 8 ]-1)*3+1,1));sum(Iline(([7 8 ]-1)*3+2,1));sum(Iline(([7 8 ]-1)*3+3,1))]-[sum(Iline(([9 10 ]-1)*3+1,1));sum(Iline(([9 10 ]-1)*3+2,1));sum(Iline(([9 10 ]-1)*3+3,1))]);
xp(56+1:56+3,1) = Sys.Cnode_3Phase(:,:,9)*([sum(Iline(([9 10 ]-1)*3+1,1));sum(Iline(([9 10 ]-1)*3+2,1));sum(Iline(([9 10 ]-1)*3+3,1))]-[sum(Iline(([11 ]-1)*3+1,1));sum(Iline(([11 ]-1)*3+2,1));sum(Iline(([11 ]-1)*3+3,1))]-[sum(Iload(([2 ]-1)*3+1,1));sum(Iload(([2 ]-1)*3+2,1));sum(Iload(([2 ]-1)*3+3,1))]);
xp(59+1:59+3,1) = Sys.Cnode_3Phase(:,:,10)*([sum(Iline(([4 11 ]-1)*3+1,1));sum(Iline(([4 11 ]-1)*3+2,1));sum(Iline(([4 11 ]-1)*3+3,1))]-[sum(Iline(([12 ]-1)*3+1,1));sum(Iline(([12 ]-1)*3+2,1));sum(Iline(([12 ]-1)*3+3,1))]);
xp(62+1:62+3,1) = Sys.Cnode_3Phase(:,:,11)*([sum(Iline(([3 12 ]-1)*3+1,1));sum(Iline(([3 12 ]-1)*3+2,1));sum(Iline(([3 12 ]-1)*3+3,1))]);
xp(65+1:65+3,1)= Sys.Lline_3Phase(:,:,1)*( [Vpp_a(1); Vpp_b(1); Vpp_c(1)]- Vbus((5-5)*3+1:(5-5)*3+3) -Sys.Gen_R(:,:,1)*Iline((1-1)*3+1:(1-1)*3+3)); % IBR voltage
xp(68+1:68+3,1)= Sys.Lline_3Phase(:,:,2)*( [Vpp_a(2); Vpp_b(2); Vpp_c(2)]- Vbus((6-5)*3+1:(6-5)*3+3) -Sys.Gen_R(:,:,2)*Iline((2-1)*3+1:(2-1)*3+3));
xp(71+1:71+3,1)= Sys.Lline_3Phase(:,:,3)*( [Vpp_a(3); Vpp_b(3); Vpp_c(3)]- Vbus((11-5)*3+1:(11-5)*3+3)-Sys.Gen_R(:,:,3)*Iline((3-1)*3+1:(3-1)*3+3));
xp(74+1:74+3,1)= Sys.Lline_3Phase(:,:,4)*( [Vpp_a(4); Vpp_b(4); Vpp_c(4)]- Vbus((10-5)*3+1:(10-5)*3+3)-Sys.Gen_R(:,:,4)*Iline((4-1)*3+1:(4-1)*3+3));
xp(77+1:77+3,1)=Sys.Lline_3Phase(:,:,5)*( Vbus((5-5)*3+1:(5-5)*3+3)-Vbus((6-5)*3+1:(6-5)*3+3)-Sys.Rline_3Phase(:,:,5) * Iline((5-1)*3+1:(5-1)*3+3) );
xp(80+1:80+3,1)=Sys.Lline_3Phase(:,:,6)*( Vbus((6-5)*3+1:(6-5)*3+3)-Vbus((7-5)*3+1:(7-5)*3+3)-Sys.Rline_3Phase(:,:,6) * Iline((6-1)*3+1:(6-1)*3+3) );
xp(83+1:83+3,1)=Sys.Lline_3Phase(:,:,7)*( Vbus((7-5)*3+1:(7-5)*3+3)-Vbus((8-5)*3+1:(8-5)*3+3)-Sys.Rline_3Phase(:,:,7) * Iline((7-1)*3+1:(7-1)*3+3) );
xp(86+1:86+3,1)=Sys.Lline_3Phase(:,:,8)*( Vbus((7-5)*3+1:(7-5)*3+3)-Vbus((8-5)*3+1:(8-5)*3+3)-Sys.Rline_3Phase(:,:,8) * Iline((8-1)*3+1:(8-1)*3+3) );
xp(89+1:89+3,1)=Sys.Lline_3Phase(:,:,9)*( Vbus((8-5)*3+1:(8-5)*3+3)-Vbus((9-5)*3+1:(9-5)*3+3)-Sys.Rline_3Phase(:,:,9) * Iline((9-1)*3+1:(9-1)*3+3) );
xp(92+1:92+3,1)=Sys.Lline_3Phase(:,:,10)*( Vbus((8-5)*3+1:(8-5)*3+3)-Vbus((9-5)*3+1:(9-5)*3+3)-Sys.Rline_3Phase(:,:,10) * Iline((10-1)*3+1:(10-1)*3+3) );
xp(95+1:95+3,1)=Sys.Lline_3Phase(:,:,11)*( Vbus((9-5)*3+1:(9-5)*3+3)-Vbus((10-5)*3+1:(10-5)*3+3)-Sys.Rline_3Phase(:,:,11) * Iline((11-1)*3+1:(11-1)*3+3) );
xp(98+1:98+3,1)=Sys.Lline_3Phase(:,:,12)*( Vbus((10-5)*3+1:(10-5)*3+3)-Vbus((11-5)*3+1:(11-5)*3+3)-Sys.Rline_3Phase(:,:,12) * Iline((12-1)*3+1:(12-1)*3+3) );
xp(101+1:101+3,1)=Sys.Lload_3Phase(:,:,1)*( Vbus((7-5)*3+1:(7-5)*3+3)-Sys.Rload_3Phase(:,:,1) * Iload((1-1)*3+1:(1-1)*3+3) );
xp(104+1:104+3,1)=Sys.Lload_3Phase(:,:,2)*( Vbus((9-5)*3+1:(9-5)*3+3)-Sys.Rload_3Phase(:,:,2) * Iload((2-1)*3+1:(2-1)*3+3) );
xp(107+1:107+4,1)=(Vt-Vtp)./Sys.T_vt;
% GLO.V1 =[GLO.V1, Vt(1)];
% GLO.V2 =[GLO.V2, Vt(2)];
% GLO.V3 =[GLO.V3, Vt(3)];
% GLO.V4 =[GLO.V4, Vt(4)];
end

97
EMT/funEMT_pre.m
Unescape View File

@ -0,0 +1,97 @@
function [xp] = funEMT_pre(t,x,Sys)
xp = zeros(size(x));
% global GLO
dta = x(1:4,1);
theta = x(4+1:2*4,1);
omg = x(2*4+1:3*4,1);
Lambda_F = x(3*4+1:4*4,1);
Lambda_D = x(4*4+1:5*4,1);
Lambda_Q1 = x(5*4+1:6*4,1);
Lambda_Q2 = x(6*4+1:7*4,1);
V1 = x(7*4+1:8*4,1);
efd = x(8*4+1:9*4,1);
% efd = min(max(efd,Sys.Emin),Sys.Emax);
P1 = x(9*4+1:10*4,1);
% P1 = min(max(P1,Sys.GVmin),Sys.GVmax);
P2 = x(10*4+1:11*4,1);
Pm = P2-Sys.Gdt.*omg;
Vbus = x(11*4+1:11*4+3*(11-4),1);
Iline = x(11*4+3*(11-4)+1:11*4+3*(11-4)+3*12,1);
Iload = x(11*4+3*(11-4)+3*12+1:11*4+3*(11-4)+3*12+3*2,1);
Vtp = x(107+1:107+4,1);
w_r = omg+1;
COS = cos(theta);
SIN = sin(theta);
Id = 2/3*( Iline(1:3:3*4).*COS + Iline(2:3:3*4).*( COS*cos(2*pi/3)+SIN*sin(2*pi/3) ) + Iline(3:3:3*4).*( COS*cos(2*pi/3)-SIN*sin(2*pi/3) ) );
Iq = 2/3*( -Iline(1:3:3*4).*SIN - Iline(2:3:3*4).*( SIN*cos(2*pi/3)-COS*sin(2*pi/3) ) - Iline(3:3:3*4).*( SIN*cos(2*pi/3)+COS*sin(2*pi/3) ) );
Lambda_ad = Sys.Lad_pp.*(-Id+Lambda_F./Sys.Lfd+Lambda_D./Sys.L1d);
Lambda_aq = Sys.Laq_pp.*(-Iq+Lambda_Q1./Sys.L1q+Lambda_Q2./Sys.L2q);
Pe = (Lambda_ad.*Iq-Lambda_aq.*Id);
Vd_pp=Sys.Vd_pp_coeff1.*Id+Sys.Vd_pp_coeff2.*Lambda_F+Sys.Vd_pp_coeff3.*Lambda_D...
-w_r.*Sys.Laq_pp.*( Lambda_Q1./Sys.L1q + Lambda_Q2./Sys.L2q ) + Sys.Lad_pp./Sys.Lfd.*efd ;
Vq_pp=Sys.Vq_pp_coeff1.*Iq+Sys.Vq_pp_coeff2.*Lambda_Q1+Sys.Vq_pp_coeff3.*Lambda_Q2...
+w_r.*Sys.Lad_pp.*( Lambda_F./Sys.Lfd + Lambda_D./Sys.L1d );
% Inv_Park(:,:,k)*[0; Vd_pp(k); Vq_pp(k)]
Vpp_a = COS.*Vd_pp -SIN.*Vq_pp;
Vpp_b = ( COS*cos(2*pi/3)+SIN*sin(2*pi/3) ).*Vd_pp -( SIN*cos(2*pi/3)-COS*sin(2*pi/3) ).*Vq_pp;
Vpp_c = ( COS*cos(2*pi/3)-SIN*sin(2*pi/3) ).*Vd_pp -( SIN*cos(2*pi/3)+COS*sin(2*pi/3) ).*Vq_pp;
Vt = sqrt(Vd_pp.^2+Vq_pp.^2); % in fact it is internal voltage
% ODEs
xp(1:4,1) = 2*pi*60*omg;
xp(4+1 : 2*4,1) = 2*pi*60*(omg+1);
xp(2*4+1: 3*4,1) = Sys.w0*(Pm - Pe - Sys.D.*omg./Sys.w0)./Sys.H/2;
xp(3*4+1: 4*4,1) = 120*pi*(efd-Sys.Rfd./Sys.Lfd.*(Lambda_F-Lambda_ad));
xp(4*4+1: 5*4,1) = -120*pi*Sys.R1d./Sys.L1d.*(Lambda_D-Lambda_ad);
xp(5*4+1: 6*4,1) = -120*pi*Sys.R1q./Sys.L1q.*(Lambda_Q1-Lambda_aq);
xp(6*4+1: 7*4,1) = -120*pi*Sys.R2q./Sys.L2q.*(Lambda_Q2-Lambda_aq);
xp(7*4+1: 8*4,1) = ( -Vtp+Sys.Vref-V1-Sys.Eta.*(Vt-Vtp)./Sys.T_vt )./Sys.Etb; % substitute d(Vt')=(Vt-Vtp)./Sys.T_vt
xp(8*4+1: 9*4,1) = ( V1.*Sys.Ek-efd.*Sys.Lad./Sys.Rfd )./Sys.Ete.*Sys.Rfd./Sys.Lad;
xp(9*4+1: 10*4,1) = ( (Sys.Pref-omg)./Sys.Gr-P1 )./Sys.Gt1;
xp(10*4+1: 11*4,1) = (P1-P2+( (Sys.Pref-omg)./Sys.Gr-P1 )./Sys.Gt1.*Sys.Gt2)./Sys.Gt3;
xp(44+1:44+3,1) = Sys.Cnode_3Phase(:,:,5)*([sum(Iline(([1 ]-1)*3+1,1));sum(Iline(([1 ]-1)*3+2,1));sum(Iline(([1 ]-1)*3+3,1))]-[sum(Iline(([5 ]-1)*3+1,1));sum(Iline(([5 ]-1)*3+2,1));sum(Iline(([5 ]-1)*3+3,1))]);
xp(47+1:47+3,1) = Sys.Cnode_3Phase(:,:,6)*([sum(Iline(([2 5 ]-1)*3+1,1));sum(Iline(([2 5 ]-1)*3+2,1));sum(Iline(([2 5 ]-1)*3+3,1))]-[sum(Iline(([6 ]-1)*3+1,1));sum(Iline(([6 ]-1)*3+2,1));sum(Iline(([6 ]-1)*3+3,1))]);
xp(50+1:50+3,1) = Sys.Cnode_3Phase(:,:,7)*([sum(Iline(([6 ]-1)*3+1,1));sum(Iline(([6 ]-1)*3+2,1));sum(Iline(([6 ]-1)*3+3,1))]-[sum(Iline(([7 8 ]-1)*3+1,1));sum(Iline(([7 8 ]-1)*3+2,1));sum(Iline(([7 8 ]-1)*3+3,1))]-[sum(Iload(([1 ]-1)*3+1,1));sum(Iload(([1 ]-1)*3+2,1));sum(Iload(([1 ]-1)*3+3,1))]);
xp(53+1:53+3,1) = Sys.Cnode_3Phase(:,:,8)*([sum(Iline(([7 8 ]-1)*3+1,1));sum(Iline(([7 8 ]-1)*3+2,1));sum(Iline(([7 8 ]-1)*3+3,1))]-[sum(Iline(([9 10 ]-1)*3+1,1));sum(Iline(([9 10 ]-1)*3+2,1));sum(Iline(([9 10 ]-1)*3+3,1))]);
xp(56+1:56+3,1) = Sys.Cnode_3Phase(:,:,9)*([sum(Iline(([9 10 ]-1)*3+1,1));sum(Iline(([9 10 ]-1)*3+2,1));sum(Iline(([9 10 ]-1)*3+3,1))]-[sum(Iline(([11 ]-1)*3+1,1));sum(Iline(([11 ]-1)*3+2,1));sum(Iline(([11 ]-1)*3+3,1))]-[sum(Iload(([2 ]-1)*3+1,1));sum(Iload(([2 ]-1)*3+2,1));sum(Iload(([2 ]-1)*3+3,1))]);
xp(59+1:59+3,1) = Sys.Cnode_3Phase(:,:,10)*([sum(Iline(([4 11 ]-1)*3+1,1));sum(Iline(([4 11 ]-1)*3+2,1));sum(Iline(([4 11 ]-1)*3+3,1))]-[sum(Iline(([12 ]-1)*3+1,1));sum(Iline(([12 ]-1)*3+2,1));sum(Iline(([12 ]-1)*3+3,1))]);
xp(62+1:62+3,1) = Sys.Cnode_3Phase(:,:,11)*([sum(Iline(([3 12 ]-1)*3+1,1));sum(Iline(([3 12 ]-1)*3+2,1));sum(Iline(([3 12 ]-1)*3+3,1))]);
xp(65+1:65+3,1)= Sys.Lline_3Phase(:,:,1)*( [Vpp_a(1); Vpp_b(1); Vpp_c(1)]- Vbus((5-5)*3+1:(5-5)*3+3) -Sys.Gen_R(:,:,1)*Iline((1-1)*3+1:(1-1)*3+3)); % IBR voltage
xp(68+1:68+3,1)= Sys.Lline_3Phase(:,:,2)*( [Vpp_a(2); Vpp_b(2); Vpp_c(2)]- Vbus((6-5)*3+1:(6-5)*3+3) -Sys.Gen_R(:,:,2)*Iline((2-1)*3+1:(2-1)*3+3));
xp(71+1:71+3,1)= Sys.Lline_3Phase(:,:,3)*( [Vpp_a(3); Vpp_b(3); Vpp_c(3)]- Vbus((11-5)*3+1:(11-5)*3+3)-Sys.Gen_R(:,:,3)*Iline((3-1)*3+1:(3-1)*3+3));
xp(74+1:74+3,1)= Sys.Lline_3Phase(:,:,4)*( [Vpp_a(4); Vpp_b(4); Vpp_c(4)]- Vbus((10-5)*3+1:(10-5)*3+3)-Sys.Gen_R(:,:,4)*Iline((4-1)*3+1:(4-1)*3+3));
xp(77+1:77+3,1)=Sys.Lline_3Phase(:,:,5)*( Vbus((5-5)*3+1:(5-5)*3+3)-Vbus((6-5)*3+1:(6-5)*3+3)-Sys.Rline_3Phase(:,:,5) * Iline((5-1)*3+1:(5-1)*3+3) );
xp(80+1:80+3,1)=Sys.Lline_3Phase(:,:,6)*( Vbus((6-5)*3+1:(6-5)*3+3)-Vbus((7-5)*3+1:(7-5)*3+3)-Sys.Rline_3Phase(:,:,6) * Iline((6-1)*3+1:(6-1)*3+3) );
xp(83+1:83+3,1)=Sys.Lline_3Phase(:,:,7)*( Vbus((7-5)*3+1:(7-5)*3+3)-Vbus((8-5)*3+1:(8-5)*3+3)-Sys.Rline_3Phase(:,:,7) * Iline((7-1)*3+1:(7-1)*3+3) );
xp(86+1:86+3,1)=Sys.Lline_3Phase(:,:,8)*( Vbus((7-5)*3+1:(7-5)*3+3)-Vbus((8-5)*3+1:(8-5)*3+3)-Sys.Rline_3Phase(:,:,8) * Iline((8-1)*3+1:(8-1)*3+3) );
xp(89+1:89+3,1)=Sys.Lline_3Phase(:,:,9)*( Vbus((8-5)*3+1:(8-5)*3+3)-Vbus((9-5)*3+1:(9-5)*3+3)-Sys.Rline_3Phase(:,:,9) * Iline((9-1)*3+1:(9-1)*3+3) );
xp(92+1:92+3,1)=Sys.Lline_3Phase(:,:,10)*( Vbus((8-5)*3+1:(8-5)*3+3)-Vbus((9-5)*3+1:(9-5)*3+3)-Sys.Rline_3Phase(:,:,10) * Iline((10-1)*3+1:(10-1)*3+3) );
xp(95+1:95+3,1)=Sys.Lline_3Phase(:,:,11)*( Vbus((9-5)*3+1:(9-5)*3+3)-Vbus((10-5)*3+1:(10-5)*3+3)-Sys.Rline_3Phase(:,:,11) * Iline((11-1)*3+1:(11-1)*3+3) );
xp(98+1:98+3,1)=Sys.Lline_3Phase(:,:,12)*( Vbus((10-5)*3+1:(10-5)*3+3)-Vbus((11-5)*3+1:(11-5)*3+3)-Sys.Rline_3Phase(:,:,12) * Iline((12-1)*3+1:(12-1)*3+3) );
xp(101+1:101+3,1)=Sys.Lload_3Phase(:,:,1)*( Vbus((7-5)*3+1:(7-5)*3+3)-Sys.Rload_3Phase(:,:,1) * Iload((1-1)*3+1:(1-1)*3+3) );
xp(104+1:104+3,1)=Sys.Lload_3Phase(:,:,2)*( Vbus((9-5)*3+1:(9-5)*3+3)-Sys.Rload_3Phase(:,:,2) * Iload((2-1)*3+1:(2-1)*3+3) );
xp(107+1:107+4,1)=(Vt-Vtp)./Sys.T_vt;
% GLO.V1 =[GLO.V1, Vt(1)];
% GLO.V2 =[GLO.V2, Vt(2)];
% GLO.V3 =[GLO.V3, Vt(3)];
% GLO.V4 =[GLO.V4, Vt(4)];
% plot(0:2*10^(-6):t,GLO.V2(1:4:end))
end

48
EMT/generator_states.m
Unescape View File

@ -0,0 +1,48 @@
function [Vt,dta,Id, Iq, Edp, Eqp, Edpp, Eqpp, Pm, Efd, Sys, omg, V1, efd, P1, P2, Epp] = generator_states(mpc, Sys, Base);
%GENERATOR_STATES calculates generator states for a given power system
%% generator inside states
Vtx = mpc.bus(:,8).*cos(mpc.bus(:,9)*pi/180);
Vty = mpc.bus(:,8).*sin(mpc.bus(:,9)*pi/180);
Vt = Vtx + 1i*Vty; % terminal voltage phasor value
Vt_gen = Vt(Sys.GenIdx,:); % terminal voltage of generator bus
GenPQ = (mpc.gen(:,2)+1i*mpc.gen(:,3))./Base;
It_gen = conj(GenPQ./Vt_gen); % current injection
Zbranch = mpc.branch(1,3)+1i*mpc.branch(1,4);
Ibr=(Vt(1,:)-Vt(5,:))/Zbranch;
EQ = Vt_gen + (Sys.Ra+1i*Sys.Xq).*It_gen;
dta = angle(EQ); % rotor angle, angle of Eq, not angle of E''
omg = zeros(length(Sys.GenIdx),1); % rotor speed
% next step, because the past calculation starts from Vt It, which are grid reference based
% thus the E' E'' are also grid state, we need to change it to d,q state
% from grid to d,q, multiply exp(-j(delta-pi/2))
Vt_gen = Vt_gen.*exp(1i*(pi/2-dta)); % d-q based
It_gen = It_gen.*exp(1i*(pi/2-dta));
Id = real(It_gen);
Iq = imag(It_gen);
% Ed''= real(E''), Eq''=imag(E'')
Ep = Vt_gen + Sys.Ra.*It_gen -Sys.Xqp.*Iq+1i*Sys.Xdp.*Id; % d-q based
Edp = real(Ep); % real of E' , d-q based
Eqp = imag(Ep); % imag of E' , d-q based
Epp = Vt_gen + Sys.Ra.*It_gen -Sys.Xqpp.*Iq+1i*Sys.Xdpp.*Id; % d-q based
Edpp = real(Epp); % real of E''
Eqpp = imag(Epp); % imag of E''
Pe = real(GenPQ) + (abs(It_gen).^2).*Sys.Ra; % electrical power
Pm = Pe + Sys.D.*omg./Sys.w0; % mechanical power
% exciter field voltage from 6th-order equation, derivative=0
Efd = (Sys.Xd-Sys.Xdpp)./(Sys.Xdp-Sys.Xdpp).*Eqp - (Sys.Xd-Sys.Xdp)./(Sys.Xdp-Sys.Xdpp).*Eqpp;
efd = Efd.*Sys.Rfd./Sys.Lad; % transfer to efd
% set Vref for exciter equation
% from transfer equation, steady state Vref-Et-V1=0
Sys.Vref = abs(Epp)+ abs(Efd)./Sys.Ek; % abs(Vt_gen) changed to abs(Epp)
V1 = abs(Efd)./Sys.Ek; % V1 is a mid-term in exciter
% from transfer equation,steady state Pref=P1*Gr; P1=P2; Pm=P2;
Sys.Pref = Pm.*Sys.Gr;
P1 = Pm;
P2 = Pm;

39
EMT/myFunction.m
Unescape View File

@ -0,0 +1,39 @@
function Sys = myFunction(mpc)
%% Generator parameters
Sys.H = [6.5 6.5 6.175 6.175]'*9; % base of 900 MW to base 100 MW
Sys.Xl = [0.2 0.2 0.2 0.2]'/9;
Sys.Xd = [1.8 1.8 1.8 1.8]'/9;
Sys.Xq = [1.7 1.7 1.7 1.7]'/9;
Sys.Xdp = [0.3 0.3 0.3 0.3]'/9;
Sys.Xqp = [0.55 0.55 0.55 0.55]'/9;
Sys.Xdpp = [0.25 0.25 0.25 0.25]'/9;
Sys.Xqpp = [0.25 0.25 0.25 0.25]'/9;
Sys.Ra = [0.0025 0.0025 0.0025 0.0025]'/9;
Sys.D = [45 45 45 45]'; % 0 on Kunder book
Sys.Td0p = [8 8 8 8]';
Sys.Tq0p = [0.4 0.4 0.4 0.4]';
Sys.Td0pp= [0.03 0.03 0.03 0.03]';
Sys.Tq0pp= [0.05 0.05 0.05 0.05]';
% add by myself
Sys.L0 = Sys.Xl*2; % L0 in VBR synchronous generator model
Sys.T_vt= 0.05; % time constant for terminal voltage measurement
%% transfer generator parameters EMT model parameters
% Paremeter_transform % provide parameters for VBR generator model
%% contollers
% exciter
Sys.Eta=1; Sys.Etb=10; Sys.Ek=100; Sys.Ete=0.1; Sys.Emax=3; Sys.Emin=0;
% TGov
Sys.Gr= 0.33; Sys.Gt1=2; Sys.Gt2=3; Sys.Gt3=15; Sys.Gdt=0.4; Sys.GVmax=10; Sys.GVmin=0;
%% Index
Sys.w0 = 1.0; % nominal frequency, p.u. useful in swing equation
Sys.GenIdx = [1 2 3 4];
Sys.NonGenIdx = [5 6 7 8 9 10 11];
Sys.LoadIdx = [7 9];
Sys.ShuntIdx = [7 9];
Sys.bus_number = mpc.bus(end,1); % total number of buses

212
EMT/plotEMTResults.m
Unescape View File

@ -0,0 +1,212 @@
function plotEMTResults(res)
close all
t_final = res.t_final;
sol_final = res.sol_final;
% %% relative rotor angles
% figure(1)
% % EMT result
% plot(t_final, (sol_final(:,1)-sol_final(:,1) )*180/pi,'r--','LineWidth',3)
% hold on
% plot(t_final, (sol_final(:,2)-sol_final(:,1) )*180/pi,'g--','LineWidth',3)
% hold on
% plot(t_final, (sol_final(:,3)-sol_final(:,1) )*180/pi,'b--','LineWidth',3)
% hold on
% plot(t_final, (sol_final(:,4)-sol_final(:,1) )*180/pi,'b--','LineWidth',3)
% hold on
% % ylim([-60,10])
% xlabel('Time (s)')
% ylabel('Relative rotor angles (degree)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Relative rotor angles')
% legend('EMT G1','EMT G2','EMT G3','EMT G4')
% %% frequency deviation
% figure(2)
% % EMT result
% plot(t_final, sol_final(:,9),'r--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,10),'g--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,11),'b--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,12),'m--','LineWidth',3)
% hold on
% % ylim([-60,10])
% % xlim([0,0.25])
% xlabel('Time (s)')
% ylabel('Frequency deviation (p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Frequency deviation')
% legend('EMT G1','EMT G2','EMT G3','EMT G4')
% %% Exciter field voltage
% figure(3)
% % EMT result
% plot(t_final, sol_final(:,33),'r--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,34),'g--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,35),'b--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,36),'m--','LineWidth',3)
% xlabel('Time (s)')
% xlim([0,0.4])
% ylabel('Exciter field voltage (p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Exciter field voltage')
% legend('EMT G1','EMT G2','EMT G3','EMT G4')
% %% P1
% figure(4)
% % EMT result
% plot(t_final, sol_final(:,37),'r--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,38),'g--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,39),'b--','LineWidth',3)
% hold on
% plot(t_final, sol_final(:,40),'m--','LineWidth',3)
%
% xlabel('Time (s)')
% ylabel('P1(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('P1')
% legend('EMT G1','EMT G2','EMT G3','EMT G4')
% %% Voltage at bus 7
% figure(5)
% % EMT result
% plot(t_final, sol_final(:,51),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,52),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,53),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Voltage at bus 7(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Voltage at bus 7')
% legend('EMT phase A','EMT phase B','EMT phase C')
%% Voltage at bus 8
figure(6)
% EMT result
plot(t_final, sol_final(:,54),'r','LineWidth',1)
hold on
plot(t_final, sol_final(:,55),'g','LineWidth',1)
hold on
plot(t_final, sol_final(:,56),'b','LineWidth',1)
xlabel('Time (s)')
ylabel('Voltage at bus 8(p.u.)')
set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
title('Voltage at bus 8')
legend('EMT phase A','EMT phase B','EMT phase C')
% %% Voltage at bus 10
% figure(7)
% % EMT result
% plot(t_final, sol_final(:,60),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,61),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,62),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Voltage at bus 10(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Voltage at bus 10')
% legend('EMT phase A','EMT phase B','EMT phase C')
% %% Voltage at bus 6
% figure(8)
% % EMT result
% plot(t_final, sol_final(:,48),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,49),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,50),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Voltage at bus 6(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Voltage at bus 6')
% legend('EMT phase A','EMT phase B','EMT phase C')
% %% current on branch 5-6
% figure(9)
% % EMT result
% plot(t_final, sol_final(:,78),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,79),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,80),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Current on branch 5-6(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Current on branch 5-6')
% legend('EMT phase A','EMT phase B','EMT phase C')
%% current on branch 6-7
figure(10)
% EMT result
plot(t_final, sol_final(:,81),'r','LineWidth',1)
hold on
plot(t_final, sol_final(:,82),'g','LineWidth',1)
hold on
plot(t_final, sol_final(:,83),'b','LineWidth',1)
xlabel('Time (s)')
ylabel('Current on branch 6-7(p.u.)')
set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
title('Current on branch 6-7')
legend('EMT phase A','EMT phase B','EMT phase C')
% %% current on branch 7-8
% figure(11)
% % EMT result
% plot(t_final, sol_final(:,84),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,85),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,86),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Current on branch 7-8(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Current on branch 7-8')
% legend('EMT phase A','EMT phase B','EMT phase C')
% %% current on branch 9-10
% figure(12)
% % EMT result
% plot(t_final, sol_final(:,96),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,97),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,98),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Current on branch 9-10(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Current on branch 9-10')
% legend('EMT phase A','EMT phase B','EMT phase C')
% %% Voltage at bus 6
% figure(13)
% % EMT result
% plot(t_final, sol_final(:,108),'r','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,109),'g','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,110),'b','LineWidth',1)
% hold on
% plot(t_final, sol_final(:,111),'b','LineWidth',1)
%
% xlabel('Time (s)')
% ylabel('Voltage at bus 6(p.u.)')
% set(gca, 'Fontname', 'Times New Roman', 'Fontsize', 12);
% title('Voltage at bus 6')
% legend('Gen1','Gen2','Gen3','Gen4')

9
EMT/setDefaultParameters.m
Unescape View File

@ -0,0 +1,9 @@
function [R_line_positive, wL_line_positive, B_line_positive, int_t, LC_fre, t_base] = setDefaultParameters()
R_line_positive = 0.0001;
wL_line_positive = 0.001;
B_line_positive = 0.00175;
int_t = 0; % could change to any number
LC_fre = 120*pi; % divide all L by w, consider time base is 1/120/pi
t_base =1/(120*pi);
end

4
EMT/setSASParameters.m
Unescape View File

@ -0,0 +1,4 @@
function [orderOfDT, timeStepOfDT] = setSASParameters()
orderOfDT=10;
timeStepOfDT=179*10^(-6);
end

15
EMT/set_simulation_parameters.m
Unescape View File

@ -0,0 +1,15 @@
function [Base, w, FtBus, time_step, time_step_on, TimeRange1, FtTime, TotalTime, TimeRange2, TimeRange3] = set_simulation_parameters()
% This function sets the simulation parameters.
%% settings
Base = 100; % baseMVA
w = 120*pi;
FtBus = 7; % grounding fault bus
time_step = 1*10^(-6); % pre and post fault
% time_step_on = 2*10^(-8); % during fault
time_step_on = 1*10^(-6); % during fault
TimeRange1 = [0,0.1]; % pre-fault simulation time (second)
FtTime = 5/60; % Fault lasting time (second) (5/60 is five cycles)
TotalTime = 0.5; %1; % Total simulation time (second)
TimeRange2 = [TimeRange1(2),TimeRange1(2)+FtTime]; % fault-on simulation period
TimeRange3 = [TimeRange2(2),TotalTime]; % post-fault simulation time

25
EMT/tline_to_matrix.m
Unescape View File

@ -0,0 +1,25 @@
function [R_line, L_line, C_line] = tline_to_matrix(R_line_positive, wL_line_positive, B_line_positive, LC_fre)
% Unit is per unit per km, pu/km
% Assume all the lines have the same per unit parameters
R_line_zero = R_line_positive*6; % from PSCAD line parameters, we can generally get this relationship
wL_line_zero = wL_line_positive*3;
L_line_positive = wL_line_positive/(LC_fre);
L_line_zero = wL_line_zero/(LC_fre);
B_line_zero = B_line_positive/1.5;
C_line_positive = 0.5*B_line_positive/(LC_fre); % multiply 0.5 because it is divided in to two parts
C_line_zero = 0.5*B_line_zero/(LC_fre);
% three phase self and mutual parameters
% Z_zero=Zs+2Zm, Z_positive=Zs-Zm
% so we get: Z_self=1/3*(Z_zero+2*Z_positive), Z_mutual=1/3*(Z_zero-Z_positive)
Rs_line=1/3*(R_line_zero+2*R_line_positive); Rm_line=1/3*(R_line_zero-R_line_positive);
Ls_line=1/3*(L_line_zero+2*L_line_positive); Lm_line=1/3*(L_line_zero-L_line_positive);
Cs_line=1/3*(C_line_zero+2*C_line_positive); Cm_line=1/3*(C_line_zero-C_line_positive);
% R L C three phase matrix
R_line=[Rs_line, Rm_line, Rm_line; Rm_line, Rs_line, Rm_line; Rm_line, Rm_line, Rs_line];
L_line=[Ls_line, Lm_line, Lm_line; Lm_line, Ls_line, Lm_line; Lm_line, Lm_line, Ls_line];
C_line=[Cs_line, Cm_line, Cm_line; Cm_line, Cs_line, Cm_line; Cm_line, Cm_line, Cs_line];
end

30
EMT/transformParameters.m
Unescape View File

@ -0,0 +1,30 @@
function [Sys,Laa0_pp,Lab0_pp,Laa2_pp] = transformParameters(Sys,t_base)
Sys.Lad=Sys.Xd-Sys.Xl;
Sys.Lfd=1./(1./(Sys.Xdp-Sys.Xl)-1./Sys.Lad);
Sys.Rfd=(Sys.Lad+Sys.Lfd)./(Sys.Td0p./t_base);
Sys.L1d=1./(1./(Sys.Xdpp-Sys.Xl)-1./Sys.Lfd-1./Sys.Lad);
Sys.R1d=(1./(1./Sys.Lad+1./Sys.Lfd)+Sys.L1d)./(Sys.Td0pp./t_base);
Sys.Laq=Sys.Xq-Sys.Xl;
Sys.L1q=1./(1./(Sys.Xqp-Sys.Xl)-1./Sys.Laq);
Sys.R1q=(Sys.Laq+Sys.L1q)./(Sys.Tq0p./t_base);
Sys.L2q=1./(1./(Sys.Xqpp-Sys.Xl)-1./Sys.L1q-1./Sys.Laq);
Sys.R2q=(1./(1./Sys.Laq+1./Sys.L1q)+Sys.L2q)./(Sys.Tq0pp./t_base);
Sys.Lad_pp=Sys.Xdpp-Sys.Xl;
Sys.Laq_pp=Sys.Xdpp-Sys.Xl;
Laa0_pp= Sys.Xl+( Sys.Lad_pp+ Sys.Laq_pp+Sys.L0-Sys.Xl )/3;
Lab0_pp= ( 2*Sys.L0-2*Sys.Xl-Sys.Lad_pp-Sys.Laq_pp )/6;
Laa2_pp= ( Sys.Lad_pp-Sys.Laq_pp )/3;
Sys.Vd_pp_coeff1=-(Sys.Rfd./(Sys.Lfd).^2+Sys.R1d./(Sys.L1d).^2).*(Sys.Lad_pp).^2;
Sys.Vd_pp_coeff2=-(Sys.Lad_pp.*Sys.Rfd./(Sys.Lfd).^2.*(1-Sys.Lad_pp./Sys.Lfd) - Sys.Lad_pp.^2.*Sys.R1d./(Sys.L1d).^2./Sys.Lfd);
Sys.Vd_pp_coeff3=(Sys.Lad_pp.^2.*Sys.Rfd./(Sys.Lfd).^2./Sys.L1d - Sys.Lad_pp.*Sys.R1d./(Sys.L1d).^2.*(1-Sys.Lad_pp./Sys.L1d));
Sys.Vq_pp_coeff1=-(Sys.R1q./(Sys.L1q).^2+Sys.R2q./(Sys.L2q).^2).*(Sys.Laq_pp).^2;
Sys.Vq_pp_coeff2=-(Sys.Laq_pp.*Sys.R1q./(Sys.L1q).^2.*(1-Sys.Laq_pp./Sys.L1q) - Sys.Laq_pp.^2.*Sys.R2q./(Sys.L2q).^2./Sys.L1q);
Sys.Vq_pp_coeff3=(Sys.Laq_pp.^2.*Sys.R1q./(Sys.L1q).^2./Sys.L2q - Sys.Laq_pp.*Sys.R2q./(Sys.L2q).^2.*(1-Sys.Laq_pp./Sys.L2q));
end

82
EMT/twoarea.m
Unescape View File

@ -0,0 +1,82 @@
function mpc = twoarea
%2_AREA_PSSE30
% April 12, 2016 13:35:48; Simulator Version 19; BuildDate 2016_2_6
%
% Converted by MATPOWER 6.0 using PSSE2MPC on 08-Apr-2018
% from '2-area_psse30.raw' using PSS/E rev 30 format.
%
% WARNINGS:
% Conversion explicitly using PSS/E revision 30
% Section label mismatch, found 'VSC DC LINE', expected 'VOLTAGE SOURCE CONVERTER'
% Section label mismatch, found 'IMPEDANCE CORRECTION TABLE', expected 'IMPEDANCE CORRECTION'
% Skipped 1 line of zone data.
% Skipped 1 line of owner data.
% Section label mismatch, found 'FACTS', expected 'FACTS CONTROL DEVICE'
% Using default voltage magnitude limits: VMIN = 0.9 p.u., VMAX = 1.1 p.u.
%
% See CASEFORMAT for details on the MATPOWER case file format.
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 100;
%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
mpc.bus = [
1 3 0 0 0 0 1 1.03 20.07 20 1 1.1 0.9;
2 2 0 0 0 0 1 1.01 10.30 20 1 1.1 0.9;
3 2 0 0 0 0 1 1.03 -7.0145 20 1 1.1 0.9;
4 2 0 0 0 0 1 1.01 -17.20 20 1 1.1 0.9;
5 1 0 0 0 0 1 1.00645 13.6071 230 1 1.1 0.9;
6 1 0 0 0 0 1 0.97812 3.5208 230 1 1.1 0.9;
7 1 967 100 0 200 1 0.9610 -4.889 230 1 1.1 0.9;
8 1 0 0 0 0 1 0.9486 -18.76 230 1 1.1 0.9;
9 1 1767 100 0 350 1 0.9714 -32.364 230 1 1.1 0.9;
10 1 0 0 0 0 1 0.9835 -23.949 230 1 1.1 0.9;
11 1 0 0 0 0 1 1.0083 -13.6407 230 1 1.1 0.9;
];
%% generator data
% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf mu_Pmax mu_Pmin mu_Qmax mu_Qmin
mpc.gen = [
1 700.10 185 300 -300 1.03 900 1 9999 -9999 0 0 0 0 0 0 0 0 0 0 0 0.0000 0.0000 0.0000 0.0000;
2 700 234.67 300 -300 1.01 900 1 9999 -9999 0 0 0 0 0 0 0 0 0 0 0 0.0000 0.0000 0.0000 0.0000;
3 719 176 300 -300 1.03 900 1 9999 -9999 0 0 0 0 0 0 0 0 0 0 0 0.0000 0.0000 0.0000 0.0000;
4 700 202.07 300 -300 1.01 900 1 9999 -9999 0 0 0 0 0 0 0 0 0 0 0 0.0000 0.0000 0.0000 0.0000;
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [
1 5 0 0.15/9 0 0 0 0 0 0 1 -360 360; % transformer. based given in Kunder's book is 900
2 6 0 0.15/9 0 0 0 0 0 0 1 -360 360; % So, using a base of 100, we need to
3 11 0 0.15/9 0 0 0 0 0 0 1 -360 360;
4 10 0 0.15/9 0 0 0 0 0 0 1 -360 360;
5 6 0.0025 0.025 0.04375 0 0 0 0 0 1 -360 360;
6 7 0.0010 0.010 0.01750 0 0 0 0 0 1 -360 360;
7 8 0.0110 0.110 0.19250 0 0 0 0 0 1 -360 360;
7 8 0.0110 0.110 0.19250 0 0 0 0 0 1 -360 360;
8 9 0.0110 0.110 0.19250 0 0 0 0 0 1 -360 360;
8 9 0.0110 0.110 0.19250 0 0 0 0 0 1 -360 360;
9 10 0.0010 0.010 0.01750 0 0 0 0 0 1 -360 360;
10 11 0.0025 0.025 0.04375 0 0 0 0 0 1 -360 360;
];
%% bus names
mpc.bus_name = {
'GEN BUS1 ';
'GEN BUS2 ';
'GEN BUS3 ';
'GEN BUS4 ';
'STA B230 ';
'STA C230 ';
'STA 2 7 ';
'STA A230 ';
'STA 3 9 ';
'STA 3 10 ';
'STA 3 11 ';
};

13
EMT/update_system_parameters.m
Unescape View File

@ -0,0 +1,13 @@
function Sys = update_system_parameters(Sys, Laa0_pp, Lab0_pp)
% This function updates the system parameters including Emax, Emin, Gen_R and Lline_3Phase.
Sys.Emax=Sys.Emax.*Sys.Rfd./Sys.Lad;
Sys.Emin=Sys.Emin.*Sys.Rfd./Sys.Lad;
for k=1:1:length(Sys.GenIdx)
Sys.Gen_R(:,:,k)= [Sys.Ra(k), 0, 0; 0, Sys.Ra(k), 0; 0, 0, Sys.Ra(k)];
L_stator(:,:,k)= [Laa0_pp(k) Lab0_pp(k) Lab0_pp(k); Lab0_pp(k) Laa0_pp(k) Lab0_pp(k); Lab0_pp(k) Lab0_pp(k) Laa0_pp(k)]/120/pi;
% because Sys.Laa0_pp is already updated
Sys.Lline_3Phase(:,:,k) = Sys.Lline_3Phase(:,:,k)+ L_stator(:,:,k);
end

1
docs/source/data_and_models.rst
Unescape View File

@ -9,6 +9,7 @@ Currently PowerSAS.m supports extended PSAT (Matlab) data format.
In addition, PSS/E data format can be easily converted to the PSAT data format
using the ``psse2mpc`` function by MATPOWER and the ``matpower2psat`` function
by PSAT. A sample code to convert a PSS/E *raw file to the PSAT format is below:
.. code:: matlab
casename = 'ieee14'

7
example/ex_emt.m
Unescape View File

@ -0,0 +1,7 @@
clc; clear; close;
addpath('EMT')
res = runPowerSAS('emt','twoarea');
plotEMTResults(res);

6
example/output/~default_log.txt
Unescape View File

@ -0,0 +1,6 @@
03:20:33[INFO]The session ID is 20230425T032033017.
03:20:33[INFO]CASE twoarea. Timestamp=20230425T032033017.
03:20:33[INFO]Finished loading data file.
03:20:41[INFO]

79
interface/data/d_mpc_temp.m
Unescape View File

@ -1,79 +0,0 @@
% 03/16/23 File data originated from Matpower data file
%
Bus.con = [ ...
1 69 1.03 0 1 1;
2 69 1.02 -0.02798 1 1;
3 69 1 -0.0601 1 1;
4 69 0.9986 -0.07472 1 1;
5 69 1.004 -0.06432 1 1;
6 138 0.9987 -0.11 2 2;
7 138 1.007 -0.08429 2 2;
8 69 1.019 -0.02434 2 2;
9 138 1.002 -0.1275 2 2;
10 138 0.9935 -0.1302 2 2;
11 138 0.9925 -0.1229 2 2;
12 138 0.9864 -0.1289 2 2;
13 138 0.984 -0.1338 2 2;
14 138 0.9906 -0.1669 2 2;
];
SW.con = [ ...
1 100 69 1.03 0 1 -0.5 1.1 0.9 0.8144 1 1 1;
];
PV.con = [ ...
2 100 69 0.4 1.03 0.15 -0.4 1.1 0.9 1 1;
3 100 69 0.4 1.01 0.15 -0.1 1.1 0.9 1 1;
6 100 138 0.3 1.03 0.1 -0.06 1.1 0.9 1 1;
8 100 69 0.35 1.03 0.1 -0.06 1.1 0.9 1 1;
];
PQ.con = [ ...
2 100 69 0.217 0.127 1.1 0.9 0 1;
3 100 69 0.5 0.25 1.1 0.9 0 1;
4 100 69 0.478 0.1 1.1 0.9 0 1;
5 100 69 0.076 0.016 1.1 0.9 0 1;
6 100 138 0.15 0.075 1.1 0.9 0 1;
9 100 138 0.295 0.166 1.1 0.9 0 1;
10 100 138 0.09 0.058 1.1 0.9 0 1;
11 100 138 0.035 0.018 1.1 0.9 0 1;
12 100 138 0.061 0.016 1.1 0.9 0 1;
13 100 138 0.135 0.058 1.1 0.9 0 1;
14 100 138 0.2 0.07 1.1 0.9 0 1;
7 100 138 0 0 1.1 0.9 0 1;
];
Shunt.con = [ ...
9 100 138 60 0 0.19 1;
14 100 138 60 0 0.15 1;
];
Line.con = [ ...
1 2 100 69 60 0 0 0.01938 0.05917 0.0528 0 0 1 1 0 1;
1 5 100 69 60 0 0 0.05403 0.223 0.0492 0 0 1 1 0 1;
2 3 100 69 60 0 0 0.04699 0.198 0.0438 0 0 1 1 0 1;
2 4 100 69 60 0 0 0.05811 0.1763 0.034 0 0 1 1 0 1;
2 5 100 69 60 0 0 0.05695 0.1739 0.0346 0 0 1 1 0 1;
3 4 100 69 60 0 0 0.06701 0.171 0.0128 0 0 1 1 0 1;
4 5 100 69 60 0 0 0.01335 0.04211 0 0 0 1 1 0 1;
6 11 100 138 60 0 0 0.09498 0.1989 0 0 0 1 1 0 1;
6 12 100 138 60 0 0 0.1229 0.2558 0 0 0 1 1 0 1;
6 13 100 138 60 0 0 0.06615 0.1303 0 0 0 1 1 0 1;
7 9 100 138 60 0 0 -0 0.11 0 0 0 1 1 0 1;
9 10 100 138 60 0 0 0.03181 0.0845 0 0 0 1 1 0 1;
9 14 100 138 60 0 0 0.1271 0.2704 0 0 0 1 1 0 1;
10 11 100 138 60 0 0 0.08205 0.1921 0 0 0 1 1 0 1;
12 13 100 138 60 0 0 0.2209 0.1999 0 0 0 1 1 0 1;
13 14 100 138 60 0 0 0.1709 0.348 0 0 0 1 1 0 1;
4 7 100 69 60 0 0.5 0 0.2091 0 0.9968 0 0.2 0.2 0 1;
4 9 100 69 60 0 0.5 0 0.5562 0 0.9968 0 0.2 0.2 0 1;
6 5 100 138 60 0 2 0 0.252 0 0.9968 0 0.5 0.5 0 1;
8 7 100 69 60 0 0.5 0 0.1762 0 0.9968 0 0.5 0.5 0 1;
];
Bus.names = {...
'BUS1'; 'BUS2'; 'BUS3'; 'BUS4'; 'BUS5';
'BUS6'; 'BUS7'; 'BUS8'; 'BUS9'; 'BUS10';
'BUS11'; 'BUS12'; 'BUS13'; 'BUS14'};

2
runPowerSAS.m
Unescape View File

@ -99,6 +99,8 @@ end
%% main body
if strcmp(simType,'pf')
res=calPF(SysData,options,caseName);
elseif strcmp(simType,'emt')
res=calEMT(data);
elseif strcmp(simType,'cpf')
res=calCPF(SysData,options,caseName,varargin{:});
elseif strcmp(simType,'tsa')

Write Preview
Loading…
Cancel
Save