More than 1 year has passed since last update.

今度こそ？使い物になるフリーの数理最適化（混合整数最適化）ソルバー（付きインターフェース） Python-MIP

Last updated at 2022-11-10Posted at 2021-09-04

(2022/11/10 追記)

Python-MIP付属のCOIN-CBCよりもかなり速い数理最適化ソルバーとして知られるSCIPが無償で商用利用可（少なくともSCIPの中身を変更しないで利用する場合はソースの開示義務なし）になったそうです。

参考：https://qiita.com/samuelladoco/items/703bf78ea66e8369c455

現在、Python-MIPからSCIPを利用することはできません。方や、PuLPからは利用できます。Python-MIPからでもSCIPを利用できるようになる予定だそうなので、それまでの間はPython-MIPではなくPuLPを使うほうがよさげです。

参考：https://github.com/ohmsha/PyOptBook/issues/15

また、最近SciPyパッケージにも最適化の機能が搭載されたようですが、整数最適化に関して言えば scipy.optimize.milp よりもSCIPが速い（SCIPはCOIN-CBCしより速いので）ですし、 scipy.optimize.milp はモデリング部分のコードが面倒になる（データを行列形式で与えなければならない？）ので、Python-MIPやPuLP（ソルバーとしてSCIPを使用）するのがよいです。
参考1：https://mirucacule.hatenablog.com/entry/2022/08/06/190000
参考2：Mittlemanのベンチマーク結果にそこら辺の比較があったはず

(2021/09/06 追記)

本記事の末尾に、大阪大学梅谷先生のツイッターのつぶやきを転記し、私の説明も追加しました。先生のおっしゃることはもっともなので、お読みください。また、記事中の一部の記述を修正しました。

主張

Pythonにて無償で商用利用可能なツールを使い、数理最適化問題の一部である混合整数最適化問題（MIP）・整数最適化問題（IP）・線形最適化問題（LP）を解きたい場合、PuLPパッケージを使うことはよく知られていますが、Python-MIPのほうが絶対よいです。今までPuLPを使っていた方はそれを捨ててPython-MIPを使ってください。これから始める方はPuLPに手を出さずにPython-MIPを無償で商用利用可能なソルバーのファーストチョイスにしてください。

なぜ

PuLPのものよりもPython-MIPに付属している混合整数最適化ソルバーのほうがよいからです
- 付属しているのは両者ともCOIN-CBCですが、Python-MIP付属のもののほうが、ビルド日が新しい・マルチスレッドに対応・ファイル経由ではなくDLLとやりとり、と言う点でPuLP付属のものより優れています
(2021/09/06 修正) ある環境においてPython-MIPの文法でコードを書いて、（有償である）Gurobi Optimizerが導入されている別の環境においてそれを動かす場合、コードを変えなくても後者の環境では自動的にソルバーとしてGurobiが使用されるという、開発上や運用上の楽さがあります
- Gurobi Optimizerの有効なライセンスファイルを配置し、pipやcondaからではなく公式インストーラーからインストールし、環境変数GUROBI_HOMEまたはPATHが設定されている場合です（インストーラーからインストールすれば自動的に設定されます）
PuLPよりもPython-MIPのほうが公式ドキュメント（APIリファレンス）がしっかり整備されています
PuLPよりもPython-MIPのほうが込み入った機能（lazy constraintsなど）が多く実装されているようです

インストール～実行（ベンチマーク）して出力を得るまで

想定環境

Python 3.8
Windows 10 64bit

インストール

pipでmipというパッケージをダウンロード＆インストールすればよいだけです。

具体的には、コマンドプロンプトなどにて

pip.exeにパスが通っている場合は pip install mip
python.exeにパスが通っている場合は python -m pip install mip
パス？なにそれという方は、ptyhon.exeがインストールされているフォルダーをShift + 右クリックし「Powershell ウィンドウをここで開く(S)」を選択し .\python.exe -m pip install mip

サンプルコード

2都市間の距離が非対称な巡回セールスマン問題を解くことにします。巡回セールスマン問題のMIP定式化はいくつかあるのですが、今回はよく知られた定式化の中では比較的速く解ける(single commodity) flow formulation定式化を採用しました。

