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 http://www.pserc.cornell.edu/matpower/ 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