Rで非線形制約条件付非線形最適化を行う方法 by Rsolnp −２

R

↓導入の記事はこちら
http://d.hatena.ne.jp/teramonagi/20091217/1261048574

簡単な例題を通して動作を確かめんとす。

library(Rsolnp)
#(x,y)=(1,2)で最小値をとるような凸関数
objectiveFunc <- function(x_)
{
  return(sum((x_-c(1,2))^2))
}
#決定変数のスタート値
x0 <- c(100,100)
#最適化
solution <- solnp(x0,fun = objectiveFunc)
#結果出力
print(solution)

これを走らせると最後のprint分から

> print(solution)
$pars
[1] 1 2
$convergence
[1] 0
$values
[1] 1.940500e+04 4.909325e-11 1.149875e-14
$lagrange
[1] 0
$hessian
          [,1]      [,2]
[1,] 1.5050760 0.4999742
[2,] 0.4999742 1.4949240
$ineqx0
NULL
$nfuneval
[1] 46
$outer.iter
[1] 2
$elapsed
Time difference of 0.01641297 secs

な感じで詳細な結果が帰ってくる。結果自体は

> str(solution)
List of 9
 $ pars       : num [1:2] 1 2
 $ convergence: num 0
 $ values     : num [1:3] 1.94e+04 4.91e-11 1.15e-14
 $ lagrange   : num 0
 $ hessian    : num [1:2, 1:2] 1.51 0.5 0.5 1.49
 $ ineqx0     : NULL
 $ nfuneval   : num 46
 $ outer.iter : num 2
 $ elapsed    :Class 'difftime'  atomic [1:1] 0.016
  .. ..- attr(*, "units")= chr "secs"

というようにList型で返ってくるのであとは適当に必要なものを取得できる。
それぞれの要素の意味は・・・（マニュアルより）

  • $ pars : 最適パラメーター（最適化でいう決定変数）
  • $ convergence: 収束したか否かを０（収束）、１（非収束）で返却
  • $ values : 各最適化ステップでの目的関数値
  • $ lagrange : ラグランジュ乗数
  • $ hessian : 最適解でのへシアン
  • $ ineqx0 : 不等式制約を等式制約になおすためのスラック変数の値
  • $ nfuneval : 目的関数の評価回数
  • $ outer.iter : 繰り返し回数の上限値
  • $ elapsed : 計算時間

となる。

これに非線形な等式制約をいれるには以下のように書く。

library(Rsolnp)
#(x,y)=(1,2)で最小値をとるような凸関数
objectiveFunc <- function(x_)
{
  return(sum((x_-c(1,2))^2))
}
#適当な非線形制約
equalityConstraint <- function(x_)
{
    return(x_[1]^4)
}

#決定係数のスタート値
x0 <- c(100,100)
#最適化：制約条件 x^4(eqfunで指定） = 16(eqBで指定）
solution <- solnp(x0,fun = objectiveFunc,eqfun = equalityConstraint,eqB = c(16))
#結果出力
print(solution)

これを実行すると(2､2)で最適値をとることがわかり、たしかに最適解になってる。
不等式制約もほぼ同様に書くことができて

library(Rsolnp)
#(x,y)=(1,2)で最小値をとるような凸関数
objectiveFunc <- function(x_)
{
  return(sum((x_-c(1,2))^2))
}
#適当な非線形制約(楕円）
inequalityConstraint <- function(x_)
{
    return(x_%*%matrix(1,2,2)%*%x)
}
#決定係数のスタート値
x0 <- c(100,100)
#最適化（楕円制約,最適解は結局(1,2)の点）
solution <- solnp(x0,fun = objectiveFunc,ineqfun = inequalityConstraint,ineqUB = c(16),ineqLB=c(0))
#結果出力
print(solution)

とすることができる。このケースの場合、結局最適解は同じで(1,2)の点になる。
もし、等式・不等式ともに複数の制約条件をいれたいときは制約条件関数の返り値とeqB,ineqLU,ineqLB引数をベクトルにすればOK。