flow formulationの参考リンク先

以下のコードは、上記の参考リンク先やテキストなどに記載されている定式化とは少し違いますが、たぶんあっていると思います。

sample.py

# Import
# -----------------------------------------------------------------------------
from __future__ import annotations
import random
import time
from typing import Tuple
# 
import mip
# -----------------------------------------------------------------------------


# Main
# -----------------------------------------------------------------------------
# 
print(f'Input generation \n')
# ----------------------------------------------------------------------
# # Vertices
n: int = 300
#
# Vertices
V: list[int] = [i for i in range(1, n + 1)]
#
# (型ヒントを使用)
VertexType = type(V[0])
ArcType = Tuple[VertexType, VertexType]
# 
# Directed arcs
A: list[ArcType] = [(i, j) for i in V for j in V if i != j]
#
# Distance function
random.seed(1)
d: dict[ArcType, int] = {a: random.randint(1, 99) for a in A}
# 
# (便宜上、最初に出発する都市を設定)
v_first: VertexType = V[0]
# 
# Time limit
time_limit_s: int = 60 * 60
# 
# Debug
# print(f'  V={V}')
# print(f'  A={A}')
# print(f'  d={d}')
# print(f'  v_first={v_first}')
# print(f'  Time limit={time_limit_s}(s) \n')
# ----------------------------------------------------------------------


print(f'Instance generation \n')
# ----------------------------------------------------------------------
formulation_time_start_s: float = time.perf_counter()
# 
# Model
#   ここでソルバーを指定する。ソルバーを指定しない場合、
#   Guorbiを(pipやcondaからではなく)インストーラーから入れ
#   環境変数GUROBI_HOMEやパスが設定されているならば、Gurobiが選ばれる。
#   そうでない場合は、Python-MIPに同こんされているCOIN-CBCが選ばれる。
m: mip.Model = mip.Model(
    name='TSPorHamiltonianPathFF', 
    # solver_name=mip.CBC, 
    # solver_name=mip.GUROBI, 
    )

#   Variables
x: dict[ArcType, mip.Var] = {}
for a in A:
    x[a] = m.add_var(name=f'Var_x({a})', var_type=mip.BINARY)
    del a
#
f: dict[ArcType, mip.Var] = {}
for a in [aa for aa in A if aa[1] is not v_first]:
    f[a] = m.add_var(name=f'Var_f({a})')
    del a

#   Objective
m.objective = mip.minimize(mip.xsum((d[a]) * x[a] for a in A))

#   Constraints
#       x
for v in V:
    m.add_constr(
        mip.xsum(x[a] for a in A if a[1] is v) == 1, 
        name=f'Con_In({v})')
    m.add_constr(
        mip.xsum(x[a] for a in A if a[0] is v) == 1, 
        name=f'Con_Out({v})')
    del v
#
#       f
for v in [vv for vv in V if vv is not v_first]:
    m.add_constr(
           mip.xsum(f[a] for a in A if a[1] is v) 
         - mip.xsum(f[a] for a in A if a[0] is v and a[1] is not v_first) 
        == 1, 
        name=f'Con_flow({v}')
    del v
#
#       f-x
for a in [aa for aa in A if aa[0] is v_first]:
    m.add_constr(f[a] == (len(V) - 1) * x[a], name=f'Con_f-x({a})')
    del a
for a in [aa for aa in A if aa[0] is not v_first and aa[1] is not v_first]:
    m.add_constr(f[a] <= (len(V) - 2) * x[a], name=f'Con_f-x({a})')
    del a
#
formulation_time: float = time.perf_counter() - formulation_time_start_s
del formulation_time_start_s
# ----------------------------------------------------------------------


