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-2006 Federico Milano

global Fl PQ DAE Bus Settings Varname

V = DAE.V(Fl.bus);

switch flag
 case 0
  for i = 1:Fl.n
    idx = findbus(PQ,Fl.bus(i));
    if isempty(idx)
      fm_disp(['No PQ load found for initializing frequency ', ...
	       'dependent load at bus ',Varname.bus{Fl.bus(i)}])
    else
      P = PQ.P0(idx)*Fl.con(i,2)/100;
      Q = PQ.Q0(idx)*Fl.con(i,5)/100;
      PQ = pqsub(PQ,idx,P,Q);
      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;
      PQ = remove(PQ,idx,'zero');
    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