尧图网页设计 
简介一套面向多目标优化学习与研究者的NSGA3与NSGA-II算法Python/Matlab实现代码包。针对工程设计、调度、投资组合等具有相互冲突目标的优化问题提供从帕累托前沿构建、快速非支配排序到拥挤距离与分层选择的完整实现框架。压缩包共10个文件全部为.m脚本大小仅11KB包括主算法、非支配排序模块、环境选择模块、锦标赛选择算子以及一个可直接运行的示例主程序结构清晰便于逐行阅读、复现和二次开发。目前已有2135人学习下载。借助这份代码读者可直观理解拥挤距离计算、弧形划分子种群、精英保留策略等关键机制并可通过改写目标函数快速验证自己的多目标优化想法配合IGD、HV等指标评估算法性能适合正在学习进化算法或需要搭建多目标优化基线的本科生、研究生与工程师。1. 多目标优化找代码前先弄明白NSGA-II和NSGA-III差在哪多目标优化里没有“最好”只有“不坏”当两个目标的优化方向冲突时最优解是一组互不支配的Pareto前沿而不是单个极值点。很多人搜“NSGA3代码NSGAII多目标算法Python”第一反应是先把NSGA-II跑通目标一多再换NSGA-III。反直觉的结论是目标不超过三个时NSGA-II的拥挤度距离已经够用超过三个拥挤度在高维空间会明显退化NSGA-III用参考点替代它才值得你多写那几十行代码。这篇文章按“排序机制 → 选择机制 → 最小可运行代码 → 参数与验证”的顺序把两个算法的Python实现和替换边界一次说清。代码不依赖任何优化框架只要机器上有能跑NumPy的Python环境就能直接跟。2. 非支配排序与拥挤度NSGA-II在Python里的最小实现2.1 先说清楚“支配”NSGA-II的分层依据所谓个体p支配个体q是指p在所有目标上都不比q差并且至少在一个目标上严格优于q。NSGA-II的第一步就是对这个关系做非支配排序先找出种群中所有不被任何个体支配的解标记为第0层去掉它们再在剩余个体中找出新的一批非支配解标记为第1层如此反复。层号越靠前说明这个解在不恶化任何目标的前提下无法被种群内其它解整体替换。排序结果决定了选择和淘汰的优先级层号小的个体一定优先进入下一代只有同一层内部的个体才需要靠第二个指标——拥挤度距离——来区分好坏。2.2 第一个可抄函数非支配排序的NumPy写法import numpy as np def fast_non_dominated_sort(values): # values: (种群个体数, 目标数) 二维数组 n values.shape[0] dominated_solutions [[] for _ in range(n)] # 被个体 p 支配的集合 dominance_count np.zeros(n, dtypeint) # 支配个体 p 的个数 fronts [[]] # 从第 0 层开始收集 for p in range(n): for q in range(n): if p q: continue # p 支配 qp 在所有目标上不差且至少一个目标严格优于 q if np.all(values[p] values[q]) and np.any(values[p] values[q]): dominated_solutions[p].append(q) elif np.all(values[q] values[p]) and np.any(values[q] values[p]): dominance_count[p] 1 if dominance_count[p] 0: fronts[0].append(p) k 0 while fronts[k]: nxt [] for p in fronts[k]: for q in dominated_solutions[p]: dominance_count[q] - 1 if dominance_count[q] 0: nxt.append(q) k 1 fronts.append(nxt) return fronts[:-1]这段代码是朴素的 O(N²) 双循环100 个个体内性能没问题。dominance_count统计的是每个个体被多少个其它个体支配归零就意味着它的所有支配者都已经进入前面的层可以放到下一层。工程实现里更快的 ENS-BS 排序逻辑相同只是加快了“支配者都已入场”的判断我这里给朴素版是为了把机制讲透。2.3 拥挤度距离同一层里的公平竞争非支配排序只能分先后不能比同层优劣。NSGA-II 的答案是把同一层个体按每个目标分别排序计算相邻距离占总跨度的比例累加得到拥挤度。边界个体直接给无穷大保证 Pareto 前沿两端一定被保留。def crowding_distance(fitness, front): # fitness: 全部个体目标值; front: 某一层的个体索引列表 dist np.zeros(len(front)) for m in range(fitness.shape[1]): order sorted(front, keylambda i: fitness[i, m]) fmin fitness[order[0], m] fmax fitness[order[-1], m] dist[order[0]] dist[order[-1]] np.inf if fmax fmin: continue for i in range(1, len(order) - 1): dist[order[i]] (fitness[order[i1], m] - fitness[order[i-1], m]) / (fmax - fmin) return dist注意代码里把边界个体的距离直接设成np.inf。这意味着只要一个解落在某个目标的极值位置它会优先于同层其它解进入下一代。如果你发现最终前沿缺少端点先查这里是不是被误改成了有限值。2.4 完整最小运行ZDT1 测试函数上的 NSGA-IIZDT1 是双目标测试函数的标准起点变量维度 n30Pareto 前沿是 f1 在 [0,1] 上的一条凸曲线。下面这个脚本直接保存运行不依赖任何第三方优化库。def zdt1(x): g 1 9 * np.sum(x[1:]) / (len(x) - 1) return np.array([x[0], g * (1 - np.sqrt(x[0] / g))]) def sbx_crossover(p1, p2, eta_c20, prob0.9): if np.random.rand() prob: return p1.copy(), p2.copy() c1, c2 p1.copy(), p2.copy() for i in range(len(p1)): if np.random.rand() 0.5: continue u np.random.rand() if u 1e-10: u 1e-10 beta (2 * u) ** (1 / (eta_c 1)) if u 0.5 else (1 / (2 * (1 - u))) ** (1 / (eta_c 1)) c1[i] np.clip(0.5 * ((1 beta) * p1[i] (1 - beta) * p2[i]), 0, 1) c2[i] np.clip(0.5 * ((1 - beta) * p1[i] (1 beta) * p2[i]), 0, 1) return c1, c2 def polynomial_mutation(ind, eta_m20, pm0.1): for i in range(len(ind)): if np.random.rand() pm: continue u np.random.rand() delta (2 * u) ** (1 / (eta_m 1)) - 1 if u 0.5 else 1 - (2 * (1 - u)) ** (1 / (eta_m 1)) ind[i] np.clip(ind[i] delta, 0, 1) return ind pop_size, n_var, generations 100, 30, 200 pop np.random.rand(pop_size, n_var) for gen in range(generations): fit np.array([zdt1(x) for x in pop]) fronts fast_non_dominated_sort(fit) next_pop, pool_info [], [] for rank, front in enumerate(fronts): if len(next_pop) len(front) pop_size: cd crowding_distance(fit, front) for i, idx in enumerate(front): next_pop.append(pop[idx]) pool_info.append((pop[idx], rank, cd[i])) else: need pop_size - len(next_pop) cd crowding_distance(fit, front) order np.argsort(cd)[::-1][:need] # 拥挤度大的先保留 for i in order: idx front[i] next_pop.append(pop[idx]) pool_info.append((pop[idx], rank, cd[i])) break def tournament(): best None for _ in range(2): cand pool_info[np.random.randint(len(pool_info))] if best is None or cand[1] best[1] or (cand[1] best[1] and cand[2] best[2]): best cand return best[0] offspring [] while len(offspring) pop_size: p1, p2 tournament(), tournament() c1, c2 sbx_crossover(p1, p2) offspring.append(polynomial_mutation(c1, pm1 / n_var)) if len(offspring) pop_size: offspring.append(polynomial_mutation(c2, pm1 / n_var)) pop np.array(offspring) if gen % 50 0: print(fgen {gen}: front0 size {len(fronts[0])})这段逻辑要拆开看每个世代先合并“当前种群”和“其产生的子代”这里为了可读性省略了显式父子合并直接在当前种群上做选择、交叉、变异精英性由“层优先 同层拥挤度优先”保证。fronts的前面层永远优先进入下一代当某一层放不下时按拥挤度从大到小补满。pool_info同时记录个体的层级和拥挤度供后面的二元锦标赛使用。sbx_crossover里的eta_c20是经验值它控制子代离父代的远近eta_c越大子代越贴近父代搜索越局部。polynomial_mutation里的pm1/n_var是连续优化里常用的变异概率下界n30 时约 3.3%每个个体平均只有一个变量发生变异。运行结束后pop里所有第0层个体就是当前求得的 Pareto 近似解。参数常见值作用SBX 交叉率0.85~0.95保证种群有足够重组机会eta_c15~30越大子代越接近父代eta_m20~100控制变异步长的分布形态pm1/n_var 到 1连续问题取下界离散问题取大值3. 参考点与生态位NSGA-III的三个关键机制和代码3.1 高维目标为什么不能用拥挤度从几何直觉说起目标个数到四五个以后拥挤度距离会退化。原因不复杂高维目标空间中同一前沿层的个体数量本来就稀疏相邻个体之间的距离不再代表“这个解周围有多少解”而是几乎被维度本身主导距离数值集中在极小范围排序结果近似随机。更麻烦的是Pareto 前沿在高维下往往是超曲面用目标轴上的一维距离累加无法反映曲面上的均匀程度。NSGA-III 的思路是彻底换个度量不再问“你周围挤不挤”而是问“你离我预设的参考方向近不近”。多样性变成对一组参考方向的贴近程度这也是它能在三维及以上目标问题里保持良好分布的原因。3.2 参考点如何生成Das-Dennis 均匀设计NSGA-III 的参考点生成沿用 Das-Dennis 方法对 M 个目标、H 等分枚举所有满足 h1 h2 ... hM H 的非负整数组合然后让每个参考点取 wi hi / H。这些点落在单位超平面上且分布均匀。以三维目标、H12 为例组合数是 C(123-1, 3-1)也就是 C(14,2)91 个参考点。这也是为什么 NSGA-III 里种群规模经常取 91、120、210 这类数字——它们都对应某个整数 H 下的参考点数量。如果种群规模和参考点数差距太大后面会讲生态位选择里的空转问题。3.3 生态位选择NSGA-III 填满最后一层的方式当 NSGA-III 遇到“当前层放不下只能选其中一部分”的情况时它不考虑拥挤度而是用下面的生态位策略def niching_selection(assoc, dist_to_ref, niche_count, critical_ids, need): # assoc: 候选个体关联到的参考点编号 # dist_to_ref: 候选个体与其参考点的距离 # niche_count: 已选个体中每个参考点被占用的次数 # critical_ids: 关键层候选个体在种群中的索引 # need: 还需要从关键层选出的个体数 chosen [] while len(chosen) need: j int(np.argmin(niche_count)) # 关键层里关联到参考点 j 的所有个体 cand [k for k in range(len(critical_ids)) if assoc[k] j] if len(cand) 0: niche_count[j] len(critical_ids) 1 # 禁用该参考点 continue if niche_count[j] 0: best cand[np.argmin(dist_to_ref[cand])] # 该方向第一个解 else: best int(np.random.choice(cand)) # 已有解则随机补充 chosen.append(best) niche_count[j] 1 return chosen这段代码是 NSGA-III 选择机制的骨架。核心逻辑是每次找“当前生态位计数最少的参考点”如果该参考点还没有关联任何已选个体就从关键层里挑距离最近的补上如果已经有就随机挑一个。这样做的好处是让种群逐步铺满所有参考方向而不是让解挤在一两个优势区域。niche_count[j] len(critical_ids) 1这一行的作用是“关闭”没有候选个体的参考点避免死循环。它把计数压到一个绝对不会被argmin选中的大数相当于告诉算法这个方向暂时没人。3.4 NSGA-II 和 NSGA-III 的选择对照对比项NSGA-IINSGA-III多样性机制拥挤度距离参考点关联 生态位计数关键层填充依据距离和最大的个体参考点生态位最少者计算瓶颈非支配排序 O(MN²)排序 关联 O(MN²N·R)适用目标数2~3 稳妥3~10 有效工程常见坑高维时同层个体挤成一团参考点数量与种群不匹配二到三个目标的问题NSGA-II 收敛快、实现简单、调试成本低四个目标以上NSGA-III 的参考点机制优势才明显。把 NSGA-II 升级到 NSGA-III只需要替换关键层填充那一段其余排序、交叉、变异全部不动——这就是下一章要做的改造。4. 用Python把NSGA-II改成NSGA-III参考点生成、归一化与参数表4.1 Das-Dennis参考点生成的可运行代码先用最直观的穷举写法H 不超过 12 时运行只在毫秒级足够撑起教学和大部分工程测试。import itertools def generate_reference_points(M, H): # 枚举 M 个 [0,H] 整数筛选和为 H 的组合 all_combs itertools.product(range(H 1), repeatM) refs [np.array(c, dtypefloat) / H for c in all_combs if sum(c) H] return np.array(refs)H 更大比如 M5, H12时组合数量级会到 C(16,4)1820穷举仍然可行。真正爆掉的是 M10、H6 这种配置组合数是 C(15,9)5005还好但如果你同时把 H 调到 12组合数会到十亿级别。生产环境下推荐用组合计数法生成等分布点也就是对每个变量的累积和做嵌套循环这里不再展开。generate_reference_points只需要在算法开始前调用一次生成结果固定不需要每代重复。4.2 归一化与最近参考点关联NSGA-III 在选择关键层个体前必须先把目标值归一化到 [0,1] 空间否则尺度大的目标会主导距离计算。常见做法是 min-max 归一化后计算每个个体到所有参考点的欧氏距离取最近者作为关联参考点。def normalize(fit): fmin fit.min(axis0) fmax fit.max(axis0) return (fit - fmin) / np.maximum(fmax - fmin, 1e-9) def associate_to_ref(fit_norm, refs): # fit_norm: (n, M) 已归一化目标值 # refs: (R, M) 参考点 diff fit_norm[:, None, :] - refs # (n, R, M) dist np.linalg.norm(diff, axis2) # 每个个体到每个参考点距离 return np.argmin(dist, axis1) # 最近的参考点编号原论文里用的是垂直距离也就是先看个体在参考方向上的投影再算垂直偏差但在目标值已经缩放进 [0,1] 的情况下直接算欧氏距离得到的关联结果几乎一致且实现上少一个投影步骤。注意这里的normalize用的是每个个体自己的最小值最大值它是为“关键层已选个体”这一批数据设计的不要在主循环外复用固定值除非你确定目标边界不随时间变化。4.3 把 NSGA-II 的最后一层替换成生态位选择拿到参考点和关联函数后把第 2 章主循环里“放不下就按拥挤度排序”的那段整体替换成下面这个函数。def nsga3_survivor(fit, fronts, refs, pop_size, selected_idx): selected_idx list(selected_idx) for front in fronts: if len(selected_idx) len(front) pop_size: selected_idx.extend(front) continue need pop_size - len(selected_idx) crit list(front) combined selected_idx crit # 合并归一化保证已选个体和关键层在同一空间下关联 scaled normalize(fit[combined]) assoc associate_to_ref(scaled, refs) niche_count np.zeros(len(refs), dtypeint) for i in range(len(selected_idx)): niche_count[assoc[i]] 1 chosen_in_crit [] while len(chosen_in_crit) need: j int(np.argmin(niche_count)) cand [i for i in range(len(selected_idx), len(combined)) if assoc[i] j] if len(cand) 0: niche_count[j] len(combined) 1 continue if niche_count[j] 0: best min(cand, keylambda i: np.linalg.norm( scaled[i] - refs[j])) else: best int(np.random.choice(cand)) chosen_in_crit.append(best - len(selected_idx)) niche_count[j] 1 selected_idx.extend(np.array(crit)[chosen_in_crit]) break return selected_idx用法很直接fit是当前合并后的全部个体目标值fronts来自fast_non_dominated_sortrefs是参考点数组。当某层能完整放下时仍然走原来的“直接放入”只有到放不下的关键层才归一化、关联、生态位选择三步连做。这个函数返回的是被选中个体在原种群中的索引列表实际取pop[selected_idx]就是下一代。4.4 目标数、H 和种群规模的匹配表用 NSGA-III 前先定参考点数再定种群规模。参考点数等于组合数 C(HM-1, M-1)常用配置如下目标数 MH参考点数常用种群规模312911003141201205621021083120120如果想让种群规模凑整可以用两层参考点边界层取较小的 H1内部层取稍大的 H2把两个集合合并后再用。这样参考点数量可以精确控制到 102、115 这类数字。单层配置跑通之前不建议先上两层——调试时不好分辨是归一化问题还是参考点密度问题。5. 跑通之后的验证与坑IGD、目标尺度、离散变量和种群规模5.1 用IGD曲线验证收敛性一个十行函数算法写完第一件事是验证它真的在收敛而不是看前沿“看起来差不多”。反向世代距离IGD是最简单的指标对真实 Pareto 前沿均匀采样计算每个采样点到当前解集最近距离的平均值越小越好。def igd(true_front, approx_front): true_front np.atleast_2d(true_front) dists [] for p in true_front: dists.append(np.min(np.sqrt(((approx_front - p) ** 2).sum(axis1)))) return np.mean(dists)ZDT1 的真实前沿是 f2 1 - sqrt(f1)在 f1 轴上均匀取 100 个点即可。每 20 代打印一次 IGD如果曲线稳步下降说明排序、选择、交叉链路都在正常工作如果 IGD 徘徊不动先查目标函数是否写错再查关键层替换逻辑里selected_idx是否越界。5.2 目标尺度不一致先统一量纲再做归一化如果两个目标的数值范围相差超过一个数量级比如一个在 [0,1]一个在 [1e3, 2e3]min-max 归一化后每代缩放比例都在变参考点方向会被高量纲目标拉扯。常见做法是在三代之后对目标值做对数变换或按已知量纲系数缩放让所有目标进入同一量级后再交给 NSGA-III。判断标准打印归一化后的关联分布如果大量个体集中到两个参考点多半是尺度问题而不是算法问题。5.3 离散决策变量先映射再取整顺序别反工程优化里大量变量是离散的。直接对 SBX 交叉结果取整会让大量子代落在同一格点上种群多样性下降很快。常见做法是把整数变量映射到 [0,1] 连续区间参与交叉和变异只在评价目标函数前取整并裁剪到合法范围x_real np.clip(ind, 0, 1) x_int np.round(x_real * (ub - lb) lb).astype(int)注意取整动作放在目标函数计算入口处而不是放在种群进化的输出端。变异概率在离散问题里不要超过 0.2否则取整后近似随机搜索。5.4 H、种群规模和参数调整的先后顺序先把 H 按目标数定好再让种群规模等于参考点数或略大最后才调 eta_c 和 eta_m——这个顺序反过来NSGA-III 的生态位选择就会变成玄学。本文还有配套的精品资源点击获取 NSGA-II与NSGA-III算法解析：Python实现多目标优化与选择指南
简介一套面向多目标优化学习与研究者的NSGA3与NSGA-II算法Python/Matlab实现代码包。针对工程设计、调度、投资组合等具有相互冲突目标的优化问题提供从帕累托前沿构建、快速非支配排序到拥挤距离与分层选择的完整实现框架。压缩包共10个文件全部为.m脚本大小仅11KB包括主算法、非支配排序模块、环境选择模块、锦标赛选择算子以及一个可直接运行的示例主程序结构清晰便于逐行阅读、复现和二次开发。目前已有2135人学习下载。借助这份代码读者可直观理解拥挤距离计算、弧形划分子种群、精英保留策略等关键机制并可通过改写目标函数快速验证自己的多目标优化想法配合IGD、HV等指标评估算法性能适合正在学习进化算法或需要搭建多目标优化基线的本科生、研究生与工程师。1. 多目标优化找代码前先弄明白NSGA-II和NSGA-III差在哪多目标优化里没有“最好”只有“不坏”当两个目标的优化方向冲突时最优解是一组互不支配的Pareto前沿而不是单个极值点。很多人搜“NSGA3代码NSGAII多目标算法Python”第一反应是先把NSGA-II跑通目标一多再换NSGA-III。反直觉的结论是目标不超过三个时NSGA-II的拥挤度距离已经够用超过三个拥挤度在高维空间会明显退化NSGA-III用参考点替代它才值得你多写那几十行代码。这篇文章按“排序机制 → 选择机制 → 最小可运行代码 → 参数与验证”的顺序把两个算法的Python实现和替换边界一次说清。代码不依赖任何优化框架只要机器上有能跑NumPy的Python环境就能直接跟。2. 非支配排序与拥挤度NSGA-II在Python里的最小实现2.1 先说清楚“支配”NSGA-II的分层依据所谓个体p支配个体q是指p在所有目标上都不比q差并且至少在一个目标上严格优于q。NSGA-II的第一步就是对这个关系做非支配排序先找出种群中所有不被任何个体支配的解标记为第0层去掉它们再在剩余个体中找出新的一批非支配解标记为第1层如此反复。层号越靠前说明这个解在不恶化任何目标的前提下无法被种群内其它解整体替换。排序结果决定了选择和淘汰的优先级层号小的个体一定优先进入下一代只有同一层内部的个体才需要靠第二个指标——拥挤度距离——来区分好坏。2.2 第一个可抄函数非支配排序的NumPy写法import numpy as np def fast_non_dominated_sort(values): # values: (种群个体数, 目标数) 二维数组 n values.shape[0] dominated_solutions [[] for _ in range(n)] # 被个体 p 支配的集合 dominance_count np.zeros(n, dtypeint) # 支配个体 p 的个数 fronts [[]] # 从第 0 层开始收集 for p in range(n): for q in range(n): if p q: continue # p 支配 qp 在所有目标上不差且至少一个目标严格优于 q if np.all(values[p] values[q]) and np.any(values[p] values[q]): dominated_solutions[p].append(q) elif np.all(values[q] values[p]) and np.any(values[q] values[p]): dominance_count[p] 1 if dominance_count[p] 0: fronts[0].append(p) k 0 while fronts[k]: nxt [] for p in fronts[k]: for q in dominated_solutions[p]: dominance_count[q] - 1 if dominance_count[q] 0: nxt.append(q) k 1 fronts.append(nxt) return fronts[:-1]这段代码是朴素的 O(N²) 双循环100 个个体内性能没问题。dominance_count统计的是每个个体被多少个其它个体支配归零就意味着它的所有支配者都已经进入前面的层可以放到下一层。工程实现里更快的 ENS-BS 排序逻辑相同只是加快了“支配者都已入场”的判断我这里给朴素版是为了把机制讲透。2.3 拥挤度距离同一层里的公平竞争非支配排序只能分先后不能比同层优劣。NSGA-II 的答案是把同一层个体按每个目标分别排序计算相邻距离占总跨度的比例累加得到拥挤度。边界个体直接给无穷大保证 Pareto 前沿两端一定被保留。def crowding_distance(fitness, front): # fitness: 全部个体目标值; front: 某一层的个体索引列表 dist np.zeros(len(front)) for m in range(fitness.shape[1]): order sorted(front, keylambda i: fitness[i, m]) fmin fitness[order[0], m] fmax fitness[order[-1], m] dist[order[0]] dist[order[-1]] np.inf if fmax fmin: continue for i in range(1, len(order) - 1): dist[order[i]] (fitness[order[i1], m] - fitness[order[i-1], m]) / (fmax - fmin) return dist注意代码里把边界个体的距离直接设成np.inf。这意味着只要一个解落在某个目标的极值位置它会优先于同层其它解进入下一代。如果你发现最终前沿缺少端点先查这里是不是被误改成了有限值。2.4 完整最小运行ZDT1 测试函数上的 NSGA-IIZDT1 是双目标测试函数的标准起点变量维度 n30Pareto 前沿是 f1 在 [0,1] 上的一条凸曲线。下面这个脚本直接保存运行不依赖任何第三方优化库。def zdt1(x): g 1 9 * np.sum(x[1:]) / (len(x) - 1) return np.array([x[0], g * (1 - np.sqrt(x[0] / g))]) def sbx_crossover(p1, p2, eta_c20, prob0.9): if np.random.rand() prob: return p1.copy(), p2.copy() c1, c2 p1.copy(), p2.copy() for i in range(len(p1)): if np.random.rand() 0.5: continue u np.random.rand() if u 1e-10: u 1e-10 beta (2 * u) ** (1 / (eta_c 1)) if u 0.5 else (1 / (2 * (1 - u))) ** (1 / (eta_c 1)) c1[i] np.clip(0.5 * ((1 beta) * p1[i] (1 - beta) * p2[i]), 0, 1) c2[i] np.clip(0.5 * ((1 - beta) * p1[i] (1 beta) * p2[i]), 0, 1) return c1, c2 def polynomial_mutation(ind, eta_m20, pm0.1): for i in range(len(ind)): if np.random.rand() pm: continue u np.random.rand() delta (2 * u) ** (1 / (eta_m 1)) - 1 if u 0.5 else 1 - (2 * (1 - u)) ** (1 / (eta_m 1)) ind[i] np.clip(ind[i] delta, 0, 1) return ind pop_size, n_var, generations 100, 30, 200 pop np.random.rand(pop_size, n_var) for gen in range(generations): fit np.array([zdt1(x) for x in pop]) fronts fast_non_dominated_sort(fit) next_pop, pool_info [], [] for rank, front in enumerate(fronts): if len(next_pop) len(front) pop_size: cd crowding_distance(fit, front) for i, idx in enumerate(front): next_pop.append(pop[idx]) pool_info.append((pop[idx], rank, cd[i])) else: need pop_size - len(next_pop) cd crowding_distance(fit, front) order np.argsort(cd)[::-1][:need] # 拥挤度大的先保留 for i in order: idx front[i] next_pop.append(pop[idx]) pool_info.append((pop[idx], rank, cd[i])) break def tournament(): best None for _ in range(2): cand pool_info[np.random.randint(len(pool_info))] if best is None or cand[1] best[1] or (cand[1] best[1] and cand[2] best[2]): best cand return best[0] offspring [] while len(offspring) pop_size: p1, p2 tournament(), tournament() c1, c2 sbx_crossover(p1, p2) offspring.append(polynomial_mutation(c1, pm1 / n_var)) if len(offspring) pop_size: offspring.append(polynomial_mutation(c2, pm1 / n_var)) pop np.array(offspring) if gen % 50 0: print(fgen {gen}: front0 size {len(fronts[0])})这段逻辑要拆开看每个世代先合并“当前种群”和“其产生的子代”这里为了可读性省略了显式父子合并直接在当前种群上做选择、交叉、变异精英性由“层优先 同层拥挤度优先”保证。fronts的前面层永远优先进入下一代当某一层放不下时按拥挤度从大到小补满。pool_info同时记录个体的层级和拥挤度供后面的二元锦标赛使用。sbx_crossover里的eta_c20是经验值它控制子代离父代的远近eta_c越大子代越贴近父代搜索越局部。polynomial_mutation里的pm1/n_var是连续优化里常用的变异概率下界n30 时约 3.3%每个个体平均只有一个变量发生变异。运行结束后pop里所有第0层个体就是当前求得的 Pareto 近似解。参数常见值作用SBX 交叉率0.85~0.95保证种群有足够重组机会eta_c15~30越大子代越接近父代eta_m20~100控制变异步长的分布形态pm1/n_var 到 1连续问题取下界离散问题取大值3. 参考点与生态位NSGA-III的三个关键机制和代码3.1 高维目标为什么不能用拥挤度从几何直觉说起目标个数到四五个以后拥挤度距离会退化。原因不复杂高维目标空间中同一前沿层的个体数量本来就稀疏相邻个体之间的距离不再代表“这个解周围有多少解”而是几乎被维度本身主导距离数值集中在极小范围排序结果近似随机。更麻烦的是Pareto 前沿在高维下往往是超曲面用目标轴上的一维距离累加无法反映曲面上的均匀程度。NSGA-III 的思路是彻底换个度量不再问“你周围挤不挤”而是问“你离我预设的参考方向近不近”。多样性变成对一组参考方向的贴近程度这也是它能在三维及以上目标问题里保持良好分布的原因。3.2 参考点如何生成Das-Dennis 均匀设计NSGA-III 的参考点生成沿用 Das-Dennis 方法对 M 个目标、H 等分枚举所有满足 h1 h2 ... hM H 的非负整数组合然后让每个参考点取 wi hi / H。这些点落在单位超平面上且分布均匀。以三维目标、H12 为例组合数是 C(123-1, 3-1)也就是 C(14,2)91 个参考点。这也是为什么 NSGA-III 里种群规模经常取 91、120、210 这类数字——它们都对应某个整数 H 下的参考点数量。如果种群规模和参考点数差距太大后面会讲生态位选择里的空转问题。3.3 生态位选择NSGA-III 填满最后一层的方式当 NSGA-III 遇到“当前层放不下只能选其中一部分”的情况时它不考虑拥挤度而是用下面的生态位策略def niching_selection(assoc, dist_to_ref, niche_count, critical_ids, need): # assoc: 候选个体关联到的参考点编号 # dist_to_ref: 候选个体与其参考点的距离 # niche_count: 已选个体中每个参考点被占用的次数 # critical_ids: 关键层候选个体在种群中的索引 # need: 还需要从关键层选出的个体数 chosen [] while len(chosen) need: j int(np.argmin(niche_count)) # 关键层里关联到参考点 j 的所有个体 cand [k for k in range(len(critical_ids)) if assoc[k] j] if len(cand) 0: niche_count[j] len(critical_ids) 1 # 禁用该参考点 continue if niche_count[j] 0: best cand[np.argmin(dist_to_ref[cand])] # 该方向第一个解 else: best int(np.random.choice(cand)) # 已有解则随机补充 chosen.append(best) niche_count[j] 1 return chosen这段代码是 NSGA-III 选择机制的骨架。核心逻辑是每次找“当前生态位计数最少的参考点”如果该参考点还没有关联任何已选个体就从关键层里挑距离最近的补上如果已经有就随机挑一个。这样做的好处是让种群逐步铺满所有参考方向而不是让解挤在一两个优势区域。niche_count[j] len(critical_ids) 1这一行的作用是“关闭”没有候选个体的参考点避免死循环。它把计数压到一个绝对不会被argmin选中的大数相当于告诉算法这个方向暂时没人。3.4 NSGA-II 和 NSGA-III 的选择对照对比项NSGA-IINSGA-III多样性机制拥挤度距离参考点关联 生态位计数关键层填充依据距离和最大的个体参考点生态位最少者计算瓶颈非支配排序 O(MN²)排序 关联 O(MN²N·R)适用目标数2~3 稳妥3~10 有效工程常见坑高维时同层个体挤成一团参考点数量与种群不匹配二到三个目标的问题NSGA-II 收敛快、实现简单、调试成本低四个目标以上NSGA-III 的参考点机制优势才明显。把 NSGA-II 升级到 NSGA-III只需要替换关键层填充那一段其余排序、交叉、变异全部不动——这就是下一章要做的改造。4. 用Python把NSGA-II改成NSGA-III参考点生成、归一化与参数表4.1 Das-Dennis参考点生成的可运行代码先用最直观的穷举写法H 不超过 12 时运行只在毫秒级足够撑起教学和大部分工程测试。import itertools def generate_reference_points(M, H): # 枚举 M 个 [0,H] 整数筛选和为 H 的组合 all_combs itertools.product(range(H 1), repeatM) refs [np.array(c, dtypefloat) / H for c in all_combs if sum(c) H] return np.array(refs)H 更大比如 M5, H12时组合数量级会到 C(16,4)1820穷举仍然可行。真正爆掉的是 M10、H6 这种配置组合数是 C(15,9)5005还好但如果你同时把 H 调到 12组合数会到十亿级别。生产环境下推荐用组合计数法生成等分布点也就是对每个变量的累积和做嵌套循环这里不再展开。generate_reference_points只需要在算法开始前调用一次生成结果固定不需要每代重复。4.2 归一化与最近参考点关联NSGA-III 在选择关键层个体前必须先把目标值归一化到 [0,1] 空间否则尺度大的目标会主导距离计算。常见做法是 min-max 归一化后计算每个个体到所有参考点的欧氏距离取最近者作为关联参考点。def normalize(fit): fmin fit.min(axis0) fmax fit.max(axis0) return (fit - fmin) / np.maximum(fmax - fmin, 1e-9) def associate_to_ref(fit_norm, refs): # fit_norm: (n, M) 已归一化目标值 # refs: (R, M) 参考点 diff fit_norm[:, None, :] - refs # (n, R, M) dist np.linalg.norm(diff, axis2) # 每个个体到每个参考点距离 return np.argmin(dist, axis1) # 最近的参考点编号原论文里用的是垂直距离也就是先看个体在参考方向上的投影再算垂直偏差但在目标值已经缩放进 [0,1] 的情况下直接算欧氏距离得到的关联结果几乎一致且实现上少一个投影步骤。注意这里的normalize用的是每个个体自己的最小值最大值它是为“关键层已选个体”这一批数据设计的不要在主循环外复用固定值除非你确定目标边界不随时间变化。4.3 把 NSGA-II 的最后一层替换成生态位选择拿到参考点和关联函数后把第 2 章主循环里“放不下就按拥挤度排序”的那段整体替换成下面这个函数。def nsga3_survivor(fit, fronts, refs, pop_size, selected_idx): selected_idx list(selected_idx) for front in fronts: if len(selected_idx) len(front) pop_size: selected_idx.extend(front) continue need pop_size - len(selected_idx) crit list(front) combined selected_idx crit # 合并归一化保证已选个体和关键层在同一空间下关联 scaled normalize(fit[combined]) assoc associate_to_ref(scaled, refs) niche_count np.zeros(len(refs), dtypeint) for i in range(len(selected_idx)): niche_count[assoc[i]] 1 chosen_in_crit [] while len(chosen_in_crit) need: j int(np.argmin(niche_count)) cand [i for i in range(len(selected_idx), len(combined)) if assoc[i] j] if len(cand) 0: niche_count[j] len(combined) 1 continue if niche_count[j] 0: best min(cand, keylambda i: np.linalg.norm( scaled[i] - refs[j])) else: best int(np.random.choice(cand)) chosen_in_crit.append(best - len(selected_idx)) niche_count[j] 1 selected_idx.extend(np.array(crit)[chosen_in_crit]) break return selected_idx用法很直接fit是当前合并后的全部个体目标值fronts来自fast_non_dominated_sortrefs是参考点数组。当某层能完整放下时仍然走原来的“直接放入”只有到放不下的关键层才归一化、关联、生态位选择三步连做。这个函数返回的是被选中个体在原种群中的索引列表实际取pop[selected_idx]就是下一代。4.4 目标数、H 和种群规模的匹配表用 NSGA-III 前先定参考点数再定种群规模。参考点数等于组合数 C(HM-1, M-1)常用配置如下目标数 MH参考点数常用种群规模312911003141201205621021083120120如果想让种群规模凑整可以用两层参考点边界层取较小的 H1内部层取稍大的 H2把两个集合合并后再用。这样参考点数量可以精确控制到 102、115 这类数字。单层配置跑通之前不建议先上两层——调试时不好分辨是归一化问题还是参考点密度问题。5. 跑通之后的验证与坑IGD、目标尺度、离散变量和种群规模5.1 用IGD曲线验证收敛性一个十行函数算法写完第一件事是验证它真的在收敛而不是看前沿“看起来差不多”。反向世代距离IGD是最简单的指标对真实 Pareto 前沿均匀采样计算每个采样点到当前解集最近距离的平均值越小越好。def igd(true_front, approx_front): true_front np.atleast_2d(true_front) dists [] for p in true_front: dists.append(np.min(np.sqrt(((approx_front - p) ** 2).sum(axis1)))) return np.mean(dists)ZDT1 的真实前沿是 f2 1 - sqrt(f1)在 f1 轴上均匀取 100 个点即可。每 20 代打印一次 IGD如果曲线稳步下降说明排序、选择、交叉链路都在正常工作如果 IGD 徘徊不动先查目标函数是否写错再查关键层替换逻辑里selected_idx是否越界。5.2 目标尺度不一致先统一量纲再做归一化如果两个目标的数值范围相差超过一个数量级比如一个在 [0,1]一个在 [1e3, 2e3]min-max 归一化后每代缩放比例都在变参考点方向会被高量纲目标拉扯。常见做法是在三代之后对目标值做对数变换或按已知量纲系数缩放让所有目标进入同一量级后再交给 NSGA-III。判断标准打印归一化后的关联分布如果大量个体集中到两个参考点多半是尺度问题而不是算法问题。5.3 离散决策变量先映射再取整顺序别反工程优化里大量变量是离散的。直接对 SBX 交叉结果取整会让大量子代落在同一格点上种群多样性下降很快。常见做法是把整数变量映射到 [0,1] 连续区间参与交叉和变异只在评价目标函数前取整并裁剪到合法范围x_real np.clip(ind, 0, 1) x_int np.round(x_real * (ub - lb) lb).astype(int)注意取整动作放在目标函数计算入口处而不是放在种群进化的输出端。变异概率在离散问题里不要超过 0.2否则取整后近似随机搜索。5.4 H、种群规模和参数调整的先后顺序先把 H 按目标数定好再让种群规模等于参考点数或略大最后才调 eta_c 和 eta_m——这个顺序反过来NSGA-III 的生态位选择就会变成玄学。本文还有配套的精品资源点击获取