fm_oxl.m

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

function  fm_oxl(flag)
% FM_OXL define Over Excitation Limiters
%
% FM_OXL(FLAG)
%       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-2006 Federico Milano

global Bus Exc Syn DAE Oxl Varname

% OXL state variable
V_oxl = DAE.x(Oxl.v);

% data
xd = Oxl.con(:,4);
xq = Oxl.con(:,5);
idx = find(Oxl.con(:,3));
xd(idx) = Syn.con(Oxl.syn(idx),8);
xq(idx) = Syn.con(Oxl.syn(idx),13);
I_lim = Oxl.con(:,6);

% voltage at generator bus
Vg = DAE.V(Oxl.bus);
% generator powers
Pg = -Syn.Pg(Oxl.syn);
Qg = -Syn.Qg(Oxl.syn);

% some useful quantities
t1 = xq.*Qg;
t2 = 1./Vg;
t4 = Vg+t1.*t2;
t5 = t4.*t4;
t6 = Pg.*Pg;
t7 = t5+t6;
t8 = sqrt(t7);
t11 = xd./xq-1.0;
t12 = xq.*t4;
t13 = Qg.*t2;
t15 = xq.*xq;
t16 = t15.*t6;
t17 = Vg.*Vg;
t18 = 1./t17;
t21 = t11.*(t12.*t13+t16.*t18);
t22 = 1./t8;

switch flag
 case 0

  % initialization and check for field current limits
  DAE.x(Oxl.v) = zeros(Oxl.n,1);
  Oxl.If = t8+t21.*t22;
  idx = find(Oxl.If > I_lim);
  for i = 1:length(idx)
    Oxl.con(i,6) = 1.2*Oxl.If(i);
    fm_disp(['Field current of Synchronous Machine #',int2str(Oxl.syn(idx(i))), ...
             ' ( -> ',Varname.bus{Oxl.bus(idx(i))},') is over its thermal limit.'])
    fm_disp(['Field current limit of OXL #',int2str(i),' has been changed to "',num2str(Oxl.con(i,6)),'"'])
  end
  idx = find(Oxl.con(:,2) == 0);
  for i = 1:length(idx)
    Oxl.con(i,2) = 10;
    fm_disp(['OXL #',int2str(i),' has null T0. The default value T0 = 10 s will be used.'])
  end
  fm_disp('Initialization of Over Excitation Limiters completed.')

 case 3

  % time derivative and check for limits
  Oxl.If = t8+t21.*t22;
  DAE.f(Oxl.v) = (Oxl.If - I_lim)./Oxl.con(:,2);
  idx1 = find(V_oxl <= 0);
  if ~isempty(idx1)
    DAE.x(Oxl.v(idx1)) = 0;
    idx2 = find(DAE.f(Oxl.v(idx1)) < 0);
    DAE.f(Oxl.v(idx1(idx2))) = 0;
  end

 case 4

  K0 = 1./Oxl.con(:,2);
  DAE.Fx = DAE.Fx + sparse(Oxl.v,Oxl.v,-1e-6,DAE.n,DAE.n);
  idx = find(V_oxl > 0);
  if ~isempty(idx)

    % other useful quantities
    t30 = t7(idx).*t7(idx);
    t32 = t8(idx)./t30;
    t36 = t22(idx).*t4(idx);
    t41 = t12(idx).*t2(idx);
    t49 = 1.0-t1(idx).*t18(idx);

    % partial derivatives of the field current with respect to Vg, Pg and Qg
    dI_dV = -K0(idx).*(t36.*t49+t11(idx).*(xq(idx).*t49.*t13(idx)-t12(idx).*Qg(idx).*t18(idx)- ...
                                           2.0*t16(idx)./t17(idx)./Vg(idx)).*t22(idx)-t21(idx).*t32.*t4(idx).*t49);
    dI_dP = -K0(idx).*(t22.*Pg(idx)+2.0*t11(idx).*t15(idx).*Pg(idx).*t18(idx).*t22(idx)- ...
                       t21(idx).*t32.*Pg(idx));
    dI_dQ = -K0(idx).*(t36.*xq(idx).*t2(idx)+t11(idx).*(t15(idx).*t18(idx).*Qg(idx)+ ...
                                                      t41).*t22(idx)-t21(idx).*t32.*t41);

    % Jacobian matrices DAE.Fx & DAE.Fy
    for i = 1:length(idx)

      h = Oxl.exc(idx(i));
      j = Oxl.bus(idx(i));
      m = Oxl.syn(idx(i));
      k = Oxl.v(idx(i));

      % update Jacobians of Synchronous Machines & AVRs
      type_exc = Exc.con(Oxl.exc(idx(i)),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

      % OXL Jacobian of algebraic variables (ineherited from synchronous machines)
      DAE.Fy(k,j) = dI_dP(i)*Syn.J11(m) + dI_dQ(i)*Syn.J21(m);
      DAE.Fy(k,j+Bus.n) = dI_dP(i)*Syn.J12(m) + dI_dQ(i)*Syn.J22(m) + dI_dV(i);

      % OXL Jacobian of state variables (inherited from synchronous machines)
      DAE.Fx(k,Syn.delta(m)) = dI_dP(i)*Syn.Gp(m,1) + dI_dQ(i)*Syn.Gq(m,1);

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

  end

 case 5

  fm_windup(Oxl.v,1e10,0);

end