NEWTONPF Solves the power flow using a full Newton's method. [V, CONVERGED, I] = NEWTONPF(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] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) 0002 %NEWTONPF Solves the power flow using a full Newton's method. 0003 % [V, CONVERGED, I] = NEWTONPF(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 if 0014 % this parameter is not given. Returns the final complex voltages, a 0015 % flag which indicates whether it converged or not, and the number of 0016 % iterations performed. 0017 % 0018 % See also RUNPF. 0019 0020 % MATPOWER 0021 % $Id: newtonpf.m,v 1.13 2010/04/26 19:45:25 ray Exp $ 0022 % by Ray Zimmerman, PSERC Cornell 0023 % Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC) 0024 % 0025 % This file is part of MATPOWER. 0026 % See http://www.pserc.cornell.edu/matpower/ for more info. 0027 % 0028 % MATPOWER is free software: you can redistribute it and/or modify 0029 % it under the terms of the GNU General Public License as published 0030 % by the Free Software Foundation, either version 3 of the License, 0031 % or (at your option) any later version. 0032 % 0033 % MATPOWER is distributed in the hope that it will be useful, 0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0036 % GNU General Public License for more details. 0037 % 0038 % You should have received a copy of the GNU General Public License 0039 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0040 % 0041 % Additional permission under GNU GPL version 3 section 7 0042 % 0043 % If you modify MATPOWER, or any covered work, to interface with 0044 % other modules (such as MATLAB code and MEX-files) available in a 0045 % MATLAB(R) or comparable environment containing parts covered 0046 % under other licensing terms, the licensors of MATPOWER grant 0047 % you additional permission to convey the resulting work. 0048 0049 %% default arguments 0050 if nargin < 7 0051 mpopt = mpoption; 0052 end 0053 0054 %% options 0055 tol = mpopt(2); 0056 max_it = mpopt(3); 0057 verbose = mpopt(31); 0058 0059 %% initialize 0060 converged = 0; 0061 i = 0; 0062 V = V0; 0063 Va = angle(V); 0064 Vm = abs(V); 0065 0066 %% set up indexing for updating V 0067 npv = length(pv); 0068 npq = length(pq); 0069 j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses 0070 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses 0071 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses 0072 0073 %% evaluate F(x0) 0074 mis = V .* conj(Ybus * V) - Sbus; 0075 F = [ real(mis([pv; pq])); 0076 imag(mis(pq)) ]; 0077 0078 %% check tolerance 0079 normF = norm(F, inf); 0080 if verbose > 0 0081 fprintf('(Newton)\n'); 0082 end 0083 if verbose > 1 0084 fprintf('\n it max P & Q mismatch (p.u.)'); 0085 fprintf('\n---- ---------------------------'); 0086 fprintf('\n%3d %10.3e', i, normF); 0087 end 0088 if normF < tol 0089 converged = 1; 0090 if verbose > 1 0091 fprintf('\nConverged!\n'); 0092 end 0093 end 0094 0095 %% do Newton iterations 0096 while (~converged && i < max_it) 0097 %% update iteration counter 0098 i = i + 1; 0099 0100 %% evaluate Jacobian 0101 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); 0102 0103 j11 = real(dSbus_dVa([pv; pq], [pv; pq])); 0104 j12 = real(dSbus_dVm([pv; pq], pq)); 0105 j21 = imag(dSbus_dVa(pq, [pv; pq])); 0106 j22 = imag(dSbus_dVm(pq, pq)); 0107 0108 J = [ j11 j12; 0109 j21 j22; ]; 0110 0111 %% compute update step 0112 dx = -(J \ F); 0113 0114 %% update voltage 0115 if npv 0116 Va(pv) = Va(pv) + dx(j1:j2); 0117 end 0118 if npq 0119 Va(pq) = Va(pq) + dx(j3:j4); 0120 Vm(pq) = Vm(pq) + dx(j5:j6); 0121 end 0122 V = Vm .* exp(1j * Va); 0123 Vm = abs(V); %% update Vm and Va again in case 0124 Va = angle(V); %% we wrapped around with a negative Vm 0125 0126 %% evalute F(x) 0127 mis = V .* conj(Ybus * V) - Sbus; 0128 F = [ real(mis(pv)); 0129 real(mis(pq)); 0130 imag(mis(pq)) ]; 0131 0132 %% check for convergence 0133 normF = norm(F, inf); 0134 if verbose > 1 0135 fprintf('\n%3d %10.3e', i, normF); 0136 end 0137 if normF < tol 0138 converged = 1; 0139 if verbose 0140 fprintf('\nNewton''s method power flow converged in %d iterations.\n', i); 0141 end 0142 end 0143 end 0144 0145 if verbose 0146 if ~converged 0147 fprintf('\nNewton''s method power did not converge in %d iterations.\n', i); 0148 end 0149 end