NEWTONPF_I_CART Solves power flow using full Newton's method (current/cartesian)
[V, CONVERGED, I] = NEWTONPF_I_CART(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
Solves for bus voltages using a full Newton-Raphson method, using nodal
current balance equations and cartesian 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_POLAR.0001 function [V, converged, i] = newtonpf_I_cart(Ybus, Sbus, V0, ref, pv, pq, mpopt) 0002 %NEWTONPF_I_CART Solves power flow using full Newton's method (current/cartesian) 0003 % [V, CONVERGED, I] = NEWTONPF_I_CART(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 cartesian 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_POLAR. 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 Vm = abs(V); 0054 Vmpv = Vm(pv); 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 - Q of pv buses 0061 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - Vr of pq buses 0062 j5 = j4 + 1; j6 = j4 + npv; %% j5:j6 - Vr of pv buses 0063 j7 = j6 + 1; j8 = j6 + npq; %% j7:j8 - Vi of pq buses 0064 j9 = j8 + 1; j10= j8 + npv; %% j9:j10- Vi of pv buses 0065 0066 %% evaluate F(x0) 0067 Sb = Sbus(Vm); 0068 Sb(pv) = real(Sb(pv)) + 1j * imag(V(pv) .* conj(Ybus(pv, :) * V)); 0069 mis = Ybus * V - conj(Sb ./ V); 0070 F = [ real(mis([pv; pq])); 0071 imag(mis([pv; pq])); 0072 V(pv) .* conj(V(pv)) - Vmpv.^2 ]; 0073 0074 %% check tolerance 0075 normF = norm(F, inf); 0076 if mpopt.verbose > 1 0077 fprintf('\n it max Ir & Ii mismatch (p.u.)'); 0078 fprintf('\n---- ---------------------------'); 0079 fprintf('\n%3d %10.3e', i, normF); 0080 end 0081 if normF < tol 0082 converged = 1; 0083 if mpopt.verbose > 1 0084 fprintf('\nConverged!\n'); 0085 end 0086 end 0087 0088 %% attempt to pick fastest linear solver, if not specified 0089 if isempty(lin_solver) 0090 nx = length(F); 0091 if nx <= 10 || have_feature('octave') 0092 lin_solver = '\'; %% default \ operator 0093 else %% MATLAB and nx > 10 or Octave and nx > 2000 0094 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0095 end 0096 end 0097 0098 %% do Newton iterations 0099 while (~converged && i < max_it) 0100 %% update iteration counter 0101 i = i + 1; 0102 0103 %% evaluate Jacobian 0104 dImis_dQ = sparse(pv, pv, 1j./conj(V(pv)), n, n); 0105 dV2_dVr = sparse(1:npv, npq+(1:npv), 2*real(V(pv)), npv, npv+npq); 0106 dV2_dVi = sparse(1:npv, npq+(1:npv), 2*imag(V(pv)), npv, npv+npq); 0107 [dImis_dVr, dImis_dVi] = dImis_dV(Sb, Ybus, V, 1); 0108 0109 %% handling of derivatives for voltage dependent loads 0110 %% (not yet implemented) goes here 0111 0112 j11 = real(dImis_dQ([pv; pq], pv)); 0113 j12 = real(dImis_dVr([pv; pq], [pq; pv])); 0114 j13 = real(dImis_dVi([pv; pq], [pq; pv])); 0115 j21 = imag(dImis_dQ([pv; pq], pv)); 0116 j22 = imag(dImis_dVr([pv; pq], [pq; pv])); 0117 j23 = imag(dImis_dVi([pv; pq], [pq; pv])); 0118 j31 = sparse(npv, npv); 0119 j32 = dV2_dVr; 0120 j33 = dV2_dVi; 0121 0122 J = [ j11 j12 j13; 0123 j21 j22 j23; 0124 j31 j32 j33; ]; 0125 0126 %% compute update step 0127 dx = mplinsolve(J, -F, lin_solver); 0128 0129 %% update voltage 0130 if npv 0131 V(pv) = V(pv) + dx(j5:j6) + 1j * dx(j9:j10); 0132 Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j1:j2)); 0133 end 0134 if npq 0135 V(pq) = V(pq) + dx(j3:j4) + 1j * dx(j7:j8); 0136 end 0137 0138 %% evalute F(x) 0139 mis = Ybus * V - conj(Sb ./ V); 0140 F = [ real(mis([pv; pq])); 0141 imag(mis([pv; pq])); 0142 V(pv) .* conj(V(pv)) - Vmpv.^2 ]; 0143 0144 %% check for convergence 0145 normF = norm(F, inf); 0146 if mpopt.verbose > 1 0147 fprintf('\n%3d %10.3e', i, normF); 0148 end 0149 if normF < tol 0150 converged = 1; 0151 if mpopt.verbose 0152 fprintf('\nNewton''s method power flow (current balance, cartesian) converged in %d iterations.\n', i); 0153 end 0154 end 0155 end 0156 0157 if mpopt.verbose 0158 if ~converged 0159 fprintf('\nNewton''s method power flow (current balance, cartesian) did not converge in %d iterations.\n', i); 0160 end 0161 end