fm_statcom.m

电力系统的psat

function  fm_statcom(flag)
% FM_STATCOM defines Statcom components
%
% FM_STATCOM(FLAG)
%       FLAG = 0 initialization
%       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
%Update:    07-Apr-2003
%Version:   1.0.1
%
%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 Statcom Bus DAE Settings

type = Statcom.con(:,2);
V = DAE.V(Statcom.bus);

Vrefac = Statcom.con(:,6);
Vrefdc = Statcom.con(:,7);
I_max = Statcom.con(:,8);
I_min = Statcom.con(:,9);
Rac = Statcom.con(:,10);
Xac = Statcom.con(:,11);
Rdc = Statcom.con(:,12);
Cdc = Statcom.con(:,13);
Ka = Statcom.con(:,14);
Kd = Statcom.con(:,15);
T1 = Statcom.con(:,16);
T2 = Statcom.con(:,17);
Kp = Statcom.con(:,18);
Ki = Statcom.con(:,19);
Kmac = Statcom.con(:,20);
Tmac = Statcom.con(:,21);
Kmdc = Statcom.con(:,22);
Tmdc = Statcom.con(:,23);

Vdc = DAE.x(Statcom.Vdc);
alpha = DAE.x(Statcom.alpha);
m = DAE.x(Statcom.m);
Vmdc = DAE.x(Statcom.Vmdc);
Vmac = DAE.x(Statcom.Vmac);
if Settings.init
  a_max = Statcom.dat(:,1);
  a_min = Statcom.dat(:,2);
else
  a_max =  ones(Statcom.n,1);
  a_min = -ones(Statcom.n,1);
end
m_max = Statcom.dat(:,3);
m_min = Statcom.dat(:,4);
Z2 = Rac.*Rac+Xac.*Xac;
G = Rac./Z2;
B = Xac./Z2;
Ca = cos(alpha);
Sa = sin(alpha);
k0 = sqrt(3/8);
k1 = k0*m.*Vdc;

