数理最適化 スタートガイド

NAG Library for Python で始める数理最適化

NAG数値計算ライブラリ > NAG数理最適化ソルバー > 数理最適化 スタートガイド

数理最適化 スタートガイド

イントロダクション

このドキュメントは、NAG Library for Python を用いて、最適化問題を解きたい方々向けのスタートガイドです。

最適化問題とは、ある制約条件の下で、与えられた目的関数の値を最小化（または最大化）するような変数の値を求める問題のことを指します。例えば、コストを最小限に抑えつつ利益を最大化するような生産計画を立てる問題などが該当します。

このドキュメントでは、NAG Library for Python のインストール、ご利用方法、ドキュメント参照方法、などをExampleコードを交えながらステップを追って説明します。

インストール

NAG Library for PythonはWindows 64-bit、Linux 64-bit、Intel Mac 64-bitのいずれかでご利用いただけます。

以下の手順に従ってインストールを行って下さい。(全プラットフォーム共通）

尚、インストールについてのより詳細な情報はNAG Library for Pythonインストールガイドをご参照下さい。

インストールコマンド

ソフトウエア本体及びドキュメント類は、コマンドラインから以下のコマンドでインストールする事ができます。

python -m pip install --extra-index-url https://www.nag.com/downloads/py/naginterfaces_mkl naginterfaces

ライセンスキーのインストール

NAG Library for Python の利用にはライセンスキーが必要です。

ライセンスキーはテキスト形式のものか、もしくはUSB キーライセンスとなります。

テキスト形式のライセンスキーは「nag.key」という名前でホームディレクトリに保存してください。

~/nag.key

USB キーライセンスの場合はご利用のマシンに USB キーを装着してください。

トライアルライセンスご希望の方

ライセンスキーをお持ちでない方で、トライアルをご希望の方は こちら からお申込み下さい。お申し込みが受け付けられるとメールでテキスト形式のトライアルライセンスファイル「nag.key」が送られてきますので、それをホームディレクトリに保存して下さい。

インストールの確認方法

インストールが上手くいったかどうかの確認は以下のコマンドで行っていただけます。

python -c "from naginterfaces import quick_check; quick_check()"

何も問題が無い場合は、各種情報が表示され、最後に以下のように ‘all ok’ と表示されます。

QUICK CHECK: all ok

何か問題があった場合（例えばライセンスが見つかれない等）はその旨のエラーメッセージが表示されます。

アンインストール

アンインストールは以下のコマンドで行って下さい。

python -m pip uninstall naginterfaces

実践チュートリアル：NAG Library for Pythonを使った最適化問題の解決

このセクションでは、NAG Library for Pythonを使って最適化問題を解くための具体的な手順を示します。

  • （１）基本的な最適化問題の解決（Rosenbrock最適化問題）：
    • まずは簡単な例を使って、NAG Libraryの基本的な使い方が学べます。
  • （２）線形計画問題（薬剤配合の最適化）：
    • 線形の最適化問題を解決する方法を紹介します。線形計画問題を解くための関数 lp_solveの使い方が学べます。
  • （３）非線形計画問題（カーボンフットプリント削減問題）：
    • 非線形の最適化問題を解決する方法を紹介します。非線形計画問題を解くための関数 nlp1_solveの使い方が学べます。
  • （４）混合整数計画問題（予算制約下仕入れ最適化問題）：
    • 整数制約を含んだ混合整数計画問題を解決する方法を紹介します。混合整数計画問題を解くための関数 handle_solve_milpの使い方が学べます。

（１）基本的な最適化問題の解決（Rosenbrock最適化問題）

まずは簡単な例を使って、NAG Libraryの基本的な使い方を学びます。 ここでは最適化ソルバーのインポート、目的関数の定義、初期値の設定、ソルバーの呼び出し、そして結果の表示方法についての基本的な流れを学びます。

Rosenbrock最適化問題

Rosenbrock関数は最適化のテストでよく用いられ、最小値を見つけるのが難しい関数です。

ここでは、NAGライブラリ for Pythonで提供される bounds_quasi_func_easy を利用して、一般化Rosenbrock関数（以下）の最小化を行います。bounds_quasi_func_easy は、変数の上限と下限の制約条件下で関数の最小値を求めるための関数です。

\(f(x)=\sum_{i=1}^{N-1}[100(x_{i+1}-x_i^2)^2+(1-x_i)^2]\)

この問題はN=4の場合、点[1,1,1,1]で最小値0に到達する事が知られています。

手順

まず、最適化ソルバー関数をインポートします。

from naginterfaces.library.opt import bounds_quasi_func_easy

次に目的関数を定義します。目的関数とは、最適化したい関数のことです。ここではラムダ式を用いてRosenbrock関数を定義していますが、通常の関数定義を用いても問題ありません。ここで定義した関数は、後でソルバーに引数として渡します。

rosen = lambda x: (
    sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1.0-x[:-1])**2.0)
)

次に初期値を与えます。初期値は、最適化アルゴリズムが解の探索を開始する点のことで、ここでは[1.2, 1.0, 1.2, 1.0]を与えています。最適化アルゴリズムは、この初期値から始めて、徐々に解を改善していきます。

x = [1.2, 1.0, 1.2, 1.0]

