编写一个函数来优化非线性函数 nleqslv

Question

我有一个我想求解的非线性方程组，但我发现很难在 nleqslv.

中编写优化函数

这是我想用数学方法做的；我想最小化:

所以 a 值是常量，我想搜索 最小化 这个总和的 x 值。

问题是对我的价值观的疯狂限制。从数学上讲，顺序是：

所以，如果 k

但是，所有 x 的总和必须等于 0:

最后，每个 xi 都可以是负数（第一个除外）并受最小值和最大值约束，这些值是它之前所有 xs 的累积和的函数:

我在 R 中设置了假值，以解决此问题的简单版本：

a <- c(20, 34, 22, 27)
scm <- 9300
finj <- function(x, inj_max){
  (1/(10 * x+1)) * inj_max
}

fwit <- function(x, wit_max){
  log( x+ 1) * wit_max
}
INJ <- 4650
WIT <- 4650

这些函数将最后一个约束转化为：

fn <- function(x){
  
  #(0 <= x1) can't be expressed - so I put it x1 + 0.0000001
  c(
    x[1] - scm,
    x[1] + 0.0000001,
    x[1] + x[2] - scm,
    x[1] + x[2] + 0.0000001,
    x[1] + x[2] + x[3] - scm,
    x[1] + x[2] + x[3] + 0.0000001,
    x[1] + x[2] + x[3] + x[4],
    #now start inj/wit constraints
    x[2] * (10 * x[1] +1) - INJ,
    x[2] + log(x[1] +1) * WIT + 0.0000001,
    x[3] * (10 * (x[1] + x[2]) +1) - INJ,
    x[3] + log(x[1] + x[2] +1) * WIT + 0.0000001,
    x[4] * (10 * (x[1] + x[2] + x[3]) +1) - INJ,
    x[4] + log(x[1] + x[2]  + x[3] + 1) * WIT+ 0.0000001
    )
  
}

nleqslv(c(4650, -4650, 4650, -4650), fn)

我也写了这个 function 并试图解决它但是我得到了错误：

Error in nleqslv(c(4650, -4650, 4650, -4650), fn) : 
  Length of fn result <> length of x!

我得到这个错误是合乎逻辑的，因为我有太多的约束，所以我不知道如何解决这个优化问题或者我如何重写约束以避免这个错误。

Answer 1

我认为您不需要 nleqslv，它适用于非线性方程组。您正在尝试最小化具有多个参数的单个函数。 optim 来自 base R 的应该可以工作。

至于约束，每个参数都有最小值和最大值，但边界是顺序相关的，这使得它有点棘手。一种方法是将输入顺序转换为允许的 space。这允许函数接受任何实际值作为输入，因为它会自动转换它们以满足约束。我使用 pnorm 进行转换。

要考虑的另一件事是该问题具有 N - 1 自由度，因为 sum(x) 必须为 0。处理该问题的方法是仅将 N - 1 参数传递给要优化的函数，然后将 x[N] 设置为 -sum(x[-N]).

下面是一些使用您的“假值”的示例代码：

scm <- 9300
INJ <- 4650
WIT <- 4650

a <- c(20, 34, 22, 27)

fT <- function(xT) {
  # transforms the input values xT into values that meet the problem constraints
  x <- numeric(length(xT) + 1)
  mini <- 0             # the minimum for parameter 1
  maxi <- min(scm, INJ) # the maximum for parameter 1
  x[1] <- (maxi - mini)*(pnorm(xT[1])) + mini # transform xT[1] to a value between mini and maxi
  xcumsum <- x[1]
  
  for (i in 2:length(xT)) {
    mini <- max(-xcumsum, -WIT*log(xcumsum + 1))     # calculate the minimum for parameter i
    maxi <- min(scm - xcumsum, INJ/(10*xcumsum + 1)) # calculate the maximum for parameter i
    x[i] <- (maxi - mini)*(pnorm(xT[i])) + mini      # transform xT[i] to a value between mini and maxi
    xcumsum <- xcumsum + x[i]
  }
  
  x[i + 1] <- -xcumsum
  return(x)
}

fn <- function(xT) {
  return(sum(a*fT(xT)))
}

# optimize fn using a vector of N - 1 zeros as the initial guess
> optim(numeric(length(a) - 1), fn)
$par
[1]  17.3 -23.2   9.2

$value
[1] -88350

$counts
function gradient 
      32       NA 

$convergence
[1] 0

$message
NULL

optim$par 返回的值是转换后的值。使用 fT:

反转转换

x <- fT(optim(numeric(length(a) - 1), fn)$par)
> x
[1]  4650 -4650  4650 -4650

将 x[1:3] 传递给 fn 给出最小化的函数值：

> fn(head(x, -1))
[1] -88350

通过手动计算 fn 结帐：

> sum(a*x)
[1] -88350

更新 1：从关于收敛到次优局部最小值的评论中，我尝试了适用于此处的 optim 的各种可用方法，并且“L-BFGS-B”方法确实找到了这种情况下的全局，但它是很难说它是否会普遍收敛到全局最小值：

cm <- 4650
INJ <- 4650
WIT <- 4650
a <- c(20, 19, 22, 27)

> optim(numeric(length(a) - 1), fn)$value
[1] -32550
> optim(numeric(length(a) - 1), fn, method = "BFGS")$value
[1] -32550
> optim(numeric(length(a) - 1), fn, method = "CG")$value
[1] -32550
> optim(numeric(length(a) - 1), fn, method = "L-BFGS-B")$value
[1] -37200
> optim(numeric(length(a) - 1), fn, method = "SANN")$value
[1] -32550.1

更新 2：为了回答有关上述代码如何处理 N - 1 参数的问题，我将指出几点：

一个N - 1长度向量被传递给optim（参见numeric(length(a) - 1)）
fT 接受一个向量 (xT) 并输出一个长度为 length(xT) + 1 的向量（参见 x <- numeric(length(xT) + 1) 和 x[i + 1] <- -xcumsum）
我从未创建对象 N，但是 i = N - 1 一旦 for 循环完成，所以 x[i + 1] <- -xcumsum 影响 x 与 [=43] 相同=] 因为 xcumsum 是滞后的累积总和

编写一个函数来优化非线性函数 nleqslv

Writing a function to optimize non-linear function nleqslv

optimization

r

nonlinear-optimization