nAG数値計算ライブラリ
> 最適化アルゴリズムExample集
> 疎な非線形計画問題の内点法による解法
非線形計画法-NLP(疎、内点法)
このExampleは、Hock と Schittkowski の問題73に基づく非線形計画問題を、内点法を用いて解いています。問題は疎な形式で定式化され、非線形の目的関数と線形・非線形の制約条件を含んでいます。このExampleの目的は、naginterfaces.library.opt.handle_solve_ipoptを用いて非線形計画問題を解く方法を示すことです。
目的関数:
タスク | 式 |
---|---|
minimize | \(24.55x_0 + 26.75x_1 + 39x_2 + 40.5x_3\) |
決定変数:
変数 | 範囲 |
---|---|
\(x_0\) | \([0, \infty)\) |
\(x_1\) | \([0, \infty)\) |
\(x_2\) | \([0, \infty)\) |
\(x_3\) | \([0, \infty)\) |
制約条件:
制約 | 式 |
---|---|
線形制約1 | \(2.3x_0 + 5.6x_1 + 11.1x_2 + 1.3x_3 \geq 5\) |
線形制約2 | \(x_0 + x_1 + x_2 + x_3 = 1\) |
非線形制約 | \(12.0x_0 + 11.9x_1 + 41.8x_2 + 52.1x_3 - 1.645 \sqrt{0.28x_0^2 + 0.19x_1^2 + 20.5x_2^2 + 0.62x_3^2} \geq 21\) |
決定変数は非負の実数で、\(x_0, x_1, x_2,
x_3\) の4変数です。
目的関数は、各変数に重み付けされた線形関数の最小化です。
線形制約は2本あり、1本目は各変数の重み付き和が5以上という不等式制約、2本目は各変数の和が1に等しいという等式制約です。
非線形制約は1本で、各変数の重み付き和から変数の2乗和の平方根にスケーリングされた値を引いたものが21以上という不等式制約です。
cb_hess関数は、非線形制約のヘッシアン行列の非ゼロ要素を計算しています。
Exampleの実行コマンド:
python -m naginterfaces.library.examples.opt.handle_solve_ipopt_ex
ソースコード表示コマンド:
python -c "import inspect; from naginterfaces.library.examples.opt import handle_solve_ipopt_ex; print(''.join(inspect.getsourcelines(handle_solve_ipopt_ex)[0]))"
出力結果例:
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.
マニュアル:
ソース:
#!/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
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:
= lambda x, inform: (
cb_objfun 24.55*x[0] + 26.75*x[1] + 39.*x[2] + 40.5*x[3], inform,
)
def cb_objgrd(x, fdx, inform):
= [24.55, 26.75, 39., 40.5]
fdx[:] return inform
= lambda x, _ncnln, inform: (
cb_confun
[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."""
= sqrt(
tmp 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."""
= hx.size
nnzh if idf == 0:
# Objective is linear.
for i in range(nnzh):
= 0.0
hx[i] return inform
= sqrt(
tmp 0.62*x[3]**2 + 20.5*x[2]**2 + 0.19*x[1]**2 + 0.28*x[0]**2
)= tmp*(
tmp 3]**2 + 33.064516129032258064*x[2]**2 +
x[0.30645161290322580645*x[1]**2 + 0.45161290322580645161*x[0]**2
)0] = (
hx[-0.4606*x[3]**2 - 15.229516129032258064*x[2]**2 -
0.14115161290322580645*x[1]*x[1])/tmp
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] = (
hx[-0.31255*x[3]**2 - 10.334314516129032258*x[2]**2 -
0.14115161290322580645*x[0]**2
/tmp
)5] = (10.334314516129032258*x[1]*x[2])/tmp
hx[6] = (0.31255*x[1]*x[3])/tmp
hx[7] = (
hx[-33.7225*x[3]**2 - 10.334314516129032258*x[1]**2 -
15.229516129032258065*x[0]**2
/tmp
)8] = (33.7225*x[2]*x[3])/tmp
hx[9] = (
hx[-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] * lamda[0]
hx[i] else:
assert idf == 1
return inform
# The initial guess:
= [1., 1., 1., 1.]
x = len(x)
nvar = 2*nvar
nnzu # Create a handle for the problem:
= opt.handle_init(nvar)
handle
# Define the bounds:
opt.handle_set_simplebounds(
handle,=[0.]*nvar, bu=[1.e20]*nvar,
bl
)
# Define the non-linear objective:
opt.handle_set_nlnobj(
handle,=[1, 2, 3, 4],
idxfd
)
# Define the linear constraints:
opt.handle_set_linconstr(
handle,=[5., 1.], bu=[1.e20, 1.],
bl=[1, 1, 1, 1, 2, 2, 2, 2],
irowb=[1, 2, 3, 4, 1, 2, 3, 4],
icolb=[2.3, 5.6, 11.1, 1.3, 1.0, 1.0, 1.0, 1.0],
b
)+= 2*2
nnzu
# Define the nonlinear constraints:
opt.handle_set_nlnconstr(
handle,=[21.], bu=[1.e20],
bl=[1, 1, 1, 1], icolgd=[1, 2, 3, 4],
irowgd
)+= 2*1
nnzu
# Define the structure of the Hessian, dense upper triangle:
= np.empty(10, dtype=int)
irowh = np.empty(10, dtype=int)
icolh = 0
idx for i in range(1, nvar + 1):
for j in range(i, nvar + 1):
= j
icolh[idx] = i
irowh[idx] += 1
idx # For the objective function:
opt.handle_set_nlnhess(
handle,=0, irowh=irowh, icolh=icolh,
idf
)# For the constraint function:
opt.handle_set_nlnhess(
handle,=1, irowh=irowh, icolh=icolh,
idf
)
opt.handle_opt_set('Stop Tolerance 1 = 2.5e-8',
handle,
)
opt.handle_opt_set('Print Level = 0',
handle,
)
# Solve the problem:
= opt.handle_solve_ipopt(
ip_soln
handle, x,=cb_objfun, objgrd=cb_objgrd, confun=cb_confun, congrd=cb_congrd,
objfun=cb_hess,
hess
)
print(
'At the solution the objective function is {:.7e}.'.format(
0]
ip_soln.rinfo[
)
)
# Destroy the handle:
opt.handle_free(handle)
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(None, verbose=True, report=False,
=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
optionflags
).failed
)