最適化アルゴリズムExample集: 疎な非線形計画問題の内点法による解法
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.

マニュアル:

handle_solve_ipoptのマニュアル

ソース:

#!/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:
    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=int)
    icolh = np.empty(10, dtype=int)
    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
    )

関連情報
Privacy Policy  /  Trademarks