数学やプログラミングの備忘録

数理最適化, Python, C++をメインに紹介するブログ。

MENU

pulpを使いこなすための備忘録

pulpは、pythonで数理最適化のモデルを記述するためのモジュール です。モデル記述後、そのまま指定した最適化ソフトウェアで解くことができます。デフォルトでは、最適化ソフトウェアとして同梱されているCBCが呼び出されます。pulpCBCは、COIN-ORにより開発されています。

pulppythonのモジュールということだけあって、記述がとても楽です解きたい問題を速攻で記述できます

インストール方法も簡単です(condaも同様です)。

pip install pulp

この記事では、そんな便利なpulpを使いこなすための備忘録です。

参考:PuLPによるモデル作成方法 — Pythonオンライン学習サービス PyQ(パイキュー)ドキュメント

目次

  1. pulpが記述できる数理最適化問題モデル
  2. 最小の記述例
  3. 使い方
  4. 高速化・便利ワザ・注意点

1. pulpが記述できる数理最適化モデル

pulpが記述できる数理最適化の問題は、

  • 目的関数: 線形関数
  • 制約条件: 線形等式制約、線形不等式制約
  • 変数: 連続変数、整数変数、01変数

を含むものとなります。2次関数はどうやら記述できないみたいです。pulpで記述できるものを表でまとめると、こんな感じ。

モデル 連続変数 整数変数 01変数
線形計画問題
整数計画問題
0-1整数計画問題
混合整数計画問題

2. 最小の記述例

ナップサック問題(0-1整数計画問題)を例にpulpで記述してみます。


\begin{array}{ll}
\max_{x} & 3 x_1 + 4 x_2 + 5 x_3,\\
\text{s.t.} & 2 x_1 + 2 x_2 + x_3  \le 3,\\
&x_{i}\in \{0,1\}\ (i = 1,2, 3).\\
\end{array}

pythonコード:

import pulp

# (1) 問題を定義
problem = pulp.LpProblem( "test", pulp.LpMaximize )

# (2) 変数を定義
x = [ pulp.LpVariable( 'x_{}'.format( i ), cat=pulp.LpBinary ) for i in range(3) ]

# (3) 係数をセット
c = [ 3, 4, 5 ]
a = [ 2, 2, 1 ]
b = 3

# (4) 目的関数をセット
assert len(c) == len(x) == len(a)
problem += pulp.lpDot( c, x )

# (5) 制約条件をセット
problem += pulp.lpDot( a, x ) <= b

# (6) 解く
result = problem.solve()

# 結果を出力
print('objective value: {}'.format(pulp.value(problem.objective)))
print('solution')
for i in range(len(x)):
    print('x_{} = {}'.format( i, pulp.value(x[i]) ))

実行結果:

objective value: 9.0
solution
x_0 = 0.0
x_1 = 1.0
x_2 = 1.0

解きたいナップサック問題の係数と、コード中の c, a, b を目で追えば直感的に理解できるかなと思います。

例のように規模の小さい問題であれば、係数をダイレクトに書いちゃうのもありですね。

#problem += pulp.lpDot( c, x )
#problem += pulp.lpDot( a, x ) <= b
#この部分は↓書き換え可能です
problem += 3 * x[0] + 4 * x[1] + 5 * x[2]
problem += 2 * x[0] + 2 * x[1] + 1 * x[2] <= 3

3. 使い方

モデルの記述方法&求解は、ずばり6ステップです:

  1. pulp.LpProblem(): 問題を定義
  2. pulp.LpVariable(): 変数を定義
  3. 係数をセット
  4. +=: 目的関数をセット
  5. +=: 制約条件をセット
  6. solve: 求解スタート

問題を定義

problem = pulp.LpProblem( "test", pulp.LpMaximize )
  • 第1引数: 問題の名前
  • 第2引数: 問題のタイプを指定します
    • 最大化 → pulp.LpMaximize
    • 最小化 → pulp.LpMinimize

変数を定義

x = [ pulp.LpVariable( 'x_{}'.format( i ), cat=pulp.LpBinary ) for i in range(3) ]

上記の例ではリスト内包表記で、変数を定義しています。dict型でもOKですので、必要に応じて使い分けることができます。

引数は他にもあって、変数の上限と下限の制約を設定することができます。設定しなければ、デフォルト値が使われます。

  • 第1引数(必須): 変数の名前
  • その他よく使う引数:
    • lowBound: 下限を設定(デフォルトは  -\infty)
    • upBound: 上限を設定(デフォルトは  \infty)
    • cat: 変数のタイプを設定
      • pulp.LpContinuous: 連続変数(デフォルト)
      • pulp.LpBinary: 0-1変数, 上下限は設定しなくてOK
      • pulp.LpInteger: 整数変数

例えば、 0 \le y \le 100, y \in \mathbb{Z} (整数変数) を定義したければ、

y = pulp.LpVariable( 'y', lowBound=0, upBound=100, cat=pulp.LpInteger )

と書けばよいです。

係数をセット

係数あるいは係数の計算に必要な数値をリスト等の配列に入れておきます

c = [ 3, 4, 5 ]
a = [ 2, 2, 1 ]
b = 3

