fm_tap.m

一个较好的MATLAB潮流程序

function fm_tap(flag)
% FM_TAP define tap changer on voltage dependent load
%
% FM_TAP(FLAG)
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       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 Tap DAE

% Tap.con:
%          1.  # bus,
%          2.  Power base [MVA],
%          3.  H: integral behaviour deviation,
%          4.  K: Integral time constant,
%          5.  an max: superior tap ratio limit,
%          6.  an min: inferior tap ratio limit,
%          7.  vrif,
%          8.  Pn: nominal active power rate,
%          9.  Qn: nominal reactive power rate,
%         10.  ap: active power exponent,
%         11.  aq: reactive power exponent.

tap = DAE.x(Tap.m);

bustap = Bus.int(Tap.con(:,1));
v = DAE.V(bustap);
hltc = Tap.con(:,4);
kltc = Tap.con(:,5);
mmax = Tap.con(:,6);
mmin = Tap.con(:,7);
vrif = Tap.con(:,8);

switch flag
 case 1

  DAE.gp(bustap) = DAE.gp(bustap) + Tap.con(:,9).*((v./tap).^Tap.con(:,11));
  DAE.gq(bustap) = DAE.gq(bustap) + Tap.con(:,10).*((v./tap).^Tap.con(:,12));

 case 2

  for i=1:Tap.n
    h = bustap(i);
    DAE.J12(h,h)= DAE.J12(h,h) + Tap.con(i,9)*Tap.con(i,11)*(v(i)^(Tap.con(i,11)-1))/(tap(i).^Tap.con(i,11));
    DAE.J22(h,h)= DAE.J22(h,h) + Tap.con(i,10)*Tap.con(i,12)*(v(i)^(Tap.con(i,12)-1))/(tap(i).^Tap.con(i,12));
  end

 case 3

  % calcolo della derivata prima degli stati

  DAE.f(Tap.m) = -hltc.*tap + kltc.*(v./tap - vrif);
  for i = 1:Tap.n
    if (tap(i) >= mmax(i))
      tap(i) = mmax(i);
      if(DAE.f(Tap.m(i)) > 0); DAE.f(Tap.m(i)) = 0;end
    end
    if (tap(i) <= mmin(i))
      tap(i) = mmin(i);
      if(DAE.f(Tap.m(i)) < 0); DAE.f(Tap.m(i)) = 0;end
    end
  end
  DAE.x(Tap.m) = tap;

 case 4

  % calcolo degli jacobiani DAE.Fx, DAE.Fy, DAE.Gx
  for i=1:Tap.n
    k=Tap.m(i);
    DAE.Fx(k,k) = -hltc(i) - v(i)*kltc(i)/(tap(i)*tap(i));
    DAE.Fy(k,bustap(i)+Bus.n)= kltc(i)/tap(i);
    DAE.Gx(bustap(i),k) = -Tap.con(i,9)*Tap.con(i,11)*(v(i)^Tap.con(i,11))/(tap(i)^(Tap.con(i,11)+1));
    DAE.Gx(bustap(i)+Bus.n,k) = -Tap.con(i,10)*Tap.con(i,12)*(v(i)^Tap.con(i,12))/(tap(i)^(Tap.con(i,12)+1));
  end

 case 5
  for i = 1:Tap.n
    if ((tap(i) >= mmax(i) | tap(i) <= mmin(i)) & DAE.f(Tap.m(i)) == 0)
      k = Tap.m(i);
      DAE.tn(k) = 0;
      DAE.Ac(:,k) = zeros(2*Bus.n+DAE.n,1);
      DAE.Ac(k,:) = zeros(1,DAE.n+2*Bus.n);
      DAE.Ac(k,k) = DAE.Fx(k,k);
    end
  end
end