使用Rsolnp求解优化问题时遇非数值参数二元运算符错误

需要在R中求解带约束的非凸优化问题：

最大化目标函数detection_cost(x)，约束条件包括：

1. W是概率分布（sum(W)=1，0≤W_i≤1）
2. 传感器与源的距离D_ij满足0<D_ij≤R
3. X/Y为二进制变量，分别表示传感器是否使用、源是否被检测
4. 已使用传感器数量≤n，已检测源数量≤m

目标函数中，检测概率通过逻辑函数计算：D_j = 1 / (1 + exp(-beta_i * sum(Alpha_i * D_ji * X_j)))

用户编写代码后运行出现错误：

Error in fun(pars, ...) - eqB : non-numeric argument to binary operator

原代码如下：

library(Rsolnp)

set.seed(123)
m <- 3  # number of sources
n <- 4  # number of sensors
D <- matrix(runif(m * n, 1, 10), nrow = m)  # initial distances from each sensor to each source
beta=runif(m,0,1)
Alpha=runif(n, 0, 1)
R=25

W <- runif(m, 0, 1)  # random weights of each source
W <- W / sum(W)  # normalize weights so that they sum to 1
Y <- rep(1, m)  # binary variables indicating which sources are detected
X=round(runif(n,0,1)) # binary variables indicating which sensors are used

param=list(W,beta,Alpha,D,X,Y)
x0 <- unlist(param)

detection_cost <- function(x) {

  W <- x[1:m]
  beta <- x[seq(m+1,2*m)]
  Alpha <- x[seq(2*m+1,2*m+n)]
  D <- t(matrix(x[seq(3*m+2,1+n*m+3*m)],ncol=m,nrow=n))
  Y=rev(x)[seq(1,m)]
  X=rev(x)[seq(m+1,m+n)]
  D_j=1 / (1 + exp(-beta * sum(Alpha * D * X) ))
  
  sum(Y*W*D_j)

}

# Define the constraints
constraints <- function(x) {
  
  # W must be a probability distribution
  
  eq1 <- sum(x[1:m]) - 1
  
  ineq1 <- x[1:m] - 0
  
  ineq2 <- x[1:m] - 1
  
  # The distances D_j between the n sensors and the m sources must satisfy D_j <= R
  
  ineq3 <- x[seq(3*m+2,1+n*m+3*m)] - R
  
  # The distances D_j between the n sensors and the m sources must satisfy D_j > 0
  
  ineq4 <- -x[seq(3*m+2,1+n*m+3*m)] 
  
  # For any target i, any sensor j, D_ji <= R
  
  ineq5 <- t(matrix(x0[seq(3*m+2,1+n*m+3*m)],ncol=m,nrow=n)) - R
  
  # The sum of used sensors should not exceed n
  
  ineq6 <- sum(rev(x)[seq(m+1,m+n)]) - n
  
  # The sum of used targets should not exceed m
  
  ineq7 <- sum(rev(x)[seq(1,m)]) - m
  
  return(list(eq1=eq1, ineq1=ineq1, ineq2=ineq2, ineq3=ineq3, ineq4=ineq4, ineq5=ineq5, ineq6=ineq6 , ineq7=ineq7))
}

# Set the initial values for the optimization variables

lb <- c(rep(0, m), rep(-10, m) , rep(-10, n) , rep(0,n*m) , rep(0, m),rep(0, n))
ub <- c(rep(1, m), rep(10, m) , rep(10, n), rep(R,n*m), rep(1, m),rep(1, n) )

# Solve the optimization problem using the L-BFGS-B method
sol <- solnp(par=x0, fun=detection_cost, eqfun=constraints, LB=lb, UB=ub, control=list(trace=1))

1. solnp约束函数返回格式错误：solnp要求eqfun必须返回数值向量，而原代码返回了包含多个元素的列表，导致内部运算时类型不匹配，触发"非数值参数"错误。
2. 变量索引计算错误：
   • 提取D矩阵时，索引seq(3*m+2,1+n*m+3*m)完全错误，正确的起始位置应为2*m + n + 1（W占m位，beta占m位，Alpha占n位，之后是D的n*m个元素）
   • 使用rev(x)提取X/Y极易出错，应该通过固定索引范围直接提取
3. 约束逻辑错误：ineq5错误引用了初始值x0，而非当前优化参数x，导致约束无法随迭代更新
4. 目标函数维度不匹配：sum(Alpha * D * X)是全局求和，无法得到每个源对应的检测概率，导致向量运算维度错误
5. 二进制变量处理问题：solnp不支持整数约束，原代码仅设置上下限0-1，后续需要专门处理，但这不是当前错误的直接原因

library(Rsolnp)

set.seed(123)
m <- 3  # 源数量
n <- 4  # 传感器数量
R <- 25

