December 3, 2011

Pythonで線形計画法を解いてみる


線形計画法(Linear Programming; LP or linear optimization)は,目的関数も制約関数も線形であるときの最適化問題を解く手法です.SciPy + GLPK + CVXOPT + OpenOpt + FuncDesignerで実装されているのを使ってみました.説明はWikipediaから抜粋.

線形計画法は以下のようにある制約関数のもとである目的関数を最大化あるいは最小化するもの.

maximize or minimize (目的関数)
{\b c}^T {\b x}
subject to (制約関数)
{\mb A}{\b x} \le {\b b}
{\b x} \ge {\b 0}
xは変数,A,b,cは係数.制約関数の形が凸多面体であるので目的関数は最適化可能.

たとえば
ある農夫がL [km^2]の広さの農地を持っているとき小麦と大麦を育てて売上を最大化させたい.農夫は肥料をF [kg]と農薬をP [kg]だけ持っている.小麦を育てるためには肥料F1 [kg/km^2]と農薬P1 [kg/km^2]が必要である.一方,大麦を育てるためには肥料F2 [kg/km^2]と農薬P2 [kg/km^2]が必要である.それぞれの売値が小麦S1[yen/kg],大麦S2[yen/kg]だとすると,売上を最大化するためにはどのように農地を振り分ければいいでしょう?

という問題を線形計画法で表すと

Maximize:
S1x1 + S2x2
Subject to:
0 ≤ x1 + x2 ≤ L
0 ≤ F1x1 + F2x2 ≤ F
0 ≤ P1x1 + P2x2 ≤ P
x1 ≥ 0, x2 ≥ 0

これを行列形式で書くと

Maximize
\left[\begin{array}{cc}S_1 & S_2 \\\end{array}\right]\left[\begin{array}{c}x_1 \\x_2\end{array}\right]
Subject to

\lef[\beg{a}{c}x_1\\x_2\end{a}\right]\ge\lef[\beg{a}{c}0\\0\end{a}\right]

となる.



準備


環境はMac OS X 10.7 (Lion).要XCode 4.2GitHomebrewsetuptools


SciPy http://fonnesbeck.github.com/ScipySuperpack/

以下を実行.
git clone git://github.com/fonnesbeck/ScipySuperpack
cd ScipySuperpack
sh install_superpack.sh
Are you installing from a repository from this machine?にはyと答える.


GLPK (GNU Linear Programming Kit) http://www.gnu.org/software/glpk/

GLPK(Gnu Linear Programming Kit)は,線形計画問題(LP)や混合整数計画問題(MIP)を解くためのソルバー.GNU GPL.Homebrewを使ってインストール.
sudo brew install glpk


CVXOPT http://abel.ee.ucla.edu/cvxopt/download/index.html

CVXOPTはPythonの凸最適化(convex optimization)のためのパッケージ.sageでも利用可.GNU GPL.左のメニューのDownloadを選び,Installation from sourceのPython 2.5+をダウンロードして./srcのsetup.pyのBUILD_GLPKを0から1に書き換え,
# Set to 1 if you are installing the glpk module.
BUILD_GLPK = 1
以下を実行.
sudo setup.py install


OpenOpt + FuncDesigner http://www.openopt.org/Install

OpenOptは数値最適化フレームワーク.自前のソルバーに加え,複数のソルバーが使える.BSD.以下を実行.
sudo easy_install openopt
sudo easy_install FuncDesigner



実装


上の問題でS1=3, S2=2, L=5, F1=1, F2=3, F=10, P1=2, P2=1, P=9とした線形計画法を解いてみる.

Maximize:
3x1 + 2x2
Subject to:
0 ≤ x1 + x2 ≤ 5
0 ≤ x1 + 3x2 ≤ 10
0 ≤ 2x1 + x2 ≤ 9
x1 ≥ 0, x2 ≥ 0

コード

結果
------------------------- OpenOpt 0.36 -------------------------
solver: glpk problem: unnamed type: LP goal: max
iter objFunVal log10(maxResidual)
0 -0.000e+00 -100.00
GLPK Simplex Optimizer, v4.47
3 rows, 2 columns, 6 non-zeros
Preprocessing...
3 rows, 2 columns, 6 non-zeros
Scaling...
A: min|aij| = 1.000e+00 max|aij| = 3.000e+00 ratio = 3.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part = 3
* 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0)
* 2: obj = -1.400000000e+01 infeas = 0.000e+00 (0)
OPTIMAL SOLUTION FOUND
1 1.400e+01 -100.00
istop: 1000 (optimal)
Solver: Time Elapsed = 0.01 CPU Time Elapsed = 0.006718
objFunValue: 14 (feasible, MaxResidual = 0)
(4.0, 1.0)
最適解はx1=4.0 x2=1.0のとき14であることがわかりました.

確認のためにGoogleを使ってグラフで描画してみます.制約条件が.目的関数が.最大化なら目的関数はなるべく右上の方に持っていきたい.ちょうど青と黄の交点が最適解(4,1)になっていることがわかります.


今回はだいぶ簡単な例題を扱いましたが,OpenOptではもっと大規模なもの扱えるようです.またOpenOptでは線形計画法のみならず,二次計画法(QP)や非線形計画法(NLP)など様々な最適化問題が解けるようです.CookbookにOpenOptやFuncDesingerで扱えるサンプルがまとまっているようです.あ,っていうかSciKitのほうを使ってみればよかったか…?