前回の記事で，SciPy==1.9.0 で新たに導入された scipy.optimize.milp に関する簡単な紹介をしました．
本記事では，上記の記事で扱った問題と同様の問題を Pyhton-MIP で求解する方法を紹介します．Python-MIP のより詳細な使用方法については公式ドキュメント [1] を参照してください．
解くべき混合整数計画問題（再掲）
次の混合整数計画問題を考えます：
\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}$ とします．
scipy.optimize.milp による求解（再掲）
必要なパッケージを import しておきます．
import numpy as np from scipy.optimize import milp from scipy.optimize import LinearConstraint from scipy.optimize import Bounds
解くべき問題の係数行列などを定義します．
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
制約条件に関わる部分を定義します．
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
Python-MIP の CBC による求解
必要なパッケージを import しておきます．
import numpy as np import mip
同様の問題を扱うため，係数行列などの定義は同様です．
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
次に，決定変数，目的関数，制約条件をそれぞれ定義して求解を行います．
model = mip.Model() # define decision variable x = {} for i in range(n): x[i] = model.add_var(lb=0, ub=mip.INF, var_type=mip.INTEGER) # define objective function model.objective = mip.minimize(mip.xsum(c[i] * x[i] for i in range(n))) # define constraints for i in range(m): model.add_constr( mip.xsum(A[i,:][j] * x[j] for j in range(n)) <= b[i] ) # solve model.optimize()
model.optimize() を実行すると以下のような出力が得られます．実際には，最適化計算の各反復における情報が出力されますが，ここでは省略します．
Cbc0032I Strong branching done 29094 times (142226 iterations), fathomed 305 nodes and fixed 796 variables Cbc0041I Maximum depth 33, 3997 variables fixed on reduced cost (complete fathoming 1136 times, 157339 nodes taking 667811 iterations) Total time (CPU seconds): 8.99 (Wallclock seconds): 9.63 <OptimizationStatus.OPTIMAL: 0>
得られた最適値と最適解を確認してみましょう．
print('Objective value: {model.objective_value:.15}'.format(**locals())) print('Solution: ', end='') for v in model.vars: if v.x > 1e-5: print('{v.name} = {v.x}'.format(**locals())) print(' ', end='')
実行すると以下のような結果が得られます．
Objective value: 7.44184450406387
Solution: var(0) = 1.0
var(2) = 1.0
var(5) = 1.0
var(6) = 1.0
var(10) = 1.0
var(13) = 4.0
var(20) = 2.0
var(21) = 1.0
var(26) = 3.0
var(27) = 1.0
var(28) = 1.0
var(30) = 1.0
var(31) = 3.0
var(34) = 1.0
var(37) = 1.0
最適値，最適解ともに scipy.optimize.milp で計算した結果と整合していることが分かります．なお，求解にかかった時間は CPU 時間だけで比較すると，milp は 17.8 秒，Python-MIP (cbc) は 8.99 秒でした．
