fm_spf.m

一个较好的MATLAB潮流程序

function fm_spf
% FM_SPF solve standard power flow by means of the NR method
%       and fast decoupled power flow (XB and BX variations)
%       with either a single or distributed slack bus model.
%
% FM_SPF
%
%see the properties of the Settings structure for options.
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    09-Jul-2003
%Version:   1.0.1
%
%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.

fm_var
fm_disp
fm_disp('Newton-Raphson Method for Power Flow Computation')
if Settings.show
  fm_disp(['Data file "',Path.data,File.data,'"'])
end
nodyn = 0;

% these computations are needed only the first time the power flow is run
if ~Settings.init
  check = fm_ncomp;     % bus type initialization
  if ~check, return, end
  % report components parameters to system bases
  if Settings.conv, fm_base, end
  fm_y                  % create Admittance Matrix
                        % create the FM_CALL FUNCTION
  if Settings.show, fm_disp('Writing file "fm_call" ...',1), end
  fm_wcall;
  Settings.refbus = SW.bus(1); % reference bus (for phase angle)
  fm_dynlf; % indicization of components used in power flow computations
end

% memory allocation for equations and Jacobians
if (DAE.n~=0)
  DAE.f = ones(DAE.n,1);  % differential equations
  DAE.x = ones(DAE.n,1);  % state variables
  fm_xfirst;
  if Settings.octave
    DAE.Fx = zeros(DAE.n,DAE.n);   % df/dx
    DAE.Fy = zeros(DAE.n,2*Bus.n); % df/dy
    DAE.Gx = zeros(2*Bus.n,DAE.n); % dg/dx
  else
    DAE.Fx = sparse(DAE.n,DAE.n);   % df/dx
    DAE.Fy = sparse(DAE.n,2*Bus.n); % df/dy
    DAE.Gx = sparse(2*Bus.n,DAE.n); % dg/dx
  end
else  % no dynamic elements
  nodyn = 1;
  DAE.n = 1;
  DAE.f = 0;
  DAE.x = 0;
  DAE.Fx = 1;
  if Settings.octave
    DAE.Fy = zeros(1,2*Bus.n);
    DAE.Gx = zeros(2*Bus.n,1);
  else
    DAE.Fy = sparse(1,2*Bus.n);
    DAE.Gx = sparse(2*Bus.n,1);
  end
end

% check PF solver method
FDPF = Settings.pfsolver;
if FDPF > 1
  if Mn.n, mn = sum(~Mn.con(:,8));  else, mn = 0; end
  if Pl.n, pl = sum(~Pl.con(:,11)); else, pl = 0; end
  ncload = mn + pl + Lines.n + (~nodyn) + Tcsc.n;
  if ncload | Settings.distrsw,
    FDPF = 1; % force using standard NR method
  else
    fm_b
    no_sw = 1:Bus.n;
    no_sw(SW.bus) = [];
    no_g = 1:Bus.n;
    no_g([SW.bus; PV.bus]) = [];
    Bp = Line.Bp(no_sw,no_sw);
    Bpp = Line.Bpp(no_g,no_g);
    [Lp, Up, Pp] = lu(Bp);
    [Lpp, Upp, Ppp] = lu(Bpp);
  end
end

switch FDPF
 case 1, fm_disp('PF solver: Newton-Raphson method')
 case 2, fm_disp('PF solver: XB fast decoupled power flow method')
 case 3, fm_disp('PF solver: BX fast decoupled power flow method')
end

if Settings.distrsw
  Jkg = sparse(DAE.n+2*Bus.n,1);
  Jrow = sparse(1,DAE.n+SW.bus,1,1,DAE.n+2*Bus.n+1);
  DAE.kg = 0;
  idx = find(SW.con(:,10) == 0);
  if length(idx) > 1
    fm_disp('More than one slack bus with zero power indication found.')
    fm_disp('Single slack bus model will be used')
    Setting.distrsw = 0;
  end
  if idx
    SW.con(1,10) = sum(PQ.con(:,4))-sum(PV.con(:,4));
    fm_disp('Slack bus with zero power direction found.')
    fm_disp('P_slack = sum(P_load)-sum(P_gen) will be used.')
    fm_disp('Only PQ loads and PV generators are used.')
    fm_disp('If there are convergence problems, use the single slack bus model.')
  end
  if SW.n,
    ksw = SW.con(:,11);
    Jkg(DAE.n+SW.bus) = -ksw.*SW.con(:,10);
  end
  if PV.n,
    kpv = PV.con(:,10);
    Jkg(DAE.n+PV.bus) = -PV.con(:,10).*PV.con(:,4);
  end
end