次に変数の上下限制約を設定します。上下限制約とは、変数の取りうる値の範囲を制限する制約条件の一種です。ここでは、各変数の下限を[0, 0, 0, 0]、上限を[2, 2, 2, 2]としています。iboundは上下限制約の種類を指定するためのパラメータで、今回は各変数の下限と上限を個別に指定しているため、ibound=0としています。

n = len(x)
bl, bu = ([0.0]*n, [2.0]*n)
ibound = 0

ソルバーを呼び出して問題を解きます。bounds_quasi_func_easy関数に、先ほど定義した目的関数や初期値、上下限制約を渡して最適化を実行します。

opt_soln = bounds_quasi_func_easy(
    ibound, rosen, bl, bu, x,
)

得られた結果を表示します。opt_soln.fは最適化後の目的関数の最小値（最適値）を表し、opt_soln.xは最適解における変数の値（最適解）を表しています。

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

最適解は(1, 1, 1, 1)、最適値は0であることが知られているため、正しく最適解が求まっていることが確認できました。

まとめ

NAG Lirbraryを使用した最適化問題の解決プロセスは、以下の基本的な要素で構成されています：

  1. 必要なソルバー関数をインポートする

  2. 以下の要素を定義する（順不同）：

    • 最適化したい目的関数
    • 変数の初期値
    • 問題に応じた制約条件（例：変数の上下限）

  3. ソルバー関数を呼び出し、定義した情報を引数として渡す

  4. ソルバーから返された結果（最適解と最適値）を表示・解析する

この基本的な構造は、NAG Libraryの他の最適化関数にも適用可能です。関数や制約条件を適切に定義し、適切なソルバーを選択することで、線形計画問題、非線形計画問題、混合整数計画問題など、様々な最適化問題に対応できます。次のセクションでは、これらの異なるタイプの最適化問題に対するNAG Libraryの適用例を見ていきます。

（２）線形計画問題（薬剤配合の最適化）

線形計画問題の実践的な応用例として、高血圧治療に用いられる3種類の薬剤A、B、Cの配合比を最適化する問題を考えてみましょう。線形計画問題は、目的関数と制約条件がすべて線形である最適化問題のことを指します。

ここでは、高血圧治療に用いられる3種類の薬剤A、B、Cの配合比を最適化する問題を考えます。 各薬剤には治療効果と副作用があり、それぞれの用量によって変化します。

  • 治療効果は、各薬剤の用量（mg）に治療効果係数を乗じたものの合計で表されます。治療効果係数は、1mgあたりの治療効果の大きさを表す定数です。

  • 副作用は、各薬剤の用量（mg）に副作用係数を乗じたものの合計で表されます。副作用係数は、1mgあたりの副作用の大きさを表す定数です。

以下は問題を解くために必要なパラメータ値です：

薬剤 治療効果係数 副作用係数 最小用量 (mg) 最大用量 (mg)
A 0.8 0.3 10 50
B 1.2 0.5 20 60
C 1.0 0.4 15 45

副作用の許容値は 60 とします。

例えば、薬剤Aを30mg、薬剤Bを40mg、薬剤Cを20mg使用する場合、以下のようになります：

  • 治療効果 = 0.8 × 30 + 1.2 × 40 + 1.0 × 20 = 92
  • 副作用 = 0.3 × 30 + 0.5 × 40 + 0.4 × 20 = 37

この場合、副作用の合計値は37であり、許容値の60以下となっています。

この問題の目的は、副作用の合計値が許容値以下となる条件の下で、治療効果を最大化するような各薬剤の最適な用量の組み合わせを見つけることです。 つまり、副作用を一定の範囲内に抑えながら、最大限の治療効果が得られる薬剤の配合比を求めることです。

目的関数:

目的 式
治療効果の最大化 \(\max \sum_{i=1}^{3} e_i x_i\)

決定変数:

変数 説明 範囲
\(x_i\) 薬剤 \(i\) の用量 (mg) \(l_i \leq x_i \leq u_i\)

制約条件:

制約 式 説明
副作用の許容値 \(\sum_{i=1}^{3} s_i x_i \leq S_{max}\) 副作用係数と各薬剤の用量の積の総和が許容値以下でなければならない
各薬剤の用量範囲 \(l_i \leq x_i \leq u_i, \quad i = 1, 2, 3\) 各薬剤の用量は最小値と最大値の範囲内でなければならない

コード例: 以下に、この線形計画問題を NAG Library for Python の LP 問題専用ソルバー関数 lp_solve を用いて解くコード例を示します。今回の問題は変数と制約条件の数が比較的少なく、小規模な線形計画問題であるため、汎用的な LP ソルバーである lp_solve を選択しました。

import numpy as np
from naginterfaces.library import opt

# パラメータ設定
# 各薬剤の治療効果係数と副作用係数
e = np.array([0.8, 1.2, 1.0])  # 治療効果係数
s = np.array([0.3, 0.5, 0.4])  # 副作用係数

# 各薬剤の用量範囲 (最小用量、最大用量)
l = np.array([10, 20, 15])  # 最小用量 (mg)
u = np.array([50, 60, 45])  # 最大用量 (mg)

# 副作用の許容最大値
S_max = 60

# 目的関数の係数（治療効果を最大化するため、負の符号を用いて最小化問題に変換）
cvec = -e

# 変数と制約の上下限
bl = np.concatenate([l, np.array([-np.inf])])  # 下限値 (x_i の下限と副作用の合計の下限)
bu = np.concatenate([u, np.array([S_max])])    # 上限値 (x_i の上限と副作用の合計の上限)

