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.
If the 'cpf.enforce_v_lims' option is set to true (default is false)
then the continuation power flow is set to terminate if any bus voltage
magnitude exceeds its minimum or maximum limit.
If the 'cpf.enforce_flow_lims' option is set to true (default is false)
then the continuation power flow is set to terminate if any line MVA
flow exceeds its rateA limit.
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
when cpf.enforce_v_lims == true
- Any bus voltage magnitude limit is reached
when cpf.enforce_flow_lims == true
- Any branch MVA flow limit is reached
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 % If the 'cpf.enforce_v_lims' option is set to true (default is false) 0094 % then the continuation power flow is set to terminate if any bus voltage 0095 % magnitude exceeds its minimum or maximum limit. 0096 % 0097 % If the 'cpf.enforce_flow_lims' option is set to true (default is false) 0098 % then the continuation power flow is set to terminate if any line MVA 0099 % flow exceeds its rateA limit. 0100 % 0101 % Possible CPF termination modes: 0102 % when cpf.stop_at == 'NOSE' 0103 % - Reached steady state loading limit 0104 % - Nose point eliminated by limit induced bifurcation 0105 % when cpf.stop_at == 'FULL' 0106 % - Traced full continuation curve 0107 % when cpf.stop_at == <target_lam_val> 0108 % - Reached desired lambda 0109 % when cpf.enforce_p_lims == true 0110 % - All generators at PMAX 0111 % when cpf.enforce_q_lims == true 0112 % - No REF or PV buses remaining 0113 % when cpf.enforce_v_lims == true 0114 % - Any bus voltage magnitude limit is reached 0115 % when cpf.enforce_flow_lims == true 0116 % - Any branch MVA flow limit is reached 0117 % other 0118 % - Base case power flow did not converge 0119 % - Base and target case have identical load and generation 0120 % - Corrector did not converge 0121 % - Could not locate <event_name> event 0122 % - Too many rollback steps triggered by callbacks 0123 % 0124 % Examples: 0125 % results = runcpf('case9', 'case9target'); 0126 % results = runcpf('case9', 'case9target', ... 0127 % mpoption('cpf.adapt_step', 1)); 0128 % results = runcpf('case9', 'case9target', ... 0129 % mpoption('cpf.enforce_q_lims', 1)); 0130 % results = runcpf('case9', 'case9target', ... 0131 % mpoption('cpf.stop_at', 'FULL')); 0132 % 0133 % See also MPOPTION, RUNPF. 0134 0135 % MATPOWER 0136 % Copyright (c) 1996-2017, Power Systems Engineering Research Center (PSERC) 0137 % by Ray Zimmerman, PSERC Cornell, 0138 % Shrirang Abhyankar, Argonne National Laboratory, 0139 % and Alexander Flueck, IIT 0140 % 0141 % This file is part of MATPOWER. 0142 % Covered by the 3-clause BSD License (see LICENSE file for details). 0143 % See https://matpower.org for more info. 0144 0145 %%----- initialize ----- 0146 %% define named indices into bus, gen, branch matrices 0147 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0148 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0149 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0150 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0151 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0152 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0153 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0154 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0155 0156 %% default arguments 0157 if nargin < 5 0158 solvedcase = ''; %% don't save solved case 0159 if nargin < 4 0160 fname = ''; %% don't print results to a file 0161 if nargin < 3 0162 mpopt = mpoption; %% use default options 0163 if nargin < 2 0164 targetcasedata = 'case9target'; 0165 if nargin < 1 0166 basecasedata = 'case9'; %% default data file is 'case9.m' 0167 end 0168 end 0169 end 0170 end 0171 end 0172 0173 %% options 0174 step = mpopt.cpf.step; %% continuation step length 0175 parm = mpopt.cpf.parameterization; %% parameterization 0176 adapt_step = mpopt.cpf.adapt_step; %% use adaptive step size? 0177 qlim = mpopt.cpf.enforce_q_lims; %% enforce reactive limits 0178 plim = mpopt.cpf.enforce_p_lims; %% enforce active limits 0179 vlim = mpopt.cpf.enforce_v_lims; %% enforce voltage magnitude limits 0180 flim = mpopt.cpf.enforce_flow_lims; %% enforce branch flow limits 0181 0182 %% register event functions (for event detection) 0183 %% and CPF callback functions (for event handling and other tasks) 0184 cpf_events = []; 0185 cpf_callbacks = []; 0186 %% to handle CPF termination 0187 if ischar(mpopt.cpf.stop_at) && strcmp(mpopt.cpf.stop_at, 'NOSE'); 0188 cpf_events = cpf_register_event(cpf_events, 'NOSE', 'cpf_nose_event', mpopt.cpf.nose_tol, 1); 0189 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_nose_event_cb', 51); 0190 else %% FULL or target lambda 0191 cpf_events = cpf_register_event(cpf_events, 'TARGET_LAM', 'cpf_target_lam_event', mpopt.cpf.target_lam_tol, 1); 0192 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_target_lam_event_cb', 50); 0193 end 0194 %% to handle branch flow limits 0195 if flim 0196 cpf_events = cpf_register_event(cpf_events, 'FLIM', 'cpf_flim_event', mpopt.cpf.flow_lims_tol, 1); 0197 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_flim_event_cb', 53); 0198 end 0199 %% to handle voltage limits 0200 if vlim 0201 cpf_events = cpf_register_event(cpf_events, 'VLIM', 'cpf_vlim_event', mpopt.cpf.v_lims_tol, 1); 0202 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_vlim_event_cb', 52); 0203 end 0204 %% to handle reactive power limits 0205 if qlim 0206 cpf_events = cpf_register_event(cpf_events, 'QLIM', 'cpf_qlim_event', mpopt.cpf.q_lims_tol, 1); 0207 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_qlim_event_cb', 41); 0208 end 0209 %% to handle active power limits 0210 if plim 0211 cpf_events = cpf_register_event(cpf_events, 'PLIM', 'cpf_plim_event', mpopt.cpf.p_lims_tol, 1); 0212 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_plim_event_cb', 40); 0213 end 0214 cpf_callbacks = cpf_register_callback(cpf_callbacks, 'cpf_default_callback', 0); 0215 %% user callbacks 0216 if ~isempty(mpopt.cpf.user_callback) 0217 %% convert to cell array, if necessary 0218 if iscell(mpopt.cpf.user_callback) 0219 user_callbacks = mpopt.cpf.user_callback; 0220 else 0221 user_callbacks = {mpopt.cpf.user_callback}; 0222 end 0223 for k = 1:length(user_callbacks) 0224 %% convert each to struct, if necessary 0225 if isstruct(user_callbacks{k}) 0226 ucb = user_callbacks{k}; 0227 else 0228 ucb = struct('fcn', user_callbacks{k}); 0229 end 0230 %% set default priority, args if not provided 0231 if ~isfield(ucb, 'priority') 0232 ucb.priority = []; 0233 end 0234 if ~isfield(ucb, 'args') 0235 ucb.args = []; 0236 end 0237 %% register the user callback 0238 cpf_callbacks = cpf_register_callback(cpf_callbacks, ... 0239 ucb.fcn, ucb.priority, ucb.args); 0240 end 0241 end 0242 nef = length(cpf_events); %% number of event functions registered 0243 ncb = length(cpf_callbacks); %% number of callback functions registered 0244 0245 %% set power flow options 0246 if mpopt.verbose > 4 0247 mpopt_pf = mpoption(mpopt, 'verbose', 2); 0248 else 0249 mpopt_pf = mpoption(mpopt, 'verbose', 0); 0250 end 0251 mpopt_pf = mpoption(mpopt_pf, 'pf.enforce_q_lims', qlim); 0252 0253 %% load base case data 0254 mpcb = loadcase(basecasedata); 0255 0256 %% clip base active generator outputs to PMAX, if necessary 0257 idx_pmax = []; %% indices of generators clipped at PMAX 0258 if plim 0259 idx_pmax = find( mpcb.gen(:, GEN_STATUS) > 0 & ... 0260 mpcb.gen(:, PG) - mpcb.gen(:, PMAX) > -mpopt.cpf.p_lims_tol); 0261 if mpopt.verbose && ~isempty(idx_pmax) 0262 fprintf('base case real power output of gen %d reduced from %g to %g MW (PMAX)\n', ... 0263 [idx_pmax mpcb.gen(idx_pmax, PG) mpcb.gen(idx_pmax, PMAX)]'); 0264 end 0265 mpcb.gen(idx_pmax, PG) = mpcb.gen(idx_pmax, PMAX); 0266 end 0267 0268 %% run base case power flow 0269 [rb, success] = runpf(mpcb, mpopt_pf); 0270 if success 0271 done = struct('flag', 0, 'msg', ''); 0272 if qlim 0273 %% find buses that were converted to PQ or REF by initial power flow 0274 b2ref = rb.bus(:, BUS_TYPE) == REF & mpcb.bus(:, BUS_TYPE) ~= REF; 0275 b2pq = rb.bus(:, BUS_TYPE) == PQ & mpcb.bus(:, BUS_TYPE) ~= PQ; 0276 end 0277 mpcb = rb; %% update base case with solved power flow 0278 else 0279 done = struct('flag', 1, 'msg', 'Base case power flow did not converge.'); 0280 results = rb; 0281 results.cpf = struct(); 0282 end 0283 0284 if ~done.flag 0285 %% read target case data 0286 mpct = loadcase(targetcasedata); 0287 if size(mpct.branch,2) < QT %% add zero columns to branch for flows if needed 0288 mpct.branch = [ mpct.branch zeros(size(mpct.branch, 1), QT-size(mpct.branch,2)) ]; 0289 end 0290 0291 %% convert both to internal indexing 0292 mpcb = ext2int(mpcb, mpopt); 0293 mpct = ext2int(mpct, mpopt); 0294 i2e_gen = mpcb.order.gen.i2e; 0295 if qlim 0296 b2ref = e2i_data(mpcb, b2ref, 'bus'); 0297 b2pq = e2i_data(mpcb, b2pq, 'bus'); 0298 end 0299 nb = size(mpcb.bus, 1); 0300 0301 %% get bus index lists of each type of bus 0302 [ref, pv, pq] = bustypes(mpcb.bus, mpcb.gen); 0303 0304 %% generator info 0305 %% find generators that are ON and at voltage-controlled buses 0306 ong = find(mpcb.gen(:, GEN_STATUS) > 0 ... 0307 & mpcb.bus(mpcb.gen(:, GEN_BUS), BUS_TYPE) ~= PQ); 0308 gbus = mpcb.gen(ong, GEN_BUS); %% what buses are they at? 0309 0310 %% make sure target case is same as base case w.r.t 0311 %% bus types, GEN_STATUS, Qg and slack Pg 0312 if qlim 0313 bb = find(b2ref | b2pq); 0314 mpct.bus(bb, BUS_TYPE) = mpcb.bus(bb, BUS_TYPE); 0315 end 0316 if any(mpcb.bus(:, BUS_TYPE) ~= mpct.bus(:, BUS_TYPE)) 0317 error('runcpf: BUS_TYPE of all buses must be the same in base and target cases'); 0318 end 0319 if any(mpcb.gen(:, GEN_STATUS) ~= mpct.gen(:, GEN_STATUS)) 0320 error('runcpf: GEN_STATUS of all generators must be the same in base and target cases'); 0321 end 0322 mpct.gen(ong, QG) = mpcb.gen(ong, QG); 0323 if qlim 0324 %% find generators that are ON and at buses that were 0325 %% converted to PQ by initial power flow 0326 g2pq = find(mpcb.gen(:, GEN_STATUS) > 0 & b2pq(mpcb.gen(:, GEN_BUS))); 0327 mpct.gen(g2pq, QG) = mpcb.gen(g2pq, QG); 0328 end 0329 for k = 1:length(ref) 0330 refgen = find(gbus == ref(k)); 0331 mpct.gen(ong(refgen), PG) = mpcb.gen(ong(refgen), PG); 0332 end 0333 0334 %% zero transfers for gens that exceed PMAX limits, if necessary 0335 if plim 0336 idx_pmax = find( mpcb.gen(:, GEN_STATUS) > 0 & ... 0337 mpcb.gen(:, PG) - mpcb.gen(:, PMAX) > -mpopt.cpf.p_lims_tol & ... 0338 mpct.gen(:, PG) - mpct.gen(:, PMAX) > -mpopt.cpf.p_lims_tol); 0339 if ~isempty(idx_pmax) 0340 if mpopt.verbose 0341 fprintf('target case real power output of gen %d reduced from %g to %g MW (PMAX)\n', ... 0342 [i2e_gen(idx_pmax) mpct.gen(idx_pmax, PG) mpcb.gen(idx_pmax, PG)]'); 0343 end 0344 mpct.gen(idx_pmax, PG) = mpcb.gen(idx_pmax, PG); 0345 end 0346 end 0347 0348 %%----- run the continuation power flow ----- 0349 t0 = tic; 0350 if mpopt.verbose 0351 v = mpver('all'); 0352 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); 0353 fprintf(' -- AC Continuation Power Flow\n'); 0354 if mpopt.verbose > 1 0355 fprintf('step %3d : lambda = %6.3f, %2d Newton steps\n', 0, 0, mpcb.iterations); 0356 end 0357 end 0358 0359 %% build admittance matrices 0360 [Ybus, Yf, Yt] = makeYbus(mpcb.baseMVA, mpcb.bus, mpcb.branch); 0361 0362 %% functions for computing base and target case V-dependent complex bus 0363 %% power injections: (generation - load) 0364 Sbusb = @(Vm)makeSbus(mpcb.baseMVA, mpcb.bus, mpcb.gen, mpopt, Vm); 0365 Sbust = @(Vm)makeSbus(mpct.baseMVA, mpct.bus, mpct.gen, mpopt, Vm); 0366 0367 %% initialize variables 0368 cont_steps = 0; 0369 lam = 0; 0370 V = mpcb.bus(:, VM) .* exp(1j * pi/180 * mpcb.bus(:, VA)); 0371 rollback = 0; %% flag to indicate that a step must be rolled back 0372 locating = 0; %% flag to indicate that an event interval was detected, 0373 %% but the event has not yet been located 0374 rb_cnt_ef = 0; %% counter for rollback steps triggered by event function intervals 0375 rb_cnt_cb = 0; %% counter for rollback steps triggered directly by callbacks 0376 0377 %% initialize tangent predictor: z = [dx;dlam] 0378 z = [zeros(2*nb, 1); 1]; 0379 direction = 1; 0380 z = cpf_tangent(V, lam, Ybus, Sbusb, Sbust, pv, pq, ... 0381 z, V, lam, parm, direction); 0382 0383 %% initialize state for current continuation step 0384 cx = struct(... %% current state 0385 'lam_hat', lam, ... %% predicted lambda 0386 'V_hat', V, ... %% predicted V 0387 'lam', lam, ... %% corrected lambda 0388 'V', V, ... %% corrected V 0389 'z', z, ... %% normalized tangent predictor 0390 'default_step', step, ... %% default step size 0391 'default_parm', parm, ... %% default parameterization 0392 'this_step', [], ... %% step size for this step only 0393 'this_parm', [], ... %% parameterization for this step only 0394 'step', step, ... %% current step size 0395 'parm', parm, ... %% current parameterization 0396 'events', [], ... %% event log 0397 'cb', struct(), ... %% user state, for callbacks (replaces cb_state) 0398 'ef', {cell(nef, 1)} ... %% event function values 0399 ); 0400 0401 %% input args for callbacks 0402 cb_data = struct( ... 0403 'mpc_base', mpcb, ... %% MATPOWER case struct - base case 0404 'mpc_target', mpct, ... %% MATPOWER case struct - target case 0405 'Sbusb', Sbusb, ... %% function for computing base bus inj 0406 'Sbust', Sbust, ... %% function for computing target bus inj 0407 'Ybus', Ybus, ... %% bus admittance matrix 0408 'Yf', Yf, ... %% branch admittance matrix, from end 0409 'Yt', Yt, ... %% branch admittance matrix, to end 0410 'ref', ref, ... %% vector of ref bus indices 0411 'pv', pv, ... %% vector of PV bus indices 0412 'pq', pq, ... %% vector of PQ bus indices 0413 'idx_pmax', idx_pmax, ... %% vector of gen indices of gens at PMAX 0414 'mpopt', mpopt ); %% MATPOWER option struct 0415 0416 %% initialize event function values 0417 for k = 1:nef 0418 cx.ef{k} = cpf_events(k).fcn(cb_data, cx); 0419 end 0420 0421 %% invoke callbacks - "initialize" context 0422 for k = 1:ncb 0423 [nx, cx, done, rollback, evnts, cb_data] = cpf_callbacks(k).fcn( ... 0424 cont_steps, cx, cx, cx, done, 0, [], cb_data, cpf_callbacks(k).args); 0425 end 0426 0427 %% check for case with no transfer 0428 if norm(Sbust(abs(cx.V)) - Sbusb(abs(cx.V))) == 0 0429 done.flag = 1; 0430 done.msg = 'Base case and target case have identical load and generation'; 0431 end 0432 0433 cont_steps = cont_steps + 1; 0434 px = cx; %% initialize state for previous continuation step 0435 while ~done.flag 0436 %% initialize next candidate with current state 0437 nx = cx; 0438 0439 %% prediction for next step 0440 [nx.V_hat, nx.lam_hat] = cpf_predictor(cx.V, cx.lam, cx.z, cx.step, cb_data.pv, cb_data.pq); 0441 0442 %% correction 0443 [nx.V, success, i, nx.lam] = cpf_corrector(Ybus, cb_data.Sbusb, nx.V_hat, cb_data.ref, cb_data.pv, cb_data.pq, ... 0444 nx.lam_hat, cb_data.Sbust, cx.V, cx.lam, cx.z, cx.step, cx.parm, mpopt_pf); 0445 if ~success %% corrector failed 0446 done.flag = 1; 0447 done.msg = sprintf('Corrector did not converge in %d iterations.', i); 0448 if mpopt.verbose 0449 fprintf('step %3d : stepsize = %-9.3g lambda = %6.3f corrector did not converge in %d iterations\n', cont_steps, cx.step, nx.lam, i); 0450 end 0451 cont_steps = max(cont_steps - 1, 1); %% go back to last step, but not to 0 0452 break; 0453 end 0454 0455 %% compute new tangent direction, based on a previous state: tx 0456 if nx.step == 0 %% if this is a re-do step, cx and nx are the same 0457 tx = px; %% so use px as the previous state 0458 else %% otherwise 0459 tx = cx; %% use cx as the previous state 0460 end 0461 nx.z = cpf_tangent(nx.V, nx.lam, Ybus, cb_data.Sbusb, cb_data.Sbust, ... 0462 cb_data.pv, cb_data.pq, ... 0463 cx.z, tx.V, tx.lam, nx.parm, direction); 0464 0465 %% detect events 0466 for k = 1:nef 0467 nx.ef{k} = cpf_events(k).fcn(cb_data, nx); %% update event functions 0468 end 0469 [rollback, evnts, nx.ef] = cpf_detect_events(cpf_events, nx.ef, cx.ef, nx.step, mpopt.verbose); 0470 0471 %% adjust step-size to locate event function zero, if necessary 0472 if rollback %% current step overshot 0473 %% rollback and initialize next step size based on rollback and previous 0474 rx = nx; %% save state we're rolling back from 0475 rx.evnts = evnts; %% and critical event info 0476 cx.this_step = evnts.step_scale * rx.step; 0477 cx.this_parm = rx.parm; %% keep same parameterization as last step 0478 locating = 1; %% enter "locating" mode (or stay in it) 0479 rb_cnt_ef = rb_cnt_ef + 1; %% increment rollback counter for ef intervals 0480 if rb_cnt_ef > 26 0481 done.flag = 1; 0482 done.msg = sprintf('Could not locate %s event!', evnts.name); 0483 end 0484 if mpopt.verbose > 3 0485 loc_msg = sprintf('OVERSHOOT : f = [%g, <<%g>>], step <-- %.4g', ... 0486 cx.ef{evnts.eidx}(evnts.idx(1)), ... 0487 rx.ef{evnts.eidx}(evnts.idx(1)), cx.this_step); 0488 end 0489 elseif locating 0490 if evnts(1).zero %% found the zero! 0491 %% step size will be reset to previously used default step size 0492 locating = 0; %% exit "locating" mode 0493 rb_cnt_ef = 0; %% reset rollback counter for ef intervals 0494 if mpopt.verbose > 3 0495 loc_msg = sprintf('ZERO! : f = %g, step <-- %.4g', ... 0496 nx.ef{rx.evnts.eidx}(rx.evnts.idx(1)), nx.default_step); 0497 end 0498 else %% prev rollback undershot 0499 %% initialize next step size based on critical event function 0500 %% values from prev rollback step and current step 0501 rx_ef = rx.ef{rx.evnts.eidx}(rx.evnts.idx(1)); 0502 cx_ef = nx.ef{rx.evnts.eidx}(rx.evnts.idx(1)); 0503 step_scale = cx_ef / (cx_ef - rx_ef); 0504 nx.this_step = step_scale * (rx.step - nx.step); 0505 rb_cnt_ef = 0; %% reset rollback counter for ef intervals 0506 if mpopt.verbose > 3 0507 loc_msg = sprintf('UNDERSHOOT : f [<<%g>>, %g], step <-- %.4g', ... 0508 cx_ef, rx_ef, nx.this_step); 0509 end 0510 end 0511 else 0512 loc_msg = ''; 0513 direction = 1; 0514 end 0515 0516 %% invoke callbacks - "iterations" context 0517 rb = rollback; 0518 for k = 1:ncb 0519 [nx, cx, done, rollback, evnts, cb_data] = cpf_callbacks(k).fcn( ... 0520 cont_steps, nx, cx, px, done, rollback, evnts, cb_data, cpf_callbacks(k).args); 0521 end 0522 if ~rb && rollback %% rollback triggered by callback (vs event function interval) 0523 rb_cnt_cb = rb_cnt_cb + 1; %% increment rollback counter for callbacks 0524 if rb_cnt_cb > 26 0525 done.flag = 1; 0526 done.msg = 'Too many rollback steps triggered by callbacks!'; 0527 end 0528 else 0529 if ~done.flag && evnts(1).zero 0530 mpce = cpf_current_mpc(cb_data.mpc_base, cb_data.mpc_target, Ybus, Yf, Yt, cb_data.ref, cb_data.pv, cb_data.pq, nx.V, nx.lam, mpopt); 0531 J = makeJac(mpce); 0532 opts.tol = 1e-3; 0533 opts.it = 2*nb; 0534 %% retain original direction for buses with negative V-Q 0535 %% sensitivity, otherwise change direction based on tangent 0536 %% direction and manifold (eigenvalue) 0537 if isempty(find(nx.z(nb+pq) > 0, 1)) 0538 direction = sign(nx.z(end)*min(real(eigs(J,1,'SR',opts)))); 0539 end 0540 end 0541 rb_cnt_cb = 0; %% reset rollback counter for callbacks 0542 end 0543 0544 %% print iteration information 0545 if mpopt.verbose > 1 0546 %% set label for rollback step counter 0547 if rb_cnt_ef 0548 sub_step = char('a' + rb_cnt_ef - 1); 0549 elseif rb_cnt_cb 0550 sub_step = char('A' + rb_cnt_cb - 1); 0551 else 0552 sub_step = ' '; 0553 end 0554 if mpopt.verbose > 4 0555 fprintf('step %3d%s : stepsize = %-9.3g lambda = %6.3f', cont_steps, sub_step, cx.step, nx.lam); 0556 else 0557 fprintf('step %3d%s : stepsize = %-9.3g lambda = %6.3f %2d corrector Newton steps', cont_steps, sub_step, cx.step, nx.lam, i); 0558 end 0559 if rollback 0560 fprintf(' ^ ROLLBACK\n'); 0561 else 0562 fprintf('\n'); 0563 end 0564 if mpopt.verbose > 3 && ~isempty(loc_msg) 0565 fprintf(' LOCATING -- %s\n', loc_msg); 0566 end 0567 end 0568 0569 %% log events 0570 for k = 1:length(evnts) 0571 if evnts(k).log 0572 e = struct( 'k', cont_steps, ... 0573 'name', evnts(k).name, ... 0574 'idx', evnts(k).idx, ... 0575 'msg', evnts(k).msg ); 0576 if isempty(nx.events) 0577 nx.events = e; 0578 else 0579 nx.events(end+1) = e; 0580 end 0581 end 0582 if (mpopt.verbose > 2 && evnts(k).log) || ... 0583 (mpopt.verbose > 3 && evnts(k).eidx) 0584 fprintf(' %s\n', evnts(k).msg); 0585 end 0586 end 0587 0588 %% adapt stepsize if requested and not terminating, locating a zero 0589 %% or re-doing a step after changing the problem data 0590 if adapt_step && ~done.flag && ~locating && ~evnts(1).zero && nx.step ~= 0 0591 cpf_error = norm([angle(nx.V( cb_data.pq)); abs(nx.V( [cb_data.pv;cb_data.pq])); nx.lam] - ... 0592 [angle(nx.V_hat(cb_data.pq)); abs(nx.V_hat([cb_data.pv;cb_data.pq])); nx.lam_hat], inf); 0593 0594 %% new nominal step size is current size * tol/err, but we reduce 0595 %% the change from the current size by a damping factor and limit 0596 %% increases to a factor of 2 0597 step_scale = min(2, 1 + mpopt.cpf.adapt_step_damping * ... 0598 (mpopt.cpf.adapt_step_tol/cpf_error - 1)); 0599 nx.default_step = nx.step * step_scale; 0600 0601 %% limit step-size 0602 if nx.default_step > mpopt.cpf.step_max 0603 nx.default_step = mpopt.cpf.step_max; 0604 end 0605 if nx.default_step < mpopt.cpf.step_min 0606 nx.default_step = mpopt.cpf.step_min; 0607 end 0608 end 0609 0610 %% if this is a normal step 0611 if ~rollback 0612 px = cx; %% save current state before update 0613 cx = nx; %% update current state to next candidate 0614 if ~done.flag 0615 cont_steps = cont_steps + 1; 0616 end 0617 end 0618 0619 %% set step size and parameterization, from one-time or defaults 0620 if isempty(cx.this_step) 0621 cx.step = cx.default_step; 0622 else 0623 cx.step = cx.this_step; 0624 cx.this_step = []; %% disable for next time 0625 end 0626 if isempty(cx.this_parm) 0627 cx.parm = cx.default_parm; 0628 else 0629 cx.parm = cx.this_parm; 0630 cx.this_parm = []; %% disable for next time 0631 end 0632 end 0633 0634 %% invoke callbacks - "final" context 0635 cpf_results = struct(); %% initialize results struct 0636 for k = 1:ncb 0637 [nx, cx, done, rollback, evnts, cb_data, cpf_results] = ... 0638 cpf_callbacks(k).fcn(-cont_steps, nx, cx, px, ... 0639 done, rollback, evnts, cb_data, cpf_callbacks(k).args, cpf_results); 0640 end 0641 cpf_results.events = cx.events; %% copy eventlog to results 0642 0643 %% update final case with solution 0644 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); 0645 mpct.et = toc(t0); 0646 mpct.success = success; 0647 0648 %%----- output results ----- 0649 %% convert back to original bus numbering & print results 0650 n = size(cpf_results.V, 2); 0651 cpf_results.V_hat = i2e_data(mpct, cpf_results.V_hat, NaN(nb,n), 'bus', 1); 0652 cpf_results.V = i2e_data(mpct, cpf_results.V, NaN(nb,n), 'bus', 1); 0653 results = int2ext(mpct); 0654 results.cpf = cpf_results; 0655 0656 %% zero out result fields of out-of-service gens & branches 0657 if ~isempty(results.order.gen.status.off) 0658 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0659 end 0660 if ~isempty(results.order.branch.status.off) 0661 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0662 end 0663 end 0664 0665 results.cpf.done_msg = done.msg; 0666 if mpopt.verbose 0667 fprintf('CPF TERMINATION: %s\n', done.msg); 0668 end 0669 0670 if fname 0671 [fd, msg] = fopen(fname, 'at'); 0672 if fd == -1 0673 error(msg); 0674 else 0675 if mpopt.out.all == 0 0676 printpf(results, fd, mpoption(mpopt, 'out.all', -1)); 0677 else 0678 printpf(results, fd, mpopt); 0679 end 0680 fclose(fd); 0681 end 0682 end 0683 printpf(results, 1, mpopt); 0684 0685 %% save solved case 0686 if solvedcase 0687 savecase(solvedcase, results); 0688 end 0689 0690 if nargout 0691 res = results; 0692 if nargout > 1 0693 suc = success; 0694 end 0695 % else %% don't define res, so it doesn't print anything 0696 end