switch flag
case 0 % initialization

 % m_max
 Statcom.dat(:,3) = (Vrefac+0.05*I_max)./Vrefdc/k0;
 m_max = Statcom.dat(:,3);

 % m_min
 Statcom.dat(:,4) = (Vrefac-0.05*I_min)./Vrefdc/k0;
 m_min = Statcom.dat(:,4);
 idx = find(type ~= 1);
 if idx, Statcom.dat(idx,4) = 0.5*0.9/k0; m_min(idx) = Statcom.dat(idx,4); end

 % a_max
 CC = Vrefac.*Vrefac.*G-Vrefdc.*Vrefdc./Rdc-Rac.*I_max.*I_max;
 BB = -k0.*m_max.*Vrefdc.*Vrefac.*G;
 AA = -k0.*m_max.*Vrefdc.*Vrefac.*B;
 d = AA.*AA+BB.*BB;
 b = BB.*CC./d;
 c = (CC.*CC-AA.*AA)./d;
 Statcom.dat(:,1) = max(acos(b+sqrt(b.*b-c)),0.1);

 % a_min
 CC = Vrefac.*Vrefac.*G-Vrefdc.*Vrefdc./Rdc-Rac.*I_min.*I_min;
 BB = -k0.*m_min.*Vrefdc.*Vrefac.*G;
 AA = -k0.*m_min.*Vrefdc.*Vrefac.*B;
 d = AA.*AA+BB.*BB;
 b = BB.*CC./d;
 c = (CC.*CC-AA.*AA)./d;
 Statcom.dat(:,2) = min(-acos(b+sqrt(b.*b-c)),-0.1);

 case 1 % algebraic equations

  Statcom.Id =  V.*G-k1.*G.*Ca+k1.*B.*Sa;
  Statcom.Iq = -V.*B+k1.*B.*Ca+k1.*G.*Sa;
  DAE.gp = DAE.gp + sparse(Statcom.bus,1,V.*Statcom.Id,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Statcom.bus,1,V.*Statcom.Iq,Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J12 = DAE.J12 + sparse(Statcom.bus,Statcom.bus, V.*G+Statcom.Id,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Statcom.bus,Statcom.bus,-V.*B+Statcom.Iq,Bus.n,Bus.n);

 case 3 % differential equations

  Id = Statcom.Id;
  Iq = Statcom.Iq;

  DAE.f(Statcom.Vdc) = V.*Id./Cdc./Vdc-Vdc./Rdc./Cdc - ...
      Rac./Cdc.*(Id.*Id+Iq.*Iq)./Vdc;
  DAE.f(Statcom.m) = (-Kd.*m+Ka.*(T1./Tmac-1).*Vmac + ...
                      Ka.*Vrefac-Ka.*T1.*Kmac./Tmac.*V)./T2;
  DAE.f(Statcom.Vmdc) = (Kmdc.*Vdc-Vmdc)./Tmdc;
  DAE.f(Statcom.Vmac) = (Kmac.*V-Vmac)./Tmac;

  jdx = find(m >= m_max & DAE.f(Statcom.m) > 0);
  if jdx, DAE.f(Statcom.m(jdx)) = 0; end
  jdx = find(m <= m_min & DAE.f(Statcom.m) < 0);
  if jdx, DAE.f(Statcom.m(jdx)) = 0; end
  DAE.x(Statcom.m) = max(DAE.x(Statcom.m),m_min);
  DAE.x(Statcom.m) = min(DAE.x(Statcom.m),m_max);

  i = find(type == 1);
  if i,
    DAE.f(Statcom.alpha(i)) = (Kp(i)./Tmdc(i)-Ki(i)).*Vmdc(i)+ ...
        Ki(i).*Vrefdc(i)-Kp(i).*Kmdc(i).*Vdc(i)./Tmdc(i);
  end
  i = find(type == 2);
  if i,
    DAE.f(Statcom.alpha(i)) = -Kp(i).*k0.*DAE.f(Statcom.m(i)) + ...
        Ki(i).*(0.9 - k0.*DAE.x(Statcom.m(i)));
  end
  i = find(type == 3);
  if i,
    DAE.f(Statcom.alpha(i)) = (Kp(i).*(0.9-k0*DAE.x(Statcom.m(i))) ...
                               - alpha(i))./Ki(i);
  end

  jdx = find(alpha >= a_max & DAE.f(Statcom.alpha) > 0);
  if jdx,
    DAE.f(Statcom.alpha(jdx)) = 0;
  end
  jdx = find(alpha <= a_min & DAE.f(Statcom.alpha) < 0);
  if jdx,
    DAE.f(Statcom.alpha(jdx)) = 0;
  end
  DAE.x(Statcom.alpha) = max(DAE.x(Statcom.alpha),a_min);
  DAE.x(Statcom.alpha) = min(DAE.x(Statcom.alpha),a_max);

 case 4 % state Jacobians

  Id = Statcom.Id;
  Iq = Statcom.Iq;
  I2 = Id.*Id+Iq.*Iq;

  P1 =  k0*(B.*Sa-G.*Ca);
  P2 = P1.*m;
  P3 = P1.*Vdc;
  Q4 = -P1.*m.*Vdc;
  Q1 =  k0*(B.*Ca+G.*Sa);
  Q2 = Q1.*m;
  Q3 = Q1.*Vdc;
  P4 = Q1.*m.*Vdc;

  M1 = Ka.*(T1./Tmac-1)./T2;
  M2 = -Ka.*T1.*Kmac./Tmac./T2;

  DAE.Fx = DAE.Fx + sparse(Statcom.Vdc,Statcom.Vdc, ...
           (V.*(P2-Id./Vdc)./Vdc-1./Rdc+ ...
            Rac.*(I2./Vdc-2*(Id.*P2+Iq.*Q2))./Vdc)./Cdc,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Statcom.m,Statcom.m,-Kd./T2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Statcom.Vmdc,Statcom.Vdc,Kmdc./Tmdc,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Statcom.Vmdc,Statcom.Vmdc,-1./Tmdc,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Statcom.Vmac,Statcom.Vmac,-1./Tmac,DAE.n,DAE.n);

  DAE.Fy = DAE.Fy + sparse(Statcom.Vdc,Statcom.bus+Bus.n, ...
           ((V.*G+Id)./Vdc-2*Rac.*(Id.*G-Iq.*B)./Vdc)./Cdc,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Statcom.Vmac,Statcom.bus+Bus.n,Kmac./Tmac,DAE.n,2*Bus.n);

  DAE.Gx = DAE.Gx + sparse(Statcom.bus,Statcom.Vdc,V.*P2,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Statcom.bus+Bus.n,Statcom.Vdc,V.*Q2,2*Bus.n,DAE.n);

  k = find(m < m_max & m > m_min);
  if k,
    DAE.Fx = DAE.Fx + sparse(Statcom.Vdc(k),Statcom.m(k), ...
             (V(k).*P3(k)./Vdc(k)-2*Rac(k).*(Id(k).*P3(k)+ ...
             Iq(k).*Q3(k))./Vdc(k))./Cdc(k),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Statcom.m(k),Statcom.Vmac(k),M1(k),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Statcom.m(k),Statcom.bus+Bus.n,M2(k),DAE.n,2*Bus.n);
    DAE.Gx = DAE.Gx + sparse(Statcom.bus(k),Statcom.m(k),V(k).*P3(k),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(Statcom.bus(k)+Bus.n,Statcom.m(k),V(k).*Q3(k),2*Bus.n,DAE.n);
  end

  k = find(alpha < a_max & alpha > a_min);
  if k,
    DAE.Fx = DAE.Fx + sparse(Statcom.Vdc(k),Statcom.alpha(k), ...
             (V(k).*P4(k)./Vdc(k)-2*Rac(k).*(Id(k).*P4(k)+ ...
             Iq(k).*Q4(k))./Vdc(k))./Cdc(k),DAE.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(Statcom.bus(k),Statcom.alpha(k),V(k).*P4(k),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(Statcom.bus(k)+Bus.n,Statcom.alpha(k),V(k).*Q4(k),2*Bus.n,DAE.n);
    j = find(type(k) == 1);
    if j,
      h = k(j);
      DAE.Fx = DAE.Fx + sparse(Statcom.alpha(h),Statcom.Vdc(h),-Kp(h).*Kmdc(h)./Tmdc(h),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Statcom.alpha(h),Statcom.Vmdc(h),Kp(h)./Tmdc(h)-Ki(h),DAE.n,DAE.n);
    end
    j = find(type(k) == 2);
    if j,
      h = k(j);
      n = find(m(h) < m_max(h) & m(h) > m_min(h));
      if n,
        DAE.Fx = DAE.Fx + sparse(Statcom.alpha(h(n)),Statcom.m(h(n)), ...
                 -Kp(h(n)).*k0.*-Kd(h(n))./T2(h(n))-k0*Ki(h(n)), ...
                                 DAE.n,DAE.n);
      end
      DAE.Fx = DAE.Fx + sparse(Statcom.alpha(h),Statcom.Vmac(h),-k0*Kp(h).*M1(h),DAE.n,DAE.n);
      DAE.Fy = DAE.Fx + sparse(Statcom.alpha(h),Statcom.bus(h)+Bus.n,-k0*Kp(h).*M2(h),DAE.n,2*Bus.n);
    end
    j = find(type(k) == 3);
    if j,
      h = k(j);
      n = find(m(h) < m_max(h) & m(h) > m_min(h));
      if n,
        DAE.Fx = DAE.Fx + sparse(Statcom.alpha(h(n)), ...
                                 Statcom.m(h(n)),-k0*Kp(h(n))./Ki(h(n)),DAE.n,DAE.n);
      end
      DAE.Fx = DAE.Fx + sparse(Statcom.alpha(h),Statcom.alpha(h),-1./Ki(h),DAE.n,DAE.n);
    end
  end
  k = find(alpha >= a_max | alpha <= a_min);
  if k,
    DAE.Fx = DAE.Fx + sparse(Statcom.alpha(k),Statcom.alpha(k),-1, ...
                             DAE.n,DAE.n);
  end

 case 5 % non-windup limiters

  jdx = find((alpha >= a_max | alpha <= a_min) &  DAE.f(Statcom.alpha) == 0);
  if jdx
    k = Statcom.alpha(jdx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = -eye(length(jdx));
    else
      DAE.Ac(k,k) = -speye(length(jdx));
    end
  end
  jdx = find((m >= m_max | m <= m_min) &  DAE.f(Statcom.m) == 0);
  if jdx
    k = Statcom.m(jdx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = -eye(length(jdx));
    else
      DAE.Ac(k,k) = -speye(length(jdx));
    end
  end

end