目的関数をセット

目的関数とは、最大化(最小化)したい関数のことです。セットの仕方として、ここではオススメ順に 4つ紹介します

lpDot を使う

problem += pulp.lpDot( c, x )

lpDot は高速で、記述が楽なので、一番オススメの方法です。problem += c[0] * x[0] + c[1] * x[1] + c[2] * x[2] という意味になります。

lpSum を使う

problem += pulp.lpSum( c[i] * x[i] for i in range(3) )

今回の例ではlpDotを使えばいいんですが、lpSumを使ってリスト内包表記みたに書くこともできます。lpSumも高速です。ifif .. else .. を使うこともできます。例えば、

problem += pulp.lpSum( 3 * x[i] if i%2 == 0 else 5 * x[i] for i in range(3))

sum を使う

problem += sum( c[i] * x[i] for i in range(3) )

lpDotlpSumに比べ極端に遅いです。数千以上の変数を扱う問題ならば、sum は避けるべきです。

直接書く

problem += 3 * x[0] + 4 * x[1] + 5 * x[2]

規模が大きいと書くの大変ですし、コードが汚くなりますね(笑)

制約条件をセット

セットできる制約条件の種類は3つです:

  •  \sum a_i x_i \le bproblem += pulp.lpDot( a, x ) <= b
  •  \sum a_i x_i = bproblem += pulp.lpDot( a, x ) == b
  •  \sum a_i x_i \ge bproblem += pulp.lpDot( a, x ) >= b

セット方法も目的関数同様にいくつかあります(オススメ順):

  • lpDot を使う
  • lpSum を使う
  • sum を使う
  • 直接書く

求解スタート

result = problem.solve()

solve() 関数で求解スタートです。環境によっては、ログが出たり出なかったりします。返り値は次のようになっています。

{-3: 'Undefined',   # 未定義
 -2: 'Unbounded',   # 非有界
 -1: 'Infeasible',  # 実行不可能
 0: 'Not Solved',   # 解けなかった
 1: 'Optimal'}      # 最適解を発見

4. 高速化・便利ワザ・注意点

sum は使うな

既に何度も書いてきましたが、目的関数や制約条件をセットするときは、lpDot or lpSum を使った方がいいです。速度が段違いです。

lpDot を使うときは長さに気をつけろ

lpDot は、内部で zip が使われるらしいので、長さが短い方に合わせられます。意図していない動作にならならいように、assert でチェックしておくといいでしょう。

assert len(c) == len(x) == len(a)
problem += pulp.lpDot( c, x )
problem += pulp.lpDot( a, x ) <= b

記述内容は表示できる

初めてpulpを使うなら、意図したモデルになっているか確認したくなりますよね。例えば、上の例なら、

print(problem)

で、記述したモデルを表示できます。

出力:

test:
MAXIMIZE
3*x_0 + 4*x_1 + 5*x_2 + 0
SUBJECT TO
_C1: 2 x_0 + 2 x_1 + x_2 <= 3

VARIABLES
0 <= x_0 <= 1 Integer
0 <= x_1 <= 1 Integer
0 <= x_2 <= 1 Integer

求解オプション(並列化, 時間制限, and ログ出力)

# オプション無し
#result = problem.solve()
# オプション有り
result = problem.solve(pulp.PULP_CBC_CMD( msg=0, threads=4, timeLimit=100 ))
  • msg: 1 →ログ有り, 0→ログ無し
  • threads: 2 以上でCBCの処理を並列計算で行います(基本的に使えるだけ使った方が性能が上がります)
  • timeLimit: 時間制限です(ただし, 実時間でなくCPU時間)

関数solve() の返り値を信じるな

以下のように、pulp.LpStatus[result] を出力すれば、結果を確認することができますが、

result = problem.solve(pulp.PULP_CBC_CMD( msg=0, threads=4, timeLimit=100 ))
print('result: {}'.format(pulp.LpStatus[result])) 

これはいつも正しい結果では無いことを確認しています。おそらくバグです。このバグの有無は、pulpのバージョンにより異なります。バグを確認している状況は、

  • MacOS/Linux
  • pulp 2.3
  • threads=8
  • timeLimit=100
  • CBCのログを見る限り、タイムオーバーだが結果はOptimal

感覚的には、「並列計算+タイムオーバー」がバグの発生条件です。

pulp.value(x[i]) の結果も信用するな

ほとんどの場合、信用していいです。しかし、ごく稀に、とても難しい問題の場合で、30分たっても制約条件を満たす解を発見できない場合があります。このような問題に対して、30分の時間制限で関数solve() を呼び出し、pulp.value(x[i]) と解を参照すると、、

その解は制約条件を満たさないけれど、参照自体はできてしまいます。お気をつけください。

LPファイルとして書き出す

problem.writeLP("test.lp")

中身:

\* test *\
Maximize
OBJ: 3 x_0 + 4 x_1 + 5 x_2
Subject To
_C1: 2 x_0 + 2 x_1 + x_2 <= 3
Binaries
x_0
x_1
x_2
End

書き出しはできますが、読み込みはできません(正直不便)。

おわり