s_simu.m

这是在MATLAB环境下进行电力系统暂态稳定分析的pst主程序。使用它可对特定故障的电力系统模型进行暂态稳定分析

% s_simu.m  
% 9:39 AM 16 December 1997
% An m.file to simulate power system transients  
% using the Matlab Power System Toolbox
% This m-file takes the dynamic and load flow data and
% calculates the response of the power system to a fault
% which is specified in a switching file
% see one of the supplied data files (data2a.m) for the 
% switching file format

% Version: 1.2
% Author:  Graham Rogers
% Date:    August 1997
% Modification : add induction generator
% Version: 1.1
% Author:  Graham Rogers
% Date:    August 1997
% Modification : add modulation for load, exciter and governor
% Version  1.0
% Author:  Graham Rogers
% Date: February 1997
% (c) Copyright: Joe Chow/ Cherry Tree Scientific Software 1991 to 1997 - All rights reserved
%
clear
clear global
close % close graphics windows
tic % set timer
plot_now=0;
jay = sqrt(-1);
pst_var % set up global variables 

% load input data from m.file
disp('non-linear simulation')
% input data file
[dfile,pathname]=uigetfile('d*.m','Select  File');
if pathname == 0
   error(' you must select a valid data file')
else
   lfile =length(dfile);
   % strip off .m and convert to lower case
   dfile = lower(dfile(1:lfile-2));
   eval(dfile);
end
% check for valid dynamic data file
if isempty(mac_con)
   error(' the selected file is not a valid data file')
end
if isempty(sw_con)
   error(' the selected file has no switching data')
end
sys_freq = input('enter the base system frequency in Hz - [60]');
if isempty(sys_freq);sys_freq = 60;end
basrad = 2*pi*sys_freq; % default system frequency is 60 Hz
basmva = input('enter system base MVA - [100]');
if isempty(basmva);basmva =100;end
syn_ref = 0 ;     % synchronous reference frame
ibus_con = []; % ignore infinite buses in transient simulation

% solve for loadflow - loadflow parameter
ans = input('Do you want to solve loadflow > (y/n)[y] ','s');
if isempty(ans);ans='y';end
if ans=='Y';ans='y';end
if ans == 'y'
   if isempty(dcsp_con)
      n_conv = 0;
      n_dcl = 0;
      tol = 1e-9;   % tolerance for convergence
      iter_max = 10; % maximum number of iterations
      acc = 1.0;   % acceleration factor
      [bus_sol,line,line_flw] = ...
         loadflow(bus,line,tol,iter_max,acc,'n',2);
      bus = bus_sol;  % solved loadflow solution needed for
      % initialization
      save sim_fle bus line n_conv n_dcl
   else
      [bus_sol,line,line_flw,rec_par,inv_par, line_par] = lfdcs(bus,line);
      bus = bus_sol;
      save sim_fle bus line rec_par  inv_par line_par
   end 
else
   load sim_fle 
end
%set machine, exciter, governor and pss indexes
% note: dc index set in dc load flow
f=mac_indx;
f=exc_indx;
f=tg_indx;
f=pss_indx;
f = svc_indx;
f = lm_indx;
f = rlm_indx;
[n_mot,dummy] = size(ind_con);
[n_ig,dummy] = size(igen_con);
if isempty(n_mot); n_mot = 0;end
if isempty(n_ig); n_ig = 0; end
ntot = n_mac+n_mot+n_ig;
ngm = n_mac + n_mot;
n_pm = n_mac;
disp(' ')
disp('Performing simulation.')
%
% construct simulation switching sequence as defined in sw_con
tswitch(1) = sw_con(1,1);
k = 1;
n_switch = length(sw_con(:,1));
k_inc = zeros(n_switch-1,1);
t_switch = zeros(n_switch,1);
for sw_count = 1:n_switch-1
   h(sw_count) = sw_con(sw_count,7);
   if h(sw_count)==0, h(sw_count) = 0.01;end
   k_inc(sw_count) = fix((sw_con(sw_count+1,1)-sw_con(sw_count,1))/h(sw_count));
   if k_inc(sw_count)==0;k_inc(sw_count)=1;end
   h(sw_count) = (sw_con(sw_count+1,1)-sw_con(sw_count,1))/k_inc(sw_count);
   t_switch(sw_count+1) =t_switch(sw_count) +  k_inc(sw_count)*h(sw_count);
   t(k:k-1+k_inc(sw_count)) = t_switch(sw_count):h(sw_count):t_switch(sw_count+1)-h(sw_count);
   k=k+k_inc(sw_count);
