PuLPによるモデル作成方法

Pythonで線形最適化を行うには、PuLPライブラリーを使用します。 PuLPは、COIN-ORプロジェクトで開発されたものです。

参考:COIN-OR

PuLPのインストール

コマンドプロンプトで、pip install -U pulpでできます(macOSではpip3 install -U pulp)。

モデル作成の手順

  1. モデル(LpProblem)の作成
  2. 変数(LpVariable)の作成
  3. 目的関数の設定
  4. 制約条件の追加

※ 1と2の順番は逆でも構いません。3と4の順番も逆でも構いません。また、3と4はなくても構いません。

問題とモデルの違い

問題には、2種類の意味があります。

  • 解決したいと思っている課題。目標と現実の違い。
  • 答えを想定している問い。試験の問題など。最短路問題のように名前がついている問題の多くは、この問題。

「数理的アプローチによる問題解決」コースにおける問題は、前者になります。 一方、モデルは、問題をPCなどで扱えるように表現したものです。 最短路問題のように名前がついている問題は、 研究者にとって「解決したいと思っている目標が、数式で表現できる」ため、モデルのように詳細に記述できます。 PuLPの開発者は研究者でもあったので、モデルのことを問題と考えたのだと思います。しかし、LpProblemは、一般的にはモデルと呼ばれるべきものです。 問題とモデルについて、変えることを考えると、わかりやすいかもしれません。

  • モデルを変えること:変数や目的関数や制約条件をとらえなおすことです。解決したいと思っていることを変えることではありません。
  • 問題を変えること:目標や現実のとらえ方を変えることです。リフレーミングなどがあります(リフレーミング - wikipedia)。

モデルの作成

モデルには、最小化最大化の2種類あります。PuLPでは、モデル作成時に指定します。

  • 最小化のモデル作成:m = LpProblem(sense=LpMinimize)
  • 最大化のモデル作成:m = LpProblem(sense=LpMaximize)

デフォルトはsense=LpMinimizeなので、引数を省略すると最小化になります。

  • LpProblemの第一引数は、モデル名です。省略すると'NoName'になります。
    • モデル名は、nameプロパティーで確認できます。
  • 最小化と最大化のどちらであるかは、senseプロパティーで確認できます。
  • Jupyter Notebookで、モデルを評価すると、下記のように、LPフォーマットで確認できます。各行は、モデル名、モデルの種類、制約条件、変数宣言になります。
NoName:
MAXIMIZE
None
VARIABLES

参考:宮代 隆平 の web ページ(整数計画法メモ)- LP ファイルについて

変数作成

変数とは、意思決定の対象です。具体例を上げてみましょう。

  • 自宅から学校への最短路を求めたい場合の変数:各経路を通るかどうか
  • 生産額を最大化したい場合の変数:製品ごとの生産量
  • 輸送費用を最小化したい場合の変数:経路ごとの輸送量

PuLPの変数の種類は、以下の3種類あります。 線形最適化では、連続変数だけ使いますLpVariableは、デフォルトで連続変数になります。

  • 連続変数:実数値を表す変数です。
  • 0-1変数:0または1のみを表す変数です。
  • 整数変数:整数値を表す変数です。

LpVariableのオプション

  • 第1引数:変数名(必須)
  • lowBound:変数の下限。デフォルトはNoneで下限なしです。
  • upBound:変数の上限。デフォルトはNoneで上限なしです。
  • cat:変数の種類。デフォルトはpulp.LpContinuousで連続変数です。

LpVariableの作成例

  • LpVariable('x')自由変数を作成します。自由変数とは、-∞から∞までの値を取ることができる変数です。
  • LpVariable('y', lowBound=0)非負変数を作成します。非負変数とは、0から∞までの値を取ることができる変数です。
  • LpVariable('y', lowBound=0, upBound=1):上下限のある変数を作成します。範囲は、lowBoundの値からupBoundの値までです。

LpVariableの注意事項

  • LpVariableは、必ず変数名を指定する必要があります。また、全ての変数は、異なる名前になっていないといけません。

実際に同じ変数名('x')を作って実行してみましょう。

m = LpProblem()
x = LpVariable('x', lowBound=0)
y = LpVariable('x', lowBound=0)  # 変数名が同じ!
m += x + y
m.solve()

下記のようにエラーになりますが、エラーメッセージを見ても原因がわかりません。同じ変数名を作成しないように注意しましょう。

PulpSolverError                           Traceback (most recent call last)
中略
PulpSolverError: Pulp: Error while executing .../cbc

目的関数の設定

目的関数とは、どうしたいかを決める関数です。具体例を上げてみましょう。

  • 自宅から学校への最短路を求めたい場合の目的関数:経路の距離
  • 生産額を最大化したい場合の目的関数:総生産額
  • 輸送費用を最小化したい場合の目的関数:総輸送費用