# 制約係数行列（副作用制約）
a = np.array([s])

# 初期解の設定
x_init = np.array([10, 20, 15])  # 初期用量

# コミュニケーションオブジェクトの初期化 ('lp_solve'関数用)
comm = opt.nlp1_init('lp_solve')

# 線形プログラミング問題を解く
_, x, _, obj, _, _ = opt.lp_solve(bl, bu, x_init, comm, a, cvec)

# 結果の表示
print(f"最適解：薬剤Aの用量 = {x[0]:.2f} mg, 薬剤Bの用量 = {x[1]:.2f} mg, 薬剤Cの用量 = {x[2]:.2f} mg")
print(f"最大治療効果：{-obj:.2f}")

# 副作用の合計値を計算
total_side_effects = np.dot(s, x)
print(f"副作用値：{total_side_effects:.2f}")

コードの詳細説明

この例題では、NAG Library for Pythonの非線形計画ソルバー lp_solve を用いて、薬剤配合問題を解いています。以下では、コードの各部分について詳しく説明します。

目的関数の設定

cvec = -e

この部分で目的関数を設定しています。

lp_solve 関数は、目的関数が以下のような線形（または1次）の形式であることを前提としています：

\(f(x) = c_1x_1 + c_2x_2 + ... + c_nx_n\)

ここで、\(c_i\) は係数、\(x_i\) は決定変数です。この形式が決まっているため、係数のみを渡せば目的関数を完全に定義できます。

今回の問題の目的関数は次のようになります：

\(\max f(x) = 0.8x_1 + 1.2x_2 + 1.0x_3\)

この式の係数ベクトルは [0.8, 1.2, 1.0] となり、これが e に相当します。cvec = -e として符合を反転させているのは、lp_solveが最小化問題と解く設計のソルバーだからです。 目的関数の係数の符号を反転させることにより、最大化問題を等価な最小化問題に変換しています。

変数と制約の上下限の設定

bl = np.concatenate([l, np.array([-np.inf])])  # 下限値 (x_i の下限と副作用の合計の下限)
bu = np.concatenate([u, np.array([S_max])])    # 上限値 (x_i の上限と副作用の合計の上限)

この部分では、変数（各薬剤の用量）と制約（副作用の合計）の上下限を設定しています。lp_solve では、変数の上下限と線形制約の上下限を一つのベクトルで指定する必要があります。このベクトルの長さは、変数の数と線形制約の数の和に等しくなります。

この問題では、変数は3つの薬剤の用量 \(x_1, x_2, x_3\) であり、線形制約は副作用の上限に関する1つの制約のみです。したがって、bl と bu のベクトルの長さは \(3 + 1 = 4\) となります。

bl は下限値を表し、各薬剤の最小用量 l と副作用の合計の下限（負の無限大）を結合したベクトルです。bu は上限値を表し、各薬剤の最大用量 u と副作用の合計の上限 S_max を結合したベクトルです。

具体的には、以下のようになります：

bl = [10, 20, 15, -∞]
bu = [50, 60, 45, 60]

制約係数行列の設定

a = np.array([s])

ここでは、線形制約条件の係数行列 a を設定しています。この問題では、副作用の制約条件 \(\sum_{i=1}^{3} s_i x_i \leq S_{max}\) のみ存在します。

重要なポイントは、lp_solve が制約条件の係数を常に2次元配列として受け取ることです。そのため、例えば今回のように制約が1つだけであっても、以下のように1行の2次元配列として定義する必要があります：

a = [[0.3, 0.5, 0.4]]

この行列は不等式制約 \(0.3x_1 + 0.5x_2 + 0.4x_3 \leq S_{max}\) の左辺の係数を表します。右辺の \(S_{max}\) と不等号 \(\leq\) の情報は、別のパラメータで指定します。

具体的に、lp_solve 関数では、すべての制約条件に対して bl (lower bound) と bu (upper bound) の両方の値を設定することで、右辺と不等号の情報を表現します。今回の \(\leq\) 制約の場合：

  • bl の対応する要素に負の無限大（実際には非常に小さな値）を設定
  • bu の対応する要素に \(S_{max}\) の値を設定

これにより、\(0.3x_1 + 0.5x_2 + 0.4x_3 \leq S_{max}\) という不等式制約全体を表現できます。

なお、複数の制約がある場合は、行列 a に複数の行が存在することになり、それぞれの制約に対応した bl と bu の要素を適切に設定する必要があります。

初期値の設定とソルバーの呼び出し

x_init = np.array([10, 20, 15])  # 初期用量
comm = opt.nlp1_init('lp_solve')
_, x, _, obj, _, _ = opt.lp_solve(bl, bu, x_init, comm, a, cvec)

ここでは、まず初期値 x_init を設定しています。初期値は、最適化アルゴリズムが解の探索を開始する点となります。

次に、opt.nlp1_init('lp_solve') を呼び出して、lp_solve 関数の必須引数であるコミュニケーションオブジェクト comm を初期化しています。 今回の例では単純に初期化しているだけですが、複雑な問題設定や他のNAGルーチンとの連携などの場合に利用されます。

最後に、opt.lp_solve 関数を呼び出して線形計画問題を解いています。この関数には、以下の引数を渡しています：

引数 説明
bl 変数と制約の下限値を表すベクトル
bu 変数と制約の上限値を表すベクトル
x_init 初期値を表すベクトル
comm コミュニケーションオブジェクト
a 線形制約の係数行列
cvec 目的関数の係数ベクトル

