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:    Hugo M. Ayres (revised by Federico Milano)
%Date:      18-Feb-2006
%Version:   1.0.0
%
%E-mail:    [e-mail]
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2006 Hugo M. Ayres & Federico Milano

global Statcom Syn Pod Bus DAE PV

V = DAE.V(Statcom.bus);
ist = DAE.x(Statcom.ist);
Kr = Statcom.con(:,5);
Tr = Statcom.con(:,6);
ist_max = Statcom.con(:,7);
ist_min = Statcom.con(:,8);

% Algebraic Equations:

switch flag
 case 0

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

  DAE.x(Statcom.ist) = Bus.Qg(Statcom.bus)./V;
  idx = find(DAE.x(Statcom.ist) > ist_max);
  if idx, statcomwarn(idx,' Ish is over its max limit.'), end
  idx = find(DAE.x(Statcom.ist) < ist_min);
  if idx, statcomwarn(idx,' Ish is under its min limit.'), end
  DAE.x(Statcom.ist) = max(DAE.x(Statcom.ist),ist_min);
  DAE.x(Statcom.ist) = min(DAE.x(Statcom.ist),ist_max);

  % reference voltages
  Statcom.Vref = DAE.x(Statcom.ist)./Kr + V;

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

 case 1 % algebraic equations

  DAE.gq = DAE.gq - sparse(Statcom.bus,1,ist.*V,Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J22 = DAE.J22 - sparse(Statcom.bus,Statcom.bus,ist,Bus.n,Bus.n);

 case 3 % differential equations

  % Updating reference voltages
  Vref = Statcom.Vref;
  if Statcom.pod
    Vref(Statcom.pod) = Vref(Statcom.pod) + DAE.x(Pod.Vs(Pod.statcom));
  end
    
  u = ~(ist >= ist_max & DAE.f(Statcom.ist) > 0) & ...
      ~(ist <= ist_min & DAE.f(Statcom.ist) < 0);
  DAE.f(Statcom.ist) = u.*(Kr.*(Vref-V)-ist)./Tr;  
  DAE.x(Statcom.ist) = min(ist_max,DAE.x(Statcom.ist));
  DAE.x(Statcom.ist) = max(ist_min,DAE.x(Statcom.ist));

 case 4 % Jacobians of state variables

  DAE.Fx = DAE.Fx + sparse(Statcom.ist,Statcom.ist,-1./Tr,DAE.n,DAE.n);     

  u = (ist <= ist_max & ist >= ist_min);
  
  DAE.Gx = DAE.Gx - sparse(Bus.n+Statcom.bus,Statcom.ist,u.*V,2*Bus.n,DAE.n);
  DAE.Fy = DAE.Fy - sparse(Statcom.ist,Bus.n+Statcom.bus,u.*Kr./Tr, ...
                           DAE.n,2*Bus.n);
  if Statcom.pod
    DAE.Fx(Statcom.ist(Statcom.pod),Pod.Vs(Pod.statcom)) = ...
        Pod.u(Pod.statcom).*u(Statcom.pod).*Kr(Statcom.pod)./Tr(Statcom.pod);
  end
  
 case 5 % non-windup limiters

  a = find((ist >= ist_max | ist <= ist_min) & DAE.f(Statcom.ist) == 0);
  if ~isempty(a)
    DAE.tn(Statcom.ist(a)) = 0;
    DAE.Ac(Statcom.ist(a),:) = 0;
    DAE.Ac(:,Statcom.ist(a)) = 0;
    DAE.Ac(Statcom.ist(a),Statcom.ist(a)) = speye(length(a));
  end

end  

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