はじめに
この記事の目的：
- SciPy==1.9.0 で導入された
scipy.optimize.milpの簡単な説明 scipy.optimize.milpの使い方の解説
この記事で扱わないこと
scipy.optimize.milpで採用されている最適化手法- 他のソルバーとの比較
[追記]
SciPy==1.9.0 のリリース
Atsushi さんによる以下の記事を参照してください．
先行記事
書籍 [2] で挙げられている問題を scipy.optimize.milp を用いて解く方法について解説されています．
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
まず，簡単な問題として Mathworks の intlinprog のドキュメントで扱われている以下のような混合整数計画問題を考えます：
\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 は整数制約を表します．他に 2 と 3 を設定することもできますが，それらの説明は 公式ドキュメント を参照してください．また，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ではじめる数理最適化: ケーススタディでモデリングのスキルを身につけよう』
