R语言nloptr包NLOPT_LN_AUGLAG算法STRING_ELT错误排查及等式约束优化问题解决咨询

问题背景 #

我正在使用R语言的nloptr包优化一个自定义函数，核心代码如下：

library('nloptr')
hn <- function(x, n) {
  hret <- 0
  if (n == 0) {
    hret <- 1
    return (hret)
  } else if (n == 1) {
    hret <- 2*x
    return (hret)
  } else {
    hn2 <- 1
    hn1 <- 2*x
    all_n <- seq(from = 2, to = n, by = 1)
    for (ni in all_n) {
      hn = (2*x*hn1/sqrt(ni)) + (2*sqrt( (ni-1)/ni)*hn2)
      #print(hn)
      hn2 = hn1
      hn1 = hn
    }
    hret <- hn
    return (hret)
  }
}
term <- function(alpha, r, theta, n) {
  beta = alpha*cosh(r) - Conj(alpha)*exp(1i*theta)*(sinh(r))
  hnterm <- beta/(sqrt(exp(1i*theta)*sinh(2*r)))
  term4 <- hn(hnterm, n)
  logterm1 <- (1/2)*log(cosh(r))
  logterm2 <- -((1/2)*(abs(alpha)^2)) + ((1/2)* (Conj(alpha)^2))*exp(1i*theta)*tanh(r)
  logterm3 <- (n/2)*( log (((1/2)*exp(1i*theta)*tanh(r)) ))
  logterm4 <- log ( term4)
  logA <- logterm1 + logterm2 + logterm3 + logterm4
  A <- exp(logA)
  retval <- c(A)
  return (A)
}
PESQ <- function(x, alpha) {
  p0 <- x[1]
  p1 <- x[2]
  beta <- x[3]
  r <- x[4]
  theta <- x[5]
  N <- 30
  NI <- seq(from = 0, to = N, by = 1)
  elements <- rep(0+1i*0, length(NI))
  elements_abs_sqr <- rep(0, length(NI))
  pr <- rep(0, length(NI))
  total <- 0 + 1i*0
  for (n in NI) {
    w <-term(2*alpha + beta, r, theta, n)
    elements[n+1] <- w
    elements_abs_sqr[n+1] <-(abs(w)^2)
  }
  total <- sum(elements_abs_sqr)
  for (n in NI) {
    pr[n+1] <- Re(elements[n+1]/sqrt(total))
    pr[n+1] <- pr[n+1]^2
  }
  p_off_given_on <- pr[1]
  elements <- rep(0+1i*0, length(NI))
  elements_abs_sqr <- rep(0, length(NI))
  pr <- rep(0, length(NI))
  total <- 0 + 1i*0
  for (n in NI) {
    w <-term(beta, r, theta, n)
    elements[n+1] <- w
    elements_abs_sqr[n+1] <-(abs(w)^2)
  }
  total <- sum(elements_abs_sqr)
  for (n in NI) {
    pr[n+1] <- Re(elements[n+1]/sqrt(total))
    pr[n+1] <- pr[n+1]^2
  }
  p_on_given_off = 1 - pr[1]
  P_e = p0*p_off_given_on + p1*p_on_given_off
  return(P_e)
}
eval_g_eq <- function(x) {
  return ( x[1] + x[2] - 1)
}
lb <- c(0, 0, -Inf, 0.001, -pi)
ub <- c(1, 1, Inf, Inf, pi)
local_opts <- list("algorithm" = "NLOPT_LD_MMA", "xtol_rel"=1.0e-18)
# Set optimization options.
opts <- list("algorithm" = "NLOPT_LN_AUGLAG", "xtol_rel" = 1.0e-18, "local_opts" = local_opts, "maxeval" = 10000)
x0 <- c(0.1,0.9, 0.1, 0.01, 0.7853982)
alpha <- 0.65
eval_g_ineq <- function(x) {
  return (c (- x[1] - x[2], x[1] + x[2] - 1) )
}
eval_f <- function(x) {
  ret = PESQ(x, alpha)
  return(ret)
}
res <- nloptr ( x0 = x0, eval_f = eval_f, eval_g_eq = eval_g_eq, eval_g_ineq = eval_g_ineq, lb = lb, ub = ub, opts = opts )
print(res)

