RUNCPF Runs a full AC continuation power flow
[RESULTS, SUCCESS] = RUNCPF(BASECASEDATA, TARGETCASEDATA, ...
MPOPT, FNAME, SOLVEDCASE)
Runs a full AC continuation power flow using a normalized tangent
predictor and selected parameterization scheme, optionally
returning a RESULTS struct and SUCCESS flag. Step size can be
fixed or adaptive.
Inputs (all are optional):
BASECASEDATA : either a MATPOWER case struct or a string containing
the name of the file with the case data defining the base loading
and generation (default is 'case9')
(see also CASEFORMAT and LOADCASE)
TARGETCASEDATA : either a MATPOWER case struct or a string
containing the name of the file with the case data defining the
target loading and generation (default is 'case9target')
MPOPT : MATPOWER options struct to override default options
can be used to specify the parameterization, output options,
termination criteria, and more (see also MPOPTION).
FNAME : name of a file to which the pretty-printed output will
be appended
SOLVEDCASE : name of file to which the solved case will be saved
in MATPOWER case format (M-file will be assumed unless the
specified name ends with '.mat')
Outputs (all are optional):
RESULTS : results struct, with the following fields:
(all fields from the input MATPOWER case, i.e. bus, branch,
gen, etc., but with solved voltages, power flows, etc.)
order - info used in external <-> internal data conversion
et - elapsed time in seconds
success - success flag, 1 = succeeded, 0 = failed
cpf - CPF output struct whose content depends on any
user callback functions, where default contains fields:
done_msg - string with message describing cause of
continuation termination
iterations - number of continuation steps performed
lam - (nsteps+1) row vector of lambda values from
correction steps
lam_hat - (nsteps+1) row vector of lambda values from
prediction steps
max_lam - maximum value of lambda in RESULTS.cpf.lam
steps - (nsteps+1) row vector of stepsizes taken at each
continuation step
V - (nb x nsteps+1) complex bus voltages from
correction steps
V_hat - (nb x nsteps+1) complex bus voltages from
prediction steps
events - an array of structs of size nevents with the
following fields:
k - continuation step number at which event was located
name - name of event
idx - index(es) of critical elements in corresponding
event function, e.g. index of generator reaching VAr
limit
msg - descriptive text detailing the event
SUCCESS : the success flag can additionally be returned as
a second output argument
Calling syntax options:
results = runcpf;
results = runcpf(basecasedata, targetcasedata);
results = runcpf(basecasedata, targetcasedata, mpopt);
results = runcpf(basecasedata, targetcasedata, mpopt, fname);
results = runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase)
[results, success] = runcpf(...);
If the 'cpf.enforce_q_lims' option is set to true (default is false) then,
if any generator reaches its reactive power limits during the AC
continuation power flow,
- the corresponding bus is converted to a PQ bus, and the problem
is modified to eliminate further reactive transfer on this bus
- the voltage magnitude at the bus will deviate from the specified
setpoint to satisfy the reactive power limit,
- if the reference bus is converted to PQ, further real power transfer
for the bus is also eliminated, and the first remaining PV bus is
selected as the new slack, resulting in the transfers at both
reference buses potentially deviating from the specified values
- if all reference and PV buses are converted to PQ, RUNCPF terminates
with an infeasibility message.
If the 'cpf.enforce_p_lims' option is set to true (default is fals) then,
if any generator reaches its maximum active power limit during the AC
continuation power flow,
- the problem is modified to eliminate further active transfer by
this generator
- if the generator was at the reference bus, it is converted to PV
and the first remaining PV bus is selected as the new slack.
Possible CPF termination modes:
when cpf.stop_at == 'NOSE'
- Reached steady state loading limit
- Nose point eliminated by limit induced bifurcation
when cpf.stop_at == 'FULL'
- Traced full continuation curve
when cpf.stop_at == <target_lam_val>
- Reached desired lambda
when cpf.enforce_p_lims == true
- All generators at PMAX
when cpf.enforce_q_lims == true
- No REF or PV buses remaining
other
- Base case power flow did not converge
- Base and target case have identical load and generation
- Corrector did not converge
- Could not locate <event_name> event
- Too many rollback steps triggered by callbacks
Examples:
results = runcpf('case9', 'case9target');
results = runcpf('case9', 'case9target', ...
mpoption('cpf.adapt_step', 1));
results = runcpf('case9', 'case9target', ...
mpoption('cpf.enforce_q_lims', 1));
results = runcpf('case9', 'case9target', ...
mpoption('cpf.stop_at', 'FULL'));
See also MPOPTION, RUNPF.0001 function [res, suc] = ... 0002 runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase) 0003 %RUNCPF Runs a full AC continuation power flow 0004 % [RESULTS, SUCCESS] = RUNCPF(BASECASEDATA, TARGETCASEDATA, ... 0005 % MPOPT, FNAME, SOLVEDCASE) 0006 % 0007 % Runs a full AC continuation power flow using a normalized tangent 0008 % predictor and selected parameterization scheme, optionally 0009 % returning a RESULTS struct and SUCCESS flag. Step size can be 0010 % fixed or adaptive. 0011 % 0012 % Inputs (all are optional): 0013 % BASECASEDATA : either a MATPOWER case struct or a string containing 0014 % the name of the file with the case data defining the base loading 0015 % and generation (default is 'case9') 0016 % (see also CASEFORMAT and LOADCASE) 0017 % TARGETCASEDATA : either a MATPOWER case struct or a string 0018 % containing the name of the file with the case data defining the 0019 % target loading and generation (default is 'case9target') 0020 % MPOPT : MATPOWER options struct to override default options 0021 % can be used to specify the parameterization, output options, 0022 % termination criteria, and more (see also MPOPTION). 0023 % FNAME : name of a file to which the pretty-printed output will 0024 % be appended 0025 % SOLVEDCASE : name of file to which the solved case will be saved 0026 % in MATPOWER case format (M-file will be assumed unless the 0027 % specified name ends with '.mat') 0028 % 0029 % Outputs (all are optional): 0030 % RESULTS : results struct, with the following fields: 0031 % (all fields from the input MATPOWER case, i.e. bus, branch, 0032 % gen, etc., but with solved voltages, power flows, etc.) 0033 % order - info used in external <-> internal data conversion 0034 % et - elapsed time in seconds 0035 % success - success flag, 1 = succeeded, 0 = failed 0036 % cpf - CPF output struct whose content depends on any 0037 % user callback functions, where default contains fields: 0038 % done_msg - string with message describing cause of 0039 % continuation termination 0040 % iterations - number of continuation steps performed 0041 % lam - (nsteps+1) row vector of lambda values from 0042 % correction steps 0043 % lam_hat - (nsteps+1) row vector of lambda values from 0044 % prediction steps 0045 % max_lam - maximum value of lambda in RESULTS.cpf.lam 0046 % steps - (nsteps+1) row vector of stepsizes taken at each 0047 % continuation step 0048 % V - (nb x nsteps+1) complex bus voltages from 0049 % correction steps 0050 % V_hat - (nb x nsteps+1) complex bus voltages from 0051 % prediction steps 0052 % events - an array of structs of size nevents with the 0053 % following fields: 0054 % k - continuation step number at which event was located 0055 % name - name of event 0056 % idx - index(es) of critical elements in corresponding 0057 % event function, e.g. index of generator reaching VAr 0058 % limit 0059 % msg - descriptive text detailing the event 0060 % SUCCESS : the success flag can additionally be returned as 0061 % a second output argument 0062 % 0063 % Calling syntax options: 0064 % results = runcpf; 0065 % results = runcpf(basecasedata, targetcasedata); 0066 % results = runcpf(basecasedata, targetcasedata, mpopt); 0067 % results = runcpf(basecasedata, targetcasedata, mpopt, fname); 0068 % results = runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase) 0069 % [results, success] = runcpf(...); 0070 % 0071 % If the 'cpf.enforce_q_lims' option is set to true (default is false) then, 0072 % if any generator reaches its reactive power limits during the AC 0073 % continuation power flow, 0074 % - the corresponding bus is converted to a PQ bus, and the problem 0075 % is modified to eliminate further reactive transfer on this bus 0076 % - the voltage magnitude at the bus will deviate from the specified 0077 % setpoint to satisfy the reactive power limit, 0078 % - if the reference bus is converted to PQ, further real power transfer 0079 % for the bus is also eliminated, and the first remaining PV bus is 0080 % selected as the new slack, resulting in the transfers at both 0081 % reference buses potentially deviating from the specified values 0082 % - if all reference and PV buses are converted to PQ, RUNCPF terminates 0083 % with an infeasibility message. 0084 % 0085 % If the 'cpf.enforce_p_lims' option is set to true (default is fals) then, 0086 % if any generator reaches its maximum active power limit during the AC 0087 % continuation power flow, 0088 % - the problem is modified to eliminate further active transfer by 0089 % this generator 0090 % - if the generator was at the reference bus, it is converted to PV 0091 % and the first remaining PV bus is selected as the new slack. 0092 % 0093 % Possible CPF termination modes: 0094 % when cpf.stop_at == 'NOSE' 0095 % - Reached steady state loading limit 0096 % - Nose point eliminated by limit induced bifurcation 0097 % when cpf.stop_at == 'FULL' 0098 % - Traced full continuation curve 0099 % when cpf.stop_at == <target_lam_val> 0100 % - Reached desired lambda 0101 % when cpf.enforce_p_lims == true 0102 % - All generators at PMAX 0103 % when cpf.enforce_q_lims == true 0104 % - No REF or PV buses remaining 0105 % other 0106 % - Base case power flow did not converge 0107 % - Base and target case have identical load and generation 0108 % - Corrector did not converge 0109 % - Could not locate <event_name> event 0110 % - Too many rollback steps triggered by callbacks 0111 % 0112 % Examples: 0113 % results = runcpf('case9', 'case9target'); 0114 % results = runcpf('case9', 'case9target', ... 0115 % mpoption('cpf.adapt_step', 1)); 0116 % results = runcpf('case9', 'case9target', ... 0117 % mpoption('cpf.enforce_q_lims', 1)); 0118 % results = runcpf('case9', 'case9target', ... 0119 % mpoption('cpf.stop_at', 'FULL')); 0120 % 0121 % See also MPOPTION, RUNPF. 0122 0123 % MATPOWER 0124 % Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) 0125 % by Ray Zimmerman, PSERC Cornell, 0126 % Shrirang Abhyankar, Argonne National Laboratory, 0127 % and Alexander Flueck, IIT 0128 % 0129 % This file is part of MATPOWER. 0130 % Covered by the 3-clause BSD License (see LICENSE file for details). 0131 % See http://www.pserc.cornell.edu/matpower/ for more info. 0132 0133 %%----- initialize ----- 0134 %% define named indices into bus, gen, branch matrices 0135 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0136 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0137 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0138 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0139 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0140 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0141 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0142 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0143 0144 %% default arguments 0145 if nargin < 5 0146 solvedcase = ''; %% don't save solved case 0147 if nargin < 4 0148 fname = ''; %% don't print results to a file 0149 if nargin < 3 0150 mpopt = mpoption; %% use default options 0151 if nargin < 2 0152 targetcasedata = 'case9target'; 0153 if nargin < 1 0154 basecasedata = 'case9'; %% default data file is 'case9.m' 0155 end 0156 end 0157 end 0158 end 0159 end 0160 0161 %% options 0162 step = mpopt.cpf.step; %% continuation step length 0163 parm = mpopt.cpf.parameterization; %% parameterization 0164 adapt_step = mpopt.cpf.adapt_step; %% use adaptive step size? 0165 qlim = mpopt.cpf.enforce_q_lims; %% enforce reactive limits 0166 plim = mpopt.cpf.enforce_p_lims; %% enforce active limits 0167 0168 %% register event functions (for event detection) 0169 %% and CPF callback functions (for event handling and other tasks) 0170 cpf_events = []; 0171 cpf_callbacks = []; 0172 %% to handle CPF termination 0173 if ischar(mpopt.cpf.stop_at) && strcmp(mpopt.cpf.stop_at, 'NOSE'); 0174 cpf_events = cpf_register_event(cpf_events, 'NOSE', 'cpf_nose_event', mpopt.cpf.nose_tol, 1); 0175 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_nose_event_cb', 51); 0176 else %% FULL or target lambda 0177 cpf_events = cpf_register_event(cpf_events, 'TARGET_LAM', 'cpf_target_lam_event', mpopt.cpf.target_lam_tol, 1); 0178 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_target_lam_event_cb', 50); 0179 end 0180 %% to handle reactive power limits 0181 if qlim 0182 cpf_events = cpf_register_event(cpf_events, 'QLIM', 'cpf_qlim_event', mpopt.cpf.q_lims_tol, 1); 0183 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_qlim_event_cb', 41); 0184 end 0185 %% to handle active power limits 0186 if plim 0187 cpf_events = cpf_register_event(cpf_events, 'PLIM', 'cpf_plim_event', mpopt.cpf.p_lims_tol, 1); 0188 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_plim_event_cb', 40); 0189 end 0190 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_default_callback', 0); 0191 %% user callbacks 0192 if ~isempty(mpopt.cpf.user_callback) 0193 %% convert to cell array, if necessary 0194 if iscell(mpopt.cpf.user_callback) 0195 user_callbacks = mpopt.cpf.user_callback; 0196 else 0197 user_callbacks = {mpopt.cpf.user_callback}; 0198 end 0199 for k = 1:length(user_callbacks) 0200 %% convert each to struct, if necessary 0201 if isstruct(user_callbacks{k}) 0202 ucb = user_callbacks{k}; 0203 else 0204 ucb = struct('fcn', user_callbacks{k}); 0205 end 0206 %% set default priority, args if not provided 0207 if ~isfield(ucb, 'priority') 0208 ucb.priority = []; 0209 end 0210 if ~isfield(ucb, 'args') 0211 ucb.args = []; 0212 end 0213 %% register the user callback 0214 cpf_callbacks = cpf_register_callback(cpf_callbacks, ... 0215 ucb.fcn, ucb.priority, ucb.args); 0216 end 0217 end 0218 nef = length(cpf_events); %% number of event functions registered 0219 ncb = length(cpf_callbacks); %% number of callback functions registered 0220 0221 %% set power flow options 0222 if mpopt.verbose > 4 0223 mpopt_pf = mpoption(mpopt, 'verbose', 2); 0224 else 0225 mpopt_pf = mpoption(mpopt, 'verbose', 0); 0226 end 0227 mpopt_pf = mpoption(mpopt_pf, 'pf.enforce_q_lims', mpopt.cpf.enforce_q_lims); 0228 0229 %% load base case data 0230 mpcb = loadcase(basecasedata); 0231 0232 %% clip base active generator outputs to PMAX, if necessary 0233 idx_pmax = []; %% indices of generators clipped at PMAX 0234 if plim 0235 idx_pmax = find( mpcb.gen(:, GEN_STATUS) > 0 & ... 0236 mpcb.gen(:, PG) - mpcb.gen(:, PMAX) > -mpopt.cpf.p_lims_tol); 0237 if mpopt.verbose && ~isempty(idx_pmax) 0238 fprintf('base case real power output of gen %d reduced from %g to %g MW (PMAX)\n', ... 0239 [idx_pmax mpcb.gen(idx_pmax, PG) mpcb.gen(idx_pmax, PMAX)]'); 0240 end 0241 mpcb.gen(idx_pmax, PG) = mpcb.gen(idx_pmax, PMAX); 0242 end 0243 0244 %% run base case power flow 0245 [mpcb, suc] = runpf(mpcb, mpopt_pf); 0246 if suc 0247 done = struct('flag', 0, 'msg', ''); 0248 else 0249 done = struct('flag', 1, 'msg', 'Base case power flow did not converge.'); 0250 results = mpcb; 0251 results.cpf = struct(); 0252 end 0253 0254 if ~done.flag 0255 %% read target case data 0256 mpct = loadcase(targetcasedata); 0257 if size(mpct.branch,2) < QT %% add zero columns to branch for flows if needed 0258 mpct.branch = [ mpct.branch zeros(size(mpct.branch, 1), QT-size(mpct.branch,2)) ]; 0259 end 0260 0261 %% convert both to internal indexing 0262 mpcb = ext2int(mpcb); 0263 mpct = ext2int(mpct); 0264 i2e_gen = mpct.order.gen.i2e; 0265 nb = size(mpcb.bus, 1); 0266 0267 %% get bus index lists of each type of bus 0268 [ref, pv, pq] = bustypes(mpcb.bus, mpcb.gen); 0269 0270 %% generator info 0271 %% find generators that are ON and at voltage-controlled buses 0272 ong = find(mpcb.gen(:, GEN_STATUS) > 0 ... 0273 & mpcb.bus(mpcb.gen(:, GEN_BUS), BUS_TYPE) ~= PQ); 0274 gbus = mpcb.gen(ong, GEN_BUS); %% what buses are they at? 0275 0276 %% make sure target case is same as base case w.r.t 0277 %% bus types, GEN_STATUS, Qg and slack Pg 0278 mpct.bus(:, BUS_TYPE) = mpcb.bus(:, BUS_TYPE); 0279 ont = find(mpct.gen(:, GEN_STATUS) > 0 ... 0280 & mpct.bus(mpct.gen(:, GEN_BUS), BUS_TYPE) ~= PQ); 0281 if length(ong) ~= length(ont) || any(ong ~= ont) 0282 error('runcpf: GEN_STATUS of all generators must be the same in base and target cases'); 0283 end 0284 mpct.gen(ong, QG) = mpcb.gen(ong, QG); 0285 for k = 1:length(ref) 0286 refgen = find(gbus == ref(k)); 0287 mpct.gen(ong(refgen), PG) = mpcb.gen(ong(refgen), PG); 0288 end 0289 0290 %% zero transfers for gens that exceed PMAX limits, if necessary 0291 if plim 0292 idx_pmax = find( mpcb.gen(:, GEN_STATUS) > 0 & ... 0293 mpcb.gen(:, PG) - mpcb.gen(:, PMAX) > -mpopt.cpf.p_lims_tol & ... 0294 mpct.gen(:, PG) - mpct.gen(:, PMAX) > -mpopt.cpf.p_lims_tol); 0295 if ~isempty(idx_pmax) 0296 if mpopt.verbose 0297 fprintf('target case real power output of gen %d reduced from %g to %g MW (PMAX)\n', ... 0298 [i2e_gen(idx_pmax) mpct.gen(idx_pmax, PG) mpcb.gen(idx_pmax, PG)]'); 0299 end 0300 mpct.gen(idx_pmax, PG) = mpcb.gen(idx_pmax, PG); 0301 end 0302 end 0303 0304 %%----- run the continuation power flow ----- 0305 t0 = clock; 0306 if mpopt.verbose 0307 v = mpver('all'); 0308 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); 0309 fprintf(' -- AC Continuation Power Flow\n'); 0310 if mpopt.verbose > 1 0311 fprintf('step %3d : lambda = %6.3f, %2d Newton steps\n', 0, 0, mpcb.iterations); 0312 end 0313 end 0314 0315 %% build admittance matrices 0316 [Ybus, Yf, Yt] = makeYbus(mpcb.baseMVA, mpcb.bus, mpcb.branch); 0317 0318 %% functions for computing base and target case V-dependent complex bus 0319 %% power injections: (generation - load) 0320 Sbusb = @(Vm)makeSbus(mpcb.baseMVA, mpcb.bus, mpcb.gen, mpopt, Vm); 0321 Sbust = @(Vm)makeSbus(mpct.baseMVA, mpct.bus, mpct.gen, mpopt, Vm); 0322 0323 %% initialize variables 0324 cont_steps = 0; 0325 lam = 0; 0326 V = mpcb.bus(:, VM) .* exp(sqrt(-1) * pi/180 * mpcb.bus(:, VA)); 0327 rollback = 0; %% flag to indicate that a step must be rolled back 0328 locating = 0; %% flag to indicate that an event has interval was detected, 0329 %% but the event has not yet been located 0330 rb_cnt_ef = 0; %% counter for rollback steps triggered by event function intervals 0331 rb_cnt_cb = 0; %% counter for rollback steps triggered directly by callbacks 0332 0333 %% initialize tangent predictor: z = [dx;dlam] 0334 z = [zeros(2*nb, 1); 1]; 0335 z = cpf_tangent(V, lam, Ybus, Sbusb, Sbust, pv, pq, z, V, lam, parm); 0336 0337 %% initialize state for current continuation step 0338 cx = struct(... %% current state 0339 'lam_hat', lam, ... %% predicted lambda 0340 'V_hat', V, ... %% predicted V 0341 'lam', lam, ... %% corrected lambda 0342 'V', V, ... %% corrected V 0343 'z', z, ... %% normalized tangent predictor 0344 'default_step', step, ... %% default step size 0345 'default_parm', parm, ... %% default parameterization 0346 'this_step', [], ... %% step size for this step only 0347 'this_parm', [], ... %% parameterization for this step only 0348 'step', step, ... %% current step size 0349 'parm', parm, ... %% current parameterization 0350 'events', [], ... %% event log 0351 'cb', struct(), ... %% user state, for callbacks (replaces cb_state) 0352 'ef', {cell(nef, 1)} ... %% event function values 0353 ); 0354 0355 %% input args for callbacks 0356 cb_data = struct( ... 0357 'mpc_base', mpcb, ... %% MATPOWER case struct - base case 0358 'mpc_target', mpct, ... %% MATPOWER case struct - target case 0359 'Sbusb', Sbusb, ... %% function for computing base bus inj 0360 'Sbust', Sbust, ... %% function for computing target bus inj 0361 'Ybus', Ybus, ... %% bus admittance matrix 0362 'Yf', Yf, ... %% branch admittance matrix, from end 0363 'Yt', Yt, ... %% branch admittance matrix, to end 0364 'ref', ref, ... %% vector of ref bus indices 0365 'pv', pv, ... %% vector of PV bus indices 0366 'pq', pq, ... %% vector of PQ bus indices 0367 'idx_pmax', idx_pmax, ... %% vector of gen indices of gens at PMAX 0368 'mpopt', mpopt ); %% MATPOWER option struct 0369 0370 %% initialize event function values 0371 for k = 1:nef 0372 cx.ef{k} = cpf_events(k).fcn(cb_data, cx); 0373 end 0374 0375 %% invoke callbacks - "initialize" context 0376 for k = 1:ncb 0377 [nx, cx, done, rollback, evnts, cb_data] = cpf_callbacks(k).fcn( ... 0378 cont_steps, cx, cx, cx, done, 0, [], cb_data, cpf_callbacks(k).args); 0379 end 0380 0381 %% check for case with no transfer 0382 if norm(Sbust(abs(cx.V)) - Sbusb(abs(cx.V))) == 0 0383 done.flag = 1; 0384 done.msg = 'Base case and target case have identical load and generation'; 0385 end 0386 0387 cont_steps = cont_steps + 1; 0388 px = cx; %% initialize state for previous continuation step 0389 while ~done.flag 0390 %% initialize next candidate with current state 0391 nx = cx; 0392 0393 %% prediction for next step 0394 [nx.V_hat, nx.lam_hat] = cpf_predictor(cx.V, cx.lam, cx.z, cx.step, cb_data.pv, cb_data.pq); 0395 0396 %% correction 0397 [nx.V, success, i, nx.lam] = cpf_corrector(Ybus, cb_data.Sbusb, nx.V_hat, cb_data.ref, cb_data.pv, cb_data.pq, ... 0398 nx.lam_hat, cb_data.Sbust, cx.V, cx.lam, cx.z, cx.step, cx.parm, mpopt_pf); 0399 if ~success 0400 done.flag = 1; 0401 done.msg = sprintf('Corrector did not converge in %d iterations.', i); 0402 if mpopt.verbose 0403 fprintf('step %3d : stepsize = %-9.3g lambda = %6.3f corrector did not converge in %d iterations\n', cont_steps, cx.step, nx.lam, i); 0404 end 0405 cont_steps = cont_steps - 1; 0406 break; 0407 end 0408 0409 %% compute new tangent direction, based on a previous state: tx 0410 if nx.step == 0 %% if this is a re-do step, cx and nx are the same 0411 tx = px; %% so use px as the previous state 0412 else %% otherwise 0413 tx = cx; %% use cx as the previous state 0414 end 0415 nx.z = cpf_tangent(nx.V, nx.lam, Ybus, cb_data.Sbusb, cb_data.Sbust, cb_data.pv, cb_data.pq, ... 0416 cx.z, tx.V, tx.lam, nx.parm); 0417 0418 %% detect events 0419 for k = 1:nef 0420 nx.ef{k} = cpf_events(k).fcn(cb_data, nx); %% update event functions 0421 end 0422 [rollback, evnts, nx.ef] = cpf_detect_events(cpf_events, nx.ef, cx.ef, nx.step, mpopt.verbose); 0423 0424 %% adjust step-size to locate event function zero, if necessary 0425 if rollback %% current step overshot 0426 %% rollback and initialize next step size based on rollback and previous 0427 rx = nx; %% save state we're rolling back from 0428 rx.evnts = evnts; %% and critical event info 0429 cx.this_step = evnts.step_scale * rx.step; 0430 cx.this_parm = rx.parm; %% keep same parameterization as last step 0431 locating = 1; %% enter "locating" mode (or stay in it) 0432 rb_cnt_ef = rb_cnt_ef + 1; %% increment rollback counter for ef intervals 0433 if rb_cnt_ef > 26 0434 done.flag = 1; 0435 done.msg = sprintf('Could not locate %s event!', evnts.name); 0436 end 0437 if mpopt.verbose > 3 0438 loc_msg = sprintf('OVERSHOOT : f = [%g, <<%g>>], step <-- %.4g', ... 0439 cx.ef{evnts.eidx}(evnts.idx(1)), ... 0440 rx.ef{evnts.eidx}(evnts.idx(1)), cx.this_step); 0441 end 0442 elseif locating 0443 if evnts(1).zero %% found the zero! 0444 %% step size will be reset to previously used default step size 0445 locating = 0; %% exit "locating" mode 0446 rb_cnt_ef = 0; %% reset rollback counter for ef intervals 0447 if mpopt.verbose > 3 0448 loc_msg = sprintf('ZERO! : f = %g, step <-- %.4g', ... 0449 nx.ef{rx.evnts.eidx}(rx.evnts.idx(1)), nx.default_step); 0450 end 0451 else %% prev rollback undershot 0452 %% initialize next step size based on critical event function 0453 %% values from prev rollback step and current step 0454 rx_ef = rx.ef{rx.evnts.eidx}(rx.evnts.idx(1)); 0455 cx_ef = nx.ef{rx.evnts.eidx}(rx.evnts.idx(1)); 0456 step_scale = cx_ef / (cx_ef - rx_ef); 0457 nx.this_step = step_scale * (rx.step - nx.step); 0458 rb_cnt_ef = 0; %% reset rollback counter for ef intervals 0459 if mpopt.verbose > 3 0460 loc_msg = sprintf('UNDERSHOOT : f [<<%g>>, %g], step <-- %.4g', ... 0461 cx_ef, rx_ef, nx.this_step); 0462 end 0463 end 0464 else 0465 loc_msg = ''; 0466 end 0467 0468 %% invoke callbacks - "iterations" context 0469 rb = rollback; 0470 for k = 1:ncb 0471 [nx, cx, done, rollback, evnts, cb_data] = cpf_callbacks(k).fcn( ... 0472 cont_steps, nx, cx, px, done, rollback, evnts, cb_data, cpf_callbacks(k).args); 0473 end 0474 if ~rb && rollback %% rollback triggered by callback (vs event function interval) 0475 rb_cnt_cb = rb_cnt_cb + 1; %% increment rollback counter for callbacks 0476 if rb_cnt_cb > 26 0477 done.flag = 1; 0478 done.msg = 'Too many rollback steps triggered by callbacks!'; 0479 end 0480 else 0481 rb_cnt_cb = 0; %% reset rollback counter for callbacks 0482 end 0483 0484 %% print iteration information 0485 if mpopt.verbose > 1 0486 %% set label for rollback step counter 0487 if rb_cnt_ef 0488 sub_step = char('a' + rb_cnt_ef - 1); 0489 elseif rb_cnt_cb 0490 sub_step = char('A' + rb_cnt_cb - 1); 0491 else 0492 sub_step = ' '; 0493 end 0494 if mpopt.verbose > 4 0495 fprintf('step %3d%s : stepsize = %-9.3g lambda = %6.3f', cont_steps, sub_step, cx.step, nx.lam); 0496 else 0497 fprintf('step %3d%s : stepsize = %-9.3g lambda = %6.3f %2d corrector Newton steps', cont_steps, sub_step, cx.step, nx.lam, i); 0498 end 0499 if rollback 0500 fprintf(' ^ ROLLBACK\n'); 0501 else 0502 fprintf('\n'); 0503 end 0504 if mpopt.verbose > 3 && ~isempty(loc_msg) 0505 fprintf(' LOCATING -- %s\n', loc_msg); 0506 end 0507 end 0508 0509 %% log events 0510 for k = 1:length(evnts) 0511 if evnts(k).log 0512 e = struct( 'k', cont_steps, ... 0513 'name', evnts(k).name, ... 0514 'idx', evnts(k).idx, ... 0515 'msg', evnts(k).msg ); 0516 if isempty(nx.events) 0517 nx.events = e; 0518 else 0519 nx.events(end+1) = e; 0520 end 0521 end 0522 if (mpopt.verbose > 2 && evnts(k).log) || ... 0523 (mpopt.verbose > 3 && evnts(k).eidx) 0524 fprintf(' %s\n', evnts(k).msg); 0525 end 0526 end 0527 0528 %% adapt stepsize if requested and not terminating, locating a zero 0529 %% or re-doing a step after changing the problem data 0530 if adapt_step && ~done.flag && ~locating && ~evnts(1).zero && nx.step ~= 0 0531 cpf_error = norm([angle(nx.V( cb_data.pq)); abs(nx.V( [cb_data.pv;cb_data.pq])); nx.lam] - ... 0532 [angle(nx.V_hat(cb_data.pq)); abs(nx.V_hat([cb_data.pv;cb_data.pq])); nx.lam_hat], inf); 0533 0534 %% new nominal step size is current size * tol/err, but we reduce 0535 %% the change from the current size by a damping factor and limit 0536 %% increases to a factor of 2 0537 step_scale = min(2, 1 + mpopt.cpf.adapt_step_damping * ... 0538 (mpopt.cpf.adapt_step_tol/cpf_error - 1)); 0539 nx.default_step = nx.step * step_scale; 0540 0541 %% limit step-size 0542 if nx.default_step > mpopt.cpf.step_max 0543 nx.default_step = mpopt.cpf.step_max; 0544 end 0545 if nx.default_step < mpopt.cpf.step_min 0546 nx.default_step = mpopt.cpf.step_min; 0547 end 0548 end 0549 0550 %% if this is a normal step 0551 if ~rollback 0552 px = cx; %% save current state before update 0553 cx = nx; %% update current state to next candidate 0554 if ~done.flag 0555 cont_steps = cont_steps + 1; 0556 end 0557 end 0558 0559 %% set step size and parameterization, from one-time or defaults 0560 if isempty(cx.this_step) 0561 cx.step = cx.default_step; 0562 else 0563 cx.step = cx.this_step; 0564 cx.this_step = []; %% disable for next time 0565 end 0566 if isempty(cx.this_parm) 0567 cx.parm = cx.default_parm; 0568 else 0569 cx.parm = cx.this_parm; 0570 cx.this_parm = []; %% disable for next time 0571 end 0572 end 0573 0574 %% invoke callbacks - "final" context 0575 cpf_results = struct(); %% initialize results struct 0576 for k = 1:ncb 0577 [nx, cx, done, rollback, evnts, cb_data, cpf_results] = ... 0578 cpf_callbacks(k).fcn(-cont_steps, nx, cx, px, ... 0579 done, rollback, evnts, cb_data, cpf_callbacks(k).args, cpf_results); 0580 end 0581 cpf_results.events = cx.events; %% copy eventlog to results 0582 0583 %% update final case with solution 0584 mpct = cpf_current_mpc(cb_data.mpc_base, cb_data.mpc_target, Ybus, Yf, Yt, cb_data.ref, cb_data.pv, cb_data.pq, cx.V, cx.lam, mpopt); 0585 mpct.et = etime(clock, t0); 0586 mpct.success = success; 0587 0588 %%----- output results ----- 0589 %% convert back to original bus numbering & print results 0590 n = cpf_results.iterations + 1; 0591 cpf_results.V_hat = i2e_data(mpct, cpf_results.V_hat, NaN(nb,n), 'bus', 1); 0592 cpf_results.V = i2e_data(mpct, cpf_results.V, NaN(nb,n), 'bus', 1); 0593 results = int2ext(mpct); 0594 results.cpf = cpf_results; 0595 0596 %% zero out result fields of out-of-service gens & branches 0597 if ~isempty(results.order.gen.status.off) 0598 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0599 end 0600 if ~isempty(results.order.branch.status.off) 0601 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0602 end 0603 end 0604 0605 results.cpf.done_msg = done.msg; 0606 if mpopt.verbose 0607 fprintf('CPF TERMINATION: %s\n', done.msg); 0608 end 0609 0610 if fname 0611 [fd, msg] = fopen(fname, 'at'); 0612 if fd == -1 0613 error(msg); 0614 else 0615 if mpopt.out.all == 0 0616 printpf(results, fd, mpoption(mpopt, 'out.all', -1)); 0617 else 0618 printpf(results, fd, mpopt); 0619 end 0620 fclose(fd); 0621 end 0622 end 0623 printpf(results, 1, mpopt); 0624 0625 %% save solved case 0626 if solvedcase 0627 savecase(solvedcase, results); 0628 end 0629 0630 if nargout 0631 res = results; 0632 if nargout > 1 0633 suc = success; 0634 end 0635 % else %% don't define res, so it doesn't print anything 0636 end