線形計画法(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]だけ持っている.小麦を育てるためには肥料F1 [kg/km^2]と農薬P1 [kg/km^2]が必要である.一方,大麦を育てるためには肥料F2 [kg/km^2]と農薬P2 [kg/km^2]が必要である.それぞれの売値が小麦S1 [yen/kg],大麦S2 [yen/kg]だとすると,売上を最大化するためにはどのように農地を振り分ければいいでしょう?
という問題を線形計画法で表すと
Maximize:
S1 x1 + S2 x2 Subject to:
0 ≤ x1 + x2 ≤ L
0 ≤ F1 x1 + F2 x2 ≤ F
0 ≤ P1 x1 + P2 x2 ≤ P
x1 ≥ 0, x2 ≥ 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と答える.
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
1 =3, S
2 =2, L=5, F
1 =1, F
2 =3, F=10, P
1 =2, P
2 =1, P=9とした線形計画法を解いてみる.
Maximize:
3x1 + 2x2 Subject to:
0 ≤ x1 + x2 ≤ 5
0 ≤ x1 + 3x2 ≤ 10
0 ≤ 2x1 + x2 ≤ 9
x1 ≥ 0, x2 ≥ 0
コード
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from FuncDesigner import *
from openopt import LP
# 変数の定義
x1,x2 = oovars(2)
# 目的関数
f = 3*x1 + 2*x2
# 初期解
startPoint = {x1:0, x2:0}
# 制約条件
constraints = [x1+x2<=5, x1+3*x2<=10, 2*x1+x2<=9]
# 式の定義
p = LP(f, startPoint, constraints=constraints, goal='max')
# 解
r = p.solve('glpk')
x1_opt,x2_opt = r(x1,x2)
print(x1_opt,x2_opt)
結果
------------------------- 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
1 =4.0 x
2 =1.0のとき14であることがわかりました.
確認のためにGoogleを使ってグラフで描画してみます.制約条件が
青 ,
赤 ,
黄 .目的関数が
緑 .最大化なら目的関数はなるべく右上の方に持っていきたい.ちょうど青と黄の交点が最適解(4,1)になっていることがわかります.
今回はだいぶ簡単な例題を扱いましたが,OpenOptではもっと大規模なもの扱えるようです.またOpenOptでは線形計画法のみならず,二次計画法(QP)や非線形計画法(NLP)など様々な最適化問題が解けるようです.
Cookbook にOpenOptやFuncDesingerで扱えるサンプルがまとまっているようです.あ,っていうか
SciKit のほうを使ってみればよかったか…?