NEWTONPF_I_POLAR Solves power flow using full Newton's method (current/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 (current/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_feature('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(pv, pv, 1j./conj(V(pv)), 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 % %% evaluate Jacobian 0118 % dImis_dQ = sparse(pv, pv, 1j./conj(V(pv)), n, n); 0119 % [dImis_dVa, dImis_dVm] = dImis_dV(Sb, Ybus, V); 0120 % 0121 % %% handling of derivatives for voltage dependent loads 0122 % %% (not yet implemented) goes here 0123 % 0124 % j11 = real(dImis_dVa([pv; pq], [pv; pq])); 0125 % j12 = real(dImis_dQ([pv; pq], pv)); 0126 % j13 = real(dImis_dVm([pv; pq], pq)); 0127 % j21 = imag(dImis_dVa([pv; pq], [pv; pq])); 0128 % j22 = imag(dImis_dQ([pv; pq], pv)); 0129 % j23 = imag(dImis_dVm([pv; pq], pq)); 0130 % 0131 % J = [ j11 j12 j13; 0132 % j21 j22 j23; ]; 0133 0134 %% compute update step 0135 dx = mplinsolve(J, -F, lin_solver); 0136 0137 %% update voltage 0138 if npv 0139 Va(pv) = Va(pv) + dx(j1:j2); 0140 Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j5:j6)); 0141 end 0142 if npq 0143 Va(pq) = Va(pq) + dx(j3:j4); 0144 Vm(pq) = Vm(pq) + dx(j7:j8); 0145 end 0146 V = Vm .* exp(1j * Va); 0147 Vm = abs(V); %% update Vm and Va again in case 0148 Va = angle(V); %% we wrapped around with a negative Vm 0149 0150 %% evalute F(x) 0151 mis = Ybus * V - conj(Sb ./ V); 0152 F = [ real(mis([pv; pq])); 0153 imag(mis([pv; pq])) ]; 0154 0155 %% check for convergence 0156 normF = norm(F, inf); 0157 if mpopt.verbose > 1 0158 fprintf('\n%3d %10.3e', i, normF); 0159 end 0160 if normF < tol 0161 converged = 1; 0162 if mpopt.verbose 0163 fprintf('\nNewton''s method power flow (current balance, polar) converged in %d iterations.\n', i); 0164 end 0165 end 0166 end 0167 0168 if mpopt.verbose 0169 if ~converged 0170 fprintf('\nNewton''s method power flow (current balance, polar) did not converge in %d iterations.\n', i); 0171 end 0172 end