# 初始化参数
W_init <- runif(m, 0, 1)
W_init <- W_init / sum(W_init)  # 归一化为概率分布
beta_init <- runif(m, 0, 1)
Alpha_init <- runif(n, 0, 1)
D_init <- matrix(runif(m * n, 1, 10), nrow = m)  # m行n列：每行对应一个源到n个传感器的距离
Y_init <- rep(1, m)
X_init <- round(runif(n, 0, 1))

# 构造初始参数向量：按W, beta, Alpha, D(按列展开), Y, X的顺序
x0 <- c(W_init, beta_init, Alpha_init, as.vector(D_init), Y_init, X_init)

# 目标函数：最大化，所以返回负数值（solnp默认最小化）
detection_cost <- function(x) {
  # 提取各变量
  W <- x[1:m]
  beta <- x[(m+1):(2*m)]
  Alpha <- x[(2*m+1):(2*m+n)]
  D <- matrix(x[(2*m+n+1):(2*m+n+m*n)], nrow = m, ncol = n)  # m行n列，对应源到传感器的距离
  Y <- x[(2*m+n+m*n+1):(2*m+n+m*n+m)]
  X <- x[(2*m+n+m*n+m+1):length(x)]
  
  # 逐源计算检测概率D_j
  D_j <- numeric(m)
  for (i in 1:m) {
    # 第i个源的检测概率：beta[i] * sum(Alpha * D[i,] * X)
    linear_term <- beta[i] * sum(Alpha * D[i,] * X)
    D_j[i] <- 1 / (1 + exp(-linear_term))
  }
  
  # 最大化目标函数，solnp默认最小化，所以返回负值
  -sum(Y * W * D_j)
}

# 等式约束函数：返回向量，满足eqfun(x) = 0
eq_constraints <- function(x) {
  # W的和为1
  sum_W <- sum(x[1:m]) - 1
  return(sum_W)
}

# 不等式约束函数：返回向量，满足ineqfun(x) <= 0
ineq_constraints <- function(x) {
  # 1. W >= 0：x[1:m] - 0 >=0 → 转换为 -(x[1:m]) <=0
  ineq_W_lower <- -x[1:m]
  # 2. W <=1：x[1:m] -1 <=0
  ineq_W_upper <- x[1:m] - 1
  # 3. D >0：-D <=0 → -x[D区域] <=0
  D_idx <- (2*m+n+1):(2*m+n+m*n)
  ineq_D_lower <- -x[D_idx]
  # 4. D <=R：x[D区域] - R <=0
  ineq_D_upper <- x[D_idx] - R
  # 5. 已使用传感器数量 <=n：sum(X) -n <=0
  X_idx <- (2*m+n+m*n+m+1):length(x)
  ineq_X_sum <- sum(x[X_idx]) - n
  # 6. 已检测源数量 <=m：sum(Y) -m <=0
  Y_idx <- (2*m+n+m*n+1):(2*m+n+m*n+m)
  ineq_Y_sum <- sum(x[Y_idx]) - m
  
  # 合并所有不等式约束
  c(ineq_W_lower, ineq_W_upper, ineq_D_lower, ineq_D_upper, ineq_X_sum, ineq_Y_sum)
}

# 设置变量上下限
lb <- c(rep(0, m),          # W: 0-1
        rep(-10, m),        # beta: 宽松范围
        rep(-10, n),        # Alpha: 宽松范围
        rep(0, m*n),        # D: 0-R
        rep(0, m),          # Y: 0-1
        rep(0, n))          # X:0-1
ub <- c(rep(1, m),
        rep(10, m),
        rep(10, n),
        rep(R, m*n),
        rep(1, m),
        rep(1, n))

# 求解优化问题
sol <- solnp(par = x0,
             fun = detection_cost,
             eqfun = eq_constraints,
             eqB = 0,  # 等式约束目标值
             ineqfun = ineq_constraints,
             ineqLB = rep(-Inf, length(ineq_constraints(x0))),  # 不等式下界：-Inf <= ineqfun(x) <=0
             ineqUB = rep(0, length(ineq_constraints(x0))),
             LB = lb,
             UB = ub,
             control = list(trace = 1))

# 查看结果
print(sol)
# 提取最优目标值（注意转换为正值，因为目标函数返回的是负值）
optimal_value <- -sol$values[length(sol$values)]
cat("最优目标值：", optimal_value, "\n")

• 约束函数拆分：将等式约束和不等式约束拆分为两个独立函数，分别返回数值向量，符合solnp的要求
• 变量索引修正：重新计算各变量在参数向量中的位置，避免使用rev(x)这类易出错的操作
• 目标函数调整：由于solnp默认最小化，将最大化问题转换为最小化负目标函数；同时修正检测概率的计算逻辑，逐源处理确保维度匹配
• 约束逻辑修正：将所有不等式约束转换为<=0的形式，并修正ineq5的参数引用错误
• 二进制变量说明：solnp不支持整数约束，这里仅设置0-1的上下限做近似处理。如果需要严格的二进制变量，建议使用nloptr搭配整数约束，或专门的混合整数优化包如ROI

内容的提问来源于stack exchange，提问作者Tou Mou