NEWTONPF_I_CART Solves power flow using full Newton's method (power/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 (power/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 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 - 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 - Vi of pv buses 0063 j7 = j6 + 1; j8 = j6 + npq; %% j7:j9 - Vi 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_dVr, dImis_dVi] = dImis_dV(Sb, Ybus, V, 1); 0104 dImis_dVi(:, pv) = dImis_dVi(:, pv) - ... 0105 dImis_dVr(:, pv) * sparse(1:npv, 1:npv, imag(V(pv))./real(V(pv)), npv, npv); 0106 dImis_dVr(:, pv) = dImis_dQ(:, pv); 0107 0108 %% handling of derivatives for voltage dependent loads 0109 %% (not yet implemented) goes here 0110 0111 j11 = real(dImis_dVr([pv; pq], [pv; pq])); 0112 j12 = real(dImis_dVi([pv; pq], [pv; pq])); 0113 j21 = imag(dImis_dVr([pv; pq], [pv; pq])); 0114 j22 = imag(dImis_dVi([pv; pq], [pv; pq])); 0115 0116 J = [ j11 j12; 0117 j21 j22; ]; 0118 0119 %% compute update step 0120 dx = mplinsolve(J, -F, lin_solver); 0121 0122 %% update voltage 0123 if npv 0124 Va(pv) = Va(pv) + dx(j5:j6) ./ real(V(pv)); 0125 Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j1:j2)); 0126 end 0127 if npq 0128 Vm(pq) = Vm(pq) + (real(V(pq))./Vm(pq)) .* dx(j3:j4) ... 0129 + (imag(V(pq))./Vm(pq)) .* dx(j7:j8); 0130 Va(pq) = Va(pq) + (real(V(pq))./(Vm(pq).^2)) .* dx(j7:j8) ... 0131 - (imag(V(pq))./(Vm(pq).^2)) .* dx(j3:j4); 0132 end 0133 V = Vm .* exp(1j * Va); 0134 Vm = abs(V); %% update Vm and Va again in case 0135 Va = angle(V); %% we wrapped around with a negative Vm 0136 0137 %% evalute F(x) 0138 mis = Ybus * V - conj(Sb ./ V); 0139 F = [ real(mis([pv; pq])); 0140 imag(mis([pv; pq])) ]; 0141 0142 %% check for convergence 0143 normF = norm(F, inf); 0144 if mpopt.verbose > 1 0145 fprintf('\n%3d %10.3e', i, normF); 0146 end 0147 if normF < tol 0148 converged = 1; 0149 if mpopt.verbose 0150 fprintf('\nNewton''s method power flow (current balance, cartesian) converged in %d iterations.\n', i); 0151 end 0152 end 0153 end 0154 0155 if mpopt.verbose 0156 if ~converged 0157 fprintf('\nNewton''s method power flow (current balance, cartesian) did not converge in %d iterations.\n', i); 0158 end 0159 end