end
%
k = sum(k_inc)+1; % k is the total number of time steps in the simulation
t(k) = sw_con(n_switch,1);
[n dummy]=size(mac_con) ;
n_bus = length(bus(:,1));
%
% create zero matrices for variables to make algorithm more efficient?
z = zeros(n,k);
z1 = zeros(1,k);
zm = zeros(1,k);if n_mot>1;zm = zeros(n_mot,k);end
zig = zeros(1,k);if n_ig>1;zig = zeros(n_ig,k);end
zdc = zeros(2,k);if n_conv>2; zdc = zeros(n_conv,k);end
zdcl = zeros(1,k);if n_dcl>1;zdcl=zeros(n_dcl,k);end
% set dc parameters  
Vdc = zdc;
i_dc = zdc;  
P_dc = zdc; dc_sig = zdc; cur_ord = zdc;
alpha = zdcl; 
gamma = zdcl;  
dc_sig = zdc;
i_dcr = zdcl; i_dci = zdcl; v_dcc = zdcl;
di_dcr = zdcl; di_dci = zdcl; dv_dcc = zdcl;
v_conr = zdcl; v_coni = zdcl; dv_conr = zdcl; dv_coni = zdcl;
if n_conv~=0
   Vdc(r_idx,1) = rec_par(:,2); Vdc(i_idx,1) = inv_par(:,2);
   i_dc(r_idx,1) = line_par; i_dc(i_idx,1) = line_par;
   i_dcr(:,1) = i_dc(r_idx,1); i_dci(:,1) = i_dc(i_idx,1);
   alpha(:,1) = rec_par(:,1)*pi/180;
   gamma(:,1) = inv_par(:,1)*pi/180;
end
v_p = z1;
theta = zeros(n_bus+1,k);bus_v = zeros(n_bus+1,k);
mac_ang = z; mac_spd = z; dmac_ang = z; dmac_spd = z;
pmech = z; pelect = z; mac_ref = z1;  sys_ref = z1; 
edprime = z; eqprime = z; dedprime = z; deqprime = z;
psikd = z; psikq = z; dpsikd = z; dpsikq = z;
pm_sig = z;
z_tg = zeros(1,k);if n_tg~=0;z_tg = zeros(n_tg,k);end
tg1 = z_tg; tg2 = z_tg; tg3 = z_tg; 
dtg1 = z_tg; dtg2 = z_tg; dtg3 = z_tg;
z_pss = zeros(1,k);if n_pss~=0;z_pss = zeros(n_pss,k);end
pss1 = z_pss; pss2 = z_pss; pss3 = z_pss;
dpss1 = z_pss; dpss2 = z_pss; dpss3 = z_pss;
curd = z; curq = z; curdg = z; curqg = z; fldcur = z;
ed = z; eq = z; eterm = z; qelect = z;
vex = z; cur_re = z; cur_im = z; psi_re = z; psi_im = z;
ze = zeros(1,k);if n_exc~=0; ze = zeros(n_exc,k);end
V_B = ze;exc_sig = ze;
V_TR = ze; V_R = ze; V_A = ze; V_As = ze; Efd = ze; R_f = ze;
dV_TR = ze; dV_R = ze; dV_As = ze; dEfd = ze; dR_f = ze;
pss_out = ze;
vdp = zm; vqp = zm; slip = zm; 
dvdp = zm; dvqp = zm; dslip = zm;
s_mot = zm; p_mot = zm; q_mot = zm;
vdpig = zig; vqpig = zig; slig = zig; 
dvdpig = zig; dvqpig = zig; dslig = zig;
s_igen = zig; pig = zig; qig = zig; tmig = zig;
if n_svc~=0
   B_cv = zeros(n_svc,k); dB_cv = zeros(n_svc,k);svc_sig = zeros(n_svc,k);
   B_con = zeros(n_svc,k); dB_con=zeros(n_svc,k);
