R中lpSolve求解线性规划崩溃问题排查求助

问题描述 #

使用lpSolve包求解大规模0-1线性规划时，仅因目标函数向量不同，出现两种截然不同的结果：

• 采用f.obj1作为目标函数时，程序可正常运行并返回结果
• 采用f.obj2作为目标函数时，直接导致R 4.2.0崩溃

推测第二个案例对应的线性规划问题不可行，但小规模不可行问题中lpSolve会正常返回“无解”提示，大规模场景下却触发崩溃。需找到避免崩溃的解决办法或lpSolve配置方案。

复现代码 #

library(lpSolve)
# 读取约束矩阵与目标函数数据
f.con <- as.matrix(read.csv('https://raw.githubusercontent.com/hdoran/Blog/master/SOData/f.con.csv'))
dat <- read.csv('https://raw.githubusercontent.com/hdoran/Blog/master/SOData/dat.csv')

constraintMat <- structure(list(
  constraint = c("total", "strand1", "strand2", "strand3", "strand4", 
                 "tminus1min", "tminus1max", "at0min", "at0max", "tplus1min", "tplus1max"),
  value = c(100L, 25L, 25L, 25L, 25L, 5L, 8L, 22L, 26L, 40L, 45L),
  direction = c("=", "=", "=", "=", "=", ">=", "<=", ">=", "<=", ">=", "<=")
), class = "data.frame", row.names = c(NA, -11L))

K <- 2000
f.dir <- c(constraintMat$direction, rep("<=", K))
f.rhs <- c(constraintMat$value, rep(1, K)) 

# 正常运行的案例
result <- lp("max", dat$f.obj1, f.con, f.dir, f.rhs, all.bin=TRUE) 
# 导致R崩溃的案例
result <- lp("max", dat$f.obj2, f.con, f.dir, f.rhs, all.bin=TRUE) 

可能原因 #

• lpSolve的分支定界求解器在处理大规模不可行0-1规划问题时，可能存在内存管理漏洞或数值稳定性问题，触发异常内存访问导致R崩溃。
• 大规模约束矩阵与不可行目标函数的组合，可能导致求解器陷入异常计算流程，无法正常返回“无解”状态。

解决办法 #

1. 提前验证问题可行性 #

在求解原目标函数前，先求解一个以全0向量为目标的可行性问题，提前判断问题是否可行，避免触发崩溃：

# 求解可行性问题：目标函数设为全0，仅验证是否存在可行解
feas_result <- lp("max", rep(0, ncol(f.con)), f.con, f.dir, f.rhs, all.bin=TRUE)
if (feas_result$status != 0) {
  cat("问题不可行，终止求解\n")
} else {
  # 问题可行时再求解原目标函数
  result <- lp("max", dat$f.obj2, f.con, f.dir, f.rhs, all.bin=TRUE)
}

2. 开启预处理优化 #

启用lpSolve的presolve参数，让求解器提前对问题进行预处理，识别不可行性或简化问题，降低崩溃概率：

result <- lp("max", dat$f.obj2, f.con, f.dir, f.rhs, all.bin=TRUE, presolve=TRUE)

3. 切换至更稳健的求解器包 #

若lpSolve的稳定性问题无法解决，可尝试其他专门的线性规划包：

• lpSolveAPI：lpSolve的API封装版本，提供更精细的控制和更可靠的错误处理：

library(lpSolveAPI)
# 创建LP模型，指定约束数与变量数
lprec <- make.lp(nrow = length(f.dir), ncol = ncol(f.con))
# 设置目标函数与最大化方向
set.objfn(lprec, dat$f.obj2)
lp.control(lprec, sense = "max")
# 批量设置约束
for (i in seq_along(f.dir)) {
  set.row(lprec, i, f.con[i, ])
  set.constr.type(lprec, i, f.dir[i])
  set.rhs(lprec, i, f.rhs[i])
}
# 设置所有变量为二进制
set.type(lprec, 1:ncol(f.con), "binary")
# 求解模型
solve_status <- solve(lprec)
if (solve_status != 0) {
  cat("问题不可行或求解失败\n")
} else {
  cat("最优目标值：", get.objective(lprec), "\n")
}

• ROI：支持多种求解器后端的统一接口，可切换至更稳定的求解器（如GLPK、CPLEX等）。

4. 检查数值稳定性 #

确认目标函数与约束矩阵的数值无极端值（过大/过小），必要时对目标函数进行缩放（如除以最大值），避免求解器出现数值溢出。

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