Home > matpower7.1 > lib > 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.

   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.

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 %   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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005