gausspf

GAUSSPF  Solves the power flow using a Gauss-Seidel method.
   [V, CONVERGED, I] = GAUSSPF(YBUS, SBUS, V0, 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, 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 struct 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] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt)
0002 %GAUSSPF  Solves the power flow using a Gauss-Seidel method.
0003 %   [V, CONVERGED, I] = GAUSSPF(YBUS, SBUS, V0, 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, and column vectors with
0007 %   the lists of bus indices for the swing bus, PV buses, and PQ buses,
0008 %   respectively. The bus voltage vector contains the set point for
0009 %   generator (including ref bus) buses, and the reference angle of the
0010 %   swing bus, as well as an initial guess for remaining magnitudes and
0011 %   angles. MPOPT is a MATPOWER options struct which can be used to
0012 %   set the termination tolerance, maximum number of iterations, and
0013 %   output options (see MPOPTION for details). Uses default options
0014 %   if this parameter is not given. Returns the final complex voltages,
0015 %   a flag which indicates whether it converged or not, and the number
0016 %   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 %   and Alberto Borghetti, University of Bologna, Italy
0024 %
0025 %   This file is part of MATPOWER.
0026 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0027 %   See https://matpower.org for more info.
0028 
0029 %% default arguments
0030 if nargin < 7
0031     mpopt = mpoption;
0032 end
0033 
0034 %% options
0035 tol     = mpopt.pf.tol;
0036 max_it  = mpopt.pf.gs.max_it;
0037 
0038 %% initialize
0039 converged = 0;
0040 i = 0;
0041 V = V0;
0042 Vm = abs(V);
0043 
0044 %% set up indexing for updating V
0045 npv = length(pv);
0046 npq = length(pq);
0047 
0048 %% evaluate F(x0)
0049 mis = V .* conj(Ybus * V) - Sbus;
0050 F = [   real(mis([pv; pq]));
0051         imag(mis(pq))   ];
0052 
0053 %% check tolerance
0054 normF = norm(F, inf);
0055 if mpopt.verbose > 1
0056     fprintf('\n it    max P & Q mismatch (p.u.)');
0057     fprintf('\n----  ---------------------------');
0058     fprintf('\n%3d        %10.3e', i, normF);
0059 end
0060 if normF < tol
0061     converged = 1;
0062     if mpopt.verbose > 1
0063         fprintf('\nConverged!\n');
0064     end
0065 end
0066 
0067 %% do Gauss-Seidel iterations
0068 while (~converged && i < max_it)
0069     %% update iteration counter
0070     i = i + 1;
0071 
0072     %% update voltage
0073     %% at PQ buses
0074     for k = pq(1:npq)'
0075         V(k) =  V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k);
0076     end
0077 
0078     %% at PV buses
0079     if npv
0080         for k = pv(1:npv)'
0081             Sbus(k) = real(Sbus(k)) + 1j * imag( V(k) .* conj(Ybus(k,:) * V));
0082             V(k) =  V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k);
0083 %           V(k) = Vm(k) * V(k) / abs(V(k));
0084         end
0085         V(pv) = Vm(pv) .* V(pv) ./ abs(V(pv));
0086     end
0087 
0088     %% evalute F(x)
0089     mis = V .* conj(Ybus * V) - Sbus;
0090     F = [   real(mis(pv));
0091             real(mis(pq));
0092             imag(mis(pq))   ];
0093 
0094     %% check for convergence
0095     normF = norm(F, inf);
0096     if mpopt.verbose > 1
0097         fprintf('\n%3d        %10.3e', i, normF);
0098     end
0099     if normF < tol
0100         converged = 1;
0101         if mpopt.verbose
0102             fprintf('\nGauss-Seidel power flow converged in %d iterations.\n', i);
0103         end
0104     end
0105 end
0106 
0107 if mpopt.verbose
0108     if ~converged
0109         fprintf('\nGauss-Seidel power flow did not converge in %d iterations.\n', i);
0110     end
0111 end