fm_int.m

这是一个很适合研究和学习用的电力系统仿真软件

function  fm_int
% FM_INT time domain integration routines:
%       1 - Forward Euler
%       2 - Trapezoidal Method
%
% FM_INT
%
%see also FM_TSTEP, FM_OUT and the Settings structure
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    16-Jan-2003
%Update:    27-Feb-2003
%Update:    01-Aug-2003
%Update:    11-Sep-2003
%Version:   1.0.4
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Fig Settings Snapshot Hdl
global Bus File DAE Theme
global SW PV PQ Fault Mn Ltc Syn Exc Oxl
global Varout Breaker Line Path clpsat

if ~autorun('Time Domain Simulation',1), return, end

tic

% initial messages
% -----------------------------------------------------------------------

fm_disp

if DAE.n == 0 & ~clpsat.init
  Settings.ok = 0;
  uiwait(fm_choice('No dynamic component is loaded. Continue anyway?',1))
  if ~Settings.ok
    fm_disp('Time domain simulation aborted.',2)
    return
  end
end

fm_disp('Time domain simulation')
switch Settings.method
 case 1,
  fm_disp('Implicit Euler integration method')
 case 2,
  fm_disp('Trapezoidal integration method')
end
fm_disp(['Data file "',Path.data,File.data,'"'])
if ~isempty(Path.pert),
  fm_disp(['Perturbation file "',Path.pert,File.pert,'"'])
end
if (strcmp(File.pert,'pert') & strcmp(Path.pert,Path.psat)) | ...
      isempty(File.pert)
  fm_disp('No perturbation file set.',1)
end

if Fig.main
  hdl = findobj(Fig.main,'Tag','PushClose');
  set(hdl,'String','Stop');
  set(Fig.main,'UserData',1);
end

if Settings.plot
  if Settings.plottype == 1 & ~DAE.n
    Settings.plottype = 2;
    fm_disp('Cannot plot state variables (dynamic order = 0).')
    fm_disp('Bus voltages will be plotted during the TD simulation.')
  end
  maxlegend = min(Bus.n,7);
  switch Settings.plottype
   case 1
    maxlegend = min(DAE.n,7);
    idx0 = 0;
   case 2
    idx0 = DAE.n;
   case 3
    idx0 = DAE.n+Bus.n;
   case 4
    idx0 = DAE.n+2*Bus.n+2*Syn.n+Exc.n+Oxl;
   case 5
    idx0 = DAE.n+3*Bus.n+2*Syn.n+Exc.n+Oxl;
   case 6
    maxlegend = 3;
    idx0 = DAE.n+4*Bus.n+2*Syn.n+Exc.n+Oxl;
  end
end

% check settings
% ------------------------------------------------------------------

iter_max = Settings.dynmit;
tol = Settings.dyntol;
Dn = 1;
if DAE.n, Dn = DAE.n; end
identica = speye(max(Dn,1));

if (Fault.n | Breaker.n) & PQ.n & ~Settings.pq2z
  if clpsat.init
    if clpsat.pq2z, Settings.pq2z = 1; end
  else
    uiwait(fm_choice(['Convert (recommended) PQ loads to constant impedances?']))
    if Settings.ok, Settings.pq2z = 1; end
  end
end

% convert PQ loads to shunt admittances (if required)
Line.Y = Line.Y + pqshunt(PQ);
if Settings.pq2z, PQ = remove(PQ,0,'onlyload'); end

% set up variables
% ----------------------------------------------------------------

DAE.t = Settings.t0;
fm_call('i');
DAE.tn = DAE.f;
if isempty(DAE.tn), DAE.tn = 0; end
DAE.g = [DAE.gp; DAE.gq];

% graphical settings
% ----------------------------------------------------------------

