Python
最適化
pulp

Python+PuLPによるタダで仕事に使える数理最適化

More than 1 year has passed since last update.

はじめに

数理最適化とか数理計画とか呼ばれる、こんな数学の問題

問題(1): 線形最適化問題(線形計画問題)
\begin{alignat}{2}
 &\mbox{最大化} 
 &\qquad x + y & \\
 &\mbox{制約} 
 & 2x + y &\leq 2 \\
 &
 & x + 2y &\leq 2 \\
 &
 & x &\geq 0 \\
 &
 & y &\geq 0 \\
\end{alignat}

a - コピー.png

や、こんな数学の問題

問題(2): 整数最適化問題(整数計画問題)
\begin{alignat}{3}
 &\mbox{Minimize} 
 &\qquad \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & 
 &\qquad & \\
 &\mbox{subject to} 
 & \sum_{j \in J} x_{ij} &\leq 1 
 & &\forall i \in I \\
 &
 &\sum_{i \in I} x_{ij} &= 1
 & &\forall j \in J \\
 &
 & x_{ij} &\in \{0,1\}
 & &\forall i \in I, \quad \forall j \in J
\end{alignat}

を解いてくれる、PythonのPuLPパッケージを紹介します。

対象ですが、目的関数と制約式が線形式(1次式)、変数が連続値または離散値またはそれらが混在したもので記述できるようなものとします。

余談ですが、問題(2)は、整数線形最適化問題(整数線形計画問題)と呼ばれることもあります。また、連続変数と離散変数が混在した問題は、混合整数最適化問題(混合整数計画問題、混合整数線形最適化問題、混合整数線形計画問題)と呼ばれます。

PuLPとは

PuLP自体は、Python上で数理最適化問題を簡単にコーディングできるモデリングAPIモデリング言語だと思ったほうがよいです。実際に数式を解いてくれるソルバーとは別物なのですが、COINというソルバーが同こんされているので、PuLPをインストールするだけで問題を解けるようになります。PuLP自体はCOIN以外にもいくつかのソルバーに対応しているので、それらのソルバーが別途インストールされていれば、PuLPから呼び出すソルバーを(簡単に)切り替えることができます。

参考

先人たちのQiitaエントリー

本家サイト

なぜPython+PuLPなのか

あくまで個人的な見解ですが…

お値段

Gurobi OptimizerとかIBM ILOG CPLEX Optimization StudioとかFICO Xpress Optimization SuiteとかNumerical Optimizerとかは高いし…
 → PuLPは無料!同こんされているソルバーCOINも無料!!

お速さ

GLPKやlp_solveは遅いし…Microsoft Excelは遅いうえに扱える変数や制約式の数も少ないし…
 → PuLPに同こんされているCOINはフリーのソルバーでは悪くないほう!ちょっとコードを変えるだけで商用ソルバーも呼べる!!

お商用利用

SCIPは商用利用したけりゃメールくださいって書いてあるし…
 → PuLP & COINは商用利用OK!組み込んだ製品のソースコードの公開義務もない!!

おモデル記述言語

C言語で有名なK&RのKさんが開発にかかわったというモデリング言語AMPLや、IBM ILOG CPLEX Optimization StudioとかFICO Xpress Optimization SuiteとかNumerical Optimizerとかにもある独自のモデリング言語、あれって数式みたいにモデルを記述できるからわかりやすいよね。でもほかじゃ使えないし、込み入ったことをやろうとするとマニュアルを見なきゃいけないし、ググっても情報少ないし…
 → PuLPも数式みたいにモデルを書ける!Pythonなのでほかのアルゴリズムとの連携も簡単だし、わからなくなったらググれば情報がたくさん!!モデル部分の記述はそのままにして呼び出すソルバーだけ変えられるので、実質PuLPが共通言語になり得る!!!

お開発環境

