NSGA-II实战：用Python手把手教你解决多目标优化问题（附完整代码）

NSGA-II实战：用Python手把手教你解决多目标优化问题（附完整代码）

在工程设计和产品研发中，我们常常面临一个令人头疼的局面：多个目标相互拉扯，难以取舍。比如，设计一款新型电动汽车，你既希望它的续航里程尽可能长，又希望制造成本尽可能低，同时还要求加速性能足够出色。这些目标往往是相互冲突的——提升续航可能需要更大的电池，从而增加成本和重量；追求极致加速性能又可能牺牲能耗。传统的单目标优化方法在这里显得力不从心，因为你无法用一个简单的“总分”来衡量方案的优劣。

这时，多目标优化算法就派上了用场。它不追求单一的“最优解”，而是致力于寻找一组“最佳权衡解”，也就是Pareto最优解集。在这组解里，你无法再改进任何一个目标而不损害其他至少一个目标。对于决策者来说，这就像获得了一份详尽的“菜单”，上面列出了所有可能的“套餐”，你可以根据实际情况（比如预算、技术路线偏好）来选择最合适的那一个。

在众多多目标进化算法中，NSGA-II（带精英策略的快速非支配排序遗传算法） 无疑是应用最广泛、最经典的算法之一。它由Kalyanmoy Deb教授等人于2002年提出，因其高效的排序机制和出色的解集分布性，迅速成为该领域的标杆。今天，我们不谈复杂的数学推导，而是聚焦于如何用Python亲手实现它，并解决一个经典的测试问题。我们将从零开始，一步步构建代码，并深入探讨关键参数如何影响最终结果，让你不仅能“会用”，更能“懂行”。

1. 环境准备与问题定义

在开始编码之前，我们需要一个清晰的“作战计划”。首先，确保你的Python环境已经就绪。我们主要会用到numpy进行数值计算，matplotlib进行结果可视化。如果你使用Anaconda，这些库通常已经预装。如果没有，可以通过以下命令快速安装：

pip install numpy matplotlib

接下来，我们需要一个具体的“战场”来测试我们的算法。这里我们选择ZDT1问题，它是多目标优化领域一个非常经典的基准测试函数。为什么选它？因为它结构清晰，Pareto前沿（所有最优解的集合在目标空间中的形状）是已知的、简单的凸曲线，非常适合用来验证算法的正确性。

ZDT1问题定义如下：

• 决策变量：x 是一个30维的向量，每个分量都在[0, 1]区间内。
• 目标函数1：f1(x) = x1。我们希望它越小越好。
• 目标函数2：g(x) = 1 + 9 * (sum(x2, x3, ..., x30) / 29)，然后 f2(x) = g(x) * (1 - sqrt(x1 / g(x)))。同样，我们希望f2越小越好。

这个问题的Pareto前沿是当g(x)=1时，即x2到x30全部为0，此时f2 = 1 - sqrt(f1)，是一条从(0,1)到(1,0)的光滑曲线。

为了在代码中清晰地表达这个问题，我们创建一个独立的函数文件 problems.py：

# problems.py
import numpy as np

def ZDT1(individual):
    """
    计算ZDT1测试问题的两个目标函数值。
    参数:
        individual: 一个一维numpy数组，长度为30，代表一个解。
    返回:
        一个包含两个目标函数值的元组 (f1, f2)。
    """
    num_variables = len(individual)
    f1 = individual[0]  # 第一个目标就是第一个决策变量
    g = 1.0 + 9.0 * np.sum(individual[1:]) / (num_variables - 1)
    f2 = g * (1.0 - np.sqrt(f1 / g))
    return np.array([f1, f2])

这个函数将成为我们算法评估每一个解的“裁判”。有了明确的问题，我们就可以开始构建算法的主体了。

2. NSGA-II核心组件实现

NSGA-II的精妙之处在于其快速非支配排序和拥挤度距离计算。这两个机制共同作用，引导种群在进化的过程中，既朝着真正的Pareto前沿收敛（收敛性），又能让解在整个前沿上均匀散布（分布性）。

2.1 快速非支配排序

