Description of toggle

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         ];
toggle_softlims

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
