使用 R 包"nloptr"最大化非线性约束问题

Question

我需要使用R包'nloptr'最大化一些问题的目标函数。我尝试了基本规则"最大化f(x(<=>最小化-f(x("，但它不起作用。我不确定使用它或有其他方法出了什么问题。

下面是一个完整的例子。当前的解决方案只是具有最小目标值的初始向量。但是，我应该得到最大化目标函数的解决方案。有人能帮我拿吗？谢谢！

library(nloptr)
X = log(rbind(c(1.350, 8.100),   
c(465.000,     423.000),
c(36.330  ,    119.500),
c(27.660   ,   115.000), 
c(1.040     ,  5.500),
c(11700.000, 50.000),  
c(2547.000  ,  4603.000),
c(187.100    , 419.000), 
c(521.000  ,   655.000), 
c(10.000    ,  115.000), 
c(3.300    ,   25.600),  
c(529.000   ,  680.000), 
c(207.000  ,   406.000), 
c(62.000    ,  1320.000),
c(6654.000 ,   5712.000),
c(9400.000  ,  70.000),  
c(6.800      , 179.000), 
c(35.000   ,   56.000),  
c(0.120  ,     1.000),   
c(0.023   ,    0.400),   
c(2.500  ,     12.100),  
c(55.500  ,    175.000), 
c(100.000  ,   157.000), 
c(52.160   ,   440.000), 
c(87000.000 ,  154.500), 
c(0.280  ,     1.900),   
c(0.122  ,     3.000),   
c(192.000 ,    180.000)))
n = nrow(X)
q = 0.5
x0 = cbind(8,4)
x01 = x0[1]
x02 = x0[2]
x1 = X[,1]
x2 = X[,2]
pInit = c(0.1614860, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 
0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 
0.0000000, 0.0000000, 0.0000000, 0.7124934, 0.0000000, 0.0000000,
0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 
0.1260206, 0.0000000, 0.0000000, 0.0000000)
eval_f0 = function(p) {
obj0 = mean((n * p ) ^ q)
grad0 = rbind(q * ((n * p) ^ (q - 1))/((mean((n * p ) ^ q))^2))
return(list("objective" = obj0, "gradient" = grad0))
}
eval_g_eq0 = function(p) {
sum0 = sum(x1 * p) - x01
sum1 =  sum(x2 * p) - x02
sum2 = sum(p) - 1
constr0 = rbind(sum0, sum1, sum2)
grad0 = rbind(x1, x2, rep(1,n))
return(list("constraints" = constr0, "jacobian" = grad0))
}
local_opts <- list( "algorithm" = "NLOPT_LD_AUGLAG",
"xtol_rel" = 1.0e-7 )
opts <- list( "algorithm" = "NLOPT_LD_AUGLAG",
"xtol_rel" = 1.0e-7,
"maxeval" = 10000,
"local_opts" = local_opts )
res1 = nloptr(x0 = c(pInit), 
eval_f = eval_f0, 
lb = c(rep(0, n)), 
ub = c(rep(Inf, n)), 
eval_g_eq = eval_g_eq0, 
opts = opts )
weight = res1$solution
fval0 = res1$objective
print(list(fval0, weight))

Answer 1

请注意，起点pInit处的梯度(和雅可比(不是有限的，这使得任何基于梯度的求解器都很难执行此任务。我将使用一个不同的起点，稍微远离边界。

无论如何，使用阿拉巴马包中的拉格朗日解算器似乎更容易找到最大值。根据您上面的定义，最高可达x1 = X[,1]; x2 = X[,2]，可能的解决方案如下：

f1 <- function(p) mean((n * p ) ^ q)
heq1 <- function(p) 
c(sum(x1 * p) - x01, sum(x2 * p) - x02, sum(p) - 1)

为了简单起见，我们让求解器计算梯度和雅可比矩阵。若要找到最大值，请将解算器应用于目标函数的负值。

sol <- alabama::auglag(rep(0.1, 28), fn=function(p) -f1(p), heq=heq1)
cat("The maximum value is:", -sol$value, 'n')
## The maximum value is: 0.7085338

满足相等条件，见

heq1(sol$par)
## [1] -1.685957e-08  3.721533e-08 -2.935964e-08

找到的解决方案是

sol$par
##  [1] 0.012186842 0.006640286 0.006706268 0.006418224 0.014501609 0.405618998  
##  [7] 0.003531462 0.005458189 0.005582029 0.005158098 0.008072278 0.005510394  
## [13] 0.005653117 0.002935642 0.003861549 0.123009564 0.004021866 0.009866779  
## [19] 0.024385229 0.027101557 0.011436010 0.006184886 0.007473135 0.004162962  
## [25] 0.245429952 0.019978294 0.010919515 0.008195238

我很想知道这对你来说是否是一个合理的解决方案！我检查了几个起点，结果总是一样的。