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

Compare commits

base: ANL-CEEESA:5b6625a1428f25e4093bb6db6271f7b7e41270cd
Branches Tags
ANL-CEEESA:release
ANL-CEEESA:develop
...
compare: ANL-CEEESA:59d62c45d2c03f4baa3013029d0af8b94a780f86
Branches Tags
ANL-CEEESA:release
ANL-CEEESA:develop

7 Commits
5b6625a142 ... 59d62c45d2

Author SHA1 Message Date
yangliu0101 59d62c45d2
Merge pull request #5 from yangliu0101/doc
2 years ago
yangliu0101 d1f3c13664
Merge pull request #4 from yangliu0101/release
2 years ago
Yang Liu 68c835ca1e enable emt in runPowerSAS.m function
2 years ago
Yang Liu e3a4af8e2c delete wrong files
2 years ago
Yang Liu 16c1547ba9 add basic version of EMT module
2 years ago
Yang Liu bb6526ed86 debug the layout of PSS/E interface
3 years ago
Yang Liu f025da30b2 add interface with PSS/E
3 years ago
34 changed files with 2080 additions and 0 deletions
Show Stats

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

13
docs/source/data_and_models.rst
Unescape View File

@ -5,6 +5,19 @@ Data and Models
~~~~~~~~~~~~~~~~~~~~~~~~~
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'
[mpc, warnings] = psse2mpc(horzcat(casename,'.raw'),horzcat(casename,'.m'));
check = matpower2psat(horzcat(casename,'.m'), pwd);
Support for other formats and data format conversion features will be
added in future versions.

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_ieee14.m
Unescape View File

@ -0,0 +1,79 @@
% 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'};

161
interface/data/d_ieee39.m
Unescape View File

