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-2015 by Power System Engineering Research Center (PSERC)
0022 %   by Ray Zimmerman, PSERC Cornell
0023 %   and Alberto Borghetti, University of Bologna, Italy
0024 %
0025 %   $Id: gausspf.m 2644 2015-03-11 19:34:22Z ray $
0026 %
0027 %   This file is part of MATPOWER.
0028 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0029 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0030 
0031 %% default arguments
0032 if nargin < 7
0033     mpopt = mpoption;
0034 end
0035 
0036 %% options
0037 tol     = mpopt.pf.tol;
0038 max_it  = mpopt.pf.gs.max_it;
0039 
0040 %% initialize
0041 converged = 0;
0042 i = 0;
0043 V = V0;
0044 Vm = abs(V);
0045 
0046 %% set up indexing for updating V
0047 npv = length(pv);
0048 npq = length(pq);
0049 
0050 %% evaluate F(x0)
0051 mis = V .* conj(Ybus * V) - Sbus;
0052 F = [   real(mis([pv; pq]));
0053         imag(mis(pq))   ];
0054 
0055 %% check tolerance
0056 normF = norm(F, inf);
0057 if mpopt.verbose > 1
0058     fprintf('\n it    max P & Q mismatch (p.u.)');
0059     fprintf('\n----  ---------------------------');
0060     fprintf('\n%3d        %10.3e', i, normF);
0061 end
0062 if normF < tol
0063     converged = 1;
0064     if mpopt.verbose > 1
0065         fprintf('\nConverged!\n');
0066     end
0067 end
0068 
0069 %% do Gauss-Seidel iterations
0070 while (~converged && i < max_it)
0071     %% update iteration counter
0072     i = i + 1;
0073 
0074     %% update voltage
0075     %% at PQ buses
0076     for k = pq(1:npq)'
0077         V(k) =  V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k);
0078     end
0079 
0080     %% at PV buses
0081     if npv
0082         for k = pv(1:npv)'
0083             Sbus(k) = real(Sbus(k)) + 1j * imag( V(k) .* conj(Ybus(k,:) * V));
0084             V(k) =  V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k);
0085 %           V(k) = Vm(k) * V(k) / abs(V(k));
0086         end
0087         V(pv) = Vm(pv) .* V(pv) ./ abs(V(pv));
0088     end
0089 
0090     %% evalute F(x)
0091     mis = V .* conj(Ybus * V) - Sbus;
0092     F = [   real(mis(pv));
0093             real(mis(pq));
0094             imag(mis(pq))   ];
0095 
0096     %% check for convergence
0097     normF = norm(F, inf);
0098     if mpopt.verbose > 1
0099         fprintf('\n%3d        %10.3e', i, normF);
0100     end
0101     if normF < tol
0102         converged = 1;
0103         if mpopt.verbose
0104             fprintf('\nGauss-Seidel power flow converged in %d iterations.\n', i);
0105         end
0106     end
0107 end
0108 
0109 if mpopt.verbose
0110     if ~converged
0111         fprintf('\nGauss-Seidel power flow did not converge in %d iterations.\n', i);
0112     end
0113 end