fm_cluster.m

电力系统的psat

function  fm_cluster(flag)
% FM_CLUSTER define Cluster Controllers
%
% FM_CLUSTER(FLAG)
%       FLAG = 0 -> initialization
%       FLAG = 1 -> algebraic equations
%       FLAG = 2 -> algebraic Jacobians
%       FLAG = 3 -> differential equations
%       FLAG = 4 -> state Jacobians
%       FLAG = 5 -> non-windup limiters)
%
%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 Cluster DAE Bus Svc Exc Syn Cac

%Data Format Cluster.con:
%          col #1: AVR or SVC number
%          col #2: type (1) AVR, (2) SVC
%          col #3: Central Area Control number
%          col #4: T Integral time constant [Second]
%          col #5: Xtg None [None]
%          col #6: Xeq None [None]
%          col #7: Qgr None [None]
%          col #8: Vs_max
%          col #9: Vs_max

Vs = DAE.x(Cluster.Vs);
T = Cluster.con(:,4);
Xtg = Cluster.con(:,5);
Xeq = Cluster.con(:,6);
Qgr = Cluster.con(:,7);
Vs_max = Cluster.con(:,8);
Vs_min = Cluster.con(:,9);

switch flag
 case 0 % initialization

  % check time constants
  idx = find(T == 0);
  if idx
    clusterwarn(idx, ['Time constant T cannot be zero. T = 0.001 s ' ...
                      'will be used.'])
  end
  Cluster.con(idx,5) = 0.001;

  % variable initialization
  DAE.x(Cluster.Vs) = 0;
  Vs = DAE.x(Cluster.Vs);

  % set reactive power references
  % (the CAC signal q is q = 1 at the steady state e.q.)
  type1 = find(Cluster.con(:,2) == 1);
  type2 = find(Cluster.con(:,2) == 2);
  type3 = find(Cluster.con(:,2) == 3);
  Qg = zeros(Cluster.n,1);
  if type1
    Qg(type1) = -Syn.Qg(Cluster.syn(type1));
  end
  if type2
    Qg(type2) = Svc.Be(Cluster.svc(type2)).* ...
        DAE.V(Cluster.bus(type2)).*DAE.V(Cluster.bus(type2));
  end
  Cluster.con(:,7) = Qg;

  % check limits
  idx = find(Vs > Vs_max);
  if idx
    clusterwarn(idx, [' State variable Vs is over its maximum ' ...
                      'limit.'])
  end
  idx = find(Vs < Vs_min);
  if idx
    clusterwarn(idx, [' State variable Vs is under its minimum ' ...
                      'limit.'])
  end
  fm_disp('Initialization of Cluster''s completed.')

 case 1 % algebraic equations


 case 2 % algebraic Jacobians


 case 3 % differential equations

  type1 = find(Cluster.con(:,2) == 1);
  type2 = find(Cluster.con(:,2) == 2);
  Qg = zeros(Cluster.n,1);
  if type1
    Qg(type1) = -Syn.Qg(Cluster.syn(type1));
  end
  if type2
    Qg(type2) = Svc.Be(Cluster.svc(type2)).* ...
        DAE.V(Cluster.bus(type2)).*DAE.V(Cluster.bus(type2));
  end
  q = Cac.q(Cluster.cac);
  DAE.f(Cluster.Vs) = ((Qgr.*q-Qg).*(Xtg+Xeq))./T;
  % non-windoup limits
  idx = find(Vs >= Vs_max & DAE.f(Cluster.Vs) > 0);
  if idx, DAE.f(Cluster.Vs(idx)) = 0; end
  DAE.x(Cluster.Vs) = min(Vs,Vs_max);
  idx = find(Vs <= Vs_min & DAE.f(Cluster.Vs) < 0);
  if idx, DAE.f(Cluster.Vs(idx)) = 0; end
  DAE.x(Cluster.Vs) = max(Vs,Vs_min);

 case 4 % state variable Jacobians

  q1 = DAE.x(Cac.q1(Cluster.cac));
  KP = Cac.con(Cluster.cac,7);
  q1_max = Cac.con(Cluster.cac,8);
  q1_min = Cac.con(Cluster.cac,9);
  idx = find((q1 >= q1_max | q1 <= q1_min) & DAE.f(Cac.q1) == 0);
  if idx
    DAE.Fx = DAE.Fx + ...
             sparse(Cluster.Vs(idx),Cac.q1(Cluster.cac(idx)), ...
                    Qgr(idx).*(Xtg(idx)+Xeq(idx))./T(idx),DAE.n,DAE.n);
  end
  DAE.Fy = DAE.Fy + sparse(Cluster.Vs,Cac.bus(Cluster.cac)+Bus.n, ...
                           -KP.*Qgr.*(Xtg+Xeq)./T,DAE.n,2*Bus.n);

  type = Cluster.con(:,2);
  dVs_dQ = (Xtg+Xeq)./T;
  idx = find(Vs <= Vs_max & Vs >= Vs_min);

  for i = 1:length(idx)

    switch type(idx(i))
     case 1 % the Cluster is connected to an AVR

      h = Cluster.exc(idx(i));
      j = Cluster.bus(idx(i));
      m = Cluster.syn(idx(i));
      k = Cluster.Vs(idx(i));

      % update Jacobians of Synchronous Machines & AVRs
      type_exc = Exc.con(h,2);
      switch type_exc
       case 1
        DAE.Fx(Exc.vr1(h),k) = -DAE.Fx(Exc.vr1(h), Exc.vm(h));
        DAE.Fx(Exc.vr2(h),k) = -DAE.Fx(Exc.vr2(h), Exc.vm(h));
       case 2
        DAE.Fx(Exc.vr1(h),k) = -DAE.Fx(Exc.vr1(h), Exc.vm(h));
       case 3
        vrmax = Exc.con(h,3);
        vrmin = Exc.con(h,4);
        if Syn.vf(m) < vrmax & Syn.vf(m) > vrmin
          Td10 = Syn.con(m,11);
          Td20 = Syn.con(m,12);
          Taa  = Syn.con(m,24);
          xd   = Syn.con(m,8);
          x1d  = Syn.con(m,9);
          xL   = Syn.con(m,6);
          ord  = Syn.con(m,5);
          vm3  = DAE.x(Exc.vm(h));
          Kr   = Exc.con(h,5);
          T2r  = Exc.con(h,6);
          T1r  = Exc.con(h,7);
          Kr1  = Kr*T1r/T2r;
          v0   = Exc.con(h,9);
          k1   = -Kr1*vm3/v0;
          switch ord
           case 3,    DAE.Fx(Syn.e1q(m),k) = -k1/Td10;
           case 4,    DAE.Fx(Syn.e1q(m),k) = -k1/Td10;
           case 5.1,  DAE.Fx(Syn.e1q(m),k) = -k1/Td10;
           case 5.3,  DAE.Fx(Syn.e1q(m),k) = -k1*(xd+xl)/(x1d+xL)/Td10;
           otherwise,
            DAE.Fx(Syn.e1q(m),k) = -k1*(1-Taa/Td10)/Td10;
            DAE.Fx(Syn.e2q(m),k) = -k1*Taa/Td10/Td20;
          end
        end
        DAE.Fx(Exc.vr3(h),k) = -DAE.Fx(Exc.vr3(h), Exc.vm(h));
      end

      % Cluster Control Jacobian of algebraic variables
      % (ineherited from synchronous machines)
      DAE.Fy(k,j) = dVs_dQ(idx(i))*Syn.J21(m);
      DAE.Fy(k,j+Bus.n) = dVs_dQ(idx(i))*Syn.J22(m);

      % Cluster Jacobian of state variables
      % (inherited from synchronous machines)
      DAE.Fx(k,Syn.delta(m)) = dVs_dQ(idx(i))*Syn.Gq(m,1);

      type_syn = Syn.con(Cluster.syn(idx(i)),5);
      switch type_syn
       case 3
        DAE.Fx(k,Syn.e1q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,3);
       case 4
        DAE.Fx(k,Syn.e1q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,3);
        DAE.Fx(k,Syn.e1d(m))  = dVs_dQ(idx(i))*Syn.Gq(m,4);
       case 5.1
        DAE.Fx(k,Syn.e1q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,3);
        DAE.Fx(k,Syn.e2d(m))  = dVs_dQ(idx(i))*Syn.Gq(m,4);
       case 5.2
        DAE.Fx(k,Syn.e2q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,3);
        DAE.Fx(k,Syn.e2d(m))  = dVs_dQ(idx(i))*Syn.Gq(m,4);
       case 5.3
        DAE.Fx(k,Syn.e1q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,3);
        DAE.Fx(k,Syn.psiq(m)) = dVs_dQ(idx(i))*Syn.Gq(m,7);
        DAE.Fx(k,Syn.psid(m)) = dVs_dQ(idx(i))*Syn.Gq(m,8);
       case 6
        DAE.Fx(k,Syn.e2q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,3);
        DAE.Fx(k,Syn.e2d(m))  = dVs_dQ(idx(i))*Syn.Gq(m,4);
       case 8
        DAE.Fx(k,Syn.e2q(m))  = dVs_dQ(idx(i))*Syn.Gq(m,5);
        DAE.Fx(k,Syn.e2d(m))  = dVs_dQ(idx(i))*Syn.Gq(m,6);
        DAE.Fx(k,Syn.psiq(m)) = dVs_dQ(idx(i))*Syn.Gq(m,7);
        DAE.Fx(k,Syn.psid(m)) = dVs_dQ(idx(i))*Syn.Gq(m,8);
      end

     case 2 % Cluster Controller  is connected to an SVC

      h = Cluster.svc(idx(i));
      j = Cluster.bus(idx(i));
      k = Cluster.Vs(idx(i));

      % Cluster Control Jacobian of algebraic variables (ineherited from SVC)
      DAE.Fy(k,j+Bus.n) = -2*dVs_dQ(idx(i))*Svc.Be(h).*DAE.V(j);

      % Cluster Jacobian of state variables (inherited from SVC)
      type_svc = Svc.con(Cluster.svc(idx(i)),5);

      switch type_svc
       case 1
        bcv = DAE.x(Svc.bcv(h));
        Tr = Svc.con(h,6);
        Kr = Svc.con(h,7);
        bcv_max = Svc.con(h,9);
        bcv_min = Svc.con(h,10);
        if (bcv <= bcv_max & bcv >= bcv_min) & DAE.f(Svc.bcv(h)) ~= 0
          DAE.Fx(Svc.bcv(h),k) = Kr/Tr;
          DAE.Fx(k,Svc.bcv(h)) = -dVs_dQ(idx(i))*DAE.V(j)^2;
        end
       case 2
        alpha = DAE.x(Svc.alpha(h));
        a_max = Svc.con(h,9);
        a_min = Svc.con(h,10);
        %Vref2 = Svc.con(h,8);
        T2 = Svc.con(h,6);
        K = Svc.con(h,7);
        Kd = Svc.con(h,11);
        T1 = Svc.con(h,12);
        Km = Svc.con(h,13);
        Tm = Svc.con(h,14);
        xl = Svc.con(h,15);
        xc = Svc.con(h,16);
        if alpha < a_max & alpha > a_min
          k1 = -K*T1/T2/Tm*Km;
          k2 = K*T1/T2/Tm-K/T2;
          k3 = -2*DAE.V(j)*DAE.V(j)*(1-cos(2*alpha))/xl/pi;
          DAE.Fx(Svc.alpha(h),k) = K/T2;
          DAE.Fy(k,Svc.alpha(h)) = -k3*dVs_dQ(idx(i));
        end
      end

    end

  end

 case 5 % non-windup limiters

  global Settings
  idx = find((Vs >= Vs_max | Vs <= Vs_min) & DAE.f(Cluster.Vs) == 0);
  if ~isempty(idx)
    k = Cluster.Vs(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = eye(length(idx));
    else
      DAE.Ac(k,k) = speye(length(idx));
    end
  end

end


% -------------------------------------------------------------------
% function for creating warning messages
function clusterwarn(idx, msg)
fm_disp(strcat('Warning: CLUSTER #',int2str(idx),msg))