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, dtype=int) # 支配个体 p 的个数 fronts = [[]] # 从第 0 层开始收集 for p in range(n): for q in range(n): if p == q: continue # p 支配 q：p 在所有目标上不差，且至少一个目标严格优于 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, key=lambda 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[i+1], m] - fitness[order[i-1], m]) / (fmax - fmin) return dist

注意代码里把边界个体的距离直接设成np.inf。这意味着只要一个解落在某个目标的极值位置，它会优先于同层其它解进入下一代。如果你发现最终前沿缺少端点，先查这里是不是被误改成了有限值。

2.4 完整最小运行：ZDT1 测试函数上的 NSGA-II

ZDT1 是双目标测试函数的标准起点：变量维度 n=30，Pareto 前沿是 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_c=20, prob=0.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_m=20, pm=0.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, pm=1 / n_var)) if len(offspring) < pop_size: offspring.append(polynomial_mutation(c2, pm=1 / n_var)) pop = np.array(offspring) if gen % 50 == 0: print(f"gen {gen}: front0 size = {len(fronts[0])}")

这段逻辑要拆开看：每个世代先合并“当前种群”和“其产生的子代”，这里为了可读性省略了显式父子合并，直接在当前种群上做选择、交叉、变异，精英性由“层优先 + 同层拥挤度优先”保证。fronts的前面层永远优先进入下一代，当某一层放不下时，按拥挤度从大到小补满。pool_info同时记录个体的层级和拥挤度，供后面的二元锦标赛使用。

sbx_crossover里的eta_c=20是经验值，它控制子代离父代的远近：eta_c越大，子代越贴近父代，搜索越局部。polynomial_mutation里的pm=1/n_var是连续优化里常用的变异概率下界，n=30 时约 3.3%，每个个体平均只有一个变量发生变异。运行结束后，pop里所有第0层个体就是当前求得的 Pareto 近似解。

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。这些点落在单位超平面上，且分布均匀。

以三维目标、H=12 为例，组合数是 C(12+3-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-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), repeat=M) refs = [np.array(c, dtype=float) / H for c in all_combs if sum(c) == H] return np.array(refs)

H 更大（比如 M=5, H=12）时，组合数量级会到 C(16,4)=1820，穷举仍然可行。真正爆掉的是 M=10、H=6 这种配置，组合数是 C(15,9)=5005，还好；但如果你同时把 H 调到 12，组合数会到十亿级别。生产环境下推荐用组合计数法生成等分布点，也就是对每个变量的累积和做嵌套循环，这里不再展开。generate_reference_points只需要在算法开始前调用一次，生成结果固定，不需要每代重复。

4.2 归一化与最近参考点关联

NSGA-III 在选择关键层个体前，必须先把目标值归一化到 [0,1] 空间，否则尺度大的目标会主导距离计算。常见做法是 min-max 归一化后，计算每个个体到所有参考点的欧氏距离，取最近者作为关联参考点。

def normalize(fit): fmin = fit.min(axis=0) fmax = fit.max(axis=0) 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, axis=2) # 每个个体到每个参考点距离 return np.argmin(dist, axis=1) # 最近的参考点编号

原论文里用的是垂直距离，也就是先看个体在参考方向上的投影，再算垂直偏差；但在目标值已经缩放进 [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), dtype=int) 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, key=lambda 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_sort，refs是参考点数组。当某层能完整放下时，仍然走原来的“直接放入”；只有到放不下的关键层，才归一化、关联、生态位选择三步连做。这个函数返回的是被选中个体在原种群中的索引列表，实际取pop[selected_idx]就是下一代。

4.4 目标数、H 和种群规模的匹配表

用 NSGA-III 前先定参考点数，再定种群规模。参考点数等于组合数 C(H+M-1, M-1)，常用配置如下：

如果想让种群规模凑整，可以用两层参考点：边界层取较小的 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(axis=1)))) 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 的生态位选择就会变成玄学。

本文还有配套的精品资源，点击获取