fm_tg.m

一个较好的MATLAB潮流程序

function fm_tg(flag)
%FM_TG define Turbine Governor
%
%FM_TG(FLAG)
%       FLAG = 0 initializations
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global Bus Tg DAE Syn

type_tg = Tg.con(:,2);
ty1 = find(type_tg == 1);
ty2 = find(type_tg == 2);

if ty1
  tg1 = DAE.x(Tg.tg1(ty1));
  tg2 = DAE.x(Tg.tg2(ty1));
  tg3 = DAE.x(Tg.tg3(ty1));
end
if ty2
  tg = DAE.x(Tg.tg(ty2));
end

switch flag
 case 0     % initialization

  Porder = Syn.pm(Tg.syn);
  M = Syn.con(:,18);

  if (~isempty(ty1))

    Porder1 = Porder(ty1);
    gain = 1./Tg.con(ty1,4);
    tmax = Tg.con(ty1,5);
    tmin = Tg.con(ty1,6);
    Ts = Tg.con(ty1,7);
    Tc = Tg.con(ty1,8);
    T3 = Tg.con(ty1,9);
    T4 = Tg.con(ty1,10);
    T5 = Tg.con(ty1,11);

    as = 1./Ts;
    ac = 1./Tc;
    a5 = 1./T5;
    K1 = T3.*ac;
    K2 = 1 - K1;
    K3 = T4.*a5;
    K4 = 1 - K3;

    Tg.dat1 = [Tg.tg1(ty1), ...              %  1
               Tg.tg2(ty1), ...              %  2
               Tg.tg3(ty1), ...              %  3
               gain, ...                     %  4
               tmax, ...                     %  5
               Porder1, ...                  %  6
               as, ...                       %  7
               ac, ...                       %  8
               a5, ...                       %  9
               K1, ...                       % 10
               K2, ...                       % 11
               K3, ...                       % 12
               K4, ...                       % 13
               Syn.omega(Tg.syn(ty1)), ...   % 14
               1./M(Tg.syn(ty1))];           % 15

    DAE.x(Tg.dat1(:,1)) = Tg.dat1(:,6);
    DAE.x(Tg.dat1(:,2)) = Tg.dat1(:,11).*Tg.dat1(:,6);
    DAE.x(Tg.dat1(:,3)) = Tg.dat1(:,13).*Tg.dat1(:,6);
    DAE.f(Tg.dat1(:,1)) = 0;
    DAE.f(Tg.dat1(:,2)) = 0;
    DAE.f(Tg.dat1(:,3)) = 0;

  end

  if (~isempty(ty2))

    Porder2 = Porder(ty2);
    gain = 1./Tg.con(ty2,4);
    tmax = Tg.con(ty2,5);
    tmin = Tg.con(ty2,6);
    T1 = Tg.con(ty2,8);
    T2 = Tg.con(ty2,7);
    a2 = 1./T2;
    K1 = gain.*T1.*a2;
    K2 = gain - K1;

    Tg.dat2 = [Tg.tg(ty2), ...               %  1
               gain, ...                     %  2
               tmax, ...                     %  3
               tmin, ...                     %  4
               Porder2, ...                  %  5
               a2, ...                       %  6
               K1, ...                       %  7
               K2, ...                       %  8
               Syn.omega(Tg.syn(ty2)), ...   %  9
               1./M(Tg.syn(ty2))];           % 10

    DAE.x(Tg.dat2(:,1)) = 0;
    DAE.f(Tg.dat2(:,1)) = 0;
  end

  maxlmt = find(Porder > Tg.con(:,5));
  minlmt = find(Porder < Tg.con(:,6));
  for i = 1:length(maxlmt)
    fm_disp(['Pmech greater than maximum limit for TG #', ...
             num2str(maxlmt(i))],1)
  end
  for i = 1:length(minlmt)
    fm_disp(['Pmech less than minimum limit for TG #', ...
             num2str(minlmt(i))],1)
  end
  if isempty(maxlmt) & isempty(minlmt)
    fm_disp('Initialization of Turbine Gorvernors completed.',1)
  else
    fm_disp('Initialization of Turbine Gorvernors failed.',2)
    return
  end

 case 3   % differential equations

  if ~isempty(ty1)
    tin = Tg.dat1(:,6) + Tg.dat1(:,4).*(Tg.con(ty1,3)-DAE.x(Tg.dat1(:,14)));
    tin = max(tin,Tg.con(ty1,6));
    tin = min(tin,Tg.con(ty1,5));
    Tg.dat1(:,16) = tin;
    DAE.f(Tg.dat1(:,1)) = Tg.dat1(:,7).*(-tg1 + tin);
    DAE.f(Tg.dat1(:,2)) = Tg.dat1(:,8).*(-tg2 + Tg.dat1(:,11).*tg1);
    DAE.f(Tg.dat1(:,3)) = Tg.dat1(:,9).*(-tg3 + Tg.dat1(:,13).*(tg2 + Tg.dat1(:,10).*tg1));
    Syn.pm(Tg.syn(ty1)) = tg3 + Tg.dat1(:,12).*(tg2 + Tg.dat1(:,10).*tg1);
  end

  if ~isempty(ty2)
    DAE.f(Tg.dat2(:,1)) = -Tg.dat2(:,6).*(tg + Tg.dat2(:,8).*(DAE.x(Tg.dat2(:,9)) - Tg.con(ty2,3)));
    Syn.pm(Tg.syn(ty2)) = tg - Tg.dat2(:,7).*(DAE.x(Tg.dat2(:,9)) - Tg.con(ty2,3)) + Tg.dat2(:,5);
    Syn.pm(Tg.syn(ty2)) = max(Syn.pm(Tg.syn(ty2)), Tg.dat2(:,4));
    Syn.pm(Tg.syn(ty2)) = min(Syn.pm(Tg.syn(ty2)), Tg.dat2(:,3));
  end

 case 4   % Jacobians

  if ~isempty(ty1)
    tin = Tg.dat1(:,16);
    idx = find(tin < Tg.con(ty1,5) & tin > Tg.con(ty1,6));  % windup limiter
    if ~isempty(idx),
      DAE.Fx = DAE.Fx + sparse(Tg.dat1(idx,1),Tg.dat1(idx,14),-Tg.dat1(idx,7).*Tg.dat1(idx,4),DAE.n,DAE.n);
    end
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,14),Tg.dat1(:,1),Tg.dat1(:,12).*Tg.dat1(:,15).*Tg.dat1(:,10),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,14),Tg.dat1(:,2),Tg.dat1(:,12).*Tg.dat1(:,15),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,14),Tg.dat1(:,3),Tg.dat1(:,15),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,1),Tg.dat1(:,1),-Tg.dat1(:,7),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,2),Tg.dat1(:,2),-Tg.dat1(:,8),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,3),Tg.dat1(:,3),-Tg.dat1(:,9),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,2),Tg.dat1(:,1),Tg.dat1(:,8).*Tg.dat1(:,11),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,3),Tg.dat1(:,2),Tg.dat1(:,9).*Tg.dat1(:,13),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat1(:,3),Tg.dat1(:,1),Tg.dat1(:,10).*Tg.dat1(:,13).*Tg.dat1(:,9),DAE.n,DAE.n);
  end

  if ~isempty(ty2)
    tin = tg - Tg.dat2(:,7).*(DAE.x(Tg.dat2(:,9)) - Tg.con(ty2,3)) + Tg.dat2(:,5);
    idx = find(tin < Tg.dat2(:,3) & tin > Tg.dat2(:,4)); % windup limiter
    if ~isempty(idx),
      DAE.Fx = DAE.Fx + sparse(Tg.dat2(idx,9),Tg.dat2(idx,9),-Tg.dat2(idx,7).*Tg.dat2(idx,10),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Tg.dat2(idx,9),Tg.dat2(idx,1),Tg.dat2(idx,10),DAE.n,DAE.n);
    end
    DAE.Fx = DAE.Fx + sparse(Tg.dat2(:,1),Tg.dat2(:,1),-Tg.dat2(:,6),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tg.dat2(:,1),Tg.dat2(:,9),-Tg.dat2(:,6).*Tg.dat2(:,8),DAE.n,DAE.n);
  end

 case 5     %  non-windup limiters

  %  void

end