fm_mn.m

电力系统的psat

function fm_mn(flag)
% FM_MN defines Monomial Loads
%
% FM_MN(FLAG)
%       FLAG = 0 initialization
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic 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 Mn PQ DAE Bus Settings

V = DAE.V(Mn.bus);

switch flag
 case 0
  i = find(Mn.con(:,8));
  if i
    for j = 1:length(i)
      k = i(j);
      idx = find(PQ.bus == Mn.bus(k));
      if isempty(idx),
	fm_disp(['No PQ load found for initializing monomial ', ...
		 'load at bus #',int2str(Mn.bus(k))])
      else
	PQ.con(idx,4) = PQ.con(idx,4) - ...
	    PQ.P0(idx)*Mn.con(k,4)/100;
	PQ.con(idx,5) = PQ.con(idx,5) - ...
	    PQ.Q0(idx)*Mn.con(k,5)/100;
	Mn.con(k,4)  = Mn.con(k,4)*PQ.P0(idx)/(V(k)^Mn.con(k,6))/100;
	Mn.con(k,5)  = Mn.con(k,5)*PQ.Q0(idx)/(V(k)^Mn.con(k,7))/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
  end
  fm_disp('Initialization of monomial loads completed.')
 case 1
  if Settings.init
    DAE.gp = DAE.gp + sparse(Mn.bus,1,Mn.con(:,4).*V.^Mn.con(:,6),Bus.n,1);
    DAE.gq = DAE.gq + sparse(Mn.bus,1,Mn.con(:,5).*V.^Mn.con(:,7),Bus.n,1);
  elseif Mn.init
    i = Mn.init;
    DAE.gp = DAE.gp + sparse(Mn.bus(i),1,Mn.con(i,4).*V(i).^Mn.con(i,6),Bus.n,1);
    DAE.gq = DAE.gq + sparse(Mn.bus(i),1,Mn.con(i,5).*V(i).^Mn.con(i,7),Bus.n,1);
  end
 case 2
  if Settings.init
    DAE.J12 = DAE.J12 + sparse(Mn.bus,Mn.bus,Mn.con(:,4).* ...
			       Mn.con(:,6).*V.^(Mn.con(:,6)-1),Bus.n,Bus.n);
    DAE.J22 = DAE.J22 + sparse(Mn.bus,Mn.bus,Mn.con(:,5).* ...
			       Mn.con(:,7).*V.^(Mn.con(:,7)-1),Bus.n,Bus.n);
  elseif Mn.init
    i = Mn.init;
    DAE.J12 = DAE.J12 + sparse(Mn.bus(i),Mn.bus(i),Mn.con(i,4).* ...
			       Mn.con(i,6).*V(i).^(Mn.con(i,6)-1),Bus.n,Bus.n);
    DAE.J22 = DAE.J22 + sparse(Mn.bus(i),Mn.bus(i),Mn.con(i,5).* ...
			       Mn.con(i,7).*V(i).^(Mn.con(i,7)-1),Bus.n,Bus.n);
  end
end