1. 问题描述与挑战
在数据分析、实验设计或样本分配等场景中，我们经常需要将一个包含m个元素的原始数据集（超集）划分为n个互不重叠、且大小预定的子集。一个常见的优化目标是使每个子集的统计特性（例如均值）尽可能地与原始超集的特性保持一致。具体来说，给定一个超集 s 及其包含的 m 个元素，以及 n 个预期的子集大小 x0, x1, ..., xn-1（其中 sum(xi) == m），目标是创建这些子集，使得每个子集的均值与超集的均值最为接近。我们通常通过最小化所有子集均值与超集均值之间绝对差异的总和来量化这一目标。
这是一个典型的组合优化问题，其挑战在于：
- 无放回抽样： 超集中的每个元素只能被分配到一个子集中，且仅使用一次。
- 固定子集大小： 每个子集必须严格满足其预设的元素数量。
- 均值优化： 这是一个全局优化目标，需要权衡不同子集之间的分配，以达到整体最优。
- 计算复杂度： 随着超集元素数量和子集数量的增加，可能的组合呈指数级增长，导致穷举法不可行。在实际应用中，算法需要在合理的时间内（例如1秒内）完成对中等规模数据的处理。
2. 数学建模：混合整数线性规划 (MILP)
这种类型的分配问题可以被归类为“集合划分问题”（Set Partitioning Problem）的一个变种，其中加入了特定的目标函数（均值优化）和额外的约束。混合整数线性规划（MILP）提供了一种强大的框架来精确解决这类问题。
2.1 目标函数
我们的目标是使每个子集的均值 mean(subset_s) 尽可能接近超集的均值 mean(superset)。这等价于使每个子集的总和 sum(subset_s) 尽可能接近 subset_size_s * mean(superset)。因此，我们可以定义目标函数为最小化所有子集总和与目标总和之间绝对差异的总和：
$$ \min \sum{s=0}^{N-1} | \sum{i \in \text{subset}_s} \text{element}_i - (\text{size}_s \times \text{mean}(\text{superset})) | $$
2.2 决策变量
我们引入二元决策变量 x_s_i： x_s_i = 1 如果超集中的第 i 个元素被分配到第 s 个子集中。 x_s_i = 0 否则。
2.3 约束条件
- 子集大小约束： 每个子集 s 必须包含预定数量的元素 size_s。 $$ \sum{i=0}^{M-1} x{s,i} = \text{size}_s \quad \forall s \in {0, \dots, N-1} $$
- 元素唯一性约束： 超集中的每个元素 i 只能被分配到一个且仅一个子集中。 $$ \sum{s=0}^{N-1} x{s,i} = 1 \quad \forall i \in {0, \dots, M-1} $$
- 绝对值线性化： 在线性规划中，通常通过引入辅助变量和不等式来处理绝对值。对于每个子集 s，我们定义其总和误差 err_s： $$ \text{err}s = \sum{i=0}^{M-1} (\text{element}i \times x{s,i}) - (\text{size}_s \times \text{mean}(\text{superset})) $$ 然后引入一个非负辅助变量 abs_err_s 来表示 |err_s|，并添加以下约束： $$ \text{abs_err}_s \ge \text{err}_s $$ $$ \text{abs_err}_s \ge -\text{err}_s $$ 最终，目标函数变为最小化 sum(abs_err_s)。
3. 使用 PuLP 进行求解
PuLP 是一个用 Python 编写的线性规划建模工具，它允许用户以直观的方式定义优化问题，并调用各种求解器（如CBC、GLPK、Gurobi等）来解决。
以下是一个使用 PuLP 解决上述问题的示例代码：
from statistics import mean
import pulp
def solve_subset_partitioning(superset_elements, subset_sizes):
"""
使用 PuLP 解决基于均值优化的数据集子集划分问题。
Args:
superset_elements (list): 超集中的所有元素列表。
subset_sizes (list): N个子集各自的目标大小列表。
Returns:
tuple: (list of lists, list of floats) 分割后的子集元素列表和每个子集的均值。
"""
N = len(subset_sizes)
M = len(superset_elements)
# 验证输入
if sum(subset_sizes) != M:
raise ValueError("所有子集大小之和必须等于超集元素总数。")
# 计算超集均值
superset_mean = mean(superset_elements)
# 创建 PuLP 问题实例
set_partitioning_model = pulp.LpProblem("Set_Partitioning_Model", pulp.LpMinimize)
# 决策变量：x_s_i = 1 如果超集中的第 i 个元素被分配到第 s 个子集中
# covering[s] 是一个列表，其中包含子集 s 的 M 个二元变量
covering = {}
for s in range(N):
vals = []
for i, v in enumerate(superset_elements):
vals.append(
pulp.LpVariable(
f"x_set_{s}_element_idx_{i:>02}_val_{v}",
lowBound=0, # 0
upBound=1, # 1
cat=pulp.LpBinary, # 二进制变量
)
)
covering[s] = vals
# 辅助变量：用于处理绝对误差
abs_sum_errs = []
for s_i in range(N):
abs_sum_errs.append(pulp.LpVariable(f"set_{s_i}_sum_error_abs", lowBound=0))
# 目标函数：最小化所有子集绝对误差之和
set_partitioning_model += pulp.lpSum(abs_sum_errs), "Minimize_Absolute_Sum_Errors"
# 添加约束
for s_i, st_vars in covering.items():
# 计算每个子集的实际总和
current_set_sum = pulp.lpSum([p * superset_elements[i] for i, p in enumerate(st_vars)])
# 计算每个子集的目标总和 (子集大小 * 超集均值)
target_set_sum = subset_sizes[s_i] * superset_mean
# 定义子集总和误差变量
set_sum_err = pulp.LpVariable(f"set_{s_i}_sum_error")
set_partitioning_model += set_sum_err == current_set_sum - target_set_sum, f"Set_{s_i}_Sum_Error_Definition"
# 绝对值线性化约束
set_partitioning_model += abs_sum_errs[s_i] >= set_sum_err, f"Abs_Error_Upper_Bound_Pos_{s_i}"
set_partitioning_model += abs_sum_errs[s_i] >= -set_sum_err, f"Abs_Error_Upper_Bound_Neg_{s_i}"
# 约束1: 每个子集的大小必须符合预设
for s_i, st_vars in enumerate(covering.values()):
set_partitioning_model += pulp.lpSum(st_vars) == subset_sizes[s_i], f"Set_{s_i}_Size_Constraint"
# 约束2: 超集中的每个元素只能被使用一次
# zip(*covering.values()) 将所有子集的变量列表转置，以便按元素索引迭代
for i, element_vars in enumerate(zip(*covering.values())):
set_partitioning_model += (
pulp.lpSum(element_vars) == 1,
f"Element_{i}_Used_Once_Constraint",
)
# 求解模型
set_partitioning_model.solve(pulp.PULP_CBC_CMD(msg=False)) # 使用默认的CBC求解器，静默模式
# 提取结果
if pulp.LpStatus[set_partitioning_model.status] == "Optimal":
result_subsets = []
result_means = []
for s_i, st_vars in covering.items():
current_subset_elements = [
superset_elements[i] for i, var in enumerate(st_vars) if var.value() == 1
]
result_subsets.append(current_subset_elements)
result_means.append(mean(current_subset_elements))
return result_subsets, result_means
else:
print(f"未能找到最优解。状态: {pulp.LpStatus[set_partitioning_model.status]}")
return [], []
# 示例 1：完美分配
print("--- 示例 1：完美分配 ---")
superset1 = [100]*5 + [101]*10 + [102]*5
subset_sizes1 = [2, 4, 14]
subsets1, means1 = solve_subset_partitioning(superset1, subset_sizes1)
print(f"超集均值: {mean(superset1)}")
for i, (subset, mean_val) in enumerate(zip(subsets1, means1)):
print(f"子集 {chr(65+i)} ({len(subset)} 元素): {subset}, 均值: {mean_val}")
# 预期输出：所有子集均值均为 101
# 示例 2：最佳拟合
print("\n--- 示例 2：最佳拟合 ---")
superset2 = [100]*5 + [103]*10 + [104]*5
subset_sizes2 = [2, 4, 14]
subsets2, means2 = solve_subset_partitioning(superset2, subset_sizes2)
print(f"超集均值: {mean(superset2)}")
for i, (subset, mean_val) in enumerate(zip(subsets2, means2)):
print(f"子集 {chr(65+i)} ({len(subset)} 元素): {subset}, 均值: {mean_val}")
# 预期输出：子集均值尽可能接近 102.5代码解析：
- 初始化： 定义超集元素、子集大小，并计算超集均值。
- LpProblem： 创建一个 PuLP 问题实例，目标是最小化 (pulp.LpMinimize)。
- 决策变量 covering： 这是一个字典，键是子集索引 s，值是一个列表，包含了 M 个 pulp.LpBinary 变量。每个变量 x_s_i 代表超集中的第 i 个元素是否属于第 s 个子集。
- 辅助变量 abs_sum_errs： 用于存储每个子集总和与目标总和之间绝对差异的辅助变量，lowBound=0 确保其非负。
- 目标函数： 设置为最小化 abs_sum_errs 中所有元素的和。
- 误差定义与绝对值约束： 遍历每个子集，计算其目标总和 (subset_size * superset_mean)。然后定义 set_sum_err 为实际总和与目标总和之差。最后，通过两个不等式 abs_sum_errs[s_i] >= set_sum_err 和 abs_sum_errs[s_i] >= -set_sum_err 来实现绝对值的线性化。
- 子集大小约束： 确保每个子集中的元素数量（即对应 x_s_i 变量之和）等于其预设的 subset_sizes[s_i]。
- 元素唯一性约束： 确保超集中的每个元素 i（通过 zip(*covering.values()) 遍历）在所有子集中的 x_s_i 变量之和为1，即每个元素仅被分配一次。
- 求解： 调用 set_partitioning_model.solve() 启动求解器。默认情况下，PuLP 会使用其自带的 CBC 求解器。
- 结果提取： 检查求解状态，如果找到最优解，则遍历决策变量，提取每个子集包含的元素，并计算其均值。
4. 启发式算法：Karmarkar-Karp (Largest Differencing Method)
当精确求解（如MILP）因问题规模过大而变得不可行时，启发式算法提供了一种快速获得近似解的方法。Karmarkar-Karp 算法（也称为最大差分法）是解决数集划分问题的一种著名启发式算法，其目标是将一个数集划分为 k 个子集，使这些子集的和尽可能接近，即最小化最大子集和与最小子集和之间的差异。
优点： 速度快，易于实现。 局限性：
- 无法指定子集大小： Karmarkar-Karp 算法的主要目的是平衡子集和，而不是严格控制子集中的元素数量。这与我们问题中“固定子集大小”的要求不符。
- 不直接优化均值： 尽管它试图使子集和接近，但这并不直接等同于使子集均值接近超集均值，尤其是在子集大小不固定的情况下。
因此，Karmarkar-Karp 算法不适用于严格满足本教程最初提出的所有约束条件。然而，作为一种通用的数集划分启发式方法，它在某些场景下仍有其价值，例如当我们只需要大致平衡子集总和，而对子集大小没有严格要求时。
以下是一个使用 numberpartitioning 库实现 Karmarkar-Karp 算法的示例：
from statistics import mean
from numberpartitioning import karmarkar_karp
# 示例 2 的超集数据
superset = [100, 100, 100, 100, 100, 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104]
print("\n--- 启发式方法：Karmarkar-Karp ---")
print("超集均值:", mean(superset))
# 使用 Karmarkar-Karp 划分成 3 个部分
# 注意：此方法不接受预设的子集大小
for p in karmarkar_karp(superset, num_parts=3).partition:
print(f"子集 ({len(p)} 元素): {p}, 均值: {mean(p)}")
从输出可以看出，Karmarkar-Karp 算法生成的子集大小不固定，且其均值与超集均值的接近程度可能不如 MILP 得到的精确解。
5. 性能考量与优化策略
尽管 MILP 可以提供最优解，但其计算复杂度是 NP-hard 的。对于大规模问题，求解时间可能过长。在实际应用中，我们需要根据具体情况权衡求解精度和计算效率。
-
问题规模：
- N (子集数量): 10-25 个子集通常是 MILP 可处理的范围，但达到 100 个子集时，问题会变得非常困难。
- M (超集元素数量): 100-1000 个元素是常见情况，10000 个唯一元素则非常庞大。
- 唯一元素数量： 如果超集中有大量重复元素，可以考虑预处理，将相同元素视为一个“类别”，并为每个类别分配一定数量的元素到子集，这可能简化问题。
-
优化策略：
-
启发式预分配 + 精确调整：
- 初步均匀分配： 按照子集大小比例，从超集中随机（或均匀）抽取元素填充子集 50%-7
-
启发式预分配 + 精确调整：
