含12个变量的双目标线性规划可行性及求解咨询

一、模型可行性判断 #

你的双目标线性规划模型完全可行，理由如下：

• 决策变量为种植面积，天然满足非负约束（x₁~x₁₂ ≥0）；
• 所有约束条件均为宽松的正数上限（如两类地块上限22666.67、烟草种植上限3333.33），存在大量可行解（比如所有变量取0，或仅种植部分作物），可行域非空；
• 目标函数与约束均为线性关系，符合线性规划的基本要求。

二、Pareto最优解求解方案 #

针对「降低养分盈余（min y₁）+提高产值（max y₂）」的冲突目标，采用ε-约束法生成Pareto前沿（即所有无法同时优化两个目标的最优解集合），以下是R语言实现代码：

1. 加载依赖包与数据 #

# 安装所需包（首次运行需执行）
install.packages(c("lpSolve", "dplyr"))
library(lpSolve)
library(dplyr)

# 导入原始数据
df <- structure(list(Crop = c("Vegetable A", "Vegetable B", "Maize", 
"Barley", "Potato", "Fruit A", "Fruit B", "Rice", "Tabacco", 
"Rape crop", "Faba bean", "Other beans"), `Nutrient surplus (kg/ha)` = c(495, 
495, 287, 269, 330, 355, 355, 226, 194, 203, 130, 137), `Output value (k yuan/ha)` = c(1234.5, 
1234.5, 260.637180923077, 160.344827586207, 798.39552631579, 
1085, 1085, 385.188901345292, 1075.6125, 216.65625, 196.511045454545, 
909), Type = c("A", "B", "A", "B", "B", "A", "B", "A", "C", "B", 
"B", "B"), Area = c("x1", "x2", "x3", "x4", "x5", "x6", "x7", 
"x8", "x9", "x10", "x11", "x12")), class = c("tbl_df", "tbl", 
"data.frame"), row.names = c(NA, -12L))

# 提取目标函数系数
y1_coef <- df$`Nutrient surplus (kg/ha)`
y2_coef <- df$`Output value (k yuan/ha)` * 100  # 转换为元/公顷，与你的y2定义一致

# 定义约束条件矩阵与边界
constraints <- list(
  # 约束1：A类+烟草地块总面积 ≤22666.67
  c1 = list(indices = c(1,3,6,8,9), values = rep(1,5), rhs = 22666.67, dir = "<="),
  # 约束2：B类地块总面积 ≤22666.67
  c2 = list(indices = c(2,4,5,7,10,11,12), values = rep(1,7), rhs = 22666.67, dir = "<="),
  # 约束3：烟草种植面积 ≤3333.33
  c3 = list(indices = 9, values = 1, rhs = 3333.33, dir = "<=")
)

# 整理为lpSolve所需格式
const_mat <- matrix(0, nrow = 3, ncol = 12)
const_dir <- c()
const_rhs <- c()

for (i in 1:length(constraints)) {
  const_mat[i, constraints[[i]]$indices] <- constraints[[i]]$values
  const_dir <- c(const_dir, constraints[[i]]$dir)
  const_rhs <- c(const_rhs, constraints[[i]]$rhs)
}

2. 生成Pareto前沿 #

通过遍历y₁的不同上限值（ε），求解对应y₂的最大值，收集所有Pareto最优解：

# 先求解两个单目标最优解，确定y1的范围
# 单目标1：最小化y1（养分盈余）
lp_min_y1 <- lp("min", objective.in = y1_coef, const.mat = const_mat, 
                const.dir = const_dir, const.rhs = const_rhs, all.int = FALSE, all.pos = TRUE)
min_y1 <- lp_min_y1$objval

# 单目标2：最大化y2（产值）
lp_max_y2 <- lp("max", objective.in = y2_coef, const.mat = const_mat, 
                const.dir = const_dir, const.rhs = const_rhs, all.int = FALSE, all.pos = TRUE)
max_y1_when_max_y2 <- sum(y1_coef * lp_max_y2$solution)

# 生成y1的候选ε值（在min_y1到max_y1_when_max_y2之间取20个点）
epsilon_values <- seq(min_y1, max_y1_when_max_y2, length.out = 20)

# 遍历ε值，求解对应Pareto解
pareto_results <- data.frame()
for (eps in epsilon_values) {
  # 添加y1 ≤ eps的约束
  new_const_mat <- rbind(const_mat, y1_coef)
  new_const_dir <- c(const_dir, "<=")
  new_const_rhs <- c(const_rhs, eps)
  
  # 最大化y2
  lp_result <- lp("max", objective.in = y2_coef, const.mat = new_const_mat, 
                  const.dir = new_const_dir, const.rhs = new_const_rhs, all.int = FALSE, all.pos = TRUE)
  
  if (lp_result$status == 0) {  # 仅保留可行解
    pareto_results <- rbind(pareto_results, data.frame(
      y1 = sum(y1_coef * lp_result$solution),
      y2 = lp_result$objval,
      x1 = lp_result$solution[1],
      x2 = lp_result$solution[2],
      x3 = lp_result$solution[3],
      x4 = lp_result$solution[4],
      x5 = lp_result$solution[5],
      x6 = lp_result$solution[6],
      x7 = lp_result$solution[7],
      x8 = lp_result$solution[8],
      x9 = lp_result$solution[9],
      x10 = lp_result$solution[10],
      x11 = lp_result$solution[11],
      x12 = lp_result$solution[12]
    ))
  }
}

# 去重（避免重复的Pareto点）
pareto_results <- distinct(pareto_results)

# 查看Pareto前沿结果
print(pareto_results)

# 绘制Pareto前沿图
plot(pareto_results$y1, pareto_results$y2, type = "b", xlab = "养分盈余y1(kg)", ylab = "产值y2(元)", 
     main = "Pareto最优前沿：养分盈余vs产值")

3. 结果说明 #

• 输出的pareto_results包含所有Pareto最优解，每行对应一组种植面积分配方案，以及对应的y₁和y₂值；
• 你可以根据实际需求（比如对养分盈余的容忍度、产值目标）从Pareto前沿中选择合适的方案；
• 若需要更精细的前沿，可以增加epsilon_values的数量。

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