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 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] = 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 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 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 % Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) 0022 % by Ray Zimmerman, PSERC Cornell 0023 % 0024 % This file is part of MATPOWER. 0025 % Covered by the 3-clause BSD License (see LICENSE file for details). 0026 % See http://www.pserc.cornell.edu/matpower/ for more info. 0027 0028 %% default arguments 0029 if nargin < 7 0030 mpopt = mpoption; 0031 end 0032 0033 %% options 0034 tol = mpopt.pf.tol; 0035 max_it = mpopt.pf.nr.max_it; 0036 0037 %% initialize 0038 converged = 0; 0039 i = 0; 0040 V = V0; 0041 Va = angle(V); 0042 Vm = abs(V); 0043 0044 %% set up indexing for updating V 0045 npv = length(pv); 0046 npq = length(pq); 0047 j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses 0048 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses 0049 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses 0050 0051 %% evaluate F(x0) 0052 mis = V .* conj(Ybus * V) - Sbus(Vm); 0053 F = [ real(mis([pv; pq])); 0054 imag(mis(pq)) ]; 0055 0056 %% check tolerance 0057 normF = norm(F, inf); 0058 if mpopt.verbose > 1 0059 fprintf('\n it max P & Q mismatch (p.u.)'); 0060 fprintf('\n---- ---------------------------'); 0061 fprintf('\n%3d %10.3e', i, normF); 0062 end 0063 if normF < tol 0064 converged = 1; 0065 if mpopt.verbose > 1 0066 fprintf('\nConverged!\n'); 0067 end 0068 end 0069 0070 %% do Newton iterations 0071 while (~converged && i < max_it) 0072 %% update iteration counter 0073 i = i + 1; 0074 0075 %% evaluate Jacobian 0076 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); 0077 [dummy, neg_dSd_dVm] = Sbus(Vm); 0078 dSbus_dVm = dSbus_dVm - neg_dSd_dVm; 0079 0080 j11 = real(dSbus_dVa([pv; pq], [pv; pq])); 0081 j12 = real(dSbus_dVm([pv; pq], pq)); 0082 j21 = imag(dSbus_dVa(pq, [pv; pq])); 0083 j22 = imag(dSbus_dVm(pq, pq)); 0084 0085 J = [ j11 j12; 0086 j21 j22; ]; 0087 0088 %% compute update step 0089 dx = -(J \ F); 0090 0091 %% update voltage 0092 if npv 0093 Va(pv) = Va(pv) + dx(j1:j2); 0094 end 0095 if npq 0096 Va(pq) = Va(pq) + dx(j3:j4); 0097 Vm(pq) = Vm(pq) + dx(j5:j6); 0098 end 0099 V = Vm .* exp(1j * Va); 0100 Vm = abs(V); %% update Vm and Va again in case 0101 Va = angle(V); %% we wrapped around with a negative Vm 0102 0103 %% evalute F(x) 0104 mis = V .* conj(Ybus * V) - Sbus(Vm); 0105 F = [ real(mis(pv)); 0106 real(mis(pq)); 0107 imag(mis(pq)) ]; 0108 0109 %% check for convergence 0110 normF = norm(F, inf); 0111 if mpopt.verbose > 1 0112 fprintf('\n%3d %10.3e', i, normF); 0113 end 0114 if normF < tol 0115 converged = 1; 0116 if mpopt.verbose 0117 fprintf('\nNewton''s method power flow converged in %d iterations.\n', i); 0118 end 0119 end 0120 end 0121 0122 if mpopt.verbose 0123 if ~converged 0124 fprintf('\nNewton''s method power flow did not converge in %d iterations.\n', i); 0125 end 0126 end