NEWTONPF_I_HYBRID Solves power flow using full Newton's method (current/hybrid)
[V, CONVERGED, I] = NEWTONPF_I_HYBRID(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
Solves for bus voltages using a full Newton-Raphson method, using nodal
current 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_S_CART, NEWTONPF_I_POLAR.0001 function [V, converged, i] = newtonpf_I_hybrid(Ybus, Sbus, V0, ref, pv, pq, mpopt) 0002 %NEWTONPF_I_HYBRID Solves power flow using full Newton's method (current/hybrid) 0003 % [V, CONVERGED, I] = NEWTONPF_I_HYBRID(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 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_S_CART, NEWTONPF_I_POLAR. 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 n = length(V0); 0057 0058 %% set up indexing for updating V 0059 npv = length(pv); 0060 npq = length(pq); 0061 j1 = 1; j2 = npv; %% j1:j2 - Q of pv buses 0062 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - Vr of pq buses 0063 j5 = j4 + 1; j6 = j4 + npv; %% j5:j6 - Vi of pv buses 0064 j7 = j6 + 1; j8 = j6 + npq; %% j7:j9 - Vi of pq 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 0073 %% check tolerance 0074 normF = norm(F, inf); 0075 if mpopt.verbose > 1 0076 fprintf('\n it max Ir & Ii mismatch (p.u.)'); 0077 fprintf('\n---- ---------------------------'); 0078 fprintf('\n%3d %10.3e', i, normF); 0079 end 0080 if normF < tol 0081 converged = 1; 0082 if mpopt.verbose > 1 0083 fprintf('\nConverged!\n'); 0084 end 0085 end 0086 0087 %% attempt to pick fastest linear solver, if not specified 0088 if isempty(lin_solver) 0089 nx = length(F); 0090 if nx <= 10 || have_feature('octave') 0091 lin_solver = '\'; %% default \ operator 0092 else %% MATLAB and nx > 10 or Octave and nx > 2000 0093 lin_solver = 'LU3'; %% LU decomp with 3 output args, AMD ordering 0094 end 0095 end 0096 0097 %% do Newton iterations 0098 while (~converged && i < max_it) 0099 %% update iteration counter 0100 i = i + 1; 0101 0102 %% evaluate Jacobian 0103 dImis_dQ = sparse(1:n, 1:n, 1j./conj(V), n, n); 0104 [dImis_dVr, dImis_dVi] = dImis_dV(Sb, Ybus, V, 1); 0105 dImis_dVi(:, pv) = ... 0106 dImis_dVi(:, pv) * sparse(1:npv, 1:npv, real(V(pv)), npv, npv) - ... 0107 dImis_dVr(:, pv) * sparse(1:npv, 1:npv, imag(V(pv)), npv, npv); 0108 dImis_dVr(:, pv) = dImis_dQ(:, pv); 0109 0110 %% handling of derivatives for voltage dependent loads 0111 %% (not yet implemented) goes here 0112 0113 j11 = real(dImis_dVr([pv; pq], [pv; pq])); 0114 j12 = real(dImis_dVi([pv; pq], [pv; pq])); 0115 j21 = imag(dImis_dVr([pv; pq], [pv; pq])); 0116 j22 = imag(dImis_dVi([pv; pq], [pv; pq])); 0117 0118 J = [ j11 j12; 0119 j21 j22; ]; 0120 0121 %% compute update step 0122 dx = mplinsolve(J, -F, lin_solver); 0123 0124 %% update voltage 0125 if npv 0126 Va(pv) = Va(pv) + dx(j5:j6); 0127 Sb(pv) = real(Sb(pv)) + 1j * (imag(Sb(pv)) + dx(j1:j2)); 0128 end 0129 if npq 0130 Vm(pq) = Vm(pq) + (real(V(pq))./Vm(pq)) .* dx(j3:j4) ... 0131 + (imag(V(pq))./Vm(pq)) .* dx(j7:j8); 0132 Va(pq) = Va(pq) + (real(V(pq))./(Vm(pq).^2)) .* dx(j7:j8) ... 0133 - (imag(V(pq))./(Vm(pq).^2)) .* dx(j3:j4); 0134 end 0135 V = Vm .* exp(1j * Va); 0136 Vm = abs(V); %% update Vm and Va again in case 0137 Va = angle(V); %% we wrapped around with a negative Vm 0138 0139 %% evalute F(x) 0140 mis = Ybus * V - conj(Sb ./ V); 0141 F = [ real(mis([pv; pq])); 0142 imag(mis([pv; pq])) ]; 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, hybrid) 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, hybrid) did not converge in %d iterations.\n', i); 0160 end 0161 end