plot_now = 0;
if ~clpsat.init | Fig.main
  if Settings.plot
    if Fig.plot
      figure(Fig.plot);
    else
      fm_plotfig;
    end
  elseif Settings.status
    fm_bar('open')
    fm_simtd('init')
    idxo = 0;
  else
    fm_disp(['t = ',num2str(Settings.t0),' s'],3)
    perc = 0;
    perc_old = 0;
  end
  drawnow
end

% ----------------------------------------------------------------
% initializations
% ----------------------------------------------------------------

t = Settings.t0;
k = 1;
h = fm_tstep(1,1,0,Settings.t0);
if Fault.n ~= 0, fm_fault(0,0), end
fm_breaker(0,0)
inc = zeros(Dn + 2*Bus.n,1);
callpert = 1;

% output initialization
% ----------------------------------------------------------------

fm_out(0,0,0);
fm_out(2,Settings.t0,k);

% time vector of snapshots, faults and breaker events
% ----------------------------------------------------------------

fixed_times = [];

n_snap = length(Snapshot);
if n_snap > 1
  snap_times = zeros(n_snap-1,1);
  for i = 2:n_snap
    snap_times(i-1,1) = Snapshot(i).time;
  end
  fixed_times = [fixed_times; snap_times];
end

if Fault.n
  idx = find(Fault.con(:,5) ~= Fault.con(:,6));
  if ~isempty(idx)
    fixed_times = [fixed_times; ...
                   Fault.con(idx,5)-1e-6; Fault.con(idx,5); ...
                   Fault.con(idx,6)-1e-6; Fault.con(idx,6)];
  end
end

if Breaker.n
  idx = find(Breaker.con(:,7) ~= Breaker.con(:,8));
  if ~isempty(idx)
    fixed_times = [fixed_times; ...
                   Breaker.con(idx,7)-1e-6; Breaker.con(idx,7); ...
                   Breaker.con(idx,8)-1e-6; Breaker.con(idx,8)];
  end
end

fixed_times = sort(fixed_times);

% compute max rotor angle difference
diff_max = anglediff;

% ================================================================
% ----------------------------------------------------------------
% Main loop
% ----------------------------------------------------------------
% ================================================================

inc = zeros(Dn+2*Bus.n,1);

