fm_pq.m

一个较好的MATLAB潮流程序

function fm_pq(flag)
% FM_PQ define PQ loads
%
% FM_PQ(FLAG)
%       FLAG  1 -> algebraic equations
%             2 -> algebraic Jacobian matrices
%
%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 PQ Bus DAE

switch flag
 case 1

  DAE.gp(PQ.bus) = PQ.con(:,4) + DAE.gp(PQ.bus);
  DAE.gq(PQ.bus) = PQ.con(:,5) + DAE.gq(PQ.bus);
  a = find(DAE.V(PQ.bus) < PQ.con(:,7) & PQ.con(:,8));
  b = find(DAE.V(PQ.bus) > PQ.con(:,6) & PQ.con(:,8));
  if ~isempty(a)
    k = PQ.bus(a);
    DAE.gp(k) = PQ.con(a,4).*DAE.V(k).*DAE.V(k)./PQ.con(a,7)./PQ.con(a,7) ...
	+ DAE.gp(k) - PQ.con(a,4);
    DAE.gq(k) = PQ.con(a,5).*DAE.V(k).*DAE.V(k)./PQ.con(a,7)./PQ.con(a,7) ...
	+ DAE.gq(k) - PQ.con(a,5);
  end
  if ~isempty(b)
    k = PQ.bus(b);
    DAE.gp(k) = PQ.con(b,4).*DAE.V(k).*DAE.V(k)./PQ.con(b,6)./PQ.con(b,6) + ...
	DAE.gp(k) - PQ.con(b,4);
    DAE.gq(k) = PQ.con(b,5).*DAE.V(k).*DAE.V(k)./PQ.con(b,6)./PQ.con(b,6) + ...
	DAE.gq(k) - PQ.con(b,5);
  end

 case 2

  a = find(DAE.V(PQ.bus) < PQ.con(:,7) & PQ.con(:,8));
  b = find(DAE.V(PQ.bus) > PQ.con(:,6) & PQ.con(:,8));
  if ~isempty(a)
    h = PQ.bus(a);
    DAE.J12(h,h)= DAE.J12(h,h) + ...
        diag(2*PQ.con(a,4).*DAE.V(h)./PQ.con(a,7)./PQ.con(a,7));
    DAE.J22(h,h)= DAE.J22(h,h) + ...
        diag(2*PQ.con(a,5).*DAE.V(h)./PQ.con(a,7)./PQ.con(a,7));
  end
  if ~isempty(b)
    h = PQ.bus(b);
    DAE.J12(h,h)= DAE.J12(h,h) + ...
        diag(2*PQ.con(b,4).*DAE.V(h)./PQ.con(b,6)./PQ.con(b,6));
    DAE.J22(h,h)= DAE.J22(h,h) + ...
        diag(2*PQ.con(b,5).*DAE.V(h)./PQ.con(b,6)./PQ.con(b,6));
  end

end