eval_nln_constraint_hess

EVAL_NLN_CONSTRAINT_HESS  Builds and returns Hessian of nonlinear constraints.
   D2G = OM.EVAL_NLN_CONSTRAINT_HESS(X, LAM, ISEQ)
   Builds the Hessian of the full set of nonlinear equality or inequality
   constraints (ISEQ equal to 1 or 0, respectively) for given values of
   the optimization vector X and dual variables LAM, based on constraints
   added by ADD_NLN_CONSTRAINT.

       g(X) = 0
       h(X) <= 0

   Example:
       d2G = om.eval_nln_constraint_hess(x, lam, 1)
       d2H = om.eval_nln_constraint_hess(x, lam, 0)

   See also OPT_MODEL, ADD_NLN_CONSTRAINT, EVAL_NLN_CONSTRAINT.

0001 function d2G = eval_nln_constraint_hess(om, x, lam, iseq)
0002 %EVAL_NLN_CONSTRAINT_HESS  Builds and returns Hessian of nonlinear constraints.
0003 %   D2G = OM.EVAL_NLN_CONSTRAINT_HESS(X, LAM, ISEQ)
0004 %   Builds the Hessian of the full set of nonlinear equality or inequality
0005 %   constraints (ISEQ equal to 1 or 0, respectively) for given values of
0006 %   the optimization vector X and dual variables LAM, based on constraints
0007 %   added by ADD_NLN_CONSTRAINT.
0008 %
0009 %       g(X) = 0
0010 %       h(X) <= 0
0011 %
0012 %   Example:
0013 %       d2G = om.eval_nln_constraint_hess(x, lam, 1)
0014 %       d2H = om.eval_nln_constraint_hess(x, lam, 0)
0015 %
0016 %   See also OPT_MODEL, ADD_NLN_CONSTRAINT, EVAL_NLN_CONSTRAINT.
0017 
0018 %   MP-Opt-Model
0019 %   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
0020 %   by Ray Zimmerman, PSERC Cornell
0021 %
0022 %   This file is part of MP-Opt-Model.
0023 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0024 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0025 
0026 %% get constraint type
0027 if iseq         %% equality constraints
0028     om_nlx = om.nle;
0029 else            %% inequality constraints
0030     om_nlx = om.nli;
0031 end
0032 
0033 %% initialize d2G (use transpose for speed on older versions of MATLAB)
0034 d2Gt = sparse(om.var.N, om.var.N);
0035 
0036 %% calls to substruct() are relatively expensive, so we pre-build the
0037 %% structs for addressing cell and numeric array fields, updating only
0038 %% the subscripts before use
0039 sc = struct('type', {'.', '{}'}, 'subs', {'', 1});  %% cell array field
0040 sn = sc; sn(2).type = '()';                         %% num array field
0041 
0042 %% fill in each piece
0043 for k = 1:om_nlx.NS
0044     name = om_nlx.order(k).name;
0045     idx  = om_nlx.order(k).idx;
0046     if isempty(idx)
0047         if ~isfield(om_nlx.data.hess, name)
0048             continue;   %% skip, there is no function handle stored here,
0049                         %% the function value for this named set was included
0050                         %% in the value computed by a previous named set
0051         end
0052         N = om_nlx.idx.N.(name);    %% number of constraint functions
0053                                     %% evaluated for this named set
0054         if isfield(om_nlx.data.include, name)
0055             N = N + sum(om_nlx.data.include.(name).N);
0056         end
0057     else
0058         % (calls to substruct() are relatively expensive ...
0059         % sn = substruct('.', name, '()', idx);
0060         % sc = substruct('.', name, '{}', idx);
0061         % ... so replace them with these more efficient lines)
0062         sn(1).subs = name;
0063         sn(2).subs = idx;
0064         sc(1).subs = name;
0065         sc(2).subs = idx;
0066         N = subsref(om_nlx.idx.N, sn);
0067     end
0068     if N                                %% non-zero number of rows
0069         if isempty(idx)
0070             d2G_fcn = om_nlx.data.hess.(name);  %% Hessian fcn for kth constraint set
0071             i1 = om_nlx.idx.i1.(name);          %% starting row index
0072             iN = i1 + N - 1;                    %% ending row index
0073             vs = om_nlx.data.vs.(name);         %% var sets
0074         else
0075             d2G_fcn = subsref(om_nlx.data.hess, sc);  %% Hessian fcn for kth constraint set
0076             i1 = subsref(om_nlx.idx.i1, sn);    %% starting row index
0077             iN = subsref(om_nlx.idx.iN, sn);    %% ending row index
0078             vs = subsref(om_nlx.data.vs, sc);   %% var sets
0079         end
0080         xx = om.varsets_x(x, vs);
0081         d2Gk = d2G_fcn(xx, lam(i1:iN));     %% evaluate kth Hessian
0082         
0083         nk = size(d2Gk, 2);
0084         if isempty(vs)          %% all rows of x
0085             if nk == om.var.N
0086                 d2Gkt_full = d2Gk';
0087             else                %% must have added vars since adding
0088                                 %% this constraint set
0089                 d2Gk_all_cols = sparse(nk, om.var.N);
0090                 d2Gk_all_cols(:, 1:nk) = d2Gk;
0091                 d2Gkt_full = sparse(om.var.N, om.var.N);
0092                 d2Gkt_full(:, 1:nk) = d2Gk_all_cols';
0093             end
0094         else                    %% selected rows of x
0095             jj = om.varsets_idx(vs);    %% indices for var set
0096             d2Gk_all_cols = sparse(nk, om.var.N);
0097             d2Gk_all_cols(:, jj) = d2Gk;
0098             d2Gkt_full = sparse(om.var.N, om.var.N);
0099             d2Gkt_full(:, jj) = d2Gk_all_cols';
0100         end
0101         d2Gt = d2Gt + d2Gkt_full;
0102     end
0103 end
0104 d2G = d2Gt';