while (t < Settings.tf) & (t + h > t) & ~diff_max
  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end
  if (t + h > Settings.tf), h = Settings.tf - t; end
  actual_time = t + h;

  % check not to jump disturbances
  index_times = find(fixed_times > t & fixed_times < t+h);
  if ~isempty(index_times);
    actual_time = min(fixed_times(index_times));
    h = actual_time - t;
  end

  % set global time
  DAE.t = actual_time;

  % backup of actual variables
  if isempty(DAE.x), DAE.x = 0; end
  xa = DAE.x;
  anga = DAE.a;
  Va = DAE.V;

  % initialize NR loop
  iterazione = 1;
  inc(1) = 1;
  if isempty(DAE.f), DAE.f = 0; end
  fn = DAE.f;

  % compute perturbations
  if ~isempty(fixed_times)
    if ~isempty(find(fixed_times == actual_time))
      if Fault.n
        fm_fault(1,actual_time)
      end
      if Breaker.n
        fm_breaker(1,actual_time)
      end
    end
  end

  if callpert
    try
      if Settings.hostver >= 6
        feval(Hdl.pert,actual_time);
      else
        if ~isempty(Path.pert)
          cd(Path.pert)
          feval(Hdl.pert,actual_time);
          cd(Path.local)
        end
      end
    catch
      fm_disp('* * Something wrong in the perturbation file:')
      fm_disp(lasterr)
      fm_disp('* * The perturbation file will be discarded.')      
      callpert = 0;
    end
  end
      
  % Newton-Raphson loop
  while max(abs(inc)) > tol
    if (iterazione > iter_max), break,  end
    drawnow
    if Fig.main
      if ~get(Fig.main,'UserData'), break, end
    end
    if isempty(Line.Y)
      DAE.gp = zeros(Bus.n,1);
      DAE.gq = zeros(Bus.n,1);
      DAE.J11 = sparse(Bus.n,Bus.n);
      DAE.J21 = sparse(Bus.n,Bus.n);
      DAE.J12 = sparse(Bus.n,Bus.n);
      DAE.J22 = sparse(Bus.n,Bus.n);
    end

    % DAE equations
    fm_call('i');
    DAE.Fy(:,SW.bus) = 0;
    DAE.Gx(SW.bus,:) = 0;
    DAE.Fy(:,Bus.n+SW.bus) = 0;
    DAE.Gx(Bus.n+SW.bus,:) = 0;
    DAE.Fy(:,Bus.n+PV.bus) = 0;
    DAE.Gx(Bus.n+PV.bus,:) = 0;

    % check for islanded buses
    if ~isempty(Bus.island)
      kkk = Bus.island;
      DAE.Jlfv(kkk,:) = 0;
      DAE.Jlfv(:,kkk) = 0;
      DAE.Jlfv(:,kkk+Bus.n) = 0;
      DAE.Jlfv(kkk+Bus.n,:) = 0;
      DAE.Jlfv(kkk,kkk) = speye(length(kkk));
      DAE.Jlfv(kkk+Bus.n,kkk+Bus.n) = speye(length(kkk));
      DAE.g(kkk) = 0;
      DAE.g(kkk+Bus.n) = 0;
      DAE.V(kkk) = 1e-6;
      DAE.a(kkk) = 0;
    end

    % complete Jacobian matrix DAE.Ac
    switch Settings.method
     case 1  % Forward Euler
      DAE.Ac = [identica - h*DAE.Fx, -h*DAE.Fy; DAE.Gx, DAE.Jlfv];
      DAE.tn = DAE.x - xa - h*DAE.f;
     case 2  % Trapezoidal Method
      DAE.Ac = [identica - h*0.5*DAE.Fx, -h*0.5*DAE.Fy; DAE.Gx, DAE.Jlfv];
      DAE.tn = DAE.x - xa - h*0.5*(DAE.f + fn);
    end

    % non-windup limiters
    fm_call('5');
    inc = -DAE.Ac\[DAE.tn; DAE.g];

    %inc = -umfpack(DAE.Ac,'\',[DAE.tn; DAE.g]);
    DAE.x = DAE.x + inc(1:Dn);
    DAE.a = DAE.a + inc(1+Dn: Bus.n+Dn);
    DAE.V = DAE.V + inc(Dn+Bus.n+1: Dn+2*Bus.n);
    iterazione = iterazione + 1;
  end

  if (iterazione > iter_max)
    h = fm_tstep(2,0,iterazione,t);
    DAE.x = xa;
    DAE.a = anga;
    DAE.V = Va;
    DAE.f = fn;
  else
    h = fm_tstep(2,1,iterazione,t);
    t = actual_time;
    k = k+1;
    % extend output stack
    if k > length(Varout.t), fm_out(1,t,k); end

    % ----------------------------------------------------------------
    % ----------------------------------------------------------------
    % update output variables and snapshots
    % ----------------------------------------------------------------
    % ----------------------------------------------------------------

    fm_out(2,t,k);

    %  plot variables & display iteration status
    % ----------------------------------------------------------------

    i_plot = 1+k-10*fix(k/10);
    perc = (t-Settings.t0)/(Settings.tf-Settings.t0);
    if i_plot == 10
      fm_disp([' > Simulation time = ',num2str(DAE.t), ...
	       ' s (',num2str(round(perc*100)),'%)'])
    end

    if ~clpsat.init | Fig.main
      if Settings.plot
        if i_plot == 10
          plot(Varout.t(1:k),Varout.vars(1:k,idx0+[1:maxlegend]));
          set(gca,'Color',Theme.color11);
          xlabel('time (s)')
          drawnow
        end
      elseif Settings.status
        idx = (t-Settings.t0)/(Settings.tf-Settings.t0);
        fm_bar([idxo,idx])
        if i_plot == 10, fm_simtd('update'), end
        idxo = idx;
      end
    end

    % fill up snapshots
    % ------------------------------------------------------------------

    if n_snap > 1
      snap_i = find(snap_times == t)+1;
      if ~isempty(snap_i)
        flag = 1;
        DAE.gp = zeros(Bus.n,1);
        DAE.gq = zeros(Bus.n,1);
        fm_call('1');
        Pl = DAE.gp;
        Ql = DAE.gq;
        P = DAE.glfp;
        Q = DAE.glfq;
        Snapshot(snap_i).V = DAE.V;
        Snapshot(snap_i).ang = DAE.a;
        Snapshot(snap_i).x = DAE.x;
        Snapshot(snap_i).Y = Line.Y;
        Snapshot(snap_i).Pg = P + Pl;
        Snapshot(snap_i).Qg = Q + Ql;
        Snapshot(snap_i).Pl = DAE.gp;
        Snapshot(snap_i).Ql = DAE.gq;
        Snapshot(snap_i).vfd = Syn.vf;
        Snapshot(snap_i).pmech = Syn.pm;
        Snapshot(snap_i).Jlf = DAE.Jlf;
        Snapshot(snap_i).Jlfv = DAE.Jlfv;
        Snapshot(snap_i).Fx = DAE.Fx;
        Snapshot(snap_i).Fy = DAE.Fy;
        Snapshot(snap_i).Gx = DAE.Gx;
      end
    end
  end

  % compute max rotor angle difference
  diff_max = anglediff;

end

if Settings.status & ~Settings.plot
  fm_bar('close')
  fm_simtd('update')
end
if ~DAE.n, DAE.x = []; DAE.f =[]; end

% final messages
% -----------------------------------------------------------------------

if Fig.main
  if diff_max & get(Fig.main,'UserData')
    fm_disp(['Rotor angle max difference is > ', ...
             num2str(Settings.deltadelta), ...
             ' deg. Simulation stopped at t = ', ...
             num2str(t), ' s'],2);
  elseif (t < Settings.tf) & get(Fig.main,'UserData')
    fm_disp(['Singularity likely. Simulation stopped at t = ', ...
             num2str(t), ' s'],2);
  elseif ~get(Fig.main,'UserData')
    fm_disp(['Dynamic Simulation interrupted at t = ',num2str(t),' s'],2)
  else
    fm_disp(['Dynamic Simulation completed in ',num2str(toc),' s']);
  end
else
  if diff_max
    fm_disp(['Rotor angle max difference is > ', ...
             num2str(Settings.deltadelta), ...
             ' deg. Simulation stopped at t = ', ...
             num2str(t), ' s'],2);

  elseif (t < Settings.tf)
    fm_disp(['Singularity likely. Simulation stopped at t = ', ...
             num2str(t), ' s'],2);
  else
    fm_disp(['Dynamic Simulation completed in ',num2str(toc),' s']);
  end
end

% resize output varibales & final settings
% -----------------------------------------------------------------------

fm_out(3,t,k);
if Settings.beep, beep, end
Settings.xlabel = 'time (s)';
if Fig.plot, fm_plotfig, end

% future simulations do not need LF computation
Settings.init = 2;
SNB.init = 0;
LIB.init = 0;
CPF.init = 0;
OPF.init = 0;
if Fig.main, set(hdl,'String','Close'); end
DAE.t = -1; % reset global time

% compute delta difference at each step
% -----------------------------------------------------------------------
function diff_max = anglediff
global Settings Syn Bus DAE

diff_max = 0;

if ~Settings.checkdelta, return, end
if ~Syn.n, return, end

delta = DAE.x(Syn.delta);
idx = [];
for h = 1:length(Bus.island)
  idx = [idx; find(Syn.bus == Bus.island(h))];
end
if ~isempty(idx)
  delta(idx) = [];
end

delta_diff = abs(delta-min(delta));

diff_max = (max(delta_diff)*180/pi) > Settings.deltadelta;