疎な二次計画問題の解法

nAG数値計算ライブラリ > 最適化アルゴリズムExample集 > 疎な二次計画問題の解法

凸二次計画問題-QP（疎）

このExampleでは、疎な二次計画問題を解いています。目的関数は二次形式であり、制約条件は線形不等式で表現されています。問題のサイズはそれほど大きくないですが、ヘッシアン行列と制約条件行列のほとんどの要素がゼロであるため、疎な（スパースな）構造を持っています。このような疎な構造を活用することで、計算効率を高めています。

目的関数:

タスク	式
Minimize	$x_{1}^{2} + x_{2}^{2} + (x_{3} + x_{4})^{2} + x_{5}^{2} + (x_{6} + x_{7})^{2}$

決定変数:

変数	範囲
$x_{1}$	$0 \leq x_{1} \leq 200$
$x_{2}$	$0 \leq x_{2} \leq 2500$
$x_{3}$	$400 \leq x_{3} \leq 800$
$x_{4}$	$100 \leq x_{4} \leq 700$
$x_{5}$	$0 \leq x_{5} \leq 1500$
$x_{6}$	$0 \leq x_{6} < \infty$
$x_{7}$	$0 \leq x_{7} < \infty$

制約条件:

制約	式
1	$0.02 x_{1} + 0.02 x_{2} + 0.03 x_{3} + 1.0 x_{4} + 0.7 x_{5} + 0.02 x_{6} + 0.15 x_{7} \leq 2000$
2	$0.06 x_{2} + 0.75 x_{3} + 0.03 x_{4} + 0.04 x_{5} + 0.05 x_{6} + 0.04 x_{7} \leq 60$
3	$0.02 x_{3} + 1.0 x_{4} + 0.01 x_{5} + 0.08 x_{6} + 0.08 x_{7} \leq 100$
4	$1.0 x_{1} + 0.12 x_{3} + 0.02 x_{4} + 0.02 x_{6} + 0.75 x_{7} \leq 40$
5	$0.01 x_{2} + 0.8 x_{3} + 0.02 x_{5} + 1.0 x_{6} + 0.02 x_{7} \leq 30$
6	$1.0 x_{1} + 0.01 x_{2} + 0.01 x_{4} + 0.97 x_{6} \geq 1500$
7	$0.97 x_{4} + 0.03 x_{7} \geq 250$

このように、7つの決定変数 $x_{1}, \dots, x_{7}$ に対して、それぞれ変数の上下限が設定されています。目的関数は各変数の二乗項の和で表される二次形式です。 7つの線形不等式制約条件が課されており、各制約式は決定変数の一次結合で表現されています。係数行列の多くの要素が0であるため、疎な構造をしていることがわかります。

Exampleの実行コマンド:

python -m naginterfaces.library.examples.opt.qpconvex2_sparse_solve_ex

ソースコード表示コマンド:

python -c "import inspect; from naginterfaces.library.examples.opt import qpconvex2_sparse_solve_ex; print(''.join(inspect.getsourcelines(qpconvex2_sparse_solve_ex)[0]))"

出力結果例:

naginterfaces.library.opt.qpconvex2_sparse_solve Python Example Results.
Sparse QP/LP.
Function value at lowest point found is -1847784.67712.
The corresponding x is:
(0.00, 349.40, 648.85, 172.85, 407.52, 271.36, 150.02).

マニュアル:

qpconvex2_sparse_solveのマニュアル

ソース:

#!/usr/bin/env python3
"``naginterfaces.library.opt.qpconvex2_sparse_solve`` Python Example."

# nAG Copyright 2017-2019.

# pylint: disable=invalid-name,too-many-locals

from naginterfaces.library import opt

