はじめに
ここでは、Pythonで最適化を始めたばかりの人が、理解を深める一助になることを目的としています。
組合せ最適化に関しては、「組合せ最適化を使おう 」をご覧ください。
モデル、変数、目的関数、制約条件という要素があることを踏まえて進めます。
対象者
- Python で最適化を学びたい人
- Excelソルバーなどで、線形最適化を実行したことがある人
- 無料でソルバーを使いたい人
- 典型的な最適化問題の雛形を見たい人
対象問題
生産最適化問題
原料に限りがある状態で製品を生産して、利益を最大化したい
詳細は、Google検索を見てください。
今回は、元となる入力データは、下記のような表で用意されているものとします。
,原料1,原料2,原料3,利益
製品1,1,4,3,5
製品2,2,4,1,4
在庫,40,80,50
Pythonコード
import pandas as pd
from pulp import *
from ortoolpy import addvars
入力表 = pd.read_csv('data.csv', encoding='cp932', index_col=0) # 諸元の読込
変数表,在庫 = 入力表.iloc[:-1,:].copy(),入力表.iloc[-1,:-1] # 2つに分ける
モデル = LpProblem(sense=LpMaximize) # 線形最適化モデル
変数表['変数'] = addvars(len(変数表)) # 生産量を表す変数を表に追加
モデル += lpDot(変数表.利益,変数表.変数) # 総利益を表す目的関数
for 原料 in 変数表.columns[:-2]: # 各製品ごとの処理
モデル += lpDot(変数表[原料],変数表.変数) <= 在庫[原料] # 原料の使用量が在庫以下
モデル.solve() # ソルバで解を求める
変数表['結果'] = 変数表.変数.apply(value) # 変数の値(結果)を表に追加
print('目的関数',value(モデル.objective)) # 目的関数の表示
print(変数表) # 変数表の表示
目的関数 95.0
原料1 原料2 原料3 利益 変数 結果
製品1 1 4 3 5.0 v0001 15.0
製品2 2 4 1 4.0 v0002 5.0
解説
パッケージ指定
import pandas as pd
from pulp import *
from ortoolpy import addvars
利用するパッケージを import します。モデルの変数は、変数表の列として表します。
- pandas: 変数表を通してモデルを扱うために使用します。
- Anaconda利用の場合 pandas はそのまま使えます。そうでない場合「pip install pandas」でインストールできます。
- pulp: 数理最適化モデルを作ることができます。
- ortoolpy: 変数生成関数(addvars)を使用します。
- このパッケージを使わなくても、pulp.LpVariable でもできます。
- この2つは「pip install pulp ortoolpy」でインストールできます。
入力表の読込
入力表 = pd.read_csv('data.csv', encoding='cp932', index_col=0) # 諸元の読込
from io import StringIO
入力表 = pd.read_table(StringIO("""
原料1 原料2 原料3 利益
製品1 1 4 3 5
製品2 2 4 1 4
在庫 40 80 50
"""),index_col=0)
必要な諸元をファイルから読込み、入力表とします。pandas.read_csvを使うと、簡単に読込むことができます。
入力表の分割
変数表,在庫 = 入力表.iloc[:-1,:].copy(),入力表.iloc[-1,:-1] # 2つに分ける
上図の青い部分を「変数表」に、緑の部分を「在庫」に代入します。入力表.iloc でスライスを使って、表から自由に取出すことができます。
スライスの詳細は、Google検索を見てください。
モデルの作成
モデル = LpProblem(sense=LpMaximize) # 最適化モデル
pulp.LpProblem を使って最適化モデル1を作成します。sense=LpMaximize を指定することによって、目的関数を最大化したいことを示します。何も指定しないと最小化になります。
変数の列の追加
変数表['変数'] = addvars(len(変数表)) # 生産量を表す変数を表に追加
pandas では、「変数表['変数'] = …」で「変数」という名前の新しい列を追加できます。
addvarsでは変数の配列を作成できます。変数表の1行が「1変数とその属性」に対応します。
目的関数の追加
モデル += lpDot(変数表.利益,変数表.変数) # 総利益を表す目的関数
pulp では、「モデル += …」で目的関数を設定できます2。
lpDot(c,x) は内積を意味し、$c^T x$になります。
制約条件の追加
for 原料 in 変数表.columns[:-2]: # 各製品ごとの処理
モデル += lpDot(変数表[原料],変数表.変数) <= 在庫[原料] # 原料の使用量が在庫以下
pulp では、「モデル += … <= 値」で制約条件を追加できます。他には、「>=」や「==」を使えますが、「!=」は使えません。
ソルバ実行
モデル.solve() # ソルバで解を求める
作成されたモデルは、ソルバと呼ばれる別のプログラムで解くことができます。
pulpでは、いろいろなソルバを指定できます。何も指定しないと無料のCBCというソルバを使って求解します。CBCは、pulpインストール時に一緒にインストールされています。
結果の追加
変数表['結果'] = 変数表.変数.apply(value) # 変数の値(結果)を表に追加
ソルバで解くと、value(XXX) で XXX の値を取出すことができます。「変数表.変数.apply(value)」とすることにより変数の列全体にvalueを適用できます。
「変数表['結果'] = …」により新しく結果の列が作られます。
結果の表示
print('目的関数',value(モデル.objective)) # 目的関数の表示
print(変数表) # 変数表の表示
目的関数 95.0
原料1 原料2 原料3 利益 変数 結果
製品1 1 4 3 5.0 v0001 15.0
製品2 2 4 1 4.0 v0002 5.0
製品1 を 15、製品2 を 5 作ると、総利益が最大の 95 になることがわかります。
その他
- 環境構築は、Anacondaが便利です。
- Jupyter Notebook を使うと試行錯誤が捗ります。
- 式の合計は、sumではなく、lpSumを使いましょう。lpSumは$\mathcal{O}(n)$ですが、sumは$\mathcal{O}(n^2)$です。
勉強するには
- 書籍を読んだりwebで調べる。
- 勉強会や入門セミナーに参加する。
- 詳しい人に聞く。
- 仕事で使うことにして自分を追い込む。
Let's enjoy Mathematical Optimization!