Home > matpower7.1 > lib > toggle_softlims.m

toggle_softlims

PURPOSE ^

TOGGLE_SOFTLIMS Relax DC optimal power flow branch limits.

SYNOPSIS ^

function mpc = toggle_softlims(mpc, on_off)

DESCRIPTION ^

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 inequality constraints for an OPF model.

   These callbacks expect to find a 'softlims' field in the input MPC,
   where MPC.softlims is a struct with fields corresponding to the
   possible limits, namely:
       VMIN, VMAX, RATE_A, PMIN, PMAX, QMIN, QMAX, ANGMAX, ANGMIN
   Each of these is itself a struct with the following fields, all of which
   are optional:
       idx     index of affected buses, branches, or generators. These are
               row indices into the respective matrices. The default is to
               include all online elements for which the constraint in
               question is not unbounded, except for generators, which also
               exclude those used to model dispatchable loads
               (i.e. those for which isload(gen) is true).
       busnum  for bus constraints, such as VMIN and VMAX, the affected
               buses can be specified by a vector of external bus numbers
               in the 'busnum' field instead of bus row indices in the 'idx'
               field. If both are present, 'idx' overrides 'busnum'.
       cost    linear marginal cost of exceeding the original limit
               The defaults are set as:
                   base_cost x 100 $/pu    for VMAX and VMIN
                   base_cost $/MW          for RATE_A, PMAX, and PMIN
                   base_cost $/MVAr        for QMAX, QMIN
                   base_cost $/deg         for ANGMAX, ANGMIN
               where base_cost is the maximum of $1000 and twice the
               maximum generator cost of all online generators.
       hl_mod  type of modification to hard limit, hl:
                   'none'    : do *not* add soft limit, no change to original
                               hard limit
                   'remove'  : add soft limit, relax hard limit by removing
                               it completely
                   'replace' : add soft limit, relax hard limit by replacing
                               original with value specified in hl_val
                   'scale'   : add soft limit, relax hard limit by scaling
                               original by value specified in hl_val
                   'shift'   : add soft limit, relax hard limit by shifting
                               original by value specified in hl_val
       hl_val  value used to modify hard limit according to hl_mod. Ignored
               for 'none' and 'remove', required for 'replace', and optional,
               with the following defaults, for 'scale' and 'shift':
                   'scale' :   2 for positive upper limits or negative lower
                               limits, 0.5 otherwise
                   'shift' :   0.25 for VMAX and VMIN, 10 otherwise

   For limits that are left unspecified in the structure, the default
   behavior is determined by the value of the mpopt.opf.softlims.default
   option. If mpopt.opf.softlims.default = 0, then the unspecified softlims
   are ignored (hl_mod = 'none', i.e. original hard limits left in place).
   If mpopt.opf.softlims.default = 1 (default), then the unspecified
   softlims are enabled with default values, which specify to 'remove'
   the hard limit, except in the case of VMIN and PMIN, whose defaults
   are set as follows:
         .VMIN
           .hl_mod = 'replace'
           .hl_val = 0
         .PMIN
           .hl_mod = 'replace'
           .hl_val = 0     for normal generators (PMIN > 0)
           .hl_val = -Inf  for for generators with PMIN < 0 AND PMAX > 0

   With mpopt.opf.softlims.default = 0, it is still possible to enable a
   softlim with default values by setting that specification to an empty
   struct. E.g. mpc.softlims.VMAX = struct() would enable a default
   softlim on VMAX.

   The 'int2ext' callback also packages up results and stores them in
   the following output fields of results.softlims.(lim), where lim is
   one of the above mentioned limits:
       overload - amount of overload, i.e. violation of hard-limit.
       ovl_cost - total cost of overload in $/hr

   The shadow prices on the soft limit constraints are also returned in the
   relevant columns of the respective matrices (MU_SF, MU_ST for RATE_A,
   MU_VMAX for VMAX, etc.)
   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 + the shadow price on any
       binding remaining hard limit.

   See also ADD_USERFCN, REMOVE_USERFCN, RUN_USERFCN, T_OPF_SOFTLIMS.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 inequality constraints for an OPF model.
