如何用R中线性与混合整数规划构建组以最小化组内亲缘关系？

核心思路 #

将最小化组内个体间亲缘关系总和的问题转化为**整数线性规划(ILP)**问题：

1. 引入二进制变量标记个体是否属于某组
2. 引入辅助变量标记一对个体是否同组，将二次目标转化为线性形式
3. 添加约束确保每组雄雌数量符合要求、每个个体仅属于一个组
4. 使用Rsymphony求解ILP得到最优分组

可复现代码 #

1. 加载依赖包并模拟数据

library(Rsymphony)
library(kinship2)

# 模拟带亲本信息的动物数据
set.seed(123)
n_females <- 8
n_males <- 4
ped <- data.frame(
  id = c(paste0("F", 1:n_females), paste0("M", 1:n_males)),
  father = c(rep(NA, 4), sample(paste0("M", 1:4), 4, replace=TRUE), rep(NA, 2), sample(paste0("M", 1:2), 2, replace=TRUE)),
  mother = c(rep(NA, 4), sample(paste0("F", 1:4), 4, replace=TRUE), rep(NA, 2), sample(paste0("F", 1:2), 2, replace=TRUE)),
  sex = c(rep("F", n_females), rep("M", n_males))
)

# 计算个体间亲缘系数矩阵
kin_mat <- kinship(ped$id, ped$father, ped$mother)
rownames(kin_mat) <- colnames(kin_mat) <- ped$id

2. 定义分组参数

G <- 2                  # 指定组数
m_per_group <- 2        # 每组雄性数量
f_per_group <- 4        # 每组雌性数量
total_ind <- nrow(ped)  # 总个体数

3. 构建ILP模型

3.1 定义变量

# x_ig: 个体i属于组g的二进制变量
x_vars <- expand.grid(id=ped$id, group=1:G)
x_vars$name <- paste0("x_", x_vars$id, "_g", x_vars$group)

# z_ijg: 个体i和j同属组g的二进制变量（仅i<j的配对）
pairs <- combn(ped$id, 2, simplify=FALSE)
z_vars <- expand.grid(pair=lapply(pairs, paste, collapse="_"), group=1:G)
z_vars$name <- paste0("z_", z_vars$pair, "_g", z_vars$group)

all_vars <- c(x_vars$name, z_vars$name)
n_vars <- length(all_vars)

3.2 目标函数（最小化组内亲缘关系总和）

obj_coeff <- rep(0, n_vars)
# 为z变量赋值对应亲缘系数
for (k in seq_len(nrow(z_vars))) {
  pair <- strsplit(z_vars$pair[k], "_")[[1]]
  i <- pair[1]
  j <- pair[2]
  var_idx <- which(all_vars == z_vars$name[k])
  obj_coeff[var_idx] <- kin_mat[i, j]
}

3.3 约束条件

# 约束1: 每个个体仅属于一个组
constraint1 <- matrix(0, nrow=total_ind, ncol=n_vars)
for (i in seq_len(total_ind)) {
  id <- ped$id[i]
  var_idx <- which(grepl(paste0("x_", id, "_g"), all_vars))
  constraint1[i, var_idx] <- 1
}
constraint1_dir <- rep("==", total_ind)
constraint1_rhs <- rep(1, total_ind)

# 约束2: 每组雄性数量符合要求
male_ids <- ped$id[ped$sex == "M"]
constraint2 <- matrix(0, nrow=G, ncol=n_vars)
for (g in 1:G) {
  var_names <- paste0("x_", male_ids, "_g", g)
  var_idx <- match(var_names, all_vars)
  constraint2[g, var_idx] <- 1
}
constraint2_dir <- rep("==", G)
constraint2_rhs <- rep(m_per_group, G)

# 约束3: 每组雌性数量符合要求
female_ids <- ped$id[ped$sex == "F"]
constraint3 <- matrix(0, nrow=G, ncol=n_vars)
for (g in 1:G) {
  var_names <- paste0("x_", female_ids, "_g", g)
  var_idx <- match(var_names, all_vars)
  constraint3[g, var_idx] <- 1
}
constraint3_dir <- rep("==", G)
constraint3_rhs <- rep(f_per_group, G)

