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

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

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

maximize or minimize （目的関数）

subject to （制約関数）

xは変数，A,b,cは係数．制約関数の形が凸多面体であるので目的関数は最適化可能．

たとえば

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

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

Maximize:

S₁x₁ + S₂x₂

Subject to:

0 ≤ x₁ + x₂ ≤ L
0 ≤ F₁x₁ + F₂x₂ ≤ F
0 ≤ P₁x₁ + P₂x₂ ≤ P
x₁ ≥ 0, x₂ ≥ 0

これを行列形式で書くと

Maximize

Subject to

となる．

準備

環境はMac OS X 10.7 (Lion)．要XCode 4.2，Git，Homebrew，setuptools．


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

実装

上の問題でS₁=3, S₂=2, L=5, F₁=1, F₂=3, F=10, P₁=2, P₂=1, P=9とした線形計画法を解いてみる．

Maximize:

3x₁ + 2x₂

Subject to:

0 ≤ x₁ + x₂ ≤ 5
0 ≤ x₁ + 3x₂ ≤ 10
0 ≤ 2x₁ + x₂ ≤ 9
x₁ ≥ 0, x₂ ≥ 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)

最適解はx₁=4.0 x₂=1.0のとき14であることがわかりました．

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

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