我有以下优化问题:

Maximisise combined benefit functions

F(x) = 2*10^9 + 160*x - 7*x*ln(x) 
F(y) = 3*10^9 + 170*y - 7.5*y*ln(y)

subject to 

x + y <= 3*10^9
x,y >= 0

我建立了一个拉格朗日

L(x,y,λ) = F(x) + F(y) + λ(3*10^9-(x+y))

然后,在用图形计算器求解之前,我对它进行了代数简化

y = 3*10^9 - x  
e^(19/15)*x^(14/15) + x - 3*10^9 = 0
y^(15/14) + e^(19/14)*y - e^(19/14)*3*10^9 = 0

x = 1.61 * 10^9
y = 1.39 * 10^9

现在,我希望能够在R中实现这个解决方案,因为实际的函数参数可能会改变.我试着改编了我在网上找到的几种不同的解决方案,但似乎都没有奏效.

如果我理解正确的话,我需要一个非线性解算器.因此,我在Rsolnp中设置了这个问题(灵感来自this answer):

library(Rsolnp)

opt_func_log <- function(x) {
  a <- x[1] 
  b <- x[2] 
  
  ben_func <-  2e9 + 160*a - 7*a*log(a) + 3e9 + 170*b - 7.5*b*log(b)
  
  -ben_func #invert to find minimum
}

equal_const <- function(x) {
  a <- x[1] 
  b <- x[2] 
  
  a + b # budget constraint formula
  
}

eps <- .Machine$double.eps*10^2 # low number, but not zero due to logs
x0 <- c(0.1, 0.1) # starting values
budget <- 3e9 # overall budget constraint value

opt_solution_log <- solnp(pars = x0,
                          fun = opt_func_log,
                          eqfun = equal_const,
                          eqB = budget,
                          LB = c(eps,eps))

不幸的是,我没有得到一个可行的解决方案

Iter: 1 fn: -5032442923.2173     Pars:  213333.43335 213333.43335
solnp-->Redundant constraints were found. Poor
solnp-->intermediate results may result.Suggest that you
solnp-->remove redundant constraints and re-OPTIMIZE

Iter: 2 fn: -5032442923.2173     Pars:  213333.43335 213333.43335
solnp--> Solution not reliable....Problem Inverting Hessian.

我做错了什么?在这个问题中,什么约束是多余的?是我错误地定义了这个问题,还是它就是不能以这种方式解决?

推荐答案

下面是两个选项来解决100.


If you work with base R

如果需要设置linear个约束,可以使用constrOptim,例如,

ui <- rbind(-c(1, 1), c(1, 0), c(0, 1))
ci <- c(-3e9, 0, 0)
theta <- runif(2)
constrOptim(theta,
    f = opt_func_log,
    grad = NULL,
    ui = ui,
    ci = ci,
    control = list(reltol = .Machine$double.eps)
)

你就会得到

$par
[1] 1609771281 1390228719

$value
[1] -40508597161

$counts
function gradient
     754       NA 

$convergence
[1] 0

$message
NULL

$outer.iterations
[1] 3

$barrier.value
[1] -6339423

If you would like to use fmincon from package pracma

library(pracma)
fmincon(runif(2),
    opt_func_log,
    A = t(c(1, 1)),
    b = 3e9,
    lb = c(0, 0),
    tol = .Machine$double.eps
)

这给了我们

$par
[1] 1609771176 1390228824

$value
[1] -40508597161

$convergence
[1] 0

$info
$info$lambda
$info$lambda$lower
     [,1]
[1,]    0
[2,]    0

$info$lambda$upper
     [,1]
[1,]    0
[2,]    0

$info$lambda$ineqlin
[1] 4.604494 0.000000 0.000000


$info$grad
          [,1]
[1,] -4.604494
[2,] -4.604494

$info$hessian
             [,1]         [,2]
[1,] 1.778707e+00 5.176693e-09
[2,] 5.176693e-09 2.693166e-09

R相关问答推荐

在数据表中呈现数学符号

如何创建构成多个独立列条目列表的收件箱框列?

使用R中相同值创建分组观测指标

为什么在ggplot2中添加geom_text这么慢?

将非重复序列高效转换为长格式

在R中为马赛克图中的每个字段着色

2个Rscript.exe可执行文件有什么区别?

R Select()可以测试不存在的子集列

如何基于两个条件从一列中提取行

如何将一列中的值拆分到R中各自的列中

在使用具有Bray-Curtis相似性的pvCluust时计算p值

Geom_arcbar()中出错:找不到函数";geom_arcbar";

将统计检验添加到GGPUBR中的盒图,在R

如果条件匹配,则使用Mariate粘贴列名

将文本批注减少到gglot的y轴上的单个值

使用其他DF中的文件名将列表中的每个元素保存到文件中

对计算变量所有唯一值的变量进行变异

带查找数据的FCT_REORDER.帧

即使使用相同的种子,mtry值也取决于TuneGrid范围

根据单个条件变异多个列