toggle_softlims

TOGGLE_SOFTLIMS Relax DC optimal power flow branch limits.
   MPC = TOGGLE_SOFTLIMS(MPC, 'on')
   MPC = TOGGLE_SOFTLIMS(MPC, 'off')
   T_F = TOGGLE_SOFTLIMS(MPC, 'status')

   Enables, disables or checks the status of a set of OPF userfcn
   callbacks to implement relaxed branch flow limits a DC optimal flow model.

   These callbacks expect to find a 'softlims' field in the input MPC,
   where MPC.softlims is a struct with the following fields:
       idx     (optional) n x 1, index vector for branches whos flow
               limits are to be relaxed, default is to use all on-line
               branches with non-zero limits specified in RATE_A
       cost    (optional) n x 1, linear marginal cost per MW of exceeding
               RATE_A. Can optionally be a scalar, in which case it is
               applied to all soft limits. Default if not specified is
               $1000/MW.

   The 'int2ext' callback also packages up results and stores them in
   the following output fields of results.softlims:
       overload - nl x 1, amount of overload of each line in MW
       ovl_cost - nl x 1, total cost of overload in $/hr

   The shadow prices on the soft limit constraints are also returned in the
   MU_SF and MU_ST columns of the branch matrix.
   Note: These shadow prices are equal to the corresponding hard limit
       shadow prices when the soft limits are not violated. When violated,
       the shadow price on a soft limit constraint is equal to the
       user-specified soft limit violation cost.

   See also ADD_USERFCN, REMOVE_USERFCN, RUN_USERFCN, T_CASE30_USERFCNS.

