R语言多机多时段生产计划优化：约束定义与求解方法咨询

我来一步步帮你解决这个生产优化的问题，包括决策变量转换、批量约束定义、solnl错误修复，以及合适的求解函数选择。

1. 把矩阵/数组决策变量转为向量 #

R中的大多数优化函数（包括solnl）要求决策变量是一维向量，而不是矩阵或数组。你可以按以下步骤处理：

假设你的决策变量是一个三维数组x，维度为机器数×天数×产品数（即3×20×15），代表每台机器每天生产每种产品的批次数量。你可以用as.vector()把它展开为向量，之后在目标函数和约束函数里再转回原结构：

# 示例：初始化一个3×20×15的决策变量数组
init_x <- array(0, dim = c(3, 20, 15))
# 转为向量
init_vec <- as.vector(init_x)

# 在函数中将向量转回数组的工具函数
vec_to_array <- function(vec) {
  array(vec, dim = c(3, 20, 15))
}

这样，在目标函数和约束函数里，你可以先把输入的向量转回数组，再进行计算，避免处理向量时的混乱。

2. 批量定义约束（避免重复写20天的规则） #

手动写20天的约束显然不现实，我们可以用循环或向量化操作批量生成约束。下面针对你的约束逐一处理：

约束1：机器生产产品范围 #

机器1只能生产1-4，机器2生产5-12，机器3生产7-15。这个约束是说，不在允许范围内的产品，批次数量必须为0：

# 定义每台机器允许生产的产品范围
allowed <- list(
  m1 = 1:4,
  m2 = 5:12,
  m3 = 7:15
)

# 约束函数：返回所有不允许生产的产品的批次数量（必须等于0）
const_prod_range <- function(vec) {
  x <- vec_to_array(vec)
  constraints <- c()
  
  # 机器1：产品5-15的批次必须为0
  constraints <- c(constraints, as.vector(x[1, , -allowed$m1]))
  # 机器2：产品1-4、13-15的批次必须为0
  constraints <- c(constraints, as.vector(x[2, , -allowed$m2]))
  # 机器3：产品1-6的批次必须为0
  constraints <- c(constraints, as.vector(x[3, , -allowed$m3]))
  
  return(constraints)
}

约束2：每日每台机器的生产产品序列递增 #

假设你这里的“序列递增”是指：每台机器当天生产的产品编号单调非递减（即生产的产品批次对应的编号不能比之前的小）。如果决策变量是批次数量，我们可以转化为：对于每台机器的允许产品，编号小的产品批次不能超过编号大的产品批次，确保序列递增：

# 约束函数：每日每台机器的产品序列递增
const_seq_incr <- function(vec) {
  x <- vec_to_array(vec)
  constraints <- c()
  
  # 遍历每台机器、每天
  for (machine in 1:3) {
    allowed_k <- allowed[[paste0("m", machine)]]
    for (day in 1:20) {
      # 遍历允许的产品对（k1 < k2）
      for (idx in 1:(length(allowed_k)-1)) {
        k1 <- allowed_k[idx]
        k2 <- allowed_k[idx+1]
        # 约束：x[machine, day, k1] ≤ x[machine, day, k2]
        constraints <- c(constraints, x[machine, day, k1] - x[machine, day, k2])
      }
    }
  }
  
  return(constraints)
}

约束3：每种产品总产量满足最低要求 #

假设你有一个长度为15的向量min_prod，代表每种产品的最低总产量。我们先计算每种产品的总批次，再乘以系数，确保总产量≥最低要求：

# 定义每种产品的单批次PC数和单PC重量（示例值，替换为你的实际数据）
pc_per_batch <- rep(745, 15)
weight_per_pc <- rep(0.1, 15)
min_prod <- c(100, 150, 200, 180, 220, 250, 300, 280, 260, 240, 220, 200, 180, 160, 140)

# 约束函数：每种产品总产量≥最低要求
const_min_prod <- function(vec) {
  x <- vec_to_array(vec)
  # 计算每种产品的总批次：跨所有机器和天数求和
  total_batches <- apply(x, 3, sum)
  # 计算总产量
  total_prod <- total_batches * pc_per_batch * weight_per_pc
  # 转化为solnl要求的f(x) ≤0形式：min_prod - total_prod ≤0 → total_prod ≥ min_prod
  return(min_prod - total_prod)
}

约束4：机器3总产量大于机器2 #

计算机器2和机器3的总产量，约束机器3的总产量≥机器2：

# 约束函数：机器3总产量≥机器2
const_m3_gt_m2 <- function(vec) {
  x <- vec_to_array(vec)
  # 计算机器2的总产量：每种允许产品的总批次×对应系数，再求和
  m2_total_batches <- apply(x[2, , allowed$m2], 2, sum)
  m2_prod <- sum(m2_total_batches * pc_per_batch[allowed$m2] * weight_per_pc[allowed$m2])
  
  # 计算机器3的总产量
  m3_total_batches <- apply(x[3, , allowed$m3], 2, sum)
  m3_prod <- sum(m3_total_batches * pc_per_batch[allowed$m3] * weight_per_pc[allowed$m3])
  
  # 转化为solnl要求的f(x) ≤0形式：m2_prod - m3_prod ≤0 → m3_prod ≥ m2_prod
  return(m2_prod - m3_prod)
}

