s_simu.m

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

         Vr1 = V_rgprf; 
         Vr2 = V_rncprf;
         bo = boprf;
         h_sol = i_simu(k,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);  
      end
      % HVDC
      f = dc_cont(0,k,bus_sim,flag);
      % network interface for control models
      f = pss(0,k,bus_sim,flag);
      f = mexc_sig(t(k),k);
      f = smpexc(0,k,bus_sim,flag);
      f = exc_st3(0,k,bus_sim,flag);
      f = exc_dc12(0,k,bus_sim,flag);
      f = mtg_sig(t(k),k);
      f = tg(0,k,bus_sim,flag);
      i_plot=k-plot_now;
      if i_plot == 10
         plot_now = k;
         v_p(1:k)=abs(bus_v(bus_idx(1),1:k));
         % plot the voltage of the faulted bus
         plot(t(1:k),v_p(1:k),'r')
         title('Voltage Magnitude at Fault Bus');
         xlabel('time (s)');
         drawnow
      end
      % step 3b: compute dynamics and integrate
      flag = 2;
      sys_freq(k) = 1.0;
	   f = mpm_sig(t(k),k);
      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);
      f = pss(0,k,bus_sim,flag);
      f = mexc_sig(t(k),k);
      f = smpexc(0,k,bus_sim,flag);
      f = exc_st3(0,k,bus_sim,flag);
      f = exc_dc12(0,k,bus_sim,flag);
      f = mtg_sig(t(k),k);
      f = tg(0,k,bus_sim,flag);
      if n_svc~=0
         v_svc = abs(bus_v(bus_int(svc_con(:,2)),k));
         f = msvc_sig(t(k),k);
         f = svc(0,k,bus_sim,flag,v_svc);
      end
      
      if n_lmod~=0
         f = ml_sig(t(k),k);
         f = lmod(0,k,bus_sim,flag);
      end
      if n_rlmod~=0
         f = rml_sig(t(k),k);
         f = rlmod(0,k,bus_sim,flag);
      end

      if n_conv~=0
         f = mdc_sig(t(k),k);
         f = dc_cont(0,k,bus_sim,flag);
         f = dc_line(0,k,bus_sim,flag);
      end
      j = k+1;
      % following statements are predictor steps
      mac_ang(:,j) = mac_ang(:,k) + h_sol*dmac_ang(:,k); 
      mac_spd(:,j) = mac_spd(:,k) + h_sol*dmac_spd(:,k);
      edprime(:,j) = edprime(:,k) + h_sol*dedprime(:,k);
      eqprime(:,j) = eqprime(:,k) + h_sol*deqprime(:,k);
      psikd(:,j) = psikd(:,k) + h_sol*dpsikd(:,k);
      psikq(:,j) = psikq(:,k) + h_sol*dpsikq(:,k);
      Efd(:,j) = Efd(:,k) + h_sol*dEfd(:,k);
      V_R(:,j) = V_R(:,k) + h_sol*dV_R(:,k);
      V_As(:,j) = V_As(:,k) + h_sol*dV_As(:,k);
      R_f(:,j) = R_f(:,k) + h_sol*dR_f(:,k);
      V_TR(:,j) = V_TR(:,k) + h_sol*dV_TR(:,k);
      pss1(:,j) = pss1(:,k) + h_sol*dpss1(:,k);
      pss2(:,j) = pss2(:,k) + h_sol*dpss2(:,k);
      pss3(:,j) = pss3(:,k) + h_sol*dpss3(:,k);
      tg1(:,j) = tg1(:,k) + h_sol*dtg1(:,k);
      tg2(:,j) = tg2(:,k) + h_sol*dtg2(:,k);
      tg3(:,j) = tg3(:,k) + h_sol*dtg3(:,k);
      vdp(:,j) = vdp(:,k) + h_sol*dvdp(:,k);
      vqp(:,j) = vqp(:,k) + h_sol*dvqp(:,k);
      slip(:,j) = slip(:,k) + h_sol*dslip(:,k);
      vdpig(:,j) = vdpig(:,k) + h_sol*dvdpig(:,k);
      vqpig(:,j) = vqpig(:,k) + h_sol*dvqpig(:,k);
      slig(:,j) = slig(:,k) + h_sol*dslig(:,k);
      B_cv(:,j) = B_cv(:,k) + h_sol*dB_cv(:,k);
      B_con(:,j) = B_con(:,k) + h_sol*dB_con(:,k);
      lmod_st(:,j) = lmod_st(:,k) + h_sol*dlmod_st(:,k);
      rlmod_st(:,j) = rlmod_st(:,k)+h_sol*drlmod_st(:,k);
      v_conr(:,j) = v_conr(:,k) + h_sol*dv_conr(:,k);
      v_coni(:,j) = v_coni(:,k) + h_sol*dv_coni(:,k);
      i_dcr(:,j) = i_dcr(:,k) + h_sol*di_dcr(:,k);
      i_dci(:,j) = i_dci(:,k) + h_sol*di_dci(:,k);
      v_dcc(:,j) = v_dcc(:,k) + h_sol*dv_dcc(:,k);
      
      
      flag = 1;
      % mach_ref(j) = mac_ang(syn_ref,j);
      mach_ref(j) = 0;
      % perform network interface calculations again with predicted states
	  f = mpm_sig(t(j),j);
      f = mac_ind(0,j,bus_sim,flag);
      f = mac_igen(0,j,bus_sim,flag);
      f = mac_sub(0,j,bus_sim,flag);
      f = mac_tra(0,j,bus_sim,flag);
      f = mac_em(0,j,bus_sim,flag);
      % assume Vdc remains unchanged for first pass through dc controls interface
      if n_conv~=0;Vdc(:,j) = Vdc(:,k);end
      f = dc_cont(0,j,bus_sim,flag);
      
      % Calculate current injections and bus voltages and angles
      if j >= 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(j,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);     
      elseif j >=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(j,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);   
      elseif j>=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(j,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;
         Vr1 = V_rgprf; 
         Vr2 = V_rncprf;
         bo = boprf;
         h_sol = i_simu(j,ks,k_inc,h,bus_sim,Y1,Y2,Y3,Y4,Vr1,Vr2,bo);     
      end
      vex(:,j)=vex(:,k);
      cur_ord(:,j) = cur_ord(:,k);
      % reapply dc interface to recalculate firing angles with new Vdc
      f = dc_cont(0,j,bus_sim,flag);   
      f = pss(0,j,bus_sim,flag);
      f = mexc_sig(t(j),j);
      f = smpexc(0,j,bus_sim,flag);
      f = exc_st3(0,j,bus_sim,flag);
      f = exc_dc12(0,j,bus_sim,flag);
      f = tg(0,j,bus_sim,flag); 
      
      
      flag = 2;
      f = mac_ind(0,j,bus_sim,flag);
      f = mac_igen(0,j,bus_sim,flag);
      f = mac_sub(0,j,bus_sim,flag);
      f = mac_tra(0,j,bus_sim,flag);
      f = mac_em(0,j,bus_sim,flag);
      f = pss(0,j,bus_sim,flag);
      f = mexc_sig(t(j),j);
      f = smpexc(0,j,bus_sim,flag);
      f = exc_st3(0,j,bus_sim,flag);
      f = exc_dc12(0,j,bus_sim,flag);
      f = mtg_sig(t(j),j);
      f = tg(0,j,bus_sim,flag);
      if n_svc~=0
         f = msvc_sig(t(j),j);
         v_svc = abs(bus_v(bus_int(svc_con(:,2)),j));
         bus_sim = svc(0,j,bus_sim,flag,v_svc);
      end
      
      if n_lmod~=0
         f = ml_sig(t(j),j);
         f = lmod(0,j,bus_sim,flag);
      end
      if n_rlmod~=0
         f = rml_sig(t(j),j);
         f = rlmod(0,j,bus_sim,flag);
      end

      if n_conv~=0
         f = mdc_sig(t(j),j);
         f = dc_cont(0,j,bus_sim,flag);
         f = dc_line(0,j,bus_sim,flag);
      end
      % following statements are corrector steps
      mac_ang(:,j) = mac_ang(:,k) +...
         h_sol*(dmac_ang(:,k)+dmac_ang(:,j))/2.;
      mac_spd(:,j) = mac_spd(:,k) +...
         h_sol*(dmac_spd(:,k)+dmac_spd(:,j))/2.;
      edprime(:,j) = edprime(:,k) +...
         h_sol*(dedprime(:,k)+dedprime(:,j))/2.;
      eqprime(:,j) = eqprime(:,k) +...
         h_sol*(deqprime(:,k)+deqprime(:,j))/2.;
      psikd(:,j) = psikd(:,k) +...
         h_sol*(dpsikd(:,k)+dpsikd(:,j))/2.;
      psikq(:,j) = psikq(:,k) +...
         h_sol*(dpsikq(:,k)+dpsikq(:,j))/2.;
      Efd(:,j) = Efd(:,k) +...
         h_sol*(dEfd(:,k)+dEfd(:,j))/2.;
      V_R(:,j) = V_R(:,k) +...
         h_sol*(dV_R(:,k)+dV_R(:,j))/2.;
      V_As(:,j) = V_As(:,k) +...
         h_sol*(dV_As(:,k)+dV_As(:,j))/2.;
      R_f(:,j) = R_f(:,k) +...
         h_sol*(dR_f(:,k)+dR_f(:,j))/2.;
      V_TR(:,j) = V_TR(:,k) +...
         h_sol*(dV_TR(:,k)+dV_TR(:,j))/2.;
      pss1(:,j) = pss1(:,k) +h_sol*(dpss1(:,k)+dpss1(:,j))/2.;
      pss2(:,j) = pss2(:,k) +h_sol*(dpss2(:,k)+dpss2(:,j))/2.;
      pss3(:,j) = pss3(:,k) +h_sol*(dpss3(:,k)+dpss3(:,j))/2.;
      tg1(:,j) = tg1(:,k) + h_sol*(dtg1(:,k) + dtg1(:,j))/2.;
      tg2(:,j) = tg2(:,k) + h_sol*(dtg2(:,k) + dtg2(:,j))/2.;
      tg3(:,j) = tg3(:,k) + h_sol*(dtg3(:,k) + dtg3(:,j))/2.;
      vdp(:,j) = vdp(:,k) + h_sol*(dvdp(:,j) + dvdp(:,k))/2.;
      vqp(:,j) = vqp(:,k) + h_sol*(dvqp(:,j) + dvqp(:,k))/2.;
      slip(:,j) = slip(:,k) + h_sol*(dslip(:,j) + dslip(:,k))/2.;
      vdpig(:,j) = vdpig(:,k) + h_sol*(dvdpig(:,j) + dvdpig(:,k))/2.;
      vqpig(:,j) = vqpig(:,k) + h_sol*(dvqpig(:,j) + dvqpig(:,k))/2.;
      slig(:,j) = slig(:,k) + h_sol*(dslig(:,j) + dslig(:,k))/2.;
      B_cv(:,j) = B_cv(:,k) + h_sol*(dB_cv(:,j) + dB_cv(:,k))/2.;
      B_con(:,j) = B_con(:,k) + h_sol*(dB_con(:,j) + dB_con(:,k))/2.;
      lmod_st(:,j) = lmod_st(:,k) + h_sol*(dlmod_st(:,j) + dlmod_st(:,k))/2.;
      rlmod_st(:,j) = rlmod_st(:,k) + h_sol*(drlmod_st(:,j) + drlmod_st(:,k))/2.;
      v_conr(:,j) = v_conr(:,k) + h_sol*(dv_conr(:,j) + dv_conr(:,k))/2.;
      v_coni(:,j) = v_coni(:,k) + h_sol*(dv_coni(:,j) + dv_coni(:,k))/2.;
      i_dcr(:,j) = i_dcr(:,k) + h_sol*(di_dcr(:,j) + di_dcr(:,k))/2.;
      i_dci(:,j) = i_dci(:,k) + h_sol*(di_dci(:,j) + di_dci(:,k))/2.;
      v_dcc(:,j) = v_dcc(:,k) + h_sol*(dv_dcc(:,j) + dv_dcc(:,k))/2.;
   end
   kt = kt + k_inc(ks);
   ks = ks+1;
end 
et = toc;
ets = num2str(et);
disp(['elapsed time = ' ets 's'])

flag = 0;
while(flag == 0)
   disp('You can examine the system response')
   disp('Type 1 to see all machine angles in 3D')
   disp('     2 to see all machine speed deviation in 3D')
   disp('     3 to see all machine turbine powers')
   disp('     4 to see all machine electrical powers')
   disp('     5 to see all field voltages')
   disp('     6 to see all bus voltage magnitude in 3D')
   disp('     7 to see the line power flows')
   disp('     0 to quit and plot your own curves')
   sel = input('enter selection >> ');
   if isempty(sel);sel = 0;end
   if sel == 1
      mesh(t,[1:1:n_mac],mac_ang*180/pi)
      title('machine angles')
      xlabel('time in seconds')
      ylabel('internal generator number')
      zlabel('angle in degrees')
      disp('paused: press any key to continue')
      pause
   elseif sel == 2
      lt = length(t);
      mesh(t,[1:1:n_mac],mac_spd- ones(n_mac,lt))
      title('machine speed deviations')
      xlabel('time in seconds')
      ylabel('internal generator number')
      zlabel('speed in pu')
      disp('paused: press any key to continue')
      pause
   elseif sel == 3
      plot(t,pmech)
      title('turbine power')
      xlabel('time in seconds')
      ylabel('power in pu on generator base')
      disp('paused: press any key to continue')
      pause
   elseif sel == 4
      plot(t,pelect)
      title('generator electrical power')
      xlabel('time in seconds')
      ylabel('power in pu on system base')
      disp('paused: press any key to continue')
      pause
   elseif sel == 5
      plot(t,Efd)
      title('Generator field voltages')
      xlabel('time in seconds')
      ylabel('voltage in pu')
      disp('paused: press any key to continue')
      pause
   elseif sel == 6
      [nbus,dum] = size(bus_v);
      mesh(t,[1:1:nbus],abs(bus_v))
      xlabel('time in seconds')
      ylabel('bus number')
      zlabel('voltage in pu')
      disp('paused: press any key to continue')
      pause
   elseif sel == 7
      nline = length(line(:,1));
      if nline<50
         V1 = bus_v(bus_int(line(:,1)),:);
         V2 = bus_v(bus_int(line(:,2)),:);
         R = line(:,3); X = line(:,4); B = line(:,5);
         tap = line(:,6); phi = line(:,7);
         [S1,S2] = line_pq(V1,V2,R,X,B,tap,phi);
         plot(t,real(S1));
         xlabel('time in seconds')
         ylabel(' power flow per unit')
         disp('paused: press any key to continue')
         pause
      else
         % ask for lines to be plotted
         line_range = input(' enter a rangle of lines less than 50, ie 1:49   >');
         if isempty(line_range);line_range = 1:49;end
         V1 = bus_v(bus_int(line(line_range,1)),:);
         V2 = bus_v(bus_int(line(line_range,2)),:);
         R = line(line_range,3); X = line(line_range,4); 
         B = line(line_range,5);
         tap = line(line_range,6); phi = line(line_range,7);
         [S1,S2] = line_pq(V1,V2,R,X,B,tap,phi);
         plot(t,real(S1));
         xlabel('time in seconds')
         ylabel(' power flow per unit')
         disp('paused: press any key to continue')
         pause
      end
   elseif sel == 0
      flag = 1;
   else
      error('invalid selection')
   end
end