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