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 Va = angle(V); 0054 Vm = abs(V); 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 - Vi of pv buses 0061 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - Vi of pq buses 0062 0063 %% evaluate F(x0) 0064 mis = V .* conj(Ybus * V) - Sbus(Vm); 0065 F = [ real(mis([pq; pv])); 0066 imag(mis(pq)) ]; 0067 0068 %% check tolerance 0069 normF = norm(F, inf); 0070 if mpopt.verbose > 1 0071 fprintf('\n it max P & Q mismatch (p.u.)'); 0072 fprintf('\n---- ---------------------------'); 0073 fprintf('\n%3d %10.3e', i, normF); 0074 end 0075 if normF < tol 0076 converged = 1; 0077 if mpopt.verbose > 1 0078 fprintf('\nConverged!\n'); 0079 end 0080 end 0081 0082 %% attempt to pick fastest linear solver, if not specified 0083 if isempty(lin_solver) 0084 nx = length(F); 0085 if nx <= 10 || have_fcn('octave') 0086 lin_solver = '\'; %% default \ operator 0087 else %% MATLAB and nx > 10 or Octave and nx > 2000 0088 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0089 end 0090 end 0091 0092 %% do Newton iterations 0093 while (~converged && i < max_it) 0094 %% update iteration counter 0095 i = i + 1; 0096 0097 %% evaluate Jacobian 0098 [dSbus_dVr, dSbus_dVi] = dSbus_dV(Ybus, V, 1); 0099 dSbus_dVi(:, pv) = dSbus_dVi(:, pv) - ... 0100 dSbus_dVr(:, pv) * sparse(1:npv, 1:npv, imag(V(pv))./real(V(pv)), npv, npv); 0101 0102 %% handling of derivatives for voltage dependent loads 0103 %% (not yet implemented) goes here 0104 0105 j11 = real(dSbus_dVr([pq; pv], pq)); 0106 j12 = real(dSbus_dVi([pq; pv], [pv; pq])); 0107 j21 = imag(dSbus_dVr(pq, pq)); 0108 j22 = imag(dSbus_dVi(pq, [pv; pq])); 0109 0110 J = [ j11 j12; 0111 j21 j22; ]; 0112 0113 %% compute update step 0114 dx = mplinsolve(J, -F, lin_solver); 0115 0116 %% update voltage 0117 if npv 0118 Va(pv) = Va(pv) + dx(j3:j4) ./ real(V(pv)); 0119 end 0120 if npq 0121 Vm(pq) = Vm(pq) + (real(V(pq))./Vm(pq)) .* dx(j1:j2) ... 0122 + (imag(V(pq))./Vm(pq)) .* dx(j5:j6); 0123 Va(pq) = Va(pq) + (real(V(pq))./(Vm(pq).^2)) .* dx(j5:j6) ... 0124 - (imag(V(pq))./(Vm(pq).^2)) .* dx(j1:j2); 0125 end 0126 V = Vm .* exp(1j * Va); 0127 Vm = abs(V); %% update Vm and Va again in case 0128 Va = angle(V); %% we wrapped around with a negative Vm 0129 0130 %% evalute F(x) 0131 mis = V .* conj(Ybus * V) - Sbus(Vm); 0132 F = [ real(mis([pq; pv])); 0133 imag(mis(pq)) ]; 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