fm_nrlf.m

电力系统的psat

function    conv = fm_nrlf(iter_max, tol)
% FM_NRLF solve power flow with locked ste variables
%
% CONV = FM_NRLF(ITERMAX,TOL)
%       ITERMAX = max number of iterations
%       TOL = convergence tolerance
%       CONV = 1 if convergence reached, 0 otherwise
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    11-Sep-2003
%Version:   1.1.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 DAE Bus Line PV SW Settings

conv = 1;
iteration = 0;
inc = ones(2*Bus.n,1);
rbus = Settings.refbus;
Vguess = ones(Bus.n,1);

% look for PV and slack generators
if PV.n | SW.n,
  Vguess([PV.bus; SW.bus]) = DAE.V([PV.bus; SW.bus]);
end

% initialize bus voltages
DAE.V = Vguess;
y2 = [DAE.a; DAE.V];

% Newton-Raphson routine with locked state variables
while max(abs(inc)) > tol & iteration < iter_max

  if isempty(Line.Y)
    DAE.gp = zeros(Bus.n,1);
    DAE.gq = zeros(Bus.n,1);
    if Settings.octave
      DAE.J11 = zeros(Bus.n,Bus.n);
      DAE.J21 = zeros(Bus.n,Bus.n);
      DAE.J12 = zeros(Bus.n,Bus.n);
      DAE.J22 = zeros(Bus.n,Bus.n);
    else
      DAE.J11 = sparse(Bus.n,Bus.n);
      DAE.J21 = sparse(Bus.n,Bus.n);
      DAE.J12 = sparse(Bus.n,Bus.n);
      DAE.J22 = sparse(Bus.n,Bus.n);
    end
  end

  fm_call('n')  % call algebraic functions

  %DAE.Jlfv(rbus,:) = 0;
  %DAE.Jlfv(:,rbus) = 0;
  %DAE.Jlfv(rbus,rbus) = 1;
  %DAE.g(rbus) = 0;

  % check for islanded buses
  if ~isempty(Bus.island)
    k = Bus.island;
    DAE.Jlfv(k,:) = 0;
    DAE.Jlfv(:,k) = 0;
    DAE.Jlfv(:,k+Bus.n) = 0;
    DAE.Jlfv(k+Bus.n,:) = 0;
    if Settings.octave
      DAE.Jlfv(k,k) = eye(length(Bus.island));
      DAE.Jlfv(k+Bus.n,k+Bus.n) = eye(length(Bus.island));
    else
      DAE.Jlfv(k,k) = speye(length(Bus.island));
      DAE.Jlfv(k+Bus.n,k+Bus.n) = speye(length(Bus.island));
    end
    DAE.g(k) = 0;
    DAE.g(k+Bus.n) = 0;
    DAE.V(k) = 1e-6;
    DAE.a(k) = 0;
  end

  inc = -DAE.Jlfv\DAE.g; 
  y2 = y2 + inc;
  DAE.a = y2(1:Bus.n);
  DAE.V = y2(Bus.n+1:2*Bus.n);
  iteration = iteration + 1;

end

DAE.a(find(DAE.V <= 1e-6)) = 0;

% unwrap voltage phases
% [only when solution is critical because of very low voltage values]
for i = 1:Bus.n
  if abs(DAE.V(i)) < 1e-4 & abs(DAE.a(i)) > 2*pi
    while abs(DAE.a(i)) > 2*pi
      if DAE.a(i) < 0
	DAE.a(i) = DAE.a(i) + 2*pi;
      else
	DAE.a(i) = DAE.a(i) - 2*pi;
      end
    end
  end
end

% message of end of operations
if iteration >= iter_max,
  fm_disp('Solution of algebraic equations failed.')
  conv = 0;
else
  fm_disp(['Solution of algebraic equations completed in ', ...
	   num2str(iteration),' iterations.'])
end

% update time derivatives of state variables
fm_call('3')