@ -0,0 +1,161 @@
% 03/16/23 File data originated from Matpower data file
%
Bus.con = [ ...
1 345 1.031 -0.2141 1 1;
2 345 1.002 -0.2502 1 1;
3 345 0.9535 -0.3343 1 1;
4 345 0.931 -0.3363 1 2;
5 345 0.9174 -0.3049 1 2;
6 345 0.9186 -0.2919 1 2;
7 345 0.8626 -0.3252 1 2;
8 345 0.88 -0.3322 1 2;
9 345 0.9788 -0.2449 1 2;
10 345 0.9562 -0.2565 1 2;
11 345 0.9429 -0.2703 1 2;
12 138 0.9449 -0.3003 1 2;
13 345 0.9523 -0.2722 1 2;
14 345 0.9465 -0.3029 2 3;
15 345 0.9565 -0.314 2 3;
16 345 0.9767 -0.2878 2 3;
17 345 0.9749 -0.3076 1 1;
18 345 0.9652 -0.3262 1 1;
19 345 1.03 -0.2038 2 4;
20 138 0.9798 -0.2297 2 4;
21 345 0.993 -0.2428 2 3;
22 345 1.029 -0.1615 2 3;
23 345 1.023 -0.1651 2 4;
24 345 0.9871 -0.2857 2 3;
25 345 1.021 -0.2356 1 1;
26 345 1.01 -0.2661 1 1;
27 345 0.9876 -0.3071 1 1;
28 345 1.02 -0.2014 2 3;
29 345 1.024 -0.1509 2 3;
30 34.5 1.005 -0.2041 1 1;
31 34.5 0.9064 -0.1583 1 2;
32 21 0.9752 -0.1068 1 2;
33 21 0.9972 -0.112 2 4;
34 15.5 1.012 -0.1387 2 4;
35 15.5 1.049 -0.07315 2 3;
36 12.5 1.063 -0.02552 2 4;
37 12.5 1.028 -0.1139 1 1;
38 34.5 1.026 -0.02585 2 3;
39 345 1.03 -0.1913 1 2;
];
SW.con = [ ...
39 1199 345 1.03 -0.1913 0.4794 -0.1445 1.05 0.95 0.6012 1 1 1;
];
PV.con = [ ...
30 275 34.5 0.9091 1.048 0.557 -0.2112 1.05 0.95 1 1;
31 836 34.5 0.6853 1.04 0.5141 -0.1467 1.05 0.95 1 1;
32 843.7 21 0.7704 0.9831 0.5295 -0.2136 1.05 0.95 1 1;
33 1175 21 0.538 0.9972 0.4666 -0.182 1.05 0.95 1 1;
34 1080 15.5 0.4703 1.012 0.5661 -0.1741 1.05 0.95 1 1;
35 1086 15.5 0.5987 1.049 0.5469 -0.2164 1.05 0.95 1 1;
36 1025 12.5 0.5462 1.063 0.5544 -0.243 1.05 0.95 1 1;
37 970.2 12.5 0.5566 1.028 0.4571 -0.2228 1.05 0.95 1 1;
38 1684 34.5 0.4928 1.026 0.4957 -0.2119 1.05 0.95 1 1;
10 900 345 0.1111 1.004 0.3287 -0.3287 1.05 0.95 1 1;
20 900 138 0.1667 1.004 0.3287 -0.3287 1.05 0.95 1 1;
2 900 345 0.2222 1.004 0.3287 -0.3287 1.05 0.95 1 1;
25 900 345 0.2778 1.004 0.3287 -0.3287 1.05 0.95 1 1;
];
PQ.con = [ ...
3 100 345 6 2.5 1.05 0.95 0 1;
4 100 345 4.5 1.84 1.05 0.95 0 1;
7 100 345 2.338 8.4 1.05 0.95 0 1;
8 100 345 5.22 1.76 1.05 0.95 0 1;
12 100 138 1.2 0.3 1.05 0.95 0 1;
15 100 345 3.2 1.53 1.05 0.95 0 1;
16 100 345 3.294 3.23 1.05 0.95 0 1;
18 100 345 1.58 0.3 1.05 0.95 0 1;
20 100 138 6.8 1.03 1.05 0.95 0 1;
21 100 345 2.74 1.15 1.05 0.95 0 1;
23 100 345 2.475 0.846 1.05 0.95 0 1;
24 100 345 3.086 -0.922 1.05 0.95 0 1;
25 100 345 2.24 0.472 1.05 0.95 0 1;
26 100 345 1.39 0.17 1.05 0.95 0 1;
27 100 345 2.81 0.755 1.05 0.95 0 1;
28 100 345 2.06 0.276 1.05 0.95 0 1;
29 100 345 2.835 1.269 1.05 0.95 0 1;
31 100 34.5 0.8 0.4 1.05 0.95 0 1;
39 100 345 4 2.5 1.05 0.95 0 1;
1 100 345 0 0 1.05 0.95 0 1;
5 100 345 0 0 1.05 0.95 0 1;
6 100 345 0 0 1.05 0.95 0 1;
9 100 345 0 0 1.05 0.95 0 1;
11 100 345 0 0 1.05 0.95 0 1;
13 100 345 0 0 1.05 0.95 0 1;
14 100 345 0 0 1.05 0.95 0 1;
17 100 345 0 0 1.05 0.95 0 1;
19 100 345 0 0 1.05 0.95 0 1;
22 100 345 0 0 1.05 0.95 0 1;
];
Shunt.con = [ ...
4 100 345 60 0 1 1;
5 100 345 60 0 2 1;
];
Line.con = [ ...
1 2 100 345 60 0 0 0.0035 0.0411 0.6987 0 0 1.75 1.92 1.92 1;
1 39 100 345 60 0 0 0.001 0.025 0.75 0 0 1.25 1.37 1.37 1;
2 3 100 345 60 0 0 0.0013 0.0151 0.2572 0 0 6.5 7.15 7.15 1;
2 25 100 345 60 0 0 0.007 0.0086 0.146 0 0 2.6 2.86 2.86 1;
3 4 100 345 60 0 0 0.0013 0.0213 0.2214 0 0 1.25 1.37 1.37 1;
3 18 100 345 60 0 0 0.0011 0.0133 0.2138 0 0 1.25 1.37 1.37 1;
4 5 100 345 60 0 0 0.0008 0.0128 0.1342 0 0 3 3.3 3.3 1;
4 14 100 345 60 0 0 0.0008 0.0129 0.1382 0 0 3 3.3 3.3 1;
5 6 100 345 60 0 0 0.0002 0.0026 0.0434 0 0 4.62 5.08 5.08 1;
5 8 100 345 60 0 0 0.0008 0.0112 0.1476 0 0 4.62 5.08 5.08 1;
6 7 100 345 60 0 0 0.0006 0.0092 0.113 0 0 7 7.7 7.7 1;
6 11 100 345 60 0 0 0.0007 0.0082 0.1389 0 0 4.66 5.12 5.12 1;
7 8 100 345 60 0 0 0.0004 0.0046 0.078 0 0 4.62 5.08 5.08 1;
8 9 100 345 60 0 0 0.0023 0.0363 0.3804 0 0 4.66 5.12 5.12 1;
9 39 100 345 60 0 0 0.001 0.025 1.2 0 0 4 4.4 4.4 1;
10 11 100 345 60 0 0 0.0004 0.0043 0.0729 0 0 4.62 5.08 5.08 1;
10 13 100 345 60 0 0 0.0004 0.0043 0.0729 0 0 4.29 4.71 4.71 1;
13 14 100 345 60 0 0 0.0009 0.0101 0.1723 0 0 3.8 4.18 4.18 1;
14 15 100 345 60 0 0 0.0018 0.0217 0.366 0 0 1.25 1.37 1.37 1;
15 16 100 345 60 0 0 0.0009 0.0094 0.171 0 0 4.29 4.71 4.71 1;
16 17 100 345 60 0 0 0.0007 0.0089 0.1342 0 0 2.6 2.86 2.86 1;
16 19 100 345 60 0 0 0.0016 0.0195 0.304 0 0 5.52 6.07 6.07 1;
16 21 100 345 60 0 0 0.0008 0.0135 0.2548 0 0 4.29 4.71 4.71 1;
16 24 100 345 60 0 0 0.0003 0.0059 0.068 0 0 2.3 2.53 2.53 1;
17 18 100 345 60 0 0 0.0007 0.0082 0.1319 0 0 3 3.3 3.3 1;
17 27 100 345 60 0 0 0.0013 0.0173 0.3216 0 0 1.25 1.37 1.37 1;
21 22 100 345 60 0 0 0.0008 0.014 0.2565 0 0 6.5 7.15 7.15 1;
22 23 100 345 60 0 0 0.0006 0.0096 0.1846 0 0 1.25 1.37 1.37 1;
23 24 100 345 60 0 0 0.0022 0.035 0.361 0 0 4.66 5.12 5.12 1;
25 26 100 345 60 0 0 0.0032 0.0323 0.513 0 0 1.25 1.37 1.37 1;
26 27 100 345 60 0 0 0.0014 0.0147 0.2396 0 0 3.8 4.18 4.18 1;
26 28 100 345 60 0 0 0.0043 0.0474 0.7802 0 0 1.75 1.92 1.92 1;
26 29 100 345 60 0 0 0.0057 0.0625 1.029 0 0 2.3 2.53 2.53 1;
28 29 100 345 60 0 0 0.0014 0.0151 0.249 0 0 4 4.4 4.4 1;
2 30 100 345 60 0 10 0 0.0181 0 1.025 0 3.8 4.18 4.18 1;
31 6 100 34.5 60 0 0.1 0 0.025 0 0.9 0 7 7.7 7.7 1;
10 32 100 345 60 0 16.43 0 0.02 0 1.07 0 8.16 8.97 8.97 1;
12 11 100 138 60 0 0.4 0.0016 0.0435 0 1.006 0 1.25 1.37 1.37 1;
12 13 100 138 60 0 0.4 0.0016 0.0435 0 1.006 0 1.25 1.37 1.37 1;
19 20 100 345 60 0 2.5 0.0007 0.0138 0 1.06 0 2.3 2.53 2.53 1;
19 33 100 345 60 0 16.43 0.0007 0.0142 0 1.07 0 7 7.7 7.7 1;
20 34 100 138 60 0 8.903 0.0009 0.018 0 1.009 0 5.6 6.16 6.16 1;
22 35 100 345 60 0 22.26 0 0.0143 0 1.025 0 7.15 7.86 7.86 1;
23 36 100 345 60 0 27.6 0.0005 0.0272 0 1 0 5.6 6.16 6.16 1;
25 37 100 345 60 0 27.6 0.0006 0.0232 0 1.025 0 5.52 6.07 6.07 1;
29 38 100 345 60 0 10 0.0008 0.0156 0 1.025 0 8.37 9.2 9.2 1;
];
Bus.names = {...
'BUS1'; 'BUS2'; 'LOAD3'; 'LOAD4'; 'BUS5';
'BUS6'; 'LOAD7'; 'LOAD8'; 'BUS9'; 'BUS10';
'BUS11'; 'LOAD12'; 'BUS13'; 'BUS14'; 'LOAD15';
'LOAD16'; 'BUS17'; 'LOAD18'; 'BUS19'; 'LOAD20';
'LOAD21'; 'BUS22'; 'LOAD23'; 'LOAD24'; 'LOAD25';
'LOAD26'; 'LOAD27'; 'LOAD28'; 'LOAD29'; 'GEN30';
'GEN31'; 'GEN32'; 'GEN33'; 'GEN34'; 'GEN35';
'GEN36'; 'GEN37'; 'GEN38'; 'GEN39'};