else
   B_cv = zeros(1,k);dB_cv = zeros(1,k); svc_sig = zeros(1,k);
   B_con = zeros(1,k);dB_con=zeros(1,k);
end
if n_lmod ~= 0
   lmod_st = zeros(n_lmod,k); dlmod_st = lmod_st; lmod_sig = lmod_st;
else
   lmod_st = zeros(1,k); dlmod_st = lmod_st; lmod_sig = lmod_st;
end
if n_rlmod ~= 0
   rlmod_st = zeros(n_rlmod,k); drlmod_st = rlmod_st; rlmod_sig = rlmod_st;
else
   rlmod_st = zeros(1,k); drlmod_st = rlmod_st; rlmod_sig = rlmod_st;
end

sys_freq = ones(1,k);
disp('constructing reduced y matrices')
% step 1: construct reduced Y matrices 

disp('initializing motor,induction generator, svc and dc control models')       
bus = mac_ind(0,1,bus,0);% initialize induction motor
bus = mac_igen(0,1,bus,0); % initialize induction generator
bus = svc(0,1,bus,0);%initialize svc
f = dc_cont(0,1,bus,0);% initialize dc controls
% this has to be done before red_ybus is used since the motor and svc 
% initialization alters the bus matrix and dc parameters are required

y_switch % calculates the reduced y matrices for the different switching conditions
disp('initializing other models')
% step 2: initialization
theta(1:n_bus,1) = bus(:,3)*pi/180;
bus_v(1:n_bus,1) = bus(:,2).*exp(jay*theta(1:n_bus,1));
if n_conv~=0
   % change dc buses from LT to HT
   Pr = bus(rec_ac_bus,6);
   Pi = bus(inv_ac_bus,6);
   Qr = bus(rec_ac_bus,7);
   Qi = bus(inv_ac_bus,7);
   VLT= bus_v(ac_bus,1);
   i_acr = (Pr-jay*Qr)./conj(VLT(r_idx));
   i_aci = (Pi - jay*Qi)./conj(VLT(i_idx));
   IHT(r_idx,1)=i_acr;
   IHT(i_idx,1)=i_aci;
   VHT(r_idx,1) = (VLT(r_idx) + jay*dcc_pot(:,2).*i_acr);
   VHT(i_idx,1) = (VLT(i_idx) + jay*dcc_pot(:,4).*i_aci);
   bus_v(ac_bus,1) = VHT;
   theta(ac_bus,1) = angle(bus_v(ac_bus,1));
   % modify the bus matrix to the HT buses
   bus(ac_bus,2) = abs(bus_v(ac_bus,1));
   bus(ac_bus,3) = theta(ac_bus,1)*180/pi;
   SHT = VHT.*conj(IHT);
   bus(ac_bus,6) = real(SHT);
   bus(ac_bus,7) = imag(SHT);
end
flag = 0;
bus_int = bus_intprf;% pre-fault system
disp('generators')
f = mac_sub(0,1,bus,flag);   
f = mac_tra(0,1,bus,flag);
f = mac_em(0,1,bus,flag);
disp('generator controls')
f = pss(0,1,bus,flag);
f = smpexc(0,1,bus,flag);
f = exc_st3(0,1,bus,flag);
f = exc_dc12(0,1,bus,flag);
f = tg(0,1,bus,flag); 
% initialize load modulation control
if ~isempty(lmod_con)
   disp('load modulation')
   f = lmod(0,1,bus,flag);
