From bb6526ed86ae090bbaf04432a71f08afc144f9ae Mon Sep 17 00:00:00 2001 From: Yang Liu Date: Mon, 20 Mar 2023 16:25:37 -0400 Subject: [PATCH 1/4] debug the layout of PSS/E interface --- docs/source/data_and_models.rst | 1 + interface/data/d_mpc_temp.m | 79 --------------------------------- 2 files changed, 1 insertion(+), 79 deletions(-) delete mode 100644 interface/data/d_mpc_temp.m diff --git a/docs/source/data_and_models.rst b/docs/source/data_and_models.rst index 328802c..9294c3d 100644 --- a/docs/source/data_and_models.rst +++ b/docs/source/data_and_models.rst @@ -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' diff --git a/interface/data/d_mpc_temp.m b/interface/data/d_mpc_temp.m deleted file mode 100644 index beb7850..0000000 --- a/interface/data/d_mpc_temp.m +++ /dev/null @@ -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'}; - From 16c1547ba961ecfd67780062b964e81760b8f01b Mon Sep 17 00:00:00 2001 From: Yang Liu Date: Tue, 25 Apr 2023 02:58:41 -0400 Subject: [PATCH 2/4] add basic version of EMT module --- EMT/DT.m | 81 +++++++ EMT/DTsolver.m | 368 +++++++++++++++++++++++++++++++ EMT/calEMT.m | 54 +++++ EMT/calculateLambda.m | 7 + EMT/calculateSystemMatrices.m | 68 ++++++ EMT/calculate_L_temp.m | 6 + EMT/convert_to_single_vector.m | 26 +++ EMT/fun1.m | 27 +++ EMT/fun2.m | 42 ++++ EMT/fun3.m | 33 +++ EMT/funEMT_on.m | 99 +++++++++ EMT/funEMT_pre.m | 97 +++++++++ EMT/generator_states.m | 48 ++++ EMT/multiwindow.m | 375 ++++++++++++++++++++++++++++++++ EMT/myFunction.m | 39 ++++ EMT/plotEMTResults.m | 216 ++++++++++++++++++ EMT/setDefaultParameters.m | 9 + EMT/setSASParameters.m | 4 + EMT/set_simulation_parameters.m | 15 ++ EMT/tline_to_matrix.m | 25 +++ EMT/transformParameters.m | 30 +++ EMT/twoarea.m | 82 +++++++ EMT/update_system_parameters.m | 13 ++ example/ex_emt.m | 10 + 24 files changed, 1774 insertions(+) create mode 100644 EMT/DT.m create mode 100644 EMT/DTsolver.m create mode 100644 EMT/calEMT.m create mode 100644 EMT/calculateLambda.m create mode 100644 EMT/calculateSystemMatrices.m create mode 100644 EMT/calculate_L_temp.m create mode 100644 EMT/convert_to_single_vector.m create mode 100644 EMT/fun1.m create mode 100644 EMT/fun2.m create mode 100644 EMT/fun3.m create mode 100644 EMT/funEMT_on.m create mode 100644 EMT/funEMT_pre.m create mode 100644 EMT/generator_states.m create mode 100644 EMT/multiwindow.m create mode 100644 EMT/myFunction.m create mode 100644 EMT/plotEMTResults.m create mode 100644 EMT/setDefaultParameters.m create mode 100644 EMT/setSASParameters.m create mode 100644 EMT/set_simulation_parameters.m create mode 100644 EMT/tline_to_matrix.m create mode 100644 EMT/transformParameters.m create mode 100644 EMT/twoarea.m create mode 100644 EMT/update_system_parameters.m create mode 100644 example/ex_emt.m diff --git a/EMT/DT.m b/EMT/DT.m new file mode 100644 index 0000000..03e7b50 --- /dev/null +++ b/EMT/DT.m @@ -0,0 +1,81 @@ +function [t_new,x_new,y_new,LTE,LTE2,LTE3] = DT(t,x,f,GLO,h,K) +% DT performs one-step integration of differential transformation +% method to solve differential equation dx/dt = f(t,x,GLO); +% +% Input +% k current order, 0,1,2,... +% x current states +% f differential equation +% GLO parameters in the differential equation +% h time step +% K The degree of power series: X(0)+X(1)*h+...X(K)*h^K +% +% Output +% t_new next time instant +% x_new next states of state variables +% y_new next states of algebraic variables +% LTE local truncation error +% +% See also rk4, Parareal +% + + + %% Function handle to compute intermidate variables + fun_y = str2func(strcat(func2str(f),'_AE')); + fun_idx = str2func(strcat(func2str(f),'_idx')); + + %% Function handle for DT model X(k+1)=FB(X(0:k)) of differential equation dx/dt = f(x) + F = str2func(strcat(func2str(f),'_DT')); + + %% compute intermidate variables y + idx = fun_idx(GLO); + y = fun_y(x,idx,GLO); + + %% Define matrix X Y to store coefficients, Xh to store terms + Nx = length(x); + X = zeros(Nx,K+1); %%why NAN + + Ny = length(y); + Y = zeros(Ny,K+1); + + Xh= zeros(Nx,K+1); + Yh= zeros(Ny,K+1); + %% Initialize the 0th order coefficients + X(:,1) = x; + Y(:,1) = y; + + Xh(:,1) = x; % Term(0) = X(0)*h^0 + Yh(:,1) = y; + %% Perform recursive calculation of high order coefficients and terms + for k = 0:K-1 +% for k = 0:K-2 + %% DT Equation: X(k+1)=F(X(0:k)) or XX(k+2)=F(XX(1:k+1)) + [X,Y] = F(k,X,Y,GLO,idx); + Xh(:,k+2) = X(:,k+2)*h^(k+1); % Term(k+1) = X(k+1)*h^(k+1) or XX(k+2)*h^(k+1) + Yh(:,k+2) = Y(:,k+2)*h^(k+1); % Term(k+1) = X(k+1)*h^(k+1) or XX(k+2)*h^(k+1) + + end + + %% sum over all columns of matrix Xh to obtain x_new as a column vector + x_new = sum(Xh(:,1:end-1),2); + + %% +% y_new = sum(Yh,2); + y_new = fun_y(x_new,idx,GLO); + %% Update time + t_new = t + h; + eta = 0.0000000000000000001; + %% Local truncation error: last term; maximum value over all variables + % khuang: we can use relative error, since different variables have + % differnet dynamics and quantities + infcof = max(abs(Xh(:,end)./(x_new+eta))); + trc = infcof^(1/K); + LTE = infcof/(1-trc); +% LTE = infcof; +% LTE = max(abs(Xh(:,end)./x_new)); + LTE3 = max(abs(Xh(:,end-1)./(x_new+eta))); + + % relative error used to change order + LTE2 = infcof/min([norm(Xh(:,end-1)./Xh(:,end),Inf), sqrt(norm(Xh(:,end-2)./Xh(:,end),Inf)), sqrt(norm(Xh(:,end-3)./Xh(:,end-1),Inf))]); +% LTE2 = min(abs(Xh(:,end)./abs(Xh(:,end-1)))) +end \ No newline at end of file diff --git a/EMT/DTsolver.m b/EMT/DTsolver.m new file mode 100644 index 0000000..4cf1ee6 --- /dev/null +++ b/EMT/DTsolver.m @@ -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; + diff --git a/EMT/calEMT.m b/EMT/calEMT.m new file mode 100644 index 0000000..74d45b9 --- /dev/null +++ b/EMT/calEMT.m @@ -0,0 +1,54 @@ + +function [t_final, sol_final] = 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'; +t_final = [t_pre; t_on(2:end); t_post(2:end)]; +sol_final = [x_pre; x_on(2:end,:); x_post(2:end,:)]; +return \ No newline at end of file diff --git a/EMT/calculateLambda.m b/EMT/calculateLambda.m new file mode 100644 index 0000000..c40d9b7 --- /dev/null +++ b/EMT/calculateLambda.m @@ -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; \ No newline at end of file diff --git a/EMT/calculateSystemMatrices.m b/EMT/calculateSystemMatrices.m new file mode 100644 index 0000000..52a4504 --- /dev/null +++ b/EMT/calculateSystemMatrices.m @@ -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 \ No newline at end of file diff --git a/EMT/calculate_L_temp.m b/EMT/calculate_L_temp.m new file mode 100644 index 0000000..6ae9863 --- /dev/null +++ b/EMT/calculate_L_temp.m @@ -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 diff --git a/EMT/convert_to_single_vector.m b/EMT/convert_to_single_vector.m new file mode 100644 index 0000000..6b5b451 --- /dev/null +++ b/EMT/convert_to_single_vector.m @@ -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 diff --git a/EMT/fun1.m b/EMT/fun1.m new file mode 100644 index 0000000..8195701 --- /dev/null +++ b/EMT/fun1.m @@ -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 \ No newline at end of file diff --git a/EMT/fun2.m b/EMT/fun2.m new file mode 100644 index 0000000..1e3ef4b --- /dev/null +++ b/EMT/fun2.m @@ -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 \ No newline at end of file diff --git a/EMT/fun3.m b/EMT/fun3.m new file mode 100644 index 0000000..4e73421 --- /dev/null +++ b/EMT/fun3.m @@ -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)); + + diff --git a/EMT/funEMT_on.m b/EMT/funEMT_on.m new file mode 100644 index 0000000..64be6d7 --- /dev/null +++ b/EMT/funEMT_on.m @@ -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 diff --git a/EMT/funEMT_pre.m b/EMT/funEMT_pre.m new file mode 100644 index 0000000..a1cb799 --- /dev/null +++ b/EMT/funEMT_pre.m @@ -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 diff --git a/EMT/generator_states.m b/EMT/generator_states.m new file mode 100644 index 0000000..dc28ae0 --- /dev/null +++ b/EMT/generator_states.m @@ -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; \ No newline at end of file diff --git a/EMT/multiwindow.m b/EMT/multiwindow.m new file mode 100644 index 0000000..a530355 --- /dev/null +++ b/EMT/multiwindow.m @@ -0,0 +1,375 @@ +function [t,x,y,hh,KK,bx,by,bz] = multiwindow(daeModelEquations,simu_T,x0,GLO,method,solverOptions,bx0) +format long +% multiwindow performs general multi-stage numerical integration. It +% applies to RK4, DT or other use-specified numerical solvers +% +% Inputs +% daeModelEquations the struct for equations related to a DAE model +% simu_T the simulation window: start, step length, end +% x0 initial values of the DAE +% GLO parameters in the DAE model +% method solution method +% solverOptions settings and parameters of the solution method +% Outputs +% t vector of time steps in the whole simulation +% x trajectory of state variables +% y trajectory of algebraic variables +% hh the time step used in each simulation step +% KK the order used in each simulation step +% +% See also Parareal, rk4, DT, getDAEModelEquations, getDAESolver + +f = daeModelEquations.fun_DAE; +fun_idx = daeModelEquations.fun_idx; + +%% +K = solverOptions.K; +% FS1_VS0 = solverOptions.FS1_VS0; +% FO1_VO0 = solverOptions.FO1_VO0; +hmax = solverOptions.hmax; +hmin = solverOptions.hmin; +Kmax = solverOptions.Kmax; +Kmin = solverOptions.Kmin; +% eps1 = solverOptions.eps1; +% eps2 = solverOptions.eps2; +% eps3 = solverOptions.eps3; +% eps4 = solverOptions.eps4; +tol = solverOptions.RBtol; +% q1 = solverOptions.q1; +% q2 = solverOptions.q2; +% dK = solverOptions.dK; +%% +a = simu_T(1); % start time of simulation +b = simu_T(2); % end time of simulation +h = simu_T(3); % time step length + +%% Judge if x0 is column vector or not. If rox vector, convert to column vector. +m = size(x0,1); % number of rows +if m == 1 + x0 = x0'; + m =size(x0,1); +end +bm = size(bx0,1); % number of rows +if bm == 1 + bx0 = bx0'; +end +%% +fun_y = str2func(strcat(func2str(f),'_AE')); +idx = fun_idx(GLO); +y0 = fun_y(x0,idx,GLO); +by0 = fun_y(bx0,idx,GLO); +%% Initial values of t and x +% mite = ceil((b-a)./hmin); +mite = 20000; +bt = zeros(1,mite); +bx = zeros(m,mite); +by = zeros(size(y0,1),mite); +t = zeros(1,mite); +x = zeros(m,mite); +y = zeros(size(y0,1),mite); +hh = zeros(1,mite); +KK = zeros(1,mite); +t(1) = a; +x(:,1) = x0; %initial conditions +bt(:,1) = a; +bx(:,1) = bx0; + +y(:,1) = y0; +by(:,1) = by0; +hh(:,1)= h; +KK(:,1)= K; + +%% Implement the integration method +ratio_b4= nan; +theta_max=2; + + + +ite = 1 ; + +% K = -(1/2)*log(tol); +M = 1; +fac1= 0.98; +fac2 = 1; + +% fac1 = 0; +% fac2 = 0; +% + +p= 1 ; +pf1 = 2*GLO.machine; +pf2 = 2*GLO.machine+4*GLO.machine*GLO.machine+idx.Nx*(idx.Nx+idx.Ny)+(idx.Nx); +pf3 = (idx.Nx); +% pf1 = 17*GLO.machine/2; +% pf2 = 27*GLO.machine/2+4*GLO.machine*GLO.machine+idx.Nx*(idx.Nx+idx.Ny)+(idx.Nx); +% pf2 = 0; +% pf3 = 0; +while t(ite)+h <= b + +switch solverOptions.solver + case 'DT' + + [t_new,x_new,y_new, LTE,LTE2,LTE3] = method(t(ite),x(:,ite),f,GLO,h,K); + ratio = LTE/h; + K_I = 0.3/K; + K_p = 0.4/K; + if isnan(ratio_b4) + else + theta = ((tol/ratio)^(K_I))*((ratio_b4/ratio)^(K_p)); + + if mod(ite,M)==0 +% Lin = ((K+p)/(K)) +% Lin = ((K+p)^2/(K)^2) + Lin = (pf1*(K+p)^2+pf2*(K+p)+pf3)/(pf1*(K)^2+pf2*(K)+pf3); +% Lde = ((K-p)/(K)) +% Lde = ((K-p)^2/(K)^2) + Lde = (pf1*(K-p)^2+pf2*(K-p)+pf3)/(pf1*(K)^2+pf2*(K)+pf3); + +% trc_in =(LTE2); +% trc_de= (LTE3); + + trc_in =(LTE2)/(1-LTE2^((K+p))); + trc_de= (LTE3)/(1-LTE3^((K-p))); + + tol_in = tol; + tol_de = tol; + + theta_in =1*(h*tol_in/trc_in)^(1/(K+p)); + theta_de =1*(h*tol_de/trc_de)^(1/(K-p)); + +% +% theta_in =1*(tol_in/trc_in)^(1/(K+p)); +% theta_de =1*(tol_de/trc_de)^(1/(K-p)); + +% % +% % theta_in =0.9*(h*tol/trc_in)^(1/(K+p)); +% % theta_de =0.9*(h*tol/trc_de)^(1/(K-p)); + + if theta> theta_max + h_can = theta_max*h; + else + h_can = theta*h; + end + + + if h_can > hmax + h_can = hmax; + + elseif h < hmin + h_can = hmin; + + end + + if theta_in> theta_max + theta_in = theta_max; + end + if theta_de> theta_max + theta_de = theta_max; + end + + + + if h*theta_in > hmax + h_inp = hmax; + + elseif h*theta_in < hmin + h_inp = hmin; + else + + h_inp = h*theta_in; + end + + + if h*theta_de > hmax + h_dep = hmax; + + elseif h*theta_de < hmin + h_dep = hmin; + else + h_dep = h*theta_de; + end + + + + + + + + + +% if min(hh(end-1:-1:end-M))>= h*theta&&Lin = h_can||KK(end-1)K) +% K=K-p; +% theta =theta_de +% +% end + + + if (Lde=K)))&&((h>= h_can||(h==hmax&&KK(1,ite+1)<=K))&&Lin =K)) + K=K-p; + theta =theta_de; + tol = tol_de; + + elseif (h>= h_can||(h_can==hmax&&KK(1,ite+1)<=K))&&Lin theta_max + h = theta_max*h; + else + h = theta*h; + end + + + end + ratio_b4 = ratio; + + + + case 'RK4' + [t_new,x_new,y_new] = method(t(ite),x(:,ite),f,GLO,h,K); + +end + + + +% bhx = bx(:,ite); +% for bht =t(ite):0.00031: t_new(end) +% if t_new(end)-bht>=0.00031 +% [bht_new,bhx_new,bhy_new] = rk4(bht,bhx,f,GLO,0.00031,8); +% bhx = bhx_new(:,end); +% else +% [bht_new,bhx_new,bhy_new] = rk4(bht,bhx,f,GLO,t_new(end)-bht,8); +% break +% % bht = t_new(end); +% end +% end + + if ite+1 > mite + mite = mite + 20000; +% bt(:,mite) = 0; +% bx(:,mite) = 0; +% by(:,mite) = 0; + t(:,mite) = 0; + x(:,mite) = 0; + y(:,mite) = 0; + hh(:,mite) = 0; + KK(:,mite) = 0; + end + + + + if h > hmax + h = hmax; + + elseif h < hmin + h = hmin; + + end + + if K > Kmax; K = Kmax; end + if K < Kmin; K = Kmin; end +% bt(:,ite+1) = bht_new(end); +% bx(:,ite+1) = bhx_new(:,end); +% by(:,ite+1) = bhy_new(:,end); + t(:,ite+1) = t_new; + x(:,ite+1) = x_new; + y(:,ite+1) = y_new; + hh(:,ite+1) = h; + KK(:,ite+1) = K; + +% r = tol^(1/(1+K)); + +% theta =(tol*h/r_norm)^(1/K); + + +% K; + ite = ite+1; + +% if FS1_VS0 == 0 +% if LTE < eps1; h = h*q1; end +% if LTE > eps2; h = h*q2; end +% if h > hmax; h = hmax; end +% if h < hmin; h = hmin; end +% end +% +% if FO1_VO0 == 0 +% if LTE < eps3 +% K = K-dK; +% end +% if LTE > eps4 +% K = K+dK; +% end +% if K > Kmax; K = Kmax; end +% if K < Kmin; K = Kmin; end +% end + if b-t(ite)0 + h = b-t(ite); + end + +end + +%% transpose: use .' to avoid conjugate transpose in case some variables are complex +t = t(:,1:ite).'; +x = x(:,1:ite).'; +y = y(:,1:ite).'; +hh = hh(:,1:ite).'; +KK = KK(:,1:ite).'; +% bx = bx(:,1:ite).'; +% by = by(:,1:ite).'; +% bz =max(abs(x-bx),[],2); +bx = x; +by = y; +bz =max(abs(x-bx),[],2); +% bz =max((x-bx)./(max(bx+10^(-19),[],1)),[],2); +end \ No newline at end of file diff --git a/EMT/myFunction.m b/EMT/myFunction.m new file mode 100644 index 0000000..7c690d1 --- /dev/null +++ b/EMT/myFunction.m @@ -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 + diff --git a/EMT/plotEMTResults.m b/EMT/plotEMTResults.m new file mode 100644 index 0000000..13b77eb --- /dev/null +++ b/EMT/plotEMTResults.m @@ -0,0 +1,216 @@ +close all + +% t_final = [t_pre; t_on(2:end); t_post(2:end)]; +% sol_final = [x_pre; x_on(2:end,:); x_post(2:end,:)]; +% % t_final = [t_pre; t_on(2:end)]; +% % sol_final = [x_pre; x_on(2:end,:)]; +% % t_final = [t_pre]; +% % sol_final = [x_pre]; + +%% 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') diff --git a/EMT/setDefaultParameters.m b/EMT/setDefaultParameters.m new file mode 100644 index 0000000..69ffe69 --- /dev/null +++ b/EMT/setDefaultParameters.m @@ -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 diff --git a/EMT/setSASParameters.m b/EMT/setSASParameters.m new file mode 100644 index 0000000..7f2927f --- /dev/null +++ b/EMT/setSASParameters.m @@ -0,0 +1,4 @@ +function [orderOfDT, timeStepOfDT] = setSASParameters() +orderOfDT=10; +timeStepOfDT=179*10^(-6); +end \ No newline at end of file diff --git a/EMT/set_simulation_parameters.m b/EMT/set_simulation_parameters.m new file mode 100644 index 0000000..454fa6f --- /dev/null +++ b/EMT/set_simulation_parameters.m @@ -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 diff --git a/EMT/tline_to_matrix.m b/EMT/tline_to_matrix.m new file mode 100644 index 0000000..af7c403 --- /dev/null +++ b/EMT/tline_to_matrix.m @@ -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 diff --git a/EMT/transformParameters.m b/EMT/transformParameters.m new file mode 100644 index 0000000..4034049 --- /dev/null +++ b/EMT/transformParameters.m @@ -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 diff --git a/EMT/twoarea.m b/EMT/twoarea.m new file mode 100644 index 0000000..f285138 --- /dev/null +++ b/EMT/twoarea.m @@ -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 '; +}; diff --git a/EMT/update_system_parameters.m b/EMT/update_system_parameters.m new file mode 100644 index 0000000..8b0e326 --- /dev/null +++ b/EMT/update_system_parameters.m @@ -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 + diff --git a/example/ex_emt.m b/example/ex_emt.m new file mode 100644 index 0000000..d64f6d9 --- /dev/null +++ b/example/ex_emt.m @@ -0,0 +1,10 @@ + +clc; clear; close; +addpath('EMT') + +testSystem = 'twoarea'; + +[t_final, sol_final] = calEMT(testSystem); + +plotEMTResults(); + From e3a4af8e2c1e881e3ea8d39071ae38b42f1de958 Mon Sep 17 00:00:00 2001 From: Yang Liu Date: Tue, 25 Apr 2023 03:01:42 -0400 Subject: [PATCH 3/4] delete wrong files --- EMT/DT.m | 81 ---------- EMT/multiwindow.m | 375 ---------------------------------------------- 2 files changed, 456 deletions(-) delete mode 100644 EMT/DT.m delete mode 100644 EMT/multiwindow.m diff --git a/EMT/DT.m b/EMT/DT.m deleted file mode 100644 index 03e7b50..0000000 --- a/EMT/DT.m +++ /dev/null @@ -1,81 +0,0 @@ -function [t_new,x_new,y_new,LTE,LTE2,LTE3] = DT(t,x,f,GLO,h,K) -% DT performs one-step integration of differential transformation -% method to solve differential equation dx/dt = f(t,x,GLO); -% -% Input -% k current order, 0,1,2,... -% x current states -% f differential equation -% GLO parameters in the differential equation -% h time step -% K The degree of power series: X(0)+X(1)*h+...X(K)*h^K -% -% Output -% t_new next time instant -% x_new next states of state variables -% y_new next states of algebraic variables -% LTE local truncation error -% -% See also rk4, Parareal -% - - - %% Function handle to compute intermidate variables - fun_y = str2func(strcat(func2str(f),'_AE')); - fun_idx = str2func(strcat(func2str(f),'_idx')); - - %% Function handle for DT model X(k+1)=FB(X(0:k)) of differential equation dx/dt = f(x) - F = str2func(strcat(func2str(f),'_DT')); - - %% compute intermidate variables y - idx = fun_idx(GLO); - y = fun_y(x,idx,GLO); - - %% Define matrix X Y to store coefficients, Xh to store terms - Nx = length(x); - X = zeros(Nx,K+1); %%why NAN - - Ny = length(y); - Y = zeros(Ny,K+1); - - Xh= zeros(Nx,K+1); - Yh= zeros(Ny,K+1); - %% Initialize the 0th order coefficients - X(:,1) = x; - Y(:,1) = y; - - Xh(:,1) = x; % Term(0) = X(0)*h^0 - Yh(:,1) = y; - %% Perform recursive calculation of high order coefficients and terms - for k = 0:K-1 -% for k = 0:K-2 - %% DT Equation: X(k+1)=F(X(0:k)) or XX(k+2)=F(XX(1:k+1)) - [X,Y] = F(k,X,Y,GLO,idx); - Xh(:,k+2) = X(:,k+2)*h^(k+1); % Term(k+1) = X(k+1)*h^(k+1) or XX(k+2)*h^(k+1) - Yh(:,k+2) = Y(:,k+2)*h^(k+1); % Term(k+1) = X(k+1)*h^(k+1) or XX(k+2)*h^(k+1) - - end - - %% sum over all columns of matrix Xh to obtain x_new as a column vector - x_new = sum(Xh(:,1:end-1),2); - - %% -% y_new = sum(Yh,2); - y_new = fun_y(x_new,idx,GLO); - %% Update time - t_new = t + h; - eta = 0.0000000000000000001; - %% Local truncation error: last term; maximum value over all variables - % khuang: we can use relative error, since different variables have - % differnet dynamics and quantities - infcof = max(abs(Xh(:,end)./(x_new+eta))); - trc = infcof^(1/K); - LTE = infcof/(1-trc); -% LTE = infcof; -% LTE = max(abs(Xh(:,end)./x_new)); - LTE3 = max(abs(Xh(:,end-1)./(x_new+eta))); - - % relative error used to change order - LTE2 = infcof/min([norm(Xh(:,end-1)./Xh(:,end),Inf), sqrt(norm(Xh(:,end-2)./Xh(:,end),Inf)), sqrt(norm(Xh(:,end-3)./Xh(:,end-1),Inf))]); -% LTE2 = min(abs(Xh(:,end)./abs(Xh(:,end-1)))) -end \ No newline at end of file diff --git a/EMT/multiwindow.m b/EMT/multiwindow.m deleted file mode 100644 index a530355..0000000 --- a/EMT/multiwindow.m +++ /dev/null @@ -1,375 +0,0 @@ -function [t,x,y,hh,KK,bx,by,bz] = multiwindow(daeModelEquations,simu_T,x0,GLO,method,solverOptions,bx0) -format long -% multiwindow performs general multi-stage numerical integration. It -% applies to RK4, DT or other use-specified numerical solvers -% -% Inputs -% daeModelEquations the struct for equations related to a DAE model -% simu_T the simulation window: start, step length, end -% x0 initial values of the DAE -% GLO parameters in the DAE model -% method solution method -% solverOptions settings and parameters of the solution method -% Outputs -% t vector of time steps in the whole simulation -% x trajectory of state variables -% y trajectory of algebraic variables -% hh the time step used in each simulation step -% KK the order used in each simulation step -% -% See also Parareal, rk4, DT, getDAEModelEquations, getDAESolver - -f = daeModelEquations.fun_DAE; -fun_idx = daeModelEquations.fun_idx; - -%% -K = solverOptions.K; -% FS1_VS0 = solverOptions.FS1_VS0; -% FO1_VO0 = solverOptions.FO1_VO0; -hmax = solverOptions.hmax; -hmin = solverOptions.hmin; -Kmax = solverOptions.Kmax; -Kmin = solverOptions.Kmin; -% eps1 = solverOptions.eps1; -% eps2 = solverOptions.eps2; -% eps3 = solverOptions.eps3; -% eps4 = solverOptions.eps4; -tol = solverOptions.RBtol; -% q1 = solverOptions.q1; -% q2 = solverOptions.q2; -% dK = solverOptions.dK; -%% -a = simu_T(1); % start time of simulation -b = simu_T(2); % end time of simulation -h = simu_T(3); % time step length - -%% Judge if x0 is column vector or not. If rox vector, convert to column vector. -m = size(x0,1); % number of rows -if m == 1 - x0 = x0'; - m =size(x0,1); -end -bm = size(bx0,1); % number of rows -if bm == 1 - bx0 = bx0'; -end -%% -fun_y = str2func(strcat(func2str(f),'_AE')); -idx = fun_idx(GLO); -y0 = fun_y(x0,idx,GLO); -by0 = fun_y(bx0,idx,GLO); -%% Initial values of t and x -% mite = ceil((b-a)./hmin); -mite = 20000; -bt = zeros(1,mite); -bx = zeros(m,mite); -by = zeros(size(y0,1),mite); -t = zeros(1,mite); -x = zeros(m,mite); -y = zeros(size(y0,1),mite); -hh = zeros(1,mite); -KK = zeros(1,mite); -t(1) = a; -x(:,1) = x0; %initial conditions -bt(:,1) = a; -bx(:,1) = bx0; - -y(:,1) = y0; -by(:,1) = by0; -hh(:,1)= h; -KK(:,1)= K; - -%% Implement the integration method -ratio_b4= nan; -theta_max=2; - - - -ite = 1 ; - -% K = -(1/2)*log(tol); -M = 1; -fac1= 0.98; -fac2 = 1; - -% fac1 = 0; -% fac2 = 0; -% - -p= 1 ; -pf1 = 2*GLO.machine; -pf2 = 2*GLO.machine+4*GLO.machine*GLO.machine+idx.Nx*(idx.Nx+idx.Ny)+(idx.Nx); -pf3 = (idx.Nx); -% pf1 = 17*GLO.machine/2; -% pf2 = 27*GLO.machine/2+4*GLO.machine*GLO.machine+idx.Nx*(idx.Nx+idx.Ny)+(idx.Nx); -% pf2 = 0; -% pf3 = 0; -while t(ite)+h <= b - -switch solverOptions.solver - case 'DT' - - [t_new,x_new,y_new, LTE,LTE2,LTE3] = method(t(ite),x(:,ite),f,GLO,h,K); - ratio = LTE/h; - K_I = 0.3/K; - K_p = 0.4/K; - if isnan(ratio_b4) - else - theta = ((tol/ratio)^(K_I))*((ratio_b4/ratio)^(K_p)); - - if mod(ite,M)==0 -% Lin = ((K+p)/(K)) -% Lin = ((K+p)^2/(K)^2) - Lin = (pf1*(K+p)^2+pf2*(K+p)+pf3)/(pf1*(K)^2+pf2*(K)+pf3); -% Lde = ((K-p)/(K)) -% Lde = ((K-p)^2/(K)^2) - Lde = (pf1*(K-p)^2+pf2*(K-p)+pf3)/(pf1*(K)^2+pf2*(K)+pf3); - -% trc_in =(LTE2); -% trc_de= (LTE3); - - trc_in =(LTE2)/(1-LTE2^((K+p))); - trc_de= (LTE3)/(1-LTE3^((K-p))); - - tol_in = tol; - tol_de = tol; - - theta_in =1*(h*tol_in/trc_in)^(1/(K+p)); - theta_de =1*(h*tol_de/trc_de)^(1/(K-p)); - -% -% theta_in =1*(tol_in/trc_in)^(1/(K+p)); -% theta_de =1*(tol_de/trc_de)^(1/(K-p)); - -% % -% % theta_in =0.9*(h*tol/trc_in)^(1/(K+p)); -% % theta_de =0.9*(h*tol/trc_de)^(1/(K-p)); - - if theta> theta_max - h_can = theta_max*h; - else - h_can = theta*h; - end - - - if h_can > hmax - h_can = hmax; - - elseif h < hmin - h_can = hmin; - - end - - if theta_in> theta_max - theta_in = theta_max; - end - if theta_de> theta_max - theta_de = theta_max; - end - - - - if h*theta_in > hmax - h_inp = hmax; - - elseif h*theta_in < hmin - h_inp = hmin; - else - - h_inp = h*theta_in; - end - - - if h*theta_de > hmax - h_dep = hmax; - - elseif h*theta_de < hmin - h_dep = hmin; - else - h_dep = h*theta_de; - end - - - - - - - - - -% if min(hh(end-1:-1:end-M))>= h*theta&&Lin = h_can||KK(end-1)K) -% K=K-p; -% theta =theta_de -% -% end - - - if (Lde=K)))&&((h>= h_can||(h==hmax&&KK(1,ite+1)<=K))&&Lin =K)) - K=K-p; - theta =theta_de; - tol = tol_de; - - elseif (h>= h_can||(h_can==hmax&&KK(1,ite+1)<=K))&&Lin theta_max - h = theta_max*h; - else - h = theta*h; - end - - - end - ratio_b4 = ratio; - - - - case 'RK4' - [t_new,x_new,y_new] = method(t(ite),x(:,ite),f,GLO,h,K); - -end - - - -% bhx = bx(:,ite); -% for bht =t(ite):0.00031: t_new(end) -% if t_new(end)-bht>=0.00031 -% [bht_new,bhx_new,bhy_new] = rk4(bht,bhx,f,GLO,0.00031,8); -% bhx = bhx_new(:,end); -% else -% [bht_new,bhx_new,bhy_new] = rk4(bht,bhx,f,GLO,t_new(end)-bht,8); -% break -% % bht = t_new(end); -% end -% end - - if ite+1 > mite - mite = mite + 20000; -% bt(:,mite) = 0; -% bx(:,mite) = 0; -% by(:,mite) = 0; - t(:,mite) = 0; - x(:,mite) = 0; - y(:,mite) = 0; - hh(:,mite) = 0; - KK(:,mite) = 0; - end - - - - if h > hmax - h = hmax; - - elseif h < hmin - h = hmin; - - end - - if K > Kmax; K = Kmax; end - if K < Kmin; K = Kmin; end -% bt(:,ite+1) = bht_new(end); -% bx(:,ite+1) = bhx_new(:,end); -% by(:,ite+1) = bhy_new(:,end); - t(:,ite+1) = t_new; - x(:,ite+1) = x_new; - y(:,ite+1) = y_new; - hh(:,ite+1) = h; - KK(:,ite+1) = K; - -% r = tol^(1/(1+K)); - -% theta =(tol*h/r_norm)^(1/K); - - -% K; - ite = ite+1; - -% if FS1_VS0 == 0 -% if LTE < eps1; h = h*q1; end -% if LTE > eps2; h = h*q2; end -% if h > hmax; h = hmax; end -% if h < hmin; h = hmin; end -% end -% -% if FO1_VO0 == 0 -% if LTE < eps3 -% K = K-dK; -% end -% if LTE > eps4 -% K = K+dK; -% end -% if K > Kmax; K = Kmax; end -% if K < Kmin; K = Kmin; end -% end - if b-t(ite)0 - h = b-t(ite); - end - -end - -%% transpose: use .' to avoid conjugate transpose in case some variables are complex -t = t(:,1:ite).'; -x = x(:,1:ite).'; -y = y(:,1:ite).'; -hh = hh(:,1:ite).'; -KK = KK(:,1:ite).'; -% bx = bx(:,1:ite).'; -% by = by(:,1:ite).'; -% bz =max(abs(x-bx),[],2); -bx = x; -by = y; -bz =max(abs(x-bx),[],2); -% bz =max((x-bx)./(max(bx+10^(-19),[],1)),[],2); -end \ No newline at end of file From 68c835ca1e79337727c3104718b1e4285595f3ce Mon Sep 17 00:00:00 2001 From: Yang Liu Date: Tue, 25 Apr 2023 03:22:42 -0400 Subject: [PATCH 4/4] enable emt in runPowerSAS.m function --- EMT/calEMT.m | 6 +- EMT/plotEMTResults.m | 362 ++++++++++++++++---------------- example/ex_emt.m | 7 +- example/output/~default_log.txt | 6 + runPowerSAS.m | 2 + 5 files changed, 192 insertions(+), 191 deletions(-) create mode 100644 example/output/~default_log.txt diff --git a/EMT/calEMT.m b/EMT/calEMT.m index 74d45b9..270afb6 100644 --- a/EMT/calEMT.m +++ b/EMT/calEMT.m @@ -1,5 +1,5 @@ -function [t_final, sol_final] = calEMT(testSystem) +function [res] = calEMT(testSystem) %% [orderOfDT, timeStepOfDT] = setSASParameters(); @@ -49,6 +49,6 @@ x_ini_post = x_on(end,:)'; t_post = t_dt; x_post = x_dt'; -t_final = [t_pre; t_on(2:end); t_post(2:end)]; -sol_final = [x_pre; x_on(2:end,:); x_post(2:end,:)]; +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 \ No newline at end of file diff --git a/EMT/plotEMTResults.m b/EMT/plotEMTResults.m index 13b77eb..aa29153 100644 --- a/EMT/plotEMTResults.m +++ b/EMT/plotEMTResults.m @@ -1,97 +1,93 @@ +function plotEMTResults(res) close all - -% t_final = [t_pre; t_on(2:end); t_post(2:end)]; -% sol_final = [x_pre; x_on(2:end,:); x_post(2:end,:)]; -% % t_final = [t_pre; t_on(2:end)]; -% % sol_final = [x_pre; x_on(2:end,:)]; -% % t_final = [t_pre]; -% % sol_final = [x_pre]; - -%% 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') +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) @@ -108,50 +104,50 @@ 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') +% %% 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) @@ -168,49 +164,49 @@ 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') +% %% 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') diff --git a/example/ex_emt.m b/example/ex_emt.m index d64f6d9..1b50696 100644 --- a/example/ex_emt.m +++ b/example/ex_emt.m @@ -1,10 +1,7 @@ - clc; clear; close; addpath('EMT') -testSystem = 'twoarea'; +res = runPowerSAS('emt','twoarea'); -[t_final, sol_final] = calEMT(testSystem); - -plotEMTResults(); +plotEMTResults(res); diff --git a/example/output/~default_log.txt b/example/output/~default_log.txt new file mode 100644 index 0000000..2da6ae8 --- /dev/null +++ b/example/output/~default_log.txt @@ -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] + diff --git a/runPowerSAS.m b/runPowerSAS.m index 80a88ef..a39652d 100644 --- a/runPowerSAS.m +++ b/runPowerSAS.m @@ -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')