Pythonから混合整数計画ソルバー(MIP solver)を使用

はじめに

数理最適化の問題を解くときに第一に挙がることが多いのが線形/整数計画法だというイメージはあるが、 いつもheuristic的な解法に頼りきりで、使えると便利だろうなと思う場面がありながら ほとんど混合整数計画法には手を出していなかったので一度は触った方がよいと思い試してみました。 以前C++でシンプレックスで分枝限定法でソルバーを作ってみたことはあるが、性能的には悲惨な状態で ツールを使った方がはやいと感じたりしたので放置状態でした。 本格的に使用したいわけではなく考察の補助として使用できればと考えています。

PythonのMIP Solver

混合整数計画法おそらく一番に挙がるのがgurobiだと思うが、ライセンスが必要で若干面倒なので 手軽に使えるCOIN CBCを使用した。さらに性能が必要だと感じたらまた考える。 COIN CBCのPythonインターフェースライブラリとしては

• PuPL
• Python-MIP

があるようなのでこれらを使用してみて 良い方を選ぶことにしました。

数独(Sudoku)ソルバーを作る

まず練習で数独ソルバーを作りました。 数独は9*9の81マスに数字を埋めるパズルですが、1マスにつき1-9の9つのバイナリ変数を設定します。 よって変数の数は9*9*9個になります。例えば、(i,j)=(2,4)のマスに5が入る場合は$X_{5,2,4}=1$となります。 ルールでは行、列、3*3の正方形ブロックにはそれぞれ1-9の数字が一つずつ入るので、

$\sum _{i=1}^9{x}_{i,j,k}=1$

$\sum _{j=1}^9{x}_{i,j,k}=1$

$\sum _{i=1}^3\sum _{j=1}^3{x}_{i,j,k}=1$

のように制約が定式化されます。三つ目の式はブロックの制約で行・列それぞれに3つ合計9つ分が必要になります。 まずメイン部分は以下になります。 PuPLのドキュメントのものとほぼ同じです。

    
#PuPL
import pulp as pl
# Python-MIP
# import mip
import time
from functools import wraps
 
def timer(func):
    @wraps(func)
    def wrapper(*args, **kargs):
        start = time.time()
        result = func(*args, **kargs)
        elapsed_time = time.time() - start
        print(f"elapsed {elapsed_time * 1000:.2f} ms")
        return result
    return wrapper
 
def main():
    with open("sudoku_data.txt") as f:
        data = [list(map(int, s.split())) for s in f.readlines()]
 
    choices = solve_sudoku(data)
 
    for r in range(1, 10):
        for c in range(1, 10):
            for v in range(1, 10):
                if pl.value(choices[v][r][c]) == 1:
                #if choices[v][r][c].x == 1:
                    print(v, end=" ")
        print()
Copy

PulPのソルバー部分は以下になります。

    
@timer
def solve_sudoku(data):
    VALS = ROWS = COLS = range(1, 10)
    Boxes = [
        [(3 * i + k + 1, 3 * j + l + 1) for k in range(3) for l in range(3)]
        for i in range(3)
        for j in range(3)
    ]
 
    prob = pl.LpProblem("Sudoku_Problem")
 
    choices = pl.LpVariable.dicts("Choice", (VALS, ROWS, COLS), cat="Binary")
 
    for (v, r, c) in data:
        prob += choices[v][r][c] == 1
 
    for r in ROWS:
        for c in COLS:
            prob += pl.lpSum([choices[v][r][c] for v in VALS]) == 1
    
    for v in VALS:
        for r in ROWS:
            prob += pl.lpSum([choices[v][r][c] for c in COLS]) == 1
        for c in COLS:
            prob += pl.lpSum([choices[v][r][c] for r in ROWS]) == 1
        for b in Boxes:
            prob += pl.lpSum([choices[v][r][c] for (r, c) in b]) == 1
 
    prob.solve(pl.PULP_CBC_CMD(msg=0))
 
    print("Status:", pl.LpStatus[prob.status])
 
    return choices
Copy

Python-MIPのソルバー部分は以下になります。 PuPLとほとんど同じように書くことができます。

    
@timer
def solve_sudoku(data):
    VALS = ROWS = COLS = range(9)
    Boxes = [
        [(3 * i + k, 3 * j + l) for k in range(3) for l in range(3)]
        for i in range(3)
        for j in range(3)
    ]
 
    model = mip.Model()
 
    x = [[
        [model.add_var(var_type=mip.BINARY) for j in COLS]
        for i in ROWS] for k in VALS]
    
    for (v, r, c) in data:
        model += x[v-1][r-1][c-1] == 1
 
    for r in ROWS:
        for c in COLS:
            model += mip.xsum(x[v][r][c] for v in VALS) == 1
    
    for v in VALS:
        for r in ROWS:
            model += mip.xsum(x[v][r][c] for c in COLS) == 1
        for c in COLS:
            model += mip.xsum(x[v][r][c] for r in ROWS) == 1
        for b in Boxes:
            model += mip.xsum(x[v][r][c] for (r, c) in b) == 1
 
    model.threads = -1
    model.verbose = 0
 
    model.optimize()
 
    print("Status:", model.status)
 
    return x