3. 修复solnl的"$ operator is invalid for atomic vectors"错误 #

这个错误通常是因为你在约束函数里把向量当成了数据框/列表来使用$操作符，或者对向量使用了二维索引（比如x[,1]）。解决方法很简单：

• 在约束函数的开头，先把输入的向量转回矩阵/数组（用之前写的vec_to_array函数）；
• 所有操作都基于转回后的数组，而不是原向量。

比如你之前的const_m1函数，修改后应该是：

const_m1_fixed <- function(vec) {
  x <- vec_to_array(vec) # 先转回数组
  f <- c()
  f[1] <- x[1, 1, 1] # 机器1第1天产品1的批次
  f[2] <- -(x[1,1,1] -5) # x[1,1,1] ≤5
  f[3] <- x[1,1,2] -6 # x[1,1,2] ≥6
  f[4] <- -(x[1,1,2] -13) # x[1,1,2] ≤13
  f[5] <- x[1,1,3] -8 # x[1,1,3] ≥8
  return(f)
}

4. 选择合适的求解函数 #

你的问题是线性目标+线性约束（目标函数是各批次的线性组合，约束大多是线性的），可以选择以下两种工具：

选项1：NlcOptim::solnl #

适合处理非线性约束，也支持线性约束。使用时需要把所有约束整合到一个函数里，区分等式和不等式约束：

library(NlcOptim)

# 目标函数：最大化总产量，solnl默认最小化，所以返回负的总产量
obj_fun <- function(vec) {
  x <- vec_to_array(vec)
  total_batches <- apply(x, 3, sum)
  total_prod <- sum(total_batches * pc_per_batch * weight_per_pc)
  return(-total_prod)
}

# 整合所有约束
all_constraints <- function(vec) {
  c(
    const_prod_range(vec), # 等式约束：=0
    const_seq_incr(vec),   # 不等式约束：≤0
    const_min_prod(vec),   # 不等式约束：≤0
    const_m3_gt_m2(vec)    # 不等式约束：≤0
  )
}

# 定义约束类型：等式约束的数量是const_prod_range返回的长度
eq_num <- length(const_prod_range(init_vec))

# 运行优化
result <- solnl(
  x = init_vec,
  objfun = obj_fun,
  confun = all_constraints,
  eqInd = 1:eq_num, # 前eq_num个约束是等式约束（=0）
  ineqInd = (eq_num+1):length(all_constraints(init_vec)) # 剩下的是不等式约束（≤0）
)

选项2：lpSolve::lp #

如果所有约束都是线性的，lpSolve包的lp函数更高效，专门处理线性规划问题，还支持整数约束（如果批次数量必须是整数）：

library(lpSolve)

# 构建目标函数系数向量：每个决策变量对应的系数是pc_per_batch[k] * weight_per_pc[k]
obj_coeff <- rep(0, length(init_vec))
for (k in 1:15) {
  # 找到所有对应产品k的决策变量位置
  positions <- which(matrix(1:15, nrow=3, ncol=20, byrow=TRUE) == k)
  obj_coeff[positions] <- pc_per_batch[k] * weight_per_pc[k]
}

# 构建约束矩阵、方向和右边项（以约束1为例，其他约束同理构建后合并）
# 约束1：不允许的产品批次必须为0
const1_mat <- matrix(0, nrow=length(const_prod_range(init_vec)), ncol=length(init_vec))
row_idx <- 1
for (machine in 1:3) {
  allowed_k <- allowed[[paste0("m", machine)]]
  for (day in 1:20) {
    for (k in setdiff(1:15, allowed_k)) {
      pos <- (machine-1)*20*15 + (day-1)*15 + k
      const1_mat[row_idx, pos] <- 1
      row_idx <- row_idx +1
    }
  }
}
const1_dir <- rep("=", nrow(const1_mat))
const1_rhs <- rep(0, nrow(const1_mat))

# 其他约束（序列递增、最低产量、机器3>机器2）同理构建后合并
# 这里省略其他约束的矩阵构建代码，你可以按相同逻辑扩展

# 运行线性规划（最大化）
lp_result <- lp(
  direction = "max",
  objective.in = obj_coeff,
  const.mat = const1_mat, # 替换为合并后的约束矩阵
  const.dir = const1_dir, # 替换为合并后的约束方向
  const.rhs = const1_rhs, # 替换为合并后的右边项
  all.int = TRUE # 如果批次数量是整数，加上这个参数
)

注意事项 #

1. 如果你的“生产序列递增”约束是逻辑约束（比如必须生产连续的产品），那么这属于整数规划问题，需要引入二进制变量，此时lpSolve或ROI包会更合适；
2. 初始值要设置为满足所有约束的向量，避免优化陷入局部最优；
3. 注意约束的方向转换：solnl默认不等式约束是f(x) ≤0，等式约束是f(x)=0，而lpSolve可以直接指定≥/≤/=。

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