def main():
    """
    Example for :func:`naginterfaces.library.opt.qpconvex2_sparse_solve`.

    Sparse QP/LP.

    >>> main()
    naginterfaces.library.opt.qpconvex2_sparse_solve Python Example Results.
    Sparse QP/LP.
    Function value at lowest point found is -1847784.67712.
    The corresponding x is:
    (0.00, 349.40, 648.85, 172.85, 407.52, 271.36, 150.02).
    """

    print(
        'naginterfaces.library.opt.qpconvex2_sparse_solve '
        'Python Example Results.'
    )
    print('Sparse QP/LP.')

    # Cold start:
    start = 'Cold'

    cb_qphx = lambda x, _nstate: (
        [
            2*x[0], 2*x[1],
            2*(x[2] + x[3]), 2*(x[2] + x[3]),
            2*x[4],
            2*(x[5] + x[6]), 2*(x[5] + x[6]),
        ]
    )

    # The sparse A stored by column:
    a = [
        0.02, 0.02, 0.03, 1.00, 0.70, 0.02, 0.15, -200.00,
        0.06, 0.75, 0.03, 0.04, 0.05, 0.04, 1.00, -2000.00,
        0.02, 1.00, 0.01, 0.08, 0.08, 0.80, -2000.00,
        1.00, 0.12, 0.02, 0.02, 0.75, 0.04, -2000.00,
        0.01, 0.80, 0.02, 1.00, 0.02, 0.06, 0.02, -2000.00,
        1.00, 0.01, 0.01, 0.97, 0.01, 400.00,
        0.97, 0.03, 1.00, 400.00,
    ]
    # The index information for the packed A:
    loca = [1, 9, 17, 24, 31, 39, 45, 49]
    inda = [
        7, 5, 3, 1, 6, 4, 2, 8, 7, 6, 5, 4, 3, 2, 1, 8, 2, 1, 4,
        3, 7, 6, 8, 1, 7, 3, 4, 6, 2, 8, 5, 6, 7, 1, 2, 3, 4, 8,
        1, 2, 3, 6, 7, 8, 7, 2, 1, 8,
    ]
    # In this example, we make the objective vector the final column of A:
    c = [0.]
    iobj = 8
    # The number of leading nonzero columns of the Hessian matrix:
    ncolh = 7
    # The bounds on x and on the linear constraints:
    m = 8
    bl = [
        0., 0., 400., 100., 0., 0., 0.,
        2000., -1.e20, -1.e20, -1.e20, -1.e20, 1500., 250., -1.e20,
    ]
    bu = [
        200., 2500., 800., 700., 1500., 1.e20, 1.e20,
        2000., 60., 100., 40., 30., 1.e20, 300., 1.e20,
    ]
    # No constant objective term:
    objadd = 0.
    # Row and column names:
    names = [
        '...X1...', '...X2...', '...X3...', '...X4...', '...X5...',
        '...X6...', '...X7...', '..ROW1..', '..ROW2..', '..ROW3..',
        '..ROW4..', '..ROW5..', '..ROW6..', '..ROW7..', '..COST..',
    ]
    # Blank problem name:
    prob = ''
    # No elastic variables:
    n = len(loca) - 1
    helast = [0]*(n+m)
    # Initial state for cold start:
    hs = [0]*(n+m)
    x = [0.]*(n+m)
    # No superbasics to specify:
    ns = 0

    # Initialize the communication structure:
    comm = opt.qpconvex2_sparse_init()

    qp_soln = opt.qpconvex2_sparse_solve(
        start, m, ncolh, iobj, objadd, prob,
        a, inda, loca, bl, bu, names, helast, hs, x, ns, comm,
        qphx=cb_qphx, c=c,
    )

    print(
        'Function value at lowest point found is {:.5f}.'.format(qp_soln.obj)
    )
    print('The corresponding x is:')
    print('(' + ', '.join(['{:.2f}']*n).format(*qp_soln.x[:n]) + ').')

if __name__ == '__main__':
    import doctest
    import sys
    sys.exit(
        doctest.testmod(
            None, verbose=True, report=False,
            optionflags=doctest.REPORT_NDIFF,
        ).failed
    )