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