IBM ILOG CPLEX Optimization StudioとかFICO Xpress Optimization Suiteについてくる統合開発環境、いろいろ情報が表示できるしマウスでポチポチできて便利だよね。でもほかじゃ使えないし…
 → PythonなのでJupyter NotebookやSpyderやVisual Studio(タダであるCommunityエディションやExpressエディションでも)上で使える!いろんなマシーンや環境で動く!!誰でもどこでも共同開発やテストができる!!!


研究者さんや学生さんへ

研究者さんや学生さんは、研究や学習を用途として、商用ソルバーが安価にあるいは無償で手に入るかもしれません。

研究者さんは、研究内容が将来的に商用実用化に向かうもののようでしたら、その際にかかるコストも考慮して、今のうちからPython+PuLPの使用を検討ください。何回も書いてますが、PuLPからコールするソルバーは簡単に切り替えられるので、商用ソルバーの深いところ(整数最適化のパラメーターのチューニングなど)まで触る必要がないなら、モデリングAPIとしてPuLPは十分だと思います。また、Python用APIを持つ商用ソルバーからの移植や両立についても、必要となるコーディング量はそれほど多くはなりません。

学生さんは、現在のバイト先あるいは将来の就職先やその客先などで数理最適化を活用することを考えて、今のうちからPython+PuLP(COIN)の使用を検討ください。もちろん案件の規模にもよりますが、商用ソルバーの購入や保守の費用を捻出させるだけの費用対効果を説明するのはけっこうハードル高いです。企業からの委託研究や共同研究などをしている場合でも同じです。企業の人が使える選択肢を考えてください。


以上、個人的な見解でした。
※ 最近でてきたフリーのソルバーMIPCLのことはよくわからないので、触れていません。すみません。

想定環境

以下、(サポート中の)WindowsPython 3の場合を想定して書きます。MacやLinuxなどでは適宜置き換えてください。Python 2と3ではprint文の書き方などが異なっているので注意してください。

インストール

事前かくにん!

自分でインストールしたPythonのほかに、過去にGurobi Optimizerやシミュレーションソフトなどをインストールしていたならば、Pythonの処理系がひっついてきている可能性があります。その場合は、コマンドプロンプトで単にpythonpipと打ったときに、どの処理系が呼ばれるか、環境変数Pathの並び順を見て確認し、必要があれば並び替えましょう(もちろん、毎回C:\Program Files (中略) \python.exeなどとフルパスで書いてもよいですが)。python -Vと打って出力されるバージョンの数値を確認すると、どの処理系が呼ばれているかのヒントになります。

よかった。 からの実作業

(場合によっては環境変数の並び替えをした後で)新規にコマンドプロンプトを起動します。Pythonが管理者権限での書き込みが必要なフォルダーにインストールされている場合は、「管理者として実行」で起動します。

企業なんかでプロキシーに邪魔される場合は、以下のようにして一時的にコマンドプロンプトにプロキシーを設定すればいけそうです。

set HTTP_PROXY=111.222.111.222:3333
set HTTPS_PROXY=111.222.111.222:3333

参考:http://d.hatena.ne.jp/showhey810/20140905/1409892787

pip.exeのあるフォルダーにパスが通っている場合は
pip install pulp
で、そうではなくてpython.exeのあるフォルダーにパスが通っている場合は
python -m pip install pulp
でインストールします。

アップデートのしかたは
pip install -U pulp
または
python -m pip install -U pulp
です。

例(1)

冒頭の問題(1)を解きます。

\begin{alignat}{2}
 &\mbox{最大化} 
 &\qquad x + y & \\
 &\mbox{制約} 
 & 2x + y &\leq 2 \\
 &
 & x + 2y &\leq 2 \\
 &
 & x &\geq 0 \\
 &
 & y &\geq 0. \\
\end{alignat}

コード

pulp_problem_1.py
# coding: UTF-8

# 線形/整数最適化問題を解くためにPuLPをインポート
import pulp
# sys.maxsizeでintの最大値を取得するためにインポート
import sys

