最適化におけるPython（Python-MIP版）

はじめに

この記事は、「最適化におけるPython（PuLP版）」で使っているモデラーのPuLPをPython-MIPに置き換えたものです。
どちらのモデラーもデフォルトのソルバーはCBCで使い勝手もほぼ同じですが、Python-MIPは下記の利点があります。

• ソルバーとのインターフェースにCFFIを用いていて高速です。
• PyPyが使えます。Python部分の実行速度も向上できる可能性があります。
• ベクトル（NumPyの多次元配列）で制約条件を書けます。PuLPではforを使わないといけない制約条件をforを使わずにシンプルに書けます。サンプル
• PuLPよりわかりやすくモデルを作成できます。
  • 目的関数の設定でobjectiveを明記してわかりやすい。
  • 最大化か最小化を目的関数で指定できてわかりやすい。
  • 目的関数の設定が=なのでわかりやすい（PuLPでは、追加でないのに+=）。
  • モデルのクラス名がModelでわかりやすい（PuLPでは、LpProblem）。
  • 多次元配列で変数を作成でき、制約条件もブロードキャストで書ける。
  • 変数名を省略できる。また、同じ変数名でもエラーにならず、ちゃんと解ける（PuLPで同じ変数名を使うと意味不明のエラーメッセージが出る）。
  • 式 != 式を制約条件にするとエラーになる（PuLPはならない）。

概要

私は、業務で、組合せ最適化技術を用いたソフトウェア開発(例えば、物流における輸送コストの最小化など)を行っています。以前は、C++やC#を用いて、最適化のモデルを作成していましたが、最近ではPythonを用いることが多いです。
ここでは、最適化におけるPythonについて紹介します。

Pythonのメリット

Pythonを利用している理由としては、以下のような点があげられます。

• わかりやすい。数式によるモデルとPythonによるモデルが近いため、より本質的な記述に専念でき、保守しやすいモデルを作成できる。

• 短い記述量で済む。C++などに比べるとプログラムのサイズは、数分の1になる。
• 学習コストが小さい。シンプルな文法で、予約語も少ない。
• Pythonで完結できる。汎用言語であるため、種々の目的の処理もほぼPythonで記述できる。
  例えば、webからデータを取得して、集計して、分析して、最適化して、可視化するなど、すべてPythonでできる。
• ライブラリが多い。パッケージコミュニティサイト（https://pypi.org/ ）で約40万ものパッケージが公開されている。
• 様々な環境で実行できる。Windows、Mac、Linuxの各種OSや、処理系もCythonやPyPyなどがある。
• 多くの最適化ソフトウェアがPythonに対応している。有料、無料含めて、多くの最適化ソフトウェアが存在しているが、Pythonから利用できるものが多い。

Pythonは、C++などのコンパイラ言語に比べると、実行速度が遅いと言われます。しかし、最適化においては、主にモデルの作成(モデリング)にPythonを用い、最適化アルゴリズムの実行はC++などで記述された専用ソフトウェア(ソルバー)を用います。このため、最適化でPythonを利用しても、実行時間はあまり問題となりません。

最適化のモデリングでは、主にPython-MIPとpandasパッケージを用いています。

• Python-MIPは、数理モデリングのパッケージであり、pandasはデータ分析のパッケージである。
• pandasは、モデルに含まれるデータの中で、表で表現できるデータを扱うのに適しており、複雑な処理をわかりやすく記述できる。
  また、Python-MIPとpandasは内部でNumPyを利用している。
• NumPyは、CやFortranで書かれた高度に最適化された線形代数ライブラリを使用しており行列計算を効率よく計算することができる。

Python-MIPについて

数理最適化問題を解くためには、以下のステップを行います。

• モデラーで数理モデルを作成します
• ソルバーをよび出して、解を得ます

ソルバーは、数理モデルを入力とし、数理モデルを解いて、変数の値（解）を出力とするソフトウェアです。

Python-MIPは、モデラーになります。
Python-MIPでは、ソルバーとしてCBC,Gurobiなどいろいろなものが使えます。
デフォルトでは、CBCが使われます。Python-MIPをインストールすると、CBCも同時にインストールされます。

• CBC: COINプロジェクトの無料ソルバー(COIN-OR Branch and Cut の略) https://projects.coin-or.org/Cbc
• Gurobi: 高性能な商用ソルバー http://www.gurobi.com/

Python-MIPで扱うことができる問題は、混合整数最適化問題です。
混合整数最適化問題は、数理最適化問題の1種で、下記の特徴があります。

