Rによる非線形計画問題

R

統計ソフトRのRsolnpパッケージを使って、2次式制約のもと2次式の最適化を行った。
制約が2次式なので、2次計画法を解くパッケージである"quadprog"は利用できない。

今回の例題はこちら。
minimize
(x-2)^2+(y-2)^2
subject to
x^2+y^2=1

Lagrangeの未定乗数法なり、円が接する状態を幾何学的に解くなりすると最適解が求まり、最適解は
(2 \sqrt2 - 1)^2 = 3.343146
であり、このとき
x=y=\frac{\sqrt2}{2}
である。

setwd("H:\\workspace3.5_R\\RsolnpTest")

library(Rsolnp)

#equality constraint
eq.value <- c(1)
equalityConstraint <- function(x){
	A <- matrix(c(1,0,0,1),2,2)
	return(t(x) %*% A %*% x)
}

#objective function
objectiveFunction <- function(x){
	A <- matrix(c(2,0,0,2),2,2)
	b <- matrix(c(-4,-4),2,1)
	c <- matrix(c(8),1,1)
	return((1/2)* t(x) %*% A %*% x + t(x) %*% b + c)
}

#execute optimization
x <- c(1,1)
solution <- solnp(x,objectiveFunction,eqfun = equalityConstraint,
		eqB = eq.value)
solution$pars
solution$values

#correct solution ((2*sqrt(2)-1)^2)
correct.solution <- (2*sqrt(2)-1)^2
correct.solution

#-- save RData
save.image("2011.05.11_RsolnpTest.RData")

結果はこちら。

Iter: 1 fn: 3.1250	 Pars:  0.75000 0.75000
Iter: 2 fn: 3.3368	 Pars:  0.70833 0.70833
Iter: 3 fn: 3.3431	 Pars:  0.70711 0.70711
Iter: 4 fn: 3.3431	 Pars:  0.70711 0.70711
Iter: 5 fn: 3.3431	 Pars:  0.70711 0.70711

solution$pars
solution$values
[1] 0.7071068 0.7071068
[1] 2.000000 3.125000 3.336806 3.343140 3.343146 3.343146

となっている。
出力結果は最適化計算の途中計算の推移を並べて出力する。

ちなみにHelpに書いてあるのだが、最適化計算を行うのに逐次2次計画法(SQP)を使っている。

これは使えそうなパッケージの予感。