“支配”是多目标优化中最核心的概念。简单来说，解A支配解B，意味着A在所有目标上都不比B差，并且至少在一个目标上严格比B好。非支配排序就是根据这个关系，把种群里的所有个体分成不同的层级（Front）。第一层（Front 0）是那些不被任何其他个体支配的“精英”，也就是当前找到的Pareto最优解；第二层是被第一层个体支配，但不受其他层支配的解，依此类推。

实现快速非支配排序的关键是避免两两比较的O(MN³)复杂度。下面是一个高效的Python实现：

# nsga2_core.py
import numpy as np

def fast_non_dominated_sort(population, fitness):
    """
    对种群进行快速非支配排序。
    参数:
        population: 种群，一个二维numpy数组，每行是一个个体。
        fitness: 适应度矩阵，一个二维numpy数组，fitness[i]对应population[i]的目标函数值。
    返回:
        fronts: 一个列表的列表，fronts[k]包含了第k层非支配层中个体的索引。
    """
    num_individuals = len(population)
    # S[p]: 被个体p支配的个体集合
    S = [[] for _ in range(num_individuals)]
    # n[p]: 支配个体p的个体数量
    n = np.zeros(num_individuals, dtype=int)
    # rank[p]: 个体p所在的层级
    rank = np.zeros(num_individuals, dtype=int)

    fronts = [[]]  # 初始化第一层

    # 第一遍遍历：计算支配关系
    for p in range(num_individuals):
        S[p] = []
        n[p] = 0
        for q in range(num_individuals):
            if p == q:
                continue
            # 判断p是否支配q
            if dominates(fitness[p], fitness[q]):
                S[p].append(q)
            elif dominates(fitness[q], fitness[p]):
                n[p] += 1
        if n[p] == 0:  # 没有个体支配p，p属于第一层
            rank[p] = 0
            fronts[0].append(p)

    # 迭代构建后续层级
    i = 0
    while fronts[i]:  # 当第i层不为空时
        next_front = []
        for p in fronts[i]:
            for q in S[p]:
                n[q] -= 1
                if n[q] == 0:
                    rank[q] = i + 1
                    next_front.append(q)
        i += 1
        if next_front:
            fronts.append(next_front)
        else:
            break
    return fronts, rank

def dominates(a, b):
    """
    判断解a是否支配解b。
    规则：对于所有目标，a <= b，且至少存在一个目标，a < b。
    """
    # 检查a是否在所有目标上都不差于b
    not_worse = np.all(a <= b)
    # 检查a是否在至少一个目标上严格优于b
    strictly_better = np.any(a < b)
    return not_worse and strictly_better

提示：dominates函数的实现是算法的基石。这里我们假设所有目标都是最小化问题。如果你的问题是最大化某些目标，需要先将其转化为最小化（例如乘以-1）。

2.2 拥挤度距离计算

仅仅分层还不够。在同一层（比如第一层）里可能有成百上千个解，我们如何从中选出最有“代表性”的、分布最均匀的那些进入下一代呢？这就需要拥挤度距离。它衡量的是一个解在目标空间中，与相邻解之间的“孤独程度”。距离越大，说明它周围越“空旷”，保留它对维持种群多样性越有利。

计算拥挤度距离的步骤是：对于同一非支配层内的所有个体，针对每一个目标函数，将它们按该目标的值从小到大排序。边界上的两个个体（最大值和最小值）被赋予无穷大的拥挤度距离（确保它们被优先保留）。对于中间的个体，其拥挤度距离是相邻个体在该目标上函数值之差的归一化累加和。

# nsga2_core.py
def crowding_distance_assignment(front, fitness):
    """
    计算指定非支配层中所有个体的拥挤度距离。
    参数:
        front: 一个列表，包含属于同一非支配层的个体索引。
        fitness: 适应度矩阵。
    返回:
        distances: 一个一维数组，distances[i]对应front[i]的拥挤度距离。
    """
    num_individuals = len(front)
    num_objectives = fitness.shape[1]
    distances = np.zeros(num_individuals)

    if num_individuals <= 2:
        # 如果层内个体数小于等于2，直接赋予最大距离
        distances[:] = np.inf
        return distances

    # 对每个目标分别计算
    for m in range(num_objectives):
        # 获取当前目标下，该层所有个体的适应度值
        obj_values = [fitness[i][m] for i in front]
        # 按该目标值排序，得到排序后的索引
        sorted_indices = np.argsort(obj_values)

        # 边界个体距离设为无穷大
        distances[sorted_indices[0]] = np.inf
        distances[sorted_indices[-1]] = np.inf

        # 获取该目标的最大最小值，用于归一化
        f_max = obj_values[sorted_indices[-1]]
        f_min = obj_values[sorted_indices[0]]
        norm = f_max - f_min
        if norm < 1e-10:  # 防止除零
            norm = 1.0

        # 计算中间个体的拥挤度距离
        for i in range(1, num_individuals - 1):
            idx = sorted_indices[i]
            idx_next = sorted_indices[i + 1]
            idx_prev = sorted_indices[i - 1]
            distances[idx] += (obj_values[idx_next] - obj_values[idx_prev]) / norm

    return distances