57
interface/data/d_kundur.m
Unescape View File

@ -0,0 +1,57 @@
% 03/16/23 File data originated from Matpower data file
%
Bus.con = [ ...
1 20 1 0.5703 1 1;
2 20 1 0.3779 1 1;
3 20 1 0.1957 2 1;
4 20 1 0.3777 2 1;
5 230 0.9834 0.4826 1 1;
6 230 0.9691 0.2935 1 1;
7 230 0.9562 0.1425 1 1;
8 230 0.954 -0.03717 2 1;
9 230 0.9686 0.1113 2 1;
10 230 0.9838 0.2933 2 1;
];
SW.con = [ ...
1 900 20 1 0.5703 0.6667 0 1.1 0.9 0.8287 1 1 1;
];
PV.con = [ ...
2 900 20 0.7778 1 0.6667 -0.6667 1.1 0.9 1 1;
3 900 20 0.7778 1 0.6667 -0.6667 1.1 0.9 1 1;
4 900 20 0.7778 1 0.6667 -0.6667 1.1 0.9 1 1;
];
PQ.con = [ ...
7 100 230 11.59 -0.735 1.1 0.9 0 1;
8 100 230 15.75 -0.899 1.1 0.9 0 1;
5 100 230 0 0 1.1 0.9 0 1;
6 100 230 0 0 1.1 0.9 0 1;
9 100 230 0 0 1.1 0.9 0 1;
10 100 230 0 0 1.1 0.9 0 1;
];
Line.con = [ ...
5 6 100 230 60 0 0 0.005 0.05 0.075 0 0 0 0 0 1;
5 6 100 230 60 0 0 0.00501 0.05001 0.075 0 0 0 0 0 1;
6 7 100 230 60 0 0 0.002 0.02 0.03 0 0 0 0 0 1;
6 7 100 230 60 0 0 0.00201 0.02001 0.03 0 0 0 0 0 1;
7 8 100 230 60 0 0 0.02201 0.22 0.33 0 0 0 0 0 1;
7 8 100 230 60 0 0 0.02202 0.22 0.33 0 0 0 0 0 1;
7 8 100 230 60 0 0 0.022 0.22 0.33 0 0 0 0 0 1;
8 9 100 230 60 0 0 0.002 0.02 0.03 0 0 0 0 0 1;
8 9 100 230 60 0 0 0.00201 0.02001 0.03 0 0 0 0 0 1;
9 10 100 230 60 0 0 0.005 0.05 0.075 0 0 0 0 0 1;
9 10 100 230 60 0 0 0.00501 0.05001 0.075 0 0 0 0 0 1;
1 5 100 20 60 0 0.08696 0.001 0.012 0 1 0 0 0 0 1;
2 6 100 20 60 0 0.08696 0.001 0.012 0 1 0 0 0 0 1;
3 9 100 20 60 0 0.08696 0.001 0.012 0 1 0 0 0 0 1;
4 10 100 20 60 0 0.08696 0.001 0.012 0 1 0 0 0 0 1;
];
Bus.names = {...
'1'; '2'; '12'; '11'; '101';
'102'; '3'; '13'; '112'; '111'};