opt.lp_solve 関数は複数の値を返しますが、ここでは x（最適解を表すベクトル）と obj（最適値）のみを利用しています。

結果の表示と副作用の合計値の計算

print(f"最適解：薬剤Aの用量 = {x[0]:.2f} mg, 薬剤Bの用量 = {x[1]:.2f} mg, 薬剤Cの用量 = {x[2]:.2f} mg")
print(f"最大治療効果：{-obj:.2f}")

total_side_effects = np.dot(s, x)
print(f"副作用値：{total_side_effects:.2f}")

最後に、最適解 x と最適値 obj を表示しています。また、最適解 x と副作用係数ベクトル s の内積を計算することで、副作用の合計値 total_side_effects を求め、表示しています。

結果例

上記のコードを実行すると、以下のような結果が得られます：

最適解：薬剤Aの用量 = 50.00 mg, 薬剤Bの用量 = 54.00 mg, 薬剤Cの用量 = 45.00 mg
最大治療効果：149.80
副作用値：60.00

この結果から、最適な薬剤の配合比は薬剤Aを50mg、薬剤Bを54mg、薬剤Cを45mg使用することであり、このときの治療効果は149.80、副作用値は許容値の60ちょうどになることがわかります。

まとめ

このセクションでは、線形計画問題を解くために以下の主要な点を学びました：

  1. NAG Library の線形計画問題用ソルバーの一つであるlp_solve の利用方法
  2. 目的関数、制約条件、および変数の上下限のコード上での表現方法

（３）非線形計画問題（カーボンフットプリント削減問題）

次に、非線形の最適化問題を取り上げます。この問題では、限られた予算内で複数のCO2削減対策に最適な投資配分を行い、CO2削減量を最大化することを目指します。

ある企業が、年間のCO2排出量を削減するために、以下の3つの対策を検討しています：

  1. 再生可能エネルギーの導入
  2. 省エネ設備への更新
  3. 森林保全への投資

それぞれの対策には、投資額に応じたCO2削減効果がありますが、その関係は非線形です。企業は、限られた予算内で、これらの対策への最適な投資配分を決定し、CO2削減量を最大化したいと考えています。

対策 投資額の範囲（百万円） CO2削減効果（トン）
再生可能エネルギー 0 - 50 \(100 \sqrt{x_1}\)
省エネ設備 0 - 30 \(50 \sqrt{x_2}\)
森林保全 0 - 20 \(80 \log(1 + x_3)\)
CO2削減対策の投資額とCO2削減効果の関係

総予算は6000万円とします。

目的関数:

目的 式
CO2削減量の最大化 \(\max\ 100 \sqrt{x_1} + 50 \sqrt{x_2} + 80 \log(1 + x_3)\)

決定変数:

変数 説明 範囲
\(x_1\) 再生可能エネルギーへの投資額（百万円） \(0 \leq x_1 \leq 50\)
\(x_2\) 省エネ設備への投資額（百万円） \(0 \leq x_2 \leq 30\)
\(x_3\) 森林保全への投資額（百万円） \(0 \leq x_3 \leq 20\)
  • \(x_1, x_2, x_3\): それぞれの対策への投資額を表す非負の連続値です。単位は百万円です。

制約条件:

制約 式 説明
予算制約 \(x_1 + x_2 + x_3 \leq B\) 総投資額が予算以内に収まる
投資額の非負制約 \(x_1, x_2, x_3 \geq 0\) 各対策への投資額は非負
投資額の上限制約 \(x_1 \leq 50, x_2 \leq 30, x_3 \leq 20\) 各対策への投資額は上限以下

制約条件は、総投資額が予算以内に収まること、各対策への投資額が非負であること、そして各対策への投資額がそれぞれの上限以下であることを表しています。

コード例:

以下に、この問題を NAG Library for Python の非線形最適化問題を解くための関数 nlp1_solve を用いて解くコード例を示します。

import numpy as np
from naginterfaces.library import opt

# パラメータ定義
budget = 60  # 総予算 (百万円), 予算制約のBです。
limits = {"renew": 50, "save": 30, "conserve": 20}  # 各対策の投資上限 (百万円), 各対策への投資額の上限制約を表します。
coefs = {"renew": 100, "save": 50, "conserve": 80}  # 各対策のCO2削減係数, CO2削減効果の非線形関係を定義します。
epsilon = 1e-8  # 微小値, 数値的安定性を保つために使用されます。

def objective(mode, x, grad, _nstate):
    """ 目的関数：CO2削減量最大化 """
    if mode in [0, 2]:
        # 目的関数の計算, 各投資額に対するCO2削減効果の合計値を計算します。
        obj_val = -(
            coefs["renew"] * np.sqrt(x[0]) +
            coefs["save"] * np.sqrt(x[1]) +
            coefs["conserve"] * np.log(1 + x[2])
        )
    else:
        obj_val = 0.0
    if mode == 2:
        # 勾配の計算, 最適化のための勾配ベクトルを計算します。
        grad[:] = [
            -coefs["renew"] / (2 * np.sqrt(x[0]) + epsilon),
            -coefs["save"] / (2 * np.sqrt(x[1]) + epsilon),
            -coefs["conserve"] / (1 + x[2] + epsilon)
        ]
    return obj_val, grad

