CPF_CORRECTOR Solves the corrector step of a continuation power flow
[V, CONVERGED, I, LAM] = CPF_CORRECTOR(YBUS, SBUSB, V_HAT, REF, PV, PQ, ...
LAM_HAT, SBUST, VPRV, LPRV, Z, ...
STEP, PARAMETERIZATION, MPOPT)
Computes the corrector step of a continuation power flow using a
full Newton method with selected parameterization scheme.
Inputs:
YBUS : complex bus admittance matrix
SBUSB : handle of function returning nb x 1 vector of complex
base case injections in p.u. and derivatives w.r.t. |V|
V_HAT : predicted complex bus voltage vector
REF : vector of indices for REF buses
PV : vector of indices of PV buses
PQ : vector of indices of PQ buses
LAM_HAT : predicted scalar lambda
SBUST : handle of function returning nb x 1 vector of complex
target case injections in p.u. and derivatives w.r.t. |V|
VPRV : complex bus voltage vector at previous solution
LAMPRV : scalar lambda value at previous solution
STEP : continuation step length
Z : normalized tangent prediction vector
STEP : continuation step size
PARAMETERIZATION : Value of cpf.parameterization option.
MPOPT : Options struct
Outputs:
V : complex bus voltage solution vector
CONVERGED : Newton iteration count
I : Newton iteration count
LAM : lambda continuation parameter
See also RUNCPF.0001 function [V, converged, i, lam] = cpf_corrector(Ybus, Sbusb, V_hat, ref, pv, pq, ... 0002 lam_hat, Sbust, Vprv, lamprv, z, step, parameterization, mpopt) 0003 %CPF_CORRECTOR Solves the corrector step of a continuation power flow 0004 % [V, CONVERGED, I, LAM] = CPF_CORRECTOR(YBUS, SBUSB, V_HAT, REF, PV, PQ, ... 0005 % LAM_HAT, SBUST, VPRV, LPRV, Z, ... 0006 % STEP, PARAMETERIZATION, MPOPT) 0007 % 0008 % Computes the corrector step of a continuation power flow using a 0009 % full Newton method with selected parameterization scheme. 0010 % 0011 % Inputs: 0012 % YBUS : complex bus admittance matrix 0013 % SBUSB : handle of function returning nb x 1 vector of complex 0014 % base case injections in p.u. and derivatives w.r.t. |V| 0015 % V_HAT : predicted complex bus voltage vector 0016 % REF : vector of indices for REF buses 0017 % PV : vector of indices of PV buses 0018 % PQ : vector of indices of PQ buses 0019 % LAM_HAT : predicted scalar lambda 0020 % SBUST : handle of function returning nb x 1 vector of complex 0021 % target case injections in p.u. and derivatives w.r.t. |V| 0022 % VPRV : complex bus voltage vector at previous solution 0023 % LAMPRV : scalar lambda value at previous solution 0024 % STEP : continuation step length 0025 % Z : normalized tangent prediction vector 0026 % STEP : continuation step size 0027 % PARAMETERIZATION : Value of cpf.parameterization option. 0028 % MPOPT : Options struct 0029 % 0030 % Outputs: 0031 % V : complex bus voltage solution vector 0032 % CONVERGED : Newton iteration count 0033 % I : Newton iteration count 0034 % LAM : lambda continuation parameter 0035 % 0036 % See also RUNCPF. 0037 0038 % MATPOWER 0039 % Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) 0040 % by Ray Zimmerman, PSERC Cornell, 0041 % Shrirang Abhyankar, Argonne National Laboratory, 0042 % and Alexander Flueck, IIT 0043 % 0044 % Modified by Alexander J. Flueck, Illinois Institute of Technology 0045 % 2001.02.22 - corrector.m (ver 1.0) based on newtonpf.m (MATPOWER 2.0) 0046 % 0047 % Modified by Shrirang Abhyankar, Argonne National Laboratory 0048 % (Updated to be compatible with MATPOWER version 4.1) 0049 % 0050 % This file is part of MATPOWER. 0051 % Covered by the 3-clause BSD License (see LICENSE file for details). 0052 % See https://matpower.org for more info. 0053 0054 %% default arguments 0055 if nargin < 14 0056 mpopt = mpoption; 0057 end 0058 0059 %% options 0060 tol = mpopt.pf.tol; 0061 max_it = mpopt.pf.nr.max_it; 0062 0063 %% initialize 0064 converged = 0; 0065 i = 0; 0066 V = V_hat; %% initialize V with predicted V 0067 Va = angle(V); 0068 Vm = abs(V); 0069 lam = lam_hat; %% initialize lam with predicted lam 0070 0071 %% set up indexing for updating V 0072 npv = length(pv); 0073 npq = length(pq); 0074 nb = length(V); %% number of buses 0075 j1 = 1; j2 = npv; %% j1:j2 - V angle of pv buses 0076 j3 = j2 + 1; j4 = j2 + npq; %% j3:j4 - V angle of pq buses 0077 j5 = j4 + 1; j6 = j4 + npq; %% j5:j6 - V mag of pq buses 0078 j7 = j6 + 1; j8 = j6 + 1; %% j7:j8 - lambda 0079 0080 %% evaluate F(x0, lam_hat), including transfer/loading 0081 Sb = Sbusb(Vm); 0082 St = Sbust(Vm); 0083 mis = V .* conj(Ybus * V) - Sb - lam * (St - Sb); 0084 F = [ real(mis([pv; pq])); 0085 imag(mis(pq)) ]; 0086 0087 %% evaluate P(x0, lambda0) 0088 P = cpf_p(parameterization, step, z, V, lam, Vprv, lamprv, pv, pq); 0089 0090 %% augment F(x,lambda) with P(x,lambda) 0091 F = [ F; 0092 P ]; 0093 0094 %% check tolerance 0095 normF = norm(F, inf); 0096 if mpopt.verbose > 1 0097 fprintf('\n it max P & Q mismatch (p.u.)'); 0098 fprintf('\n---- ---------------------------'); 0099 fprintf('\n%3d %10.3e', i, normF); 0100 end 0101 if normF < tol 0102 converged = 1; 0103 if mpopt.verbose > 1 0104 fprintf('\nConverged!\n'); 0105 end 0106 end 0107 0108 %% do Newton iterations 0109 while (~converged && i < max_it) 0110 %% update iteration counter 0111 i = i + 1; 0112 0113 %% evaluate Jacobian 0114 [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V); 0115 [dummy, neg_dSdb_dVm] = Sbusb(Vm); 0116 [dummy, neg_dSdt_dVm] = Sbust(Vm); 0117 dSbus_dVm = dSbus_dVm - neg_dSdb_dVm - lam * (neg_dSdt_dVm - neg_dSdb_dVm); 0118 0119 j11 = real(dSbus_dVa([pv; pq], [pv; pq])); 0120 j12 = real(dSbus_dVm([pv; pq], pq)); 0121 j21 = imag(dSbus_dVa(pq, [pv; pq])); 0122 j22 = imag(dSbus_dVm(pq, pq)); 0123 0124 J = [ j11 j12; 0125 j21 j22; ]; 0126 0127 Sxf = St - Sb; 0128 dF_dlam = -[real(Sxf([pv; pq])); imag(Sxf(pq))]; 0129 [dP_dV, dP_dlam] = cpf_p_jac(parameterization, z, V, lam, Vprv, lamprv, pv, pq); 0130 0131 %% augment J with real/imag -Sxfr and z^T 0132 J = [ J dF_dlam; 0133 dP_dV dP_dlam ]; 0134 0135 %% compute update step 0136 dx = -(J \ F); 0137 0138 %% update voltage 0139 if npv 0140 Va(pv) = Va(pv) + dx(j1:j2); 0141 end 0142 if npq 0143 Va(pq) = Va(pq) + dx(j3:j4); 0144 Vm(pq) = Vm(pq) + dx(j5:j6); 0145 end 0146 V = Vm .* exp(1j * Va); 0147 Vm = abs(V); %% update Vm and Va again in case 0148 Va = angle(V); %% we wrapped around with a negative Vm 0149 0150 %% update lambda 0151 lam = lam + dx(j7:j8); 0152 0153 %% evalute F(x, lam) 0154 Sb = Sbusb(Vm); 0155 St = Sbust(Vm); 0156 mis = V .* conj(Ybus * V) - Sb - lam * (St - Sb); 0157 F = [ real(mis(pv)); 0158 real(mis(pq)); 0159 imag(mis(pq)) ]; 0160 0161 %% evaluate P(x, lambda) 0162 P = cpf_p(parameterization, step, z, V, lam, Vprv, lamprv, pv, pq); 0163 0164 %% augment F(x,lambda) with P(x,lambda) 0165 F = [ F; 0166 P ]; 0167 0168 %% check for convergence 0169 normF = norm(F, inf); 0170 if mpopt.verbose > 1 0171 fprintf('\n%3d %10.3e', i, normF); 0172 end 0173 if normF < tol 0174 converged = 1; 0175 if mpopt.verbose 0176 fprintf('\nNewton''s method corrector converged in %d iterations.\n', i); 0177 end 0178 end 0179 end 0180 0181 if mpopt.verbose 0182 if ~converged 0183 fprintf('\nNewton''s method corrector did not converge in %d iterations.\n', i); 0184 end 0185 end