fm_svc.m

电力系统的psat

function  fm_svc(flag)
% FM_SVC define Static Var Systems
%
% FM_SVC(FLAG)
%       FLAG = 0 initializations
%       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 Svc Bus DAE PV Cluster

type = Svc.con(:,5);
ty1 = find(type == 1);
ty2 = find(type == 2);
V = DAE.V(Svc.bus);

if ty1
  bcv = DAE.x(Svc.bcv);
  Tr = Svc.con(ty1,6);
  Kr = Svc.con(ty1,7);
  %Vref1 = Svc.con(ty1,8);
  bcv_max = Svc.con(ty1,9);
  bcv_min = Svc.con(ty1,10);
  Svc.Be(ty1) = DAE.x(Svc.bcv);
end
if ty2
  alpha = DAE.x(Svc.alpha);
  vm = DAE.x(Svc.vm);
  a_max = Svc.con(ty2,9);
  a_min = Svc.con(ty2,10);
  %Vref2 = Svc.con(ty2,8);
  T2 = Svc.con(ty2,6);
  K = Svc.con(ty2,7);
  Kd = Svc.con(ty2,11);
  T1 = Svc.con(ty2,12);
  Km = Svc.con(ty2,13);
  Tm = Svc.con(ty2,14);
  xl = Svc.con(ty2,15);
  xc = Svc.con(ty2,16);
  Svc.Be(ty2) = (2*DAE.x(Svc.alpha) - sin(2*DAE.x(Svc.alpha)) ...
                 - pi*(2-xl./xc))./(pi*xl);
end

% Algebraic Equations:

