Help us understand the problem. What is going on with this article?

ついに使い物になるフリーの数理最適化ソルバーが登場? MIPCL

はじめに

数理最適化問題(数理計画問題)、特に整数最適化問題(整数計画問題)の式を解いてくれる無償で商用利用可能なソルバーとしては、PuLP(デフォルトのソルバーとしてCOIN-CBCが付属)が比較的有名だと思います(使用例)。また、Google OR-Tools(こちらもデフォルトのソルバーとしてCOIN-CBCが付属)も知られています。ただし、これらのソルバーは、有償のソルバー(Gurobi, CPLEX, XPRESSなど)と比べると、計算速度が大きく劣るため、実務で解きたい問題が許容時間内に解けない場合があります。
ここでは、上記の無償ソルバーよりも速く問題を解いてくれそうなMIPCL」という新しめのソルバーを紹介します。

MIPCLとは

MIPCL
公式サイト|http://www.mipcl-cpp.appspot.com/index.html
2年くらい前でしょうか、ソルバーを比較するベンチマークサイトに突如としてMIPCLというソルバーが出現し、フリーの中ではSCIPに次ぐ性能という評価になりました。ただしSCIPは商用利用したい場合は作者への相談が必要なので、そういう敷居がない中ではMIPCLは最高性能になりました。
MIPCLのよさげなところは、特に何の設定もなしにPythonインターフェースから呼んでも、マルチスレッドで動作してくれるところにあります。(こちらのブログエントリーによると、PuLPについても、同こんのソルバーCOIN-CBCを呼ぶ際にマルチスレッドで動かすことができるようなのですが、なぜか私の環境ではうまく動きません…)

MIPCLは、整数最適化のほかにも、凸二次最適化問題(凸二次計画問題)や二次錐最適化問題(二次錐計画問題)を扱えるそうです。

ライセンス形態ですが、公式ページには「GNU Lesser Public License (GLPL)」と書いていますが、おそらく誤植です。ダウンロードしたファイルには「GNU LESSER GENERAL PUBLIC LICENSE Version 3」と書いていまして、おそらくこちらのほうが正しい表記だと思います。いわゆるLGPLですね。なので、いにしえのGPLとは違い、MIPCLを外部から呼び出すという通常の使い方をする場合は、自分のコードを公開しなくてよさそうです(よね?)。

実験環境

  • Windows 10 64bit
  • Python 3.7.1(Anacondaをインストールしてアップデート)
  • MIPCL-PY 2.1.3
  • Visual Studio Code 1.33.1

インストール

  1. 公式サイト → Download → OSを選択 → Download(個人情報は入力しなくてもよい)
  2. setup.exe → VC14ランタイムのインストールが始まる → OSが勝手に再起動する
    • すでに(類似の)ランタイムがインストールされていれば、MIPCL同こんのものはインストールしなくてもよい模様
  3. MIPCL-PY.msi → インストール対象を選択 → インストールフォルダーを選択
    • インストール対象をJust MeAll Usersにするかは、おそらく環境変数PYTHONPATHがユーザー環境変数になるかシステム環境変数になるかに影響するのではないか
      • 今どきPYTHONPATH使うパッケージって…
      • Pythonなのに.msiのインストーラーって…
    • インストールフォルダーは、デフォルトではC:\Program Files\PNN\MIPCL-PY_64\になっていますが、もしなんかいやな予感がする方はC:\Users\(ユーザー名)\Documents\Software\PNN\MIPCL-PY_64\あたりに入れると精神衛生上よさそう

最適化例(1)

簡単な例を紹介します。
問題例はこちらのエントリーの例(1)と同じものを使用します。

コード

p1.py
# coding: utf-8

# Dummy Nihongo: あ(私の使っているテキストエディタにUTF-8と認識させるため)


# Import
# -----------------------------------------------------------------------------
# MIPCL-PY
import mipcl_py.mipshell.mipshell as mp
# -----------------------------------------------------------------------------


# Main
# -----------------------------------------------------------------------------
print(f'-' * 40)
print(f'[宣言]')

# 数理最適化問題を宣言(名前をつける必要あり)
prob = mp.Problem('Problem-1')

# 変数を宣言(連続変数)
# 引数を全部指定する場合
x = mp.Var(name='x', type=mp.REAL, lb=0.0, ub=mp.VAR_INF, priority=0)
# 引数をけっこう省略する場合(type以降のデフォルト値は↑と同じ)
y = mp.Var(name='y')

# 目的関数を宣言(最大化)
mp.maximize(x + y)

# 制約条件を宣言
2 * x + y <= 2
x + 2 * y <= 2


print(f'-' * 40)
print(f'[計算]')

# 計算(計算ログを出力)
prob.optimize(silent=False)

print(f'-' * 40)
print(f'[出力]')
# テキスト出力
prob.printSolution()

print(f'-' * 32)
# 目的関数値や変数の値を個別に取得
print(f'目的関数値 = {prob.getObjVal()}')
print(f'変数 x の値 = {x.val}')
print(f'変数 y の値 = {y.val}')

print(f'-' * 40)
# -----------------------------------------------------------------------------

実行結果

----------------------------------------
[宣言]
----------------------------------------
[計算]
Scaling...
Min exp = 0,  Max exp = 1

Preprocessing Time: 0.000
Solving...

Building initial basis...
     0.001          1          0      D/S                  0.0000
     0.003          3          0      P/S                  1.3333