print(f'Optimization')
# ----------------------------------------------------------------------
optimization_time_start_s: float = time.perf_counter()
#
#   使用スレッド数: 
#       0 = お任せ(デフォルト, 複数スレ使ってくれている?)
#       -1 = 計算機の全スレ(SMTスレも使うので物理コア数だけ使うより非効率な場合も?)
#       1以上 = 指定した数のスレ(指定したぶんだけちゃんと使ってくれている?)
m.threads = -1
status: mip.OptimizationStatus = m.optimize(max_seconds=time_limit_s)
obj_val: float = m.objective_value
is_optimal: bool = (status == mip.OptimizationStatus.OPTIMAL)
#
optimization_time: float = time.perf_counter() - optimization_time_start_s
del optimization_time_start_s
print(f'')
# ----------------------------------------------------------------------


print(f'Solution check')
# ----------------------------------------------------------------------
# x and f
#   x[a] = 1 である枝を集める
xa_1s: list[ArcType] = [
    a for a in A if a in x.keys() and x[a].x > 1.0 - 0.01]
#   f[a]の値の大きい順( = 訪問順)に枝並べ替えするが、
#   最後の枝だけは f[a] が未定義なので、一度抜いて並べ替えて後で最後に差し込む
xa_lasts: list[ArcType] = [xa for xa in xa_1s if xa not in f.keys()]
assert len(xa_lasts) == 1, f'x[a] = 1 で f[a] が未定義の a が 1 つではありません'
xa_1s.remove(xa_lasts[0])
xa_1s.sort(key=lambda a: f[a].x, reverse=True)
xa_1s.append(xa_lasts[0])
del xa_lasts
assert len(xa_1s) == len(V), f'{xa_1s} と 頂点数 {len(V)} が一致しません'
# Debug
# for a in xa_ones:
#     print(
#           f'  arc={a}: ' 
#         + f'x={x[a].x}, f={f[a].x if a in f.keys() else "Undefined"}, ' 
#         + 'd={d[a]}')
#     del a
# print(f'  Sum of d[a]: {sum([d[a] for a in xa_ones])} \n')
#
# 並べ替えた枝が巡回路を構成しているか確かめる
prev_xa: ArcType = xa_1s[0]
for xa in xa_1s:
    # 最初の要素は xa is prev_xa なので検査しない
    if xa is xa_1s[0]:
        continue
    # 次の要素からは検査
    assert xa[0] is prev_xa[1], f'{prev_xa}-{xa} がつながっていません'
    prev_xa = xa
    # 最後の要素は追加検査
    if xa is xa_1s[-1]:
        assert xa[1] is xa_1s[0][0], f'{xa_1s[-1]}-{xa_1s[0]} がつながっていません'
    del xa
del prev_xa
# 
A_sol: list[ArcType] = [a for a in xa_1s]
del xa_1s
V_sol: list[VertexType] = [a[0] for a in A_sol]
V_sol.append(A_sol[-1][1])
assert len(V_sol) == len(V) + 1, f'{V_sol} と 頂点数 {len(V)} + 1 が一致しません'
assert len(V_sol) == len(set(V_sol)) + 1, f'{V_sol} に始点以外で重複要素があります'
sum_d_A_sol: int = sum(d[a] for a in A_sol)
del A_sol
assert abs(obj_val - sum_d_A_sol) < 1.0, f'{obj_val} != {sum_d_A_sol} です'
del obj_val
print(f'')
# ----------------------------------------------------------------------


print(f'Solution')
# ----------------------------------------------------------------------
print(f'  Tour                = {V_sol}')
print(f'  Sum of distance     = {sum_d_A_sol}')
print(f'  Optimality          = {is_optimal}')
print(f'  Time (formulation)  = {formulation_time:6.1f} (s)')
print(f'       (optimization) = {optimization_time:6.1f} (s)')
# ----------------------------------------------------------------------
# -----------------------------------------------------------------------------

実行（ベンチマーク）結果

上記のコードのとおり、完全有向グラフ・都市数300・2都市間の距離が非対称な問題例を解きました。2都市間の距離は適当にランダムで発生させたので、三角不等式を満たしているとは限りません。

実験環境は、CPU: i9-9900K, RAM: 64GB（本計算には、そんなにRAMは使われていません）, OS: Windows 10 64bit版です。

以下、いずれもpipで最新版をダウンロードして実行した結果です。