m += のようにして、モデルに目的関数を設定します。

式の書き方

  • 変数や定数を使った足し算や引き算が使えます。掛け算は2 * xのような定数と変数の掛け算だけ使えます。1 / xx * xx * yのような式は使えません。
  • xを変数のリストとします。リストの合計はlpSum(x)のように記述します。「x[0] + x[1] + + x[-1]」と同じですが、効率的に計算してくれます。
  • aを定数のリストとします。axの内積(a[0] * x[0] + a[1] * x[1] + + a[-1] * x[-1])は、lpDot(a, x)のように記述します。

m + = について

PuLPで最も間違えやすいと思われるのが、この式でしょう。 「モデルに式を足す」表現で「モデルに目的関数を設定する」ことを意味します。 追加ではなく設定です。下記を実行すると、目的関数は「式1+式2」ではなく、「式2」だけになります。ただし、目的関数を2度以上設定すると警告が出ます。

m += 式1
m += 式2  # 警告が出ます

目的関数の設定は、m.setObjective(式)でもできます。

lpSumについて

lpSumの代わりにsumを使っても、同様に合計できます。しかし、sumは非効率的なのでlpSumを使うようにしましょう。

lpDotについて

lpDot(係数リスト, 変数リスト)のように使えます。 内部では、zipを使っています。係数リストと変数リストの長さは同じになるようにしましょう。 長さが異なっていても、はみ出た分が無視されてエラーにはなりません。バグの元になるので気をつけましょう。

制約条件の追加

制約条件とは、守らなければならないことを指定します。具体例を上げてみましょう。

  • 自宅から学校への最短路を求めたい場合の制約条件:選んだ経路が、自宅から学校まで接続していること。
  • 生産額を最大化したい場合の制約条件:使用する原料が、原料の在庫以下であること。
  • 倉庫から工場への輸送費用を最小化したい場合の制約条件:倉庫から出荷する総量が供給量以下であること。工場への入荷する総量が需要量以上であること。

制約条件の書き方 制約条件は、以下の3種類の書き方ができます。式1と式2の片方は定数でも構いません。

m += 式1 >= 式2
m += 式1 == 式2
m += 式1 <= 式2

下記の2種類の書き方は実行エラーになります。「=」を含んだ形式にしてください。

m += 式1 > 式2
m += 式1 < 式2

下記の書き方はエラーになりません。ただし、「式 != 式2」は、TrueまたはFalseになります。結果として、下記の式は、目的関数に定数を設定します。このような制約条件は、数理モデルにはありませんが、PuLPで書いてもエラーにならないので注意しましょう。

m += 式1 != 式2

3種類の制約条件のいずれの書き方でも、意図通りの結果が得られています。

制約条件の書き方は非常にシンプルなルールですが、実に多彩なことがらを表現できます。 PyQでは、色々な書き方を通して学ぶ予定です。

PuLPでのソルバーの実行

solve()でソルバーの実行します。 内部で以下の処理が行われています。

  • モデルをMPSフォーマットとよばれる形式でファイルに出力します。
  • そのファイルを入力として、デフォルトのソルバーCBCを実行します。
  • CBCを実行した結果のファイルを読み込み、モデルのステータスと、変数の結果を設定します。
  • CBCの入力ファイルと出力ファイルを削除します。

参考:MPS(format)

ソルバーの指定方法

単にPuLPをインストールしただけですと、使えるソルバーはCBCだけです。 別途、GLPK(無料)、Gurobi(有料)、CPLEX(有料)などをインストールすると、以下のようにソルバーを変更できます。

from pulp import (LpProblem, PulpSolverError,
                  GLPK_CMD, GUROBI_CMD, CPLEX_CMD)
solver1 = GLPK_CMD()
solver2 = GUROBI_CMD()
solver3 = CPLEX_CMD()
for solver in [solver1, solver2, solver3]:
    try:
        m = LpProblem()
        m.solve(solver)
    except PulpSolverError:
        print(solver.path, 'はインストールされていません')

PuLPで使えるソルバーは、pulp.pulpTestAll()で確認できます。

参考

ステータスの確認

実行した結果のステータスコードは、m.statusで数値で確認できます。文字列で確認するときは、pulp.LpStatus[m.status]とします。

ステータスの全種類は、pulp.LpStatusで確認でき、以下の通りです。

{-3: 'Undefined',
 -2: 'Unbounded',
 -1: 'Infeasible',
 0: 'Not Solved',
 1: 'Optimal'}

結果を取得する前に、m.status == 1であることを確認しましょう。1Optimalなので、「最適解が得られた」ことになります。

変数や式や目的関数の値の確認

値は、value(オブジェクト)で確認できます。

print(value(x))  # 変数の値
print(value(x + 1)) # 式の値
print(value(m.objective))  # 目的関数の値

数理的アプローチによる問題解決