switch Settings.distrsw
 case 1
  fm_disp('Distributed slack bus model')
 case 0
  fm_disp('Single slack bus model')
end

if FDPF > 1, fm_call('fdpf'), end

iter_max = Settings.lfmit;
if iter_max < 2, iter_max = 2; end
iteration = 0;
tol = Settings.lftol;
err_max = tol+1;

% Graphical settings
fm_status('pf','init',iter_max,{'r'},{'-'},{Theme.color11})

%  Newton-Raphson routine
t0 = clock;

while (err_max > tol) & (iteration <= iter_max)

  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end

  if FDPF == 1   % Standard NR technique

    if Settings.octave
      DAE.gp = zeros(Bus.n,1);
      DAE.gq = zeros(Bus.n,1);
      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.gp = sparse(Bus.n,1);
      DAE.gq = sparse(Bus.n,1);
      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

    % component calls
    if Settings.distrsw
      fm_call('kgpf');
      if SW.n
        DAE.g(SW.bus) = DAE.g(SW.bus)-(1+DAE.kg*ksw).*SW.con(:,10);
        DAE.g(Bus.n+SW.bus)=0;
      end
      Gbus = Bus.n+[SW.bus;PV.bus];
      DAE.Jlfv(Gbus,:) = 0;
      DAE.Jlfv(:,Gbus) = 0;
      DAE.Jlfv(:,SW.bus) = 0;
      if Settings.octave
        DAE.Jlfv(Gbus,Gbus) = eye(PV.n+SW.n);
      else
        DAE.Jlfv(Gbus,Gbus) = speye(PV.n+SW.n);
      end
      DAE.Fy(:,SW.bus) = 0;
      DAE.Gx(SW.bus,:) = 0;
      if nodyn, DAE.Fx = 1; end
      inc = -[[DAE.Fx, DAE.Fy; DAE.Gx, DAE.Jlfv],Jkg;Jrow]\[DAE.f; DAE.g; 0];
      DAE.x = DAE.x + inc(1:DAE.n);
      DAE.a = DAE.a + inc(1+DAE.n:Bus.n+DAE.n);
      DAE.V = DAE.V + inc(DAE.n+Bus.n+1:DAE.n+2*Bus.n);
      DAE.kg = DAE.kg + inc(end);
    else
      fm_call('l');
      % complete Jacobian matrix
      DAE.Fy(:,SW.bus) = 0;
      DAE.Gx(SW.bus,:) = 0;
      DAE.Fy(:,Bus.n+SW.bus) = 0;
      DAE.Gx(Bus.n+SW.bus,:) = 0;
      DAE.Fy(:,Bus.n+PV.bus) = 0;
      DAE.Gx(Bus.n+PV.bus,:) = 0;
      if nodyn, DAE.Fx = 1; end

      inc = -[DAE.Fx, DAE.Fy; DAE.Gx, DAE.Jlfv]\[DAE.f; DAE.g];
      DAE.x = DAE.x + inc(1:DAE.n);
      DAE.a = DAE.a + inc(1+DAE.n:Bus.n+DAE.n);
      DAE.V = DAE.V + inc(DAE.n+Bus.n+1:DAE.n+2*Bus.n);
    end

  else % Fast Decoupled Power Flow

    % P-theta
    da = -(Up\(Lp\(Pp*DAE.gp(no_sw))));
    DAE.a(no_sw) = DAE.a(no_sw) + da;
    Vc = DAE.V.*exp(j*DAE.a);
    fm_call('fdpf')
    normP = norm(DAE.gp,inf);
    normQ = norm(DAE.gq,inf);
    if normP < tol & normQ < tol, break, end

    % Q-V
    dV = -(Upp\(Lpp\(Ppp*DAE.gq(no_g))));
    DAE.V(no_g) = DAE.V(no_g)+dV;
    Vc = DAE.V.*exp(j*DAE.a);
    fm_call('fdpf')
    normP = norm(DAE.gp,inf);
    normQ = norm(DAE.gq,inf);
    if normP < tol & normQ < tol, break, end

    inc = [normP; normQ];

    % recompute Bpp if some PV bus has been switched to PQ bus
    if Settings.pv2pq & strmatch('Switch',History.text{end})
      fm_disp('Recomputing Bpp matrix for FDPF')
      no_g = 1:Bus.n;
      no_g([SW.bus; PV.bus]) = [];
      Bpp = Line.Bpp(no_g,no_g);
      [Lpp, Upp, Ppp] = lu(Bpp);
    end
  end

  iteration = iteration + 1;
  err_max = max(abs(inc));

  fm_status('pf','update',[iteration, err_max],iteration)
  if Settings.show
    if err_max == Inf, err_max = 1e3; end 
    fm_disp(['Iteration = ', num2str(iteration), ...
             '     Maximum Convergency Error = ', ...
             num2str(err_max)],1)
  end

