nlp_hessfcn

NLP_HESSFCN  Evaluates Hessian of Lagrangian.
   LXX = NLP_HESSFCN(OM, X, LAMBDA, COST_MULT)
   LXX = NLP_HESSFCN(OM, X, LAMBDA, COST_MULT, HS)

   Hessian evaluation function, suitable for use with MIPS, FMINCON, etc.

   Inputs:
     OM : Opt-Model object
     X : optimization vector
     LAMBDA (struct)
       .eqnonlin : Lagrange multipliers on equality constraints
       .ineqnonlin : Kuhn-Tucker multipliers on inequality constraints
     COST_MULT : (optional) Scale factor to be applied to the cost
          (default = 1).
     HS : (optional) sparse matrix with tiny non-zero values specifying
          the fixed sparsity structure that the resulting LXX should match

   Outputs:
     LXX : Hessian of the Lagrangian.

   Examples:
       Lxx = nlp_hessfcn(om, x, lambda, cost_mult);
       Lxx = nlp_hessfcn(om, x, lambda, cost_mult, Hs);

   See also NLP_COSTFCN, NLP_CONSFCN.

0001 function Lxx = nlp_hessfcn(om, x, lambda, cost_mult, Hs)
0002 %NLP_HESSFCN  Evaluates Hessian of Lagrangian.
0003 %   LXX = NLP_HESSFCN(OM, X, LAMBDA, COST_MULT)
0004 %   LXX = NLP_HESSFCN(OM, X, LAMBDA, COST_MULT, HS)
0005 %
0006 %   Hessian evaluation function, suitable for use with MIPS, FMINCON, etc.
0007 %
0008 %   Inputs:
0009 %     OM : Opt-Model object
0010 %     X : optimization vector
0011 %     LAMBDA (struct)
0012 %       .eqnonlin : Lagrange multipliers on equality constraints
0013 %       .ineqnonlin : Kuhn-Tucker multipliers on inequality constraints
0014 %     COST_MULT : (optional) Scale factor to be applied to the cost
0015 %          (default = 1).
0016 %     HS : (optional) sparse matrix with tiny non-zero values specifying
0017 %          the fixed sparsity structure that the resulting LXX should match
0018 %
0019 %   Outputs:
0020 %     LXX : Hessian of the Lagrangian.
0021 %
0022 %   Examples:
0023 %       Lxx = nlp_hessfcn(om, x, lambda, cost_mult);
0024 %       Lxx = nlp_hessfcn(om, x, lambda, cost_mult, Hs);
0025 %
0026 %   See also NLP_COSTFCN, NLP_CONSFCN.
0027 
0028 %   MP-Opt-Model
0029 %   Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC)
0030 %   by Ray Zimmerman, PSERC Cornell
0031 %
0032 %   This file is part of MP-Opt-Model.
0033 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0034 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0035 
0036 %% ----- evaluate d2f -----
0037 [f, df, d2f] = nlp_costfcn(om, x);
0038 d2f = d2f * cost_mult;
0039 
0040 %%----- evaluate Hessian of equality constraints -----
0041 d2G = om.eval_nln_constraint_hess(x, lambda.eqnonlin, 1);
0042 
0043 %%----- evaluate Hessian of inequality constraints -----
0044 d2H = om.eval_nln_constraint_hess(x, lambda.ineqnonlin, 0);
0045 
0046 %%-----  do numerical check using (central) finite differences  -----
0047 if 0
0048     nx = length(x);
0049     step = 1e-5;
0050     num_d2f = sparse(nx, nx);
0051     num_d2G = sparse(nx, nx);
0052     num_d2H = sparse(nx, nx);
0053     for i = 1:nx
0054         xp = x;
0055         xm = x;
0056         xp(i) = x(i) + step/2;
0057         xm(i) = x(i) - step/2;
0058         % evaluate cost & gradients
0059         [fp, dfp] = nlp_costfcn(om, xp);
0060         [fm, dfm] = nlp_costfcn(om, xm);
0061         % evaluate constraints & gradients
0062         [Hp, Gp, dHp, dGp] = nlp_consfcn(om, xp);
0063         [Hm, Gm, dHm, dGm] = nlp_consfcn(om, xm);
0064         num_d2f(:, i) = cost_mult * (dfp - dfm) / step;
0065         num_d2G(:, i) = (dGp - dGm) * lambda.eqnonlin   / step;
0066         num_d2H(:, i) = (dHp - dHm) * lambda.ineqnonlin / step;
0067     end
0068     d2f_err = full(max(max(abs(d2f - num_d2f))));
0069     d2G_err = full(max(max(abs(d2G - num_d2G))));
0070     d2H_err = full(max(max(abs(d2H - num_d2H))));
0071     if d2f_err > 1e-6
0072         fprintf('Max difference in d2f: %g\n', d2f_err);
0073     end
0074     if d2G_err > 1e-5
0075         fprintf('Max difference in d2G: %g\n', d2G_err);
0076     end
0077     if d2H_err > 1e-6
0078         fprintf('Max difference in d2H: %g\n', d2H_err);
0079     end
0080 end
0081 
0082 Lxx = d2f + d2G + d2H;
0083 
0084 %% force specified sparsity structure
0085 if nargin > 4
0086     %% add sparse structure (with tiny values) to current matrices to
0087     %% ensure that sparsity structure matches that supplied
0088     Lxx = Lxx + Hs;
0089 
0090 %     %% check sparsity structure against that supplied
0091 %     if nnz(Lxx) ~= nnz(Hs)
0092 %         fprintf('=====> nnz(Lxx) is %d, expected %d <=====\n', nnz(Lxx), nnz(Hs));
0093 %     else
0094 %         [iHs, jHs] = find(Hs);
0095 %         [iH, jH] = find(Lxx);
0096 %         if any(iH ~= iHs) || any(jH ~= jHs)
0097 %             fprintf('=====> structure of Lxx is not as expected <=====\n');
0098 %         end
0099 %     end
0100 end