运行时遇到如下错误：

Error in nloptr(x0 = x0, eval_f = eval_f, eval_g_ineq = eval_g_ineq, eval_g_eq = eval_g_eq, : 
  STRING_ELT() can only be applied to a 'character vector', not a 'NULL'
Calls: ... withCallingHandlers -> withVisible -> eval -> eval -> nloptr
Execution halted

当改用"algorithm"="NLOPT_LN_COBYLA"并移除等式约束时代码能正常运行，但我的需求必须保留等式约束，请问该如何修复？

问题分析与解决方案 #

这个错误的核心原因有两个：算法组合不兼容和冗余约束导致的冲突，下面是具体的修复步骤：

1. 更换兼容的局部优化算法

你使用的全局算法NLOPT_LN_AUGLAG是**无导数（LN=Local No-derivative）的增广拉格朗日算法，但搭配的局部算法NLOPT_LD_MMA是有导数（LD=Local Derivative-based）**的算法，这种跨类型的组合会导致nloptr内部处理出错，触发STRING_ELT()错误。

需要把局部算法换成无导数类型的（以LN开头），比如NLOPT_LN_NELDERMEAD或NLOPT_LN_COBYLA：

# 修改局部优化器为无导数类型
local_opts <- list("algorithm" = "NLOPT_LN_NELDERMEAD", "xtol_rel"=1.0e-18)

2. 移除冗余的不等式约束

你的eval_g_ineq函数中定义的约束c(-x[1]-x[2], x[1]+x[2]-1)和等式约束x[1]+x[2]-1=0完全冲突：

• -x[1]-x[2] ≤ 0 等价于 x[1]+x[2] ≥ 0（虽然这个约束本身没问题，但和等式约束重复）
• x[1]+x[2]-1 ≥ 0 加上等式约束x[1]+x[2]=1，其实没有额外意义，反而会让优化器的可行域判断混乱。

如果没有其他需要的不等式约束，直接删除eval_g_ineq相关代码，并在nloptr调用中移除这个参数：

# 删除冗余的eval_g_ineq函数
# eval_g_ineq <- function(x) {
#   return (c (- x[1] - x[2], x[1] + x[2] - 1) )
# }

# 修改nloptr调用，移除eval_g_ineq参数
res <- nloptr ( x0 = x0, eval_f = eval_f, eval_g_eq = eval_g_eq, lb = lb, ub = ub, opts = opts )

3. 确认初始值满足等式约束

你的初始值x0 <- c(0.1,0.9, 0.1, 0.01, 0.7853982)中x0[1]+x0[2]=1，刚好满足等式约束，这一点做得很好，能帮助优化器更快收敛，不需要修改。

修改后的完整代码片段 #

library('nloptr')
# 保留你的hn、term、PESQ、eval_g_eq、eval_f函数不变...

lb <- c(0, 0, -Inf, 0.001, -pi)
ub <- c(1, 1, Inf, Inf, pi)
# 更换为无导数的局部优化算法
local_opts <- list("algorithm" = "NLOPT_LN_NELDERMEAD", "xtol_rel"=1.0e-18)
opts <- list("algorithm" = "NLOPT_LN_AUGLAG", "xtol_rel" = 1.0e-18, "local_opts" = local_opts, "maxeval" = 10000)
x0 <- c(0.1,0.9, 0.1, 0.01, 0.7853982)
alpha <- 0.65

# 移除冗余的不等式约束
res <- nloptr ( x0 = x0, eval_f = eval_f, eval_g_eq = eval_g_eq, lb = lb, ub = ub, opts = opts )
print(res)

这样修改后，代码应该能正常运行并处理等式约束了。

内容的提问来源于stack exchange，提问作者Astra Uvarova - Saturn's star