有了这两个核心武器，我们就可以构建NSGA-II的选择机制了：在选择个体进入下一代时，优先选择非支配层级更靠前（rank值更小）的个体；如果层级相同，则优先选择拥挤度距离更大的个体。这种机制完美平衡了“趋优”（收敛）和“扩散”（多样性）两个进化方向。

3. 遗传算子与主循环构建

遗传算法模仿生物进化，离不开选择、交叉、变异这三个基本操作。NSGA-II的进化流程可以概括为以下几步：

1. 初始化：随机生成一个大小为N的父代种群。
2. 评估：计算种群中每个个体的目标函数值（适应度）。
3. 主循环（迭代）： a. 选择：通过二元锦标赛选择（基于非支配层级和拥挤度距离）从父代中选出N个个体形成交配池。 b. 交叉与变异：对交配池中的个体进行交叉和变异操作，生成同样大小为N的子代种群。 c. 合并：将父代种群和子代种群合并，形成一个大小为2N的临时种群。 d. 环境选择：对合并后的2N个个体进行快速非支配排序和拥挤度距离计算。然后，按照层级从低到高、同层内按拥挤度从大到小的顺序，选出前N个个体作为新的父代。
4. 终止：重复步骤3，直到达到预设的最大迭代次数。

下面，我们用代码将这个流程具体化。首先定义遗传算子：

# operators.py
import numpy as np

def simulated_binary_crossover(parent1, parent2, eta=20):
    """
    模拟二进制交叉（SBX）。
    参数:
        parent1, parent2: 父代个体。
        eta: 分布指数，控制子代与父代的相似程度，值越大子代越接近父代。
    返回:
        child1, child2: 两个子代个体。
    """
    # 确保父代是numpy数组
    parent1 = np.array(parent1)
    parent2 = np.array(parent2)
    u = np.random.rand(len(parent1))
    beta = np.empty_like(u)
    
    # SBX的核心公式
    mask = u <= 0.5
    beta[mask] = (2 * u[mask]) ** (1.0 / (eta + 1))
    beta[~mask] = (1.0 / (2 * (1 - u[~mask]))) ** (1.0 / (eta + 1))
    
    child1 = 0.5 * ((1 + beta) * parent1 + (1 - beta) * parent2)
    child2 = 0.5 * ((1 - beta) * parent1 + (1 + beta) * parent2)
    
    # 确保子代在变量边界内（对于ZDT1是[0,1]）
    child1 = np.clip(child1, 0, 1)
    child2 = np.clip(child2, 0, 1)
    return child1, child2

def polynomial_mutation(individual, eta=20, mutation_prob=0.1):
    """
    多项式变异。
    参数:
        individual: 待变异的个体。
        eta: 分布指数。
        mutation_prob: 每个基因发生变异的概率。
    返回:
        mutated_individual: 变异后的个体。
    """
    mutated = np.array(individual, copy=True)
    for i in range(len(mutated)):
        if np.random.rand() < mutation_prob:
            u = np.random.rand()
            delta = 0.0
            if u < 0.5:
                delta = (2 * u) ** (1.0 / (eta + 1)) - 1
            else:
                delta = 1 - (2 * (1 - u)) ** (1.0 / (eta + 1))
            mutated[i] += delta
            # 边界处理
            if mutated[i] < 0:
                mutated[i] = 0
            elif mutated[i] > 1:
                mutated[i] = 1
    return mutated

