NEWTONPF_S_HYBRID Solves power flow using full Newton's method (power/hybrid)
[V, CONVERGED, I] = NEWTONPF_S_HYBRID(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
Solves for bus voltages using a full Newton-Raphson method, using nodal
power balance equations and a hybrid representation of voltages, where
a polar update is computed using a cartesian Jacobian, 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_hybrid(Ybus, Sbus, V0, ref, pv, pq, mpopt) 0002 %NEWTONPF_S_HYBRID Solves power flow using full Newton's method (power/hybrid) 0003 % [V, CONVERGED, I] = NEWTONPF_S_HYBRID(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 a hybrid representation of voltages, where 0007 % a polar update is computed using a cartesian Jacobian, given the 0008 % following inputs: 0009 % YBUS - full system admittance matrix (for all buses) 0010 % SBUS - handle to function that returns the complex bus power 0011 % injection vector (for all buses), given the bus voltage 0012 % magnitude vector (for all buses) 0013 % V0 - initial vector of complex bus voltages 0014 % REF - bus index of reference bus (voltage ang reference & gen slack) 0015 % PV - vector of bus indices for PV buses 0016 % PQ - vector of bus indices for PQ buses 0017 % MPOPT - (optional) MATPOWER option struct, used to set the 0018 % termination tolerance, maximum number of iterations, and 0019 % output options (see MPOPTION for details). 0020 % 0021 % The bus voltage vector contains the set point for generator 0022 % (including ref bus) buses, and the reference angle of the swing 0023 % bus, as well as an initial guess for remaining magnitudes and 0024 % angles. 0025 % 0026 % Returns the final complex voltages, a flag which indicates whether it 0027 % converged or not, and the number of iterations performed. 0028 % 0029 % See also RUNPF, NEWTONPF, NEWTONPF_I_POLAR, NEWTONPF_I_CART. 0030 0031 % MATPOWER 0032 % Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC) 0033 % by Ray Zimmerman, PSERC Cornell 0034 % and Baljinnyam Sereeter, Delft University of Technology 0035 % 0036 % This file is part of MATPOWER. 0037 % Covered by the 3-clause BSD License (see LICENSE file for details). 0038 % See https://matpower.org for more info. 0039 0040 %% default arguments 0041 if nargin < 7 0042 mpopt = mpoption; 0043 end 0044 0045 %% options 0046 tol = mpopt.pf.tol; 0047 max_it = mpopt.pf.nr.max_it; 0048 lin_solver = mpopt.pf.nr.lin_solver; 0049 0050 %% initialize 0051 converged = 0; 0052 i = 0; 0053 V = V0; 0054 Va = angle(V); 0055 Vm = abs(V); 0056 0057 %% set up indexing for updating V 0058 npv = length(pv); 0059 npq = length(pq); 0060 j1 = 1; j2 = npq; %% j1:j2 - Vr of pq buses 0061 j3 = j2 + 1; j4 = j2 + npv; %% j3:j4 - Vi of pv buses 0062 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - Vi of pq 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 0069 %% check tolerance 0070 normF = norm(F, inf); 0071 if mpopt.verbose > 1 0072 fprintf('\n it max P & Q mismatch (p.u.)'); 0073 fprintf('\n---- ---------------------------'); 0074 fprintf('\n%3d %10.3e', i, normF); 0075 end 0076 if normF < tol 0077 converged = 1; 0078 if mpopt.verbose > 1 0079 fprintf('\nConverged!\n'); 0080 end 0081 end 0082 0083 %% attempt to pick fastest linear solver, if not specified 0084 if isempty(lin_solver) 0085 nx = length(F); 0086 if nx <= 10 || have_feature('octave') 0087 lin_solver = '\'; %% default \ operator 0088 else %% MATLAB and nx > 10 or Octave and nx > 2000 0089 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0090 end 0091 end 0092 0093 %% do Newton iterations 0094 while (~converged && i < max_it) 0095 %% update iteration counter 0096 i = i + 1; 0097 0098 %% evaluate Jacobian 0099 [dSbus_dVr, dSbus_dVi] = dSbus_dV(Ybus, V, 1); 0100 dSbus_dVi(:, pv) = ... 0101 dSbus_dVi(:, pv) * sparse(1:npv, 1:npv, real(V(pv)), npv, npv) - ... 0102 dSbus_dVr(:, pv) * sparse(1:npv, 1:npv, imag(V(pv)), npv, npv); 0103 0104 %% handling of derivatives for voltage dependent loads 0105 %% (not yet implemented) goes here 0106 0107 j11 = real(dSbus_dVr([pq; pv], pq)); 0108 j12 = real(dSbus_dVi([pq; pv], [pv; pq])); 0109 j21 = imag(dSbus_dVr(pq, pq)); 0110 j22 = imag(dSbus_dVi(pq, [pv; pq])); 0111 0112 J = [ j11 j12; 0113 j21 j22; ]; 0114 0115 %% compute update step 0116 dx = mplinsolve(J, -F, lin_solver); 0117 0118 %% update voltage 0119 if npv 0120 Va(pv) = Va(pv) + dx(j3:j4); 0121 end 0122 if npq 0123 Vm(pq) = Vm(pq) + (real(V(pq))./Vm(pq)) .* dx(j1:j2) ... 0124 + (imag(V(pq))./Vm(pq)) .* dx(j5:j6); 0125 Va(pq) = Va(pq) + (real(V(pq))./(Vm(pq).^2)) .* dx(j5:j6) ... 0126 - (imag(V(pq))./(Vm(pq).^2)) .* dx(j1:j2); 0127 end 0128 V = Vm .* exp(1j * Va); 0129 Vm = abs(V); %% update Vm and Va again in case 0130 Va = angle(V); %% we wrapped around with a negative Vm 0131 0132 %% evalute F(x) 0133 mis = V .* conj(Ybus * V) - Sbus(Vm); 0134 F = [ real(mis([pq; pv])); 0135 imag(mis(pq)) ]; 0136 0137 %% check for convergence 0138 normF = norm(F, inf); 0139 if mpopt.verbose > 1 0140 fprintf('\n%3d %10.3e', i, normF); 0141 end 0142 if normF < tol 0143 converged = 1; 0144 if mpopt.verbose 0145 fprintf('\nNewton''s method power flow (power balance, hybrid) converged in %d iterations.\n', i); 0146 end 0147 end 0148 end 0149 0150 if mpopt.verbose 0151 if ~converged 0152 fprintf('\nNewton''s method power flow (power balance, hybrid) did not converge in %d iterations.\n', i); 0153 end 0154 end