带多不等式约束的R语言非线性优化代码错误排查求助

问题背景 #

数据集包含以下字段：

• Phs Code：药店唯一标识
• Score：药店关联评分
• MKT：药店市场价值

优化问题定义：

• 决策变量：X(i) = 第i家药店的投资额
• 目标函数：最大化 $\sum Y(i)$，其中 $Y(i) = 0.272*\log(X(i)) + Score(i)$
• 约束条件：
  1. $\exp(0.0272*\log(X(i)) + Score(i)) < 0.8*MKT(i)$
  2. $\exp(0.0272*\log(X(i)) + Score(i)) - \exp(0.0272*\log(X(i)+1) + Score(i)) \geq 1$
  3. $\sum X(i) < 18000000$（总预算限制）
  4. $X(i) > 0$（所有药店投资额为正）

当前编写代码 #

# Objective function
objfun <- function(x, score) 0.272 * log(x) + my_data$score

# Define the constraints
con <- function(x, Score, MKT) {
  con1 <- exp(0.0272 * log(x) + Score) - 0.8 * MKT
  con2 <- exp(0.0272 * log(x) + Score) - exp(0.0272 * log(x + 1) + Score) - 1
  con3 <- sum(x) - 18000000
}

# Set up the optimization problem
x0 <- rep(1, nrow(my_data))
lb <- rep(0, nrow(my_data))
ub <- rep(Inf, nrow(my_data))

problem <- list(
  objective = objfun,
  lb = lb,
  ub = ub,
  eval_g_ineq = con )

result <- nloptr(x0 = x0, 
                 eval_f = problem$objective, 
                 lb = problem$lb, 
                 ub = problem$ub, 
                 eval_g_ineq = problem$eval_g_ineq, 
                 opts = list("algorithm" = "NLOPT_LN_COBYLA"))

代码问题排查与修正 #

1. 目标函数未求和且参数传递错误 #

目标函数要求最大化所有Y(i)的和，当前代码仅返回每个Y(i)的向量，未执行求和操作；同时直接引用全局变量my_data$score而非函数传入的参数，会导致参数传递混乱。

修正：

objfun <- function(x, score) {
  sum(0.272 * log(x) + score)  # 对所有Y(i)求和，使用传入的score参数
}

2. 约束条件符号与nloptr要求不匹配 #

nloptr的eval_g_ineq要求返回的约束满足 $g(x) \leq 0$，当前代码的约束符号不符合要求：

• 约束2要求 $\exp(...) - \exp(...) \geq 1$，需转换为 $1 - (\exp(...) - \exp(...)) \leq 0$，当前代码的con2符号相反，会导致约束失效
• 需将所有约束合并为一个向量返回（逐元素约束+全局约束）

修正约束函数：

con <- function(x, Score, MKT) {
  # 约束1: exp(...) - 0.8*MKT <= 0
  con1 <- exp(0.0272 * log(x) + Score) - 0.8 * MKT
  # 约束2: 1 - [exp(...) - exp(...)] <= 0
  con2 <- 1 - (exp(0.0272 * log(x) + Score) - exp(0.0272 * log(x + 1) + Score))
  # 约束3: sum(x) - 18000000 <= 0
  con3 <- sum(x) - 18000000
  # 返回所有约束的向量
  c(con1, con2, con3)
}

3. 变量下界设置错误 #

约束4要求X(i) > 0，当前lb设为0会导致log(0)无意义，需设置极小正数下界避免计算错误。

修正：

lb <- rep(1e-5, nrow(my_data))  # 替代0，避免log(0)报错

4. nloptr函数参数传递缺失 #

调用nloptr时，eval_f和eval_g_ineq所需的额外参数（score、Score、MKT）需通过opts的eval_f_args和eval_g_ineq_args传入，否则函数无法找到这些参数。

5. 算法参数调整 #

NLOPT_LN_COBYLA是无导数局部优化算法，可增加最大迭代次数、设置收敛精度以提升优化效果。

修正后完整代码 #

# 加载nloptr包
library(nloptr)

# 假设my_data已加载，包含Score和MKT列
# my_data <- read.csv("your_data.csv")  # 示例加载数据

# 目标函数：最大化所有Y(i)的和
objfun <- function(x, score) {
  sum(0.272 * log(x) + score)
}

# 约束函数：返回g(x) <= 0的约束向量
con <- function(x, Score, MKT) {
  # 约束1：每家药店的exp值不超过0.8*MKT
  con1 <- exp(0.0272 * log(x) + Score) - 0.8 * MKT
  # 约束2：边际收益递减满足>=1的要求，转换为<=0形式
  con2 <- 1 - (exp(0.0272 * log(x) + Score) - exp(0.0272 * log(x + 1) + Score))
  # 约束3：总投资不超过预算
  con3 <- sum(x) - 18000000
  # 合并所有约束
  c(con1, con2, con3)
}

# 初始值设置：每家药店初始投资1
x0 <- rep(1, nrow(my_data))
# 下界：极小正数，避免log(0)错误
lb <- rep(1e-5, nrow(my_data))
# 上界：无限制
ub <- rep(Inf, nrow(my_data))

# 调用nloptr，传入额外参数
result <- nloptr(
  x0 = x0,
  eval_f = objfun,
  lb = lb,
  ub = ub,
  eval_g_ineq = con,
  opts = list(
    "algorithm" = "NLOPT_LN_COBYLA",
    "maxeval" = 10000,  # 增加最大迭代次数
    "xtol_rel" = 1e-6,  # 设置收敛精度
    "eval_f_args" = list(score = my_data$Score),  # 传递目标函数的额外参数
    "eval_g_ineq_args" = list(Score = my_data$Score, MKT = my_data$MKT)  # 传递约束函数的额外参数
  )
)

# 查看优化结果
print(result)

内容的提问来源于stack exchange，提问作者Chiara Abergo