fm_jimma.m

电力系统的psat

function  fm_jimma(flag)
% FM_JIMMA defines Mixed Loads
%
% FM_JIMMA(FLAG)
%       FLAG = 0 initialization
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%
%Author:    Federico Milano
%Date:      23-May-2005
%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 Jimma DAE Bus Settings PQ

x = DAE.x(Jimma.x);
V1 = DAE.V(Jimma.bus);
Tf = Jimma.con(:,5);
Plz = Jimma.con(:,6);
Pli = Jimma.con(:,7);
Plp = Jimma.con(:,8);
Qlz = Jimma.con(:,9);
Qli = Jimma.con(:,10);
Qlp = Jimma.con(:,11);
Kv = Jimma.con(:,12);
V1_0 = Jimma.dat(:,1);

switch flag
 case 0 % initialization

  for k = 1:length(Jimma.n)
    idx = find(PQ.bus == Jimma.bus(k));
    if isempty(idx)
      fm_disp(['No PQ load found for initializing Jimma''s load ', ...
               'at bus #',int2str(Jimma.bus(k))])
      Jimma.con(k,6) = 0;
      Jimma.con(k,7) = 0;
      Jimma.con(k,8) = 0;
      Jimma.con(k,9) = 0;
      Jimma.con(k,10) = 0;
      Jimma.con(k,11) = 0;      
    else
      PQ.con(idx,4) = PQ.con(idx,4) - ...
          PQ.P0(idx)*sum(Jimma.con(k,[6:8]))/100;
      PQ.con(idx,5) = PQ.con(idx,5) - ...
          PQ.Q0(idx)*sum(Jimma.con(k,[9:11]))/100;
      Jimma.con(k,6) = Jimma.con(k,6)*PQ.P0(idx)/100;
      Jimma.con(k,7) = Jimma.con(k,7)*PQ.P0(idx)/100;
      Jimma.con(k,8) = Jimma.con(k,8)*PQ.P0(idx)/100;
      Jimma.con(k,9) = Jimma.con(k,9)*PQ.Q0(idx)/100;
      Jimma.con(k,10) = Jimma.con(k,10)*PQ.Q0(idx)/100;
      Jimma.con(k,11) = Jimma.con(k,11)*PQ.Q0(idx)/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

  %check time constants
  idx = find(Tf == 0);
  if idx
    jimmawarn(idx, 'Time constant Tf cannot be zero. Tf = 0.001 s will be used.'),
  end
  Jimma.con(idx,5) = 0.001;

  %variable initialization
  DAE.x(Jimma.x) = -V1./Tf;
  x = DAE.x(Jimma.x);
  Jimma.dat(:,1) = V1;

  %check limits
  fm_disp('Initialization of Jimma''s loads completed.')

 case 1 % algebraic equations

  DAE.gp = DAE.gp + sparse(Jimma.bus,1,Plz.*(V1./V1_0).^2 + ...
                           Pli.*(V1./V1_0)+Plp,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Jimma.bus,1,Qlz.*(V1./V1_0).^2 + ...
                           Qli.*(V1./V1_0)+Qlp+Kv.*(x+V1./Tf),Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J12 = DAE.J12 + sparse(Jimma.bus,Jimma.bus, ...
                             2.*Plz.*V1./V1_0.^2+Pli./V1_0,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Jimma.bus,Jimma.bus, ...
                             2.*Qlz.*V1./V1_0.^2+Qli./V1_0+Kv./Tf,Bus.n,Bus.n);

 case 3 % differential equations

  DAE.f(Jimma.x) = (-V1./Tf-x)./Tf;

 case 4 % state variable Jacobians

  % DAE.Fx
  DAE.Fx = DAE.Fx + sparse(Jimma.x,Jimma.x,-1./Tf,DAE.n,DAE.n);

  % DAE.Fy
  DAE.Fy = DAE.Fy + sparse(Jimma.x,Jimma.bus+Bus.n,-1./Tf./Tf,DAE.n,2*Bus.n);

  % DAE.Gx
  DAE.Gx = DAE.Gx + sparse(Jimma.bus+Bus.n,Jimma.x,Kv,2*Bus.n,DAE.n);

end


% -------------------------------------------------------------------
% function for creating warning messages
function jimmawarn(idx, msg)
fm_disp(strcat('Warning: Jimma''s Load #',int2str(idx),msg))