Zenn
🔢

Python の数理最適化ソルバーを試してみた

2021/04/25に公開
Python
ソルバー
pulp
CVXPY
tech

python の数理最適化ソルバー、pulpとcxvpyを使って試してみたので、メモとして残します。

コードはこちら。

今回のサンプル問題は下記サイトから参照しています。

http://www.nct9.ne.jp/m_hiroi/light/pulp01.html

工場 (x, y) から商品を店 (a, b, c) に配送。
供給量、需要量、輸送コストが下表で与えられているとき、総輸送コストが最小となる条件を求める。

表 : 輸送コスト

店 a 店 b 店 c 供給量
工場x 10 6 16 8
工場y 8 8 4 16
需要量 12 4 8

上記表をDataFrameで使用します。

pulp

最適化問題の名称の設定と、最小化(最大化)問題の設定を行います。

import pulp

prob = pulp.LpProblem('sample', # 問題の名称
                      pulp.LpMinimize) # 最小化。なおデフォルトは最小化

最大化問題の場合は pulp.LpMaximize を引数に入れます。

店名(店a、店b、店c)をリスト化します。

item_list = ["a", "b", "c"]

変数を設定します。

# 工場Xの供給量
fact_x = [pulp.LpVariable(f'fact_x_{i}', 
                    lowBound = 0, # 下限値。デフォルト値はNone
                    upBound=None, # 上限値。デフォルト値はNone
                    cat=pulp.LpContinuous
                         # 変数の種類。デフォルトは"Continuous"(連続値)
                         # 他には pulp.LpInteger 整数、pulp.LpBinary バイナリ値
                    ) for i in shop_list]

工場xから店aに送付する商品の量をxa、というように1個ずつ変数を設定してもいいのですが、面倒なので内包表記を使います。

一部省略して抜粋すると、

[pulp.LpVariable(
    f"x_{i}" # 変数名。
    ) for i in shop_list]

となります。

変数を指定する pulp.LpVariable 内の変数名に、フォーマット済み文字列を使い、for文でリスト shop_list 内の["a", "b", "c"]を入れて一度にx_a、x_b、x_c という変数を作成します。

工場yの供給量も同様です。

# 工場yの供給量
fact_y = [pulp.LpVariable(f'fact_y_{i}', 
                    lowBound = 0, # 下限値。デフォルト値はNone
                    upBound=None, # 上限値。デフォルト値はNone
                    cat=pulp.LpContinuous
                         # 変数の種類。デフォルトは"Continuous"(連続値)
                         # 他には pulp.LpInteger 整数、pulp.LpBinary バイナリ値
                    ) for i in shop_list]

目的変数を設定します。

prob += pulp.lpSum(fact_x * df.iloc[0, :3]) + pulp.lpSum(fact_y * df.iloc[1, :3])

制約条件を設定します。

需要量

for i in range(3):
    prob += fact_x[i] + fact_y[i] >= df.iloc[2, i]

供給量

prob += pulp.lpSum(fact_x) <= df["供給量"][0]
prob += pulp.lpSum(fact_y) <= df["供給量"][1]

実行

status = prob.solve()
print("Status", pulp.LpStatus[status])

Status Optimal

実行結果のステータス 意味
'Undefined' 存在する可能性があるが実行可能な解は発見されなかった
'Unbounded' 非有界。値を無限にとれるので最適解が存在しない
'Infeasible' 実行不能。制約条件を満たす解が存在しない
'Not Solved' 未解決
'Optimal' 最適解が存在する

https://stackoverflow.com/questions/24167958/what-does-pulp-lpstatus-undefined-actually-mean

結果をDataFrameで表示

result = pd.DataFrame(
    [
        [fact_x[i].value() for i in range(3)],
        [fact_y[i].value() for i in range(3)],
    ],
    columns=["店a", "店b", "店c"],
    index=["工場x", "工場y"]
)

最小化した目的関数

prob.objective.value()

160.0

cxvpy

cxvpyはpulpと違い、変数の要素数を指定することができます。

import cvxpy as cp

fact_x = cp.Variable(
                    len(shop_list), # 要素数。2次元の場合は(x, y)と設定する
                    integer=False, # 整数の場合はTrueとする
                    boolean=False, # Binary値の場合はTrueとする
                    pos=True # 正の数
                    )

cp.Variable(3) とすれば要素数3の変数を指定できるし、cp.Variable((3, 4))とすれば 3, 4 のように多次元の要素数をもつ変数を設定できます。

fact_y = cp.Variable(
                    len(shop_list), # 要素数。2次元の場合は(x, y)と設定する
                    integer=False, # 整数の場合はTrueとする
                     boolean=False, # Binary値の場合はTrueとする
                    pos=True # 正の数
                    )

目的変数の設定と最小化か最大化の設定

exp = cp.sum(fact_x * df.iloc[:1, :3].squeeze() + fact_y * df.iloc[1:2, :3].squeeze())
obj = cp.Minimize(exp) # 最大化のときはcp.Maximize()

需要量

cxvpxの制約条件設定は、

1回目 変数 = [(イコール) 左式 >= 右式]
2回目以降 変数 += [(プラスイコール) 左式 >= 右式]

とする必要があるようです。

今回のように、for文を使っているため、1回目～3回目がまとまってしまうときは少し面倒ですが、下記のように変数=[]とするとよさそうです。

const = []

for i in range(3):
    const += [fact_x[i] + fact_y[i] >= df.iloc[2, i]]

供給量

const += [cp.sum(fact_x) <= df["供給量"][0]]
const += [cp.sum(fact_y) <= df["供給量"][1]]

問題(目的変数・制約条件)の設定と実行

prob = cp.Problem(obj, const)
status = prob.solve(verbose=True)

変数の表示

少数点以下の細かい数字が表示されるので、小数点以下15桁くらいまで表示するよう設定しました。

pd.options.display.precision = 15

result = pd.DataFrame(
    [
        [fact_x[i].value for i in range(3)],
        [fact_y[i].value for i in range(3)],
    ],
    columns=["店a", "店b", "店c"],
    index=["工場x", "工場y"]
)

最小化した目的変数の出力

prob.value

pulp と cxvpy の違い

pulp
https://coin-or.github.io/pulp/index.html

  • 線形計画のみを求めることができる。
    cxvpy
    https://www.cvxpy.org/index.html
  • 線形計画だけでなく、非線形計画も求めることができる。

数理最適化を学べるサイト・ドキュメント

以上になります、最後までお読みいただきありがとうございました。

Discussion