Solution Time: 0.6
----------------------------------------
[出力]
Optimal objective value is 1.3333
'x' = 0.6667
'y' = 0.6667
--------------------------------
目的関数値 = 1.3333333333333335
変数 x の値 = 0.6666666666666667
変数 y の値 = 0.6666666666666666
----------------------------------------

最適化例(2)

$\sum$を使った整数最適化問題の例を紹介します。
問題例はこちらのエントリーのサンプル2と同じものを使用します。ただし、クラス周りの記述はやめて、シンプルなコードにします。

コード

p2.py
# coding: utf-8

# Dummy Nihongo: あ(私の使っているテキストエディタにUTF-8と認識させるため)


# Import
# -----------------------------------------------------------------------------
# MIPCL-PY
import mipcl_py.mipshell.mipshell as mp
# -----------------------------------------------------------------------------


# Main
# -----------------------------------------------------------------------------
print(f'-' * 40)
print(f'[集合の定義]')
I = ['Aさん', 'Bさん', 'Cさん']
J = ['東京', '大阪']
c = {
    ('Aさん', '東京'):  7, ('Aさん', '大阪'):  2, 
    ('Bさん', '東京'): 10, ('Bさん', '大阪'): 13, 
    ('Cさん', '東京'):  3, ('Cさん', '大阪'):  6
}

print(f'-' * 40)
print(f'[宣言]')

# 数理最適化問題を宣言(名前をつける必要あり)
prob = mp.Problem('Problem-2')

# 変数を宣言(0-1変数)
# x[社員名, 仕事名] を gp.Var クラスのオブジェクト を表す辞書として定義
x = {}
for i in I:
    for j in J:
        # 引数のうち、 priority を省略
        x[i, j] = mp.Var(name=f'x({i},{j})', type=mp.BIN)


# 目的関数を宣言(最小化)
mp.minimize(mp.sum_(c[i, j] * x[i, j] for i in I for j in J))

# 制約条件を宣言
for i in I:
    mp.sum_(x[i, j] for j in J) <= 1
for j in J:
    mp.sum_(x[i, j] for i in I) >= 1

print(f'-' * 40)
print(f'[計算]')

# 計算(計算ログを出力)
prob.optimize(silent=False)

print(f'-' * 40)
print(f'[出力]')
# テキスト出力
prob.printSolution()

print(f'-' * 32)
# 目的関数値や変数の値を個別に取得
print(f'目的関数値 = {prob.getObjVal()}')
for i in I:
        for j in J:
            if x[i, j].val > 0.98:
                print(f"{i}に{j}を割り当て")

print(f'-' * 40)
# -----------------------------------------------------------------------------

実行結果

----------------------------------------
[集合の定義]
----------------------------------------
[宣言]
----------------------------------------
[計算]
Start preprocessing: Rows - 5, Cols - 6 (Int - 6, Bin - 6), NZs - 12
==================================== Probing =====================================
     Time  Round  Probed vars  Fixed vars  Mod. entries  New var bds  Implications
    0.002      1            6           0             0            0             3
==================================================================================
After preprocessing: Rows - 5, Cols - 6 (Int - 6, Bin - 6), NZs - 12
Scaling...
Min exp = 0,  Max exp = 0

Preprocessing Time: 0.005
Optimizing...
Building initial basis...
     0.006          1          0      D/S                  5.0000
Generating cuts...
===========================================
MIPCL version 2.1.3
Solution time: 0.009
Branch-and-Cut nodes: 1
Objective value: 5.0000 - optimality proven
----------------------------------------
[出力]
Optimal objective value is 5.0000
'x(Aさん,東京)' = 0.0000
'x(Aさん,大阪)' = 1.0000
'x(Bさん,東京)' = 0.0000
'x(Bさん,大阪)' = 0.0000
'x(Cさん,東京)' = 1.0000
'x(Cさん,大阪)' = 0.0000
--------------------------------
目的関数値 = 5.0
Aさんに大阪を割り当て
Cさんに東京を割り当て
----------------------------------------

ベンチマーク

グラフクラスタリング(与えられた無向グラフを、いくつかの密な部分グラフに分割する問題)に関するこちらの論文の定式化を、いくつかのソルバーで解いてみました。環境はCore i9-9900Kです。
この論文の定式化に基づきインスタンスを解いた印象としては、最適解を見つけるのにも少し時間がかかるものの、それ以上に、見つかった解が最適であることを保証するために多くの時間がかかるような定式化であるということが挙げられます。

インスタンス Gurobi 8.1.0 MIPCL 2.1.3 PuLP 1.6.9 (COIN-CBC, 1スレ)
1 0.6 0.9 21.4
2 0.4 1.2 11.1
3 0.6 2.1 39.4
4 0.3 1.1 16.0
5 9.8 126.8 1,589.4

(単位:秒)
(おおむね、インスタンスの数字が大きいほど、インスタンスのサイズが大きい)

残念ながら私の環境ではPuLPに同こんされているCOIN-CBCがマルチスレッドで動かなかったためフェアな比較にはなっていませんが、MIPCLは、PuLPよりも速く、ある程度のサイズまでは有償ソルバーであるGurobiとそんなに変わりない時間でインスタンスを解きました。

おわりに

MIPCLは、インストール方法だったりAPIが、ほかのソルバーと大きく違って直感的ではないため、もにょるところはありますが(例:Problemクラスにクラス変数curProbがあって、それが現在選択されている問題例を示していること)、フリーで商用利用可能なことと性能のバランスがとれていることから、試しに動かしてみる価値はあると思います。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away