• 連続(実数)変数と整数変数を使って表現される
• 目的関数と制約条件が1次式である

さらに詳細について調べたい場合は、参考サイトを参考にしてください。

Python-MIPの使い方

下記の問題を考えてみましょう。

問題

材料AとBから合成できる化学製品XとYをたくさん作成したい。
Xを1kg作るのに、Aが1kg、Bが3kg必要である。
Yを1kg作るのに、Aが2kg、Bが1kg必要である。
また、XもYも1kg当りの価格は100円である。
材料Aは16kg、Bは18kgしかないときに、XとYの価格の合計が最大になるようにするには、
XとYをどれだけ作成すればよいか求めよ。

問題を数理モデルであらわすと下記のようになります。数理モデルを式で表現することを定式化するといいます。

これをPython-MIPでモデル化して解いてみます。

python

from mip import Model, maximize, minimize, xsum

m = Model()  # 数理モデル
# 変数
x = m.add_var("x")
y = m.add_var("y")
# 目的関数
m.objective = maximize(100 * x + 100 * y)
# 制約条件
m += x + 2 * y <= 16  # 材料Aの上限
m += 3 * x + y <= 18  # 材料Bの上限
m.optimize()  # ソルバーの実行
print(x.x, y.x)
>>>
4.0 6.0

以下、順番に簡単に説明します。

パッケージのインポート

from mip import Model, maximize, minimize, xsum

数理モデルの作成

m = Model()

変数の作成

x = m.add_var(変数名)

主な変数の作り方です。

• 非負の連続変数：m.add_var(変数名)
• 自由変数：m.add_var(変数名, lb=-np.inf)
• 0-1変数：m.add_var(変数名, var_type="B")
• 非負の整数変数：m.add_var(変数名, var_type="I")
• n個の非負の連続変数のベクトル：m.add_var_tensor((n,), 変数名)
• 非負の連続変数の多次元配列：m.add_var_tensor(形状, 変数名)
• n個の0-1変数のベクトル：m.add_var_tensor((n,), 変数名, var_type="B")
• 0-1変数の多次元配列：m.add_var_tensor(形状, 変数名, var_type="B")

※ 形状は、多次元配列のshapeを意味します。1次元（ベクトル）の形状は、(n,)のように書きます。

目的関数の設定

# 最小化
m.objective = minimize(式)

# 最大化
m.objective = maximize(式)

制約条件の追加

m += 式 == 式
m += 式 <= 式
m += 式 >= 式

式の例

式は、数式のように書けます。また、add_var_tensor()の戻り値は、NumPyの多次元配列の派生クラス（LinExprTensor）なので、NumPyと同様に記述できます。
和や内積は、sum()も使えますが、xsum()の方が効率的です（参考「数理モデルにおける変数の和」）。

2 * x + 3 * y - 5

# 和
xsum(変数のベクトル)

# 内積
xsum(係数のベクトル * 変数のベクトル)

ソルバーの実行

m.optimize()

ソルバーを実行すると、ソルバーのログが出力されますが、下記をするとログ出力を抑制できます。

m.verbose = 0

変数や式や目的関数の値

ソルバーを実行して、 m.status.value == 0 であれば、下記のようにして変数などの値を取得できます。

x.x  # 変数xの値
y.x  # 変数yの値
(2 * x + 3 * y - 5).x  # 式の値
m.objective.x  # 目的関数の値
v.astype(float)  # 変数の多次元配列vの値

Python-MIPとpandasの組合せについて

Python-MIPとpandasを組合せて、pandasの表(DataFrame)で変数(Var)を管理すると、シンプルでわかりやすくモデルを作成できます。

輸送最適化問題を例にしてみてみましょう。

輸送最適化問題

倉庫群から工場群へ部品を搬送したい。輸送費が最小となる計画を求めたい。

• 倉庫群から工場群への輸送量を決めたい → 変数
• 輸送コストを最小化したい → 目的関数
• 各倉庫からの搬出は、供給可能量以下 → 制約
• 各工場への搬入は、需要量以上 → 制約

輸送費		組み立て工場
		F1	F2	F3	F4	供給
倉庫	W1	10	11	18	16	47
	W2	19	15	16	19	42
	W3	17	16	15	15	40
	需要	25	26	20	21

パラメータの設定

必要なパラメータを設定します。(数字は前表と同じ)

import numpy as np
import pandas as pd
from mip import Model, minimize, xsum

