CPF_CORRECTOR Solves the corrector step of a continuation power flow using a
full Newton method with selected parameterization scheme.
[V, CONVERGED, I, LAM] = CPF_CORRECTOR(YBUS, SBUS, V0, REF, PV, PQ, ...
LAM0, SXFR, VPRV, LPRV, Z, STEP, PARAMETERIZATION, MPOPT)
solves for bus voltages and lambda given the full system admittance
matrix (for all buses), the complex bus power injection vector (for
all buses), the initial vector of complex bus voltages, and column
vectors with the lists of bus indices for the swing bus, PV buses, and
PQ buses, respectively. 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. MPOPT is a MATPOWER options struct which can be used to
set the termination tolerance, maximum number of iterations, and
output options (see MPOPTION for details). Uses default options if
this parameter is not given. Returns the final complex voltages, a
flag which indicates whether it converged or not, the number
of iterations performed, and the final lambda.
The extra continuation inputs are LAM0 (initial predicted lambda),
SXFR ([delP+j*delQ] transfer/loading vector for all buses), VPRV
(final complex V corrector solution from previous continuation step),
LAMPRV (final lambda corrector solution from previous continuation step),
Z (normalized predictor for all buses), and STEP (continuation step size).
The extra continuation output is LAM (final corrector lambda).
See also RUNCPF.0001 function [V, converged, i, lam] = cpf_corrector(Ybus, Sbus, V0, ref, pv, pq, ... 0002 lam0, Sxfr, Vprv, lamprv, z, step, parameterization, mpopt) 0003 %CPF_CORRECTOR Solves the corrector step of a continuation power flow using a 0004 % full Newton method with selected parameterization scheme. 0005 % [V, CONVERGED, I, LAM] = CPF_CORRECTOR(YBUS, SBUS, V0, REF, PV, PQ, ... 0006 % LAM0, SXFR, VPRV, LPRV, Z, STEP, PARAMETERIZATION, MPOPT) 0007 % solves for bus voltages and lambda given the full system admittance 0008 % matrix (for all buses), the complex bus power injection vector (for 0009 % all buses), the initial vector of complex bus voltages, and column 0010 % vectors with the lists of bus indices for the swing bus, PV buses, and 0011 % PQ buses, respectively. The bus voltage vector contains the set point 0012 % for generator (including ref bus) buses, and the reference angle of the 0013 % swing bus, as well as an initial guess for remaining magnitudes and 0014 % angles. MPOPT is a MATPOWER options struct which can be used to 0015 % set the termination tolerance, maximum number of iterations, and 0016 % output options (see MPOPTION for details). Uses default options if 0017 % this parameter is not given. Returns the final complex voltages, a 0018 % flag which indicates whether it converged or not, the number 0019 % of iterations performed, and the final lambda. 0020 % 0021 % The extra continuation inputs are LAM0 (initial predicted lambda), 0022 % SXFR ([delP+j*delQ] transfer/loading vector for all buses), VPRV 0023 % (final complex V corrector solution from previous continuation step), 0024 % LAMPRV (final lambda corrector solution from previous continuation step), 0025 % Z (normalized predictor for all buses), and STEP (continuation step size). 0026 % The extra continuation output is LAM (final corrector lambda). 0027 % 0028 % See also RUNCPF. 0029 0030 % MATPOWER 0031 % Copyright (c) 1996-2015 by Power System Engineering Research Center (PSERC) 0032 % by Ray Zimmerman, PSERC Cornell, 0033 % Shrirang Abhyankar, Argonne National Laboratory, 0034 % and Alexander Flueck, IIT 0035 % 0036 % Modified by Alexander J. Flueck, Illinois Institute of Technology 0037 % 2001.02.22 - corrector.m (ver 1.0) based on newtonpf.m (MATPOWER 2.0) 0038 % 0039 % Modified by Shrirang Abhyankar, Argonne National Laboratory 0040 % (Updated to be compatible with MATPOWER version 4.1) 0041 % 0042 % $Id: cpf_corrector.m 2644 2015-03-11 19:34:22Z ray $ 0043 % 0044 % This file is part of MATPOWER. 0045 % Covered by the 3-clause BSD License (see LICENSE file for details). 0046 % See http://www.pserc.cornell.edu/matpower/ for more info. 0047 0048 %% default arguments 0049 if nargin < 14 0050 mpopt = mpoption; 0051 end 0052 0053 %% options 0054 tol = mpopt.pf.tol; 0055 max_it = mpopt.pf.nr.max_it; 0056 0057 %% initialize 0058 converged = 0; 0059 i = 0; 0060 V = V0; 0061 Va = angle(V); 0062 Vm = abs(V); 0063 lam = lam0; %% set lam to initial lam0 0064 0065 %% set up indexing for updating V 0066 npv = length(pv); 0067 npq = length(pq); 0068 nb = length(V); %% number of buses 0069 j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses 0070 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses 0071 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses 0072 j7 = j6 + 1; j8 = j6 + 1; %% j7:j8 - lambda 0073 0074 %% evaluate F(x0, lam0), including Sxfr transfer/loading 0075 mis = V .* conj(Ybus * V) - Sbus - lam*Sxfr; 0076 F = [ real(mis([pv; pq])); 0077 imag(mis(pq)) ]; 0078 0079 %% evaluate P(x0, lambda0) 0080 P = cpf_p(parameterization, step, z, V, lam, Vprv, lamprv, pv, pq); 0081 0082 %% augment F(x,lambda) with P(x,lambda) 0083 F = [ F; 0084 P ]; 0085 0086 %% check tolerance 0087 normF = norm(F, inf); 0088 if mpopt.verbose > 1 0089 fprintf('\n it max P & Q mismatch (p.u.)'); 0090 fprintf('\n---- ---------------------------'); 0091 fprintf('\n%3d %10.3e', i, normF); 0092 end 0093 if normF < tol 0094 converged = 1; 0095 if mpopt.verbose > 1 0096 fprintf('\nConverged!\n'); 0097 end 0098 end 0099 0100 %% do Newton iterations 0101 while (~converged && i < max_it) 0102 %% update iteration counter 0103 i = i + 1; 0104 0105 %% evaluate Jacobian 0106 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); 0107 0108 j11 = real(dSbus_dVa([pv; pq], [pv; pq])); 0109 j12 = real(dSbus_dVm([pv; pq], pq)); 0110 j21 = imag(dSbus_dVa(pq, [pv; pq])); 0111 j22 = imag(dSbus_dVm(pq, pq)); 0112 0113 J = [ j11 j12; 0114 j21 j22; ]; 0115 0116 dF_dlam = -[real(Sxfr([pv; pq])); imag(Sxfr(pq))]; 0117 [dP_dV, dP_dlam] = cpf_p_jac(parameterization, z, V, lam, Vprv, lamprv, pv, pq); 0118 0119 %% augment J with real/imag -Sxfr and z^T 0120 J = [ J dF_dlam; 0121 dP_dV dP_dlam ]; 0122 0123 %% compute update step 0124 dx = -(J \ F); 0125 0126 %% update voltage 0127 if npv 0128 Va(pv) = Va(pv) + dx(j1:j2); 0129 end 0130 if npq 0131 Va(pq) = Va(pq) + dx(j3:j4); 0132 Vm(pq) = Vm(pq) + dx(j5:j6); 0133 end 0134 V = Vm .* exp(1j * Va); 0135 Vm = abs(V); %% update Vm and Va again in case 0136 Va = angle(V); %% we wrapped around with a negative Vm 0137 0138 %% update lambda 0139 lam = lam + dx(j7:j8); 0140 0141 %% evalute F(x, lam) 0142 mis = V .* conj(Ybus * V) - Sbus - lam*Sxfr; 0143 F = [ real(mis(pv)); 0144 real(mis(pq)); 0145 imag(mis(pq)) ]; 0146 0147 %% evaluate P(x, lambda) 0148 P = cpf_p(parameterization, step, z, V, lam, Vprv, lamprv, pv, pq); 0149 0150 %% augment F(x,lambda) with P(x,lambda) 0151 F = [ F; 0152 P ]; 0153 0154 %% check for convergence 0155 normF = norm(F, inf); 0156 if mpopt.verbose > 1 0157 fprintf('\n%3d %10.3e', i, normF); 0158 end 0159 if normF < tol 0160 converged = 1; 0161 if mpopt.verbose 0162 fprintf('\nNewton''s method corrector converged in %d iterations.\n', i); 0163 end 0164 end 0165 end 0166 0167 if mpopt.verbose 0168 if ~converged 0169 fprintf('\nNewton''s method corrector did not converge in %d iterations.\n', i); 0170 end 0171 end