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 solution algorithm, output options
termination tolerances, 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. Default contains fields:
V_p - (nb x nsteps+1) complex bus voltages from
predictor steps
lam_p - (nsteps+1) row vector of lambda values from
predictor steps
V_c - (nb x nsteps+1) complex bus voltages from
corrector steps
lam_c - (nsteps+1) row vector of lambda values from
corrector steps
max_lam - maximum value of lambda in lam_c
iterations - number of continuation steps performed
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(...);
Alternatively, for compatibility with previous versions of MATPOWER,
some of the results can be returned as individual output arguments:
[baseMVA, bus, gen, branch, success, et] = runcpf(...);
Generator reactive power limits are ignored for the continuation power flow.
Examples:
results = runcpf('case9', 'case9target');
results = runcpf('case9', 'case9target', ...
mpoption('cpf.adapt_step', 1));
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 solution algorithm, output options 0022 % termination tolerances, 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. Default contains fields: 0038 % V_p - (nb x nsteps+1) complex bus voltages from 0039 % predictor steps 0040 % lam_p - (nsteps+1) row vector of lambda values from 0041 % predictor steps 0042 % V_c - (nb x nsteps+1) complex bus voltages from 0043 % corrector steps 0044 % lam_c - (nsteps+1) row vector of lambda values from 0045 % corrector steps 0046 % max_lam - maximum value of lambda in lam_c 0047 % iterations - number of continuation steps performed 0048 % SUCCESS : the success flag can additionally be returned as 0049 % a second output argument 0050 % 0051 % Calling syntax options: 0052 % results = runcpf; 0053 % results = runcpf(basecasedata, targetcasedata); 0054 % results = runcpf(basecasedata, targetcasedata, mpopt); 0055 % results = runcpf(basecasedata, targetcasedata, mpopt, fname); 0056 % results = runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase) 0057 % [results, success] = runcpf(...); 0058 % 0059 % Alternatively, for compatibility with previous versions of MATPOWER, 0060 % some of the results can be returned as individual output arguments: 0061 % 0062 % [baseMVA, bus, gen, branch, success, et] = runcpf(...); 0063 % 0064 % Generator reactive power limits are ignored for the continuation power flow. 0065 % 0066 % Examples: 0067 % results = runcpf('case9', 'case9target'); 0068 % results = runcpf('case9', 'case9target', ... 0069 % mpoption('cpf.adapt_step', 1)); 0070 % 0071 % See also MPOPTION, RUNPF. 0072 0073 % MATPOWER 0074 % $Id: runcpf.m 2411 2014-11-04 20:13:56Z ray $ 0075 % by Ray Zimmerman, PSERC Cornell, 0076 % Shrirang Abhyankar, Argonne National Laboratory, 0077 % and Alexander Flueck, IIT 0078 % Copyright (c) 1996-2013 by Power System Engineering Research Center (PSERC) 0079 % 0080 % This file is part of MATPOWER. 0081 % See http://www.pserc.cornell.edu/matpower/ for more info. 0082 % 0083 % MATPOWER is free software: you can redistribute it and/or modify 0084 % it under the terms of the GNU General Public License as published 0085 % by the Free Software Foundation, either version 3 of the License, 0086 % or (at your option) any later version. 0087 % 0088 % MATPOWER is distributed in the hope that it will be useful, 0089 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0090 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0091 % GNU General Public License for more details. 0092 % 0093 % You should have received a copy of the GNU General Public License 0094 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0095 % 0096 % Additional permission under GNU GPL version 3 section 7 0097 % 0098 % If you modify MATPOWER, or any covered work, to interface with 0099 % other modules (such as MATLAB code and MEX-files) available in a 0100 % MATLAB(R) or comparable environment containing parts covered 0101 % under other licensing terms, the licensors of MATPOWER grant 0102 % you additional permission to convey the resulting work. 0103 0104 %%----- initialize ----- 0105 %% define named indices into bus, gen, branch matrices 0106 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0107 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0108 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0109 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0110 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0111 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0112 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0113 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0114 0115 0116 %% default arguments 0117 if nargin < 5 0118 solvedcase = ''; %% don't save solved case 0119 if nargin < 4 0120 fname = ''; %% don't print results to a file 0121 if nargin < 3 0122 mpopt = mpoption; %% use default options 0123 if nargin < 2 0124 targetcasedata = 'case9target'; 0125 if nargin < 1 0126 basecasedata = 'case9'; %% default data file is 'case9.m' 0127 end 0128 end 0129 end 0130 end 0131 end 0132 0133 %% options 0134 step = mpopt.cpf.step; %% continuation step length 0135 parameterization = mpopt.cpf.parameterization; %% parameterization 0136 adapt_step = mpopt.cpf.adapt_step; %% use adaptive step size? 0137 cb_args = mpopt.cpf.user_callback_args; 0138 0139 %% set up callbacks 0140 callback_names = {'cpf_default_callback'}; 0141 if ~isempty(mpopt.cpf.user_callback) 0142 if iscell(mpopt.cpf.user_callback) 0143 callback_names = {callback_names{:}, mpopt.cpf.user_callback{:}}; 0144 else 0145 callback_names = {callback_names{:}, mpopt.cpf.user_callback}; 0146 end 0147 end 0148 callbacks = cellfun(@str2func, callback_names, 'UniformOutput', false); 0149 0150 %% read base case data 0151 mpcbase = loadcase(basecasedata); 0152 nb = size(mpcbase.bus, 1); 0153 0154 %% add zero columns to branch for flows if needed 0155 if size(mpcbase.branch,2) < QT 0156 mpcbase.branch = [ mpcbase.branch zeros(size(mpcbase.branch, 1), QT-size(mpcbase.branch,2)) ]; 0157 end 0158 0159 %% convert to internal indexing 0160 mpcbase = ext2int(mpcbase); 0161 [baseMVAb, busb, genb, branchb] = deal(mpcbase.baseMVA, mpcbase.bus, mpcbase.gen, mpcbase.branch); 0162 0163 %% get bus index lists of each type of bus 0164 [ref, pv, pq] = bustypes(busb, genb); 0165 0166 %% generator info 0167 onb = find(genb(:, GEN_STATUS) > 0); %% which generators are on? 0168 gbusb = genb(onb, GEN_BUS); %% what buses are they at? 0169 0170 %% read target case data 0171 mpctarget = loadcase(targetcasedata); 0172 0173 %% add zero columns to branch for flows if needed 0174 if size(mpctarget.branch,2) < QT 0175 mpctarget.branch = [ mpctarget.branch zeros(size(mpctarget.branch, 1), QT-size(mpctarget.branch,2)) ]; 0176 end 0177 0178 %% convert to internal indexing 0179 mpctarget = ext2int(mpctarget); 0180 [baseMVAt, bust, gent, brancht] = deal(mpctarget.baseMVA, mpctarget.bus, mpctarget.gen, mpctarget.branch); 0181 0182 %% get bus index lists of each type of bus 0183 %[ref, pv, pq] = bustypes(bust, gent); 0184 0185 %% generator info 0186 ont = find(gent(:, GEN_STATUS) > 0); %% which generators are on? 0187 gbust = gent(ont, GEN_BUS); %% what buses are they at? 0188 0189 %%----- run the power flow ----- 0190 t0 = clock; 0191 if mpopt.verbose > 0 0192 v = mpver('all'); 0193 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); 0194 fprintf(' -- AC Continuation Power Flow\n'); 0195 end 0196 %% initial state 0197 %V0 = ones(size(bus, 1), 1); %% flat start 0198 V0 = busb(:, VM) .* exp(sqrt(-1) * pi/180 * busb(:, VA)); 0199 vcb = ones(size(V0)); %% create mask of voltage-controlled buses 0200 vcb(pq) = 0; %% exclude PQ buses 0201 k = find(vcb(gbusb)); %% in-service gens at v-c buses 0202 V0(gbusb(k)) = genb(onb(k), VG) ./ abs(V0(gbusb(k))).* V0(gbusb(k)); 0203 0204 %% build admittance matrices 0205 [Ybus, Yf, Yt] = makeYbus(baseMVAb, busb, branchb); 0206 0207 %% compute base case complex bus power injections (generation - load) 0208 Sbusb = makeSbus(baseMVAb, busb, genb); 0209 %% compute target case complex bus power injections (generation - load) 0210 Sbust = makeSbus(baseMVAt, bust, gent); 0211 0212 %% scheduled transfer 0213 Sxfr = Sbust - Sbusb; 0214 0215 %% Run the base case power flow solution 0216 if mpopt.verbose > 2 0217 mpopt_pf = mpoption(mpopt, 'verbose', max(0, mpopt.verbose-1)); 0218 else 0219 mpopt_pf = mpoption(mpopt, 'verbose', max(0, mpopt.verbose-2)); 0220 end 0221 lam = 0; 0222 [V, success, iterations] = newtonpf(Ybus, Sbusb, V0, ref, pv, pq, mpopt_pf); 0223 if mpopt.verbose > 2 0224 fprintf('step %3d : lambda = %6.3f\n', 0, 0); 0225 elseif mpopt.verbose > 1 0226 fprintf('step %3d : lambda = %6.3f, %2d Newton steps\n', 0, 0, iterations); 0227 end 0228 0229 lamprv = lam; %% lam at previous step 0230 Vprv = V; %% V at previous step 0231 continuation = 1; 0232 cont_steps = 0; 0233 0234 %% input args for callbacks 0235 cb_data = struct( ... 0236 'mpc_base', mpcbase, ... 0237 'mpc_target', mpctarget, ... 0238 'Sxfr', Sxfr, ... 0239 'Ybus', Ybus, ... 0240 'Yf', Yf, ... 0241 'Yt', Yt, ... 0242 'ref', ref, ... 0243 'pv', pv, ... 0244 'pq', pq, ... 0245 'mpopt', mpopt ); 0246 cb_state = struct(); 0247 0248 %% invoke callbacks 0249 for k = 1:length(callbacks) 0250 cb_state = callbacks{k}(cont_steps, V, lam, V, lam, ... 0251 cb_data, cb_state, cb_args); 0252 end 0253 0254 %% tangent predictor z = [dx;dlam] 0255 z = zeros(2*length(V)+1,1); 0256 z(end,1) = 1.0; 0257 while(continuation) 0258 cont_steps = cont_steps + 1; 0259 %% prediction for next step 0260 [V0, lam0, z] = cpf_predictor(V, lam, Ybus, Sxfr, pv, pq, step, z, ... 0261 Vprv, lamprv, parameterization); 0262 0263 %% save previous voltage, lambda before updating 0264 Vprv = V; 0265 lamprv = lam; 0266 0267 %% correction 0268 [V, success, i, lam] = cpf_corrector(Ybus, Sbusb, V0, ref, pv, pq, ... 0269 lam0, Sxfr, Vprv, lamprv, z, step, parameterization, mpopt_pf); 0270 if ~success 0271 continuation = 0; 0272 if mpopt.verbose 0273 fprintf('step %3d : lambda = %6.3f, corrector did not converge in %d iterations\n', cont_steps, lam, i); 0274 end 0275 break; 0276 end 0277 if mpopt.verbose > 2 0278 fprintf('step %3d : lambda = %6.3f\n', cont_steps, lam); 0279 elseif mpopt.verbose > 1 0280 fprintf('step %3d : lambda = %6.3f, %2d corrector Newton steps\n', cont_steps, lam, i); 0281 end 0282 0283 %% invoke callbacks 0284 for k = 1:length(callbacks) 0285 cb_state = callbacks{k}(cont_steps, V, lam, V0, lam0, ... 0286 cb_data, cb_state, cb_args); 0287 end 0288 0289 if ischar(mpopt.cpf.stop_at) 0290 if strcmp(upper(mpopt.cpf.stop_at), 'FULL') 0291 if abs(lam) < 1e-8 %% traced the full continuation curve 0292 if mpopt.verbose 0293 fprintf('\nTraced full continuation curve in %d continuation steps\n',cont_steps); 0294 end 0295 continuation = 0; 0296 elseif lam < lamprv && lam - step < 0 %% next step will overshoot 0297 step = lam; %% modify step-size 0298 parameterization = 1; %% change to natural parameterization 0299 adapt_step = 0; %% disable step-adaptivity 0300 end 0301 else %% == 'NOSE' 0302 if lam < lamprv %% reached the nose point 0303 if mpopt.verbose 0304 fprintf('\nReached steady state loading limit in %d continuation steps\n',cont_steps); 0305 end 0306 continuation = 0; 0307 end 0308 end 0309 else 0310 if lam < lamprv %% reached the nose point 0311 if mpopt.verbose 0312 fprintf('\nReached steady state loading limit in %d continuation steps\n', cont_steps); 0313 end 0314 continuation = 0; 0315 elseif abs(mpopt.cpf.stop_at - lam) < 1e-8 %% reached desired lambda 0316 if mpopt.verbose 0317 fprintf('\nReached desired lambda %3.2f in %d continuation steps\n', ... 0318 mpopt.cpf.stop_at, cont_steps); 0319 end 0320 continuation = 0; 0321 elseif lam + step > mpopt.cpf.stop_at %% will reach desired lambda in next step 0322 step = mpopt.cpf.stop_at - lam; %% modify step-size 0323 parameterization = 1; %% change to natural parameterization 0324 adapt_step = 0; %% disable step-adaptivity 0325 end 0326 end 0327 0328 if adapt_step && continuation 0329 %% Adapt stepsize 0330 cpf_error = norm([angle(V(pq));abs(V([pv;pq]));lam] - [angle(V0(pq));abs(V0([pv;pq]));lam0],inf); 0331 if cpf_error < mpopt.cpf.error_tol 0332 %% Increase stepsize 0333 step = step*mpopt.cpf.error_tol/cpf_error; 0334 if step > mpopt.cpf.step_max 0335 step = mpopt.cpf.step_max; 0336 end 0337 else 0338 %% decrese stepsize 0339 step = step*mpopt.cpf.error_tol/cpf_error; 0340 if step < mpopt.cpf.step_min 0341 step = mpopt.cpf.step_min; 0342 end 0343 end 0344 end 0345 end 0346 0347 %% invoke callbacks 0348 if success 0349 cpf_results = struct(); 0350 for k = 1:length(callbacks) 0351 [cb_state, cpf_results] = callbacks{k}(cont_steps, V, lam, V0, lam0, ... 0352 cb_data, cb_state, cb_args, cpf_results); 0353 end 0354 else 0355 cpf_results.iterations = i; 0356 end 0357 0358 %% update bus and gen matrices to reflect the loading and generation 0359 %% at the noise point 0360 bust(:,PD) = busb(:,PD) + lam*(bust(:,PD) - busb(:,PD)); 0361 bust(:,QD) = busb(:,QD) + lam*(bust(:,QD) - busb(:,QD)); 0362 gent(:,PG) = genb(:,PG) + lam*(gent(:,PG) - genb(:,PG)); 0363 0364 %% update data matrices with solution 0365 [bust, gent, brancht] = pfsoln(baseMVAt, bust, gent, brancht, Ybus, Yf, Yt, V, ref, pv, pq); 0366 0367 mpctarget.et = etime(clock, t0); 0368 mpctarget.success = success; 0369 0370 %%----- output results ----- 0371 %% convert back to original bus numbering & print results 0372 [mpctarget.bus, mpctarget.gen, mpctarget.branch] = deal(bust, gent, brancht); 0373 if success 0374 n = cpf_results.iterations + 1; 0375 cpf_results.V_p = i2e_data(mpctarget, cpf_results.V_p, NaN(nb,n), 'bus', 1); 0376 cpf_results.V_c = i2e_data(mpctarget, cpf_results.V_c, NaN(nb,n), 'bus', 1); 0377 end 0378 results = int2ext(mpctarget); 0379 results.cpf = cpf_results; 0380 0381 %% zero out result fields of out-of-service gens & branches 0382 if ~isempty(results.order.gen.status.off) 0383 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0384 end 0385 if ~isempty(results.order.branch.status.off) 0386 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0387 end 0388 0389 if fname 0390 [fd, msg] = fopen(fname, 'at'); 0391 if fd == -1 0392 error(msg); 0393 else 0394 if mpopt.out.all == 0 0395 printpf(results, fd, mpoption(mpopt, 'out.all', -1)); 0396 else 0397 printpf(results, fd, mpopt); 0398 end 0399 fclose(fd); 0400 end 0401 end 0402 printpf(results, 1, mpopt); 0403 0404 %% save solved case 0405 if solvedcase 0406 savecase(solvedcase, results); 0407 end 0408 0409 if nargout 0410 res = results; 0411 if nargout > 1 0412 suc = success; 0413 end 0414 % else %% don't define res, so it doesn't print anything 0415 end