Hello! I use the following code to build a model via pyomo and check hessian
via IPOPT.
import pyomo.environ as pyo
import pyomo.core.base.external as external
import pyomo.common.gsl as gsl
dll = gsl.find_GSL()
model = pyo.ConcreteModel()
model.test_indices = pyo.RangeSet(0, 2)
x0 = [-0.2, 0.3, 0.1]
model.x = pyo.Var(model.test_indices, domain=pyo.Reals, initialize=lambda m, i:
x0[i])
model.gsl_hypot3 = external.ExternalFunction(library=dll, function="gsl_hypot3")
model.cons = pyo.Constraint(rule=lambda m: (m.gsl_hypot3(m.x[0], m.x[1],
m.x[2]) <= 100))
model.obj = pyo.Objective(rule=lambda m: m.x[0] + m.x[1], m.x[2],
sense=pyo.minimize)
optimizer = pyo.SolverFactory("ipopt")
options["derivative_test"]="second-order"
results = optimizer.solve(model, tee=True)
Then, IPOPT will report warnings like:
Starting derivative checker for second derivatives with obj_factor or lambda[i]
set to 1.5.
\* 1-th constr_hess[ 1, 3] = 1.2351616171153007e-01 v ~ -5.4190095373535634e-02
[ 3.279e+00]
\* 1-th constr_hess[ 2, 2] = -5.4190094682647579e-02 v ~ 1.2351616144287625e-01
[ 1.439e+00]
\* 1-th constr_hess[ 3, 1] = 1.2351616171153007e-01 v ~ -5.4190087606698210e-02
[ 3.279e+00]
After checking the source code of amplgsl_hypot3 function in amplgsl.cc, I
found that the original hessian order is incorrect, which does not accord with
the instruction in funcadd.h: In this case, the function should set al->hes[i +
j*(j+1)/2] to the partial derivative of the function with respect to al->ra[i]
and al->ra[j] for all 0 <= i <= j < al->nr. In amplgsl_hypot3, the terms in
hes[2] and hes[3] should be swaped. I also submit new code at
https://github.com/ampl/gsl/issues/76. May the code be merged?
I'm currently using GSL 2.7 on Windows 10 and obtain it from github. AMPL GSL
is compiled with VS 2022 and default options. The code I use for testing
hessian is writen in Python 3.8, Pyomo 6.6.2, IPOPT 3.14.12. Could you review
the submission and fix the mistake in the Hessian order of amplgsl_hypot3?
Thanks!
