如何在R中分析3变量2目标的多目标优化问题？

需要对3个变量（Total A、Total B、Total C）和2个目标（最低Nutrient surplus、最高Output value）进行分析，现有数据如下：

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 (yuan/ha)` = c(123450, 
                                                                                                                                                                                          123450, 26063.7180923077, 16034.4827586207, 79839.552631579, 
                                                                                                                                                                                          108500, 108500, 38518.8901345292, 107561.25, 21665.625, 19651.1045454545, 
                                                                                                                                                                                          90900), Type = c("A", "B", "A", "B", "B", "A", "B", "A", "C", 
                                                                                                                                                                                                           "B", "B", "B")), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
                                                                                                                                                                                                                                                                                    -12L))

约束条件：

The total area for A+C<= 22666.67
The total area for B <= 22666.67
The total area for C <= 3333.33
unit:ha (hectare)

尝试使用gMOIP包但发现不支持3变量，请问可用哪些R包实现该分析？期望得到类似帕累托前沿图的结果。

1. lpSolveAPI + 手动生成帕累托前沿 #

lpSolveAPI是线性规划基础工具，支持通过加权法手动求解多目标问题：将两个目标线性组合，遍历不同权重得到最优解集合，最终生成帕累托前沿。

实现步骤： #

（1）预处理数据，确定各类型的目标系数

先按作物类型分组，提取每个类型中营养盈余最低、产出值最高的单位系数（可根据需求调整选择逻辑）：

library(dplyr)
type_stats <- df %>%
  group_by(Type) %>%
  summarise(
    min_nutrient = min(`Nutrient surplus (kg/ha)`),
    max_output = max(`Output value (yuan/ha)`),
    .groups = "drop"
  )

# 提取目标系数向量
coef_nutrient <- c(
  type_stats$min_nutrient[type_stats$Type=="A"],
  type_stats$min_nutrient[type_stats$Type=="B"],
  type_stats$min_nutrient[type_stats$Type=="C"]
)
coef_output <- c(
  type_stats$max_output[type_stats$Type=="A"],
  type_stats$max_output[type_stats$Type=="B"],
  type_stats$max_output[type_stats$Type=="C"]
)

（2）构建线性规划模型并遍历权重

library(lpSolveAPI)

# 创建含3个变量的LP模型
lp_model <- make.lp(0, 3)

# 添加约束条件
add.constraint(lp_model, c(1, 0, 1), "<=", 22666.67)  # A+C ≤ 22666.67
add.constraint(lp_model, c(0, 1, 0), "<=", 22666.67)  # B ≤ 22666.67
add.constraint(lp_model, c(0, 0, 1), "<=", 3333.33)   # C ≤ 3333.33
set.bounds(lp_model, lower = c(0, 0, 0))  # 面积非负

# 遍历权重生成帕累托点
weights <- seq(0, 1, by = 0.05)
pareto_points <- data.frame()

for(w in weights) {
  # 组合目标：w*营养盈余（最小化） + (1-w)*(-产出值)（转为最小化问题）
  set.objfn(lp_model, w*coef_nutrient + (1-w)*(-coef_output))
  solve(lp_model)
  
  # 提取最优解与目标值
  x <- get.variables(lp_model)
  total_nutrient <- sum(x * coef_nutrient)
  total_output <- sum(x * coef_output)
  
  pareto_points <- rbind(pareto_points, data.frame(total_nutrient, total_output))
}

# 去重冗余帕累托点
pareto_points <- pareto_points[!duplicated(pareto_points), ]

（3）绘制帕累托前沿图

plot(pareto_points$total_nutrient, pareto_points$total_output, 
     type = "l", xlab = "Total Nutrient Surplus (kg)", 
     ylab = "Total Output Value (yuan)", main = "Pareto Front")
points(pareto_points$total_nutrient, pareto_points$total_output, pch = 16)

2. mco包（多目标遗传算法） #

mco专门针对多目标优化，内置NSGA-II等启发式算法，无需手动遍历权重，可直接生成帕累托前沿，适合线性/非线性多目标场景。

实现代码： #

library(mco)

# 定义目标函数：返回需最小化的两个值（营养盈余、-产出值）
obj_fun <- function(x) {
  total_nutrient <- sum(x * coef_nutrient)
  total_output <- sum(x * coef_output)
  return(c(total_nutrient, -total_output))
}

# 定义约束函数：返回值需≤0
const_fun <- function(x) {
  c(
    x[1] + x[3] - 22666.67,
    x[2] - 22666.67,
    x[3] - 3333.33
  )
}

# 设置变量边界
lower_bounds <- c(0, 0, 0)
upper_bounds <- c(22666.67, 22666.67, 3333.33)

# 运行NSGA-II算法
result <- nsga2(obj_fun, idim = 3, odim = 2, 
                lower.bounds = lower_bounds, upper.bounds = upper_bounds,
                constraints = const_fun, cdim = 3,
                popsize = 100, generations = 50)

# 提取帕累托点并绘图
pareto_front <- result$value
pareto_front[,2] <- -pareto_front[,2]  # 转回原始产出值

plot(pareto_front[,1], pareto_front[,2], 
     xlab = "Total Nutrient Surplus (kg)", ylab = "Total Output Value (yuan)",
     main = "Pareto Front (NSGA-II)", pch = 16)

3. ROI + ROI.plugin.glpk #

ROI是统一的优化接口，配合glpk插件可处理线性规划，语法更规范，支持切换多种求解器，同样通过加权法生成帕累托前沿。

实现代码片段： #

library(ROI)
library(ROI.plugin.glpk)

# 定义目标与约束
objective_L1 <- L_objective(coef_nutrient)
objective_L2 <- L_objective(-coef_output)
constraints <- L_constraint(
  matrix(c(1,0,1, 0,1,0, 0,0,1), nrow=3, byrow=TRUE),
  dir = c("<=", "<=", "<="),
  rhs = c(22666.67, 22666.67, 3333.33)
)
bounds <- V_bound(lower = c(0,0,0))

# 遍历权重求解
pareto_points <- data.frame()
for(w in seq(0,1,0.05)){
  obj <- w*objective_L1 + (1-w)*objective_L2
  opt <- OP(obj, constraints, bounds)
  sol <- ROI_solve(opt, solver = "glpk")
  x <- solution(sol)
  total_nutrient <- sum(x*coef_nutrient)
  total_output <- sum(x*coef_output)
  pareto_points <- rbind(pareto_points, data.frame(total_nutrient, total_output))
}

# 绘图
plot(pareto_points$total_nutrient, pareto_points$total_output, type="l", pch=16,
     xlab="Total Nutrient Surplus", ylab="Total Output Value")

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