Home > matpower6.0 > runcpf.m

runcpf

PURPOSE ^

RUNCPF Runs a full AC continuation power flow

SYNOPSIS ^

function [res, suc] =runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase)

DESCRIPTION ^

RUNCPF  Runs a full AC continuation power flow
   [RESULTS, SUCCESS] = RUNCPF(BASECASEDATA, TARGETCASEDATA, ...
                                   MPOPT, FNAME, SOLVEDCASE)

   Runs a full AC continuation power flow using a normalized tangent
   predictor and selected parameterization scheme, optionally
   returning a RESULTS struct and SUCCESS flag. Step size can be
   fixed or adaptive.

   Inputs (all are optional):
       BASECASEDATA : either a MATPOWER case struct or a string containing
           the name of the file with the case data defining the base loading
           and generation (default is 'case9')
           (see also CASEFORMAT and LOADCASE)
       TARGETCASEDATA : either a MATPOWER case struct or a string
           containing the name of the file with the case data defining the
           target loading and generation (default is 'case9target')
       MPOPT : MATPOWER options struct to override default options
           can be used to specify the parameterization, output options,
           termination criteria, and more (see also MPOPTION).
       FNAME : name of a file to which the pretty-printed output will
           be appended
       SOLVEDCASE : name of file to which the solved case will be saved
           in MATPOWER case format (M-file will be assumed unless the
           specified name ends with '.mat')

   Outputs (all are optional):
       RESULTS : results struct, with the following fields:
           (all fields from the input MATPOWER case, i.e. bus, branch,
               gen, etc., but with solved voltages, power flows, etc.)
           order - info used in external <-> internal data conversion
           et - elapsed time in seconds
           success - success flag, 1 = succeeded, 0 = failed
           cpf - CPF output struct whose content depends on any
                   user callback functions, where default contains fields:
               done_msg - string with message describing cause of
                       continuation termination
               iterations - number of continuation steps performed
               lam - (nsteps+1) row vector of lambda values from
                       correction steps
               lam_hat - (nsteps+1) row vector of lambda values from
                       prediction steps
               max_lam - maximum value of lambda in RESULTS.cpf.lam
               steps - (nsteps+1) row vector of stepsizes taken at each
                       continuation step
               V - (nb x nsteps+1) complex bus voltages from
                       correction steps
               V_hat - (nb x nsteps+1) complex bus voltages from
                       prediction steps
               events - an array of structs of size nevents with the
                   following fields:
                    k - continuation step number at which event was located
                    name - name of event
                    idx - index(es) of critical elements in corresponding
                       event function, e.g. index of generator reaching VAr
                       limit
                    msg - descriptive text detailing the event
       SUCCESS : the success flag can additionally be returned as
           a second output argument

   Calling syntax options:
       results = runcpf;
       results = runcpf(basecasedata, targetcasedata);
       results = runcpf(basecasedata, targetcasedata, mpopt);
       results = runcpf(basecasedata, targetcasedata, mpopt, fname);
       results = runcpf(basecasedata, targetcasedata, mpopt, fname, solvedcase)
       [results, success] = runcpf(...);

   If the 'cpf.enforce_q_lims' option is set to true (default is false) then,
   if any generator reaches its reactive power limits during the AC
   continuation power flow,
       - the corresponding bus is converted to a PQ bus, and the problem
         is modified to eliminate further reactive transfer on this bus
       - the voltage magnitude at the bus will deviate from the specified
         setpoint to satisfy the reactive power limit,
       - if the reference bus is converted to PQ, further real power transfer
         for the bus is also eliminated, and the first remaining PV bus is
         selected as the new slack, resulting in the transfers at both
         reference buses potentially deviating from the specified values
       - if all reference and PV buses are converted to PQ, RUNCPF terminates
         with an infeasibility message.

   If the 'cpf.enforce_p_lims' option is set to true (default is fals) then,
   if any generator reaches its maximum active power limit during the AC
   continuation power flow,
       - the problem is modified to eliminate further active transfer by
         this generator
       - if the generator was at the reference bus, it is converted to PV
         and the first remaining PV bus is selected as the new slack.

   Possible CPF termination modes:
       when cpf.stop_at == 'NOSE'
           - Reached steady state loading limit
           - Nose point eliminated by limit induced bifurcation
       when cpf.stop_at == 'FULL'
           - Traced full continuation curve
       when cpf.stop_at == <target_lam_val>
           - Reached desired lambda
       when cpf.enforce_p_lims == true
           - All generators at PMAX
       when cpf.enforce_q_lims == true
           - No REF or PV buses remaining
       other
           - Base case power flow did not converge
           - Base and target case have identical load and generation
           - Corrector did not converge
           - Could not locate <event_name> event
           - Too many rollback steps triggered by callbacks

   Examples:
       results = runcpf('case9', 'case9target');
       results = runcpf('case9', 'case9target', ...
                           mpoption('cpf.adapt_step', 1));
       results = runcpf('case9', 'case9target', ...
                           mpoption('cpf.enforce_q_lims', 1));
       results = runcpf('case9', 'case9target', ...
                           mpoption('cpf.stop_at', 'FULL'));

   See also MPOPTION, RUNPF.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005