# 投資範囲の設定
lower_bounds_var = [0, 0, 0]  # 各対策の投資額の下限です。
upper_bounds_var = [limits["renew"], limits["save"], limits["conserve"]]  # 上限制約です。

# 線形制約
a = np.full((1, 3), 1.0, dtype=np.float64)
lower_bounds_con = [ -1.0e20 ]
upper_bounds_con = [ budget ]

# nlp1_solveの仕様に合わせて変数の上下限と制約の上下限を結合
lower_bounds = np.concatenate((lower_bounds_var, lower_bounds_con))
upper_bounds = np.concatenate((upper_bounds_var, upper_bounds_con))

# 各対策への投資額の初期値
x_init = [15, 15, 15]  

# オプション設定
comm = opt.nlp1_init('nlp1_solve')
opt.nlp1_option_string('Print Level = 0', comm)

# 最適化実行
result = opt.nlp1_solve(
    a, lower_bounds, upper_bounds, objective, x_init, comm
)

# 結果表示
print('最適な投資額:')
print(f'再生可能エネルギーへの投資 (renew): {result.x[0]:.2f} 百万円')
print(f'省エネルギーへの投資 (save): {result.x[1]:.2f} 百万円') 
print(f'森林保全への投資 (conserve): {result.x[2]:.2f} 百万円')
print(f'達成されたCO2削減量の最大値: {-result.objf:.2f} トン')

コードの詳細説明

この例題では、NAG Library for Pythonの非線形計画ソルバー nlp1_solve を用いて、カーボンフットプリント削減問題を解いています。以下では、コードの各部分について詳しく説明します。

目的関数の設定

nlp1_solveでは、ユーザーが目的関数をコールバック関数(objective)として定義します。この関数は、modeに応じて目的関数の値 (obj_val) とその勾配 (grad) を計算し、返すように記述します。mode=0または2の場合は目的関数値を、mode=2の場合は勾配も計算します。

def objective(mode, x, grad, _nstate):
    """ 目的関数：CO2削減量最大化 """
    if mode in [0, 2]:
        # 目的関数の計算, 各投資額に対するCO2削減効果の合計値を計算します。
        obj_val = -(
            coefs["renew"] * np.sqrt(x[0]) +
            coefs["save"] * np.sqrt(x[1]) +
            coefs["conserve"] * np.log(1 + x[2])
        )
    else:
        obj_val = 0.0
    if mode == 2:
        # 勾配の計算, 最適化のための勾配ベクトルを計算します。
        grad[:] = [
            -coefs["renew"] / (2 * np.sqrt(x[0]) + epsilon),
            -coefs["save"] / (2 * np.sqrt(x[1]) + epsilon),
            -coefs["conserve"] / (1 + x[2] + epsilon)
        ]
    return obj_val, grad

今回の目的は、各対策への投資額に応じたCO2削減効果の合計の最大化です。

つまり、目的関数は次のようになります：

最大化：100 * sqrt(x1) + 50 * sqrt(x2) + 80 * log(1 + x3)

nlp1_solve関数は最小化問題を解くものであるため、最大化を行いたい場合は最小化問題に変換する必要があります。それを行うには目的関数の符号を反転させます。つまり、次のようにします：

最小化：-(100 * sqrt(x1) + 50 * sqrt(x2) + 80 * log(1 + x3))

符号を反転させて最小化することにより、結果的に元の目的関数を最大化することになります。

この関数では目的関数の勾配（偏微分）も計算させる必要があります。

今回の目的関数の勾配（偏微分）は以下のようになります。

  • d(100 * sqrt(x1)) / dx1（100 * sqrt(x1)のx1に関する微分）= 50 / sqrt(x1)
  • d(50 * sqrt(x2)) / dx2（50 * sqrt(x2)のx2に関する微分）= 25 / sqrt(x2)
  • d(80 * log(1 + x3)) / dx3（80 * log(1 + x3)のx3に関する微分）= 80 / (1 + x3)

勾配も符号の反転が必要なので、以下のように符合を反転させてgradに指定しています：

  • grad[0] = -100 / (2 * sqrt(x[0]) + epsilon)
  • grad[1] = -50 / (2 * sqrt(x[1]) + epsilon)
  • grad[2] = -80 / (1 + x[2] + epsilon)

（ここでepsilon を足しているのは、ゼロ除算や無限大になることを防ぐためです。これにより、数値計算が安定します。）

このように、目的関数とその勾配をコード内で直接定義することで、ソルバーに数式を渡すことができます。

制約条件の設定（線形制約）

# 投資範囲の設定
lower_bounds_var = [0, 0, 0]  # 各対策の投資額の下限です。
upper_bounds_var = [limits["renew"], limits["save"], limits["conserve"]]  # 上限制約です。

# 総予算の線形制約
a = np.full((1, 3), 1.0, dtype=np.float64)
lower_bounds_con = [-1.0e20]
upper_bounds_con = [budget]

# nlp1_solveの仕様に合わせて変数の上下限と制約の上下限を結合
lower_bounds = np.concatenate((lower_bounds_var, lower_bounds_con))
upper_bounds = np.concatenate((upper_bounds_var, upper_bounds_con))

