基于R的整数线性规划：员工对数工资建模的整数系数求解

问题背景 #

• 建模目标：基于154名员工的工资对数，以154×63的二进制能力矩阵（行=员工，列=工作能力，值为1/0表示是否具备该能力）为预测变量，求解63个整数系数向量x，要求x≥1
• 线性规划定义：
  • 目标函数：$\text{min } \sum(n1 + n2)$（n1为预测值超额向量，n2为预测值不足向量）
  • 约束条件：
    1. 主约束：$Ax + En1 - En2 = y\lambda + b$（E为单位矩阵，对应每个员工的残差项；y为工资对数向量；$\lambda$、b为外生参数）
    2. x必须为整数
    3. x ≥ 1，且n1、n2 ≥ 0（正差值要求）

原代码存在变量维度错误、约束矩阵构造逻辑偏差、目标函数定义不全等问题，导致结果异常，以下是修正后的完整实现：

修正后R代码与解释 #

# 加载依赖包
library(lpSolve)

# ----------------------
# 替换为你的实际输入数据
# ----------------------
# X: 154×63的二进制预测变量矩阵
# grades_and_salaries: 包含salary列的员工数据框
# lambda.min: 来自cv.glmnet的外生参数
# offset_test: 外生参数b的向量

# 1. 定义基础参数
A <- as.matrix(X)
n_employees <- nrow(A)  # 154
n_features <- ncol(A)   # 63
y <- log(grades_and_salaries$salary)
lambda <- lambda.min
b <- as.vector(offset_test)

# 2. 构造目标函数：仅最小化n1和n2的和，x的系数为0
obj <- c(rep(0, n_features),  # x的系数
         rep(1, n_employees), # n1的系数
         rep(1, n_employees))# n2的系数

# 3. 构造主约束矩阵（154行，对应每个员工）
# 主约束：A*x + I*n1 - I*n2 = y*lambda + b（I为单位矩阵）
I_mat <- diag(n_employees)
main_constraint <- cbind(A, I_mat, -I_mat)
main_rhs <- y * lambda + b
main_dir <- rep("==", n_employees)

# 4. 构造x的约束：x ≥ 1（63行）
x_constraint <- diag(n_features, nrow = n_features, ncol = length(obj))
x_rhs <- rep(1, n_features)
x_dir <- rep(">=", n_features)

# 5. 构造n1、n2的非负约束（154+154行，确保残差为正）
n1_constraint <- cbind(matrix(0, n_employees, n_features), 
                       diag(n_employees), 
                       matrix(0, n_employees, n_employees))
n1_rhs <- rep(0, n_employees)
n1_dir <- rep(">=", n_employees)

n2_constraint <- cbind(matrix(0, n_employees, n_features + n_employees), 
                       diag(n_employees))
n2_rhs <- rep(0, n_employees)
n2_dir <- rep(">=", n_employees)

# 6. 合并所有约束
constraint_mat <- rbind(main_constraint, x_constraint, n1_constraint, n2_constraint)
constraint_dir <- c(main_dir, x_dir, n1_dir, n2_dir)
constraint_rhs <- c(main_rhs, x_rhs, n1_rhs, n2_rhs)

# 7. 指定整数变量：仅x为整数（n1/n2为连续值，若需整数可修改int.vec）
int_vars <- 1:n_features

# 8. 求解线性规划
lp_model <- lp(direction = "min",
               objective.in = obj,
               const.mat = constraint_mat,
               const.dir = constraint_dir,
               const.rhs = constraint_rhs,
               int.vec = int_vars)

# 9. 提取结果
if(lp_model$status == 0) {
  x_coeffs <- lp_model$solution[1:n_features]
  cat("成功获取整数系数：\n")
  print(x_coeffs)
} else {
  cat("求解失败，状态码：", lp_model$status, "\n")
}

关键修正点 #

• 目标函数维度修正：原代码仅设置2个元素，实际需覆盖63个x、154个n1、154个n2，共371个变量
• E矩阵修正：将错误的全1矩阵改为单位矩阵，对应每个员工的残差项
• 约束完整性：补充n1、n2的非负约束（正差值的核心要求）
• 整数变量精准指定：用int.vec仅标记x为整数，避免误将n1/n2强制设为整数
• 主约束等式整理：将原约束转换为可直接求解的等式形式，确保右侧计算逻辑正确

常见失败排查 #

• 若求解状态码非0：检查$\lambda$、b的取值是否合理，是否导致约束矛盾
• 若需n1/n2也为整数：将int.vec改为1:length(obj)
• 若求解超时：添加timeout参数限制求解时长，例如lp(..., timeout = 600)

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