简介:这份 PDF 面向运筹优化、智能算法方向的研究生与工程技术人员,围绕多约束组合优化这一难题展开,选题覆盖金融投资、资源分配、生产调度等实际场景。资源以约束转目标为核心思路,将多约束问题转化为多目标优化问题,并给出改进的多目标粒子群算法 IMaOPSO:用违反约束度维护外部档案,以拥挤度和个体到理想点距离筛选全局最优,同时引入扰动变异算子并让参与变异的粒子数随迭代递减,兼顾开发能力与跳出局部最优。文中还进一步提出基于多种群协同进化的多目标粒子群算法,改进了速度更新机制并设计替换算子以维持种群多样性,最后以不同规模的多背包问题验证算法在约束满足与解质量上的有效性。整包仅 1 个 PDF 文件,约 51KB,体量轻便,便于通读与检索算法细节。目前已有 516 人学习,适合需要复现算法、撰写论文或寻找多目标优化建模思路的读者参考。
1. 多约束组合优化为什么绕不开多目标粒子群算法
互联网业务里的资源分配经常写成这种形式:200 个任务塞进 8 台机器,每台机器的 CPU、内存、带宽都是硬约束,目标还要同时压时延和成本。这就是多约束组合优化的典型形态,也是多背包问题的现实映射——物品有重量和体积,背包有多维容量,要在容量约束下把总价值做大。规模一上来,分支定界和动态规划都很难在可接受时间内收敛。多目标粒子群算法(MOPSO)的价值不在“更聪明”,而在于用种群并行搜索加外部档案保存非支配解,把约束满足和目标优化放进同一条迭代流程。这篇拆的是 IMaOPSO 这套改进:约束转目标、违反约束度维护档案、拥挤度加理想点距离选全局最优、带衰减的扰动变异,以及多种群协同进化版本。适合已经写过单目标 PSO、准备在多约束场景落地的同学。
2. 从单目标 PSO 到多目标 PSO:约束转目标与外部档案
单目标 PSO 只有一个 gbest,粒子朝它飞就行。多约束多目标问题没有单一最优解,解与解之间可能是非支配关系:A 的价值高但违反容量约束,B 完全可行但价值低,C 的价值居中且违反量极小。把其中一个强行定为 gbest,会直接把种群拉向某个偏好方向。多目标粒子群算法要额外解决三件事:怎么表示约束违反、怎么存非支配解、怎么从档案里选全局最优。IMaOPSO 的处理顺序是先把约束转成目标,再用违反约束度维护外部档案,最后用拥挤度和理想点距离选 gbest。
2.1 多背包问题的粒子编码与目标向量
多背包问题(MKP)的定义:n 个物品,m 个背包,物品 j 放进背包 i 需要消耗重量 w_ij,获得价值 p_j,每个物品至多放进一个背包,每个背包 i 的容量上限是 C_i。目标可以设为最大化总价值,约束是每个背包不超容。粒子位置用 (m, n) 的连续矩阵;解码时用 sigmoid 加 argmax,保证每个物品只进一个背包。
import numpy as np def decode(position, m, n): """把连续位置映射成 0-1 分配矩阵。 position: shape=(m, n),m 个背包、n 个物品 return: assign,shape=(m, n),每列恰好一个 1 """ prob = 1.0 / (1.0 + np.exp(-position)) # sigmoid:连续值 -> (0,1) 概率 assign = np.zeros((m, n), dtype=int) idx = np.argmax(prob, axis=0) # 每个物品选概率最大的背包 assign[idx, np.arange(n)] = 1 return assign逻辑说明:sigmoid 把速度域映射到概率域,argmax 保证每列只有一个 1,天然满足“一个物品最多进一个背包”的结构约束。如果问题允许物品不被选中,就在 argmax 之前加一层阈值:max(prob) < 0.5 的列全部置 0。参数上,m 和 n 要与权重矩阵 shape 一致;position 的初始化范围一般取 [-2, 2],太大时 sigmoid 会饱和,梯度信息几乎消失。
2.2 违反约束度:把硬约束变成可比较的目标
约束转目标的核心是把每个容量约束的超出量写成非负值,然后累加。公式可以写成 CV(x) = Σ_i max(0, Σ_j w_ij x_ij - C_i)。总价值取负号后和 CV 一起构成目标向量 F(x) = [-profit(x), CV(x)],都按最小化方向处理。这样不需要为惩罚系数反复调参——惩罚系数法最麻烦的地方是:系数小,约束被无视;系数大,可行域边缘不可探索。
def evaluate(assign, weights, profit, cap): """计算目标向量和违反约束度。 assign: (m, n) 0-1 分配矩阵 weights: (m, n) 每个物品放进每个背包的重量 profit: (n,) 每个物品的价值 cap: (m,) 每个背包的容量上限 """ total_profit = float((assign * profit[None, :]).sum()) used = (assign * weights).sum(axis=1) # 每个背包的实际载重 overload = np.maximum(0.0, used - cap) # 超容量部分 cv = float(overload.sum()) obj = np.array([-total_profit]) # 最大化价值 -> 取负最小化 return obj, cv逻辑说明:used 是各背包的载重向量,overload 只保留正值,cv 是所有超容量之和。目标向量中 -total_profit 与 cv 方向一致,便于后续支配比较。参数 weights 和 cap 要保证量纲一致,如果不同背包的重量尺度差几个数量级,先做归一化,否则 CV 会被大量纲背包主导。
| 约束类型 | 数学表达 | 常用处理 | 注意点 | | 背包容量 | Σ_j w_ij x_ij ≤ C_i | 违反约束度 CV | 量纲统一 | | 物品互斥 | Σ_i x_ij ≤ 1 | 编码阶段解决 | argmax 保证 | | 背包启用 | 未启用则不能放 | 额外 CV 项 | 和容量约束分开算 | | 数量上限 | 单背包物品数 ≤ K | 惩罚或 CV | 惩罚系数要小 |
2.3 外部档案维护:以违反约束度做准入
外部档案存非支配解。支配关系要同时看目标向量和 CV:当 A 的目标都不差于 B,且 A 的 CV 也不大于 B,并且至少有一项严格更优,就称 A 支配 B。更新档案时,先把被新解支配的旧解删掉,再追加新解;超过档案规模就按拥挤度截断。
def dominates(a_obj, a_cv, b_obj, b_cv): """判断解 a 是否支配解 b,目标统一为最小化。""" no_worse = np.all(a_obj <= b_obj) and a_cv <= b_cv strictly = np.any(a_obj < b_obj) or a_cv < b_cv return bool(no_worse and strictly) def update_archive(archive, cand_obj, cand_cv, cand_pos, max_size=100): """archive: list[(obj, cv, pos)]""" for o, c, _ in archive: if dominates(o, c, cand_obj, cand_cv): return archive # 被旧解支配,直接丢弃 archive = [(o, c, p) for o, c, p in archive if not dominates(cand_obj, cand_cv, o, c)] archive.append((cand_obj, cand_cv, cand_pos)) if len(archive) > max_size: archive.sort(key=lambda t: (t[1], np.linalg.norm(t[0]))) archive = archive[:max_size] # 简化截断,生产环境换拥挤度 return archive逻辑说明:第一个循环做快速拒绝,避免无效追加;列表推导删除被新解支配的旧解;截断时先按 CV 升序,再按目标范数升序,等价于“先可行、后收敛”。参数 max_size 取 100 到 300 比较常见:太小则前沿中段容易丢,太大则每代支配比较的平方复杂度会拖慢速度。
提示:存档里的 obj 必须统一方向,如果一部分是 -profit、另一部分是 profit,支配判断会完全反过来。
3. IMaOPSO 落地:全局最优选择、拥挤度与扰动变异
外部档案建好之后,真正决定种群往哪飞的是 gbest 的选择。单目标 PSO 直接从 pbest 和 gbest 取差值;多目标里,档案里每个解都是非支配的,随便挑一个都会带来偏好。IMaOPSO 用两个指标做选择:拥挤度距离衡量解的稀疏程度,个体与理想点的距离衡量解的收敛程度。前者保证前沿铺得开,后者保证不飘在远处。
3.1 拥挤度距离与理想点距离的双指标
拥挤度距离的计算方式:对档案里所有解的第 k 个目标排序,边界解给无穷大,中间解累加相邻目标差与目标范围的比值。理想点 z* = 每个目标在档案里的最小值。个体到理想点的距离用欧氏距离。两个指标量纲不同,先归一化,再加权打分。
def crowding_distance(objs): """objs: (N, M) 全部按最小化方向。返回 (N,) 拥挤度。""" n, m = objs.shape dist = np.zeros(n) for k in range(m): order = np.argsort(objs[:, k]) dist[order[0]] = dist[order[-1]] = np.inf fmin, fmax = objs[order[0], k], objs[order[-1], k] if fmax - fmin < 1e-12: continue for i in range(1, n - 1): nxt = objs[order[i + 1], k] prv = objs[order[i - 1], k] dist[order[i]] += (nxt - prv) / (fmax - fmin) return dist def select_gbest(archive_objs, archive_pos, alpha=0.5): """双指标打分选 gbest:拥挤度越大越好,离理想点越近越好。""" cd = crowding_distance(archive_objs) z = archive_objs.min(axis=0) # 理想点 d = np.linalg.norm(archive_objs - z, axis=1) finite = np.isfinite(cd) cd_n = np.zeros_like(cd) cd_n[finite] = cd[finite] / (cd[finite].max() + 1e-12) cd_n[~finite] = 1.0 # 边界解给满分 d_n = d / (d.max() + 1e-12) score = alpha * cd_n - (1 - alpha) * d_n return archive_pos[int(np.argmax(score))]逻辑说明:cd_n 把无穷大的边界解置为 1,中间解的拥挤度归一化;d_n 越小说明越靠近理想点。alpha 取 0.5 表示探索和开发各占一半,前期可以取 0.7 偏探索,后期取 0.3 偏收敛。参数上,每次迭代都要重新算理想点,不要在第一代算完就锁死,否则后期档案更新后理想点会失真。
3.2 扰动变异算子的衰减策略
IMaOPSO 的扰动变异是为了扩大搜索区域,同时让参与变异的粒子数随迭代减少:前期多变异探索,后期少变异开发。变异个数用 N_mut = max(1, int(N * (1 - t/T)^2)),t 是当前代数,T 是最大代数。变异方式用高斯扰动叠加到位置向量上。
def mutate(pos, t, T, sigma=0.2): """带衰减的扰动变异,返回变异后的位置和被变异的索引。""" N = pos.shape[0] k = max(1, int(N * (1 - t / T) ** 2)) # 变异粒子数随代数递减 idx = np.random.choice(N, size=k, replace=False) pos[idx] += np.random.normal(0.0, sigma, pos[idx].shape) return pos, idx逻辑说明:k 从接近 N 降到 1,前期种群能跳出局部区域,后期只保留少量扰动,避免破坏已经收敛的优异个体。参数 sigma 控制扰动幅度,常取 0.1 到 0.3,超过 0.5 会大量产生不可行解,CV 很难降下来。变异之后必须重新 decode、evaluate 和 update_archive,否则变异只改了位置,没同步到档案。
3.3 一轮完整迭代的拼接
把速度更新、位置更新、pbest 更新和档案更新串起来。速度公式里惯性权重 w 从 0.9 线性降到 0.4,c1、c2 常见各取 1.5 到 2.0,vmax 取位置范围的 10% 到 20%。
def step(swarm, pbest, archive, weights, profit, cap, t, T, c1=1.8, c2=1.8, vmax=0.2): pos, vel = swarm gbest_pos = select_gbest(np.array([o for o, _, _ in archive]), np.array([p for _, _, p in archive])) r1 = np.random.rand(*pos.shape) r2 = np.random.rand(*pos.shape) w = 0.9 - 0.5 * t / T vel = w * vel + c1 * r1 * (pbest[0] - pos) + c2 * r2 * (gbest_pos[None, :, :] - pos) vel = np.clip(vel, -vmax, vmax) pos = pos + vel pos, idx = mutate(pos, t, T) for i in range(pos.shape[0]): assign = decode(pos[i], m=weights.shape[0], n=weights.shape[1]) obj, cv = evaluate(assign, weights, profit, cap) if dominates(obj, cv, pbest[1][i], pbest[2][i]): pbest[0][i], pbest[1][i], pbest[2][i] = pos[i].copy(), obj, cv archive = update_archive(archive, obj, cv, pos[i].copy()) return (pos, vel), pbest, archive逻辑说明:gbest 从外部档案中选,pbest 也用支配关系更新,而不是比较单目标适应度。参数 w 的线性递减让前期速度大、后期速度小,vmax 防止速度爆炸。代码里把 gbest 广播成 (1, m, n) 再加到速度上,保证维度一致。
| 参数 | 常用取值 | 作用 | 调参方向 | | 种群 N | 50-200 | 搜索并行度 | 规模大加大 | | 档案 A | 100-300 | 非支配解容量 | 目标数多时加大 | | w | 0.9 -> 0.4 | 探索到开发 | 后期收敛差就降更慢 | | c1, c2 | 1.5-2.0 | 个体/全局牵引 | 早熟时降 c2 | | vmax | 0.1-0.2 | 速度上限 | 震荡时调小 | | sigma | 0.1-0.3 | 变异幅度 | CV 高时调小 | | alpha | 0.3-0.7 | 拥挤度权重 | 前沿集中时加大 |
注意:pbest 的更新必须比较目标向量和 CV,不能把多目标压成加权和再去比大小,否则档案里的非支配性会被破坏。
4. 多种群协同进化:速度更新与替换算子的工程实现
IMaOPSO 在中小规模还能撑住,但目标空间复杂、问题规模大时,单个种群容易只覆盖前沿的一段。多种群协同进化的思路是把种群拆成 K 个子群,各自负责不同区域,再定期共享信息。IMaOPSO 的多种群版本还改进了速度更新机制,并加入替换算子防止子群停滞。
4.1 多种群划分与协同信息交换
划分方式:按目标空间的偏好向量划分,比如给每个子群一个权重向量,让它偏向某个目标;或者用随机初始化加决策空间聚类。信息交换每隔 R 代做一次:从全局档案里抽非支配解,替换各子群中最差的个体。
def split_swarm(pos, vel, K): """把大种群均分成 K 个子群。""" idx = np.array_split(np.arange(pos.shape[0]), K) return [(pos[i].copy(), vel[i].copy()) for i in idx] def exchange(sub_swarms, global_archive, rng, ratio=0.1): """每 R 代从全局档案抽样,替换子群中最差个体。""" for s in range(len(sub_swarms)): pos, vel = sub_swarms[s] n_replace = max(1, int(pos.shape[0] * ratio)) chosen = rng.choice(len(global_archive), size=n_replace, replace=False) for j, ci in enumerate(chosen): pos[j] = global_archive[ci][2].copy() # 档案中的位置 vel[j] = 0.0 # 替换后速度清零 sub_swarms[s] = (pos, vel) return sub_swarms逻辑说明:拆分后每个子群独立迭代,交换时用全局档案中的解替换最差个体,速度清零是为了避免旧速度把新注入的解立刻带偏。参数 K 取 3 到 5 比较常见,太多会让每个子群样本不足;ratio 取 0.05 到 0.2,过大等于频繁重启,过小起不到信息交换作用。
4.2 速度更新机制改进
单一种群的速度更新只有 pbest 和 gbest 两项。多种群版本加上子群最优 subgbest,形成三层牵引:v = wv + c1r1*(pbest-x) + c2r2(subgbest-x) + c3r3(globalgbest-x)。c1 偏个体经验,c2 偏子群局部经验,c3 偏全局档案。前期 c1 大、c3 小,让子群各自探索;后期反过来,加速收敛。
def update_velocity(pos, vel, pbest, subgbest, globalgbest, t, T, c1=2.0, c2=1.5, c3=1.0, vmax=0.2): r1 = np.random.rand(*pos.shape) r2 = np.random.rand(*pos.shape) r3 = np.random.rand(*pos.shape) w = 0.9 - 0.5 * t / T vel = (w * vel + c1 * r1 * (pbest - pos) + c2 * r2 * (subgbest[None, :, :] - pos) + c3 * r3 * (globalgbest[None, :, :] - pos)) return np.clip(vel, -vmax, vmax)逻辑说明:三项速度分别对应个体记忆、子群局部最优和全局非支配解。c1 到 c3 的取值不需要严格归一化,但总和过大时速度容易震荡,常见做法是先固定 c1=2.0、c2=1.5、c3=1.0,再根据收敛曲线微调。vmax 仍然要卡住,尤其是决策变量做了归一化之后。
4.3 替换算子的触发条件
替换算子解决的是子群停滞:连续若干代没有向外部档案贡献新解,就判定这个子群陷入局部最优。触发后,用全局档案中的随机解替换子群最差个体,或者对档案解做小扰动后注入。
def replace_if_stagnated(sub_pos, sub_vel, sub_stagn, archive_pos, rng, limit=5, sigma=0.1): """连续 limit 代无更新时,注入档案解并重置速度。""" for s in range(len(sub_pos)): if sub_stagn[s] >= limit: worst = np.random.randint(sub_pos[s].shape[0]) donor = archive_pos[rng.randint(len(archive_pos))].copy() donor += rng.normal(0, sigma, donor.shape) # 小扰动,保持多样性 sub_pos[s][worst] = donor sub_vel[s][worst] = 0.0 sub_stagn[s] = 0 return sub_pos, sub_vel, sub_stagn逻辑说明:stagn 计数器在子群产生非支配解时清零,否则每代加一。limit 取 5 到 10,太小会把正常探索误判为停滞,太大则子群已经浪费很多代。sigma 只做小扰动,目的是让注入解不完全等同于档案解,避免多个子群迅速同质化。
| 规模 | 物品数 n | 背包数 m | 种群 N | 档案 A | 迭代 T | | 小 | 50 | 5 | 50 | 100 | 300 | | 中 | 100 | 10 | 100 | 150 | 500 | | 大 | 250 | 20 | 200 | 300 | 800 |
python main.py --items 100 --knapsacks 10 --pop 100 --archive 150 --max_iter 500 --seed 42参数说明:--items 和 --knapsacks 控制问题规模;--pop 是每个子群的粒子数还是总种群数,要看实现,我一般用总种群数均分;--archive 是外部档案上限;--seed 固定随机种子,便于复现。多跑几个种子再比均值,单次结果很容易被初始化运气带偏。
5. 多背包算例验证与调参排错实战
5.1 评价指标:HV、IGD 与约束违反率
多目标算法不能只看“最好解的价值”。常用 HV(超体积):档案解集在目标空间围出的面积/体积,越大越好;IGD 需要参考前沿,衡量解集到参考前沿的平均距离,越小越好;约束违反率是可行解占比。多约束问题里,先过滤 CV < 1e-6 的可行解,再算 HV,否则不可行解会把指标撑大。
import numpy as np, csv def export_and_hv(archive, ref=None, path="front.csv"): """导出外部档案,并在二维目标上近似计算 HV。""" feasible = [(o, c) for o, c, _ in archive if c < 1e-6] objs = np.array([o for o, _ in feasible]) if feasible else np.empty((0, 1)) with open(path, "w", newline="") as f: writer = csv.writer(f) writer.writerow(["profit", "cv"]) for o, c in feasible: writer.writerow([-o[0], c]) # 把 -profit 还原成价值 if objs.shape[0] < 2: return 0.0 if ref is None: ref = objs.max(axis=0) + 1.0 # 参考点要严格劣于所有解 order = np.argsort(objs[:, 0]) hv = 0.0 prev = ref[1] for i in order: hv += max(0.0, ref[1] - objs[i, 1]) * max(0.0, objs[i, 0] - (objs[order[i - 1], 0] if i > 0 else ref[0])) return hv逻辑说明:feasible 只保留满足约束的解,CSV 里把 -profit 还原成真实价值,方便后续画图。ref 是参考点,必须比所有解都差,否则 HV 会偏小。二维 HV 可以用排序法近似,目标数超过三个时换蒙特卡洛采样或专用库。
5.2 常见失败模式与排查顺序
| 现象 | 优先检查 | 处理 | | 外部档案一直只有 1 个解 | 支配比较方向、CV 是否全为 0 | 统一最小化方向,检查重复解码 | | 解全不可行 | sigma 过大、替换算子注入太频繁 | 降 sigma,提高档案中可行解比例 | | 前沿只覆盖一端 | 拥挤度距离失效、alpha 偏小 | 检查 fmax-fmin,调大 alpha | | 后期目标不再改善 | w 太小、变异比例已衰减到 1 | 把 w 下限提到 0.5,保留最小变异 | | 各子群解几乎相同 | 信息交换过于频繁 | 把交换间隔调大,替换 ratio 降到 0.05 |
调试时先把迭代次数压到 30 代,打印每代的档案规模、可行解占比和 CV 最小值。如果 30 代内档案规模不涨,基本可以断定支配比较或档案更新逻辑写错了,而不是算法参数问题。
5.3 参数小技巧
真正要对比不同规模多背包算例时,先用 5 个随机种子跑出可行前沿,按 CV 过滤,再比较 HV 的中位数;只看一次运行的最优价值,很容易被初始化和随机算子带偏。
本文还有配套的精品资源,点击获取