冷めたコーヒー

Weniger, aber besser

SciPy==1.9.0 で使用可能になった混合整数計画問題のソルバー scipy.optimize.milp に関する備忘録

はじめに

この記事の目的:

  • SciPy==1.9.0 で導入された scipy.optimize.milp の簡単な説明
  • scipy.optimize.milp の使い方の解説

この記事で扱わないこと

  • scipy.optimize.milp で採用されている最適化手法
  • 他のソルバーとの比較

[追記]

  • 公式ドキュメントにナップザック問題を解くチュートリアルが掲載されたので参考 [3] にリンクを追記しました (2022/09/08)

SciPy==1.9.0 のリリース

Atsushi さんによる以下の記事を参照してください.

myenigma.hatenablog.com

先行記事

書籍 [2] で挙げられている問題を scipy.optimize.milp を用いて解く方法について解説されています.

github.com

scipy.optimize.milp で解ける最適化問題

scipy.optimize.milp では以下の形式の混合整数計画問題を解くことができる:

\begin{align} \min_{\mathbf{x} \in \mathbb{R}^n} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{b}_l \leqslant\mathbf{A}\mathbf{x} \leqslant \mathbf{b}_u,\\ &\mathbf{l} \leqslant \mathbf{x} \leqslant \mathbf{u},\\ &x_i \in \mathbb{Z},\, i \in \mathcal{I} \subseteq [n]. \end{align}

ここで,$\mathcal{I}$ は添字集合 $[n] := \{1,2,\dots,n\}$ の部分集合 (e.g., $\mathcal{I} = \{2, 3, 5 \}$) を表す.

残念ながら,scipy.optimize.milp を用いて上記の形式でない混合整数計画問題を解くことは今のところできません.頑張って上記の形式に書き直すか,他の最適化モデリング言語などを用いましょう.

scipy.optimize.milp の使い方

導入

以下で scipy==1.9.0 のインストールをしましょう:

python -m pip install scipy==1.9.0

以下では,scipy.optimize に含まれるパッケージを使用します.ついでに numpy も使用するので import しておきましょう.

import numpy as np
from scipy.optimize import milp
from scipy.optimize import LinearConstraint
from scipy.optimize import Bounds

例 1

まず,簡単な問題として Mathworksintlinprog のドキュメントで扱われている以下のような混合整数計画問題を考えます:

\begin{align} \min_{x_1,x_2} \quad &8x_1 + x_2 \\ \text{subject to} \quad&x_1 + 2x_2 \geqslant -14,\\ &-4x_1 - x_2 \leqslant -33,\\ &2x_1 + x_2 \leqslant 20,\\ &x_2 \in \mathbb{Z}. \end{align}

この問題はベクトルと行列を用いて次のように表すことができます:

\begin{align} \min_{x_1,x_2} \quad &\begin{bmatrix}8 \\ 1\end{bmatrix}^\top \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \\ \text{subject to} \quad&\begin{bmatrix}-1 & -2 \\ -4 & -1 \\2 & 1\end{bmatrix}\begin{bmatrix}x_1 \\ x_2\end{bmatrix} \leqslant \begin{bmatrix}14 \\ -33 \\ 20\end{bmatrix},\\ &x_2 \in \mathbb{Z}. \end{align}

この問題を scipy.minimize.milp を用いて解けるように係数ベクトルや行列を以下のように定めます.なお,この問題は $(x_1,x_2)^\top$ に関する上下限制約はないので特別指定する必要はありません.また,$\mathbf{Ax}$ に関する下限制約もないので,以下のように b_l の各要素に -np.inf を設定します.整数制約は $x_2$ に対してのみ課されるため,以下のように integrality = np.array([0, 1]) と設定します.なお,0 は連続変数,1 は整数制約を表します.他に 23 を設定することもできますが,それらの説明は 公式ドキュメント を参照してください.また,integrality を指定しなかった場合,デフォルトでは 0 が指定されます.

