NEWTONPF Solves power flow using full Newton's method (power/polar)
[V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
Solves for bus voltages using a full Newton-Raphson method, using nodal
power balance equations and polar coordinate representation of
voltages, given the following inputs:
YBUS - full system admittance matrix (for all buses)
SBUS - handle to function that returns the complex bus power
injection vector (for all buses), given the bus voltage
magnitude vector (for all buses)
V0 - initial vector of complex bus voltages
REF - bus index of reference bus (voltage ang reference & gen slack)
PV - vector of bus indices for PV buses
PQ - vector of bus indices for PQ buses
MPOPT - (optional) MATPOWER option struct, used to set the
termination tolerance, maximum number of iterations, and
output options (see MPOPTION for details).
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.
Returns the final complex voltages, a flag which indicates whether it
converged or not, and the number of iterations performed.
See also RUNPF, NEWTONPF_S_CART, NEWTONPF_I_POLAR, NEWTONPF_I_CART.0001 function [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) 0002 %NEWTONPF Solves power flow using full Newton's method (power/polar) 0003 % [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT) 0004 % 0005 % Solves for bus voltages using a full Newton-Raphson method, using nodal 0006 % power balance equations and polar coordinate representation of 0007 % voltages, given the following inputs: 0008 % YBUS - full system admittance matrix (for all buses) 0009 % SBUS - handle to function that returns the complex bus power 0010 % injection vector (for all buses), given the bus voltage 0011 % magnitude vector (for all buses) 0012 % V0 - initial vector of complex bus voltages 0013 % REF - bus index of reference bus (voltage ang reference & gen slack) 0014 % PV - vector of bus indices for PV buses 0015 % PQ - vector of bus indices for PQ buses 0016 % MPOPT - (optional) MATPOWER option struct, used to set the 0017 % termination tolerance, maximum number of iterations, and 0018 % output options (see MPOPTION for details). 0019 % 0020 % The bus voltage vector contains the set point for generator 0021 % (including ref bus) buses, and the reference angle of the swing 0022 % bus, as well as an initial guess for remaining magnitudes and 0023 % angles. 0024 % 0025 % Returns the final complex voltages, a flag which indicates whether it 0026 % converged or not, and the number of iterations performed. 0027 % 0028 % See also RUNPF, NEWTONPF_S_CART, NEWTONPF_I_POLAR, NEWTONPF_I_CART. 0029 0030 % MATPOWER 0031 % Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC) 0032 % by Ray Zimmerman, PSERC Cornell 0033 % 0034 % This file is part of MATPOWER. 0035 % Covered by the 3-clause BSD License (see LICENSE file for details). 0036 % See https://matpower.org for more info. 0037 0038 %% default arguments 0039 if nargin < 7 0040 mpopt = mpoption; 0041 end 0042 0043 %% options 0044 tol = mpopt.pf.tol; 0045 max_it = mpopt.pf.nr.max_it; 0046 lin_solver = mpopt.pf.nr.lin_solver; 0047 0048 %% initialize 0049 converged = 0; 0050 i = 0; 0051 V = V0; 0052 Va = angle(V); 0053 Vm = abs(V); 0054 0055 %% set up indexing for updating V 0056 npv = length(pv); 0057 npq = length(pq); 0058 j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses 0059 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses 0060 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses 0061 0062 %% evaluate F(x0) 0063 mis = V .* conj(Ybus * V) - Sbus(Vm); 0064 F = [ real(mis([pv; pq])); 0065 imag(mis(pq)) ]; 0066 0067 %% check tolerance 0068 normF = norm(F, inf); 0069 if mpopt.verbose > 1 0070 fprintf('\n it max P & Q mismatch (p.u.)'); 0071 fprintf('\n---- ---------------------------'); 0072 fprintf('\n%3d %10.3e', i, normF); 0073 end 0074 if normF < tol 0075 converged = 1; 0076 if mpopt.verbose > 1 0077 fprintf('\nConverged!\n'); 0078 end 0079 end 0080 0081 %% attempt to pick fastest linear solver, if not specified 0082 if isempty(lin_solver) 0083 nx = length(F); 0084 if nx <= 10 || have_fcn('octave') 0085 lin_solver = '\'; %% default \ operator 0086 else %% MATLAB and nx > 10 or Octave and nx > 2000 0087 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0088 end 0089 end 0090 0091 %% do Newton iterations 0092 while (~converged && i < max_it) 0093 %% update iteration counter 0094 i = i + 1; 0095 0096 %% evaluate Jacobian 0097 [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V); 0098 [dummy, neg_dSd_dVm] = Sbus(Vm); 0099 dSbus_dVm = dSbus_dVm - neg_dSd_dVm; 0100 0101 j11 = real(dSbus_dVa([pv; pq], [pv; pq])); 0102 j12 = real(dSbus_dVm([pv; pq], pq)); 0103 j21 = imag(dSbus_dVa(pq, [pv; pq])); 0104 j22 = imag(dSbus_dVm(pq, pq)); 0105 0106 J = [ j11 j12; 0107 j21 j22; ]; 0108 0109 %% compute update step 0110 dx = mplinsolve(J, -F, lin_solver); 0111 0112 %% update voltage 0113 if npv 0114 Va(pv) = Va(pv) + dx(j1:j2); 0115 end 0116 if npq 0117 Va(pq) = Va(pq) + dx(j3:j4); 0118 Vm(pq) = Vm(pq) + dx(j5:j6); 0119 end 0120 V = Vm .* exp(1j * Va); 0121 Vm = abs(V); %% update Vm and Va again in case 0122 Va = angle(V); %% we wrapped around with a negative Vm 0123 0124 %% evalute F(x) 0125 mis = V .* conj(Ybus * V) - Sbus(Vm); 0126 F = [ real(mis([pv; pq])); 0127 imag(mis(pq)) ]; 0128 0129 %% check for convergence 0130 normF = norm(F, inf); 0131 if mpopt.verbose > 1 0132 fprintf('\n%3d %10.3e', i, normF); 0133 end 0134 if normF < tol 0135 converged = 1; 0136 if mpopt.verbose 0137 fprintf('\nNewton''s method power flow (power balance, polar) converged in %d iterations.\n', i); 0138 end 0139 end 0140 end 0141 0142 if mpopt.verbose 0143 if ~converged 0144 fprintf('\nNewton''s method power flow (power balance, polar) did not converge in %d iterations.\n', i); 0145 end 0146 end