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 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] = 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 vector 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 % $Id: gausspf.m,v 1.10 2010/04/26 19:45:25 ray Exp $ 0022 % by Ray Zimmerman, PSERC Cornell 0023 % and Alberto Borghetti, University of Bologna, Italy 0024 % Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC) 0025 % 0026 % This file is part of MATPOWER. 0027 % See http://www.pserc.cornell.edu/matpower/ for more info. 0028 % 0029 % MATPOWER is free software: you can redistribute it and/or modify 0030 % it under the terms of the GNU General Public License as published 0031 % by the Free Software Foundation, either version 3 of the License, 0032 % or (at your option) any later version. 0033 % 0034 % MATPOWER is distributed in the hope that it will be useful, 0035 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0037 % GNU General Public License for more details. 0038 % 0039 % You should have received a copy of the GNU General Public License 0040 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0041 % 0042 % Additional permission under GNU GPL version 3 section 7 0043 % 0044 % If you modify MATPOWER, or any covered work, to interface with 0045 % other modules (such as MATLAB code and MEX-files) available in a 0046 % MATLAB(R) or comparable environment containing parts covered 0047 % under other licensing terms, the licensors of MATPOWER grant 0048 % you additional permission to convey the resulting work. 0049 0050 %% default arguments 0051 if nargin < 7 0052 mpopt = mpoption; 0053 end 0054 0055 %% options 0056 tol = mpopt(2); 0057 max_it = mpopt(5); 0058 verbose = mpopt(31); 0059 0060 %% initialize 0061 converged = 0; 0062 i = 0; 0063 V = V0; 0064 Vm = abs(V); 0065 0066 %% set up indexing for updating V 0067 npv = length(pv); 0068 npq = length(pq); 0069 0070 %% evaluate F(x0) 0071 mis = V .* conj(Ybus * V) - Sbus; 0072 F = [ real(mis([pv; pq])); 0073 imag(mis(pq)) ]; 0074 0075 %% check tolerance 0076 normF = norm(F, inf); 0077 if verbose > 0 0078 fprintf('(Gauss-Seidel)\n'); 0079 end 0080 if verbose > 1 0081 fprintf('\n it max P & Q mismatch (p.u.)'); 0082 fprintf('\n---- ---------------------------'); 0083 fprintf('\n%3d %10.3e', i, normF); 0084 end 0085 if normF < tol 0086 converged = 1; 0087 if verbose > 1 0088 fprintf('\nConverged!\n'); 0089 end 0090 end 0091 0092 %% do Gauss-Seidel iterations 0093 while (~converged && i < max_it) 0094 %% update iteration counter 0095 i = i + 1; 0096 0097 %% update voltage 0098 %% at PQ buses 0099 for k = pq(1:npq)' 0100 V(k) = V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k); 0101 end 0102 0103 %% at PV buses 0104 if npv 0105 for k = pv(1:npv)' 0106 Sbus(k) = real(Sbus(k)) + 1j * imag( V(k) .* conj(Ybus(k,:) * V)); 0107 V(k) = V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k); 0108 % V(k) = Vm(k) * V(k) / abs(V(k)); 0109 end 0110 V(pv) = Vm(pv) .* V(pv) ./ abs(V(pv)); 0111 end 0112 0113 %% evalute F(x) 0114 mis = V .* conj(Ybus * V) - Sbus; 0115 F = [ real(mis(pv)); 0116 real(mis(pq)); 0117 imag(mis(pq)) ]; 0118 0119 %% check for convergence 0120 normF = norm(F, inf); 0121 if verbose > 1 0122 fprintf('\n%3d %10.3e', i, normF); 0123 end 0124 if normF < tol 0125 converged = 1; 0126 if verbose 0127 fprintf('\nGauss-Seidel power flow converged in %d iterations.\n', i); 0128 end 0129 end 0130 end 0131 0132 if verbose 0133 if ~converged 0134 fprintf('\nGauss-Seidel power did not converge in %d iterations.\n', i); 0135 end 0136 end