c = np.array([8, 1])
A = np.array([[-1, -2], [-4, -1], [2, 1]])
b_u = np.array([14, -33, 20])
b_l = np.full_like(b_u, -np.inf)
constraints = LinearConstraint(A, b_l, b_u) $b_l <= Ax <= b_u
integrality = np.array([0, 1]) # x_2 is integer

求解は以下のように行います.

%%time
res = milp(
    c=c,
    constraints=constraints,
    integrality=integrality,
)
print(res)

結果は以下のようになります.

            fun: 59.0
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 59.0
        mip_gap: 0.0
 mip_node_count: 1
         status: 0
        success: True
              x: array([6.5, 7. ])
CPU times: user 5.49 ms, sys: 0 ns, total: 5.49 ms
Wall time: 3.24 ms

Mathworksドキュメント に掲載されている結果と同じ結果が得られていることを確認できます.

例 2

次の混合整数計画問題を考えます:

\begin{align} \min_{\mathbf{x}} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{A}\mathbf{x} \leqslant \mathbf{b},\\ &\mathbf{x} \in \mathbb{Z}^n,\,\mathbf{x} \geqslant \mathbf{0}. \end{align}

ここで,$\mathbf{c}\in\mathbb{R}^n, \, \mathbf{A}\in\mathbb{R}^{m\times n}, \, \mathbf{b}\in\mathbb{R}^{m}$ とします.この問題の実行可能領域 $\mathcal{F}=\left\{ \mathbf{x} \in \mathbb{Z}^n \mid \mathbf{A}\mathbf{x} \leqslant \mathbf{b},\, \mathbf{x} \geqslant \mathbf{0} \right\}$ に実行可能解が存在することを保証するために,$\overline{\mathbf{x}} \in \mathbf{Z}^n$ をランダムに生成し,その $\overline{\mathbf{x}} \in \mathbf{Z}^n$ が実行可能領域 $\mathcal{F}$ に含まれるように $\mathbf{A}$ と $\mathbf{b}$ を設定します.

np.random.seed(42)
n, m = 40, 80
c = np.random.rand(n)
x = np.random.randint(0, 5, n)
A = np.random.randn(m, n)
s = 10 * np.ones_like(m)
b = A @ x + s

決定変数 $\mathbf{x}$ に対する下限制約が存在するため,scipy.optimize.Bounds を使用します.決定変数 $\mathbf{x}$ に対する上限制約は存在しないので,以下のように u の各要素に np.inf を設定します.また,今回の問題では決定変数のすべての要素は整数であることを課しているため,integrality のすべての要素を 1 とします.

l = np.zeros_like(c)
u = np.full_like(c, np.inf)
b_u = b
b_l = np.full_like(b_u, -np.inf)

bounds = Bounds(l, u)
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.ones_like(c)

求解します.

%%time
res = milp(
    c=c,
    bounds=bounds,
    constraints=constraints,
    integrality=integrality,
)
print(res)

結果は以下の通りです.

            fun: 7.441844504063867
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 7.44111388697805
        mip_gap: 9.8176881473139e-05
 mip_node_count: 49538
         status: 0
        success: True
              x: array([1., 0., 1., 0., 0., 1., 1., 0., 0., 0., 1., 0., 0., 4., 0., 0., 0.,
       0., 0., 0., 2., 1., 0., 0., 0., 0., 3., 1., 1., 0., 1., 3., 0., 0.,
       1., 0., 0., 1., 0., 0.])
CPU times: user 17.8 s, sys: 665 ms, total: 18.4 s
Wall time: 17.6 s

決定変数の次元 $n=40$ の問題を解いてみました.およそ 20 秒弱の時間で最適解を計算することができました.

おわりに

この記事では SciPy==1.9.0 で新たに導入された scipy.optimize.milp を用いて混合整数計画問題を解く方法について解説をしました.

参考

[1] scipy.optimize.milp
[2] 『Pythonではじめる数理最適化: ケーススタディモデリングのスキルを身につけよう』

[3] Knapsack problem example