diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..dfe0770 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,2 @@ +# Auto detect text files and perform LF normalization +* text=auto diff --git a/Australian14gen_following_inertia.slx b/Australian14gen_following_inertia.slx new file mode 100644 index 0000000..1e1457a Binary files /dev/null and b/Australian14gen_following_inertia.slx differ diff --git a/Australian14gen_forming_inertia.slx b/Australian14gen_forming_inertia.slx new file mode 100644 index 0000000..76d1472 Binary files /dev/null and b/Australian14gen_forming_inertia.slx differ diff --git a/Australian14gen_original.slx b/Australian14gen_original.slx new file mode 100644 index 0000000..eeb69e5 Binary files /dev/null and b/Australian14gen_original.slx differ diff --git a/Library.slx b/Library.slx new file mode 100644 index 0000000..e335f32 Binary files /dev/null and b/Library.slx differ diff --git a/ReadMe.txt b/ReadMe.txt new file mode 100644 index 0000000..f819fde --- /dev/null +++ b/ReadMe.txt @@ -0,0 +1,33 @@ + +This repository contains the SIMULINK models and MATLAB codes used as a benchmark system in the paper: +B.K.Poolla, D. Groß, and F. Dörfler, "Placement and Implementation of Grid-Forming and Grid-Following Virtual Inertia" + +dataaustralian14gen.m generates the underlying data, load model, for the grid-following and grid-forming Virtual Inertia (VI) implementations. +Australian14gen_original.slx contains the low-inertia model without any VI devices. +Library.slx contains the custom models for the VI devices in the different implementations. + +Australian14gen_following_inertia.slx contains the low-inertia model of the South-East Australian grid along with grid-following VI devices. +main_following.m contains the relevant code for the non-linear simulation of the grid-following VI implementation. + +Australian14gen_forming_inertia.slx contains the low-inertia model of the South-East Australian grid along with grid-forming VI devices. +main_forming.m contains the relevant code for the non-linear simulation of the grid-following VI implementation. + +The scripts 'main_following.m' and 'main_forming.m' call the following functions: +linearizesystem.m generates the linearized system models (A,B,C,D,G matrices) for analysis in either implementation. +modelreduction.m removes the zero eigenvalue due to the system Laplacian. +tdsim.m simulates the linearized models for the VI implementations. +plotvi.m generates the relevant time-domain plots for the metrics of interest. + + +The models are based original SIMULINK model available at https://ch.mathworks.com/matlabcentral/fileexchange/51177-australian-simplified-14-generators-ieee-benchmark + +% This source code is distributed in the hope that it will be useful, but without any warranty. +% We do request that publications in which this testbed is adopted, explicitly acknowledge that fact by citing the above mentioned paper. + + + + + + + + diff --git a/dataaustralian14gen.m b/dataaustralian14gen.m new file mode 100644 index 0000000..6a3345b --- /dev/null +++ b/dataaustralian14gen.m @@ -0,0 +1,333 @@ +PSSMODEL=1; +%bus p(MW) q (MVar) +load_case4=[... +102 270 30 +205 235 25 +206 80 10 +207 1130 120 +208 125 15 +211 1060 110 +212 1000 110 +215 290 30 +216 1105 120 +217 750 80 +306 900 90 +307 470 50 +308 620 100 +309 140 15 +312 92 10 +313 1625 165 +314 180 20 +405 730 75 +406 540 55 +407 0 0 +408 110 10 +409 190 20 +410 390 40 +411 420 45 +412 922 100 +504 180 20 +507 640 65 +508 490 50 +509 122 15]; +load_case4(20,2)=1; +load_mag=load_case4; + +%bus MW(each unit) KV H Xa Xd Xq Xd' Tdo' Xd" Tdo" Xq' Tqo' Xq" Tqo" r +%1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 +G_data=[... +101 333.3 15 3.6 0.14 1.1 0.65 0.25 8.5 0.25 0.05 0 0 0.25 0.2 0 +201 666.7 20 3.2 0.2 1.8 1.75 0.3 8.5 0.21 0.04 0.7 2 0.21 0.25 0 +202 555.6 20 2.8 0.17 2.2 2.1 0.3 4.5 0.2 0.04 0.5 1.5 0.21 0.06 0 +203 555.6 20 2.6 0.2 2.3 1.7 0.3 5 0.25 0.03 0.4 2 0.25 0.25 0 +204 666.7 20 3.2 0.2 1.8 1.75 0.3 8.5 0.21 0.04 0.7 2 0.21 0.25 0 +301 666.7 20 2.8 0.2 2.7 1.5 0.3 7.5 0.25 0.04 0.85 0.85 0.25 0.12 0 +302 444.4 20 3.5 0.15 2 1.8 0.25 7.5 0.2 0.04 0 0 0.2 0.25 0 +401 444.4 20 2.6 0.2 2.3 1.7 0.3 5 0.25 0.03 0.4 2 0.25 0.25 0 +402 333.3 20 3 0.2 1.9 1.8 0.3 6.5 0.26 0.035 0.55 1.4 0.26 0.04 0 +403 444.4 20 2.6 0.2 2.3 1.7 0.3 5 0.25 0.03 0.4 2 0.25 0.25 0 +404 333.3 20 4 0.18 2.2 1.4 0.32 9 0.24 0.04 0.75 1.4 0.24 0.13 0 +501 333.3 20 3.5 0.15 2.2 1.7 0.3 7.5 0.24 0.025 0.8 1.5 0.24 0.1 0 +502 250 15 4 0.2 2 1.5 0.3 7.5 0.22 0.04 0.8 3 0.22 0.2 0 +503 166.7 15 7.5 0.15 2.3 2 0.25 5 0.17 0.022 0.35 1 0.17 0.035 0]; +G_data(1,13)=1; +G_data(7,13)=1; +active_gen_case4=[2 4 3 3 4 6 2 3 2 3 3 2 3 1]; +G_data(:,2)=G_data(:,2).*active_gen_case4'; + +gentestcase=[2:8,11,12,14]; +genrem=setdiff(1:14,gentestcase); +Hinv=1/2*diag(G_data(gentestcase,4))^-1; + + +% 1 2 3 4 5 6 7 8 9 10 11 12 +% bus Tr KA TA TB TC KE TE KF TF TB1(s) TC1(s) +Exc=[... +101 0 200 0.1 13.25 2.5 0 0 0 0 0 0 +%201 0 300 0.01 0.7 0.35 0 0 0 0 0 0 +201 0 400 0.02 1.12 0.5 0 0 0 0 0 0 +202 0 400 0.02 0 0 1 1 0.029 1 0 0 +203 0 300 0.01 0.7 0.35 0 0 0 0 0 0 +%204 0 300 0.01 0.7 0.35 0 0 0 0 0 0 +204 0 400 0.02 1.12 0.5 0 0 0 0 0 0 +301 0 400 0.05 6.42 1.14 0 0 0 0 0 0 +302 0 200 0.05 0 0 1 1.333 0.02 0.8 0 0 +401 0 300 0.1 40 4 0 0 0 0 0 0 +402 0.02 300 0.05 9.8 1.52 0 0 0 0 0 0 +403 0 300 0.01 0.7 0.35 0 0 0 0 0 0 +404 0 250 0.2 0.0232 0.136 0 0 0 0 0 0 +501 0 1000 0.04 0 0 1 0.87 0.004 0.27 0 0 +502 0 400 0.5 16 1.4 0 0 0 0 0.05 0.6 +503 0 300 0.01 0.8 0.2 0 0 0 0 0 0 +]; + Exc(:,2)=0.02; + Exc(13,9)=.001; + Exc(13,10)=1; +% bus rating pg (MW) Heavy load case + + +% bus MVA(per unit) p(per unit) q(per unit) active units +p_case4=[... +101 333.3 .001 -97.4 2 +201 666.7 540 -30.8 4 +202 555.6 460 -2.5 3 +203 555.6 470 9.4 3 +204 666.7 399.3 -43.6 4 +301 666.7 555 16.6 6 +302 444.4 380 -9.3 2 +401 444.4 320 -21.9 3 +402 333.3 290 -2.4 2 +403 444.4 320 14.2 3 +404 333.3 217 -3.5 3 +501 333.3 280 -52.5 2 +502 250 180 -1.8 3 +503 166.7 150 2.2 1]; +p=p_case4; +p(:,3)=p(:,3).*p(:,5); +p(:,4)=p(:,4).*p(:,5); +pref=p(:,3)./(p(:,2).*p(:,5)); + +bus_v=[... +101 15 +102 330 +201 20 +202 20 +203 20 +204 20 +205 330 +206 330 +207 330 +208 330 +209 330 +210 500 +211 330 +212 330 +213 500 +214 330 +215 330 +216 330 +217 330 +301 20 +302 20 +303 500 +304 500 +305 500 +306 500 +307 500 +308 500 +309 330 +310 330 +311 330 +312 220 +313 220 +314 220 +315 275 +401 20 +402 20 +403 20 +404 20 +405 275 +406 275 +407 275 +408 275 +409 275 +410 275 +411 275 +412 275 +413 275 +414 330 +415 330 +416 330 +501 20 +502 15 +503 15 +504 275 +505 275 +506 275 +507 275 +508 275 +509 275]; +% No from To v1(kv) v2(kv) rating x%(on rating) total number of units +%1 2 3 4 5 6 7 8 +trans=[... +1 101 102 15 330 333.3 12 12 +2 201 206 20 330 666.7 16 6 +3 202 209 20 330 555.6 16 5 +4 203 208 20 330 555.6 17 4 +5 204 215 20 330 666.7 16 6 +6 209 210 330 500 625 17 4 +7 213 214 500 330 625 17 4 +8 301 303 20 500 666.7 16 8 +9 302 312 20 220 444.4 15 4 +10 304 313 500 220 500 16 2 +11 305 311 500 330 500 12 2 +12 305 314 500 220 700 17 2 +13 308 315 500 275 370 10 2 +14 401 410 20 275 444.4 15 4 +15 402 408 20 275 333.3 17 3 +16 403 407 20 275 444.4 15 4 +17 404 405 20 275 333.3 17 6 +18 413 414 275 330 750 6 3 +19 501 504 20 275 333.3 17 2 +20 502 505 15 275 166.7 16 4 +21 503 506 15 275 250 16.7 6]; +active_trans_case4=[2,4,3,3,4,4,4,6,2,2,2,2,2,3,2,3,3,3,2,3,1]; + +trans(:,6)=trans(:,6).*active_trans_case4'; +trans(:,7)=trans(:,7)/100; + +%from to r X b(pu) v(kv) No +line=[... +102 217 0.0084 0.0667 0.817 330 1 +102 217 0.0078 0.062 0.76 330 2 +102 309 0.0045 0.0356 0.437 330 3 +102 309 0.0109 0.0868 0.76 330 4 +205 206 0.0096 0.076 0.931 330 5 +205 416 0.0037 0.046 0.73 330 6 +206 207 0.0045 0.0356 0.437 330 7 +206 212 0.0066 0.0527 0.646 330 8 +206 215 0.0066 0.0527 0.646 330 9 +207 208 0.0018 0.014 0.171 330 10 +207 209 0.0008 0.0062 0.076 330 11 +208 211 0.0031 0.0248 0.304 330 12 +209 212 0.0045 0.0356 0.437 330 13 +210 213 0.001 0.0145 1.54 500 14 +211 212 0.0014 0.0108 0.133 330 15 +211 214 0.0019 0.0155 0.19 330 16 +212 217 0.007 0.0558 0.684 330 17 +214 216 0.001 0.0077 0.095 330 18 +214 217 0.0049 0.0388 0.475 330 19 +215 216 0.0051 0.0403 0.494 330 20 +215 217 0.0072 0.0574 0.703 330 21 +216 217 0.0051 0.0403 0.494 330 22 +303 304 0.001 0.014 1.48 500 23 +303 305 0.0011 0.016 1.7 500 24 +304 305 0.0003 0.004 0.424 500 25 +305 306 0.0002 0.003 0.32 500 26 +305 307 0.0003 0.0045 0.447 500 27 +306 307 0.0001 0.0012 0.127 500 28 +307 308 0.0023 0.0325 3.445 500 29 +309 310 0.009 0.0713 0.874 330 30 +310 311 0 -0.0337 0 330 31 +312 313 0.002 0.015 0.9 220 32 +313 314 0.0005 0.005 0.52 220 33 +315 509 0.007 0.05 0.19 275 34 +405 406 0.0039 0.0475 0.381 275 35 +405 408 0.0054 0.05 0.189 275 36 +405 409 0.018 0.122 0.79 275 37 +406 407 0.0006 0.0076 0.062 275 38 +407 408 0.0042 0.0513 0.412 275 39 +408 410 0.011 0.128 1.01 275 40 +409 411 0.0103 0.0709 0.46 275 41 +410 411 0.0043 0.0532 0.427 275 42 +410 412 0.0043 0.0532 0.427 275 43 +410 413 0.004 0.0494 0.4 275 44 +411 412 0.0012 0.0152 0.122 275 45 +414 415 0.002 0.025 0.39 330 46 +415 416 0.0037 0.046 0.73 330 47 +504 507 0.023 0.15 0.56 275 48 +504 508 0.026 0.019 0.87 275 49 +505 507 0.0008 0.0085 0.06 275 50 +505 508 0.0025 0.028 0.17 275 51 +506 507 0.0008 0.0085 0.06 275 52 +506 508 0.003 0.028 0.14 275 53 +507 508 0.002 0.019 0.09 275 54 +507 509 0.03 0.22 0.9 275 55 +]; + zbase=(line(:,6).^2)/100; + line(:,3)=line(:,3).*zbase; + line(:,4)=line(:,4).*zbase/(2*pi*50); + line(:,5)=line(:,5)./zbase/(2*pi*50); + +% No bus QL Qc +comp=[... +1 211 0 0 +2 212 0 400 +3 216 0 300 +4 409 0 60 +5 411 0 30 +6 414 30 0 +7 415 60 0 +8 416 60 0 +9 504 0 0 +]; +SVC=[... +1 205 1.045 -39.3 +2 313 1.015 86.7 +3 412 1 -52.2 +4 507 1.01 -4 +5 509 1.03 -109.3]; + + PSSMODEL=1; + G=1; FL=.1; KL=10; FI=1; KI=10; FH=13; KH=120; + Tr=.001; + Ka=200; Ta=.001; + Ke=1; Te=0; + Tb=0; Tc=0; + Kf=0; Tf=0; + Efmin=-5; Efmax=5; + +alpha(1:15)=0*1e11*ones(15,1); +alpha(16:30)=0*1e8*ones(15,1); +nred=90; + +%PLL parameters +% Kp=.01; %original 180 +% Ki=.007; %original 3200 +PLL_Ki=2; +PLL_Kp=500; +visat=40*1e6;% original 100 MW +Isel=eye(2*15); +Isel_forming=eye(3*15); + +thetavsc0(1)=6.6166; +thetavsc0(2)=21.2703; +thetavsc0(3)=15.5222; +thetavsc0(4)=28.8703; +thetavsc0(5)=14.0249; +thetavsc0(6)=-6.3379; +thetavsc0(7)=4.7928; +thetavsc0(8)=3.6837; +thetavsc0(9)=-1.0986; +thetavsc0(10)=23.5812; +thetavsc0(11)=10.5944; +thetavsc0(12)=2.8249; +thetavsc0(13)=-12.7055; +thetavsc0(14)=-12.4459; +thetavsc0(15)=-14.7643; + +Kvsc(1,:)=1*ones(1,15); +Kvsc(2,:)=1*ones(1,15); +Kvsc=0*Kvsc; + +dist.inp=zeros(6,1); +dist.sc=ones(6,1); +dist.time=50; + + +r_fil=5e-1*ones(15,1); % 5e-1 +l_fil=5e-3*ones(15,1); % 5e-3 +c_fil=1e-5*ones(15,1); % 1e-5 + +% Grid forming +Pinit=[0.405553120017514,0.262363656270440,0.245287066765904,0.511573619715646,0.238910511748616,0.266289198874107,0.405451933295266,0.301363467877759,0.299904641553078,0.238601214125575,0.201576196829353,0.166515071179716,0.281316548159335,0.249279775700722,0.436119451429005]; diff --git a/linearizesystem.m b/linearizesystem.m new file mode 100644 index 0000000..6d69f92 --- /dev/null +++ b/linearizesystem.m @@ -0,0 +1,159 @@ +scomp=tmp.model; + +tmp.loadflow=power_loadflow('-v2',tmp.model,'NoUpdate'); + +io=getlinio(tmp.model); +set_param(tmp.model,'IgnoredZcDiagnostic','none'); +opt = linearizeOptions('UseFullBlockNameLabels','on');%'SampleTime',0 + +disp('linearizing') +tic + +switch scomp + case 'Australian14gen_following_inertia' + tlin=50; + case 'Australian14gen_forming_inertia' + tlin=30; +end + +tmp.lin=linearize(tmp.model,io,tlin); +tc=toc; +disp(['linearization completed in ',num2str(tc),' seconds']) + + + +switch scomp + + case 'Australian14gen_following_inertia' + disp('sanity checks on linearization') + + if norm(tmp.lin.D,'fro')>1e-6 + error('Feedthrough matrix is not zero'); + end + + etol=2e-3; + if any(real(eig(tmp.lin.A))>etol) + error(['the linearized system has an eigenvalue with real part larger than ',num2str(etol)]) + end + + + % Sanity Check + [tmp.Ve,tmp.De]=eig(tmp.lin.A); + [tmp.ev.magsorted,tmp.ev.idx]=sort(abs(real(diag(tmp.De)))); + + tmp.idx.zeroev=tmp.ev.idx(1); + + if norm(tmp.lin.C*tmp.Ve(:,tmp.idx.zeroev))>etol + error('The uncontrollable mode is not unobsevable') + end + + + disp('The linearized system passed the sanity checks') + + tmp.n_x=size(tmp.lin.A,2); + tmp.n_y=size(tmp.lin.C,1); + tmp.n_u=size(tmp.lin.B,2); + + tmp.mstring=['lin']; + + nstr = min(numel(tmp.model),numel(scomp)); + + + if tmp.model(1:nstr)==scomp(1:nstr) + tmp.n_vi=15; + tmp.n_dist=6; + tmp.n_xr=nred; + + tmp.idx_gen=1:10; + tmp.idx_wgen=1:10; + tmp.idx_ddtwgen=41:50; + tmp.idx_pvi=11:25; + tmp.idx_pgen=51:60; + tmp.idx_qvi=26:40; + + tmp.n_gen=length(tmp.idx_gen); + tmp.idx_pllw=61:2:90; + tmp.idx_plldw=62:2:90; + + %sparsity pattern + Ftmp=zeros(30,90); + Ftmp(1:tmp.n_vi,tmp.idx_pllw)=eye(tmp.n_vi); + Ftmp(tmp.n_vi+1:2*tmp.n_vi,tmp.idx_plldw)=2*eye(tmp.n_vi); + + tmp.SP.d=(Ftmp==1); + tmp.SP.m=(Ftmp==2); + clear Ftmp + + tmp.opt.x0=0*ones(30,1); + end + + case 'Australian14gen_forming_inertia' + + + tmp.n_vi=15; + tmp.n_dist=6; + tmp.n_xr=nred; + + tmp.idx_gen=1:10; + tmp.idx_wgen=1:10; + tmp.idx_ddtwgen=41:50; + tmp.idx_pgen=51:60; + + tmp.n_gen=length(tmp.idx_gen); + tmp.idx_p=11:25; + tmp.idx_pvi=11:25; + tmp.idx_qvi=26:40; + + tmp.idx_w=61:75; + tmp.idx_wvi=61:75; + + disp('sanity checks on linearization') + + if norm(tmp.lin.D,'fro')>1e-6 + error('Feedthrough matrix is not zero'); + end + + %temporary control gain for sanity checks + tmp.F=kron([1 1],eye(15)); + + %temporary closed loop for sanity checks + tmp.Acl=tmp.lin.A-tmp.lin.B(:,tmp.n_dist+1:end)*tmp.F*tmp.lin.C([tmp.idx_p,tmp.idx_w],:); + + + etol=1e-5; + if any(real(eig(tmp.Acl))>etol) + error(['the linearized system has an eigenvalue with real part larger than ',num2str(etol)]) + end + + + % Sanity Check + [tmp.Ve,tmp.De]=eig(tmp.Acl); + [tmp.ev.magsorted,tmp.ev.idx]=sort(abs(real(diag(tmp.De)))); + + tmp.idx.zeroev=tmp.ev.idx(1); + + if norm(tmp.lin.C*tmp.Ve(:,tmp.idx.zeroev))>etol + error('The uncontrollable mode is not unobsevable') + end + + + disp('The linearized system passed the sanity checks') + + tmp.n_x=size(tmp.lin.A,2); + tmp.n_y=size(tmp.lin.C,1); + tmp.n_u=size(tmp.lin.B,2); + + tmp.mstring=['lin']; + + %sparsity pattern + Ftmp=zeros(tmp.n_vi,tmp.n_y); + Ftmp(:,tmp.idx_pvi)=eye(tmp.n_vi); + Ftmp(:,tmp.idx_wvi)=2*eye(tmp.n_vi); + + tmp.SP.p=(Ftmp==1); + tmp.SP.w=(Ftmp==2); + clear Ftmp + + tmp.opt.x0=-1*ones(30,1); + tmp.opt.x0(16:30)=-1*ones(15,1); +end diff --git a/main_following.m b/main_following.m new file mode 100644 index 0000000..260f9f2 --- /dev/null +++ b/main_following.m @@ -0,0 +1,89 @@ + +close all +clear all +bdclose all + +%Nonlinear Model +sys.follin.model='Australian14gen_following_inertia'; + +% Disturbance inputs +sys.follin.dist.inp(1:6,1)=0; +sys.follin.dist.inp(4)=-250e6; +sys.follin.dist.time=20; + +%load data +dataaustralian14gen +open(sys.follin.model) + +%set up filter for the RoCoF +s=tf('s'); +Gf=eye(10)*(s/(.2/3*s+1)); +sys.follin.filt=ss(Gf); +clear s Gf + +%Linearize the nonlinear system +tmp=sys.follin; +tmp.GenPU=diag(G_data(gentestcase,2)); +linearizesystem; +sys.follin=tmp; + +sys.follin.opt.m0=sum(2*G_data(genrem,4).*G_data(genrem,2)./(2*pi*50))/sys.follin.n_vi; %in MW / (rad/s^2) +sys.follin.opt.d0=sum(G_data(genrem,2)*1/(2*pi*50*0.05))/sys.follin.n_vi; %in MW / (rad/s) +sys.follin.opt.x0=-[ones(sys.follin.n_vi,1)*sys.follin.opt.d0;ones(sys.follin.n_vi,1)*sys.follin.opt.m0]./(1e2/(2*pi*50)); %to 100MW, 50 Hz base + +%Remove the zero eigenvalue +sys.follin=modelreduction(sys.follin); + + +% Gains for non-linear simulation, input in the optimization is in +% 100 MW base, input in the nonlinear model is in W. +sys.follin.opt.alpha=100 *[1.5752;3.0736;2.4685;2.5607;2.3273;1.3412;1.4971;2.4483;1.8394;0.7652;1.8701; 7.2121; 3.3551; 3.8545; 2.3868; 0.0599; 0.1003; 0.0148; 0.2209; 0.0755; 0.0075; 0.0394; 0.3280; 0.0210; 0.0312; 0.1314; 0.0914; 1.7089; 1.6426; 1.6138]; + +% Gains for linear simulation +sys.follin.opt.d=-sys.follin.opt.alpha(1:sys.follin.n_vi,1); +sys.follin.opt.m=-sys.follin.opt.alpha(sys.follin.n_vi+1:2*sys.follin.n_vi,1); + +% Form the feedback matrix for the linearized system +F=zeros(30,90); +F(sys.follin.SP.d==1)=sys.follin.opt.d; +F(sys.follin.SP.m==1)=sys.follin.opt.m; +sys.follin.opt.F=F; + +sys.follin.opt.B=[sys.follin.red.B(:,sys.follin.n_dist+1:end),sys.follin.red.B(:,sys.follin.n_dist+1:end)]; +sys.follin.opt.G=sys.follin.red.B(:,1:sys.follin.n_dist); + +sys.follin.tend=30; +sys.follin=tdsim(sys.follin); + +% Non linear simulation without VI +dist=sys.follin.dist; +disp('simulating nonlinear system with disturbance and without VI') +tic +SimOut = sim(sys.follin.model,'StopTime',int2str(sys.follin.dist.time+sys.follin.tend)); +tc=toc; + +idx.disttime=find(SimOut.ysim.Time>sys.follin.dist.time,1,'first'); +sys.follin.tdsim.nl.Ts=sys.follin.tdsim.red.Ts; + +sys.follin.tdsim.nl.Ys=interp1(SimOut.ysim.Time(idx.disttime:end)-sys.follin.dist.time,SimOut.ysim.Data(idx.disttime:end,:)-SimOut.ysim.Data(idx.disttime,:),sys.follin.tdsim.nl.Ts); + +disp(['nonlinear simulation completed in ',num2str(tc),' seconds']) + +% Non linear simulation with VI +disp('simulating nonlinear system with disturbance and VI') +alpha=-sys.follin.opt.alpha*1e8; +sys.follin.dist.time=10; +dist.time=sys.follin.dist.time; +tic +SimOut = sim(sys.follin.model,'StopTime',int2str(sys.follin.dist.time+sys.follin.tend)); +tc=toc; + +idx.disttime=find(SimOut.ysim.Time>sys.follin.dist.time,1,'first'); +sys.follin.tdsim.nlcl.Ts=sys.follin.tdsim.red.Ts; +sys.follin.tdsim.nlcl.Ys=interp1(SimOut.ysim.Time(idx.disttime:end)-sys.follin.dist.time,SimOut.ysim.Data(idx.disttime:end,:)-SimOut.ysim.Data(idx.disttime,:),sys.follin.tdsim.nlcl.Ts); + +disp(['nonlinear simulation completed in ',num2str(tc),' seconds']) + + +%plot the results +plotvi(sys.follin); \ No newline at end of file diff --git a/main_forming.m b/main_forming.m new file mode 100644 index 0000000..727a8fd --- /dev/null +++ b/main_forming.m @@ -0,0 +1,78 @@ + +close all +clear all +bdclose all + +% Nonlinear system +sys.formin.model='Australian14gen_forming_inertia'; +open(sys.formin.model) + + +%load data +dataaustralian14gen +open(sys.formin.model) + +%set up filter for the RoCoF +s=tf('s'); +Gf=eye(10)*(s/(.2/3*s+1)); +sys.formin.filt=ss(Gf); +clear s Gf + +%Linearize the nonlinear system +tmp=sys.formin; +tmp.GenPU=diag(G_data(gentestcase,2)); +linearizesystem; +sys.formin=tmp; + +% Disturbance inputs +sys.formin.dist.inp(1:6,1)=0; +sys.formin.dist.inp(4)=-250e6; +sys.formin.dist.time=15; + +%Remove the zero eigenvalue +sys.formin=modelreduction(sys.formin); + + +% create extended Input matrices +sys.formin.opt.B=sys.formin.red.B(:,sys.formin.n_dist+1:end); +sys.formin.opt.C=sys.formin.red.C; +sys.formin.opt.G=sys.formin.red.B(:,1:sys.formin.n_dist); +sys.formin.opt.Kv=sys.formin.opt.x0; + +% Set the gains for non-linear simulation +sys.formin.Kv=zeros(30,1); +sys.formin.Kv(1:2:30)=sys.formin.opt.Kv(1:sys.formin.n_vi); +sys.formin.Kv(2:2:30)=sys.formin.opt.Kv(sys.formin.n_vi+1:2*sys.formin.n_vi)*1e3/(2*pi*50); + +% Form the Feedback matrix for the linear simulation +F=zeros(sys.formin.n_vi,sys.formin.n_y); +F(sys.formin.SP.p==1)=sys.formin.opt.Kv(1:sys.formin.n_vi,1); +F(sys.formin.SP.w==1)=sys.formin.opt.Kv(sys.formin.n_vi+1:2*sys.formin.n_vi,1); + +sys.formin.opt.F=F; + +%linear simulation +sys.formin.tend=30; +sys.formin=tdsim(sys.formin); + + +% Non linear simulation using the gains obtained via the H2 optimization +dist=sys.formin.dist; +disp('simulating nonlinear system with disturbance and VI') +Kvsc(2,1:1:15)=sys.formin.Kv(1:2:30)'; +Kvsc(1,1:1:15)=sys.formin.Kv(2:2:30)'; +sys.formin.dist.time=15; +dist.time=sys.formin.dist.time; + +tic +SimOut = sim(sys.formin.model,'StopTime',int2str(sys.formin.dist.time+sys.formin.tend)); +tc=toc; + +idx.disttime=find(SimOut.ysim.Time>sys.formin.dist.time,1,'first'); +sys.formin.tdsim.nlcl.Ts=sys.formin.tdsim.red.Ts; +sys.formin.tdsim.nlcl.Ys=interp1(SimOut.ysim.Time(idx.disttime:end)-sys.formin.dist.time,SimOut.ysim.Data(idx.disttime:end,:)-SimOut.ysim.Data(idx.disttime,:),sys.formin.tdsim.nlcl.Ts); + +disp(['nonlinear simulation completed in ',num2str(tc),' seconds']) + +%plot the results +plotvi(sys.formin); \ No newline at end of file diff --git a/modelreduction.m b/modelreduction.m new file mode 100644 index 0000000..05fecf9 --- /dev/null +++ b/modelreduction.m @@ -0,0 +1,64 @@ +function sys=modelreduction(sys) + + +disp('zero eigenvalue mode corresponds to the following states:') +sys.lin.StateName(abs(sys.Ve(:,sys.idx.zeroev))>1e-3) +sys.vistrings=sys.lin.StateName(abs(sys.Ve(:,sys.idx.zeroev))>1e-3); + + switch sys.model + + case 'Australian14gen_following_inertia' + +%extract buses of vi devices +sys.vistrings=cell2mat(sys.vistrings(end-sys.n_vi+1:end)); + +%extract bus name +sys.vistrings = sys.vistrings(:,1:5); + +%Transformed system n=372 +[U,S,V]=svd(sys.lin.A); + +%shift the zero EV +sys.shift.A=U(:,1:end-1)*S(1:end-1,1:end-1)*V(:,1:end-1)'; +sys.shift.B=sys.lin.B; +sys.shift.C=sys.lin.C; + +%remove it +sys.rem.A=V(:,1:end-1)'*sys.shift.A*V(:,1:end-1); +sys.rem.B=V(:,1:end-1)'*sys.shift.B; +sys.rem.C=sys.shift.C*V(:,1:end-1); + +% Balanced trunction may be used to reduce the system size and speed up the +% computation +% +% sys.n_xr=382; +% sys.red=balred(ss(sys.rem.A,sys.rem.B*1e8,sys.rem.C,[]),sys.n_xr,'StateElimMethod','Truncate'); + +sys.red=ss(sys.rem.A,sys.rem.B*1e8,sys.rem.C,[]); +sys.mstring=[sys.mstring,'rem']; + + case 'Australian14gen_forming_inertia' + + +%extract buses of vi devices +sys.vistrings=cell2mat(sys.vistrings(end-sys.n_vi+1:end)); + +%extract bus name +sys.vistrings = sys.vistrings(:,1:6); + + +%not much to do here. Rescaling the system +scal=eye(75,75); +scal(sys.idx_wgen, sys.idx_wgen)=1e3*eye(10); +scal(sys.idx_ddtwgen, sys.idx_ddtwgen)=1e3*eye(10); +scal(sys.idx_wvi, sys.idx_wvi)=1e3*eye(15); + +sys.red.A=sys.lin.A; +sys.red.B=sys.lin.B*diag([1e8*ones(6,1);ones(15,1)]); +sys.red.C=scal*sys.lin.C; + + +sys.mstring=[sys.mstring,'rem']; + + + end \ No newline at end of file diff --git a/plotvi.m b/plotvi.m new file mode 100644 index 0000000..51b80ce --- /dev/null +++ b/plotvi.m @@ -0,0 +1,280 @@ +function plotvi(sys) + +switch sys.model + case 'Australian14gen_following_inertia' + for i=1:15 + vilabel{i}=sys.vistrings(i,3:end); + end + + cy=[255, 170, 51]/255; + cb=[39, 58, 95]/255; + + figure(10011) + b=bar([-sys.opt.d,-sys.opt.m]*1e2/(2*pi*50),'grouped'); % The units were in MW/pu + b(2).FaceColor=cy; + b(1).FaceColor=cb; + legend('damping[MW s/rad]','inertia[MW s^2/rad]','orientation','horizontal','location','northwest') + set(gca,'xticklabel',vilabel) + grid on + box on + + figure(10021),clf + subplot(4,5,1) + plot(sys.tdsim.red.Ts,50*sys.tdsim.red.Ys(:,sys.idx_wgen)) + grid on + box on + xlabel('t [s]') + ylabel('\omega [Hz]') + title('linear, reduced, open loop') + + subplot(4,5,2) + plot(sys.tdsim.red.Ts,50*sys.tdsim.red.Ys(:,sys.idx_ddtwgen)) + grid on + box on + xlabel('t [s]') + ylabel('\dot \omega [Hz/s]') + title('linear, reduced, open loop') + + subplot(4,5,3) + plot(sys.tdsim.red.Ts,sys.GenPU*sys.tdsim.red.Ys(:,sys.idx_pgen)') + grid on + box on + xlabel('t [s]') + ylabel('P_{gen} [MW]') + title('linear, reduced, open loop') + + subplot(4,5,4) + plot(sys.tdsim.red.Ts,100*sys.tdsim.red.Ys(:,sys.idx_pvi)) + grid on + box on + xlabel('t [s]') + ylabel('P_{vi} [MW]') + title('linear, reduced, open loop') + + subplot(4,5,5) + plot(sys.tdsim.red.Ts,100*sys.tdsim.red.Ys(:,sys.idx_qvi)) + grid on + box on + xlabel('t [s]') + ylabel('Q_{vi} [MW]') + title('linear, reduced, open loop') + + + subplot(4,5,6) + plot(sys.tdsim.nl.Ts,50*sys.tdsim.nl.Ys(:,sys.idx_wgen)) + grid on + box on + xlabel('t [s]') + ylabel('\omega [Hz]') + title('nonlinear, open loop') + + subplot(4,5,7) + plot(sys.tdsim.nl.Ts,50*sys.tdsim.nl.Ys(:,sys.idx_ddtwgen)) + grid on + box on + xlabel('t [s]') + ylabel('\dot \omega [Hz/s]') + title('nonlinear, open loop') + + subplot(4,5,8) + plot(sys.tdsim.nl.Ts,sys.GenPU*sys.tdsim.nl.Ys(:,sys.idx_pgen)') + grid on + box on + xlabel('t [s]') + ylabel('P_{gen} [MW]') + title('nonlinear, open loop') + + subplot(4,5,9) + plot(sys.tdsim.nl.Ts,100*sys.tdsim.nl.Ys(:,sys.idx_pvi)) + grid on + box on + xlabel('t [s]') + ylabel('P_{vi} [MW]') + title('nonlinear, open loop') + + subplot(4,5,10) + plot(sys.tdsim.nl.Ts,100*sys.tdsim.nl.Ys(:,sys.idx_qvi)) + grid on + box on + xlabel('t [s]') + ylabel('Q_{vi} [MW]') + title('nonlinear, open loop') + + subplot(4,5,11) + plot(sys.tdsim.redcl.Ts,50*sys.tdsim.redcl.Ys(:,sys.idx_wgen)) + grid on + box on + xlabel('t [s]') + ylabel('\omega [Hz]') + title('linear, reduced, closed loop') + + subplot(4,5,12) + plot(sys.tdsim.redcl.Ts,50*sys.tdsim.redcl.Ys(:,sys.idx_ddtwgen)) + grid on + box on + xlabel('t [s]') + ylabel('\dot \omega [Hz/s]') + title('linear, reduced, closed loop') + + subplot(4,5,13) + plot(sys.tdsim.redcl.Ts,sys.GenPU*sys.tdsim.redcl.Ys(:,sys.idx_pgen)') + grid on + box on + xlabel('t [s]') + ylabel('P_{gen} [MW]') + title('linear, reduced, closed loop') + + subplot(4,5,14) + plot(sys.tdsim.redcl.Ts,100*sys.tdsim.redcl.Ys(:,sys.idx_pvi)) + grid on + box on + xlabel('t [s]') + ylabel('P_{vi} [MW]') + title('linear, reduced, closed loop') + + subplot(4,5,15) + plot(sys.tdsim.redcl.Ts,100*sys.tdsim.redcl.Ys(:,sys.idx_qvi)) + grid on + box on + xlabel('t [s]') + ylabel('Q_{vi} [MW]') + title('linear, reduced, closed loop') + + subplot(4,5,16) + plot(sys.tdsim.nlcl.Ts,50*sys.tdsim.nlcl.Ys(:,sys.idx_wgen)) + grid on + box on + xlabel('t [s]') + ylabel('\omega [Hz]') + title('nonlinear, closed loop') + + subplot(4,5,17) + plot(sys.tdsim.nlcl.Ts,50*sys.tdsim.nlcl.Ys(:,sys.idx_ddtwgen)) + grid on + box on + xlabel('t [s]') + ylabel('\dot \omega [Hz/s]') + title('nonlinear, closed loop') + + subplot(4,5,18) + plot(sys.tdsim.nlcl.Ts,sys.GenPU*sys.tdsim.nlcl.Ys(:,sys.idx_pgen)') + grid on + box on + xlabel('t [s]') + ylabel('P_{gen} [MW]') + title('nonlinear, closed loop') + + subplot(4,5,19) + plot(sys.tdsim.nlcl.Ts,100*sys.tdsim.nlcl.Ys(:,sys.idx_pvi)) + grid on + box on + xlabel('t [s]') + ylabel('P_{vi} [MW]') + title('nonlinear, closed loop') + + subplot(4,5,20) + plot(sys.tdsim.nlcl.Ts,100*sys.tdsim.nlcl.Ys(:,sys.idx_qvi)) + grid on + box on + xlabel('t [s]') + ylabel('Q_{vi} [MW]') + title('nonlinear, closed loop') + + + + + + case 'Australian14gen_forming_inertia' + + for i=1:15 + vilabel{i}=sys.vistrings(i,3:end); + end + + cy=[255, 170, 51]/255; + cb=[39, 58, 95]/255; + + m_opt=-sys.Kv(1:2:30).^-1*100; + d_opt=sys.Kv(2:2:30)./sys.Kv(1:2:30)*100; + m_opt(m_opt=='inf')=0; + d_opt(d_opt=='inf')=0; + + figure(1001) + b=bar([d_opt,m_opt],'grouped'); + b(2).FaceColor=cy; + b(1).FaceColor=cb; + legend('damping [MW s/rad]','inertia [MW s^2/rad]','orientation','horizontal','location','northwest') + set(gca,'xticklabel',vilabel) + grid on + box on + + + + figure(1002),clf + + + subplot(2,4,1) + plot(sys.tdsim.redcl.Ts,50*sys.tdsim.redcl.Ys(:,sys.idx_wgen)/1000) + grid on + box on + xlabel('t [s]') + ylabel('\omega [Hz]') + title('linear, reduced, closed loop') + + subplot(2,4,2) + plot(sys.tdsim.redcl.Ts,50*sys.tdsim.redcl.Ys(:,sys.idx_ddtwgen)/1000) + grid on + box on + xlabel('t [s]') + ylabel('\dot \omega [Hz/s]') + title('linear, reduced, closed loop') + + subplot(2,4,3) + plot(sys.tdsim.redcl.Ts,sys.GenPU*sys.tdsim.redcl.Ys(:,sys.idx_pgen)') + grid on + box on + xlabel('t [s]') + ylabel('P_{gen} [MW]') + title('linear, reduced, closed loop') + + subplot(2,4,4) + plot(sys.tdsim.redcl.Ts,100*sys.tdsim.redcl.Ys(:,sys.idx_pvi)) + grid on + box on + xlabel('t [s]') + ylabel('P_{vi} [MW]') + title('linear, reduced, closed loop') + + + subplot(2,4,5) + plot(sys.tdsim.nlcl.Ts,50*sys.tdsim.nlcl.Ys(:,sys.idx_wgen)) + grid on + box on + xlabel('t [s]') + ylabel('\omega [Hz]') + title('nonlinear, closed loop') + + subplot(2,4,6) + plot(sys.tdsim.nlcl.Ts,50*sys.tdsim.nlcl.Ys(:,sys.idx_ddtwgen)) + grid on + box on + xlabel('t [s]') + ylabel('\dot \omega [Hz/s]') + title('nonlinear, closed loop') + + subplot(2,4,7) + plot(sys.tdsim.nlcl.Ts,sys.GenPU*sys.tdsim.nlcl.Ys(:,sys.idx_pgen)') + grid on + box on + xlabel('t [s]') + ylabel('P_{gen} [MW]') + title('nonlinear, closed loop') + + subplot(2,4,8) + plot(sys.tdsim.nlcl.Ts,100*sys.tdsim.nlcl.Ys(:,sys.idx_pvi)) + grid on + box on + xlabel('t [s]') + ylabel('P_{vi} [MW]') + title('nonlinear, closed loop') + +end diff --git a/tdsim.m b/tdsim.m new file mode 100644 index 0000000..0a04928 --- /dev/null +++ b/tdsim.m @@ -0,0 +1,56 @@ +function sys=tdsim(sys) + +switch sys.model + + case 'Australian14gen_following_inertia' + + %simulate linearized system with reduced order + disp('simulating linearized system with disturbance and no VI') + tic + + [sys.tdsim.red.Ys,sys.tdsim.red.Ts]=step(ss(sys.red.A,sys.opt.G(:,1:sys.n_dist)*sys.dist.inp/1e8,sys.red.C,[]),sys.tend); + + tc=toc; + disp(['linear simulation completed in ',num2str(tc),' seconds']) + + %compute damping ratios + sys.tdsim.red.damp_ratio=min(-real(eig(sys.red.A))./(abs(eig(sys.red.A)))); + + %simulate linearized system with reduced order and VI + disp('simulating linearized system with disturbance and VI') + tic + + [sys.tdsim.redcl.Ys,sys.tdsim.redcl.Ts]=step(ss(sys.red.A+sys.opt.B*sys.opt.F*sys.red.C,sys.opt.G(:,1:sys.n_dist)*sys.dist.inp/1e8,sys.red.C,[]),sys.tdsim.red.Ts); + tc=toc; + disp(['linear simulation completed in ',num2str(tc),' seconds']) + + %compute damping ratios + sys.tdsim.redcl.damp_ratio=min(-real(eig(sys.red.A+sys.opt.B*sys.opt.F*sys.red.C))./(abs(eig(sys.red.A+sys.opt.B*sys.opt.F*sys.red.C)))); + + + + case 'Australian14gen_forming_inertia' + + %simulate linearized system with disturbance and pre-stabilising VI + disp('simulating linearized system with disturbance and pre-stabilizing VI') + tic + [sys.tdsim.red.Ys,sys.tdsim.red.Ts]=step(ss(sys.Acl,sys.opt.G(:,1:sys.n_dist)*sys.dist.inp/1e8,sys.red.C,[]),sys.tend); + tc=toc; + disp(['linear simulation completed in ',num2str(tc),' seconds']) + + %compute damping ratios + sys.tdsim.red.damp_ratio=min(-real(eig(sys.Acl))./(abs(eig(sys.Acl)))); + + + %simulate linearized system with disturbance and VI + disp('simulating linearized system with disturbance and VI') + tic + [sys.tdsim.redcl.Ys,sys.tdsim.redcl.Ts]=step(ss(sys.red.A+sys.opt.B*sys.opt.F*sys.opt.C,sys.opt.G(:,1:sys.n_dist)*sys.dist.inp/1e8,sys.red.C,[]),sys.tdsim.red.Ts); + tc=toc; + disp(['linear simulation completed in ',num2str(tc),' seconds']) + + %compute damping ratios + sys.tdsim.redcl.damp_ratio=min(-real(eig(sys.red.A+sys.opt.B*sys.opt.F*sys.opt.C))./(abs(eig(sys.red.A+sys.opt.B*sys.opt.F*sys.opt.C)))); + + +end