Pythonで線形最適化問題を解く

概要

本記事では線形最適化問題に対する Python での実装例を二つ解説します。

はじめに

線形最適化問題とは、数理モデルの目的関数、制約条件が線形かつ変数が連続値をとる問題です。

実装例１

問題

以下の問題を考えてみましょう。

ある工場では製品pとqを製造しています。製品pとqを製造するには原料m、nが必要で、以下のことが分かっています。
・製品pを1kg製造するには原料mが1kg、原料nが2kg必要である。
・製品qを1kg製造するには原料mが3kg、原料nが1kg必要である。
・原料mの在庫は30kg、原料nの在庫は40kgである。
単位量あたりの利得は、製品pは1万円、製品qは2万円である。
この条件で利得を最大にするためには製品pとqをそれぞれ何kg製造すればよいでしょうか。

この問題を定式化すると

変数: x,y \geq 0 \\ 目的関数: x+2y \\ 制約条件: x+3y \leq 30 \\ \hspace{50pt} 2x+y \leq 40 \\

のようになります。

実装

上記の数理モデルをコードに落とし込むと

from pulp import *

#最大最適化問題
problem = LpProblem("LP", LpMaximize)

x = LpVariable("x", cat="Continues")
y = LpVariable("y", cat="Continues")

#条件式
problem +=  x >= 0
problem +=  y >= 0
problem += x+3*y <= 30
problem += 2*x+y <= 40
problem += x+2*y


#数理モデルを解く。
status = problem.solve()
print(LpStatus[status])

print("x = ", x.value(),"y = ", y.value(),"obj", problem.objective.value())

のようになります。コードを順に解説していきます。

解説

problem = LpProblem("LP", LpMaximize)

数理モデルを定義しています。第一引数には数理モデルの名前を任意で定義します。ここでは線形計画問題（Linear Programming）にちなんで LP としています。第二引数では LpMaximize で最大化問題を解くことを指定しています。

x = LpVariable("x", cat="Continues")
y = LpVariable("y", cat="Continues")

変数を定義しています。第一引数の名前は任意で指定します。ここでは分かりやすさのために x、y としています。第 2 引数の cat = "Continues"はその変数が連続値をとること、つまり実数値をとることを指定しています。

#条件式
problem += x >= 0
problem += y >= 0
problem += x+3*y <= 30
problem += 2*x+y <= 40

#目的関数
problem += x+2*y

制約式と目的関数を追加しています。 += の右辺に <= が付いているのが制約式で付いていないのが目的関数です。

#数理モデルを解く。
status = problem.solve()
print(LpStatus[status])

数理モデルを解いています。解けたか解けなかったかの情報が status に入ります。

print("x = ", x.value(),"y = ", y.value(),"obj = ", problem.objective.value())

最適化した変数 x、y とそれを目的関数に代入したときの値を出力しています。プログラムを実行すると x = 18.0、y = 4.0、obj = 26.0（18.0+2*4.0 = 26.0）を出力するはずです。

実装例２

より実務に近い実装例をご紹介します。本格的に実務で数理最適化をする場合、変数や制約条件がかなり多くなってきます。その場合、変数や制約条件をファイル等に保持しておいてプログラムから呼び出す処理にするとコード行数の削減とプログラムの汎用性に貢献します。

問題

以下、3 つの CSV ファイルがあるとします

requires.csv

p,m,require
p1,m1,2
p1,m2,0
p1,m3,1
p2,m1,1
p2,m2,3
p2,m3,2
p3,m1,0
p3,m2,1
p3,m3,2
p4,m1,2
p4,m2,1
p4,m3,3

これは製品ごとの製造するのに必要な材料とその個数のデータを保持しているファイルです。p 列は製品、m 列は材料、requirere 列は個数です。

stock.csv

m,stock
m1,35
m2,22
m3,27

これは材料の在庫数の情報を保持しています。m 列は材料、stock 列が在庫数です。

gain.csv

p,gain
p1,3
p2,4
p3,4
p4,5

これは製品ごとの利得を保持しています。ｐ列が製品名、gain 列が利得です。

以下、問題です。