现在，我们可以组装完整的主循环了。为了便于管理和测试不同参数，我们将算法封装成一个类：

# nsga2.py
import numpy as np
from .nsga2_core import fast_non_dominated_sort, crowding_distance_assignment
from .operators import simulated_binary_crossover, polynomial_mutation

class NSGA2:
    def __init__(self, problem_func, num_variables=30, pop_size=100, 
                 crossover_prob=0.9, mutation_prob=0.1, eta_crossover=20, eta_mutation=20, max_generations=250):
        """
        初始化NSGA-II算法。
        参数:
            problem_func: 问题函数，输入个体，输出目标值数组。
            num_variables: 决策变量维度。
            pop_size: 种群大小。
            crossover_prob: 交叉概率。
            mutation_prob: 变异概率。
            eta_crossover: 交叉分布指数。
            eta_mutation: 变异分布指数。
            max_generations: 最大进化代数。
        """
        self.problem = problem_func
        self.num_vars = num_variables
        self.pop_size = pop_size
        self.cx_prob = crossover_prob
        self.mut_prob = mutation_prob
        self.eta_cx = eta_crossover
        self.eta_mut = eta_mutation
        self.max_gen = max_generations
        
        # 边界（ZDT1问题变量范围是[0,1]）
        self.lower_bound = np.zeros(self.num_vars)
        self.upper_bound = np.ones(self.num_vars)
        
        # 记录进化过程
        self.history = {'population': [], 'fitness': []}
        
    def _initialize_population(self):
        """随机初始化种群。"""
        return np.random.rand(self.pop_size, self.num_vars)
    
    def _evaluate_population(self, population):
        """评估整个种群的目标函数值。"""
        fitness = np.array([self.problem(ind) for ind in population])
        return fitness
    
    def _selection(self, population, fitness, fronts, crowding_dist):
        """
        二元锦标赛选择，用于生成交配池。
        选择标准：先看非支配层级（rank），层级低者胜出；层级相同则看拥挤度距离，距离大者胜出。
        """
        mating_pool = []
        for _ in range(self.pop_size):
            # 随机选择两个参赛者
            idx1, idx2 = np.random.choice(len(population), 2, replace=False)
            # 比较
            if fronts[idx1] < fronts[idx2]:
                winner = idx1
            elif fronts[idx1] > fronts[idx2]:
                winner = idx2
            else:  # 同一层级，比较拥挤度距离
                if crowding_dist[idx1] > crowding_dist[idx2]:
                    winner = idx1
                else:
                    winner = idx2
            mating_pool.append(population[winner])
        return np.array(mating_pool)
    
    def run(self):
        """执行NSGA-II主循环。"""
        # 1. 初始化
        population = self._initialize_population()
        fitness = self._evaluate_population(population)
        
        for gen in range(self.max_gen):
            # 2. 非支配排序与拥挤度计算（针对父代）
            fronts_list, rank = fast_non_dominated_sort(population, fitness)
            # 将fronts列表转换为每个个体所属层级的数组
            fronts = np.zeros(len(population), dtype=int)
            for level, front in enumerate(fronts_list):
                for idx in front:
                    fronts[idx] = level
            
            crowding_dist = np.zeros(len(population))
            for front in fronts_list:
                if len(front) > 0:
                    front_dist = crowding_distance_assignment(front, fitness)
                    for i, idx in enumerate(front):
                        crowding_dist[idx] = front_dist[i]
            
            # 记录当前代信息（可选，用于分析）
            self.history['population'].append(population.copy())
            self.history['fitness'].append(fitness.copy())
            
            # 3. 选择（生成交配池）
            mating_pool = self._selection(population, fitness, fronts, crowding_dist)
            
            # 4. 交叉与变异（生成子代）
            offspring = []
            for i in range(0, self.pop_size, 2):
                parent1 = mating_pool[i]
                parent2 = mating_pool[i+1] if i+1 < self.pop_size else mating_pool[0]
                
                # 交叉
                if np.random.rand() < self.cx_prob:
                    child1, child2 = simulated_binary_crossover(parent1, parent2, self.eta_cx)
                else:
                    child1, child2 = parent1.copy(), parent2.copy()
                
                # 变异
                child1 = polynomial_mutation(child1, self.eta_mut, self.mut_prob)
                child2 = polynomial_mutation(child2, self.eta_mut, self.mut_prob)
                
                offspring.extend([child1, child2])
            offspring = np.array(offspring[:self.pop_size])  # 确保子代数量与种群大小一致
            
            # 5. 评估子代
            offspring_fitness = self._evaluate_population(offspring)
            
            # 6. 合并父代与子代
            combined_pop = np.vstack([population, offspring])
            combined_fit = np.vstack([fitness, offspring_fitness])
            
            # 7. 环境选择（从合并的2N个个体中选出N个作为新父代）
            new_fronts_list, _ = fast_non_dominated_sort(combined_pop, combined_fit)
            
            next_population = []
            next_fitness = []
            remaining = self.pop_size
            current_front_idx = 0
            
            while remaining > 0 and current_front_idx < len(new_fronts_list):
                current_front = new_fronts_list[current_front_idx]
                if len(current_front) <= remaining:
                    # 如果当前层所有个体都能加入，则全部加入
                    for idx in current_front:
                        next_population.append(combined_pop[idx])
                        next_fitness.append(combined_fit[idx])
                    remaining -= len(current_front)
                else:
                    # 如果当前层个体数超过剩余名额，则根据拥挤度距离筛选
                    front_fitness = combined_fit[current_front]
                    crowding_dist_front = crowding_distance_assignment(current_front, combined_fit)
                    # 按拥挤度距离降序排序当前层的个体索引
                    sorted_by_crowding = [current_front[i] for i in np.argsort(-crowding_dist_front)]
                    # 选取前remaining个
                    for i in range(remaining):
                        idx = sorted_by_crowding[i]
                        next_population.append(combined_pop[idx])
                        next_fitness.append(combined_fit[idx])
                    remaining = 0
                current_front_idx += 1
            
            population = np.array(next_population)
            fitness = np.array(next_fitness)
            
            # 打印进度
            if (gen + 1) % 50 == 0:
                print(f"Generation {gen+1}/{self.max_gen} completed.")
        
        # 最终的非支配排序，获取第一层（Pareto前沿）
        final_fronts_list, _ = fast_non_dominated_sort(population, fitness)
        pareto_front_indices = final_fronts_list[0]
        pareto_population = population[pareto_front_indices]
        pareto_fitness = fitness[pareto_front_indices]
        
        return pareto_population, pareto_fitness, self.history