# 约束4: 关联z变量与x变量（确保z_ijg=1当且仅当i和j同属组g）
constraint4 <- matrix(0, nrow=3*length(pairs)*G, ncol=n_vars)
constraint4_dir <- character(3*length(pairs)*G)
constraint4_rhs <- numeric(3*length(pairs)*G)
row_idx <- 1

for (pair in pairs) {
  i <- pair[1]
  j <- pair[2]
  for (g in 1:G) {
    z_var <- paste0("z_", paste(pair, collapse="_"), "_g", g)
    x_i_var <- paste0("x_", i, "_g", g)
    x_j_var <- paste0("x_", j, "_g", g)
    z_idx <- match(z_var, all_vars)
    x_i_idx <- match(x_i_var, all_vars)
    x_j_idx <- match(x_j_var, all_vars)
    
    # z_ijg ≤ x_ig
    constraint4[row_idx, c(z_idx, x_i_idx)] <- c(1, -1)
    constraint4_dir[row_idx] <- "<="
    constraint4_rhs[row_idx] <- 0
    row_idx <- row_idx + 1
    
    # z_ijg ≤ x_jg
    constraint4[row_idx, c(z_idx, x_j_idx)] <- c(1, -1)
    constraint4_dir[row_idx] <- "<="
    constraint4_rhs[row_idx] <- 0
    row_idx <- row_idx + 1
    
    # z_ijg ≥ x_ig + x_jg - 1
    constraint4[row_idx, c(z_idx, x_i_idx, x_j_idx)] <- c(1, -1, -1)
    constraint4_dir[row_idx] <- ">="
    constraint4_rhs[row_idx] <- -1
    row_idx <- row_idx + 1
  }
}

# 合并所有约束
all_constraints <- rbind(constraint1, constraint2, constraint3, constraint4)
all_dir <- c(constraint1_dir, constraint2_dir, constraint3_dir, constraint4_dir)
all_rhs <- c(constraint1_rhs, constraint2_rhs, constraint3_rhs, constraint4_rhs)

4. 求解并输出结果

# 定义变量类型：所有变量均为二进制(0/1)
var_types <- rep("B", n_vars)

# 调用Rsymphony求解ILP
result <- Rsymphony_solve_LP(
  obj = obj_coeff,
  mat = all_constraints,
  dir = all_dir,
  rhs = all_rhs,
  types = var_types,
  max = FALSE  # 最小化目标函数
)

# 提取分组结果
x_sol <- result$solution[grepl("x_", all_vars)]
names(x_sol) <- gsub("x_", "", names(x_sol))

groups <- lapply(1:G, function(g) {
  ids <- names(x_sol)[grepl(paste0("_g", g), names(x_sol)) & x_sol == 1]
  gsub(paste0("_g", g), "", ids)
})
names(groups) <- paste0("Group", 1:G)

# 按性别输出分组
for (g in names(groups)) {
  cat(g, ":\n")
  cat("  雄性:", paste(groups[[g]][groups[[g]] %in% male_ids], collapse=", "), "\n")
  cat("  雌性:", paste(groups[[g]][groups[[g]] %in% female_ids], collapse=", "), "\n\n")
}

关键说明 #

1. 亲缘系数计算：使用kinship2包基于亲本数据计算标准亲缘系数，若已有亲缘矩阵可直接替换kin_mat
2. 模型转化：通过引入z_ijg变量将二次目标（组内亲缘总和）转化为线性形式，适配Rsymphony的线性规划求解能力
3. 扩展性：可通过修改G、m_per_group、f_per_group参数调整分组规则
4. 性能提示：当个体数量较多时，变量和约束数会显著增加，可考虑使用更高效的ILP求解器（如Gurobi、CPLEX）或简化模型

内容的提问来源于stack exchange，提问作者Fernando Brito Lopes