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