0001 function mpc = toggle_softlims(mpc, on_off)
0002 %TOGGLE_SOFTLIMS Relax DC optimal power flow branch limits.
0003 %   MPC = TOGGLE_SOFTLIMS(MPC, 'on')
0004 %   MPC = TOGGLE_SOFTLIMS(MPC, 'off')
0005 %   T_F = TOGGLE_SOFTLIMS(MPC, 'status')
0006 %
0007 %   Enables, disables or checks the status of a set of OPF userfcn
0008 %   callbacks to implement relaxed branch flow limits a DC optimal flow model.
0009 %
0010 %   These callbacks expect to find a 'softlims' field in the input MPC,
0011 %   where MPC.softlims is a struct with the following fields:
0012 %       idx     (optional) n x 1, index vector for branches whos flow
0013 %               limits are to be relaxed, default is to use all on-line
0014 %               branches with non-zero limits specified in RATE_A
0015 %       cost    (optional) n x 1, linear marginal cost per MW of exceeding
0016 %               RATE_A. Can optionally be a scalar, in which case it is
0017 %               applied to all soft limits. Default if not specified is
0018 %               $1000/MW.
0019 %
0020 %   The 'int2ext' callback also packages up results and stores them in
0021 %   the following output fields of results.softlims:
0022 %       overload - nl x 1, amount of overload of each line in MW
0023 %       ovl_cost - nl x 1, total cost of overload in $/hr
0024 %
0025 %   The shadow prices on the soft limit constraints are also returned in the
0026 %   MU_SF and MU_ST columns of the branch matrix.
0027 %   Note: These shadow prices are equal to the corresponding hard limit
0028 %       shadow prices when the soft limits are not violated. When violated,
0029 %       the shadow price on a soft limit constraint is equal to the
0030 %       user-specified soft limit violation cost.
0031 %
0032 %   See also ADD_USERFCN, REMOVE_USERFCN, RUN_USERFCN, T_CASE30_USERFCNS.
0033 
0034 %   To do for future versions:
0035 %       - make softlims input field optional, convert all lines if missing
0036 %       Inputs:
0037 %       cost    n x 3, linear marginal cost per MW of exceeding each of
0038 %               RATE_A, RATE_B and RATE_C. Columns 2 and 3 are optional.
0039 %       brkpts  n x npts, allow to specify arbitrary breakpoints at which
0040 %               cost increases, defined as percentages above RATE_A.
0041 %       base_flow   n x 1, arbitrary baseline (other than RATE_A)
0042 
0043 %   MATPOWER
0044 %   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
0045 %   by Ray Zimmerman, PSERC Cornell
0046 %
0047 %   This file is part of MATPOWER.
0048 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0049 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0050 
0051 if strcmp(upper(on_off), 'ON')
0052     %% check for proper softlims inputs
0053     %% (inputs are optional, defaults handled in ext2int callback)
0054     
0055     %% add callback functions
0056     %% note: assumes all necessary data included in 1st arg (mpc, om, results)
0057     %%       so, no additional explicit args are needed
0058     mpc = add_userfcn(mpc, 'ext2int', @userfcn_softlims_ext2int);
0059     mpc = add_userfcn(mpc, 'formulation', @userfcn_softlims_formulation);
0060     mpc = add_userfcn(mpc, 'int2ext', @userfcn_softlims_int2ext);
0061     mpc = add_userfcn(mpc, 'printpf', @userfcn_softlims_printpf);
0062     mpc = add_userfcn(mpc, 'savecase', @userfcn_softlims_savecase);
0063     mpc.userfcn.status.softlims = 1;
0064 elseif strcmp(upper(on_off), 'OFF')
0065     mpc = remove_userfcn(mpc, 'savecase', @userfcn_softlims_savecase);
0066     mpc = remove_userfcn(mpc, 'printpf', @userfcn_softlims_printpf);
0067     mpc = remove_userfcn(mpc, 'int2ext', @userfcn_softlims_int2ext);
0068     mpc = remove_userfcn(mpc, 'formulation', @userfcn_softlims_formulation);
0069     mpc = remove_userfcn(mpc, 'ext2int', @userfcn_softlims_ext2int);
0070     mpc.userfcn.status.softlims = 0;
0071 elseif strcmp(upper(on_off), 'STATUS')
0072     if isfield(mpc, 'userfcn') && isfield(mpc.userfcn, 'status') && ...
0073             isfield(mpc.userfcn.status, 'softlims')
0074         mpc = mpc.userfcn.status.softlims;
0075     else
0076         mpc = 0;
0077     end
0078 else
0079     error('toggle_softlims: 2nd argument must be ''on'', ''off'' or ''status''');
0080 end
0081 
0082 
0083 %%-----  ext2int  ------------------------------------------------------
0084 function mpc = userfcn_softlims_ext2int(mpc, args)
0085 %
0086 %   mpc = userfcn_softlims_ext2int(mpc, args)
0087 %
0088 %   This is the 'ext2int' stage userfcn callback that prepares the input
0089 %   data for the formulation stage. It expects to find a 'softlims' field in
0090 %   mpc as described above. The optional args are not currently used.
0091 
0092 %% define named indices into data matrices
0093 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0094     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0095     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0096 
0097 %% check for proper softlims inputs
0098 default_cost = 1000;    %% used if cost is not specified
0099 if isfield(mpc, 'softlims')
0100     if ~isfield(mpc.softlims, 'idx')
0101         mpc.softlims.idx = [];
0102     end
0103     if ~isfield(mpc.softlims, 'cost')
0104         mpc.softlims.cost = default_cost;
0105     end
0106 else
0107     mpc.softlims.idx = [];
0108     mpc.softlims.cost = default_cost;
0109 end
0110 
0111 %% initialize some things
0112 s = mpc.softlims;
0113 o = mpc.order;
0114 nl0 = size(o.ext.branch, 1);    %% original number of branches
0115 nl = size(mpc.branch, 1);       %% number of on-line branches
0116 
0117 %% save softlims struct for external indexing
0118 mpc.order.ext.softlims = s;
0119 
0120 %%-----  convert stuff to internal indexing  -----
0121 s = softlims_defaults(s, o.ext.branch);     %% get defaults
0122 e2i = zeros(nl0, 1);
0123 e2i(o.branch.status.on) = (1:nl)';  %% ext->int branch index mapping
0124 s.idx = e2i(s.idx);
0125 k = find(s.idx == 0);   %% find idxes corresponding to off-line branches
0126 s.idx(k) = [];          %% delete them
0127 s.cost(k, :) = [];
0128 k = find(mpc.branch(s.idx, RATE_A) <= 0);   %% find branches w/o limits
0129 s.idx(k) = [];          %% delete them
0130 s.cost(k, :) = [];
0131 
0132 %%-----  remove hard limits on branches with soft limits  -----
0133 s.Pfmax = mpc.branch(s.idx, RATE_A) / mpc.baseMVA;  %% save limit first
0134 mpc.branch(s.idx, RATE_A) = 0;                      %% then remove it
0135 
0136 mpc.softlims = s;
0137 mpc.order.int.softlims = s;
0138 
0139 
0140 %%-----  formulation  --------------------------------------------------
0141 function om = userfcn_softlims_formulation(om, args)
0142 %
0143 %   om = userfcn_softlims_formulation(om, args)
0144 %
0145 %   This is the 'formulation' stage userfcn callback that defines the
0146 %   user costs and constraints for interface flow limits. It expects to
0147 %   find a 'softlims' field in the mpc stored in om, as described above. The
0148 %   optional args are not currently used.
0149 
0150 %% define named indices into data matrices
0151 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0152     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0153     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0154 
0155 %% initialize some things
0156 mpc = get_mpc(om);
0157 [baseMVA, bus, branch] = deal(mpc.baseMVA, mpc.bus, mpc.branch);
0158 s = mpc.softlims;
0159 
0160 %% form B matrices for DC model
0161 [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0162 n = size(Bf, 2);        %% dim of theta
0163 
0164 %% form constraints (flv is flow limit violation variable)
0165 %%   -Bf * Va - flv <=  Pfinj + Pfmax
0166 %%    Bf * Va - flv <= -Pfinj + Pfmax
0167 ns = length(s.idx);     %% number of soft limits
0168 I = speye(ns, ns);
0169 As = [-Bf(s.idx, :) -I; Bf(s.idx, :) -I];
0170 ls = -Inf(2*ns, 1);
0171 us = [   Pfinj(s.idx) + s.Pfmax;
0172         -Pfinj(s.idx) + s.Pfmax ];
0173 
0174 %% costs on flv variable
0175 Cw = s.cost(:, 1) * mpc.baseMVA;
0176 
0177 %% add vars, costs, constraints
0178 om = add_vars(om, 'flv', ns, zeros(ns, 1), zeros(ns, 1), Inf(ns, 1));
0179 om = add_costs(om, 'vc', struct('N', I, 'Cw', Cw), {'flv'});
0180 om = add_constraints(om, 'softlims',  As, ls, us, {'Va', 'flv'});   %% 2*ns
0181 
0182 
0183 %%-----  int2ext  ------------------------------------------------------
0184 function results = userfcn_softlims_int2ext(results, args)
0185 %
0186 %   results = userfcn_softlims_int2ext(results, args)
0187 %
0188 %   This is the 'int2ext' stage userfcn callback that converts everything
0189 %   back to external indexing and packages up the results. It expects to
0190 %   find a 'softlims' field in the results struct as described for mpc above.
0191 %   It also expects the results to contain solved branch flows and linear
0192 %   constraints named 'softlims' which are used to populate output fields
0193 %   in results.softlims. The optional args are not currently used.
0194 
0195 %% define named indices into data matrices
0196 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0197     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0198     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0199 
0200 %% get internal softlims struct
0201 s = results.softlims;
0202 
0203 %%-----  convert stuff back to external indexing  -----
0204 o = results.order;
0205 nl0 = size(o.ext.branch, 1);    %% original number of branches
0206 nl = size(results.branch, 1);   %% number of on-line branches
0207 o.branch.status.on;
0208 results.softlims = results.order.ext.softlims;
0209 
0210 %%-----  restore hard limits  -----
0211 results.branch(s.idx, RATE_A) = s.Pfmax * results.baseMVA;
0212 
0213 %%-----  results post-processing  -----
0214 %% get shadow prices
0215 n = size(results.lin.mu.u.softlims, 1) / 2;
0216 results.branch(s.idx, MU_ST) = results.lin.mu.u.softlims(1:n) / results.baseMVA;
0217 results.branch(s.idx, MU_SF) = results.lin.mu.u.softlims(n+1:end) / results.baseMVA;
0218 %% get overloads and overload costs
0219 results.softlims.overload = zeros(nl0, 1);
0220 k = find(results.branch(:, RATE_A) & ...
0221          abs(results.branch(:, PF)) > results.branch(:, RATE_A) );
0222 results.softlims.overload(o.branch.status.on(k)) = ...
0223         abs(results.branch(k, PF)) - results.branch(k, RATE_A);
0224 results.softlims.ovl_cost = zeros(size(results.softlims.overload));
0225 results.softlims.ovl_cost(o.branch.status.on(s.idx)) = ...
0226     results.softlims.overload(o.branch.status.on(s.idx)) .* s.cost(:, 1);
0227 
0228 
0229 %%-----  printpf  ------------------------------------------------------
0230 function results = userfcn_softlims_printpf(results, fd, mpopt, args)
0231 %
0232 %   results = userfcn_softlims_printpf(results, fd, mpopt, args)
0233 %
0234 %   This is the 'printpf' stage userfcn callback that pretty-prints the
0235 %   results. It expects a results struct, a file descriptor and a MATPOWER
0236 %   options struct. The optional args are not currently used.
0237 
0238 %% define named indices into data matrices
0239 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0240     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0241     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0242 
0243 %%-----  print results  -----
0244 ptol = 1e-6;        %% tolerance for displaying shadow prices
0245 isOPF           = isfield(results, 'f') && ~isempty(results.f);
0246 SUPPRESS        = mpopt.out.suppress_detail;
0247 if SUPPRESS == -1
0248     if size(results.bus, 1) > 500
0249         SUPPRESS = 1;
0250     else
0251         SUPPRESS = 0;
0252     end
0253 end
0254 OUT_ALL         = mpopt.out.all;
0255 OUT_FORCE       = mpopt.out.force;
0256 OUT_BRANCH      = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.branch);
0257 
0258 if isOPF && OUT_BRANCH && (results.success || OUT_FORCE)
0259     s = softlims_defaults(results.softlims, results.branch);   %% get defaults
0260     k = find(s.overload(s.idx) | sum(results.branch(s.idx, MU_SF:MU_ST), 2) > ptol);
0261 
0262     fprintf(fd, '\n================================================================================');
0263     fprintf(fd, '\n|     Soft Flow Limits                                                         |');
0264     fprintf(fd, '\n================================================================================');
0265     fprintf(fd, '\nBrnch   From   To     Flow      Limit    Overload     mu');
0266     fprintf(fd, '\n  #     Bus    Bus    (MW)      (MW)       (MW)     ($/MW)');
0267     fprintf(fd, '\n-----  -----  -----  --------  --------  --------  ---------');
0268     fprintf(fd, '\n%4d%7d%7d%10.2f%10.2f%10.2f%11.3f', ...
0269             [   s.idx(k), results.branch(s.idx(k), [F_BUS, T_BUS]), ...
0270                 results.branch(s.idx(k), [PF, RATE_A]), ...
0271                 s.overload(s.idx(k)), ...
0272                 sum(results.branch(s.idx(k), [MU_SF:MU_ST]), 2) ...
0273             ]');
0274     fprintf(fd, '\n                                         --------');
0275     fprintf(fd, '\n                                Total:%10.2f', ...
0276             sum(s.overload(s.idx(k))));
0277     fprintf(fd, '\n');
0278 end
0279 
0280 
0281 %%-----  savecase  -----------------------------------------------------
0282 function mpc = userfcn_softlims_savecase(mpc, fd, prefix, args)
0283 %
0284 %   mpc = userfcn_softlims_savecase(mpc, fd, mpopt, args)
0285 %
0286 %   This is the 'savecase' stage userfcn callback that prints the M-file
0287 %   code to save the 'softlims' field in the case file. It expects a
0288 %   MATPOWER case struct (mpc), a file descriptor and variable prefix
0289 %   (usually 'mpc.'). The optional args are not currently used.
0290 
0291 if isfield(mpc, 'softlims')
0292     s = mpc.softlims;
0293 
0294     fprintf(fd, '\n%%%%-----  Soft Flow Limit Data  -----%%%%\n');
0295 
0296     if isfield(s, 'idx')
0297         fprintf(fd, '%%%% branch indexes\n');
0298         fprintf(fd, '%%\tbranchidx\n');
0299         if isempty(s.idx)
0300             fprintf(fd, '%ssoftlims.idx = [];\n\n', prefix);
0301         else
0302             fprintf(fd, '%ssoftlims.idx = [\n', prefix);
0303             fprintf(fd, '\t%d;\n', s.idx);
0304             fprintf(fd, '];\n\n');
0305         end
0306     end
0307 
0308     fprintf(fd, '%%%% violation cost coefficients\n');
0309     fprintf(fd, '%%\trate_a_cost\n');
0310     fprintf(fd, '%ssoftlims.cost = [\n', prefix);
0311     fprintf(fd, '\t%g;\n', s.cost);
0312     fprintf(fd, '];\n');
0313 
0314     %% save output fields for solved case
0315     if isfield(mpc.softlims, 'overload')
0316         fprintf(fd, '\n%%%% overloads\n');
0317         fprintf(fd, '%%\toverload\n');
0318         fprintf(fd, '%ssoftlims.overload = [\n', prefix);
0319         fprintf(fd, '\t%g;\n', s.overload);
0320         fprintf(fd, '];\n');
0321 
0322         fprintf(fd, '\n%%%% overload costs\n');
0323         fprintf(fd, '%%\toverload_costs\n');
0324         fprintf(fd, '%ssoftlims.ovl_cost = [\n', prefix);
0325         fprintf(fd, '\t%g;\n', s.ovl_cost);
0326         fprintf(fd, '];\n');
0327     end
0328 end
0329 
0330 %%-----  softlims_defaults  --------------------------------------------
0331 function s = softlims_defaults(s, branch)
0332 %
0333 %   s = softlims_defaults(s, branch)
0334 %
0335 %   Takes a softlims struct that could have an empty 'idx' field or a
0336 %   scalar 'cost' field and fills them out with the defaults, where the
0337 %   default for 'idx' includes all on-line branches with non-zero RATE_A,
0338 %   and the default for the cost is to apply the scalar to each soft limit
0339 %   violation.
0340 
0341 %% define named indices into data matrices
0342 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0343     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0344     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0345 
0346 if isempty(s.idx)
0347     s.idx = find(branch(:, BR_STATUS) > 0 & branch(:, RATE_A) > 0);
0348 end
0349 if length(s.cost) == 1 && length(s.idx) > 1
0350     s.cost = s.cost * ones(size(s.idx));
0351 end