ここでは、各対策への投資額に対する上下限と総予算の線形制約を設定しています。

  • lower_bounds_varとupper_bounds_varは、それぞれの変数（各対策の投資額）の範囲を設定しています。lower_bounds_var = [0, 0, 0]は各投資額の下限を0に設定し、投資が0以上であることを示しています。upper_bounds_var = [50, 30, 20]は各投資額の上限を設定し、再生可能エネルギーへの投資が50、省エネ設備への投資が30、森林保全への投資が20を超えないことを示しています。
  • 線形制約は、行列aで係数を指定し、lower_bounds_conとupper_bounds_conで上下限を設定します。

今回の総予算の制約式は次のようになっています：

x1 + x2 + x3 <= 60

この式の各係数（今回のケースでは全て1）を係数行列aに以下のように指定しています。

a = np.full((1, 3), 1.0)

ここでlp_solveは複数の制約式の係数を二次元配列aで与えるようになっているため、その仕様に合わせ１ｘ３の行列としています。（今回は1つの制約式があり、３つの変数があるため）

そして、制約式の上限はupper_bounds_con = [60]、下限はlower_bounds_con = [-1.0e20]と与え、各対策の投資額の合計が60を超えないことを指定しています。

変数と制約の上下限の指定方法

nlp1_solveでは、変数と制約の下限と上限を、lower_boundsとupper_boundsの引数にまとめて指定します。これらの配列の大きさは、変数の数 + 線形制約の数 + 非線形制約の数であり、この順番で値を並べる必要があります。

今回の問題では、変数（投資額）が3つ、線形制約（予算制約）が1つあるので、lower_boundsとupper_boundsの配列の大きさは3 + 1 = 4となります。具体的には、以下のように変数の上下限と線形制約の上下限を結合してそれを行っています。

lower_bounds = np.concatenate((lower_bounds_var, lower_bounds_con))
upper_bounds = np.concatenate((upper_bounds_var, upper_bounds_con))

ソルバーの呼び出しと結果の取得

以上で定義した目的関数(objective)、制約式(a)、変数と制約の上下限(lower_bounds, upper_bounds)などを引数としてnlp1_solveを呼び出すことで、最適化計算を実行しています。その他、変数の初期値(x_init)や通信用のオブジェクト(comm)なども引数として渡します。

result = opt.nlp1_solve(
    a, lower_bounds, upper_bounds, objective, x_init, comm
)

nlp1_solve関数は最適化の結果を含むオブジェクトを返します。このオブジェクトから最適解を取得しています。

結果例:

上記のコードを実行すると、以下のような出力が得られます。

最適な投資額:
再生可能エネルギーへの投資 (renew): 40.64 百万円
省エネルギーへの投資 (save): 10.16 百万円
森林保全への投資 (conserve): 9.20 百万円
達成されたCO2削減量の最大値: 982.66 トン
CO2削減対策の投資額とCO2削減効果の関係

上記の結果から、提示された予算60百万円の範囲内で、再生可能エネルギーへの投資に40.64百万円、省エネルギーへの投資に10.16百万円、森林保全への投資に9.20百万円を配分することが最適であることがわかります。この投資配分により、CO2削減量は982.66トンに達します。

まとめ

このセクションでは、以下の主要な点を学びました：

  1. NAG Library の nlp1_solve 関数を用いた非線形最適化問題の解き方
  2. nlp1_solve 関数における目的関数と勾配のコールバック関数としての実装方法
  3. nlp1_solve 関数における変数の上下限と線形制約条件の設定方法

（４）混合整数計画問題（予算制約下の仕入れ最適化問題）

次に、整数制約を含む混合整数計画問題を取り上げます。

ある物流企業の在庫管理者が、予算10万円で以下の製品から最大の利益を得られるように仕入れを行いたいと考えています。各製品の仕入れ価格と利益は以下の通りです。最大の利益を得るために、どの製品を何個仕入れるべきか決定してください。

製品 仕入れ価格（円） 利益（円）
A 1,000 700
B 6,000 5,000
C 11,000 9,000
D 15,000 13,000

パラメータ:

パラメータ 説明 値
\(n\) 製品の種類数 4
\(p_i\) 製品 \(i\) の仕入れ価格 \(p_1 = 1000, p_2 = 6000, p_3 = 11000, p_4 = 15000\)
\(v_i\) 製品 \(i\) の利益 \(v_1 = 700, v_2 = 5000, v_3 = 9000, v_4 = 13000\)
\(W\) 予算 100000 円

決定変数:

変数 説明 範囲
\(x_i\) 製品 \(i\) の仕入れ数量（非負の整数） \(0 \leq x_i \leq \lfloor W/p_i \rfloor, x_i \in \mathbb{Z}\)

目的関数:

目的 式
利益最大化 \(\text{maximize} \sum_{i=1}^n v_i x_i\)

限られた予算内で総利益を最大化するような仕入れ計画を立てます。

制約条件:

制約 式 説明
予算制約 \(\sum_{i=1}^n p_i x_i \leq W\) 仕入れる製品の合計金額が予算以下でなければならない
非負制約 \(x_i \geq 0, \forall i\) 仕入れ数量は非負でなければならない
整数制約 \(x_i \in \mathbb{Z}, \forall i\) 仕入れ数量は整数でなければならない

コード例:

以下に、この最適化問題を NAG Library for Python の MILP 問題専用ソルバー関数 handle_solve_milp を用いて解くコード例を示します。

import numpy as np
from naginterfaces.library import opt
from naginterfaces.library import mip

