読者です 読者をやめる 読者になる 読者になる

J's blog

趣味で統計•データ解析をしています

RでLPを解く

R

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