nAG数値計算ライブラリ
> 最適化アルゴリズムExample集
> リバースコミュニケーション法を用いた密な非線形計画問題の解法
非線形計画-NLP（密、Reverse Communication）
このExampleは、Hock and Schittkowski Problem 71と呼ばれる非線形計画問題をリバースコミュニケーション法を用いて解いています。
目的関数:
| タスク | 式 |
|---|---|
| minimize | \(x_0 x_3 (x_0 + x_1 + x_2) + x_2\) |
- 目的関数は、\(x_0\), \(x_1\), \(x_2\), \(x_3\)の4つの決定変数からなる非線形式で表されています。
決定変数:
| 変数 | 範囲 |
|---|---|
| \(x_0\) | \(1 \leq x_0 \leq 5\) |
| \(x_1\) | \(1 \leq x_1 \leq 5\) |
| \(x_2\) | \(1 \leq x_2 \leq 5\) |
| \(x_3\) | \(1 \leq x_3 \leq 5\) |
制約条件:
| 制約 | 式 |
|---|---|
| 制約1 | \(x_0^2 + x_1^2 + x_2^2 + x_3^2 \leq 40\) |
| 制約2 | \(x_0 x_1 x_2 x_3 \geq 25\) |
- 制約条件は2つあり、いずれも決定変数の非線形の不等式で表現されています。
- 制約1は、各決定変数の2乗の和が40以下であることを要求しています。
- 制約2は、4つの決定変数の積が25以上であることを要求しています。
上記の定式化に基づき、本コードはnAGライブラリのnlp1_rcommルーチンを用いて、与えられた初期値から最適解を求めています。目的関数と制約条件の勾配は解析的に計算され、最適化アルゴリズムに渡されます。最終的に得られた目的関数値は1.7014017e+01となります。
Exampleの実行コマンド:
python -m naginterfaces.library.examples.opt.nlp1_rcomm_exソースコード表示コマンド:
python -c "import inspect; from naginterfaces.library.examples.opt import nlp1_rcomm_ex; print(''.join(inspect.getsourcelines(nlp1_rcomm_ex)[0]))"出力結果例:
naginterfaces.library.opt.nlp1_rcomm Python Example Results.
Solve Hock and Schittkowski Problem 71.
Final objective value is 1.7014017e+01
マニュアル:
ソース:
#!/usr/bin/env python3
"``naginterfaces.library.opt.nlp1_rcomm`` Python Example."
# nAG Copyright 2017-2019.
# pylint: disable=invalid-name,too-many-locals
import numpy as np
from naginterfaces.library import opt
def main():
"""
Example for :func:`naginterfaces.library.opt.nlp1_rcomm`.
Dense NLP.
>>> main()
naginterfaces.library.opt.nlp1_rcomm Python Example Results.
Solve Hock and Schittkowski Problem 71.
Final objective value is 1.7014017e+01
"""
print(
'naginterfaces.library.opt.nlp1_rcomm Python Example Results.'
)
print('Solve Hock and Schittkowski Problem 71.')
# Initialize the solver:
comm = opt.nlp1_init('nlp1_rcomm')
# The initial guess:
x = np.array([1., 5., 5., 1.])
# The linear constraints:
a = np.array([[1.]*len(x)])
# There are two nonlinear constraints defined by cb_confun:
ncnln = 2
# The bounds:
bl = [1., 1., 1., 1., -1.0E+25, -1.0E+25, 25.]
bu = [5., 5., 5., 5., 20., 40., 1.0E+25]
n = len(x)
istate = np.zeros(n + 1 + ncnln, dtype=int)
c = np.empty(max(1, ncnln))
cjac = np.zeros((ncnln, n))
clamda = np.zeros(n + 1 + ncnln)
objgrd = np.empty(n)
r = np.zeros((n, n))
irevcm = 0
itera = 0
objf = 0.
needc = np.empty(ncnln, dtype=int)
while True:
irevcm, itera, objf, needc = opt.nlp1_rcomm(
irevcm, a.shape[0], a, bl, bu, itera,
istate, c, cjac, clamda, objf, objgrd, r, x,
comm,
)
if irevcm == 0:
break
if irevcm in [1, 3]:
objf = x[0]*x[3]*(x[0] + x[1] + x[2]) + x[2]
if irevcm in [2, 3]:
objgrd[:] = [
x[3]*(2*x[0] + x[1] + x[2]),
x[0]*x[3],
x[0]*x[3] + 1.0,
x[0]*(x[0] + x[1] + x[2]),
]
if irevcm in [4, 6]:
if needc[0] > 0:
c[0] = (x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2)
if needc[1] > 0:
c[1] = x[0]*x[1]*x[2]*x[3]
if irevcm in [5, 6]:
if needc[0] > 0:
cjac[0, :] = 2*x[:]
if needc[1] > 0:
cjac[1, :] = [
x[1]*x[2]*x[3],
x[0]*x[2]*x[3],
x[0]*x[1]*x[3],
x[0]*x[1]*x[2],
]
print('Final objective value is {:.7e}'.format(objf))
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(
None, verbose=True, report=False,
optionflags=doctest.REPORT_NDIFF,
).failed
)