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-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_varfm_dispfm_disp('Newton-Raphson Method for Power Flow Computation')if Settings.show  fm_disp(['Data file "',Path.data,File.data,'"'])endnodyn = 0;% these computations are needed only the first time the power flow is runif ~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 computationsend% memory allocation for equations and Jacobiansif (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  endelse  % 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);  endend% check PF solver methodFDPF = 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);  endendswitch 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')endif 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);  endendswitch Settings.distrsw case 1  fm_disp('Distributed slack bus model') case 0  fm_disp('Single slack bus model')endif FDPF > 1, fm_call('fdpf'), enditer_max = Settings.lfmit;if iter_max < 2, iter_max = 2; enditeration = 0;tol = Settings.lftol;err_max = tol+1;% Graphical settingsfm_status('pf','init',iter_max,{'r'},{'-'},{Theme.color11})%  Newton-Raphson routinet0 = 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)  endendSettings.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 busesif 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 computationif Line.n, fm_lf(1), endDAE.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 indicizationif nodyn == 1; DAE.x = []; DAE.f = []; DAE.n = 0; enddynordold = 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  endend%  build cell arrays of variable namesif (Bus.n>0), fm_idx(1), endif (DAE.n>0), fm_idx(2), endif (Settings.vs == 1), fm_idx(3), end% initializations of state variables and componentsif Settings.static  fm_disp('* * * Dynamic components are not initialized * * *')endSettings.init = 1;fm_synitfm_excinfm_call('0');% start up of induction machinesif Mot.n > 0, fm_mot(5), end% power flow result visualizationfm_disp(['Power Flow completed in ',num2str(Settings.lftime),' s'])if Settings.showlf == 1 | Fig.stat  fm_stat;else  if Settings.beep, beep, endend% initialization of all equations & JacobiansDAE.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);endfm_call('i');% build structure "Snapshot"if isempty(Settings.t0) & Fig.main  hdl = findobj(Fig.main,'Tag','EditText3');  Settings.t0 = str2num(get(hdl,'String'));endif ~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);endfm_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