Ml-tips

PuLP使ってみた

はじめに

PuLPは、線形計画（LP）や整数計画（MILP）をPythonで書いて解くためのライブラリです。数式に近い形で目的関数と制約を書くと、あとはソルバーが最適解を返してくれます。この記事では、簡単な生産計画を例に、モデルの書き方から解の読み方までを見ます。

install

pip install pulp

例題：生産計画

製品Aと製品Bを作る工場を考えます。各製品を1個作るには「作業時間」と「原料」を消費し、それぞれ使える総量が決まっています。1個あたりの利益と消費量を表にすると次の通りです。

利益（円） 作業時間（時間） 原料（kg）
製品A 1個あたり 500 2 4
製品B 1個あたり 400 3 2
使える総量 180 200

たとえば作業時間の列を見ると、製品Aは1個2時間・製品Bは1個3時間かかり、全体で180時間しか使えません。作る量をそれぞれA, B個とすると、消費する作業時間は2A + 3Bで、これが180以下でなければなりません。原料も同じように書けます。

  • 作業時間: 2A + 3B <= 180
  • 原料: 4A + 2B <= 200

利益（500A + 400B）を最大にする生産量A, Bを求めます。制約式の係数は「その製品1個あたりの資源消費量」、右辺は「使える総量」、と読むと分かりやすいです。

モデルを書いて解く

LpProblemに最大化/最小化を指定し、LpVariableで変数を作ります。目的関数と制約は+=でproblemに足していきます。lowBound=0で変数を非負（生産量がマイナスにならない）にしています。

import pulp

prob = pulp.LpProblem("production", pulp.LpMaximize)
x = pulp.LpVariable("A", lowBound=0)   # 製品Aの生産量
y = pulp.LpVariable("B", lowBound=0)   # 製品Bの生産量

prob += 500*x + 400*y        # 目的関数: 利益
prob += 2*x + 3*y <= 180     # 作業時間
prob += 4*x + 2*y <= 200     # 原料

status = prob.solve(pulp.PULP_CBC_CMD(msg=0))

probをそのままprintすると、組み立てた数理モデルを確認できます。

print(prob)

production:
MAXIMIZE
500*A + 400*B + 0
SUBJECT TO
_C1: 2 A + 3 B <= 180

_C2: 4 A + 2 B <= 200

VARIABLES
A Continuous
B Continuous

解を読む

解のステータスはLpStatusで、最適解が見つかればOptimalです。各変数の値は.value()、目的関数の値はpulp.value(prob.objective)で取り出します。

print("status :", pulp.LpStatus[prob.status])
print("A =", x.value())
print("B =", y.value())
print("profit =", pulp.value(prob.objective))

status : Optimal
A = 30.0
B = 40.0
profit = 31000.0

製品Aを30個、製品Bを40個作るのが最適で、そのときの利益は31000でした。このとき作業時間（2*30 + 3*40 = 180）も原料（4*30 + 2*40 = 200）もちょうど使い切っています。両方の資源に余りがなく、利益を増やすにはどちらかの総量を増やすしかない、と読み取れます。

図で見る

2変数なので、実行可能領域（制約を全部満たすAとBの範囲）と最適解を平面に描けます。水色の領域が制約を満たす範囲、赤の破線が利益500A + 400Bの等高線です。利益を上げる方向へ線を平行移動していくと、実行可能領域の頂点で最後に接します。そこが最適解（A=30, B=40）で、作業時間の線と原料の線が交わる角にあたります。

線形計画の最適解は、こうして必ず実行可能領域の頂点に来ます。ソルバーはこの頂点をしらみつぶしに探すのではなく、効率的にたどって最適解に到達します。

注意点

  • pulp.LpMaximize / pulp.LpMinimize で最大化・最小化を切り替えます。
  • 変数は既定で連続値です。整数にしたいときはLpVariable("A", lowBound=0, cat="Integer")、0/1ならcat="Binary"にします（別記事で扱う予定）。
  • lowBound / upBound を指定しないと変数は下限なし（-∞〜+∞）になります。非負変数はlowBound=0を忘れないようにします。
  • ステータスはOptimalのほか、制約が矛盾するとInfeasible、目的関数が発散するとUnboundedになります。解いたら必ずLpStatusを確認します。
  • solve()は引数なしでも既定のCBCで解けます。PULP_CBC_CMD(msg=0)はソルバーのログを抑制するための指定です。

まとめ

  • pulp楽しい

参考リンク