switch flag
 case 0

  % eliminate PV components used for initializing SVC's
  global Syn
  for i = 1:Svc.n
    idxg = find(Syn.bus == Svc.bus(i));
    if ~isempty(idxg)
      svcwarn(i,[' SVC cannot be connected at the same bus as ' ...
                 'synchronous machines.'])
      continue
    end
    idx = find(PV.bus == Svc.bus(i));
    if isempty(idx)
      svcwarn(i,' no PV generator found at the bus.')
    else
      PV.con(idx,:) = [];
      PV.bus(idx) = [];
      PV.n = PV.n-1;
    end
  end

  if ty1
    DAE.x(Svc.bcv) = Bus.Qg(Svc.bus(ty1))./V(ty1)./V(ty1);
    Svc.con(ty1,8) = DAE.x(Svc.bcv)./Kr + V(ty1);
    idx = find(DAE.x(Svc.bcv) > bcv_max);
    if idx, svcwarn(ty1(idx),' b_svc is over its max limit.'), end
    idx = find(DAE.x(Svc.bcv) < bcv_min);
    if idx, svcwarn(ty1(idx),' b_svc is under its min limit.'), end
    DAE.x(Svc.bcv) = max(DAE.x(Svc.bcv),bcv_min);
    DAE.x(Svc.bcv) = min(DAE.x(Svc.bcv),bcv_max);
    Svc.Be(ty1) = DAE.x(Svc.bcv);
  end

  if ty2
    DAE.x(Svc.vm) = V(ty2)./Km;
    b = pi*(2-xl./xc)+pi*xl.*Bus.Qg(Svc.bus(ty2))./V(ty2)./V(ty2);
    % numeric solution for alpha
    for i = 1:length(ty2)
      err = 1;
      % first guess is the expansion of a 3rd order Taylor series
      a = sign(b(i))*((6*abs(b(i)))^(1/3))/2;
      iter = 0;
      while abs(err) > 1e-8
	if iter > 20,
	  svcwarn(ty2(i),' convergence not reached in computing alpha')
	  break,
	end
	ga = 2*a - sin(2*a) - b(i);
	ja = 2*(1-cos(2*a));
	err = -ga/ja;
	a = a + err;
	iter = iter + 1;
      end
      DAE.x(Svc.alpha(i)) = a;
    end
    Svc.con(ty2,8) = DAE.x(Svc.vm) + Kd./K.*DAE.x(Svc.alpha);
    idx = find(DAE.x(Svc.alpha) > a_max);
    if idx, svcwarn(ty2(idx),' alpha is over its max limit.'), end
    idx = find(DAE.x(Svc.alpha) < a_min);
    if idx, svcwarn(ty2(idx),' alpha is under its min limit.'), end
    DAE.x(Svc.alpha) = max(DAE.x(Svc.alpha),a_min);
    DAE.x(Svc.alpha) = min(DAE.x(Svc.alpha),a_max);
    Svc.Be(ty2) = (2*DAE.x(Svc.alpha) - sin(2*DAE.x(Svc.alpha)) - ...
                   pi*(2-xl./xc))./(pi*xl);
  end

  % reference voltages
  Svc.Vref = Svc.con(:,8);

  fm_disp('Initialization of Static Var Compensators completed.')

 case 1 % algebraic equations

  DAE.gq = DAE.gq + sparse(Svc.bus,1,-Svc.Be.*V.*V,Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J22 = DAE.J22 + sparse(Svc.bus,Svc.bus,-2*Svc.Be.*V,Bus.n,Bus.n);

 case 3 % differential equations

  %Vref = Svc.con(:,8);
  Vref = Svc.Vref;
  if Cluster.n,
    type = find(Cluster.con(:,2) == 2);
    if ~isempty(type)
      Vref(Cluster.svc(type)) = Vref(Cluster.svc(type)) + ...
          DAE.x(Cluster.Vs(type));
    end
  end

  if ~isempty(ty1)
    DAE.f(Svc.bcv) = (Kr.*(Vref(ty1)-V(ty1))-bcv)./Tr;
    a = find(bcv >= bcv_max & DAE.f(Svc.bcv) > 0);
    if ~isempty(a), DAE.f(Svc.bcv(a)) = 0; end
    a = find(bcv <= bcv_min & DAE.f(Svc.bcv) < 0);
    if ~isempty(a), DAE.f(Svc.bcv(a)) = 0; end
    DAE.x(Svc.bcv) = min(bcv_max,DAE.x(Svc.bcv));
    DAE.x(Svc.bcv) = max(bcv_min,DAE.x(Svc.bcv));
    %Svc.Be(ty1) = DAE.x(Svc.bcv);
  end
  if ~isempty(ty2)
    DAE.f(Svc.vm) = (Km.*V(ty2)-vm)./Tm;
    DAE.f(Svc.alpha) = -Kd./T2.*alpha + K.*T1./T2./Tm.*(vm-Km.*V(ty2)) ...
        + K./T2.*(Vref(ty2)-vm);
    a = find(alpha >= a_max & DAE.f(Svc.alpha) > 0);
    if ~isempty(a), DAE.f(Svc.alpha(a)) = 0; end
    a = find(alpha <= a_min & DAE.f(Svc.alpha) < 0);
    if ~isempty(a), DAE.f(Svc.alpha(a)) = 0; end
    DAE.x(Svc.alpha) = min(a_max,DAE.x(Svc.alpha));
    DAE.x(Svc.alpha) = max(a_min,DAE.x(Svc.alpha));
    %Svc.Be(ty2) = (2*DAE.x(Svc.alpha) - sin(2*DAE.x(Svc.alpha)) - ...
    %               pi*(2-xl./xc))./(pi*xl);
  end

 case 4 % Jacobians of state variables

  if ~isempty(ty1)
    DAE.Fx = DAE.Fx + sparse(Svc.bcv,Svc.bcv,-1./Tr,DAE.n,DAE.n);
    a = find((bcv <= bcv_max & bcv >= bcv_min) & DAE.f(Svc.bcv) ~= 0);
    if ~isempty(a)
      DAE.Gx = DAE.Gx + sparse(Bus.n+Svc.bus(ty1(a)),Svc.bcv(a), ...
                               -V(ty1(a)).^2,2*Bus.n,DAE.n);
      DAE.Fy = DAE.Fy + sparse(Svc.bcv(a),Bus.n+Svc.bus(ty1(a)), ...
                               -Kr(a)./Tr(a),DAE.n,2*Bus.n);
    end
  end
  if ~isempty(ty2)
    DAE.Fx = DAE.Fx + sparse(Svc.vm,Svc.vm,-1./Tm,DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Svc.vm,Bus.n+Svc.bus(ty2),Km./Tm,DAE.n,2*Bus.n);
    DAE.Fx = DAE.Fx + sparse(Svc.alpha,Svc.alpha,-Kd./T2,DAE.n,DAE.n);
    a = find(alpha < a_max & alpha > a_min);
    if ~isempty(a)
      k1 = -K(a).*T1(a)./T2(a)./Tm(a).*Km(a);
      k2 = K(a).*T1(a)./T2(a)./Tm(a) - K(a)./T2(a);
      k3 = -2*V(ty2(a)).*V(ty2(a)).*(1-cos(2*alpha(a)))./xl(a)/pi;
      DAE.Fy = DAE.Fy + sparse(Svc.alpha(a),Svc.bus(ty2(a))+Bus.n,k1,DAE.n,2*Bus.n);
      DAE.Fx = DAE.Fx + sparse(Svc.alpha(a),Svc.vm(a),k2,DAE.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Bus.n+Svc.bus(ty2(a)),Svc.alpha(a),k3,2*Bus.n,DAE.n);
    end
  end

 case 5 % non-windup limiters

  global Settings

  if ~isempty(ty1)
    a = find((bcv >= bcv_max | bcv <= bcv_min) & DAE.f(Svc.bcv) == 0);
    if ~isempty(a)
      DAE.tn(Svc.bcv(a)) = 0;
      DAE.Ac(Svc.bcv(a),:) = 0;
      DAE.Ac(:,Svc.bcv(a)) = 0;
      if Settings.octave
        DAE.Ac(Svc.bcv(a),Svc.bcv(a)) = eye(length(a));
      else
        DAE.Ac(Svc.bcv(a),Svc.bcv(a)) = speye(length(a));
      end
    end
  end
  if ~isempty(ty2)
    a = find((alpha >= a_max | alpha <= a_min) & DAE.f(Svc.alpha) == 0);
    if ~isempty(a)
      DAE.tn(Svc.alpha(a)) = 0;
      DAE.Ac(Svc.alpha(a),:) = 0;
      DAE.Ac(:,Svc.alpha(a)) = 0;
      if Settings.octave
        DAE.Ac(Svc.alpha(a),Svc.alpha(a)) = eye(length(a));
      else
        DAE.Ac(Svc.alpha(a),Svc.alpha(a)) = speye(length(a));
      end
    end
  end

end

% function for creating warning messages
function svcwarn(idx, msg)
global Svc Varname
fm_disp(strcat('Warning: SVC #',int2str(idx), ...
	       ' at bus <',Varname.bus(Svc.bus(idx)), ...
	       '>: ',msg))