Rによる2次計画法

今回は久しぶりにRの話題。
2次計画問題を解くのに、Rの"quadprog"パッケージを使うことにした。
例題は東工大の水野先生のサイト↓
http://www.me.titech.ac.jp/~mizu_lab/text/pdf-file/LP9-QP-problem.pdf

例題では制約が等式制約である。だがquadprogで扱える制約は不等式制約である。
等式制約を不等式制約に変えるのは難しくなく、
 a_1 x_1 + a_2 x_2 = c
ならば、
 a_1 x_1 + a_2 x_2 \geq c
 -a_1 x_1 - a_2 x_2 \geq -c
という2つの制約に落とし込めばよい。

# TODO: 狭義2次計画問題を解く
###############################################################################
setwd("H:\\RW\\QPTest")

library(quadprog)
##
## Assume we want to minimize: -(8 12) %*% b + 1/2 b^T %*% D %*% b
## D = (2  1)
##     (1  4)
## under the constraints:      A^T b >= b0
## with b0 = (2 -2)^T
## and  A = ( 2  1)
##          (-2 -1)
## we can use solve.QP as follows:
##
Dmat       <- matrix(c(2,1,1,4),2,2)
dvec       <- c(8,12)
Amat       <- matrix(c(2,1,-2,-1),2,2)
bvec       <- c(2,-2)
residue    <- 34
result <- solve.QP(Dmat,dvec,Amat,bvec=bvec)

result$solution
result$value + residue
save.image("2011.05.05_QPTest2.RData")

定数項は使えないので、residueとして残しておいて、結果に加える。

結果はこちら。

> result$solution
[1] -0.1428571  2.2857143
> result$value + residue
[1] 17.85714

となり、例題にあった解析解である
 x_1 = -\frac{1}{7},x_2 = \frac{16}{7}
とその時の最適値の
 f = \frac{125}{7}
が得られていることが分かる。