具有等式和不等式约束的 R 优化

R optimization with equality and inequality constraints

我正在尝试寻找函数的局部最小值,并且参数的总和是固定的。例如,

Fx = 10 - 5x1 + 2x2 - x3

条件如下,

x1 + x2 + x3 = 15

(x1,x2,x3) >= 0

其中 x1、x2 和 x3 的总和具有已知值,并且它们都大于零。在 R 中,它看起来像这样,

Fx = function(x) {10 - (5*x[1] + 2*x[2] + x[3])}
opt = optim(c(1,1,1), Fx, method = "L-BFGS-B", lower=c(0,0,0), upper=c(15,15,15))

我还尝试使用 constrOptim 的不等式来强制固定总和。我仍然认为这可能是一个可行的解决方法,但我无法让它发挥作用。这是实际问题的简化示例,但我们将不胜感激。

在这种情况下 optim 显然不会起作用,因为你有等式约束。 constrOptim 出于同样的原因也不起作用(我尝试将等式转换为两个不等式,即大于和小于 15,但这不适用于 constrOptim)。

但是,有一个包专门解决这种问题,那就是Rsolnp

您的使用方式如下:

#specify your function
opt_func <- function(x) {
  10 - 5*x[1] + 2 * x[2] - x[3]
}

#specify the equality function. The number 15 (to which the function is equal)
#is specified as an additional argument
equal <- function(x) {
  x[1] + x[2] + x[3] 
}

#the optimiser - minimises by default
solnp(c(5,5,5), #starting values (random - obviously need to be positive and sum to 15)
      opt_func, #function to optimise
      eqfun=equal, #equality function 
      eqB=15,   #the equality constraint
      LB=c(0,0,0), #lower bound for parameters i.e. greater than zero
      UB=c(100,100,100)) #upper bound for parameters (I just chose 100 randomly)

输出:

> solnp(c(5,5,5),
+       opt_func,
+       eqfun=equal,
+       eqB=15,
+       LB=c(0,0,0),
+       UB=c(100,100,100))

Iter: 1 fn: -65.0000     Pars:  14.99999993134  0.00000002235  0.00000004632
Iter: 2 fn: -65.0000     Pars:  14.999999973563  0.000000005745  0.000000020692
solnp--> Completed in 2 iterations
$pars
[1] 1.500000e+01 5.745236e-09 2.069192e-08

$convergence
[1] 0

$values
[1] -10 -65 -65

$lagrange
     [,1]
[1,]   -5

$hessian
          [,1]      [,2]      [,3]
[1,] 121313076 121313076 121313076
[2,] 121313076 121313076 121313076
[3,] 121313076 121313076 121313076

$ineqx0
NULL

$nfuneval
[1] 126

$outer.iter
[1] 2

$elapsed
Time difference of 0.1770101 secs

$vscale
[1] 6.5e+01 1.0e-08 1.0e+00 1.0e+00 1.0e+00

因此得到的最优值为:

$pars
[1] 1.500000e+01 5.745236e-09 2.069192e-08

表示第一个参数为15,其余为零和零。这确实是函数中的全局最小值,因为 x2 正在添加到函数中,并且 5 * x1 对结果的影响比 x3 大得多(负面)。 15,0,0的选择是根据约束对函数的解和全局最小值。

功能很好用!

这实际上是一个线性规划问题,因此自然的方法是使用线性规划求解器,例如 lpSolve 程序包。您需要提供一个 objective 函数和一个约束矩阵,求解器将完成剩下的工作:

library(lpSolve)
mod <- lp("min", c(-5, 2, -1), matrix(c(1, 1, 1), nrow=1), "=", 15)

然后可以得到最优解和objective值(加上常数项10,不提供给求解器):

mod$solution
# [1] 15  0  0
mod$objval + 10
# [1] -65

线性规划求解器应该比一般的非线性优化求解器快得多,并且返回精确的最优解应该没有问题(而不是可能会出现舍入误差的附近点)。

重写

fx = function(x) {10 - (5*x[1] + 2*x[2] + (15 - x[3] - x[2]))}

并添加额外的约束

x[1] + x[2] <= 15

您可以稍微不同地定义问题,而不是使用外部库:

x1 + x2 + x3 = 15 等同于 1/15 (x1+x2+x3) = 1,我将其表示为 y1+y2+y3=1

所以x1等于y1的15倍。

现在,您可以将 y 定义为 softmax 函数的输出 - 这可确保 y 为非负数且总和为 1。因此 y = softmax(z)

要获得完整的转换,我们可以这样写:

softmax.15 = function(z) {
  y = exp(z)/sum(exp(z))
  x = 15*y
  return(x)
}

现在您可以运行使用新变量对问题进行无约束优化:

Fx = function(z) {
  x = softmax.15(z)
  10 - (5*x[1] + 2*x[2] + x[3])
}
opt = optim(c(1,1,1), Fx, method = "L-BFGS-B")

这会输出最佳的 z,因此请记住对输出再次进行转换以获得最佳的 x

softmax.15(opt$par)
# [1] 1.500000e+01 2.445345e-08 1.147868e-08