fm_exload.m

电力系统的psat

function  fm_exload(flag)
% FM_EXLOAD defines exponential recovery load
%
% FM_EXLOAD(FLAG)
%       FLAG = 0 -> initialization
%       FLAG = 1 -> algebraic equations
%       FLAG = 2 -> algebraic Jacobians
%       FLAG = 3 -> differential equations
%       FLAG = 4 -> state Jacobians
%
%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 Exload DAE Bus PQ

%Data Structure: Exload.dat:
%          col #1: P0
%          col #2: Q0
%          col #3: V0

xp = DAE.x(Exload.xp);
xq = DAE.x(Exload.xq);
V  = DAE.V(Exload.bus);
Tp = Exload.con(:,5);
Tq = Exload.con(:,6);
as = Exload.con(:,7);
at = Exload.con(:,8);
bs = Exload.con(:,9);
bt = Exload.con(:,10);
P0 = Exload.dat(:,1);
Q0 = Exload.dat(:,2);
V0 = Exload.dat(:,3);

switch flag
 case 0 % initialization

  % parameter initialization
  Exload.dat(:,3) = V;
  for i = 1:Exload.n
    idx = find(PQ.bus == Exload.bus(i));
    if isempty(idx)
      fm_disp(['No PQ load found for initializing Exponential ', ...
	       'Recovery Load at bus #',int2str(Exload.bus(i))])
    else
      Exload.dat(i,1) = PQ.con(idx,4);
      Exload.dat(i,2) = PQ.con(idx,5);
      PQ.con(idx,:) = [];
      PQ.bus(idx) = [];
      PQ.n = PQ.n - 1;
      PQ.P0(idx) = [];
      PQ.Q0(idx) = [];
    end
  end

  % state variable initialization
  DAE.x(Exload.xp) = 0;
  xp = DAE.x(Exload.xp);
  DAE.x(Exload.xq) = 0;
  xq = DAE.x(Exload.xq);

  % message
  fm_disp('Initialization of exponential recovery loads completed.')

 case 1 % algebraic equations

  DAE.gp = DAE.gp + sparse(Exload.bus,1,xp./Tp + ...
			   P0.*(V./V0).^at,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Exload.bus,1,xq./Tq + ...
			   Q0.*(V./V0).^bt,Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J12 = DAE.J12 + sparse(Exload.bus,Exload.bus,P0.* ...
			     (V./V0).^(at-1).*at./V0,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Exload.bus,Exload.bus,Q0.* ...
			     (V./V0).^(bt-1).*bt./V0,Bus.n,Bus.n);

 case 3 % differential equations

  DAE.f(Exload.xp) = -xp./Tp + P0.*(V./V0).^as - P0.*(V./V0).^at;
  DAE.f(Exload.xq) = -xq./Tq + Q0.*(V./V0).^bs - Q0.*(V./V0).^bt;

 case 4 % state variable Jacobians

  % DAE.Fx
  DAE.Fx = DAE.Fx + sparse(Exload.xp,Exload.xp,-1./Tp,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Exload.xq,Exload.xq,-1./Tq,DAE.n,DAE.n);

  % DAE.Fy
  DAE.Fy = DAE.Fy + sparse(Exload.xp,Exload.bus+Bus.n, ...
           P0.*(V./V0).^(as-1).*as./V0-P0.*(V./V0).^(at-1).*at./V0, ...
			   DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Exload.xq,Exload.bus+Bus.n, ...
	   Q0.*(V./V0).^(bs-1).*bs./V0-Q0.*(V./V0).^(bt-1).*bt./V0, ...
			   DAE.n,2*Bus.n);

  % DAE.Gx
  DAE.Gx = DAE.Gx + sparse(Exload.bus,Exload.xp,1./Tp,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Exload.bus+Bus.n,Exload.xq,1./Tq,2*Bus.n,DAE.n);

end