はじめに¶
この記事では，
Python製のモデリングツールPyomoを使って標準的な線形計画問題を解く方法
を紹介します．
線形計画問題の数学的な話，解法アルゴリズム，ソルバー(gurobiとかcplexとか)の詳細な仕様については書きません．
Pythonで線形計画問題を解く方法は数多くあるのですが，その多さとPython界隈の開発のスピードのおかげで(せいで)，参入障壁が高くなっているように感じます．
そもそもPythonで数値計算を行うことを目的としている人々にとっては，写経でもまず動くコードが欲しいのだと思います．
(自分もそうです．CSガチ勢でもなく，あくまでツールとして使いこなすことを目的にPythonを使っています，)
そこで，今回は自分が混乱したところを中心に，線形計画問題を解く即戦力コードを紹介したいと思います．
(*) 参考記事，サイトは数は最後にまとめます．
ソルバーとモデラー¶
最初の混乱ポイントががソルバーとモデラーの話です．
まず，ソルバーは汎用的な数理問題を解くツールです，
代表的なものに，Gurobi,cplex, ipoptなどがあります．
これらのソルバーは厳格な入力フォーマットを要求し，仕様も初心者に難しいので直接使うのは困難です．
さらに，大規模な問題でも高速に解くことを目的にしているため，PythonではなくCやC++で書かれています．
ここで，登場するのがモデラーです．モデリングツールとも呼ばれます．
これはインタプリンタ言語であるPythonとソルバーのつなぎ役を担います．
具体的には，ソルバーに渡すための入力フォーマットを作成し，ソルバーが解いてくれたらPythonの変数として受け取るなどの仕事をします．
これらは当然Pythonで書かれており，直感的なPython記法で記述が可能で，初心者にも一応優しい作りとなっています．
つまり，
モデラーを使ってPythonで問題(モデル)を記述し，モデラーがソルバーに投げて，解をモデラーが受け取る
というようなシステムのイメージです．
(*)文献によっては両者の違いを区別していない場合もあります，(明確な線引きは個別のツール依存で難しいため)
ソルバーの種類¶
正直，自分もあまり詳しくないのですが有名どころは以下です．
Scipy
- 純Pythonのソルバー
- 書きやすいが，大規模問題には無力(100変数ぐらいでストレスになる計算時間)
CBC
- C#で書かれたソルバー
- COIN-OR Projectの一貫
- 混合整数計画問題に対する分枝切除法(branch and cut method)を実装している
glpk
- 純C言語の無料ソルバー
- 実装アルゴリズム(LP向け)
- primal and dual simplex methods(主問題ベースの単体法・双対問題ベースの単体法)
- primal-dual interior-point method(主双対内点法)
- branch-and-cut method(分枝切除法)
- 選択可能オプション GLPK/Using GLPSOL – Wikibooks, open books for an open world
- GLPK – GNU Project – Free Software Foundation (FSF)
ipopt
- C++で主双対内点法を実装
- 滑らか？な解がでる印象
Gurobi
- 商用ソルバー，アカデミックライセンスあり
- 選択可能なアルゴリズム Gurobi Method
- 単体法
- 内点法
- バリアー法
cplex
- IBM製のソルバー，世界最速との呼び声高い
- アカデミックライセンスで無料で利用可能(購入だとすごい高い)
- 記述言語は不明．アセンブリレベルでチューニングされているとか
- IBM Knowledge Center – Overview of LP optimizers
- ILOGの最適化ソリューション 最初の一歩
- マルチスレッドも簡単に記述できる
モデラーの種類¶
線形計画問題を扱うPythonのモデラーも複数あります．代表的な3つと僕個人の印象をまとめます，
Scipy
- Python製の最適化ツール(モデラー・ソルバー)
- アルゴリズム部分(前述の区分で言うソルバー部分)もPythonなのでひたすら遅い
- 単体法，内点法を選択可能
- 書きやすいため，規模の小さい問題には適している
PuLP
- Python製のモデラー
-
- COIN-OR Projectの一貫
- 呼出可能ソルバー
- CBC
- glpk
- 2017年夏頃より開発が止まっているらしい
- 日本語の解説本も存在するが，英語含めドキュメントは貧弱
- myPuLPというPuLP自体のモデラーも存在しややこしい
Pyomo
- Python製のモデラー
- 呼出可能ソルバー
- Gurobi
- glpk
- ipopt
- cplex
- Googleコミュニティが活発
- ドキュメントは貧弱だがなんとか頑張れる
今回はPyomoを使った方法を紹介します
線形計画問題の基本形¶
以下に線形計画問題の基本形[P]を示します．
どんな線形計画問題もこの形に落とし込めるので，この解き方をマスターすればあとは微調整でなんとかなります．
[P]
$$
\max_{x} \qquad b^{\mathsf{T}} x
\\
\mbox{s.t.} \qquad A x \leq c
\\
\qquad \qquad x \geq 0
$$
ここで，
- $x \in \mathbb{R}^{N}$ : 決定変数ベクトル
- $b \in \mathbb{R}^{N}$ : 係数ベクトル
- $A \in \mathbb{R}^{M \times N}$ : 係数行列
- $c \in \mathbb{R}^{M}$ : 定数ベクトル
です．
今回は，より具体的に，$N=2, M=2$とし，
$$
b =
\begin{bmatrix}
2 \\
3
\end{bmatrix}
\\
A =
\begin{bmatrix}
4 & 3 \\
1 & 2
\end{bmatrix}
\\
c =
\begin{bmatrix}
120 \\
60
\end{bmatrix}
$$
として計算します．
(*) ちなみにこの問題は，和歌山大学のレジュメの2ページにあったチョコレート工場の生産問題です，
環境構築¶
Pyomoをinstallします．いずれかの方法でOKです．(*)Pyomo GitHub
pip install pyomoconda install -c https://conda.anaconda.org/conda-forge pyomo
Pyomoは，各種ソルバーを自在に呼び出すことが可能ですが，ソルバーを使うためにはローカル環境にソルバーがinstallされている必要があります．
Gurobiやipoptは導入が少し手間なので，今回はgplkを使います．
glpkは以下のいずれかで導入可能です．
pip install glpkconda install -c conda-forge glpk
実装¶
長くなりましたが，実装していきます．
import pyomo.environ as pyo
from pyomo.opt import SolverFactory
import numpy as np
# 使うソルバーの確認
SolverFactory("glpk").available() == True
#変数・パラメータの設定
#係数行列の次元
N=2
M=2
#係数ベクトル，係数行列の作成
b = np.array([2,3], dtype="float")
A = np.array([[4,3],[1,2]], dtype="float")
c = np.array([120,60], dtype="float")
そして謎の初期化関数もセットで用意します，(理由はよくわかりません)
#決定変数の設定
x = np.zeros(shape=(N), dtype="float")
#決定変数初期化用関数---------------------
xx = {i:x[i] for i in range(len(x))}
def init_x(model,i):
return xx[i]
# -----------------------------------
#index list
N_index_list = list(range(N))
M_index_list = list(range(M))
ここら辺はこう書くものと割り切っています．ここで．
index_listが活躍します．(おそらくPyomoは一般的な配列のindexを扱えないのだと思います，そのため単純な
index_listによって要素指定をします．)
#不等式制約
def const(model, j):
tmp = sum(A[j,i]*model.x[i] for i in N_index_list)
return tmp <= c[j]
#目的関数
def func(model):
return sum(b[i] * model.x[i] for i in N_index_list)
#モデル
model = pyo.AbstractModel("LP-sample")
#決定変数をモデルに組み込み
model.x = pyo.Var(N_index_list, domain=pyo.NonNegativeReals, initialize=init_x)
#目的関数をモデルに組み込む
model.value = pyo.Objective(rule=func, sense=pyo.maximize)
#制約をモデルに組み込み
model.const = pyo.Constraint(M_index_list, rule=const)
デフォルトでもいい感じに解いてくれます．
#ソルバーの選択とオプション
solver_name = "glpk"
opt = SolverFactory(solver_name)
instance = model.create_instance()
instance.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) #双対変数取得用に必要
results = opt.solve(instance, tee=True)#tee=Trueでログ出力
instance.load(results)
results.write()
#解の概要を見る
instance.display()
#最適変数値だけ見る
print(instance.x[0].value)
print(instance.x[1].value)
#最適目的関数値を見る
instance.value.expr()
#双対変数の一覧を見る
instance.dual.display()
#最適双対変数値だけ見る
print(instance.dual[instance.const[1]])
print(instance.dual[instance.const[0]])
おわりに¶
おおむね即戦力になるコードの紹介ができたと思います．
ただ，Pyomoもドキュメントが少し怪しかったり，美しい唯一の書き方がなかったりするので自分なりに試行錯誤は必要です．
違うソルバーの導入方法や使い方，オプションの設定方法についてはまた記事にしたいと思います，
参考サイト¶
- carldlaird/pyomo: The main repository for the Pyomo optimization modeling software.
- GitHubページ
- Python + Pyomoによる(非線形)数値最適化 – Easy to type
- 日本語で一番わかりやすい記事でした
- Pyomo Meets Fantasy Football | The Data Barista
- 英語だけどわかりやすい
- Accessing Pyomo variable values and objective function value – hselab.org
追記
* 2019.03.22 : hikima氏の指摘により主問題[P]の誤植修正．コード内のコメントにおける，最適値と最適変数値の混同を修正．