# 製品の種類とそれぞれの仕入れ価格および利益
products = ['A', 'B', 'C', 'D']
purchase_prices = np.array([1000, 6000, 11000, 15000])  # 仕入れ価格（円）
profits = np.array([700, 5000, 9000, 13000])  # 利益（円）
total_budget = 100000  # 予算（円）

# 製品の種類数
product_count = len(products)

# 最適化問題のハンドルを初期化
handle = opt.handle_init(product_count)

# 目的関数を設定（利益の最大化）
opt.handle_set_linobj(handle, profits)

# 各製品について予算内で購入可能な最大数を計算
max_purchases = np.floor(total_budget / purchase_prices)

# 変数の上限と下限を設定
opt.handle_set_simplebounds(handle, np.zeros(product_count), max_purchases)

# 予算制約を設定
opt.handle_set_linconstr(handle, bl=[-np.inf], bu=[total_budget], irowb=[1] * product_count, icolb=list(range(1, product_count + 1)), b=purchase_prices, idlc=0)

# 各変数を整数として扱う設定
variable_indices = np.array([i for i in range(1, product_count + 1)])
opt.handle_set_property(handle, 'INTEGER', variable_indices)

# ソルバーのオプションを指定
opt.handle_opt_set(handle, 'Print Level = 0')
opt.handle_opt_set(handle, 'Task = Maximize')  # 最大化問題

# 問題を解く
purchase_quantities = mip.handle_solve_milp(handle)[0]

# 結果の表示
total_spent = np.dot(purchase_prices, purchase_quantities)  # 総仕入れ額を計算
print("最適解と仕入れ額:")
for i in range(product_count):
    print(f"{products[i]}: {round(purchase_quantities[i])} 個, 仕入れ額: {round(purchase_quantities[i] * purchase_prices[i])}円")
print(f"---------------------------\n総利益: {round(np.dot(profits, purchase_quantities))} 円")
print(f"総仕入れ額: {int(total_spent)} 円")

# ハンドルの削除
opt.handle_free(handle)

コードの詳細説明

この例題では、NAG Library for Pythonの混合整数計画ソルバー handle_solve_milp.html を用いて予算制約下の仕入れ最適化問題を解いています。以下では、コードの各部分について詳しく説明します。

目的関数の設定

この問題の目的関数は、総利益の最大化です。総利益は、各製品の利益と仕入れ数量の積の和で表されます。数式で表すと、以下のようになります。

maximize 700 * x1 + 5000 * x2 + 9000 * x3 + 13000 * x4

ここで、x1, x2, x3, x4 は、それぞれ製品 A, B, C, D の仕入れ数量を表す決定変数です。また、700, 5000, 9000, 13000 は、それぞれ製品 A, B, C, D の利益を表す係数です。この式は線形結合であり、係数だけを渡すことで目的関数を完全に定義できます。

具体的には、目的関数を設定するために、その係数をopt.handle_set_linobj 関数の第2引数に渡します。

opt.handle_set_linobj(handle, profits)

変数の上下限の設定

次に、変数の上下限を設定します。この部分では、変数の下限 bl を、非負制約を満たすためにすべて0に設定しています。また、変数の上限 bu の各要素は、予算内で購入可能なそれぞれの製品の最大数に設定します。

max_purchases = np.floor(total_budget / purchase_prices)
opt.handle_set_simplebounds(handle, np.zeros(product_count), max_purchases)

opt.handle_set_simplebounds 関数は、決定変数の上下限を設定するための関数で、第1引数は問題のハンドル、第2引数は変数の下限を表す配列、第3引数は変数の上限を表す配列です。これらの配列の長さは、handle_init 関数で指定した決定変数の数（nvar）と一致している必要があります。

予算制約の設定

この問題には予算制約があります。予算制約は以下のような線形制約式として表現できます。

[1000x_1 + 6000x_2 + 11000x_3 + 15000x_4 ]

ここで、(x_1, x_2, x_3, x_4) は、それぞれ製品 A, B, C, D の仕入れ数量を表す決定変数で、係数 1000, 6000, 11000, 15000 は、それぞれ製品 A, B, C, D の仕入れ価格です。そして右辺の 100000 が総予算です。

この線形制約を opt.handle_set_linconstr 関数を用いて以下のように指定しています。

opt.handle_set_linconstr(handle, bl=[-np.inf], bu=[total_budget], irowb=[1] * product_count, icolb=list(range(1, product_count + 1)), b=purchase_prices, idlc=0)

関数の引数 irowb、icolb、b は、制約式の係数と決定変数の対応関係を表現します。

  • irowb：各係数が属する制約式のインデックス（1から始まる）を指定します。今回は制約式が1つだけなので、全ての要素が1になります。
    • irowb = [1, 1, 1, 1]
  • icolb：各係数が対応する決定変数のインデックス（1から始まる）を指定します。今回は (x_1), (x_2), (x_3), (x_4) に対応するインデックスを指定します。
    • icolb = [1, 2, 3, 4]
  • b：制約式の係数の値を指定します。今回は (x_1), (x_2), (x_3), (x_4) の係数である仕入れ価格を指定します。
    • b = [1000, 6000, 11000, 15000]

制約式の右辺は、bl と bu の引数で指定されます。bl は下限、bu は上限を示し、これによって不等号（<=, >=）や等号（=）を設定できます。下限が負の無限大（-∞）の場合、上限だけが適用されて「<=」の制約となり、上限が正の無限大（∞）の場合、下限だけが適用されて「>=」の制約となります。両方が同じ値の場合、「=」の制約となります。