利得を最大にするには製品p1、p2、p3、p4をそれぞれ何kg製造すればよいか。

この問題を定式化すると

変数: x_{p1},x_{p2},x_{p3},x_{p4} \\ 目的関数: 3x_{p1}+4x_{p2}+4x_{p3}+5x_{p4} \\ 制約条件: x_{p1},x_{p2},x_{p3},x_{p4} \geq 0 \\ \hspace{50pt}　2x_{p1}+3x_{p2}+2x_{p4} \leq 35 \\ \hspace{60pt} 2x_{p2}+2x_{p3}+2x_{p4} \leq 22 \\ \hspace{85pt} x_{p1}+2x_{p2}+2x_{p3}+3x_{p4} \leq 27 \\

のようになります。

実装

上記の数理モデルをコードに落とし込むと

import pandas as pd
import pulp

# データの取得
require_df = pd.read_csv("./requires.csv")
gain_df = pd.read_csv("./gains.csv")
stock_df = pd.read_csv("./stocks.csv")

# リストの定義
P = gain_df["p"].to_list()
M = stock_df["m"].to_list()

# 定数の定義
stock = {row.m:row.stock for row in stock_df.itertuples()}
gain = {row.p:row.gain for row in gain_df.itertuples()}
require = {(row.p,row.m):row.require for row in require_df.itertuples()}

# 数理最適化モデルを定義
problem = pulp.LpProblem("LP2",pulp.LpMaximize)

# 変数の定義。
x = pulp.LpVariable.dicts("x", P, cat = "Continuous")

# 制約式の定義
for p in P:
    problem += x[p] >= 0
for m in M:
    problem += pulp.lpSum([require[p,m] *x[p] for p in P]) <= stock[m]

# 目的関数の定義
problem += pulp.lpSum([gain[p] * x[p] for p in P])

# 求解
status = problem.solve()
print("status", pulp.LpStatus[status])
# 計算結果の表示
for p in P:
    print(p, x[p].value())

のようになります。コードを順に解説していきます。

解説

# データの取得
require_df = pd.read_csv("./requires.csv")
gain_df = pd.read_csv("./gains.csv")
stock_df = pd.read_csv("./stocks.csv")

CSV からデータを取得しています。 CSV や Excel 形式でデータが保持されている場合は pandas で用いるのが良いかと思います。

# リストの定義
P = gain_df["p"].to_list()
M = stock_df["m"].to_list()

製品名と材料名のリストを定義しています。

stock = {row.m:row.stock for row in stock_df.itertuples()}
gain = {row.p:row.gain for row in gain_df.itertuples()}
require = {(row.p,row.m):row.require for row in require_df.itertuples()}

定数を定義しています。 CSV からタプル型でデータを抽出して辞書型で格納しています。

# 数理最適化モデルを定義
problem = pulp.LpProblem("LP2",pulp.LpMaximize)

数理モデルを定義しています。

# 変数の定義。
x = pulp.LpVariable.dicts("x", P, cat="Continuous")

変数を辞書型で定義しています。 dicts()は第一引数に指定した変数名と第二引数で指定したリストを組み合わせて辞書型で変数の一覧を生成します。x を print すると

{'p1': x_p1, 'p2': x_p2, 'p3': x_p3, 'p4': x_p4}

のようになります。

# 制約式の定義
for p in P:
    problem += x[p] >= 0
for m in M:
    problem += pulp.lpSum([require[p,m] *x[p] for p in P]) <=  stock[m]

# 目的関数の定義
problem += pulp.lpSum([gain[p] * x[p] for p in P])

制約式、目的関数を定義しています。lpSum()は引数に指定したリストの総和を生成します。

# 求解
status = problem.solve()
print("status", pulp.LpStatus[status])
# 計算結果の表示
for p in P:
    print(p, x[p].value())

数理モデルを解いています。プログラムを実行すると p1 = 17.5、p3 = 4.75、p2,p4 = 0.0 を出力するはずです。

参考文献

備考

Hakky ではエンジニアを募集中です！まずは話してみたいなどでも構いませんので、ぜひお気軽に採用ページからお問い合わせくださいませ。