代码看起来有点长，但逻辑是清晰的。关键在于第7步的环境选择，它完美体现了NSGA-II的精英策略：父代和子代同台竞技，只有最优秀的个体（层级低且分布好）才能存活下来。这种机制保证了好的解不会被轻易丢弃，加速了收敛。

4. 运行实验与结果可视化

理论说得再多，不如跑一遍看看效果。我们写一个主程序来调用上面实现的NSGA-II类，解决ZDT1问题，并将结果可视化。

# main.py
import numpy as np
import matplotlib.pyplot as plt
from problems import ZDT1
from nsga2 import NSGA2

def main():
    # 设置随机种子，确保结果可复现
    np.random.seed(42)
    
    # 初始化NSGA-II求解器
    solver = NSGA2(
        problem_func=ZDT1,
        num_variables=30,
        pop_size=100,
        crossover_prob=0.9,
        mutation_prob=0.1,
        eta_crossover=20,
        eta_mutation=20,
        max_generations=250
    )
    
    print("开始运行NSGA-II求解ZDT1问题...")
    pareto_solutions, pareto_front, history = solver.run()
    print(f"优化完成！找到 {len(pareto_front)} 个Pareto最优解。")
    
    # 绘制最终的Pareto前沿
    plt.figure(figsize=(10, 6))
    # 绘制找到的Pareto前沿点
    plt.scatter(pareto_front[:, 0], pareto_front[:, 1], c='red', s=50, edgecolors='k', label='NSGA-II Pareto Front', zorder=5)
    
    # 绘制真实的Pareto前沿（用于对比）
    true_f1 = np.linspace(0, 1, 300)
    true_f2 = 1 - np.sqrt(true_f1)
    plt.plot(true_f1, true_f2, 'b--', linewidth=2, label='True Pareto Front')
    
    plt.xlabel('Objective 1 (f1)', fontsize=12)
    plt.ylabel('Objective 2 (f2)', fontsize=12)
    plt.title('NSGA-II on ZDT1 Problem', fontsize=14)
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()
    
    # 可选：绘制进化过程中种群分布的变化（动态图或最后几代的静态图）
    # 例如，查看初始种群和最终种群的对比
    initial_pop = history['population'][0]
    initial_fit = history['fitness'][0]
    final_pop = history['population'][-1]
    final_fit = history['fitness'][-1]
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    axes[0].scatter(initial_fit[:, 0], initial_fit[:, 1], c='gray', alpha=0.6, s=30, label='Initial Population')
    axes[0].plot(true_f1, true_f2, 'b--', linewidth=2, label='True Pareto Front')
    axes[0].set_xlabel('Objective 1 (f1)')
    axes[0].set_ylabel('Objective 2 (f2)')
    axes[0].set_title('Generation 0')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    axes[1].scatter(final_fit[:, 0], final_fit[:, 1], c='green', alpha=0.6, s=30, label='Final Population')
    axes[1].plot(true_f1, true_f2, 'b--', linewidth=2, label='True Pareto Front')
    axes[1].set_xlabel('Objective 1 (f1)')
    axes[1].set_ylabel('Objective 2 (f2)')
    axes[1].set_title(f'Generation {solver.max_gen}')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()