rng = np.random.default_rng(0)
nw = 3  # 倉庫数
nf = 4  # 工場数
供給 = rng.integers(30, 50, nw)
需要 = rng.integers(20, 40, nf)
輸送費 = rng.integers(10, 20, (nw, nf))

pandasを使わない数理モデル

変数に添字でアクセスします。

m1 = Model()
v1 = m1.add_var_tensor((nw, nf), "v1")
expr = xsum(xsum(v) for v in 輸送費 * v1)
m1.objective = minimize(expr)
for i in range(nw):
    m1 += xsum(v1[i]) <= 供給[i]
for j in range(nf):
    m1 += xsum(v1[:, j]) >= 需要[j]
m1.verbose = 0
m1.optimize()
v1.astype(float)
>>>
LinExprTensor([[25., 22.,  0.,  0.],
               [ 0.,  4.,  1.,  0.],
               [ 0.,  0., 19., 21.]])

pandasを使った数理モデル

変数を表に持つことで、変数を表の属性でアクセスできます。まずは、表を作成しましょう。

dfw = pd.DataFrame({"倉庫": ["W1", "W2", "W3"], "供給": 供給})
dff = pd.DataFrame({"工場": ["F1", "F2", "F3", "F4"], "需要": 需要})
df = pd.merge(dfw, dff, "cross").assign(輸送費=輸送費.flatten())
df

	倉庫	供給	工場	需要	輸送費
0	W1	47	F1	25	10
1	W1	47	F2	26	11
2	W1	47	F3	20	18
3	W1	47	F4	21	16
4	W2	42	F1	25	19
...	...	...	...	...	...

同様に数理モデルを作って解いてみましょう。

m2 = Model()
df["Var"] = m2.add_var_tensor((len(df),), "v2")
m2.objective = minimize(xsum(df.輸送費 * df.Var))
for _, gr in df.groupby('倉庫'):
    m2 += xsum(gr.Var) <= gr.供給.iloc[0]
for _, gr in df.groupby('工場'):
    m2 += xsum(gr.Var) >= gr.需要.iloc[0]
m2.verbose = 0
m2.optimize()
df["Val"] = df.Var.astype(float)
df[df.Val > 0]

	倉庫	工場	輸送費	供給	需要	Var	Val
0	W0	F0	10	47	25	v2_0	25
1	W0	F1	11	47	26	v2_1	22
5	W1	F1	15	42	26	v2_5	4
6	W1	F2	16	42	20	v2_6	1
10	W2	F2	15	40	20	v2_10	19
11	W2	F3	15	40	21	v2_11	21

添え字を使った表現は、添え字が何を表しているか覚えていないといけませんでした。しかし、Python-MIPとpandasを組合せることによって、下記のように、数理モデルが理解しやすくなります。

• 実務ではスパース（疎）なデータが多い。表にすることでスパースなデータでも扱いやすくなる。
• 変数が持つ属性が行を見ればわかる。
• 単なるiとかではなく、倉庫などの列名が使える。
• pandasの条件式を使って、数式を組み立てられる。(参考 組合せ最適化でN Queen問題を解く)
• pandasの便利な関数(groupbyやmergeやjoinなど)が使える。

参考サイト

• Qiitaの記事
  • 組合せ最適化を使おう
  • 数理モデルにおける変数の和
  • 組合せ最適化ソルバーの威力
  • Python-MIPのLinExprTensor.astype()をDataFrameの列に追加したときの考察
  • 双対問題を調べる
  • 今度こそ？使い物になるフリーの数理最適化（混合整数最適化）ソルバー（付きインターフェース） Python-MIP
• Qiita以外の記事
  • Python-MIPによるモデル作成方法
  • 組合せ最適化(松井先生)(PDF 2ページ)
  • 組合せ最適化による問題解決の実践的なアプローチ(梅谷先生)(PDF 5ページ)
• 書籍
  • 「Pythonで学ぶ数理最適化による問題解決入門」
  • 「データ分析ライブラリーを用いた最適化モデルの作り方」
  • 「今日から使える!組合せ最適化」
  • 「Python言語によるビジネスアナリティクス」
  • 「モデリングの諸相 (シリーズ:最適化モデリング)」
• ソルバー関連
  • Python-MIP ドキュメント
  • 整数計画法メモ(宮代先生)
  • 整数計画法による定式化入門
  • 整数計画ソルバー入門
  • Gurobi Optimizer：ソルバーエンジン

---

数理最適化に関するコンサルティングや開発をしています。詳しくは、下記を参照してください。

以上