Pythonパッケージ	Gurobi (gurobipy 9.1.2)	Python-MIP (mip 1.13.0)	PuLP (pulp 2.5.0)
使用ソルバー(パッケージ付属)	Gurobi Optimizer 9.1.2	COIN-CBC 2020/11/15版	COIN-CBC 2019/12/15版
マルチスレッド	対応	対応版	非対応版
定式化にかかった時間（秒）	5.5	7.3	7.3
計算にかかった時間（秒）	257.6	609.4	874.0

今回COIN-CBCを動かしてみて気づいたのですが、解いた問題の性質とCOIN-CBCの特性上？、Python-MIPではマルチスレッドで動かせる分枝限定法が適用される場面は少なく、切除平面法で最適解を求めようとあがいている割合が多いようです。なので、分枝限定法が多く適用される問題に対しては、Python-MIPとCOIN-CBCの付属ソルバーの速度差はより開くと予想します。

サンプルコードの出力

Input generation 

Instance generation 

Optimization
Welcome to the CBC MILP Solver 
Version: devel 
Build Date: Nov 15 2020 

Starting solution of the Linear programming relaxation problem using Primal Simplex

Coin0506I Presolve 90001 (-299) rows, 178802 (-299) columns and 536107 (-598) elements
Clp0030I 13 infeas 3936.1919, obj 0.73453941 - mu 111.08889, its 52, 82471 interior
Clp0030I 23 infeas 336.4839, obj 321.96194 - mu 37.025927, its 52, 89822 interior
Clp0030I 33 infeas 280.74914, obj 331.43035 - mu 4.1131691, its 52, 90301 interior
Clp0030I 43 infeas 269.95626, obj 328.92366 - mu 0.45692739, its 52, 90337 interior
Clp0030I 49 infeas 218.14118, obj 333.89375 - mu 0.1522939, its 105, 89929 interior
Clp0030I 55 infeas 191.28678, obj 335.95861 - mu 0.050759557, its 105, 90927 interior
Clp0030I 60 infeas 114.15094, obj 350.39721 - mu 0.050759557, its 105, 89391 interior
Clp0030I 64 infeas 94.151937, obj 354.25841 - mu 0.01691816, its 105, 93987 interior
Clp0030I 68 infeas 67.847049, obj 361.15958 - mu 0.0056388228, its 105, 95848 interior
Clp0030I 72 infeas 27.07786, obj 383.67043 - mu 0.0056388228, its 105, 91564 interior
Clp0030I 76 infeas 20.252627, obj 384.47393 - mu 0.0018794197, its 105, 100345 interior
Clp0030I 80 infeas 14.024998, obj 390.04209 - mu 0.00062641057, its 105, 103160 interior
Clp1000I sum of infeasibilities 294.952 - average 0.00326636, 76280 fixed columns
Coin0506I Presolve 15018 (-75282) rows, 28340 (-150761) columns and 84972 (-451733) elements
Clp0006I 0  Obj 387.00444 Primal inf 11.549002 (9851) Dual inf 2.9823596e+12 (26490)
Clp0029I End of values pass after 28340 iterations
Clp0014I Perturbing problem by 0.001% of 17.417343 - largest nonzero change 0.00017739902 ( 0.00097137468%) - largest zero change 2.9996244e-05
Clp0006I 28552  Obj 334.1268 Dual inf 7307.6343 (9667)
Clp0000I Optimal - objective value 332
Clp0000I Optimal - objective value 332
Coin0511I After Postsolve, objective 332, infeasibilities - dual 0 (0), primal 0 (0)
Clp0006I 0  Obj 332
Clp0000I Optimal - objective value 332
Clp0000I Optimal - objective value 332
Clp0000I Optimal - objective value 332
Clp0032I Optimal objective 332 - 0 iterations time 11.292, Idiot 11.05