end
if ~isempty(rlmod_con)
   disp('reactive load modulation')
   f = rlmod(0,1,bus,flag);
end

% initialize non-linear loads
if ~isempty(load_con)
   disp('non-linear loads')
   vnc = nc_load(bus,flag,Y_ncprf,Y_ncgprf);
else
   nload = 0;
end
if ~isempty(dcsp_con)
   bus_sim = bus;
   bus_int = bus_intprf;
   Y1 = Y_gprf;
   Y2 = Y_gncprf;
   Y3 = Y_ncgprf;
   Y4 = Y_ncprf;
   Vr1 = V_rgprf; 
   Vr2 = V_rncprf;
   bo = boprf;
   h_sol = i_simu(1,1,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);
   % reinitialize dc controls
   f = dc_cont(0,1,bus,flag);
   % initialize dc line
   f = dc_line(0,1,bus,flag);  
end
H_sum = sum(mac_con(:,16)./mac_pot(:,1));
% step 3: perform a predictor-corrector integration 
kt = 0;
ks = 1;
k_tot = sum(k_inc);
lswitch = length(k_inc);
ktmax = k_tot-k_inc(lswitch);
bus_sim = bus;
plot_now = 0;
while (kt<=ktmax)
   k_start = kt+1;
   if kt==ktmax
      k_end = kt + k_inc(ks);
   else
      k_end = kt + k_inc(ks) + 1;
   end
   for k = k_start:k_end
      % step 3a: network solution
      % mach_ref(k) = mac_ang(syn_ref,k);
      mach_ref(k) = 0;
      pmech(:,k+1) = pmech(:,k);
      tmig(:,k+1) = tmig(:,k);
      
      if n_conv~=0;cur_ord(:,k+1) = cur_ord(:,k);end
      
      flag = 1;
      timestep = int2str(k);
      % network-machine interface
      f = mac_ind(0,k,bus_sim,flag);
      f = mac_igen(0,k,bus_sim,flag);
      f = mac_sub(0,k,bus_sim,flag); 
      f = mac_tra(0,k,bus_sim,flag);
      f = mac_em(0,k,bus_sim,flag);
      % Calculate current injections and bus voltages and angles
      if k >= sum(k_inc(1:3))+1
         % fault cleared
         bus_sim = bus_pf2;
         bus_int = bus_intpf2;
         Y1 = Y_gpf2;
         Y2 = Y_gncpf2;
         Y3 = Y_ncgpf2;
         Y4 = Y_ncpf2;
         Vr1 = V_rgpf2; 
         Vr2 = V_rncpf2;
         bo = bopf2;
         h_sol = i_simu(k,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);     
      elseif k >=sum(k_inc(1:2))+1
         % near bus cleared
         bus_sim = bus_pf1;
         bus_int = bus_intpf1;
         Y1 = Y_gpf1;
         Y2 = Y_gncpf1;
         Y3 = Y_ncgpf1;
         Y4 = Y_ncpf1;
         Vr1 = V_rgpf1; 
         Vr2 = V_rncpf1;
         bo = bopf1;
         h_sol = i_simu(k,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);   
      elseif k>=k_inc(1)+1
         % fault applied
         bus_sim = bus_f;
         bus_int = bus_intf;
         Y1 = Y_gf;
         Y2 = Y_gncf;
         Y3 = Y_ncgf;
         Y4 = Y_ncf;
         Vr1 = V_rgf; 
         Vr2 = V_rncf;
         bo = bof;
         h_sol = i_simu(k,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);     
      elseif k<k_inc(1)+1
         % pre fault
         bus_sim = bus;
         bus_int = bus_intprf;
         Y1 = Y_gprf;
         Y2 = Y_gncprf;
         Y3 = Y_ncgprf;
         Y4 = Y_ncprf;