fm_fl.m

电力系统的psat

function fm_fl(flag)
% FM_FL define General Exponential Voltage Frequency Dependent Loads
%
%FM_FL(FLAG)
%       FLAG = 0 initialization
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations (frequency filter)
%       FLAG = 4 state Jacobians (frequency filter)
%
%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 Fl PQ DAE Bus Settings

V = DAE.V(Fl.bus);

switch flag
 case 0
  for i = 1:Fl.n
    idx = find(PQ.bus == Fl.bus(i));
    if isempty(idx)
      fm_disp(['No PQ load found for initializing frequency ', ...
	       'dependent load at bus #',int2str(Fl.bus(i))])
    else
      PQ.con(idx,4) = PQ.con(idx,4) - ...
	  PQ.P0(idx)*Fl.con(i,2)/100;
      PQ.con(idx,5) = PQ.con(idx,5) - ...
	  PQ.Q0(idx)*Fl.con(i,5)/100;
      Fl.con(i,2)  = Fl.con(i,2)*PQ.P0(idx)/(V(i)^Fl.con(i,3))/100;
      Fl.con(i,5)  = Fl.con(i,5)*PQ.Q0(idx)/(V(i)^Fl.con(i,6))/100;
      if PQ.con(idx,4) == 0 & PQ.con(idx,5) == 0,
	PQ.con(idx,:) = [];
	PQ.bus(idx) = [];
	PQ.n = PQ.n - 1;
	PQ.P0(idx) = [];
	PQ.Q0(idx) = [];
      end
    end
  end
  DAE.x(Fl.x) = 0;
  Fl.a0 = DAE.a(Fl.bus);
  Fl.Dw = zeros(Fl.n,1);
  fm_disp('Initialization of frequency dependent loads completed.')
 case 1
  DAE.gp = DAE.gp + sparse(Fl.bus,1,Fl.con(:,2).*(V.^Fl.con(:,3)).* ...
			   (1+Fl.Dw).^Fl.con(:,4),Bus.n,1);
  DAE.gq = DAE.gq + sparse(Fl.bus,1,Fl.con(:,5).*(V.^Fl.con(:,6)).* ...
			   (1+Fl.Dw).^Fl.con(:,7),Bus.n,1);
 case 2
  P0 = Fl.con(:,2);
  Q0 = Fl.con(:,5);
  ap = Fl.con(:,3);
  aq = Fl.con(:,6);
  bp = Fl.con(:,4);
  bq = Fl.con(:,7);
  dw_da = 1./Fl.con(:,8)/Settings.rad;
  DAE.J11 = DAE.J11 + sparse(Fl.bus,Fl.bus,P0.*V.^ap.*bp.*dw_da.* ...
			     (1+Fl.Dw).^(bp-1),Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(Fl.bus,Fl.bus,P0.*ap.*V.^(ap-1).* ...
			     (1+Fl.Dw).^bp,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Fl.bus,Fl.bus,Q0.*aq.*V.^(aq-1).* ...
			     (1+Fl.Dw).^bq,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Fl.bus,Fl.bus,Q0.*V.^aq.*bq.*dw_da.* ...
			     (1+Fl.Dw).^(bq-1),Bus.n,Bus.n);
 case 3
  a = DAE.a(Fl.bus);
  k = (a-Fl.a0)./Fl.con(:,8)/Settings.rad;
  Fl.Dw = DAE.x(Fl.x) + k;
  DAE.f(Fl.x) = -Fl.Dw./Fl.con(:,8);
 case 4
  P0 = Fl.con(:,2);
  Q0 = Fl.con(:,5);
  ap = Fl.con(:,3);
  aq = Fl.con(:,6);
  bp = Fl.con(:,4);
  bq = Fl.con(:,7);
  DAE.Fx = DAE.Fx + sparse(Fl.x,Fl.x,-1./Fl.con(:,8),DAE.n,DAE.n);
  DAE.Fy = DAE.Fy + sparse(Fl.x,Fl.bus,-1./Fl.con(:,8)./ ...
			   Fl.con(:,8)/Settings.rad,DAE.n,2*Bus.n);
  DAE.Gx = DAE.Gx + sparse(Fl.bus,Fl.x,P0.*V.^ap.*bp.* ...
			   (1+Fl.Dw).^(bp-1),2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Fl.bus+Bus.n,Fl.x,Q0.*V.^aq.*bq.* ...
			   (1+Fl.Dw).^(bq-1),2*Bus.n,DAE.n);
end