Home > matpower5.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 solution algorithm, output options
           termination tolerances, and more (see also MPOPTION).
       FNAME : name of a file to which the pretty-printed output will
           be appended
       SOLVEDCASE : name of file to which the solved case will be saved
           in MATPOWER case format (M-file will be assumed unless the
           specified name ends with '.mat')

   Outputs (all are optional):
       RESULTS : results struct, with the following fields:
           (all fields from the input MATPOWER case, i.e. bus, branch,
               gen, etc., but with solved voltages, power flows, etc.)
           order - info used in external <-> internal data conversion
           et - elapsed time in seconds
           success - success flag, 1 = succeeded, 0 = failed
           cpf - CPF output struct whose content depends on any
               user callback functions. Default contains fields:
               V_p - (nb x nsteps+1) complex bus voltages from
                       predictor steps
               lam_p - (nsteps+1) row vector of lambda values from
                       predictor steps
               V_c - (nb x nsteps+1) complex bus voltages from
                       corrector steps
               lam_c - (nsteps+1) row vector of lambda values from
                       corrector steps
               max_lam - maximum value of lambda in lam_c
               iterations - number of continuation steps performed
       SUCCESS : the success flag can additionally be returned as
           a second output argument

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

       Alternatively, for compatibility with previous versions of MATPOWER,
       some of the results can be returned as individual output arguments:

       [baseMVA, bus, gen, branch, success, et] = runcpf(...);

   Generator reactive power limits are ignored for the continuation power flow.

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

   See also MPOPTION, RUNPF.

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

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005