※更新履歴

2022/08/15 Githubのソースコード修正、それに伴う記事内のサンプルコード修正
- @nariaki3551 さんに最適化結果をpulpにデコードする関数を実装していただきました。簡潔でとても便利です。ありがとうございます！
- @nariaki3551 さんにReadMeのサンプルコードのバグを修正していただきました。ありがとうございます！
2022/8/17 ソースコード修正(v0.1.2)、改善点の更新
- pypiにアップロードしたので、$ pip install pulp2matでインストールできるようになりました。
- 最大化問題に対応、ただし目的関数の値の正負が反転します。
2022/8/21 ソースコード修正(v0.1.3)、それに伴う記事内のサンプルコード修正
- @nariaki3551 さんにpull requestいただきました。以前は変数を辞書形式で定義する必要がありましたが、pulp.LpProblemのオブジェクトからとってこれるようになり、使いやすくなりました。

参考にしたもの

背景

SciPyのMILPは使えるのか？？（SciPy==1.9.0の新機能の検証）によると、ScipyのMILPソルバーが高速だが、使いにくいとのこと。確かに変数の種類が複数あると、行列表現するのは大変そう。だが、

実務MILPを解く場合は、GurobiやCPLEXを使うしかない。

こう言われてしまうと実務家として頑張りたくなってしまったので、pulpの定式をscipyの入力に変換する関数を自作することにした。引用元ではpyomoを使っているが、pyomoは使ったことがないので今回はpulpを用いた。

実装したもの

Github:

Pypi:

使い方

$ pip install pulp2matでインストールできる。

~~ReadMeに書いてあるが、使い方の注意点として、すべての変数を辞書形式で定義する必要がある。~~
※ 2022/8/21 修正: v0.1.3から、変数の持ち方は特に問わないような実装に変更した。

例えば1次元ビンパッキング問題の場合、以下$x$, $y$を定義する。

$x_{ij} \in 0,1$: アイテム$i$をビン$j$に詰め込むか否か

$y_{j} \in 0,1$: ビン$j$にアイテムが入っているか否か

pulpでの実装はこのようになる。

import pulp as pl
import numpy as np

item_sizes = np.array([7, 3, 3, 1, 6, 8, 4, 9, 5, 2])
num_items = len(item_sizes)
num_bins = len(item_sizes)
bin_size = 10

# Variables * must be defined as dictionaries
x = {
    (i, j): pl.LpVariable("x_{}_{}".format(i, j), cat=pl.LpBinary)
    for i in range(num_items)
    for j in range(num_bins)
}
y = {
    j: pl.LpVariable("y_{}".format(j), cat=pl.LpBinary)
    for j in range(num_bins)
}

problem = pl.LpProblem()

# Bin size constraint for each bin
for j in range(num_bins):
    problem += (
        pl.lpSum(
            x[i, j] * item_sizes[i] for i in range(num_items)
        )
        <= bin_size * y[j]
    )
# One-hot constraint for each item
for i in range(num_items):
    problem += pl.lpSum(x[i, j] for j in range(num_bins)) == 1

# Objective: minimize number of bins used.
problem += pl.lpSum(y[j] for j in range(num_bins))

pulp.LpProblem オブジェクトをscipy.optimize.milpへの入力に変換するために、pulp2mat.convert_all()関数を用いる。

import pulp2mat
from scipy.optimize import milp
c, integrality, constraints, bounds = pulp2mat.convert_all(problem)
result = milp(c, integrality=integrality, constraints=constraints, bounds=bounds)

改善点

~~まだ最大化問題には対応していない。最大化問題の場合は、目的関数の正負を反転させる必要がある。~~ 2022/8/17対応済み
問題サイズが大きい場合や反復実行する場合は、変換の速度が問題になりそう。
~~pip installできるようにする~~ 2022/8/17対応済み
まだビンパッキング問題でしかテストしていないので、バグなどのご指摘を歓迎します。