# 数理最適化問題(最大化)を宣言
problem = pulp.LpProblem("Problem-1", pulp.LpMaximize)

# 変数(連続)を宣言
x = pulp.LpVariable("x", 0, 999, pulp.LpContinuous)
y = pulp.LpVariable("y", 0, sys.maxsize, pulp.LpContinuous)

# 目的関数を宣言
problem += ( x + y, "Objective" )
# かっこは不要だけど、わかりやすさのため記載

# 制約条件を宣言
problem += ( 2 * x + y <= 2 , "Constraint_1" )
problem += ( x + 2 * y <= 2 , "Constraint_2" )

# 問題の式全部を表示
print("問題の式")
print("--------")
print(problem)
print("--------")

# 計算
result_status = problem.solve()

# (解が得られていれば)目的関数値や解を表示
print("")
print("計算結果")
print("********")
print("最適性 = {}".format(pulp.LpStatus[result_status]))
print("目的関数値 = {}".format(pulp.value(problem.objective)))
print("解 x = {}".format(pulp.value(x)))
print("  y = {}".format(pulp.value(y)))
print("********")

出力

問題の式
--------
Problem-1:
MAXIMIZE
1*x + 1*y + 0
SUBJECT TO
Constraint_1: 2 x + y <= 2

Constraint_2: x + 2 y <= 2

VARIABLES
x <= 999 Continuous
y <= 9.22337203685e+18 Continuous

--------

計算結果
********
最適性 = Optimal
目的関数値 = 1.33333334
解 x = 0.66666667
  y = 0.66666667
********

余談(1)

演算子のオーバーロードが使える言語は、数理最適化の式をラクに記述できますね。そうでない言語は…Javaさんは……


例(2)

冒頭の問題(2)を解きます。

\begin{alignat}{3}
 &\mbox{Minimize} 
 &\qquad \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & 
 &\qquad & \\
 &\mbox{subject to} 
 & \sum_{j \in J} x_{ij} &\leq 1 
 & &\forall i \in I \\
 &
 &\sum_{i \in I} x_{ij} &= 1
 & &\forall j \in J \\
 &
 & x_{ij} &\in \{0,1\}
 & &\forall i \in I, \quad \forall j \in J.
\end{alignat}

余談(2)

数理最適化の式って、$\sum$を頻繁に使いますよね。数理最適化のモデリング言語は、$\sum$をどう記述できるかが使いやすさの基本になってきます。その点、PuLPは安心です。
ところで、数理最適化の入門テキストで、この問題でいうところの集合$I$や$J$を、$I:=\{1, \ldots, m\}$や$J:=\{1, \ldots, n\}$などと、むりやり1から番号づけして、例えば制約条件を

\sum_{j = 1}^{n} x_{ij} \leq 1 \quad \mbox{for} \;\; i = 1, \ldots, m

などと記述しているものがありますが、$I$や$J$の要素が何を表しているのか、数字の$2$は$I$と$J$のどっちの要素なのか、混乱してくるので、集合を正の整数の列に置き換えず、集合の要素の名前も元の名前のまま使うことをお勧めします。PuLPでは、Pythonの辞書やタプルといった言語仕様を使うことで、自然にモデリングできます。

あと、この問題の$\LaTeX$コマンドで書いた数式のように、各式の縦をそろえると、見やすくなっていいですよね。

あとあと、この問題は、0-1変数$x_{ij} \in \{0,1\}$を連続変数$0 \leq x_{ij} \leq 1$などに変更しても、問題の全ユニモジュラ性(total unimodularity)のおかげで、整数の最適解が得られるのはよく知られていますが、実際に仕事でモデル化すると、どんどん制約条件が増えていって、全ユニモジュラ性が崩れちゃうのがほとんどなんですよね。


コード

pulp_problem_2.py
# coding: UTF-8

