Description of runcpf

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