end

Settings.lftime = etime(clock,t0);
Settings.iter = iteration;
if iteration > iter_max
  fm_disp(['Reached Maximum number of iteration for NR Routine without ' ...
           'Convergence'],2)
end

% Total power injection and adsorption at network buses

% restore PV buses switched to PQ buses
if Settings.pv2pq & PV.pq.n
  PV.n = PV.n + PV.pq.n;
  PV.bus = [PV.bus; PV.pq.bus];
  PV.pq.con(:,5) = DAE.V(PV.pq.bus);
  PV.con = [PV.con; PV.pq.con(:,[1:end-1])];
  PV.pq.n = 0;
  PV.pq.bus = [];
  PV.pq.con = [];
end

% Pl and Ql computation (shunts only)
DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
fm_call('pq');
Bus.Pl = DAE.gp;
Bus.Ql = DAE.gq;
%Pg and Qg computation
if Line.n, fm_lf(1), end
DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
fm_call('1');
Bus.Pg = DAE.gp + DAE.glfp;
Bus.Qg = DAE.gq + DAE.glfq;
% reset of global variables DAE.gp and DAE.gq
%DAE.gp = zeros(Bus.n,1);
%DAE.gq = zeros(Bus.n,1);

% memory allocation for dynamic variables & state variables indicization
if nodyn == 1; DAE.x = []; DAE.f = []; DAE.n = 0; end
dynordold = DAE.n;
DAE.npf = DAE.n;
fm_dynidx;
if (DAE.n~=0)
  DAE.f = [DAE.f; ones(DAE.n-dynordold,1)];  % differential equations
  DAE.x = [DAE.x; ones(DAE.n-dynordold,1)];  % state variables
  if Settings.octave
    DAE.Fx = zeros(DAE.n,DAE.n); % state Jacobian df/dx
    DAE.Fy = zeros(DAE.n,2*Bus.n); % state Jacobian df/dy
    DAE.Gx = zeros(2*Bus.n,DAE.n); % algebraic Jacobian dg/dx
  else
    DAE.Fx = sparse(DAE.n,DAE.n); % state Jacobian df/dx
    DAE.Fy = sparse(DAE.n,2*Bus.n); % state Jacobian df/dy
    DAE.Gx = sparse(2*Bus.n,DAE.n); % algebraic Jacobian dg/dx
  end
end

%  build cell arrays of variable names
if (Bus.n>0), fm_idx(1), end
if (DAE.n>0), fm_idx(2), end
if (Settings.vs == 1), fm_idx(3), end

% initializations of state variables and components
if Settings.static
  fm_disp('* * * Dynamic components are not initialized * * *')
end
Settings.init = 1;
fm_synit
fm_excin
fm_call('0');
% start up of induction machines
if Mot.n > 0, fm_mot(5), end

% power flow result visualization
fm_disp(['Power Flow completed in ',num2str(Settings.lftime),' s'])
if Settings.showlf == 1 | Fig.stat
  fm_stat;
else
  if Settings.beep, beep, end
end

% initialization of all equations & Jacobians
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

fm_call('i');

% build structure "Snapshot"
if isempty(Settings.t0) & Fig.main
  hdl = findobj(Fig.main,'Tag','EditText3');
  Settings.t0 = str2num(get(hdl,'String'));
end

if ~Settings.locksnap
  Snapshot = struct('name','Power Flow Results', ...
                    'time',Settings.t0, ...
                    'V',DAE.V, ...
                    'ang',DAE.a, ...
                    'x',DAE.x,...
                    'Y',Line.Y, ...
                    'Pg',Bus.Pg, ...
                    'Qg',Bus.Qg, ...
                    'Pl',Bus.Pl, ...
                    'Ql',Bus.Ql, ...
                    'vfd',Syn.vf, ...
                    'pmech',Syn.pm, ...
                    'Jlf',DAE.Jlf, ...
                    'Jlfv',DAE.Jlfv, ...
                    'Fx',DAE.Fx, ...
                    'Fy',DAE.Fy, ...
                    'Gx',DAE.Gx, ...
                    'Ploss',sum(DAE.glfp), ...
                    'Qloss',sum(DAE.glfq), ...
                    'it',1);
end

fm_status('pf','close')

LIB.selbus = min(LIB.selbus,Bus.n);
if Fig.lib,
  set(findobj(Fig.lib,'Tag','Listbox1'), ...
      'String',Varname.bus, ...
      'Value',LIB.selbus);
end