RでLPを解く - J's blog

RでもLP程度なら「不等式制約付き最適化問題」であっても
数値計算で求められることがわかりました。
用いる関数は、optim()のラッパー関数であるconstrOptim()を使います。

例えば、以下のような(最大化)LPを考えます。

目的関数:
　 $3x_1 + 2x_2$
制約条件:
　 $3x_1 + x_2 \le 9$
　 $2.5x_1 + 2x_2 \le 12.5$
　 $x_1 + 2x_2 \le 8$
　 $x_1, x_2 \ge 0$

(解答は, $(x_1, x_2)=(2, 3)$ のとき, 最適値12)

わかりやすさのためにまず、constrOptim()の引数のために目的関数f、制約条件uiとciを定義しておきます。fに関しては、optim()と同様に扱えます。

> # 目的関数
> f <- function(par) {
+   return(3*par[1] + 2*par[2])
+ }

uiとciに関しては、「 $ui$ %*% $par \ge ci$ 」を満たすように構成します。
parは推定するパラメータです。

> # 制約条件その1
> (ui <- matrix(c(-3, -1, -2.5, -2, -1, -2, 1, 0, 0, 1), nc=2, byrow=TRUE))
     [,1] [,2]
[1,] -3.0   -1
[2,] -2.5   -2
[3,] -1.0   -2
[4,]  1.0    0
[5,]  0.0    1
> # 制約条件その2
> (ci <- c(-9, -12.5, -8, 0, 0))
[1]  -9.0 -12.5  -8.0   0.0   0.0

さて、目的関数と制約条件の定義が済んだので、実行に入ります。

引数thetaには推定パラメータの初期値を与え、
またこの関数はデフォルトで目的関数の最小化になっているので、control=list(fnscale=-1)と加えることで、最大化にしています。
引数gradはfのグラディエント関数で、これにNULLを与えるとNelder-Mead法、関数を与えると準ニュートン法による計算となります。

> # 実行（ Nelder-Mead法 ）
> # constrOptim(theta=初期値, f=目的関数, grad=fのグラディエント関数, ui, ci, control=optim()へ渡すパラメータ)
> constrOptim(c(1, 1), f, grad=NULL, ui=ui, ci=ci, control=list(fnscale=-1))
$par
[1] 2 3

$value
[1] 12

$counts
function gradient 
     538       NA 

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 5

$barrier.value
[1] 0.002829033

結果は、 $(x_1, x_2) = (2, 3)$ が最適解であり、最適値12が得られました。

目的関数fの引数はベクトルでなければならない点と、
初期値が0だとうまくいかない点(理由はよくわからない)には注意が必要です。

参考: 線形不等式制約付きの最適化関数 constrOptim - RjpWiki