运行这段代码，你会看到两张图。第一张图展示了算法找到的Pareto前沿（红点）与真实前沿（蓝色虚线）的对比。理想情况下，红点应该紧密地分布在蓝色虚线附近。第二张图对比了初始随机种群和最终种群的分布，你可以清晰地看到进化过程如何将散乱的点驱赶到最优前沿附近。

5. 关键参数调优与影响分析

NSGA-II的性能很大程度上依赖于其参数设置。盲目使用默认值可能无法得到理想的结果。下面我们通过一系列对比实验，来看看几个核心参数如何影响算法的表现。我们将主要调整种群大小（pop_size）、交叉分布指数（eta_crossover）和变异概率（mutation_prob）。

为了系统地进行比较，我们设计一个实验，固定其他参数，只改变其中一个，观察最终Pareto前沿的收敛性（离真实前沿有多近）和分布性（解是否均匀覆盖整个前沿）。我们可以用两个简单的指标来量化：

• 收敛性指标（GD）：计算找到的解到真实Pareto前沿的最小距离的平均值。值越小越好。
• 分布性指标（Spacing）：计算找到的解彼此之间的距离标准差。值越小，说明解分布越均匀。

我们修改一下主程序，进行批量实验并汇总结果：

# parameter_study.py
import numpy as np
import matplotlib.pyplot as plt
from problems import ZDT1
from nsga2 import NSGA2

def calculate_gd(pf_found, pf_true):
    """计算世代距离（Generational Distance），衡量收敛性。"""
    min_distances = []
    for point in pf_found:
        # 计算该点到真实前沿上所有点的欧氏距离，取最小值
        distances = np.sqrt(np.sum((pf_true - point) ** 2, axis=1))
        min_distances.append(np.min(distances))
    return np.mean(min_distances)

def calculate_spacing(pf):
    """计算间距（Spacing），衡量分布均匀性。"""
    if len(pf) < 2:
        return 0.0
    distances = []
    for i in range(len(pf)):
        # 计算该点到其他所有点的最小距离（除了自己）
        d = np.sqrt(np.sum((pf - pf[i]) ** 2, axis=1))
        d[i] = np.inf  # 排除自身
        distances.append(np.min(d))
    d_mean = np.mean(distances)
    spacing = np.sqrt(np.sum((np.array(distances) - d_mean) ** 2) / (len(pf) - 1))
    return spacing

# 定义参数研究范围
param_studies = {
    'pop_size': [50, 100, 200],
    'eta_crossover': [5, 20, 50],
    'mutation_prob': [0.01, 0.1, 0.3]
}

# 生成真实的Pareto前沿参考点（用于计算GD）
true_f1 = np.linspace(0, 1, 500)
true_f2 = 1 - np.sqrt(true_f1)
pf_true = np.column_stack([true_f1, true_f2])

results = {}
default_params = {
    'num_variables': 30,
    'crossover_prob': 0.9,
    'max_generations': 250
}

