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