0009 %
0010 %   These callbacks expect to find a 'softlims' field in the input MPC,
0011 %   where MPC.softlims is a struct with fields corresponding to the
0012 %   possible limits, namely:
0013 %       VMIN, VMAX, RATE_A, PMIN, PMAX, QMIN, QMAX, ANGMAX, ANGMIN
0014 %   Each of these is itself a struct with the following fields, all of which
0015 %   are optional:
0016 %       idx     index of affected buses, branches, or generators. These are
0017 %               row indices into the respective matrices. The default is to
0018 %               include all online elements for which the constraint in
0019 %               question is not unbounded, except for generators, which also
0020 %               exclude those used to model dispatchable loads
0021 %               (i.e. those for which isload(gen) is true).
0022 %       busnum  for bus constraints, such as VMIN and VMAX, the affected
0023 %               buses can be specified by a vector of external bus numbers
0024 %               in the 'busnum' field instead of bus row indices in the 'idx'
0025 %               field. If both are present, 'idx' overrides 'busnum'.
0026 %       cost    linear marginal cost of exceeding the original limit
0027 %               The defaults are set as:
0028 %                   base_cost x 100 $/pu    for VMAX and VMIN
0029 %                   base_cost $/MW          for RATE_A, PMAX, and PMIN
0030 %                   base_cost $/MVAr        for QMAX, QMIN
0031 %                   base_cost $/deg         for ANGMAX, ANGMIN
0032 %               where base_cost is the maximum of $1000 and twice the
0033 %               maximum generator cost of all online generators.
0034 %       hl_mod  type of modification to hard limit, hl:
0035 %                   'none'    : do *not* add soft limit, no change to original
0036 %                               hard limit
0037 %                   'remove'  : add soft limit, relax hard limit by removing
0038 %                               it completely
0039 %                   'replace' : add soft limit, relax hard limit by replacing
0040 %                               original with value specified in hl_val
0041 %                   'scale'   : add soft limit, relax hard limit by scaling
0042 %                               original by value specified in hl_val
0043 %                   'shift'   : add soft limit, relax hard limit by shifting
0044 %                               original by value specified in hl_val
0045 %       hl_val  value used to modify hard limit according to hl_mod. Ignored
0046 %               for 'none' and 'remove', required for 'replace', and optional,
0047 %               with the following defaults, for 'scale' and 'shift':
0048 %                   'scale' :   2 for positive upper limits or negative lower
0049 %                               limits, 0.5 otherwise
0050 %                   'shift' :   0.25 for VMAX and VMIN, 10 otherwise
0051 %
0052 %   For limits that are left unspecified in the structure, the default
0053 %   behavior is determined by the value of the mpopt.opf.softlims.default
0054 %   option. If mpopt.opf.softlims.default = 0, then the unspecified softlims
0055 %   are ignored (hl_mod = 'none', i.e. original hard limits left in place).
0056 %   If mpopt.opf.softlims.default = 1 (default), then the unspecified
0057 %   softlims are enabled with default values, which specify to 'remove'
0058 %   the hard limit, except in the case of VMIN and PMIN, whose defaults
0059 %   are set as follows:
0060 %         .VMIN
0061 %           .hl_mod = 'replace'
0062 %           .hl_val = 0
0063 %         .PMIN
0064 %           .hl_mod = 'replace'
0065 %           .hl_val = 0     for normal generators (PMIN > 0)
0066 %           .hl_val = -Inf  for for generators with PMIN < 0 AND PMAX > 0
0067 %
0068 %   With mpopt.opf.softlims.default = 0, it is still possible to enable a
0069 %   softlim with default values by setting that specification to an empty
0070 %   struct. E.g. mpc.softlims.VMAX = struct() would enable a default
0071 %   softlim on VMAX.
0072 %
0073 %   The 'int2ext' callback also packages up results and stores them in
0074 %   the following output fields of results.softlims.(lim), where lim is
0075 %   one of the above mentioned limits:
0076 %       overload - amount of overload, i.e. violation of hard-limit.
0077 %       ovl_cost - total cost of overload in $/hr
0078 %
0079 %   The shadow prices on the soft limit constraints are also returned in the
0080 %   relevant columns of the respective matrices (MU_SF, MU_ST for RATE_A,
0081 %   MU_VMAX for VMAX, etc.)
0082 %   Note: These shadow prices are equal to the corresponding hard limit
0083 %       shadow prices when the soft limits are not violated. When violated,
0084 %       the shadow price on a soft limit constraint is equal to the
0085 %       user-specified soft limit violation cost + the shadow price on any
0086 %       binding remaining hard limit.
0087 %
0088 %   See also ADD_USERFCN, REMOVE_USERFCN, RUN_USERFCN, T_OPF_SOFTLIMS.
0089 
0090 %   To do for future versions:
0091 %       Inputs:
0092 %       cost    n x 3, linear marginal cost per MW of exceeding each of
0093 %               RATE_A, RATE_B and RATE_C. Columns 2 and 3 are optional.
0094 %       brkpts  n x npts, allow to specify arbitrary breakpoints at which
0095 %               cost increases, defined as percentages above RATE_A.
0096 %       base_flow   n x 1, arbitrary baseline (other than RATE_A)
0097 
0098 %   MATPOWER
0099 %   Copyright (c) 2009-2018, Power Systems Engineering Research Center (PSERC)
0100 %   by Ray Zimmerman, PSERC Cornell
0101 %   and Eran Schweitzer, Arizona State University
0102 %
0103 %   This file is part of MATPOWER.
0104 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0105 %   See https://matpower.org for more info.
0106 
0107 if strcmp(upper(on_off), 'ON')
0108     %% check for proper softlims inputs
0109     %% (inputs are optional, defaults handled in ext2int callback)
0110     
0111     %% add callback functions
0112     %% note: assumes all necessary data included in 1st arg (mpc, om, results)
0113     %%       so, no additional explicit args are needed
0114     mpc = add_userfcn(mpc, 'ext2int', @userfcn_softlims_ext2int);
0115     mpc = add_userfcn(mpc, 'formulation', @userfcn_softlims_formulation);
0116     mpc = add_userfcn(mpc, 'int2ext', @userfcn_softlims_int2ext);
0117     mpc = add_userfcn(mpc, 'printpf', @userfcn_softlims_printpf);
0118     mpc = add_userfcn(mpc, 'savecase', @userfcn_softlims_savecase);
0119     mpc.userfcn.status.softlims = 1;
0120 elseif strcmp(upper(on_off), 'OFF')
0121     mpc = remove_userfcn(mpc, 'savecase', @userfcn_softlims_savecase);
0122     mpc = remove_userfcn(mpc, 'printpf', @userfcn_softlims_printpf);
0123     mpc = remove_userfcn(mpc, 'int2ext', @userfcn_softlims_int2ext);
0124     mpc = remove_userfcn(mpc, 'formulation', @userfcn_softlims_formulation);
0125     mpc = remove_userfcn(mpc, 'ext2int', @userfcn_softlims_ext2int);
0126     mpc.userfcn.status.softlims = 0;
0127 elseif strcmp(upper(on_off), 'STATUS')
0128     if isfield(mpc, 'userfcn') && isfield(mpc.userfcn, 'status') && ...
0129             isfield(mpc.userfcn.status, 'softlims')
0130         mpc = mpc.userfcn.status.softlims;
0131     else
0132         mpc = 0;
0133     end
0134 else
0135     error('toggle_softlims: 2nd argument must be ''on'', ''off'' or ''status''');
0136 end
0137 
0138 
0139 %%-----  ext2int  ------------------------------------------------------
0140 function mpc = userfcn_softlims_ext2int(mpc, mpopt, args)
0141 %
0142 %   mpc = userfcn_softlims_ext2int(mpc, mpopt, args)
0143 %
0144 %   This is the 'ext2int' stage userfcn callback that prepares the input
0145 %   data for the formulation stage. It expects to find a 'softlims' field in
0146 %   mpc as described above. The optional args are not currently used.
0147 
0148 %% define named indices into data matrices
0149 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0150     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0151 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0152     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0153     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0154 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0155     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0156     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0157 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0158 
0159 %% check for cartesian voltage coordinates
0160 if ~isempty(mpopt) && isfield(mpopt, 'opf') && mpopt.opf.v_cartesian
0161     error('userfcn_softlims_ext2int: TOGGLE_SOFTLIMS is not implemented for OPF with cartesian voltages. Please set the MPOPT.opf.v_cartesian option to 0 for use with TOGGLE_SOFTLIMS.');
0162 end
0163 
0164 %% structures used to index into softlims in a loop
0165 lims = softlims_lim2mat(mpopt);
0166 mat2lims = softlims_mat2lims(mpopt);
0167 
0168 %% set up softlims defaults
0169 mpc.softlims = softlims_defaults(mpc, mpopt);
0170 
0171 %% initialize some things
0172 slims = softlims_init(mpc, mpopt);
0173 o = mpc.order;
0174 
0175 %% save softlims struct with external indexing
0176 mpc.order.ext.softlims = slims;
0177 
0178 %%-----  convert stuff to internal indexing  -----
0179 for m = fieldnames(mat2lims).'
0180     mat_name = m{:};
0181     n0  = size(o.ext.(mat_name), 1);    %% original number
0182     n   = size(mpc.(mat_name), 1);      %% on-line number
0183     e2i = zeros(n0, 1);                 %% ext->int index mapping
0184     if strcmp(mat_name, 'gen')
0185         %% for gens, account for permutation of gen matrix
0186         e2i(o.gen.status.on(o.gen.i2e)) = (1:n)';
0187     else
0188         e2i(o.(mat_name).status.on) = (1:n)';
0189     end
0190     for lm = mat2lims.(mat_name)
0191         lim = lm{:};
0192         s = slims.(lim);
0193         if ~strcmp( s.hl_mod, 'none')
0194             s.idx = e2i(s.idx);
0195             k = find(s.idx == 0);    %% find idxs corresponding to off-line elements
0196             s.idx(k)     = [];       %% delete them
0197             s.cost(k, :) = [];
0198             s.sav(k)     = [];
0199             s.ub(k)      = [];
0200             if isfield(s, 'hl_val') && ~isscalar(s.hl_val)
0201                 s.hl_val(k)  = [];
0202             end
0203             slims.(lim) = s;
0204         end
0205     end
0206 end
0207 
0208 %%-----  remove hard limits on elements with soft limits  -----
0209 %% Note: any new hard limits will be implemented as upper bounds
0210 %%       on the soft limit violation variable
0211 for lm = fieldnames(lims).'
0212     lim = lm{:};
0213     s = slims.(lim);
0214     if ~strcmp(s.hl_mod, 'none')
0215         mat_name = lims.(lim).mat_name; %% mpc sub-matrix name
0216         col = lims.(lim).col;           %% mpc sub-matrix column
0217         mpc.(mat_name)(s.idx, col) = s.rval;
0218     end
0219 end
0220 
0221 mpc.softlims = slims;
0222 mpc.order.int.softlims = slims;
0223 
0224 
0225 %%-----  formulation  --------------------------------------------------
0226 function om = userfcn_softlims_formulation(om, mpopt, args)
0227 %
0228 %   om = userfcn_softlims_formulation(om, mpopt, args)
0229 %
0230 %   This is the 'formulation' stage userfcn callback that defines the
0231 %   user costs and constraints for interface flow limits. It expects to
0232 %   find a 'softlims' field in the mpc stored in om, as described above. The
0233 %   optional args are not currently used.
0234 
0235 %% define named indices into data matrices
0236 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0237     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0238     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0239 
0240 %% initialize some things
0241 mpc = om.get_mpc();
0242 [baseMVA, bus, branch] = deal(mpc.baseMVA, mpc.bus, mpc.branch);
0243 nb = size(bus, 1);
0244 ng = size(mpc.gen, 1);
0245 
0246 %% structures used to index into softlims in a loop
0247 lims = softlims_lim2mat(mpopt);
0248 
0249 %% temporarily save mpopt in om for use by int2ext callback
0250 om.userdata.mpopt = mpopt;
0251 
0252 %%-----  add variables, costs, and constraints  -----
0253 for lm = fieldnames(lims).'
0254     lim = lm{:};
0255     s = mpc.softlims.(lim);
0256     if strcmp(s.hl_mod, 'none')
0257         continue;
0258     end
0259     varname = ['s_', lower(lim)];
0260     cstname = ['cs_', lower(lim)];
0261     ns = length(s.idx);  %% number of softlims
0262 
0263     %% variables and costs
0264     switch lim
0265         case {'VMIN', 'VMAX'}
0266             ub = s.ub;                      %% bound already in pu
0267             Cw = s.cost(:, 1);              %% cost already in $/pu
0268         case {'RATE_A', 'PMIN', 'PMAX', 'QMIN', 'QMAX'}
0269             ub = s.ub / baseMVA;            %% bound MVA -> pu
0270             Cw = s.cost(:, 1) * baseMVA;    %% cost in $/MVA -> $/pu
0271         case {'ANGMIN', 'ANGMAX'}
0272             ub = s.ub * pi/180;             %% bound deg -> rad
0273             Cw = s.cost(:, 1) * 180/pi;     %% cost in $/deg -> $/rad
0274         otherwise
0275             error('userfcn_soflims_formulation: limit %s is unknown ', lim)
0276     end
0277     om.add_var(varname, ns, zeros(ns, 1), zeros(ns, 1), ub);
0278     om.add_quad_cost(cstname, [], Cw, 0, {varname});
0279 
0280     %% linear constraints that are identical for DC and AC formulations
0281     switch lim
0282         case 'ANGMIN'
0283             %% theta_f - theta_t + s_angmin >= s.sav
0284             ns = length(s.idx);
0285             Av = sparse([1:ns,1:ns].', ...
0286                     [branch(s.idx, F_BUS); branch(s.idx, T_BUS)], ...
0287                     [ones(ns,1);-ones(ns,1)], ns, nb);
0288             As = speye(ns);
0289             lb = s.sav * pi/180;
0290             ub = Inf(ns,1);
0291 
0292             om.add_lin_constraint('soft_angmin', [Av As], lb, ub, ...
0293                 {'Va', 's_angmin'});
0294         case 'ANGMAX'
0295             %% theta_f - theta_t - s_angmax <= s.sav
0296             ns = length(s.idx);
0297             Av = sparse([1:ns,1:ns].', ...
0298                     [branch(s.idx, F_BUS); branch(s.idx, T_BUS)], ...
0299                     [ones(ns,1);-ones(ns,1)], ns, nb);
0300             As = speye(ns);
0301             lb = -Inf(ns,1);
0302             ub = s.sav * pi/180;
0303 
0304             om.add_lin_constraint('soft_angmax', [Av -As], lb, ub, ...
0305                 {'Va', 's_angmax'});
0306         case 'PMIN'
0307             %% Pg + s_pmin >= s.sav
0308             ns = length(s.idx);
0309             Av = sparse(1:ns, s.idx, 1, ns, ng);
0310             As = speye(ns);
0311             lb = s.sav / baseMVA;
0312             ub = Inf(ns,1);
0313 
0314             om.add_lin_constraint('soft_pmin', [Av As], lb, ub, ...
0315                 {'Pg', 's_pmin'});
0316         case 'PMAX'
0317             %% Pg - s_pmax <= s.sav
0318             ns = length(s.idx);
0319             Av = sparse(1:ns, s.idx, 1, ns, ng);
0320             As = speye(ns);
0321             lb = -Inf(ns,1);
0322             ub = s.sav / baseMVA;
0323 
0324             om.add_lin_constraint('soft_pmax', [Av -As], lb, ub, ...
0325                 {'Pg', 's_pmax'});
0326     end
0327 end
0328 
0329 %% RATE_A (different for DC and AC) and other AC-only constraints
0330 isDC = strcmp(mpopt.model, 'DC');
0331 lims = {'RATE_A'};
0332 if ~isDC
0333     lims = {'VMIN', 'VMAX', 'QMIN', 'QMAX', lims{:}};
0334 end
0335 for lm = lims
0336     lim = lm{:};
0337     s = mpc.softlims.(lim);
0338     if strcmp(s.hl_mod, 'none') %% skip
0339         continue;
0340     end
0341     ns = length(s.idx);         %% number of softlims
0342 
0343     switch lim
0344         case 'RATE_A'
0345             if isDC
0346                 %% fetch Bf matrix for DC model
0347                 Bf = om.get_userdata('Bf');
0348                 Pfinj = om.get_userdata('Pfinj');
0349 
0350                 %% form constraints
0351                 %%    Bf * Va - s_rate_a <= -Pfinj + Pfmax
0352                 %%   -Bf * Va - s_rate_a <=  Pfinj + Pfmax
0353                 I = speye(ns, ns);
0354                 Asf = [ Bf(s.idx, :) -I];
0355                 Ast = [-Bf(s.idx, :) -I];
0356                 lsf = -Inf(ns, 1);
0357                 lst = lsf;
0358                 usf = -Pfinj(s.idx) + s.sav/baseMVA;
0359                 ust =  Pfinj(s.idx) + s.sav/baseMVA;
0360                 vs = {'Va', 's_rate_a'};
0361 
0362                 om.add_lin_constraint('softPf', Asf, lsf, usf, vs);    %% ns
0363                 om.add_lin_constraint('softPt', Ast, lst, ust, vs);    %% ns
0364             else
0365                 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0366 
0367                 fcn = @(x)softlims_fcn(x, mpc, Yf(s.idx, :), Yt(s.idx, :), ...
0368                                        s.idx, mpopt, s.sav/baseMVA);
0369                 hess = @(x, lam)softlims_hess(x, lam, mpc, Yf(s.idx, :), ...
0370                                               Yt(s.idx, :), s.idx, mpopt);
0371                 om.add_nln_constraint({'softSf', 'softSt'}, [ns;ns], 0, ...
0372                                         fcn, hess, {'Va', 'Vm', 's_rate_a'});
0373             end
0374         case 'VMIN'
0375             %% Vm + s_vmin >= s.sav
0376             Av  = sparse(1:ns, s.idx, 1, ns, nb);
0377             As  = speye(ns);
0378             lb  = s.sav;
0379             ub  = Inf(ns, 1);
0380 
0381             om.add_lin_constraint('soft_vmin', [Av As], lb, ub, ...
0382                                     {'Vm', 's_vmin'});
0383         case 'VMAX'
0384             %% Vm - s_vmax <= s.sav
0385             Av  = sparse(1:ns, s.idx, 1, ns, nb);
0386             As  = speye(ns);
0387             lb  = -Inf(ns, 1);
0388             ub  = s.sav;
0389 
0390             om.add_lin_constraint('soft_vmax', [Av -As], lb, ub, ...
0391                                     {'Vm', 's_vmax'});
0392         case 'QMIN'
0393             %% Qg + s_pmin >= s.max
0394             Av  = sparse(1:ns, s.idx, 1, ns, ng);
0395             As  = speye(ns);
0396             lb  = s.sav / baseMVA;
0397             ub  = Inf(ns, 1);
0398 
0399             om.add_lin_constraint('soft_qmin', [Av As], lb, ub, ...
0400                                     {'Qg', 's_qmin'});
0401         case 'QMAX'
0402             %% Qg - s_pmax <= s.sav
0403             Av  = sparse(1:ns, s.idx, 1, ns, ng);
0404             As  = speye(ns);
0405             lb  = -Inf(ns, 1);
0406             ub  = s.sav / baseMVA;
0407 
0408             om.add_lin_constraint('soft_qmax', [Av -As], lb, ub, ...
0409                                     {'Qg', 's_qmax'});
0410     end
0411 end
0412 
0413 
0414 %%-----  int2ext  ------------------------------------------------------
0415 function results = userfcn_softlims_int2ext(results, mpopt, args)
0416 %
0417 %   results = userfcn_softlims_int2ext(results, mpopt, args)
0418 %
0419 %   This is the 'int2ext' stage userfcn callback that converts everything
0420 %   back to external indexing and packages up the results. It expects to
0421 %   find a 'softlims' field in the results struct as described for mpc above.
0422 %   It also expects the results to contain values in bus, gen and branch
0423 %   and constraints for the specified softlims, with the following possible
0424 %   names, which are used to populate output fields in results.softlims:
0425 %       soft_Pf, softPt (DC only)
0426 %       soft_Sf, softSt (AC only)
0427 %       soft_angmin (both)
0428 %       soft_angmax (both)
0429 %       soft_pmin (both)
0430 %       soft_pmax (both)
0431 %       soft_vmin (AC only)
0432 %       soft_vmax (AC only)
0433 %       soft_qmin (AC only)
0434 %       soft_qmax (AC only)
0435 %   The optional args are not currently used.
0436 
0437 %% define named indices into data matrices
0438 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0439     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0440 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0441     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0442     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0443 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0444     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0445     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0446 
0447 %% get internal softlims struct and mpopt
0448 slims = results.softlims;
0449 isOPF = isfield(results, 'f') && ~isempty(results.f);
0450 if isOPF
0451     mpopt = results.om.get_userdata('mpopt');   %% extract and remove mpopt from om
0452     results.om.userdata = rmfield(results.om.userdata, 'mpopt');
0453     isDC = strcmp(mpopt.model, 'DC');
0454 end
0455 
0456 %% structures used to index into softlims in a loop
0457 lims = softlims_lim2mat(mpopt);
0458 
0459 %%-----  convert stuff back to external indexing  -----
0460 o = results.order;
0461 results.softlims = o.ext.softlims;
0462 
0463 %%-----  restore hard limits  -----
0464 for lm = fieldnames(lims).'
0465     lim = lm{:};
0466     s = slims.(lim);
0467     if strcmp(s.hl_mod, 'none')
0468         continue;
0469     end
0470     mat_name = lims.(lim).mat_name; %% mpc sub-matrix name
0471     col = lims.(lim).col;           %% mpc sub-matrix column
0472     results.(mat_name)(s.idx, col) = s.sav;
0473 end
0474 
0475 %%-----  remove rval, sav and ub fields  -----
0476 for lm = fieldnames(lims).'
0477     lim = lm{:};
0478     if strcmp(slims.(lim).hl_mod, 'none')
0479         continue;
0480     end
0481     results.softlims.(lim) = rmfield(results.softlims.(lim), 'sav');
0482     results.softlims.(lim) = rmfield(results.softlims.(lim), 'rval');
0483     results.softlims.(lim) = rmfield(results.softlims.(lim), 'ub');
0484 end
0485 
0486 %%-----  results post-processing  -----
0487 %% get overloads and overload costs
0488 tol = 1e-8;
0489 if isOPF
0490     for lm = fieldnames(lims).'
0491         lim = lm{:};
0492         s = slims.(lim);
0493         if strcmp(s.hl_mod, 'none')
0494             continue;
0495         end
0496 
0497         mat_name  = lims.(lim).mat_name;    %% mpc sub-matrix name
0498         n0 = size(o.ext.(mat_name), 1);     %% original number
0499         n  = size(results.(mat_name), 1);   %% number on-line
0500         results.softlims.(lim).overload = zeros(n0, 1);
0501         results.softlims.(lim).ovl_cost = zeros(n0, 1);
0502         varname = ['s_', lower(lim)];  %% variable name
0503         if isfield(results.var.val, varname)
0504             var = results.var.val.(varname);    %% violation in p.u.
0505             var(var < tol) = 0;                 %% zero out violations < tol
0506         end
0507         switch lim
0508             case {'VMIN', 'VMAX'}               %% p.u. (no change)
0509                 if isDC
0510                     continue;
0511                 end
0512             case {'RATE_A', 'PMIN', 'PMAX'}     %% p.u. -> MVA|MW
0513                 var = var * results.baseMVA;
0514             case {'QMIN', 'QMAX'}               %% p.u. -> MVAr
0515                 if isDC
0516                     continue;
0517                 end
0518                 var = var * results.baseMVA;
0519             case {'ANGMIN', 'ANGMAX'}           %% rad -> deg
0520                 var = var * 180/pi;
0521             otherwise
0522                 error('userfcn_soflims_formulation: limit %s is unknown ', lim)
0523         end
0524 
0525         %% NOTE: o.(mat_name).status.on is an nx1 vector where n is the INTERNAL
0526         %% number of elements. Entries are the EXTERNAL locations (row numbers).
0527         if strcmp(mat_name, 'gen')
0528             results.softlims.(lim).overload(o.gen.status.on(o.gen.i2e(s.idx))) = var;
0529             results.softlims.(lim).ovl_cost(o.gen.status.on(o.gen.i2e(s.idx))) = var .* s.cost(:, 1);
0530         else
0531             results.softlims.(lim).overload(o.(mat_name).status.on(s.idx)) = var;
0532             results.softlims.(lim).ovl_cost(o.(mat_name).status.on(s.idx)) = var .* s.cost(:, 1);
0533         end
0534 
0535         cname = ['soft_' lower(lim)];
0536         switch lim
0537             case {'ANGMAX', 'PMAX', 'VMAX', 'QMAX'}
0538                 mu = results.lin.mu.u.(cname);
0539             case {'ANGMIN', 'PMIN', 'VMIN', 'QMIN'}
0540                 mu = results.lin.mu.l.(cname);
0541             case 'RATE_A'
0542                 if isDC
0543                     muF = results.lin.mu.u.softPf;
0544                     muT = results.lin.mu.u.softPt;
0545                 else
0546                     muF = results.nli.mu.softSf;
0547                     muT = results.nli.mu.softSt;
0548                 end
0549         end
0550         switch lim
0551             case 'ANGMAX'
0552                 results.branch(slims.ANGMAX.idx, MU_ANGMAX) = mu * pi/180;
0553             case 'ANGMIN'
0554                 results.branch(slims.ANGMIN.idx, MU_ANGMIN) = mu * pi/180;
0555             case 'PMAX'
0556                 results.gen(slims.PMAX.idx, MU_PMAX) = mu / results.baseMVA;
0557             case 'PMIN'
0558                 results.gen(slims.PMIN.idx, MU_PMIN) = mu / results.baseMVA;
0559             case 'RATE_A'
0560                 results.branch(slims.RATE_A.idx, MU_SF) = muF / results.baseMVA;
0561                 results.branch(slims.RATE_A.idx, MU_ST) = muT / results.baseMVA;
0562         end
0563         if ~isDC    %% AC
0564             switch lim
0565                 case 'RATE_A'
0566                     if upper(mpopt.opf.flow_lim(1)) ~= 'P'
0567                         s = slims.RATE_A;
0568                         var = results.softlims.RATE_A.overload(o.branch.status.on(s.idx));
0569                         %% conversion factor for squared constraints (2*F)
0570                         cf = 2 * (s.sav + var) / results.baseMVA;
0571                         results.branch(s.idx, MU_SF) = ...
0572                             results.branch(s.idx, MU_SF) .* cf;
0573                         results.branch(s.idx, MU_ST) = ...
0574                             results.branch(s.idx, MU_ST) .* cf;
0575                     end
0576                 case 'VMAX'
0577                     results.bus(slims.VMAX.idx, MU_VMAX) = mu;
0578                 case 'VMIN'
0579                     results.bus(slims.VMIN.idx, MU_VMIN) = mu;
0580                 case 'QMAX'
0581                     results.gen(slims.QMAX.idx, MU_QMAX) = mu / results.baseMVA;
0582                 case 'QMIN'
0583                     results.gen(slims.QMIN.idx, MU_QMIN) = mu / results.baseMVA;
0584             end
0585         end
0586     end
0587 end
0588 
0589 
0590 %%-----  printpf  ------------------------------------------------------
0591 function results = userfcn_softlims_printpf(results, fd, mpopt, args)
0592 %
0593 %   results = userfcn_softlims_printpf(results, fd, mpopt, args)
0594 %
0595 %   This is the 'printpf' stage userfcn callback that pretty-prints the
0596 %   results. It expects a results struct, a file descriptor and a MATPOWER
0597 %   options struct. The optional args are not currently used.
0598 
0599 %% define named indices into data matrices
0600 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0601     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0602 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0603     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0604     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0605 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0606     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0607     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0608 
0609 %%-----  print results  -----
0610 ptol = 1e-6;        %% tolerance for displaying shadow prices
0611 isOPF           = isfield(results, 'f') && ~isempty(results.f);
0612 SUPPRESS        = mpopt.out.suppress_detail;
0613 if SUPPRESS == -1
0614     if size(results.bus, 1) > 500
0615         SUPPRESS = 1;
0616     else
0617         SUPPRESS = 0;
0618     end
0619 end
0620 OUT_ALL         = mpopt.out.all;
0621 OUT_FORCE       = mpopt.out.force;
0622 OUT_V_LIM       = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.v);
0623 OUT_LINE_LIM    = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.line);
0624 OUT_PG_LIM      = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.pg);
0625 OUT_QG_LIM      = OUT_ALL == 1 || (OUT_ALL == -1 && ~SUPPRESS && mpopt.out.lim.qg);
0626 
0627 if isOPF && (results.success || OUT_FORCE)
0628     [bus, gen, branch] = deal(results.bus, results.gen, results.branch);
0629     isAC = strcmp(mpopt.model, 'AC');
0630     if isAC
0631         s = results.softlims.VMAX;
0632         if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0633             fprintf(fd, '\n================================================================================');
0634             fprintf(fd, '\n|     Soft Voltage Upper Bounds                                                |');
0635             fprintf(fd, '\n================================================================================');
0636             k = find(s.overload(s.idx) | bus(s.idx, MU_VMAX) > ptol);
0637             if isempty(k)
0638                 fprintf(fd,'\nNo violations.\n');
0639             else
0640                 fprintf(fd, '\nBus    Voltage   Limit   Overload    mu');
0641                 fprintf(fd, '\n  #    Mag(pu)   (pu)     (pu)     ($/pu)');
0642                 fprintf(fd, '\n-----  -------  -------  -------  ---------');
0643                 fprintf(fd, '\n%5d%8.3f%9.3f%9.3f%11.3f',...
0644                     [ bus(s.idx(k), BUS_I), bus(s.idx(k),[VM, VMAX]),...
0645                       s.overload(s.idx(k)), ...
0646                       bus(s.idx(k), MU_VMAX)...
0647                     ]');
0648                 fprintf(fd, '\n                        --------');
0649                 fprintf(fd, '\n               Total:%10.2f', ...
0650                         sum(s.overload(s.idx(k))));
0651                 fprintf(fd, '\n');
0652             end
0653         end
0654         s = results.softlims.VMIN;
0655         if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0656             fprintf(fd, '\n================================================================================');
0657             fprintf(fd, '\n|     Soft Voltage Lower Bounds                                                |');
0658             fprintf(fd, '\n================================================================================');
0659             k = find(s.overload(s.idx) | bus(s.idx, MU_VMIN) > ptol);
0660             if isempty(k)
0661                 fprintf(fd,'\nNo violations.\n');
0662             else
0663                 fprintf(fd, '\n----------------------------------------');
0664                 fprintf(fd, '\n Bus   Voltage   Limit   Overload    mu');
0665                 fprintf(fd, '\n  #    Mag(pu)   (pu)     (pu)     ($/pu)');
0666                 fprintf(fd, '\n-----  -------  -------  -------  ---------');
0667                 fprintf(fd, '\n%5d%8.3f%9.3f%9.3f%11.3f',...
0668                     [ bus(s.idx(k), BUS_I), bus(s.idx(k),[VM, VMIN]),...
0669                       s.overload(s.idx(k)), ...
0670                       bus(s.idx(k), MU_VMIN)...
0671                     ]');
0672                 fprintf(fd, '\n                        --------');
0673                 fprintf(fd, '\n               Total:%10.2f', ...
0674                         sum(s.overload(s.idx(k))));
0675                 fprintf(fd, '\n');
0676             end
0677         end
0678     end
0679     s = results.softlims.PMAX;
0680     if OUT_PG_LIM && ~strcmp(s.hl_mod,'none')
0681         k = find(s.overload(s.idx) | gen(s.idx, MU_PMAX) > ptol);
0682         fprintf(fd, '\n================================================================================');
0683         fprintf(fd, '\n|     Soft Generator Active Power Upper Bounds                                 |');
0684         fprintf(fd, '\n================================================================================');
0685         if isempty(k)
0686             fprintf(fd,'\nNo violations.\n');
0687         else
0688             fprintf(fd, '\nGen     Bus  Generation  Limit   Overload    mu');
0689             fprintf(fd, '\n  #      #     P (MW)    (MW)     (MW)     ($/MW)');
0690             fprintf(fd, '\n-----  -----  --------  -------  -------  ---------');
0691             fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0692                 [ s.idx(k), gen(s.idx(k),[GEN_BUS, PG, PMAX]),...
0693                   s.overload(s.idx(k)), ...
0694                   gen(s.idx(k), MU_PMAX)...
0695                 ]');
0696             fprintf(fd, '\n                                --------');
0697             fprintf(fd, '\n                       Total:%10.2f', ...
0698                     sum(s.overload(s.idx(k))));
0699             fprintf(fd, '\n');
0700         end
0701     end
0702     s = results.softlims.PMIN;
0703     if OUT_PG_LIM && ~strcmp(s.hl_mod,'none')
0704         k = find(s.overload(s.idx) | gen(s.idx, MU_PMIN) > ptol);
0705         fprintf(fd, '\n================================================================================');
0706         fprintf(fd, '\n|     Soft Generator Active Power Lower Bounds                                 |');
0707         fprintf(fd, '\n================================================================================');
0708         if isempty(k)
0709             fprintf(fd,'\nNo violations.\n');
0710         else
0711             fprintf(fd, '\nGen     Bus  Generation  Limit   Overload    mu');
0712             fprintf(fd, '\n  #      #     P (MW)    (MW)     (MW)     ($/MW)');
0713             fprintf(fd, '\n-----  -----  --------  -------  -------  ---------');
0714             fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0715                 [ s.idx(k), gen(s.idx(k),[GEN_BUS, PG, PMIN]),...
0716                   s.overload(s.idx(k)), ...
0717                   gen(s.idx(k), MU_PMIN)...
0718                 ]');
0719             fprintf(fd, '\n                                --------');
0720             fprintf(fd, '\n                       Total:%10.2f', ...
0721                     sum(s.overload(s.idx(k))));
0722             fprintf(fd, '\n');
0723         end
0724     end
0725     if isAC
0726         s = results.softlims.QMAX;
0727         if OUT_QG_LIM && ~strcmp(s.hl_mod,'none')
0728             k = find(s.overload(s.idx) | gen(s.idx, MU_QMAX) > ptol);
0729             fprintf(fd, '\n================================================================================');
0730             fprintf(fd, '\n|     Soft Generator Reactive Power Upper Bounds                               |');
0731             fprintf(fd, '\n================================================================================');
0732             if isempty(k)
0733                 fprintf(fd,'\nNo violations.\n');
0734             else
0735                 fprintf(fd, '\nGen     Bus  Generation  Limit   Overload    mu');
0736                 fprintf(fd, '\n  #      #    Q (MVAr)  Q (MVAr)  (MVAr)   ($/MVAr)');
0737                 fprintf(fd, '\n-----  -----  --------  -------  -------  ---------');
0738                 fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0739                     [ s.idx(k), gen(s.idx(k),[GEN_BUS, QG, QMAX]),...
0740                       s.overload(s.idx(k)), ...
0741                       gen(s.idx(k), MU_QMAX)...
0742                     ]');
0743                 fprintf(fd, '\n                                --------');
0744                 fprintf(fd, '\n                       Total:%10.2f', ...
0745                         sum(s.overload(s.idx(k))));
0746                 fprintf(fd, '\n');
0747             end
0748         end
0749         s = results.softlims.QMIN;
0750         if OUT_QG_LIM && ~strcmp(s.hl_mod,'none')
0751             k = find(s.overload(s.idx) | gen(s.idx, MU_QMIN) > ptol);
0752             fprintf(fd, '\n================================================================================');
0753             fprintf(fd, '\n|     Soft Generator Reactive Power Lower Bounds                               |');
0754             fprintf(fd, '\n================================================================================');
0755             if isempty(k)
0756                 fprintf(fd,'\nNo violations.\n');
0757             else
0758                 fprintf(fd, '\nGen     Bus  Generation  Limit   Overload    mu');
0759                 fprintf(fd, '\n  #      #    Q (MVAr)  Q (MVAr)  (MVAr)   ($/MVAr)');
0760                 fprintf(fd, '\n-----  -----  --------  -------  -------  ---------');
0761                 fprintf(fd, '\n%5d%7d%9.2f%9.2f%9.2f%11.3f',...
0762                     [ s.idx(k), gen(s.idx(k),[GEN_BUS, QG, QMIN]),...
0763                       s.overload(s.idx(k)), ...
0764                       gen(s.idx(k), MU_QMIN)...
0765                     ]');
0766                 fprintf(fd, '\n                                --------');
0767                 fprintf(fd, '\n                       Total:%10.2f', ...
0768                         sum(s.overload(s.idx(k))));
0769                 fprintf(fd, '\n');
0770             end
0771         end
0772     end
0773     s = results.softlims.RATE_A;
0774     if OUT_LINE_LIM && ~strcmp(s.hl_mod, 'none')
0775         fprintf(fd, '\n================================================================================');
0776         fprintf(fd, '\n|     Soft Branch Flow Limits                                                  |');
0777         fprintf(fd, '\n================================================================================');
0778         k = find(s.overload(s.idx) | sum(branch(s.idx, MU_SF:MU_ST), 2) > ptol);
0779         if isempty(k)
0780             fprintf(fd,'\nNo violations.\n');
0781         else
0782             %% line flow constraints
0783             lim_type = upper(mpopt.opf.flow_lim(1));
0784             if ~isAC || lim_type == 'P' || lim_type == '2'  %% |P| limit
0785                 Ff = branch(s.idx(k), PF);
0786                 Ft = branch(s.idx(k), PT);
0787 %               str = '\n  #     Bus    Bus     (MW)      (MW)      (MW)     ($/MW)';
0788                 str = '\n  #     Bus    Bus        (MW)         (MW)         (MW)        ($/MW)';
0789             elseif lim_type == 'I'                          %% |I| limit
0790                 %% create map of external bus numbers to bus indices
0791                 i2e = bus(:, BUS_I);
0792                 e2i = sparse(max(i2e), 1);
0793                 e2i(i2e) = (1:size(bus, 1))';
0794 
0795                 V = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));
0796                 Ff = abs( (branch(s.idx(k), PF) + 1j * branch(s.idx(k), QF)) ./ V(e2i(branch(s.idx(k), F_BUS))) );
0797                 Ft = abs( (branch(s.idx(k), PT) + 1j * branch(s.idx(k), QT)) ./ V(e2i(branch(s.idx(k), T_BUS))) );
0798 %               str = '\n  #     Bus    Bus     (MA)      (MA)      (MA)     ($/MA)';
0799                 str = '\n  #     Bus    Bus   (kA*basekV)  (kA*basekV)  (kA*basekV)  ($/kA/basekV)';
0800             else                                            %% |S| limit
0801                 Ff = abs(branch(s.idx(k), PF) + 1j * branch(s.idx(k), QF));
0802                 Ft = abs(branch(s.idx(k), PT) + 1j * branch(s.idx(k), QT));
0803 %               str = '\n  #     Bus    Bus     (MVA)     (MVA)     (MVA)    ($/MVA)';
0804                 str = '\n  #     Bus    Bus        (MVA)       (MVA)        (MVA)        ($/MVA)';
0805             end
0806             flow = max(abs(Ff), abs(Ft));
0807 
0808 %           fprintf(fd, '\nBrnch   From   To      Flow      Limit   Overload     mu');
0809             fprintf(fd, '\nBrnch   From   To         Flow        Limit      Overload         mu');
0810             fprintf(fd, str);
0811 %           fprintf(fd, '\n-----  -----  -----  --------  --------  --------  ---------');
0812             fprintf(fd, '\n-----  -----  -----  -----------  -----------  -----------  -------------');
0813 %           fprintf(fd, '\n%4d%7d%7d%10.2f%10.2f%10.2f%11.3f', ...
0814             fprintf(fd, '\n%4d%7d%7d%13.2f%13.2f%13.2f%14.3f', ...
0815                     [   s.idx(k), branch(s.idx(k), [F_BUS, T_BUS]), ...
0816                         flow, branch(s.idx(k), RATE_A), ...
0817                         s.overload(s.idx(k)), ...
0818                         sum(branch(s.idx(k), MU_SF:MU_ST), 2) ...
0819                     ]');
0820             fprintf(fd, '\n                                               -----------');
0821             fprintf(fd, '\n                                      Total:   %10.2f', ...
0822                     sum(s.overload(s.idx(k))));
0823             fprintf(fd, '\n');
0824         end
0825     end
0826     s = results.softlims.ANGMAX;
0827     if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0828         k = find(s.overload(s.idx) | branch(s.idx, MU_ANGMAX) > ptol);
0829         delta = calc_branch_angle(results);
0830         fprintf(fd, '\n================================================================================');
0831         fprintf(fd, '\n|     Soft Maximum Angle Difference Limits                                     |');
0832         fprintf(fd, '\n================================================================================');
0833         if isempty(k)
0834             fprintf(fd,'\nNo violations.\n');
0835         else
0836             fprintf(fd, '\nBrnch   From   To     Angle     Limit    Overload     mu');
0837             fprintf(fd, '\n  #     Bus    Bus    (deg)     (deg)     (deg)     ($/MW)');
0838             fprintf(fd, '\n-----  -----  -----  --------  --------  --------  ---------');
0839             fprintf(fd, '\n%4d%7d%7d%10.3f%10.3f%10.3f%11.3f', ...
0840                 [ s.idx(k), branch(s.idx(k), [F_BUS, T_BUS]), ...
0841                   delta(s.idx(k)), branch(s.idx(k), ANGMAX),...
0842                   s.overload(s.idx(k)),...
0843                   branch(s.idx(k), MU_ANGMAX)...
0844                 ]');
0845             fprintf(fd, '\n                                         --------');
0846             fprintf(fd, '\n                                Total:%10.2f', ...
0847                     sum(s.overload(s.idx(k))));
0848             fprintf(fd, '\n');
0849         end
0850     end
0851     s = results.softlims.ANGMIN;
0852     if OUT_V_LIM && ~strcmp(s.hl_mod,'none')
0853         k = find(s.overload(s.idx) | branch(s.idx, MU_ANGMIN) > ptol);
0854         delta = calc_branch_angle(results);
0855         fprintf(fd, '\n================================================================================');
0856         fprintf(fd, '\n|     Soft Minimum Angle Difference Limits                                     |');
0857         fprintf(fd, '\n================================================================================');
0858         if isempty(k)
0859             fprintf(fd,'\nNo violations.\n');
0860         else
0861             fprintf(fd, '\nBrnch   From   To     Angle     Limit    Overload     mu');
0862             fprintf(fd, '\n  #     Bus    Bus    (deg)     (deg)     (deg)     ($/MW)');
0863             fprintf(fd, '\n-----  -----  -----  --------  --------  --------  ---------');
0864             fprintf(fd, '\n%4d%7d%7d%10.3f%10.3f%10.3f%11.3f', ...
0865                 [ s.idx(k), branch(s.idx(k), [F_BUS, T_BUS]), ...
0866                   delta(s.idx(k)), branch(s.idx(k), ANGMIN),...
0867                   s.overload(s.idx(k)),...
0868                   branch(s.idx(k), MU_ANGMIN)...
0869                 ]');
0870             fprintf(fd, '\n                                         --------');
0871             fprintf(fd, '\n                                Total:%10.2f', ...
0872                     sum(s.overload(s.idx(k))));
0873             fprintf(fd, '\n');
0874         end
0875     end
0876 end
0877 
0878 
0879 %%-----  savecase  -----------------------------------------------------
0880 function mpc = userfcn_softlims_savecase(mpc, fd, prefix, args)
0881 %
0882 %   mpc = userfcn_softlims_savecase(mpc, fd, prefix, args)
0883 %
0884 %   This is the 'savecase' stage userfcn callback that prints the M-file
0885 %   code to save the 'softlims' field in the case file. It expects a
0886 %   MATPOWER case struct (mpc), a file descriptor and variable prefix
0887 %   (usually 'mpc.'). The optional args are not currently used.
0888 
0889 %% structures used to index into softlims in a loop
0890 lims = softlims_lim2mat();
0891 
0892 %% convenience structure for the different fields
0893 slfieldnames = { 'hl_mod', 'hl_val', 'idx', 'busnum', 'cost', ...
0894     'overload', 'ovl_cost' };
0895 fields = struct( ...
0896     'hl_mod',   struct('desc', 'type of hard limit modification', 'tok', '%s'),...
0897     'hl_val',   struct('desc', 'value(s) used to set new hard limit', 'tok', '%g'),...
0898     'busnum',   struct('desc', 'bus numbers for soft voltage limit', 'tok', '%d'),...
0899     'idx',      struct('desc', '%s matrix row indices', 'tok', '%d'),...
0900     'cost',     struct('desc', 'violation cost coefficient', 'tok', '%g'),...
0901     'overload', struct('desc', 'violation of original limit', 'tok', '%0.10g'),...
0902     'ovl_cost', struct('desc', 'overload cost', 'tok', '%.10g') ...
0903 );
0904 
0905 if isfield(mpc, 'softlims')
0906     fprintf(fd, '\n%%%%-----  OPF Soft Limit Data  -----%%%%\n');
0907     limits = fieldnames(lims).';
0908     for lm = limits
0909         lim = lm{:};
0910         s = mpc.softlims.(lim);
0911         if ~strcmp(lim, limits{1})
0912             fprintf(fd, '\n');
0913         end
0914         fprintf(fd,'%%%% %s soft limit data\n', lim);
0915         for f = slfieldnames
0916             f = f{1};
0917             if isfield(s, f)
0918                 if strcmp(f, 'idx')
0919                     if ismember(lim, {'VMIN', 'VMAX'})
0920                         desc = fields.busnum.desc;
0921                     else
0922                         desc = sprintf(fields.idx.desc, lims.(lim).mat_name);
0923                     end
0924                 else
0925                     desc = fields.(f).desc;
0926                 end
0927                 if strcmp(f, 'hl_mod')
0928                     str = sprintf('%ssoftlims.%s.hl_mod = ''%s'';', prefix, lim, s.hl_mod);
0929                     pad = repmat(' ', 1, 40-length(str));
0930                     fprintf(fd, '%s%s%%%% %s\n', str, pad, desc);
0931                 else
0932                     str = sprintf('%ssoftlims.%s.%s = [', prefix, lim, f);
0933                     pad = repmat(' ', 1, 40-length(str));
0934                     fprintf(fd, '%s%s%%%% %s\n', str, pad, desc);
0935                     fprintf(fd, ['\t', fields.(f).tok, ';\n'], s.(f));
0936                     fprintf(fd, '];\n');
0937                 end
0938             end
0939         end
0940     end
0941 end
0942 
0943 
0944 %%-----  softlims_lim2mat  --------------------------------------------
0945 function lim2mat = softlims_lim2mat(mpopt)
0946 
0947 %% define named indices into data matrices
0948 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0949     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0950 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0951     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0952     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0953 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0954     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0955     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0956 
0957 if nargin < 1
0958     mpopt = [];
0959 end
0960 if ~isempty(mpopt) && isfield(mpopt, 'model') && strcmp(mpopt.model, 'DC')
0961     lim2mat = struct(...
0962         'PMAX',   struct('mat_name', 'gen',    'col', PMAX,   'dir',  1 ), ...
0963         'PMIN',   struct('mat_name', 'gen',    'col', PMIN,   'dir', -1 ), ...
0964         'RATE_A', struct('mat_name', 'branch', 'col', RATE_A, 'dir',  1 ), ...
0965         'ANGMAX', struct('mat_name', 'branch', 'col', ANGMAX, 'dir',  1 ), ...
0966         'ANGMIN', struct('mat_name', 'branch', 'col', ANGMIN, 'dir', -1 ) ...
0967     );
0968 else
0969     lim2mat = struct(...
0970         'VMAX',   struct('mat_name', 'bus',    'col', VMAX,   'dir',  1 ), ...
0971         'VMIN',   struct('mat_name', 'bus',    'col', VMIN,   'dir', -1 ), ...
0972         'PMAX',   struct('mat_name', 'gen',    'col', PMAX,   'dir',  1 ), ...
0973         'PMIN',   struct('mat_name', 'gen',    'col', PMIN,   'dir', -1 ), ...
0974         'QMAX',   struct('mat_name', 'gen',    'col', QMAX,   'dir',  1 ), ...
0975         'QMIN',   struct('mat_name', 'gen',    'col', QMIN,   'dir', -1 ), ...
0976         'RATE_A', struct('mat_name', 'branch', 'col', RATE_A, 'dir',  1 ), ...
0977         'ANGMAX', struct('mat_name', 'branch', 'col', ANGMAX, 'dir',  1 ), ...
0978         'ANGMIN', struct('mat_name', 'branch', 'col', ANGMIN, 'dir', -1 ) ...
0979     );
0980 end
0981 
0982 
0983 %%-----  softlims_mat2lims  --------------------------------------------
0984 function mat2lims = softlims_mat2lims(mpopt)
0985 if nargin < 1
0986     mpopt = [];
0987 end
0988 if ~isempty(mpopt) && isfield(mpopt, 'model') && strcmp(mpopt.model, 'DC')
0989     mat2lims = struct(...
0990         'branch', {{'ANGMAX', 'ANGMIN', 'RATE_A'}}, ...
0991         'gen', {{'PMAX', 'PMIN'}} ...
0992     );
0993 else
0994     mat2lims = struct(...
0995         'bus', {{'VMAX', 'VMIN'}}, ...
0996         'branch', {{'ANGMAX', 'ANGMIN', 'RATE_A'}}, ...
0997         'gen', {{'PMAX', 'PMIN', 'QMAX', 'QMIN'}} ...
0998     );
0999 end
1000 
1001 
1002 %%-----  softlims_defaults  --------------------------------------------
1003 function slims = softlims_defaults(mpc, mpopt)
1004 %
1005 %   slims = softlims_defaults(mpc, mpopt)
1006 %
1007 %   Returns a softlims field for mpc in which any missing inputs have
1008 %   been filled in with defaults, so that each limit type has the all
1009 %   4 input fields (idx, cost, hl_mod, hl_cost) and idx has been expanded
1010 %   to a vector, unless hl_mod == 'none', in which case the others need
1011 %   not be present (and are, in any case, ignored).
1012 
1013 %% define named indices into data matrices
1014 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
1015     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
1016 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
1017     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
1018     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
1019 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
1020     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
1021     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
1022 
1023 %% initialization
1024 lims = softlims_lim2mat(mpopt);
1025 if isfield(mpc, 'softlims')
1026     slims = mpc.softlims;
1027 else
1028     slims = struct();
1029 end
1030 if isempty(mpopt)
1031     warning('softlims_defaults: Assuming ''mpopt.opf.softlims.default'' = 1, since mpopt was not provided.');
1032     use_default = 1;
1033 else
1034     use_default = mpopt.opf.softlims.default;
1035 end
1036 
1037 %% get maximum generator marginal cost (of online generators)
1038 max_gen_cost = max(margcost(mpc.gencost, mpc.gen(:, PMAX)));
1039 
1040 %% set defaults for each element
1041 for lm = fieldnames(lims).'
1042     lim = lm{:};
1043     mat  = mpc.order.ext.(lims.(lim).mat_name); %% mpc sub-matrix
1044     col = lims.(lim).col;                       %% mpc sub-matrix column
1045     specified = isfield(slims, lim);
1046     if ~specified && ~use_default
1047         slims.(lim).hl_mod = 'none';
1048     else
1049         if specified
1050             s = slims.(lim);
1051 
1052             %% ignore if hl_mod = 'none'
1053             if isfield(s, 'hl_mod') && strcmp(s.hl_mod, 'none')
1054                 continue;
1055             end
1056         else
1057             s = struct();
1058         end
1059 
1060         %% default base cost is the maximum between a predefined value and
1061         %% 2 times the maximum generation marginal cost. A scaling is applied
1062         %% to account for the fact that we are essentially equating 1MW with
1063         %% 0.01 p.u. voltage or 1 deg angle difference.
1064         default_base_cost = max(1000, 2 * max_gen_cost);
1065         cost_scale = struct( ...
1066             'VMAX', 100, 'VMIN', 100, ...   %% $/pu
1067             'ANGMAX', 1, 'ANGMIN', 1, ...   %% $/deg
1068             'RATE_A', 1, ...                %% $/MW
1069             'PMIN', 1, 'PMAX', 1, ...       %% $/MW
1070             'QMIN', 1, 'QMAX', 1 );         %% $/MVar
1071         default_cost = cost_scale.(lim) * default_base_cost;
1072 
1073         %% idx: indices of bounds to relax
1074         %% idxfull is full list of candidate index values for given constraint
1075         idxfull = softlims_default_idx(mat, lim, col);
1076 
1077         %% convert bus numbers to bus row indices, if necessary, for VMIN/VMAX
1078         if ismember(lim, {'VMAX', 'VMIN'}) && ...
1079                     (~isfield(s, 'idx') || isempty(s.idx)) && ...
1080                     isfield(s, 'busnum') && ~isempty(s.busnum)
1081             s.idx = find(ismember(mat(:, BUS_I), s.busnum));
1082         end
1083 
1084         %% check that idx is a scalar or vector
1085         if isfield(s, 'idx')
1086             if size(s.idx, 2) > 1
1087                 s.idx = s.idx.';        %% make row vector into column vector
1088                 if size(s.idx, 2) > 1
1089                     error('softlim_defaults: mpc.softlims.%s.idx must be a scalar or vector', lim);
1090                 end
1091             end
1092         else    %% set default value for s.idx to full set of constraints
1093             s.idx = idxfull;
1094         end
1095 
1096         %% if there are no constraints to relax, set hl_mod to 'none' and skip
1097         if isempty(s.idx)
1098             s.hl_mod = 'none';
1099             slims.(lim) = s;
1100             continue;
1101         end
1102 
1103         %% cost: cost of violating softlim
1104         if isfield(s, 'cost') && ~isempty(s.cost)
1105             if isscalar(s.cost)     %% expand to vector if specified as scalar
1106                 s.cost = s.cost * ones(size(s.idx));
1107             end
1108         else    %% not specified, use default cost
1109             s.cost = default_cost * ones(size(s.idx));
1110         end
1111 
1112         %% type of limit
1113         ubsign = struct('VMAX', 1 , 'VMIN', -1, 'RATE_A', 1, ...
1114                         'PMAX', 1, 'PMIN', -1, 'QMAX', 1, 'QMIN', -1, ...
1115                         'ANGMAX', 1, 'ANGMIN', -1);
1116         shift_defaults = struct('VMAX', 0.25 , 'VMIN', 0.25, 'RATE_A', 10, ...
1117                         'PMAX', 10, 'PMIN', 10, 'QMAX', 10, 'QMIN', 10, ...
1118                         'ANGMAX', 10, 'ANGMIN', 10);
1119         if isfield(s, 'hl_mod') %% use specified soft limits
1120             switch s.hl_mod
1121                 case 'remove'   %% new hard limit is infinite (hl_val unused)
1122                 case 'scale'    %% new hard limit is original limit * hl_val
1123                     if ~isfield(s, 'hl_val')
1124                         %% use 2 if ubsign and original constraint have
1125                         %% same sign, otherwise use 1/2
1126                         orig_lim = mat(s.idx, col);
1127                         s.hl_val = 2 * ones(size(s.idx));   %% scale up by 2
1128                         k = find(ubsign.(lim) * orig_lim < 0);  %% unless opp. sign
1129                         s.hl_val(k) = 1/2;                  %% then scale down by 2
1130                     end
1131                 case 'shift'    %% new hard limit is original limit + ubsign * hl_val
1132                     if ~isfield(s, 'hl_val')
1133                         s.hl_val = shift_defaults.(lim);
1134                     end
1135                 case 'replace'  %% new hard limit is hl_val (no default value)
1136                 otherwise
1137                     error('softlims_element_defaults: unknown hard limit modification %s', s.hl_mod);
1138             end
1139         else                    %% use default soft limits
1140             switch lim
1141                 case 'PMIN'
1142                     %% for normal gens (PMIN > 0) PMIN is relaxed to 0, not removed
1143                     %% for gens with negative PMIN, it is relaxed to -Inf
1144                     s.hl_mod = 'replace';
1145                     s.hl_val = zeros(size(s.idx));
1146                     s.hl_val(mat(s.idx, PMIN) < 0) = -Inf;
1147                 case 'VMIN'
1148                     %% VMIN is relaxed to zero, not removed
1149                     s.hl_mod = 'replace';
1150                     s.hl_val = 0;
1151                 case {'QMIN', 'ANGMIN'}
1152                     s.hl_mod = 'remove';
1153                     s.hl_val = -Inf;
1154                 otherwise
1155                     s.hl_mod = 'remove';
1156                     s.hl_val = Inf;
1157             end
1158         end
1159 
1160         slims.(lim) = s;
1161     end
1162 end
1163 
1164 
1165 %%-----  softlims_default_idx  --------------------------------------------
1166 function idx = softlims_default_idx(mat, lim, col)
1167 %
1168 %   idx = softlims_default_idx(mat, lim, col)
1169 %
1170 %   Returns the full IDX vector used as the default for the constraint
1171 %   specified by lim.
1172 
1173 %% define named indices into data matrices
1174 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
1175     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
1176 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
1177     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
1178     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
1179 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
1180     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
1181     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
1182 
1183 %% idxfull is the full list of candidate values for idx for given constraint
1184 %% type, all are row indices into EXTERNAL matrix
1185 switch lim
1186     case {'VMAX', 'VMIN'}
1187         %% all buses
1188         idx = [1:size(mat, 1)]';
1189     case {'ANGMAX', 'ANGMIN'}
1190         %% all active branches with meaninful limit (not 0, +360 or -360)
1191         idx = find(mat(:, BR_STATUS) > 0 & mat(:, col) & abs(mat(:, col)) < 360 );
1192     case 'RATE_A'
1193         %% all active branches with flow limit (RATE_A not 0)
1194         idx = find(mat(:, BR_STATUS) > 0 & mat(:, RATE_A) > 0);
1195     case {'PMAX', 'QMAX', 'QMIN'}
1196         %% all active generators (excluding dispatchable loads) w/finite limits
1197         idx = find( (mat(:, GEN_STATUS) > 0) & ~isload(mat) & ~isinf(mat(:, col)) );
1198     case 'PMIN'
1199         %% all active generators (excluding dispatchable loads)
1200         idx = find( (mat(:, GEN_STATUS) > 0) & ~isload(mat) );
1201 end
1202 
1203 
1204 %%-----  softlims_init  --------------------------------------------
1205 function slims = softlims_init(mpc, mpopt)
1206 %
1207 %   slims = softlims_init(mpc, mpopt)
1208 %
1209 %   Returns a softlims field for mpc in which any missing inputs have
1210 %   been filled in with defaults, so that each limit type has the all
1211 %   4 input fields (idx, cost, hl_mod, hl_cost) and idx has been expanded
1212 %   to a vector, unless hl_mod == 'none', in which case the others need
1213 %   not be present (and are, in any case, ignored).
1214 
1215 %% define named indices into data matrices
1216 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
1217     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
1218 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
1219     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
1220     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
1221 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
1222     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
1223     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
1224 
1225 %% initialization
1226 lims = softlims_lim2mat(mpopt);
1227 slims = mpc.softlims;
1228 
1229 %% set defaults for each element
1230 for lm = fieldnames(lims).'
1231     lim = lm{:};
1232     s = slims.(lim);
1233 
1234     %% ignore if hl_mod = 'none'
1235     if isfield(s, 'hl_mod') && strcmp(s.hl_mod, 'none')
1236         continue;
1237     end
1238 
1239     mat  = mpc.order.ext.(lims.(lim).mat_name); %% mpc sub-matrix
1240     col = lims.(lim).col;                       %% mpc sub-matrix column
1241     d = lims.(lim).dir;                         %% constraint direction
1242 
1243     %% idx: indices of bounds to relax
1244     %% idxfull is full list of candidate index values for given constraint
1245     idxfull = softlims_default_idx(mat, lim, col);
1246 
1247     %% idxmask is a boolean vector the size of s.idx, where
1248     %% idxmask(i) is 1 if s.idx(i) is in idxfull and 0 otherwise
1249     idxmask = ismember(s.idx, idxfull);
1250     s.idx = s.idx(idxmask); %% remove possibly irrelevant entries entered by user
1251 
1252     %% if there are no constraints to relax, skip
1253     if isempty(s.idx)
1254         slims.(lim) = s;
1255         continue;
1256     end
1257 
1258     %% save original bounds in s.sav
1259     s.sav = mat(s.idx, col);
1260 
1261     %% cost: cost of violating softlim
1262     %% filter cost to include only valid constraints
1263     try
1264         s.cost = s.cost(idxmask);
1265     catch ME
1266         warning('softlims_init: something went wrong when handling the cost for limit %s. Perhaps the size of the ''cost'' vector didn''t match the size of the ''idx'' vector?', lim);
1267         rethrow(ME);
1268     end
1269 
1270     %% upper bound on constraint violation variable (always a vector)
1271     %% d = 1 for upper bounds, -1 for lower bounds
1272     %% ub = d * ( new limit - original limit )
1273     if isfield(s, 'hl_mod')
1274         switch s.hl_mod
1275             case 'remove'   %% upper bound is infinite
1276                 s.ub = Inf(size(s.idx));
1277             case 'scale'
1278                 %% new limit = original limit * hl_val
1279                 %% ub = d * original limit * (hl_val - 1)
1280                 switch vectorcheck(s, idxmask)
1281                     case 0
1282                         error('softlims_init: provided hl_val vector does not conform to idx vector. When specifying hl_val as a vector idx must also be explicitly given');
1283                     case 1          %% vector
1284                         s.ub = d * s.sav .* (s.hl_val(idxmask) - 1);
1285                     case 2          %% scalar
1286                         s.ub = d * s.sav  * (s.hl_val - 1);
1287                     case {1, 3}     %% default, which is a vector
1288                         s.ub = d * s.sav .* (s.hl_val - 1);
1289                 end
1290             case 'shift'
1291                 %% new limit = original limit + d * hl_val
1292                 %% ub = hl_val
1293                 switch vectorcheck(s, idxmask)
1294                     case 0
1295                         error('softlims_init: provided hl_val vector does not conform to idx vector. When specifying hl_val as a vector idx must also be explicitly given');
1296                     case 1          %% vector
1297                         s.ub = s.hl_val(idxmask);
1298                     case {2, 3}     %% scalar (including default which is a scalar)
1299                         s.ub = s.hl_val * ones(size(s.idx));
1300                 end
1301             case 'replace'
1302                 %% new limit = hl_val
1303                 %% ub = d * ( hl_val - original limit )
1304                 switch vectorcheck(s, idxmask)
1305                     case 0
1306                         error('softlims_init: provided hl_val vector does not conform to idx vector. When specifying hl_val as a vector idx must also be explicitly given');
1307                     case 1          %% vector
1308                         s.ub = d * ( s.hl_val(idxmask) - s.sav );
1309                     case 2          %% scalar
1310                         s.ub = d * ( s.hl_val - s.sav );
1311                     case 3          %% no default for 'replace'
1312                         error('softlims_init: for hard limit ''replace'' modification, replacement value hl_val must be specified');
1313                 end
1314         end
1315         %% check that all ub are non negative (all 'hl_val' were valid)
1316         if any(s.ub < 0)
1317             error('softlims_init: some hl_val for %s results in more restrictive, rather than relaxed, hard constraint.', lim);
1318         end
1319     end
1320 
1321     %% rval: replacement value used to eliminate original hard constraint
1322     switch lim
1323         case {'ANGMAX', 'ANGMIN', 'RATE_A'}
1324             s.rval = 0;
1325         case {'VMAX', 'PMAX', 'QMAX'}
1326             s.rval = Inf;
1327         case {'QMIN', 'PMIN','VMIN'}
1328             s.rval = -Inf;
1329         otherwise
1330             error('softlims_init: limit %s does not have an ''rval'' assigned', lim);
1331     end
1332     slims.(lim) = s;
1333 end
1334 
1335 
1336 %%----- vector check ---------------------------------------------------
1337 function out = vectorcheck(s, idx)
1338 % Utility function for softlims_init(). Checks whether the hl_val
1339 % field in the softlims element struct is a scalar or a vector whose size
1340 % conforms with the idx vector (prior to filtering with idxmask)
1341 % OUTPUT:
1342 %           out      0  hl_val is not scalar and does NOT conform to idx
1343 %                    1  hl_val is not scalar and DOES conform to idx
1344 %                    2  hl_val is a scalar
1345 %                    3  hl_val is not given (use default)
1346 
1347 if ~isfield(s, 'hl_val')
1348     out = 3;
1349 elseif isscalar(s.hl_val)
1350     out = 2;
1351 else
1352     out = all(size(s.hl_val) == size(idx));
1353 end
1354 
1355 
1356 %%-----  softlims_fcn  -------------------------------------------------
1357 function [h, dh] = softlims_fcn(x, mpc, Yf, Yt, il, mpopt, Fmax)
1358 %
1359 %   Evaluates AC branch flow soft limit constraints and Jacobian.
1360 
1361 %% options
1362 lim_type = upper(mpopt.opf.flow_lim(1));
1363 
1364 %% form constraints
1365 %% compute flow (opf.flow_lim = 'P') or square of flow ('S', 'I', '2')
1366 %% note that the Fmax used by opf_branch_flow_fcn() from branch matrix is zero
1367 if nargout == 1
1368     h = opf_branch_flow_fcn(x(1:2), mpc, Yf, Yt, il, mpopt);
1369 else
1370     [h, dh] = opf_branch_flow_fcn(x(1:2), mpc, Yf, Yt, il, mpopt);
1371 end
1372 flv = x{3};     %% flow limit violation variable
1373 if lim_type == 'P'
1374     %%   Ff(Va,Vm) - flv <=  Fmax ===> Ff(Va,Vm) - flv - Fmax <= 0
1375     %%   Ft(Va,Vm) - flv <=  Fmax ===> Ff(Va,Vm) - flv - Fmax <= 0
1376     tmp2 = Fmax + flv;
1377 else    %% lim_type == 'S', 'I', '2'
1378     %%   |Ff(Va,Vm)| - flv <=  Fmax ===> |Ff(Va,Vm)|.^2 - (flv + Fmax).^2 <= 0
1379     %%   |Ft(Va,Vm)| - flv <=  Fmax ===> |Ff(Va,Vm)|.^2 - (flv + Fmax).^2 <= 0
1380     tmp1 = Fmax + flv;
1381     tmp2 = tmp1.^2;
1382 end
1383 h = h - [tmp2; tmp2];
1384 if nargout == 2
1385     ns = length(il);        %% number of soft limits
1386     if lim_type == 'P'
1387         tmp3 = -speye(ns, ns);
1388     else    %% lim_type == 'S', 'I', '2'
1389         tmp3 = spdiags(-2*tmp1, 0, ns, ns);
1390     end
1391     dh = [ dh [tmp3; tmp3] ];
1392 end
1393 
1394 
1395 %%-----  softlims_hess  ------------------------------------------------
1396 function d2H = softlims_hess(x, lambda, mpc, Yf, Yt, il, mpopt)
1397 %
1398 %   Evaluates AC branch flow soft limit constraint Hessian.
1399 
1400 %% options
1401 lim_type = upper(mpopt.opf.flow_lim(1));
1402 
1403 %% form Hessian
1404 d2H = opf_branch_flow_hess(x(1:2), lambda, mpc, Yf, Yt, il, mpopt);
1405 nh = size(d2H, 1);
1406 ns = length(il);        %% number of soft limits
1407 
1408 if lim_type == 'P'
1409     tmp = sparse(ns, ns);
1410 else    %% lim_type == 'S', 'I', '2'
1411     tmp = spdiags(-2*lambda, 0, ns, ns);
1412 end
1413 d2H = [     d2H         sparse(nh, ns);
1414         sparse(ns, nh)      tmp         ];

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