fm_cluster.m

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

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-2006 Federico Milano

global Cluster DAE Bus Exc Svc Syn Cac

%Data Format Cluster.con:
%          col #1: Central Area Control number
%          col #2: AVR or SVC number
%          col #3: type (1) AVR, (2) SVC
%          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.)
  Qg = zeros(Cluster.n,1);
  if Cluster.exc 
    Qg(Cluster.exc) = -Syn.Qg(Syn.cluster);
  end
  if Cluster.svc
    Qg(Cluster.svc) = getq(Svc,Svc.cluster);
  end
  Cluster.con(:,7) = Qg;
  Cluster.dVsdQ = (Xtg+Xeq)./T;
  Cluster.u = ones(Cluster.n,1);

  % 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

  Qg = zeros(Cluster.n,1);
  if Cluster.exc 
    Qg(Cluster.exc) = -Syn.Qg(Syn.cluster);
  end
  if Cluster.svc
    Qg(Cluster.svc) = getq(Svc,Svc.cluster);
  end
  q = Cac.q(Cluster.cac);
  DAE.f(Cluster.Vs) = ((Qgr.*q-Qg).*(Xtg+Xeq))./T;
  % non-windup 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);
  Cluster.u = Vs <= 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(:,3);
  idx = find(Cluster.u);
  
  for i = 1:length(idx)

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

      j = Cluster.bus(idx(i));
      k = Cluster.Vs(idx(i));
      h = Exc.cluster(idx(i));
      m = Syn.cluster(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

      alpha = Cluster.dVsdQ(idx(i));

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

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

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

    end
    
  end

 case 5 % non-windup limiters

  fm_windup(Cluster.Vs,Cluster.con(:,8),Cluster.con(:,9))
  
end


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