# 線形/整数最適化問題を解くためにPuLPをインポート
import pulp
# 計算時間を計るのに time をインポート
import time



# 作業員の集合(便宜上、リストを用いる)
I = ["Aさん", "Bさん", "Cさん"]

print("作業員の集合 I = {:}".format(I))


# タスクの集合(便宜上、リストを用いる)
J = ["仕事イ", "仕事ロ", "仕事ハ"]

print("タスクの集合 J = {:}".format(J))


# 作業員 i を タスク j に割り当てたときのコストの集合(一時的なリスト)
cc = [
      [ 1,  2,  3],
      [ 4,  6,  8],
      [10, 13, 16],
     ]

# cc はリストであり、添え字が数値なので、
# 辞書 c を定義し、cc[0][0] は c["Aさん","仕事イ"] でアクセスできるようにする
c = {} # 空の辞書
for i in I:
    for j in J:
        c[i,j] = cc[I.index(i)][J.index(j)]

print("コスト c[i,j]: ")
for i in I:
    for j in J:
        print("c[{:},{:}] = {:2d},  ".format(i, j, c[i,j]), end = "")
    print("")
print("")



# 数理最適化問題(最小化)を宣言
problem = pulp.LpProblem("Problem-2", pulp.LpMinimize)
# pulp.LpMinimize : 最小化 
# pulp.LpMaximize : 最大化


# 変数集合を表す辞書
x = {} # 空の辞書
       # x[i,j] または x[(i,j)] で、(i,j) というタプルをキーにしてバリューを読み書き

# 0-1変数を宣言
for i in I:
    for j in J:
        x[i,j] = pulp.LpVariable("x({:},{:})".format(i,j), 0, 1, pulp.LpInteger)
        # 変数ラベルに '[' や ']' や '-' を入れても、なぜか '_' に変わる…?

# 内包表記も使える
#x_suffixes = [(i,j) for i in I for j in J]
#x = pulp.LpVariable.dicts("x", x_suffixes, cat = pulp.LpBinary) 
# lowBound, upBound を指定しないと、それぞれ -無限大, +無限大 になる

# pulp.LpContinuous : 連続変数
# pulp.LpInteger    : 整数変数
# pulp.LpBinary     : 0-1変数


# 目的関数を宣言
problem += pulp.lpSum(c[i,j] * x[i,j] for i in I for j in J), "TotalCost"
#problem += sum(c[i,j] * x[i,j] for i in I for j in J)
# としてもOK


# 制約条件を宣言
# 各作業員 i について、割り当ててよいタスク数は1つ以下
for i in I:
    problem += sum(x[i,j] for j in J) <= 1, "Constraint_leq_{:}".format(i)
    # 制約条件ラベルに '[' や ']' や '-' を入れても、なぜか '_' に変わる…?

# 各タスク j について、割り当てられる作業員数はちょうど1人
for j in J:
    problem += sum(x[i,j] for i in I) == 1, "Constraint_eq_{:}".format(j)


# 問題の式全部を表示
print("問題の式")
print("--------")
print(problem)
print("--------")
print("")



# 計算
# ソルバー指定
solver = pulp.solvers.PULP_CBC_CMD()
# pulp.solvers.PULP_CBC_CMD() : PuLP付属のCoin-CBC
# pulp.solvers.GUROBI_CMD()   : Gurobiをコマンドラインから起動 (.lpファイルを一時生成)
# pulp.solvers.GUROBI()       : Gurobiをライブラリーから起動 (ライブラリーの場所指定が必要)
# ほかにもいくつかのソルバーに対応
# (使用例)
#if pulp.solvers.GUROBI_CMD().available():
#    solver = pulp.solvers.GUROBI_CMD()

# 時間計測開始
time_start = time.clock()

result_status = problem.solve(solver)
# solve()の()内でソルバーを指定できる
# 何も指定しない場合は pulp.solvers.PULP_CBC_CMD() 

# 時間計測終了
time_stop = time.clock()



