NEWTONPF_I_POLAR Solves power flow using full Newton's method (power/cartesian)
[V, CONVERGED, I] = NEWTONPF_I_POLAR(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
Solves for bus voltages using a full Newton-Raphson method, using nodal
current 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, NEWTONPF_S_CART, NEWTONPF_I_CART.0001 function [V, converged, i] = newtonpf_I_polar(Ybus, Sbus, V0, ref, pv, pq, mpopt) 0002 %NEWTONPF_I_POLAR Solves power flow using full Newton's method (power/cartesian) 0003 % [V, CONVERGED, I] = NEWTONPF_I_POLAR(YBUS, SBUS, V0, REF, PV, PQ, MPOPT) 0004 % 0005 % Solves for bus voltages using a full Newton-Raphson method, using nodal 0006 % current 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, NEWTONPF_S_CART, NEWTONPF_I_CART. 0029 0030 % MATPOWER 0031 % Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC) 0032 % by Ray Zimmerman, PSERC Cornell 0033 % and Baljinnyam Sereeter, Delft University of Technology 0034 % 0035 % This file is part of MATPOWER. 0036 % Covered by the 3-clause BSD License (see LICENSE file for details). 0037 % See https://matpower.org for more info. 0038 0039 %% default arguments 0040 if nargin < 7 0041 mpopt = mpoption; 0042 end 0043 0044 %% options 0045 tol = mpopt.pf.tol; 0046 max_it = mpopt.pf.nr.max_it; 0047 lin_solver = mpopt.pf.nr.lin_solver; 0048 0049 %% initialize 0050 converged = 0; 0051 i = 0; 0052 V = V0; 0053 Va = angle(V); 0054 Vm = abs(V); 0055 n = length(V0); 0056 0057 %% set up indexing for updating V 0058 npv = length(pv); 0059 npq = length(pq); 0060 j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses 0061 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses 0062 j5 = j4 + 1; j6 = j4 + npv; %% j5:j6 - Q of pv buses 0063 j7 = j6 + 1; j8 = j6 + npq; %% j7:j8 - V mag of pq buses 0064 0065 %% evaluate F(x0) 0066 Sb = Sbus(Vm); 0067 Sb(pv) = real(Sb(pv)) + 1j * imag(V(pv) .* conj(Ybus(pv, :) * V)); 0068 mis = Ybus * V - conj(Sb ./ V); 0069 F = [ real(mis([pv; pq])); 0070 imag(mis([pv; pq])) ]; 0071 0072 %% check tolerance 0073 normF = norm(F, inf); 0074 if mpopt.verbose > 1 0075 fprintf('\n it max Ir & Ii mismatch (p.u.)'); 0076 fprintf('\n---- ---------------------------'); 0077 fprintf('\n%3d %10.3e', i, normF); 0078 end 0079 if normF < tol 0080 converged = 1; 0081 if mpopt.verbose > 1 0082 fprintf('\nConverged!\n'); 0083 end 0084 end 0085 0086 %% attempt to pick fastest linear solver, if not specified 0087 if isempty(lin_solver) 0088 nx = length(F); 0089 if nx <= 10 || have_fcn('octave') 0090 lin_solver = '\'; %% default \ operator 0091 else %% MATLAB and nx > 10 or Octave and nx > 2000 0092 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0093 end 0094 end 0095 0096 %% do Newton iterations 0097 while (~converged && i < max_it) 0098 %% update iteration counter 0099 i = i + 1; 0100 0101 %% evaluate Jacobian 0102 dImis_dQ = sparse(1:n, 1:n, 1j./conj(V), n, n); 0103 [dImis_dVa, dImis_dVm] = dImis_dV(Sb, Ybus, V); 0104 dImis_dVm(:, pv) = dImis_dQ(:, pv); 0105 0106 %% handling of derivatives for voltage dependent loads 0107 %% (not yet implemented) goes here 0108 0109 j11 = real(dImis_dVa([pv; pq], [pv; pq])); 0110 j12 = real(dImis_dVm([pv; pq], [pv; pq])); 0111 j21 = imag(dImis_dVa([pv; pq], [pv; pq])); 0112 j22 = imag(dImis_dVm([pv; pq], [pv; pq])); 0113 0114 J = [ j11 j12; 0115 j21 j22; ]; 0116 0117 %% compute update step 0118 dx = mplinsolve(J, -F, lin_solver); 0119 0120 %% update voltage 0121 if npv 0122 Va(pv) = Va(pv) + dx(j1:j2); 0123 Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j5:j6)); 0124 end 0125 if npq 0126 Va(pq) = Va(pq) + dx(j3:j4); 0127 Vm(pq) = Vm(pq) + dx(j7:j8); 0128 end 0129 V = Vm .* exp(1j * Va); 0130 Vm = abs(V); %% update Vm and Va again in case 0131 Va = angle(V); %% we wrapped around with a negative Vm 0132 0133 %% evalute F(x) 0134 mis = Ybus * V - conj(Sb ./ V); 0135 F = [ real(mis([pv; pq])); 0136 imag(mis([pv; pq])) ]; 0137 0138 %% check for convergence 0139 normF = norm(F, inf); 0140 if mpopt.verbose > 1 0141 fprintf('\n%3d %10.3e', i, normF); 0142 end 0143 if normF < tol 0144 converged = 1; 0145 if mpopt.verbose 0146 fprintf('\nNewton''s method power flow (current balance, polar) converged in %d iterations.\n', i); 0147 end 0148 end 0149 end 0150 0151 if mpopt.verbose 0152 if ~converged 0153 fprintf('\nNewton''s method power flow (current balance, polar) did not converge in %d iterations.\n', i); 0154 end 0155 end