params_quad_cost

PARAMS_QUAD_COST  Returns the cost parameters for quadratic costs.
   [Q, C] = OM.PARAMS_QUAD_COST()
   [Q, C] = OM.PARAMS_QUAD_COST(NAME)
   [Q, C] = OM.PARAMS_QUAD_COST(NAME, IDX_LIST)
   [Q, C, K] = OM.PARAMS_QUAD_COST(...)
   [Q, C, K, VS] = OM.PARAMS_QUAD_COST(...)

   With no input parameters, it assembles and returns the parameters
   for the aggregate quadratic cost from all quadratic cost sets added
   using ADD_QUAD_COST. The values of these parameters are cached
   for subsequent calls. The parameters are Q, C, and optionally K,
   where the quadratic cost is of the form
       F(X) = 1/2 * X'*Q*X + C'*X + K

   If a NAME is provided then it simply returns the parameters for the
   corresponding named set. Likewise for indexed named sets specified
   by NAME and IDX_LIST. In this case, Q and K may be vectors, corresponding
   to a cost function of the form
       F(X) = 1/2 * Q .* X.^2 + C .* X + K

   An optional 4th output argument VS indicates the variable sets used by
   this cost set. The size of Q and C will be consistent with VS.

   See also OPT_MODEL, ADD_QUAD_COST.