# (解が得られていれば)目的関数値や解を表示
print("計算結果")
print("********")
print("最適性 = {:}, 目的関数値 = {:}, 計算時間 = {:} (秒)"
      .format(pulp.LpStatus[result_status], pulp.value(problem.objective), 
              time_stop - time_start))
print("解 x[i,j]: ")
for i in I:
    for j in J:
        print("{:} = {:},  "
              .format(x[i,j].name, x[i,j].value()), end="")
    print("")
print("********")

出力

作業員の集合 I = ['Aさん', 'Bさん', 'Cさん']
タスクの集合 J = ['仕事イ', '仕事ロ', '仕事ハ']
コスト c[i,j]: 
c[Aさん,仕事イ] =  1,  c[Aさん,仕事ロ] =  2,  c[Aさん,仕事ハ] =  3,  
c[Bさん,仕事イ] =  4,  c[Bさん,仕事ロ] =  6,  c[Bさん,仕事ハ] =  8,  
c[Cさん,仕事イ] = 10,  c[Cさん,仕事ロ] = 13,  c[Cさん,仕事ハ] = 16,  

問題の式
--------
Problem-2:
MINIMIZE
1*x(Aさん,仕事イ) + 3*x(Aさん,仕事ハ) + 2*x(Aさん,仕事ロ) + 4*x(Bさん,仕事イ) + 8*x(Bさん,仕事ハ) + 6*x(Bさん,仕事ロ) + 10*x(Cさん,仕事イ) + 16*x(Cさん,仕事ハ) + 13*x(Cさん,仕事ロ) + 0
SUBJECT TO
Constraint_leq_Aさん: x(Aさん,仕事イ) + x(Aさん,仕事ハ) + x(Aさん,仕事ロ) <= 1

Constraint_leq_Bさん: x(Bさん,仕事イ) + x(Bさん,仕事ハ) + x(Bさん,仕事ロ) <= 1

Constraint_leq_Cさん: x(Cさん,仕事イ) + x(Cさん,仕事ハ) + x(Cさん,仕事ロ) <= 1

Constraint_eq_仕事イ: x(Aさん,仕事イ) + x(Bさん,仕事イ) + x(Cさん,仕事イ) = 1

Constraint_eq_仕事ロ: x(Aさん,仕事ロ) + x(Bさん,仕事ロ) + x(Cさん,仕事ロ) = 1

Constraint_eq_仕事ハ: x(Aさん,仕事ハ) + x(Bさん,仕事ハ) + x(Cさん,仕事ハ) = 1

VARIABLES
0 <= x(Aさん,仕事イ) <= 1 Integer
0 <= x(Aさん,仕事ハ) <= 1 Integer
0 <= x(Aさん,仕事ロ) <= 1 Integer
0 <= x(Bさん,仕事イ) <= 1 Integer
0 <= x(Bさん,仕事ハ) <= 1 Integer
0 <= x(Bさん,仕事ロ) <= 1 Integer
0 <= x(Cさん,仕事イ) <= 1 Integer
0 <= x(Cさん,仕事ハ) <= 1 Integer
0 <= x(Cさん,仕事ロ) <= 1 Integer

--------

計算結果
********
最適性 = Optimal, 目的関数値 = 19.0, 計算時間 = 0.029293706568961887 (秒)
解 x[i,j]: 
x(Aさん,仕事イ) = 0.0,  x(Aさん,仕事ロ) = 0.0,  x(Aさん,仕事ハ) = 1.0,  
x(Bさん,仕事イ) = 0.0,  x(Bさん,仕事ロ) = 1.0,  x(Bさん,仕事ハ) = 0.0,  
x(Cさん,仕事イ) = 1.0,  x(Cさん,仕事ロ) = 0.0,  x(Cさん,仕事ハ) = 0.0,  
********

おわりに

ボブの絵画教室よりも簡単だと思うので、ぜひPython+PuLPに触れてみてください。