今回の問題では線形制約が１つですが、線形制約が複数ある場合の例を以下に示します。

例1： [1000x_1 + 6000x_2 ] [11000x_3 + 15000x_4 ]

この場合、irowb、icolb、b、bl、bu は以下のように設定します。

irowb = [1, 1, 2, 2]  # 第1制約式に属する係数のインデックスを1、第2制約式に属する係数のインデックスを2とする
icolb = [1, 2, 3, 4]  # 各係数が対応する変数のインデックス
b = [1000, 6000, 11000, 15000]  # 各係数の値
bl = [50000, 80000]  # 各制約式の下限値
bu = [np.inf, np.inf]  # 各制約式の上限値（ここでは両方とも正の無限大）

例2： [1000x_1 + 6000x_2 ] [15000x_2 + 11000x_3 + 8000x_4 = 30000]

この場合、irowb、icolb、b、bl、bu は以下のように設定します。

irowb = [1, 1, 2, 2, 2]  # 第1制約式に属する係数のインデックスを1、第2制約式に属する係数のインデックスを2とする
icolb = [1, 2, 2, 3, 4]  # 各係数が対応する変数のインデックス
b = [1000, 6000, 15000, 11000, 8000]  # 各係数の値
bl = [-np.inf, 30000]  # 各制約式の下限値（第1制約式の下限値は負の無限大、第2制約式の下限値は30000）
bu = [50000, 30000]  # 各制約式の上限値（第1制約式の上限値は50000、第2制約式の上限値は30000）

変数の整数性の設定

ここでは、opt.handle_set_property を使用して、特定の変数が整数値をとるように設定します。これは、製品の仕入れ数量が整数でなければならないためです。

variable_indices = np.array([i for i in range(1, product_count + 1)])
opt.handle_set_property(handle, 'INTEGER', variable_indices)

variable_indices は、1から始まる変数のインデックス配列です。例えば、product_count が4の場合、variable_indices は [1, 2, 3, 4] となりますが、これらのインデックスに対応する変数が整数として扱われるようになります。

ソルバーのオプションの指定

ここでは、opt.handle_opt_set 関数を使用してソルバーのオプションを設定します。

opt.handle_opt_set(handle, 'Print Level = 0')
opt.handle_opt_set(handle, 'Task = Maximize')  # 最大化問題
  • 'Print Level = 0' は、出力を最小限に抑えるオプションです。
  • 'Task = Maximize' は、今回の問題が最大化問題であることを示します。デフォルトではソルバーは最小化問題として設定されているため、最大化問題にするためにこのオプションを指定します。

問題の解法と結果の取得

mip.handle_solve_milp 関数を使用して、MILP 問題を解きます。この関数は最適解を含む配列や追加情報を返します。ここでは、最適解のみを取得しています。

purchase_quantities = mip.handle_solve_milp(handle)[0]

結果の表示

得られた最適解を使用して、各製品の仕入れ数量と仕入れ額、総利益、総仕入れ額を計算し、結果を表示します。

total_spent = np.dot(purchase_prices, purchase_quantities)  # 総仕入れ額を計算
print("最適解と仕入れ額:")
for i in range(product_count):
    print(f"{products[i]}: {round(purchase_quantities[i])} 個, 仕入れ額: {round(purchase_quantities[i] * purchase_prices[i])}円")
print(f"---------------------------\n総利益: {round(np.dot(profits, purchase_quantities))} 円")
print(f"総仕入れ額: {int(total_spent)} 円")

ハンドルの削除

opt.handle_free(handle)

opt.handle_free 関数を使用して、最適化問題のハンドルを削除しています。これにより、ハンドルに関連付けられたメモリが解放されます。

まとめ

このセクションでは、以下の主要な点を学びました：

  1. NAG Library の混合整数計画問題を解く関数 handle_solve_milp の使い方
  2. 目的関数の設定方法
  3. 変数の非負制約を含む上下限の設定方法
  4. 線形制約条件の設定方法（スパース行列形式を使用し、インデックスで変数を指定）
  5. 変数の整数制約の指定方法 (opt.handle_set_property)
  6. ソルバーのオプション設定方法 (opt.handle_opt_set)
  7. 結果の取得方法

関連資料

NAGライブラリ最適化機能の詳細

  • URL: https://www.nag-j.co.jp/naglib/optimization/index.htm
  • 説明: NAGライブラリが提供する数多くの最適化関数についての情報がまとめられています。本ガイドで取り上げた関数以外にも、線形計画法、非線形計画法、二次計画法など、様々な最適化問題に対応する手法が網羅されています。

分野別の数理最適化サンプルコード

  • URL: https://www.nag-j.co.jp/naglib/optimization/python/opt-examples/
  • 説明: NAG Library for Pythonを用いた数理最適化のサンプルコード集。本ガイドで取り上げた例以外にも、様々な最適化問題で実際に使われそうな問題の具体的な実装例を通して理解を深めることができます。

問題の種類別サンプルコード

NAGライブラリの開発環境

  • URL: https://www.nag-j.co.jp/naglib/environment.htm
  • 説明: NAGライブラリはPython以外の言語でも利用可能です。さまざまなプログラミング環境や開発環境での利用方法について情報が提供されています。

NAG Library for Python マニュアル

NAG Library for Python トップページ


NAG数理最適化ソルバーを試す
関連情報