Starting MIP optimization
Cgl0004I processed model has 90001 rows, 178802 columns (89700 integer (89700 of which binary)) and 536107 elements
Coin3009W Conflict graph built in 0.153 seconds, density: 0.042%
Cgl0015I Clique Strengthening extended 0 cliques, 0 were dominated
Cbc0038I Initial state - 464 integers unsatisfied sum - 50.3197
Cbc0038I Pass   1: (21.88 seconds) suminf.    3.42356 (12) obj. 542.247 iterations 50630
Cbc0038I Pass   2: (32.69 seconds) suminf.    2.33557 (10) obj. 510.711 iterations 27483
Cbc0038I Pass   3: (42.73 seconds) suminf.    2.40940 (9) obj. 545.681 iterations 25885
Cbc0038I Pass   4: (45.66 seconds) suminf.    1.71812 (19) obj. 563.975 iterations 6496
Cbc0038I Pass   5: (56.74 seconds) suminf.    0.00000 (0) obj. 814 iterations 23791
Cbc0038I Solution found of 814
Cbc0038I Relaxing continuous gives 814
Cbc0038I Before mini branch and bound, 89204 integers at bound fixed and 88647 continuous
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 50 rows 44 columns
Cbc0038I Mini branch and bound did not improve solution (67.53 seconds)
Cbc0038I Round again with cutoff of 764.9
Cbc0038I Pass   6: (93.79 seconds) suminf.    3.42356 (12) obj. 542.247 iterations 4056
Cbc0038I Pass   7: (104.80 seconds) suminf.    2.33557 (10) obj. 510.711 iterations 31477
Cbc0038I Pass   8: (117.48 seconds) suminf.    2.40940 (9) obj. 545.681 iterations 32174
Cbc0038I Pass   9: (120.56 seconds) suminf.    1.71812 (19) obj. 563.975 iterations 5295
Cbc0038I Pass  10: (122.42 seconds) suminf.    0.50336 (6) obj. 738.809 iterations 689
Cbc0038I Pass  11: (124.08 seconds) suminf.    0.33557 (8) obj. 736.082 iterations 872
Cbc0038I Pass  12: (132.11 seconds) suminf.    2.04027 (7) obj. 764.9 iterations 14826
Cbc0038I Pass  13: (133.34 seconds) suminf.    1.91946 (4) obj. 723.56 iterations 993
Cbc0038I Pass  14: (134.55 seconds) suminf.    0.00000 (0) obj. 725 iterations 1637
Cbc0038I Solution found of 725
Cbc0038I Relaxing continuous gives 725
Cbc0038I Before mini branch and bound, 89200 integers at bound fixed and 88643 continuous
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 53 rows 50 columns
Cbc0038I Mini branch and bound improved solution from 725 to 682 (140.53 seconds)
Cbc0038I Round again with cutoff of 611.2
Cbc0038I Pass  15: (141.02 seconds) suminf.    3.42356 (12) obj. 542.247 iterations 1840
Cbc0038I Pass  16: (151.80 seconds) suminf.    2.33557 (10) obj. 510.711 iterations 32081
Cbc0038I Pass  17: (163.26 seconds) suminf.    2.40940 (9) obj. 545.681 iterations 28799
Cbc0038I Pass  18: (166.54 seconds) suminf.    1.71812 (19) obj. 563.975 iterations 5539
Cbc0038I Pass  19: (170.35 seconds) suminf.    0.81250 (10) obj. 611.2 iterations 5445
Cbc0038I Pass  20: (172.17 seconds) suminf.    0.00000 (0) obj. 593 iterations 1111
Cbc0038I Solution found of 593
Cbc0038I Relaxing continuous gives 593
Cbc0038I Before mini branch and bound, 89202 integers at bound fixed and 88644 continuous
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 77 rows 71 columns
Cbc0038I Mini branch and bound did not improve solution (177.99 seconds)
Cbc0038I Round again with cutoff of 514
Cbc0038I Pass  21: (263.81 seconds) suminf.    2.80027 (15) obj. 514 iterations 18446
Cbc0038I Pass  22: (266.53 seconds) suminf.    2.67114 (23) obj. 495.935 iterations 4203
Cbc0038I Pass  23: (272.54 seconds) suminf.    1.76230 (8) obj. 514 iterations 7706
Cbc0038I Pass  24: (278.98 seconds) suminf.    1.03438 (7) obj. 514 iterations 13838
Cbc0038I Pass  25: (289.76 seconds) suminf.    0.00000 (0) obj. 456 iterations 17375
Cbc0038I Solution found of 456
Cbc0038I Relaxing continuous gives 456
Cbc0038I Before mini branch and bound, 89205 integers at bound fixed and 88648 continuous
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 90 rows 73 columns
Cbc0038I Mini branch and bound did not improve solution (295.60 seconds)
Cbc0038I Round again with cutoff of 418.1
Cbc0038I Reduced cost fixing fixed 10604 variables on major pass 5
Cbc0038I Pass  26: (335.51 seconds) suminf.    2.78498 (17) obj. 418.1 iterations 4154
Cbc0038I Pass  27: (342.22 seconds) suminf.    1.28859 (7) obj. 373.06 iterations 10968
Cbc0038I Pass  28: (345.87 seconds) suminf.    1.28859 (16) obj. 362.161 iterations 5454
Cbc0038I Pass  29: (358.78 seconds) suminf.    2.23334 (6) obj. 418.1 iterations 24590
Cbc0038I Pass  30: (364.74 seconds) suminf.    0.00000 (0) obj. 398 iterations 4126
Cbc0038I Solution found of 398
Cbc0038I Relaxing continuous gives 398
Cbc0038I Before mini branch and bound, 89207 integers at bound fixed and 88648 continuous
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 56 rows 45 columns
Cbc0038I Mini branch and bound did not improve solution (372.36 seconds)
Cbc0038I Round again with cutoff of 371
Cbc0038I Reduced cost fixing fixed 53425 variables on major pass 6
Cbc0038I Pass  31: (415.19 seconds) suminf.    1.28859 (7) obj. 363.916 iterations 10287
Cbc0038I Pass  32: (417.92 seconds) suminf.    1.28859 (12) obj. 360.777 iterations 4378
Cbc0038I Pass  33: (420.21 seconds) suminf.    4.09837 (10) obj. 371 iterations 1833
Cbc0038I Pass  34: (421.85 seconds) suminf.    3.93427 (13) obj. 371 iterations 3637
Cbc0038I Pass  35: (424.69 seconds) suminf.    1.69112 (18) obj. 371 iterations 4147
Cbc0038I Pass  36: (425.08 seconds) suminf.    1.50335 (17) obj. 371 iterations 226
Cbc0038I Pass  37: (432.28 seconds) suminf.    0.00000 (0) obj. 335 iterations 13365
Cbc0038I Solution found of 335
Cbc0038I Relaxing continuous gives 335
Cbc0038I Before mini branch and bound, 89207 integers at bound fixed and 88648 continuous
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 43 rows 34 columns
Cbc0038I Mini branch and bound did not improve solution (438.21 seconds)
Cbc0038I After 462.79 seconds - Feasibility pump exiting with objective of 335 - took 456.39 seconds
Cbc0012I Integer solution of 335 found by feasibility pump after 0 iterations and 0 nodes (498.97 seconds)
Cbc0038I Full problem 90001 rows 178802 columns, reduced to 854 rows 743 columns
Cbc0012I Integer solution of 332 found by RINS after 0 iterations and 0 nodes (543.55 seconds)
Cbc0030I Thread 0 used 0 times,  waiting to start 3.5513699, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 1 used 0 times,  waiting to start 3.3475089, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 2 used 0 times,  waiting to start 3.11816, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 3 used 0 times,  waiting to start 2.8999031, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 4 used 0 times,  waiting to start 2.6695158, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 5 used 0 times,  waiting to start 2.439141, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 6 used 0 times,  waiting to start 2.208766, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 7 used 0 times,  waiting to start 1.99402, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 8 used 0 times,  waiting to start 1.7480149, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 9 used 0 times,  waiting to start 1.5166428, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 10 used 0 times,  waiting to start 1.2862649, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 11 used 0 times,  waiting to start 1.071516, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 12 used 0 times,  waiting to start 0.84114099, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 13 used 0 times,  waiting to start 0.61176586, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 14 used 0 times,  waiting to start 0.39601684, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Thread 15 used 0 times,  waiting to start 0.1510241, 0 cpu time, 0 locks, 0 locked, 0 waiting for locks
Cbc0030I Main thread 0 waiting for threads,  1 locks, 0 locked, 0 waiting for locks
Cbc0001I Search completed - best objective 332, took 0 iterations and 0 nodes (548.24 seconds)
Cbc0035I Maximum depth 0, 75886 variables fixed on reduced cost
Total time (CPU seconds):       597.80   (Wallclock seconds):       597.80