94
interface/data/ieee14.m
Unescape View File

@ -0,0 +1,94 @@
function mpc = ieee14
%IEEE14
% PSS(R)E 32 RAW created by rawd32 SUN, MAY 17 2020 20:38
% IEEE 14 BUS TEST CASE
% 03/06/14 CONTO 100.0 1962 W
%
% Converted by MATPOWER 7.1 using PSSE2MPC on 16-Mar-2023
% from 'ieee14.raw' using PSS/E rev 32 format.
%
% WARNINGS:
% Skipped 2 lines of zone data.
% Skipped 2 lines of owner data.
% 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 0 69 1 1.1 0.9;
2 2 21.7 12.7 0 0 1 1.0197 -1.6032 69 1 1.1 0.9;
3 2 50 25 0 0 1 1.00042 -3.4433 69 1 1.1 0.9;
4 1 47.8 10 0 0 1 0.99858 -4.2812 69 1 1.1 0.9;
5 1 7.6 1.6 0 0 1 1.00443 -3.685 69 1 1.1 0.9;
6 2 15 7.5 0 0 2 0.99871 -6.3024 138 2 1.1 0.9;
7 1 0 0 0 0 2 1.00682 -4.8292 138 2 1.1 0.9;
8 2 0 0 0 0 2 1.01895 -1.3945 69 2 1.1 0.9;
9 1 29.5 16.6 0 19 2 1.00193 -7.3053 138 2 1.1 0.9;
10 1 9 5.8 0 0 2 0.99351 -7.46 138 2 1.1 0.9;
11 1 3.5 1.8 0 0 2 0.99245 -7.0444 138 2 1.1 0.9;
12 1 6.1 1.6 0 0 2 0.98639 -7.3874 138 2 1.1 0.9;
13 1 13.5 5.8 0 0 2 0.98403 -7.6654 138 2 1.1 0.9;
14 1 20 7 0 15 2 0.99063 -9.5636 138 2 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
mpc.gen = [
1 81.442 1.962 100 -50 1.03 100 1 200 50 0 0 0 0 0 0 0 0 0 0 0;
2 40 15 15 -40 1.03 100 1 50 10 0 0 0 0 0 0 0 0 0 0 0;
3 40 15 15 -10 1.01 100 1 50 10 0 0 0 0 0 0 0 0 0 0 0;
6 30 10 10 -6 1.03 100 1 50 10 0 0 0 0 0 0 0 0 0 0 0;
8 35 10 10 -6 1.03 100 1 50 10 0 0 0 0 0 0 0 0 0 0 0;
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [
1 2 0.01938 0.05917 0.0528 100 100 0 0 0 1 -360 360;
1 5 0.05403 0.22304 0.0492 100 100 0 0 0 1 -360 360;
2 3 0.04699 0.19797 0.0438 100 100 0 0 0 1 -360 360;
2 4 0.05811 0.17632 0.034 100 100 0 0 0 1 -360 360;
2 5 0.05695 0.17388 0.0346 100 100 0 0 0 1 -360 360;
3 4 0.06701 0.17103 0.0128 100 100 0 0 0 1 -360 360;
4 5 0.01335 0.04211 0 100 100 0 0 0 1 -360 360;
6 11 0.09498 0.1989 0 100 100 0 0 0 1 -360 360;
6 12 0.12291 0.25581 0 100 100 0 0 0 1 -360 360;
6 13 0.06615 0.13027 0 100 100 0 0 0 1 -360 360;
7 9 -0 0.11001 0 100 100 0 0 0 1 -360 360;
9 10 0.03181 0.0845 0 100 100 0 0 0 1 -360 360;
9 14 0.12711 0.27038 0 100 100 0 0 0 1 -360 360;
10 11 0.08205 0.19207 0 100 100 0 0 0 1 -360 360;
12 13 0.22092 0.19988 0 100 100 0 0 0 1 -360 360;
13 14 0.17093 0.34802 0 100 100 0 0 0 1 -360 360;
4 7 0 0.20912 0 20 20 0 0.99677 0 1 -360 360;
4 9 0 0.55618 0 20 20 0 0.99677 0 1 -360 360;
6 5 0 0.25202 0 50 50 0 0.99677 0 1 -360 360;
8 7 0 0.17615 0 50 50 0 0.99677 0 1 -360 360;
];
%% bus names
mpc.bus_name = {
'BUS1 ';
'BUS2 ';
'BUS3 ';
'BUS4 ';
'BUS5 ';
'BUS6 ';
'BUS7 ';
'BUS8 ';
'BUS9 ';
'BUS10 ';
'BUS11 ';
'BUS12 ';
'BUS13 ';
'BUS14 ';
};

181
interface/data/ieee39.m
Unescape View File

@ -0,0 +1,181 @@
function mpc = ieee39
%IEEE39
% PSS(R)E-33.11 FRI, APR 13 2018 10:22
% NEW ENGLAND TEST SYSTEM, 39 BUSES, 9 GENERATORS, 2 Areas, 4 Zones
% UPDATED BY JOSE CONTO
%
% Converted by MATPOWER 7.1 using PSSE2MPC on 16-Mar-2023
% from 'ieee39.raw' using PSS/E rev 33 format.
%
% WARNINGS:
% Section label mismatch, found 'VSC DC LINE', expected 'VOLTAGE SOURCE CONVERTER'
% Skipped 4 lines of zone data.
% Skipped 2 lines of owner data.
% Section label mismatch, found 'FACTS DEVICE', expected 'FACTS CONTROL DEVICE'
% Section label mismatch, found 'GNE', expected 'GNE DEVICE'
%
% 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 1 0 0 0 0 1 1.03145 -12.2668 345 1 1.05 0.95;
2 2 0 0 0 0 1 1.00247 -14.3363 345 1 1.05 0.95;
3 1 600 250 0 0 1 0.95348 -19.1513 345 1 1.05 0.95;
4 1 450 184 0 100 1 0.93103 -19.2658 345 2 1.05 0.95;
5 1 0 0 0 200 1 0.91744 -17.4681 345 2 1.05 0.95;
6 1 0 0 0 0 1 0.91858 -16.7257 345 2 1.05 0.95;
7 1 233.8 840 0 0 1 0.86255 -18.6341 345 2 1.05 0.95;
8 1 522 176 0 0 1 0.87995 -19.032 345 2 1.05 0.95;
9 1 0 0 0 0 1 0.9788 -14.0303 345 2 1.05 0.95;
10 2 0 0 0 0 1 0.95616 -14.6962 345 2 1.05 0.95;
11 1 0 0 0 0 1 0.94295 -15.4844 345 2 1.05 0.95;
12 1 120 30 0 0 1 0.94491 -17.2039 138 2 1.05 0.95;
13 1 0 0 0 0 1 0.95234 -15.5933 345 2 1.05 0.95;
14 1 0 0 0 0 2 0.94651 -17.355 345 3 1.05 0.95;
15 1 320 153 0 0 2 0.95647 -17.9924 345 3 1.05 0.95;
16 1 329.4 323 0 0 2 0.97673 -16.4894 345 3 1.05 0.95;
17 1 0 0 0 0 1 0.97495 -17.6268 345 1 1.05 0.95;
18 1 158 30 0 0 1 0.9652 -18.6881 345 1 1.05 0.95;
19 1 0 0 0 0 2 1.02958 -11.6791 345 4 1.05 0.95;
20 2 680 103 0 0 2 0.97981 -13.1611 138 4 1.05 0.95;
21 1 274 115 0 0 2 0.99298 -13.9095 345 3 1.05 0.95;
22 1 0 0 0 0 2 1.02902 -9.2532 345 3 1.05 0.95;
23 1 247.5 84.6 0 0 2 1.02329 -9.4586 345 4 1.05 0.95;
24 1 308.6 -92.2 0 0 2 0.98714 -16.3712 345 3 1.05 0.95;
25 2 224 47.2 0 0 1 1.02124 -13.4991 345 1 1.05 0.95;
26 1 139 17 0 0 1 1.01034 -15.2439 345 1 1.05 0.95;
27 1 281 75.5 0 0 1 0.98758 -17.5932 345 1 1.05 0.95;
28 1 206 27.6 0 0 2 1.02015 -11.5383 345 3 1.05 0.95;
29 1 283.5 126.9 0 0 2 1.02428 -8.6481 345 3 1.05 0.95;
30 2 0 0 0 0 1 1.00458 -11.6965 34.5 1 1.05 0.95;
31 2 80 40 0 0 1 0.90644 -9.0711 34.5 2 1.05 0.95;
32 2 0 0 0 0 1 0.97522 -6.1171 21 2 1.05 0.95;
33 2 0 0 0 0 2 0.9972 -6.4145 21 4 1.05 0.95;
34 2 0 0 0 0 2 1.0123 -7.9444 15.5 4 1.05 0.95;
35 2 0 0 0 0 2 1.0493 -4.1911 15.5 3 1.05 0.95;
36 2 0 0 0 0 2 1.0635 -1.4621 12.5 4 1.05 0.95;
37 2 0 0 0 0 1 1.0278 -6.5254 12.5 1 1.05 0.95;
38 2 0 0 0 0 2 1.0265 -1.4812 34.5 3 1.05 0.95;
39 3 400 250 0 0 1 1.03 -10.96 345 2 1.05 0.95;
];
%% 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
mpc.gen = [
30 250 153.18 153.18 -58.084 1.0475 275 1 250 0 0 0 0 0 0 0 0 0 0 0 0;
31 572.93 429.804 429.804 -122.665 1.04 836 1 760 0 0 0 0 0 0 0 0 0 0 0 0;
32 650 446.722 446.722 -180.218 0.9831 843.7 1 767 0 0 0 0 0 0 0 0 0 0 0 0;
33 632 242.94 548.208 -213.78 0.9972 1174.8 1 1068 0 0 0 0 0 0 0 0 0 0 0 0;
34 508 229.094 611.495 -188.028 1.0123 1080.2 1 982 0 0 0 0 0 0 0 0 0 0 0 0;
35 650 361.684 593.788 -234.972 1.0493 1085.7 1 987 0 0 0 0 0 0 0 0 0 0 0 0;
36 560 185.829 568.372 -249.132 1.0635 1025.2 1 932 0 0 0 0 0 0 0 0 0 0 0 0;
37 540 158.09 443.468 -216.122 1.0278 970.2 1 882 0 0 0 0 0 0 0 0 0 0 0 0;
38 830 187.815 834.775 -356.889 1.0265 1684.1 1 1531 0 0 0 0 0 0 0 0 0 0 0 0;
39 720.876 345.605 574.85 -173.261 1.03 1199 1 1090 0 0 0 0 0 0 0 0 0 0 0 0;
10 100 27.003 295.816 -295.816 1.0037 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
20 150 27.003 295.816 -295.816 1.0037 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
2 200 27.003 295.816 -295.816 1.0037 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
25 250 27.003 295.816 -295.816 1.0037 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [
1 2 0.0035 0.0411 0.6987 175 192 192 0 0 1 -360 360;
1 39 0.001 0.025 0.75 125 137 137 0 0 1 -360 360;
2 3 0.0013 0.0151 0.2572 650 715 715 0 0 1 -360 360;
2 25 0.007 0.0086 0.146 260 286 286 0 0 1 -360 360;
3 4 0.0013 0.0213 0.2214 125 137 137 0 0 1 -360 360;
3 18 0.0011 0.0133 0.2138 125 137 137 0 0 1 -360 360;
4 5 0.0008 0.0128 0.1342 300 330 330 0 0 1 -360 360;
4 14 0.0008 0.0129 0.1382 300 330 330 0 0 1 -360 360;
5 6 0.0002 0.0026 0.0434 462 508 508 0 0 1 -360 360;
5 8 0.0008 0.0112 0.1476 462 508 508 0 0 1 -360 360;
6 7 0.0006 0.0092 0.113 700 770 770 0 0 1 -360 360;
6 11 0.0007 0.0082 0.1389 466 512 512 0 0 1 -360 360;
7 8 0.0004 0.0046 0.078 462 508 508 0 0 1 -360 360;
8 9 0.0023 0.0363 0.3804 466 512 512 0 0 1 -360 360;
9 39 0.001 0.025 1.2 400 440 440 0 0 1 -360 360;
10 11 0.0004 0.0043 0.0729 462 508 508 0 0 1 -360 360;
10 13 0.0004 0.0043 0.0729 429 471 471 0 0 1 -360 360;
13 14 0.0009 0.0101 0.1723 380 418 418 0 0 1 -360 360;
14 15 0.0018 0.0217 0.366 125 137 137 0 0 1 -360 360;
15 16 0.0009 0.0094 0.171 429 471 471 0 0 1 -360 360;
16 17 0.0007 0.0089 0.1342 260 286 286 0 0 1 -360 360;
16 19 0.0016 0.0195 0.304 552 607 607 0 0 1 -360 360;
16 21 0.0008 0.0135 0.2548 429 471 471 0 0 1 -360 360;
16 24 0.0003 0.0059 0.068 230 253 253 0 0 1 -360 360;
17 18 0.0007 0.0082 0.1319 300 330 330 0 0 1 -360 360;
17 27 0.0013 0.0173 0.3216 125 137 137 0 0 1 -360 360;
21 22 0.0008 0.014 0.2565 650 715 715 0 0 1 -360 360;
22 23 0.0006 0.0096 0.1846 125 137 137 0 0 1 -360 360;
23 24 0.0022 0.035 0.361 466 512 512 0 0 1 -360 360;
25 26 0.0032 0.0323 0.513 125 137 137 0 0 1 -360 360;
26 27 0.0014 0.0147 0.2396 380 418 418 0 0 1 -360 360;
26 28 0.0043 0.0474 0.7802 175 192 192 0 0 1 -360 360;
26 29 0.0057 0.0625 1.029 230 253 253 0 0 1 -360 360;
28 29 0.0014 0.0151 0.249 400 440 440 0 0 1 -360 360;
2 30 0 0.0181 0 380 418 418 1.025 0 1 -360 360;
31 6 0 0.025 0 700 770 770 0.9 0 1 -360 360;
10 32 0 0.02 0 816 897 897 1.07 0 1 -360 360;
12 11 0.0016 0.0435 0 125 137 137 1.006 0 1 -360 360;
12 13 0.0016 0.0435 0 125 137 137 1.006 0 1 -360 360;
19 20 0.0007 0.0138 0 230 253 253 1.06 0 1 -360 360;
19 33 0.0007 0.0142 0 700 770 770 1.07 0 1 -360 360;
20 34 0.0009 0.018 0 560 616 616 1.009 0 1 -360 360;
22 35 0 0.0143 0 715 786 786 1.025 0 1 -360 360;
23 36 0.0005 0.0272 0 560 616 616 1 0 1 -360 360;
25 37 0.0006 0.0232 0 552 607 607 1.025 0 1 -360 360;
29 38 0.0008 0.0156 0 837 920 920 1.025 0 1 -360 360;
];
%% bus names
mpc.bus_name = {
'BUS1 ';
'BUS2 ';
'LOAD3 ';
'LOAD4 ';
'BUS5 ';
'BUS6 ';
'LOAD7 ';
'LOAD8 ';
'BUS9 ';
'BUS10 ';
'BUS11 ';
'LOAD12 ';
'BUS13 ';
'BUS14 ';
'LOAD15 ';
'LOAD16 ';
'BUS17 ';
'LOAD18 ';
'BUS19 ';
'LOAD20 ';
'LOAD21 ';
'BUS22 ';
'LOAD23 ';
'LOAD24 ';
'LOAD25 ';
'LOAD26 ';
'LOAD27 ';
'LOAD28 ';
'LOAD29 ';
'GEN30 ';
'GEN31 ';
'GEN32 ';
'GEN33 ';
'GEN34 ';
'GEN35 ';
'GEN36 ';
'GEN37 ';
'GEN38 ';
'GEN39 ';
};

80
interface/data/kundur.m
Unescape View File

@ -0,0 +1,80 @@
function mpc = kundur
%KUNDUR
% PSS(R)E 32 RAW created by rawd32 THU, JAN 30 2020 10:51
% MODIFIED KUNDUR'S TWO-AREA TEST SYSTEM, DISTRIBUTED WITH ANDES
% SEE THE BOOK "POWER SYSTEM STABILITY AND CONTROL" FOR ORIGINAL DATA
%
% Converted by MATPOWER 7.1 using PSSE2MPC on 16-Mar-2023
% from 'kundur.raw' using PSS/E rev 32 format.
%
% WARNINGS:
% Skipped 1 line of zone data.
% Skipped 1 line of owner data.
% 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 32.6732 20 1 1.1 0.9;
2 2 0 0 0 0 1 1 21.6548 20 1 1.1 0.9;
3 2 0 0 0 0 2 1 11.2148 20 1 1.1 0.9;
4 2 0 0 0 0 2 1 21.6398 20 1 1.1 0.9;
5 1 0 0 0 0 1 0.98337 27.6488 230 1 1.1 0.9;
6 1 0 0 0 0 1 0.96908 16.8176 230 1 1.1 0.9;
7 1 1159 -73.5 0 0 1 0.95621 8.1662 230 1 1.1 0.9;
8 1 1575 -89.9 0 0 2 0.954 -2.1295 230 1 1.1 0.9;
9 1 0 0 0 0 2 0.96856 6.3774 230 1 1.1 0.9;
10 1 0 0 0 0 2 0.98377 16.8036 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
mpc.gen = [
1 745.861 143.612 600 0 1 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
2 700 300 600 -600 1 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
3 700 550 600 -600 1 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
4 700 -100 600 -600 1 900 1 900 0 0 0 0 0 0 0 0 0 0 0 0;
];
%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [
5 6 0.005 0.05 0.075 0 0 0 0 0 1 -360 360;
5 6 0.00501 0.05001 0.075 0 0 0 0 0 1 -360 360;
6 7 0.002 0.02 0.03 0 0 0 0 0 1 -360 360;
6 7 0.00201 0.02001 0.03 0 0 0 0 0 1 -360 360;
7 8 0.02201 0.22001 0.33 0 0 0 0 0 1 -360 360;
7 8 0.02202 0.22002 0.33 0 0 0 0 0 1 -360 360;
7 8 0.022 0.22 0.33 0 0 0 0 0 1 -360 360;
8 9 0.002 0.02 0.03 0 0 0 0 0 1 -360 360;
8 9 0.00201 0.02001 0.03 0 0 0 0 0 1 -360 360;
9 10 0.005 0.05 0.075 0 0 0 0 0 1 -360 360;
9 10 0.00501 0.05001 0.075 0 0 0 0 0 1 -360 360;
1 5 0.001 0.012 0 0 0 0 1 0 1 -360 360;
2 6 0.001 0.012 0 0 0 0 1 0 1 -360 360;
3 9 0.001 0.012 0 0 0 0 1 0 1 -360 360;
4 10 0.001 0.012 0 0 0 0 1 0 1 -360 360;
];
%% bus names
mpc.bus_name = {
'1 ';
'2 ';
'12 ';
'11 ';
'101 ';
'102 ';
'3 ';
'13 ';
'112 ';
'111 ';
};

7
interface/data/psse2psat.m
Unescape View File

@ -0,0 +1,7 @@
casename = 'ieee14';
%
[mpc, warnings] = psse2mpc(horzcat(casename,'.raw'),horzcat(casename,'.m'));
check = matpower2psat(horzcat(casename,'.m'), pwd);
return

51
interface/ex_cpf.m
Unescape View File

@ -0,0 +1,51 @@
% EXAMPLE: continuation power flow
%
% in Software PowerSAS.m
% Author: Rui Yao <ruiyao@ieee.org>
%
% Copyright (C) 2021, UChicago Argonne, LLC. All rights reserved.
%
% OPEN SOURCE LICENSE
%
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
%
%
% ******************************************************************************************************
% DISCLAIMER
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
% OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% ***************************************************************************************************
clear;
clc;
global Global_Settings
Global_Settings.logLevel='INFO';
res_npcc140_zip=runPowerSAS('pf','d_ieee14');
res_npcc140_zip_cpf=runPowerSAS('cpf','d_ieee14',[],res_npcc140_zip.SysDataBase.zip,[],res_npcc140_zip.snapshot);
plotCurves(1,res_npcc140_zip_cpf.t,res_npcc140_zip_cpf.stateCurve,res_npcc140_zip_cpf.SysDataBase,'v');
%%
res_npcc140_zip=runPowerSAS('pf','d_ieee39');
res_npcc140_zip_cpf=runPowerSAS('cpf','d_ieee39',[],res_npcc140_zip.SysDataBase.zip,[],res_npcc140_zip.snapshot);
plotCurves(1,res_npcc140_zip_cpf.t,res_npcc140_zip_cpf.stateCurve,res_npcc140_zip_cpf.SysDataBase,'v');
%%
res_npcc140_zip=runPowerSAS('pf','d_kundur');
res_npcc140_zip_cpf=runPowerSAS('cpf','d_kundur',[],res_npcc140_zip.SysDataBase.zip,[],res_npcc140_zip.snapshot);
plotCurves(1,res_npcc140_zip_cpf.t,res_npcc140_zip_cpf.stateCurve,res_npcc140_zip_cpf.SysDataBase,'v');
return

38
interface/ex_pf.m
Unescape View File

@ -0,0 +1,38 @@
% EXAMPLE: steady-state
%
% in Software PowerSAS.m
% Author: Rui Yao <ruiyao@ieee.org>
%
% Copyright (C) 2021, UChicago Argonne, LLC. All rights reserved.
%
% OPEN SOURCE LICENSE
%
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
%
%
% ******************************************************************************************************
% DISCLAIMER
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
% OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% ***************************************************************************************************
clear;
clc;
global Global_Settings
Global_Settings.logLevel='INFO';
res_ieee14=runPowerSAS('pf','d_ieee14.m');
res_ieee39=runPowerSAS('pf','d_ieee39.m');
res_kundur=runPowerSAS('pf','d_kundur.m');

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