np.random.seed(42)  # 固定种子，使不同实验的初始种群一致

for param_name, param_values in param_studies.items():
    print(f"\n=== 研究参数: {param_name} ===")
    gd_list = []
    spacing_list = []
    
    for value in param_values:
        print(f"  正在运行 {param_name} = {value}...")
        kwargs = default_params.copy()
        kwargs[param_name] = value
        # 设置其他参数为默认值
        if param_name != 'pop_size':
            kwargs['pop_size'] = 100
        if param_name != 'eta_crossover':
            kwargs['eta_crossover'] = 20
        if param_name != 'mutation_prob':
            kwargs['mutation_prob'] = 0.1
            
        solver = NSGA2(problem_func=ZDT1, **kwargs)
        _, pareto_front, _ = solver.run()
        
        gd = calculate_gd(pareto_front, pf_true)
        spacing = calculate_spacing(pareto_front)
        
        gd_list.append(gd)
        spacing_list.append(spacing)
        print(f"    世代距离(GD): {gd:.6f}, 间距(Spacing): {spacing:.6f}")
    
    results[param_name] = {
        'values': param_values,
        'GD': gd_list,
        'Spacing': spacing_list
    }

# 可视化结果
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
param_names = list(param_studies.keys())

for idx, param in enumerate(param_names):
    ax = axes[idx]
    values = results[param]['values']
    gd = results[param]['GD']
    spacing = results[param]['Spacing']
    
    ax.plot(values, gd, 'o-', linewidth=2, markersize=8, label='Generational Distance (GD)', color='tab:blue')
    ax.set_xlabel(param, fontsize=11)
    ax.set_ylabel('GD (lower is better)', color='tab:blue', fontsize=11)
    ax.tick_params(axis='y', labelcolor='tab:blue')
    ax.grid(True, alpha=0.3)
    
    ax2 = ax.twinx()
    ax2.plot(values, spacing, 's--', linewidth=2, markersize=8, label='Spacing (lower is better)', color='tab:red')
    ax2.set_ylabel('Spacing', color='tab:red', fontsize=11)
    ax2.tick_params(axis='y', labelcolor='tab:red')
    
    # 合并图例
    lines1, labels1 = ax.get_legend_handles_labels()
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax.legend(lines1 + lines2, labels1 + labels2, loc='upper center')
    ax.set_title(f'Impact of {param}')

plt.tight_layout()
plt.show()

运行这个参数研究脚本，你会得到三张趋势图。从这些图中，我们可以总结出一些实用的调参经验：

在实际项目中，我通常的做法是：先用中等规模的种群（如100）和默认参数（交叉概率0.9，变异概率0.1，分布指数20）跑一遍，看看大致结果。如果发现收敛性不好（解离真实前沿远），可以尝试增大种群规模或代数。如果发现解都挤在前沿的某一段（分布性差），可以适当提高变异概率或降低交叉分布指数，增加探索能力。这是一个需要根据具体问题反复微调的过程。

6. 进阶：从ZDT1到真实世界问题

掌握了在ZDT1这个“操场”上的训练后，是时候将NSGA-II应用到更真实的场景了。其核心框架不变，你需要调整的只是问题定义部分。我们来看一个简化版的工程案例：天线阵列设计优化。

假设我们需要设计一个由N个阵元组成的线性天线阵列。我们希望：

1. 目标1：主瓣宽度尽可能窄（提高方向性）。
2. 目标2：最高旁瓣电平尽可能低（减少干扰）。
3. 决策变量：每个阵元的激励幅度（Amplitude）和相位（Phase）。

这是一个典型的多目标优化问题，两个目标相互冲突（通常，更窄的主瓣会导致更高的旁瓣）。我们可以这样定义问题函数：

# real_world_problem.py
import numpy as np