Solution check

Solution
  Tour                = [1, 161, 84, 200, 242, 147, 87, 93, 121, 229, 173, 184, 50, 280, 180, 76, 217, 98, 48, 11, 254, 155, 248, 192, 16, 95, 46, 68, 255, 277, 206, 34, 169, 263, 274, 246, 80, 129, 168, 160, 125, 131, 162, 52, 191, 175, 176, 
113, 8, 149, 203, 23, 72, 186, 4, 114, 148, 105, 269, 279, 201, 209, 65, 110, 188, 36, 86, 56, 70, 196, 38, 61, 137, 177, 27, 212, 39, 150, 145, 2, 281, 223, 247, 51, 74, 132, 55, 179, 31, 20, 287, 25, 264, 172, 189, 91, 30, 166, 285, 215, 244, 238, 140, 290, 17, 232, 152, 130, 53, 60, 117, 90, 193, 216, 126, 252, 241, 120, 47, 268, 224, 42, 111, 96, 41, 154, 289, 82, 37, 75, 214, 204, 265, 272, 79, 159, 15, 165, 276, 83, 40, 211, 6, 99, 237, 108, 181, 67, 221, 135, 258, 187, 164, 
134, 251, 78, 296, 115, 28, 153, 22, 295, 94, 133, 167, 112, 190, 102, 259, 245, 202, 250, 142, 205, 3, 178, 7, 104, 88, 151, 92, 9, 136, 10, 109, 43, 122, 144, 292, 26, 197, 124, 194, 283, 5, 236, 89, 128, 174, 199, 101, 240, 271, 170, 219, 73, 267, 286, 146, 298, 294, 54, 300, 208, 100, 24, 49, 118, 123, 262, 256, 228, 158, 57, 85, 218, 156, 81, 58, 64, 243, 249, 77, 266, 233, 275, 288, 210, 198, 185, 213, 222, 29, 299, 278, 139, 63, 143, 71, 282, 119, 195, 141, 226, 284, 270, 35, 227, 291, 261, 171, 97, 293, 183, 138, 260, 231, 69, 106, 225, 116, 103, 21, 59, 19, 13, 66, 157, 33, 253, 239, 235, 12, 32, 182, 230, 163, 107, 297, 234, 220, 62, 14, 44, 273, 127, 18, 257, 207, 45, 1]
  Sum of distance     = 332
  Optimality          = True
  Time (formulation)  =    7.3 (s)
       (optimization) =  609.4 (s)