0001 function [Q, c, K, vs] = params_quad_cost(om, name, idx)
0002 %PARAMS_QUAD_COST  Returns the cost parameters for quadratic costs.
0003 %   [Q, C] = OM.PARAMS_QUAD_COST()
0004 %   [Q, C] = OM.PARAMS_QUAD_COST(NAME)
0005 %   [Q, C] = OM.PARAMS_QUAD_COST(NAME, IDX_LIST)
0006 %   [Q, C, K] = OM.PARAMS_QUAD_COST(...)
0007 %   [Q, C, K, VS] = OM.PARAMS_QUAD_COST(...)
0008 %
0009 %   With no input parameters, it assembles and returns the parameters
0010 %   for the aggregate quadratic cost from all quadratic cost sets added
0011 %   using ADD_QUAD_COST. The values of these parameters are cached
0012 %   for subsequent calls. The parameters are Q, C, and optionally K,
0013 %   where the quadratic cost is of the form
0014 %       F(X) = 1/2 * X'*Q*X + C'*X + K
0015 %
0016 %   If a NAME is provided then it simply returns the parameters for the
0017 %   corresponding named set. Likewise for indexed named sets specified
0018 %   by NAME and IDX_LIST. In this case, Q and K may be vectors, corresponding
0019 %   to a cost function of the form
0020 %       F(X) = 1/2 * Q .* X.^2 + C .* X + K
0021 %
0022 %   An optional 4th output argument VS indicates the variable sets used by
0023 %   this cost set. The size of Q and C will be consistent with VS.
0024 %
0025 %   See also OPT_MODEL, ADD_QUAD_COST.
0026 
0027 %   MP-Opt-Model
0028 %   Copyright (c) 2017-2020, Power Systems Engineering Research Center (PSERC)
0029 %   by Ray Zimmerman, PSERC Cornell
0030 %
0031 %   This file is part of MP-Opt-Model.
0032 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0033 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0034 
0035 if nargin > 1       %% individual set
0036     if nargin < 3
0037         idx = {};
0038     end
0039     if isempty(idx)                 %% name, no index provided
0040         if numel(om.qdc.idx.i1.(name)) == 1     %% simple named set
0041             Q = om.qdc.data.Q.(name);
0042             c = om.qdc.data.c.(name);
0043             K = om.qdc.data.k.(name);
0044             if nargout > 3
0045                 vs = om.qdc.data.vs.(name);
0046             end
0047         else                                    %% indexing required
0048             error('@opt_model/params_quad_cost: quadratic cost set ''%s'' requires an IDX_LIST arg', name);
0049         end
0050     else                            %% indexed named set
0051         % (calls to substruct() are relatively expensive ...
0052         % s = substruct('.', name, '{}', idx);
0053         % ... so replace it with these more efficient lines)
0054         sc = struct('type', {'.', '{}'}, 'subs', {name, idx});
0055         Q = subsref(om.qdc.data.Q, sc);
0056         c = subsref(om.qdc.data.c, sc);
0057         K = subsref(om.qdc.data.k, sc);
0058         if nargout > 3
0059             vs = subsref(om.qdc.data.vs, sc);
0060         end
0061     end
0062 else                %% aggregate
0063     cache = om.qdc.params;
0064     if isempty(cache)       %% build the aggregate
0065         nx = om.var.N;          %% number of variables
0066         if om.qdc.NS < 25 || om.qdc.NS < 100 && nx < 300
0067             %% METHOD 1: Add sparse matrices (original method)
0068             Qt = sparse(nx, nx);    %% transpose of quadratic coefficients
0069             c = zeros(nx, 1);       %% linear coefficients
0070             K = 0;                  %% constant term
0071             for k = 1:om.qdc.NS
0072                 name = om.qdc.order(k).name;
0073                 idx  = om.qdc.order(k).idx;
0074                 N = om.getN('qdc', name, idx);
0075                 [Qk, ck, kk, vs] = om.params_quad_cost(name, idx);
0076                 haveQ = ~isempty(Qk);
0077                 havec = ~isempty(ck);
0078                 nk = max(size(Qk, 1), size(ck, 1));     %% size of Qk and/or ck
0079                 if isempty(vs)
0080                     if nk == nx     %% full size
0081                         if size(Qk, 2) == 1     %% Qk is a column vector
0082                             Qkt_full = spdiags(Qk, 0, nx, nx);
0083                         elseif haveQ            %% Qk is a matrix
0084                             Qkt_full = Qk';
0085                         end
0086                         if havec
0087                             ck_full = ck;
0088                         end
0089                     else            %% vars added since adding this cost set
0090                         if size(Qk, 2) == 1     %% Qk is a column vector
0091                             Qkt_full = sparse(1:nk, 1:nk, Qk, nx, nx);
0092                         elseif haveQ            %% Qk is a matrix
0093                             Qk_all_cols = sparse(nk, nx);
0094                             Qk_all_cols(:, 1:nk) = Qk;
0095                             Qkt_full(:, 1:nk) = Qk_all_cols';
0096                         end
0097                         if havec
0098                             ck_full = zeros(nx, 1);
0099                             ck_full(1:nk) = ck;
0100                         end
0101                     end
0102                 else
0103                     jj = om.varsets_idx(vs);    %% indices for var set
0104                     if size(Qk, 2) == 1     %% Qk is a column vector
0105                         Qkt_full = sparse(jj, jj, Qk, nx, nx);
0106                     elseif haveQ            %% Qk is a matrix
0107                         Qk_all_cols = sparse(nk, nx);
0108                         Qk_all_cols(:, jj) = Qk;
0109                         Qkt_full = sparse(nx, nx);
0110                         Qkt_full(:, jj) = Qk_all_cols';
0111                     end
0112                     if havec
0113                         ck_full = zeros(nx, 1);
0114                         ck_full(jj) = ck;
0115                     end
0116                 end
0117                 if haveQ
0118                     Qt = Qt + Qkt_full;
0119                 end
0120                 if havec
0121                     c = c + ck_full;
0122                 end
0123                 if length(kk) == 1
0124                     K = K + N * kk;     %% N handles case where k is expanded to vector cost
0125                 else
0126                     K = K + sum(kk);
0127                 end
0128             end
0129             Q = Qt';
0130        else
0131             %% METHOD 2: construct using single call to sparse()
0132             Q_ijv = cell(om.qdc.NS, 3); %% indices/values to construct Q
0133             c_ijv = cell(om.qdc.NS, 3); %% indices/values to construct c
0134             K = 0;                  %% constant term
0135             for k = 1:om.qdc.NS
0136                 name = om.qdc.order(k).name;
0137                 idx  = om.qdc.order(k).idx;
0138                 N = om.getN('qdc', name, idx);
0139                 [Qk, ck, kk, vs] = om.params_quad_cost(name, idx);
0140                 haveQ = ~isempty(Qk);
0141                 havec = ~isempty(ck);
0142                 if haveQ
0143                     [i, j, v] = find(Qk);
0144                 end
0145                 if havec
0146                     [ic, jc, vc] = find(ck);
0147                 end
0148                 if isempty(vs)
0149                     if size(Qk, 2) == 1     %% Qk is a column vector
0150                         Q_ijv(k, :) = {i, i, v};
0151                     elseif haveQ            %% Qk is a matrix
0152                         Q_ijv(k, :) = {i, j, v};
0153                     end
0154                     if havec
0155                         c_ijv(k, :) = {ic, jc, vc};
0156                     end
0157                 else
0158                     jj = om.varsets_idx(vs)';    %% indices for var set
0159                     if size(Qk, 2) == 1     %% Qk is a column vector
0160                         Q_ijv(k, :) = {jj(i), jj(i), v};
0161                     elseif haveQ            %% Qk is a matrix
0162                         Q_ijv(k, :) = {jj(i), jj(j), v};
0163                     end
0164                     if havec
0165                         c_ijv(k, :) = {jj(ic), jc, vc};
0166                     end
0167                 end
0168                 if length(kk) == 1
0169                     K = K + N * kk;     %% N handles case where k is expanded to vector cost
0170                 else
0171                     K = K + sum(kk);
0172                 end
0173             end
0174             Q = sparse( vertcat(Q_ijv{:,1}), ...
0175                         vertcat(Q_ijv{:,2}), ...
0176                         vertcat(Q_ijv{:,3}), nx, nx);
0177             c = accumarray(vertcat(c_ijv{:,1}), vertcat(c_ijv{:,3}), [nx 1]);
0178         end
0179 
0180         %% cache aggregated parameters
0181         om.qdc.params = struct('Q', Q, 'c', c, 'k', K);
0182     else                    %% return cached values
0183         Q = cache.Q;
0184         c = cache.c;
0185         K = cache.k;
0186     end
0187     if nargout > 3
0188         vs = {};
0189     end
0190 end