[R言語] Rのお部屋::あーるのおへや[別館]

[R] [R言語] [2次計画法] [2次計画問題] [非線形計画法] [数理計画法] [最適化]

Rのお部屋で自分用のメモとして書いておいた記事（２次計画法（２次計画問題）のメモ）を、少しわかりやすくまとめてみることにしました．

自分用のメモのため、見に来てくれた方にはわかりにくいものとなっていましたので．

Rで二次計画法を行うには「quadprog」というライブラリを使います．

quadprogの説明（別ウインドウでPDFを開きます）にノーテーションを合わせます

min(-d^Tb + 1/2 b^TDb)
with constraints A^Tb ≧ b₀

ここで

b ： 求めたい変数のベクトル n×1の行列（要素n個のベクトル）
D ： n×nの対角行列
d ： n×1の行列（要素n個のベクトル）
A ： 制約条件の左辺の係数
b₀ ： 右辺の定数項

1/2 b^TDb

-d^Tb

制約条件は≧としかなっていませんが、等式制約と不等式制約両方扱えます（あたりまえ）．
ここでquadprogの制約条件は行列Aと制約条件の右辺のb₀で設定しますが
等式制約 → 不等式制約
の順でAとb₀を設定しなくてはなりません．これは後で説明します．

関数solve.QPのサンプルの例を実行してみましょう．最初のコメント部分に以下のように書いてあります．

##
## Assume we want to minimize: -(0 5 0) %*% b + 1/2 b^T b
## under the constraints: A^T b >= b0
## with b0 = (-8,2,0)^T
## and
##     (-4 2 0)
## A = (-3 1 -2)
##     ( 0 0 1)
## we can use solve.QP as follows:
##

これはまず求めたい変数の次元が３で、目的関数が

-(0 5 0) %*% b + 1/2 b^T b

制約条件はbを（b1、b2、b3）とすると、 b0 = (-8,2,0)^Tが右辺、Aが左辺ですから

-4×b1 ＋ -3×b2 ＋       ≧ -8
 2×b1 ＋  1×b2 ＋       ≧  2
          -2×b2 ＋ 1×b3 ≧  0

ここで2次の項を構成する対角行列のDがありませんが

1  0  0
0  1  0
0  0  1

Rのソースを記述します．DmatがD行列を表しています．Dは単位行列ですので以下でDmatを設定できます．

Dmat <- matrix(0,3,3)
diag(Dmat) <- 1
Dmat
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

次に目的関数の1次の項のdvecを設定します．ノーテーションではn×1の行列と書きましたが必ず1列なのでRでいうところのただのベクトルで設定します．

dvec <- c(0,5,0)

Amatは上記制約条件のとおりですが、以下の並びになるように設定します．制約条件の不等式の係数の並びとは行と列が逆になっていますので注意が必要です．

     (-4 2 0)
 A = (-3 1 -2)
     ( 0 0 1)

以下でAmatを設定します．

Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
Amat
     [,1] [,2] [,3]
[1,]   -4    2    0
[2,]   -3    1   -2
[3,]    0    0    1

bvecは制約条件の右辺ですので

bvec <- c(-8,2,0)

実行して

> solve.QP(Dmat,dvec,Amat,bvec=bvec)
$solution
[1] 0.4761905 1.0476190 2.0952381

$value
[1] -2.380952

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.0000000 0.2380952 2.0952381

$iact
[1] 3 2

等式制約がある場合
Rによる二次計画法 例