(2021/09/06 追記) 大阪大学梅谷先生のツイッターのコメントの転記

ソルバー(=アルゴリズム)とインターフェース(=モデリング)を区別しないと読者を混乱させるので気をつける必要があります．COIN-CBCやGurobiがソルバー，PuLPやPython-MIPがインターフェースです． https://t.co/bN98kRCPDt
— Umepon (@shunji_umetani) September 5, 2021

はい、そのとおりです。Python言語の場合を例に話をしますと、pipでインストールできるgurobipy, mip (Python-MIPのパッケージ名), pulpがモデリング言語（モデリングインタフェース、ソルバーとのインターフェース、ソルバーのAPI）にあたり、それらのモデリング言語を用いて記述された数理最適化問題を受け取り解こうとしてくれる計算エンジンがソルバーになります（optimize()やsolve()といったメソッドを呼ぶと何やらログが吐き出されながらCPUがぶん回りますが、それを起こしているものです）。ただし、ダウンロードしたものにモデリング言語機能だけが入っていて肝心の最適化の計算ができずソルバーは別途インストールが必要という話だと不便なので、gurobipyにはGurobi Optimizer（有償ライセンスを契約していなければ小さいサイズしか解けない試用版ライセンスが適用される）、mipやpulpには（バージョンなどが異なる）COIN-CBCというように、何らかのソルバーが付属しているモデリング言語が多いです。特に指定をしなければ、デフォルトで付属しているソルバーがモデリング言語により記述された数理最適化問題を受け取り計算を開始します。

