在R中求解带二次约束的线性目标函数问题

针对你将Excel Solver（GRG2）迁移至R的需求，这里提供两种可行的局部最优解求解方案，适配Q非正定的非凸约束场景，同时支持对偶问题复用：

方案1：使用nloptr包（匹配GRG2的非线性优化逻辑） #

nloptr支持带约束的非线性局部优化，可指定梯度类或无梯度类算法，完美适配你的问题场景。

原问题求解（最小化线性目标，带二次约束+上下界） #

library(nloptr)

# 定义核心参数
Q <- matrix(c(0, 3897.6, 3897.6, 0), 2)
q <- c(2, 0)
P <- c(1.24, 10)
x_cur <- c(0.1, 0.01)
ref_val <- t(x_cur) %*% Q %*% x_cur + sum(q * x_cur) 

# 线性目标函数
obj_fun <- function(x) {
  sum(P * x)
}

# 二次约束转换为nloptr要求的 g(x) <= 0 形式
constraint_fun <- function(x) {
  ref_val - (t(x) %*% Q %*% x + sum(q * x))
}

# 优化配置：选用带约束的梯度类算法，贴近GRG2逻辑
opts <- list(
  algorithm = "LD_LBFGS",
  xtol_rel = 1e-6,
  maxeval = 1000,
  print_level = 0
)

# 执行求解
result <- nloptr(
  x0 = x_cur,
  eval_f = obj_fun,
  eval_g_ineq = constraint_fun,
  lb = rep(0, 2),
  ub = rep(1, 2),
  opts = opts
)

# 查看结果与约束验证
result$solution
result$objective
t(result$solution) %*% Q %*% result$solution + sum(q * result$solution) >= ref_val - 1e-6

对偶问题求解（最大化二次型，带线性约束） #

假设对偶问题为max x'Qx + q'x，约束为P'x <= 给定常数+0<=x<=1，复用框架如下：

# 对偶目标函数（转为最小化问题）
dual_obj_fun <- function(x) {
  -(t(x) %*% Q %*% x + sum(q * x))
}

# 线性约束转换为 g(x) <= 0 形式
const_val <- t(P) %*% x_cur  # 示例常数
dual_constraint_fun <- function(x) {
  sum(P * x) - const_val
}

# 求解
dual_result <- nloptr(
  x0 = x_cur,
  eval_f = dual_obj_fun,
  eval_g_ineq = dual_constraint_fun,
  lb = rep(0, 2),
  ub = rep(1, 2),
  opts = opts
)

# 还原最大化目标值与解
-(dual_result$objective)
dual_result$solution

方案2：optim配合标准化惩罚函数（解决尺度问题） #

通过变量标准化统一尺度，结合optim的L-BFGS-B算法（自带上下界约束），规避惩罚项尺度失衡问题：

# 变量标准化：x = 0.5*(y+1)，将x∈[0,1]转换为y∈[-1,1]
obj_penalty <- function(y) {
  x <- 0.5*(y + 1)
  current_q <- t(x) %*% Q %*% x + sum(q * x)
  # 约束不满足时添加惩罚项
  if (current_q < ref_val - 1e-6) {
    return(sum(P * x) + 1e6 * (ref_val - current_q)^2)
  } else {
    return(sum(P * x))
  }
}

# 执行优化
optim_result <- optim(
  par = 2*x_cur - 1,  # 转换为标准化初始值
  fn = obj_penalty,
  method = "L-BFGS-B",
  lower = rep(-1, 2),
  upper = rep(1, 2),
  control = list(maxit = 1000, reltol = 1e-6)
)

# 还原原始变量并验证
optim_x <- 0.5*(optim_result$par + 1)
optim_x
t(optim_x) %*% Q %*% optim_x + sum(q * optim_x) >= ref_val - 1e-6

关键说明 #

• nloptr的LD_LBFGS算法基于梯度优化，与GRG2核心逻辑一致，收敛效率高，适合非凸约束场景；
• 若无需梯度计算，可替换为LN_COBYLA算法（无梯度约束优化），兼容性更强；
• 标准化惩罚函数方法通过统一变量尺度解决惩罚项失衡问题，适合快速适配简单场景；
• 两种方案的代码框架可直接复用，仅需调整目标函数与约束定义即可处理对偶问题。

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