fdpf

FDPF  Solves the power flow using a fast decoupled method.
   [V, CONVERGED, I] = FDPF(YBUS, SBUS, V0, BP, BPP, REF, PV, PQ, MPOPT)
   solves for bus voltages given the full system admittance matrix (for
   all buses), the complex bus power injection vector (for all buses),
   the initial vector of complex bus voltages, the FDPF matrices B prime
   and B double prime, and column vectors with the lists of bus indices
   for the swing bus, PV buses, and PQ buses, respectively. The bus voltage
   vector contains the set point for generator (including ref bus)
   buses, and the reference angle of the swing bus, as well as an initial
   guess for remaining magnitudes and angles. MPOPT is a MATPOWER options
   vector which can be used to set the termination tolerance, maximum
   number of iterations, and output options (see MPOPTION for details).
   Uses default options if this parameter is not given. Returns the
   final complex voltages, a flag which indicates whether it converged
   or not, and the number of iterations performed.

   See also RUNPF.

0001 function [V, converged, i] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt)
0002 %FDPF  Solves the power flow using a fast decoupled method.
0003 %   [V, CONVERGED, I] = FDPF(YBUS, SBUS, V0, BP, BPP, REF, PV, PQ, MPOPT)
0004 %   solves for bus voltages given the full system admittance matrix (for
0005 %   all buses), the complex bus power injection vector (for all buses),
0006 %   the initial vector of complex bus voltages, the FDPF matrices B prime
0007 %   and B double prime, and column vectors with the lists of bus indices
0008 %   for the swing bus, PV buses, and PQ buses, respectively. The bus voltage
0009 %   vector contains the set point for generator (including ref bus)
0010 %   buses, and the reference angle of the swing bus, as well as an initial
0011 %   guess for remaining magnitudes and angles. MPOPT is a MATPOWER options
0012 %   vector which can be used to set the termination tolerance, maximum
0013 %   number of iterations, and output options (see MPOPTION for details).
0014 %   Uses default options if this parameter is not given. Returns the
0015 %   final complex voltages, a flag which indicates whether it converged
0016 %   or not, and the number of iterations performed.
0017 %
0018 %   See also RUNPF.
0019 
0020 %   MATPOWER
0021 %   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
0022 %   by Ray Zimmerman, PSERC Cornell
0023 %
0024 %   This file is part of MATPOWER.
0025 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0026 %   See https://matpower.org for more info.
0027 
0028 %% default arguments
0029 if nargin < 7
0030     mpopt = mpoption;
0031 end
0032 
0033 %% options
0034 tol     = mpopt.pf.tol;
0035 max_it  = mpopt.pf.fd.max_it;
0036 lu_vec  = have_feature('lu_vec');
0037 
0038 %% initialize
0039 converged = 0;
0040 i = 0;
0041 V = V0;
0042 Va = angle(V);
0043 Vm = abs(V);
0044 
0045 %% set up indexing for updating V
0046 npv = length(pv);
0047 npq = length(pq);
0048 
0049 %% evaluate initial mismatch
0050 mis = (V .* conj(Ybus * V) - Sbus(Vm)) ./ Vm;
0051 P = real(mis([pv; pq]));
0052 Q = imag(mis(pq));
0053 
0054 %% check tolerance
0055 normP = norm(P, inf);
0056 normQ = norm(Q, inf);
0057 if mpopt.verbose > 1
0058     fprintf('\niteration     max mismatch (p.u.)  ');
0059     fprintf('\ntype   #        P            Q     ');
0060     fprintf('\n---- ----  -----------  -----------');
0061     fprintf('\n  -  %3d   %10.3e   %10.3e', i, normP, normQ);
0062 end
0063 if normP < tol && normQ < tol
0064     converged = 1;
0065     if mpopt.verbose > 1
0066         fprintf('\nConverged!\n');
0067     end
0068 end
0069 
0070 %% reduce B matrices
0071 Bp = Bp([pv; pq], [pv; pq]);
0072 Bpp = Bpp(pq, pq);
0073 
0074 %% factor B matrices
0075 if lu_vec
0076     [Lp,  Up,  pp,  qp ] = lu(Bp,  'vector');
0077     [Lpp, Upp, ppp, qpp] = lu(Bpp, 'vector');
0078     [junk, iqp ] = sort(qp);
0079     [junk, iqpp] = sort(qpp);
0080     % [~, iqp ] = sort(qp);
0081     % [~, iqpp] = sort(qpp);
0082 else
0083     [Lp, Up, Pp] = lu(Bp);
0084     [Lpp, Upp, Ppp] = lu(Bpp);
0085 end
0086 
0087 %% do P and Q iterations
0088 while (~converged && i < max_it)
0089     %% update iteration counter
0090     i = i + 1;
0091 
0092     %%-----  do P iteration, update Va  -----
0093     if lu_vec
0094         dVa = -( Up \  (Lp \ P(pp)) );
0095         dVa = dVa(iqp);
0096     else
0097         dVa = -( Up \  (Lp \ (Pp * P)));
0098     end
0099 
0100     %% update voltage
0101     Va([pv; pq]) = Va([pv; pq]) + dVa;
0102     V = Vm .* exp(1j * Va);
0103 
0104     %% evalute mismatch
0105     mis = (V .* conj(Ybus * V) - Sbus(Vm)) ./ Vm;
0106     P = real(mis([pv; pq]));
0107     Q = imag(mis(pq));
0108     
0109     %% check tolerance
0110     normP = norm(P, inf);
0111     normQ = norm(Q, inf);
0112     if mpopt.verbose > 1
0113         fprintf('\n  P  %3d   %10.3e   %10.3e', i, normP, normQ);
0114     end
0115     if normP < tol && normQ < tol
0116         converged = 1;
0117         if mpopt.verbose
0118             fprintf('\nFast-decoupled power flow converged in %d P-iterations and %d Q-iterations.\n', i, i-1);
0119         end
0120         break;
0121     end
0122 
0123     %%-----  do Q iteration, update Vm  -----
0124     if lu_vec
0125         dVm = -( Upp \ (Lpp \ Q(ppp)) );
0126         dVm = dVm(iqpp);
0127     else
0128         dVm = -( Upp \ (Lpp \ (Ppp * Q)) );
0129     end
0130 
0131     %% update voltage
0132     Vm(pq) = Vm(pq) + dVm;
0133     V = Vm .* exp(1j * Va);
0134 
0135     %% evalute mismatch
0136     mis = (V .* conj(Ybus * V) - Sbus(Vm)) ./ Vm;
0137     P = real(mis([pv; pq]));
0138     Q = imag(mis(pq));
0139     
0140     %% check tolerance
0141     normP = norm(P, inf);
0142     normQ = norm(Q, inf);
0143     if mpopt.verbose > 1
0144         fprintf('\n  Q  %3d   %10.3e   %10.3e', i, normP, normQ);
0145     end
0146     if normP < tol && normQ < tol
0147         converged = 1;
0148         if mpopt.verbose
0149             fprintf('\nFast-decoupled power flow converged in %d P-iterations and %d Q-iterations.\n', i, i);
0150         end
0151         break;
0152     end
0153 end
0154 
0155 if mpopt.verbose
0156     if ~converged
0157         fprintf('\nFast-decoupled power flow did not converge in %d iterations.\n', i);
0158     end
0159 end