#!/usr/bin/env python3
"``naginterfaces.library.opt.handle_solve_ipopt`` Python Example."
# NAG Copyright 2018-2020.
# pylint: disable=invalid-name,too-many-arguments,too-many-locals,too-many-statements
from math import sqrt
import numpy as np
from naginterfaces.library import opt
from naginterfaces.base import utils as b_utils
[docs]def main():
"""
Example for :func:`naginterfaces.library.opt.handle_solve_ipopt`.
Interior-point solver for sparse NLP.
>>> main()
naginterfaces.library.opt.handle_solve_ipopt Python Example Results.
Solving a problem based on Hock and Schittkowski Problem 73.
Solving with a nonlinear objective.
At the solution the objective function is 2.9894378e+01.
"""
print(
'naginterfaces.library.opt.handle_solve_ipopt Python Example Results.'
)
print('Solving a problem based on Hock and Schittkowski Problem 73.')
print('Solving with a nonlinear objective.')
# The 'non'linear objective:
cb_objfun = lambda x, inform: (
24.55*x[0] + 26.75*x[1] + 39.*x[2] + 40.5*x[3], inform,
)
def cb_objgrd(x, fdx, inform):
fdx[:] = [24.55, 26.75, 39., 40.5]
return inform
cb_confun = lambda x, _ncnln, inform: (
[
12.0*x[0] + 11.9*x[1] + 41.8*x[2] + 52.1*x[3] - 1.645 *
sqrt(
0.28*x[0]**2 + 0.19*x[1]**2 + 20.5*x[2]**2 + 0.62*x[3]**2
)
],
inform,
)
def cb_congrd(x, gdx, inform):
"""The Jacobian of the nonlinear constraints."""
tmp = sqrt(
0.62*x[3]**2 + 20.5*x[2]**2 + 0.19*x[1]**2 + 0.28*x[0]**2
)
gdx[:] = [
(12.0*tmp-0.4606*x[0])/tmp,
(11.9*tmp-0.31255*x[1])/tmp,
(41.8*tmp-33.7225*x[2])/tmp,
(52.1*tmp-1.0199*x[3])/tmp,
]
return inform
def cb_hess(x, idf, sigma, lamda, hx, inform):
"""Hessian function."""
nnzh = hx.size
if idf == 0:
# Objective is linear.
for i in range(nnzh):
hx[i] = 0.0
return inform
tmp = sqrt(
0.62*x[3]**2 + 20.5*x[2]**2 + 0.19*x[1]**2 + 0.28*x[0]**2
)
tmp = tmp*(
x[3]**2 + 33.064516129032258064*x[2]**2 +
0.30645161290322580645*x[1]**2 + 0.45161290322580645161*x[0]**2
)
hx[0] = (
-0.4606*x[3]**2 - 15.229516129032258064*x[2]**2 -
0.14115161290322580645*x[1]*x[1])/tmp
hx[1] = (0.14115161290322580645*x[0]*x[1])/tmp
hx[2] = (15.229516129032258064*x[0]*x[2])/tmp
hx[3] = (0.4606*x[0]*x[3])/tmp
hx[4] = (
-0.31255*x[3]**2 - 10.334314516129032258*x[2]**2 -
0.14115161290322580645*x[0]**2
)/tmp
hx[5] = (10.334314516129032258*x[1]*x[2])/tmp
hx[6] = (0.31255*x[1]*x[3])/tmp
hx[7] = (
-33.7225*x[3]**2 - 10.334314516129032258*x[1]**2 -
15.229516129032258065*x[0]**2
)/tmp
hx[8] = (33.7225*x[2]*x[3])/tmp
hx[9] = (
-33.7225*x[2]**2 - 0.31255*x[1]**2 - 0.4606*x[0]**2
)/tmp
if idf == -1:
for i in range(nnzh):
hx[i] = hx[i] * lamda[0]
else:
assert idf == 1
return inform
# The initial guess:
x = [1., 1., 1., 1.]
nvar = len(x)
nnzu = 2*nvar
# Create a handle for the problem:
handle = opt.handle_init(nvar)
# Define the bounds:
opt.handle_set_simplebounds(
handle,
bl=[0.]*nvar, bu=[1.e20]*nvar,
)
# Define the non-linear objective:
opt.handle_set_nlnobj(
handle,
idxfd=[1, 2, 3, 4],
)
# Define the linear constraints:
opt.handle_set_linconstr(
handle,
bl=[5., 1.], bu=[1.e20, 1.],
irowb=[1, 1, 1, 1, 2, 2, 2, 2],
icolb=[1, 2, 3, 4, 1, 2, 3, 4],
b=[2.3, 5.6, 11.1, 1.3, 1.0, 1.0, 1.0, 1.0],
)
nnzu += 2*2
# Define the nonlinear constraints:
opt.handle_set_nlnconstr(
handle,
bl=[21.], bu=[1.e20],
irowgd=[1, 1, 1, 1], icolgd=[1, 2, 3, 4],
)
nnzu += 2*1
# Define the structure of the Hessian, dense upper triangle:
irowh = np.empty(10, dtype=b_utils.EngineIntNumPyType)
icolh = np.empty(10, dtype=b_utils.EngineIntNumPyType)
idx = 0
for i in range(1, nvar + 1):
for j in range(i, nvar + 1):
icolh[idx] = j
irowh[idx] = i
idx += 1
# For the objective function:
opt.handle_set_nlnhess(
handle,
idf=0, irowh=irowh, icolh=icolh,
)
# For the constraint function:
opt.handle_set_nlnhess(
handle,
idf=1, irowh=irowh, icolh=icolh,
)
opt.handle_opt_set(
handle, 'Stop Tolerance 1 = 2.5e-8',
)
opt.handle_opt_set(
handle, 'Print Level = 0',
)
# Solve the problem:
ip_soln = opt.handle_solve_ipopt(
handle, x,
objfun=cb_objfun, objgrd=cb_objgrd, confun=cb_confun, congrd=cb_congrd,
hess=cb_hess,
)
print(
'At the solution the objective function is {:.7e}.'.format(
ip_soln.rinfo[0]
)
)
# Destroy the handle:
opt.handle_free(handle)
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(
None, verbose=True, report=False,
optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
).failed
)
