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:    fmilano@thunderbox.uwaterloo.ca
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

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
  % bus type initialization
  fm_ncomp
  if ~Settings.ok, return, end
  % report components parameters to system bases
  if Settings.conv, fm_base, end
  % create the admittance matrix
  fm_y
  % 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;
  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
else  % no dynamic elements
  nodyn = 1;
  DAE.n = 1;
  DAE.f = 0;
  DAE.x = 0;
  DAE.Fx = 1;
  DAE.Fy = sparse(1,2*Bus.n);
  DAE.Gx = sparse(2*Bus.n,1);
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);
  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
  DAE.Fk = sparse(DAE.n,1);
  DAE.Gk = sparse(2*Bus.n,1);
  Jrow = sparse(1,DAE.n+SW.bus,1,1,DAE.n+2*Bus.n+1);
  DAE.kg = 0;
  if ~swgamma(SW,'sum')
    fm_disp('Slack buses have zero power loss participation factor.')
    fm_disp('Single slack bus model will be used')
    Setting.distrsw = 0;
  end
  if totp(SW) == 0
    SW = setpg(totp(PQ)-totp(PV),1);
    fm_disp('Slack buses have zero generated power.')
    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
  Gkcall(SW);
  Gkcall(PV);
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;
err_vec = [];
alfatry = 1;
alfa = 1;
safety = 0.9;
pgrow = -0.2;
pshrnk = -0.25;
try
  errcon = (5/safety)^(1/pgrow);
catch
  errcon = 1.89e-4;
end

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

%  Newton-Raphson routine
t0 = clock;

while (err_max > tol) & (iteration <= iter_max) & (alfa > 1e-5)

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

  if FDPF == 1   % Standard NR technique

    refreshJac;

    % component calls
    if Settings.distrsw
      fm_call('kgpf');
      glambda(SW,1,DAE.kg);
      DAE.g(Bus.n+SW.bus) = 0;

      Gbus = Bus.n+[SW.bus;PV.bus];
      DAE.Jlfv(Gbus,:) = 0;
      DAE.Jlfv(:,Gbus) = 0;
      DAE.Jlfv(:,SW.bus) = 0;
      DAE.Jlfv(Gbus,Gbus) = speye(PV.n+SW.n);
      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],[DAE.Fk;DAE.Gk];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(1+DAE.n+Bus.n:DAE.n+2*Bus.n);
      DAE.kg = DAE.kg + inc(end);

    elseif Settings.robust

      xold = DAE.x;
      aold = DAE.a;
      Vold = DAE.V;

      while 1

        refreshJac;
        fm_call('l');
        refreshGen(nodyn);
        k1 = -alfa*[DAE.Fx, DAE.Fy; DAE.Gx, DAE.Jlfv]\[DAE.f; DAE.g];

        DAE.x = xold + 0.5*k1(1:DAE.n);
        DAE.a = aold + 0.5*k1(1+DAE.n:Bus.n+DAE.n);
        DAE.V = Vold + 0.5*k1(1+DAE.n+Bus.n:DAE.n+2*Bus.n);
        refreshJac;
        fm_call('l');
        refreshGen(nodyn);
        k2 = -alfa*[DAE.Fx, DAE.Fy; DAE.Gx, DAE.Jlfv]\[DAE.f; DAE.g];

        DAE.x = xold + 0.5*k2(1:DAE.n);
        DAE.a = aold + 0.5*k2(1+DAE.n:Bus.n+DAE.n);
        DAE.V = Vold + 0.5*k2(1+DAE.n+Bus.n:DAE.n+2*Bus.n);
        refreshJac;
        fm_call('l');
        refreshGen(nodyn);
        k3 = -alfa*[DAE.Fx, DAE.Fy; DAE.Gx, DAE.Jlfv]\[DAE.f; DAE.g];

        DAE.x = xold + k3(1:DAE.n);
        DAE.a = aold + k3(1+DAE.n:Bus.n+DAE.n);
        DAE.V = Vold + k3(1+DAE.n+Bus.n:DAE.n+2*Bus.n);
        refreshJac;
        fm_call('l');
        refreshGen(nodyn);
        k4 = -alfa*[DAE.Fx, DAE.Fy; DAE.Gx, DAE.Jlfv]\[DAE.f; DAE.g];

        % compute RK4 increment of variables
        inc = (k1+2*(k2+k3)+k4)/6;

        % to estimate RK error, use the RK2:Dy and RK4:Dy.
        yerr = max(abs(abs(k2)-abs(inc)));
        if yerr > 0.01
          alfa = max(0.985*alfa,0.9);
        else
          alfa = min(1.015*alfa,1.1);
        end

        DAE.x = xold + inc(1:DAE.n);
        DAE.a = aold + inc(1+DAE.n:Bus.n+DAE.n);
        DAE.V = Vold + inc(1+DAE.n+Bus.n:DAE.n+2*Bus.n);
        break

      end
      %disp(['alpha = ',num2str(alfa)])

    else % standard Newton-Raphson method

      % complete Jacobian matrix
      fm_call('l');
      refreshGen(nodyn)
      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(1+DAE.n+Bus.n: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));
  err_vec(iteration) = err_max;
  
  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

  % stop if the error increases too much 
  if iteration > 4
    if err_vec(iteration) > 1000*err_vec(1)
      fm_disp('The error is increasing too much.')
      fm_disp('Convergence is likely not reachable.')
      break
    end
  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 injections and consumptions at network buses

% 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;
% take away PQ generator powers
rmpqgen(PQ)

%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;
SW = setpg(SW,'all',Bus.Pg(SW.bus));
% add PQ generator powers
addpqgen(PQ)

% memory allocation for dynamic variables & state variables indicization
if nodyn == 1; DAE.x = []; DAE.f = []; DAE.n = 0; end
DAE.npf = DAE.n;
fm_dynidx;
if (DAE.n~=0)
  DAE.f = [DAE.f; ones(DAE.n-DAE.npf,1)];  % differential equations
  DAE.x = [DAE.x; ones(DAE.n-DAE.npf,1)];  % state variables
  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

%  build cell arrays of variable names
fm_idx(1)
if (Settings.vs == 1), fm_idx(2), 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
% refresh PQ powers
% PQ = pqactual(PQ,'all');

% 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
refreshJac
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 & Bus.n <= 5000 & DAE.n < 5000
  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);
else
  if Bus.n > 5000
    fm_disp(['Snapshots are disabled for networks with more than 5000 ' ...
             'buses.'])
  end
  if DAE.n > 5000
    fm_disp(['Snapshots are disabled for networks with more than 5000 ' ...
             'state variables.'])
  end  
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


% -----------------------------------------------
function refreshJac

global DAE Bus

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);

% -----------------------------------------------
function refreshGen(nodyn)

global DAE Bus SW PV

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
 
% 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;
  DAE.Jlfv(k,k) = speye(length(k));
  DAE.Jlfv(k+Bus.n,k+Bus.n) = speye(length(k));
  DAE.g(k) = 0;
  DAE.g(k+Bus.n) = 0;
  DAE.V(k) = 1e-6;
  DAE.a(k) = 0;
end