「Gurobi Optimizerがインストールされている環境では、コードを変えなくても自動的にソルバーとしてGurobiが使用される」が最大の利点ならば，Gurobiは商用ソルバーなので，タイトルの「フリーの数理最適化ソルバー」はミスリーディングを引き起こします．
— Umepon (@shunji_umetani) September 5, 2021

記述を「(2021/09/06 修正) ある環境においてPython-MIPの文法でコードを書いて、（有償である）Gurobi Optimizerが導入されている別の環境においてそれを動かす場合、コードを変えなくても後者の環境では自動的にソルバーとしてGurobiが使用されるという、開発上や運用上の楽さがあります」に修正しました。意図を表すには言葉足らずでした。

想定として、まず開発において、（最安である）マシン固定ライセンスでGurobiが使用可能な状態にある重い計算用や本番想定用のマシン（CPUがよいPCやサーバー）と、開発者（複数人もあり得る）がコーディング～簡単なテストまでを行うマシン（Gurobiはインストールされていない）が異なる場合を念頭に置いてました。運用においては、確かに、開発会社（SIerやコンサル会社）の開発機と運用会社（事業会社）の運用機にてそれぞれGurobiが導入されているケースが多そうですが、請負契約でも客先にて客が保有するGurobiライセンスを使用して開発をした後、瑕疵対応に備えてソルバー以外は自社で同様の環境を構築しておくことや、同一会社内で開発から運用をする場合に1ライセンスだけGurobiを保有しライセンスが割り当てられたマシンを開発機から本番機に移行した後も開発機で最適化問題が解ける状態を維持することが考えられます（本番機が障害を起こしその原因がソルバーの違いに起因する場合、原因の特定が困難になるというリスクはありますが）。

ちなみに，GurobiはもともとPythonをインターフェースとして開発しており，GurobiユーザにとってPuLPやPython-MIPを使うことでモデリングが楽になることはありません（おそらくサブセット）．利点は，ソルバーを切り替える手間が省けることぐらいです．
— Umepon (@shunji_umetani) September 5, 2021

おっしゃるとおりです。

実務では事例ごとにモデルを作成するので，コピペだけで他の事例に対応できることは少ないです．開発時はGurobiで運用時はCBCという使い分けはしませんから，Gurobiが利用できる状況でわざわざPuLPやPython-MIPを使う動機は多くはないと思います．
— Umepon (@shunji_umetani) September 5, 2021

「Gurobiが利用できる状況でわざわざPuLPやPython-MIPを使う動機は多くはない」のは、おっしゃるとおりです。2つ上の記述と重複しますが、コーディング～簡単なテストまでを手元端末のCOIN-CBCで、重い計算や総合テスト・それ以降の計算を行う別環境でGurobiを使う場合を想定しています。

ソルバーにGurobiではなくCOIN-CBCを使うならば，PuLPよりもPython-MIPの方が使い勝手が良いのはその通りだと思います．
— Umepon (@shunji_umetani) September 5, 2021

この記事でいいたかったのはまさにこれです。

正確にはGurobiなどのソルバーはインターフェースも備えています．主旨は「Gurobiのインターフェースは使い勝手が悪くないので，他に動機がなければGurobiユーザがサードパーティのインターフェースを使う必要はなさそう」です． https://t.co/MHMcoRhSsV
— Umepon (@shunji_umetani) September 5, 2021

おっしゃるとおりです。（まだ開発版ですが）今年の4月にGurobi開発者の手によるgurobipy-stubsパッケージが登場したことで（pip install gurobipy-stubsでインストール）、Visual Studio Codeなどの統合開発環境（IDE）でGuorbi OptimizerのPython APIのコード補完が効くようになったので、利便性はさらに増しました（他の言語はもともとコード補完が効いていたように記憶しています）。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up