Copy

実際に実行したところ、 https://sudoku.com/ のEvilレベルを100ms以下で解くことができました。 普通の枝狩りDFSでもこれほどは速くは解けないのでかなり速いです。 この程度ならPuPL、Python-MIPどちらでも速度や使用感の差はありません。

TSPソルバー

続いて二つのソルバーの比較と使用感を確認するためTSPソルバーを作りました。 こちらのソースコードを書いている途中、PuLPの付属CBCバイナリにはマルチスレッドが対応していないことに気づきました。 coinbrewという独自のビルダーでビルドをしなくてはならないようなので、こちらのマルチスレッド化はあきらめました。 Python-MIPの方は確認したところマルチスレッドにも対応していました。 無向グラフの場合の定式化は以下になります。 まず目的関数は以下を最小化することです。dijは頂点間距離、xijは辺ijを使用するかどうかのバイナリ変数です。

$\sum _{i,j}{d}_{ij}{x}_{ij}$

各頂点の入次数と出次数は1なので

$\sum _{j\in V/\{i\}}{x}_{i,j}=1$

$\sum _{i\in V/\{j\}}{x}_{i,j}=1$

また複数のループが発生しないように頂点0(始点/終点)以外の頂点にポテンシャルを導入して たどる順番にポテンシャルが増加するようにします。

$v_i-(n+1)x_{ij}\ge v_j-n$

PuLPの場合ソースコードは以下になります。

    
@timer
def solve_tsp(dists):
    n = len(dists)
 
    prob = pl.LpProblem("TSP", pl.LpMinimize)
    
    x = pl.LpVariable.dicts("vertex",(range(n), range(n)), cat="Binary")
    y = pl.LpVariable.dicts("potential", range(n), 0, None)
 
    prob += pl.lpSum(dists[i][j] * x[i][j] for i in range(n) for j in range(n))
 
    for i in range(n):
        prob += pl.lpSum(x[i][j] for j in range(n) if i!=j) == 1
 
    for i in range(n):
        prob += pl.lpSum(x[j][i] for j in range(n) if i!=j) == 1
 
    for i in range(n):
        for j in range(n):
            if i!=0 and j!=0 and i!=j:
                prob += y[i] - (n+1)*x[i][j] >= y[j]-n
 
    prob.solve(pl.PULP_CBC_CMD(msg=0))
 
    print("Status:", pl.LpStatus[prob.status])
    print("objective", pl.value(prob.objective))
 
    return x
Copy

Python-MIPの場合は以下になります。 model.threadsを-1としてmultithread化しています。

    
@timer
def solve_tsp(dists):
    n = len(dists)
 
    model = mip.Model()
 
    x = [
        [model.add_var(var_type=mip.BINARY) for j in range(n)]
        for i in range(n)]
    y = [model.add_var() for i in range(n)]
 
    model.objective = mip.minimize(mip.xsum(dists[i][j]*x[i][j] for i in range(n) for j in range(n)))
 
    for i in range(n):
        model += mip.xsum(x[i][j] for j in range(n) if i!=j) == 1
 
    for i in range(n):
        model += mip.xsum(x[j][i] for j in range(n) if i!=j) == 1
 
    for i in range(n):
        for j in range(n):
            if i!=0 and j!=0 and i!=j:
                model += y[i] - (n+1)*x[i][j] >= y[j]-n
 
    model.threads = -1
    model.verbose = 0
 
    model.optimize()
 
    print("Status:", model.status)
    print("objective", model.objective_value)
 
    return x
Copy

最初は頂点100辺数100*99で試しましたが短時間では完了しそうになかったので、頂点数を30に減らして行いました。 頂点の生成にはまずx,y座標をrandom整数で0-300の範囲で作成し、それぞれの頂点間距離をユークリッド距離で算出し小数点以下を切り捨てて整数にして辺のコストとしました。 以下は上記ソースコードで得られたTSPの解を 自作で作ったRust+wasmのTSPソルバー兼ビジュアライザーによって図示したものです。

インターフェースライブラリ比較

以下に30頂点870辺TSPの実行時間による比較を示します。 シングルスレッドで差が出ているのはcbcバイナリのバージョンの差と思われます。 Python-MIPはマルチスレッドの効果を確認できました。

TSPをMIPソルバーで解くのは向いていないようでしたが、今回はお試しと比較が目的なので問題はありません。 変数の数と式の数は辺数オーダーになるので仕方ありません。 どちらかというと近似解を求めてからMIPソルバーにかけるという方法の方が良いように思います。

その他機能等について

双方のライブラリにもwarm start、time limit、下界(上界)とのgap差条件による停止、 モデルの保存、特定の変数の固定機能は使えるようでどちらもほぼ同じ感じで使用可能です。 ただofficialのドキュメントはPython-MIPの方が詳しく記載されていました。 マルチスレッドの機能の面も考慮するとPuPLよりはPython-MIPの方が使用しやすそうです。 ただ1.13.0よりバージョンが更新されていないことだけは気になります。 しばらくはPython-MIPをメインに使用することにします。

Reference