def antenna_array_design(individual, num_elements=10, wavelength=1.0, element_spacing=0.5):
    """
    简化天线阵列设计问题。
    个体编码：[A1, A2, ..., AN, Phi1, Phi2, ..., PhiN]，前半部分是幅度，后半部分是相位。
    目标1：主瓣宽度（Beamwidth）。
    目标2：最大旁瓣电平（Sidelobe Level, SLL）。
    """
    N = num_elements
    amplitudes = individual[:N]  # 幅度，假设范围[0, 1]
    phases = individual[N:]      # 相位，假设范围[-pi, pi]
    
    # 计算阵列因子（Array Factor, AF）
    theta = np.linspace(-np.pi, np.pi, 361)  # 角度范围
    AF = np.zeros_like(theta, dtype=complex)
    for n in range(N):
        AF += amplitudes[n] * np.exp(1j * (phases[n] + 2 * np.pi * element_spacing * n * np.sin(theta) / wavelength))
    AF_magnitude = np.abs(AF)
    AF_normalized = AF_magnitude / np.max(AF_magnitude)  # 归一化
    AF_db = 20 * np.log10(AF_normalized + 1e-10)  # 转换为分贝
    
    # 目标1：主瓣宽度（-3dB点之间的宽度）
    main_lobe_mask = (theta >= -np.pi/18) & (theta <= np.pi/18)  # 假设主瓣在±10度内寻找
    main_lobe_pattern = AF_db[main_lobe_mask]
    peak_index = np.argmax(main_lobe_pattern)
    peak_power = main_lobe_pattern[peak_index]
    half_power = peak_power - 3  # -3dB点
    # 简化计算：找到主瓣两侧-3dB点的角度差（这里需要更精细的插值，仅为示例）
    # 此处简化为主瓣内低于-3dB的点的角度范围
    below_3db = main_lobe_pattern < half_power
    if np.any(below_3db):
        beamwidth = np.sum(below_3db) * (theta[1] - theta[0]) * 180 / np.pi  # 转换为度
    else:
        beamwidth = 10.0  # 默认值
    
    # 目标2：最大旁瓣电平（排除主瓣区域后的最大值）
    sidelobe_mask = ~main_lobe_mask
    max_sidelobe = np.max(AF_db[sidelobe_mask]) if np.any(sidelobe_mask) else -30
    
    # 我们希望主瓣宽度越小越好，旁瓣电平越低越好（所以取负值，转化为最小化问题）
    f1 = beamwidth
    f2 = -max_sidelobe  # 因为max_sidelobe是负的分贝值，取负后变成正数，最小化它相当于最大化旁瓣抑制
    
    return np.array([f1, f2])

将这个函数替换掉之前ZDT1，作为NSGA2类的problem_func参数，并相应地调整决策变量的边界（幅度[0,1]，相位[-π, π]）和维度（2*N），你就可以用同样的NSGA-II框架来优化天线设计了。最终得到的Pareto前沿会展示一系列“主瓣宽度-旁瓣电平”的权衡方案，供射频工程师根据实际系统要求进行选择。

这个例子展示了NSGA-II的强大通用性。无论是机械结构设计（轻量化 vs 强度）、投资组合优化（收益 vs 风险）、还是机器学习超参数调优（精度 vs 模型复杂度），只要你能将问题建模为多个相互冲突的目标，并定义好决策变量和评估函数，NSGA-II就能为你探索那片广阔的权衡空间。

代码写到这里，我们已经完成了一个从理论到实践、从标准测试函数到工程问题的完整闭环。过程中最深的体会是，理解算法每个步骤背后的意图，远比死记硬背公式更重要。比如，拥挤度距离的计算，本质上是在同一优先级的解中，奖励那些“特立独行”的个体，避免种群扎堆，这是维持多样性的关键。而快速非支配排序中的精英保留策略，则确保了搜索不会“开倒车”，好的发现不会被遗忘。

在实际应用中，你可能会遇到计算耗时太长的问题。这时可以考虑对目标函数评估进行向量化加速，或者对非支配排序等核心步骤尝试用numba进行JIT编译。如果问题变量维度极高（成百上千维），标准的NSGA-II可能会遇到“维度灾难”，这时可以关注其改进版本NSGA-III，它通过引入参考点来更好地处理高维目标空间。

多目标优化没有唯一的“正确答案”，它提供的是一组“可能的最佳答案”。NSGA-II的价值，就在于高效、自动地为你勾勒出这片可能性的疆域，将复杂的多目标决策，转化为一张清晰的可视化地图。剩下的，就是结合你的领域知识，在那条Pareto前沿上，做出最适合当前情境的抉择。