基于NSGA-II实现DLTZ1帕累托最优前沿收敛的不稳定问题排查

问题概况 #

• 实现的NSGA-II算法逻辑自认为正确，但稳定性极差
• 针对DLTZ1测试函数：
  • 使用标准g(x)计算逻辑时，仅约1/10的运行能收敛到预期的3D帕累托前沿
  • 手动固定g(x)=0.5时，算法表现完全正常
• 算法在原论文的基础测试函数上可得到预期结果
• 提升种群规模至300可解决收敛问题，但当前使用的是文献推荐的参数

DLTZ1公式 #

DLTZ1的3目标场景定义（$n$为变量总数，$k=2$）：
$$
\begin{align*}
g(x_{m:}) &= 100\left( (n-m+1) + \sum_{i=m}^{n} (x_i - 0.5)^2 - \cos(20\pi(x_i - 0.5)) \right) \
f_1(x) &= \frac{1}{2}x_1(1+g(x_{m:})) \
f_2(x) &= \frac{1}{2}(1-x_1)x_2(1+g(x_{m:})) \
f_3(x) &= \frac{1}{2}(1-x_1)(1-x_2)(1+g(x_{m:}))
\end{align*}
$$

可能的原因分析 #

1. g(x)的数值稳定性问题：
   DLTZ1的g(x)依赖后$n-k$个变量的求和计算，若初始种群在这些变量上分布极端（比如大量个体的$x_i$接近0或1），会导致g(x)值异常偏大，直接扭曲目标函数的比例关系，使算法搜索方向偏离帕累托前沿。手动固定g(x)=0.5时，目标函数比例关系稳定，因此算法能正常收敛。

2. 遗传操作的适配性不足：
   若交叉/变异算子在处理DLTZ1变量空间时无法维持种群多样性，会导致大部分运行陷入局部最优。比如实数编码的变异步长过大，可能使g(x)突变幅度过大；交叉概率过低则无法充分交换优秀个体的基因信息。

3. 非支配排序/拥挤度计算的细节偏差：
   3目标场景下的非支配排序或拥挤度计算若存在精度问题（比如浮点数比较阈值设置不当），可能导致真正的帕累托前沿个体被错误淘汰，进而影响算法收敛稳定性。

调试建议 #

• 检查g(x)计算与种群初始化：
  输出初始种群及每一代的g(x)值，确认是否存在大量异常值。若有，可将后$n-k$个变量的初始化范围限制在[0.4, 0.6]（缩小初始波动），观察稳定性是否提升。
• 微调遗传操作参数：
  • 尝试将SBX交叉概率从默认的0.9提升至0.95，增大优秀基因的交换概率
  • 调整变异步长（比如将多项式变异的分布指数从20调至10），降低g(x)的突变幅度
• 验证核心模块正确性：
  拿少量已知的帕累托前沿个体代入非支配排序和拥挤度计算模块，对比标准实现的结果，排查逻辑错误
• 临时调整精英保留策略：
  在选择操作中优先保留拥挤度高的个体，确保帕累托前沿附近的个体不被过早淘汰

附可运行的NSGA-II实现代码 #

import numpy as np

def dltz1(x, m=3):
    n = len(x)
    k = n - m + 1
    x_m = x[m-1:]
    g = 100 * (k + np.sum((x_m - 0.5)**2 - np.cos(20 * np.pi * (x_m - 0.5))))
    f = np.zeros(m)
    f[0] = 0.5 * x[0] * (1 + g)
    for i in range(1, m-1):
        f[i] = 0.5 * np.prod(x[:i]) * (1 - x[i]) * (1 + g)
    f[-1] = 0.5 * np.prod(x[:m-1]) * (1 + g)
    return f

class NSGAII:
    def __init__(self, pop_size=100, n_vars=10, m_obj=3, max_gen=200):
        self.pop_size = pop_size
        self.n_vars = n_vars
        self.m_obj = m_obj
        self.max_gen = max_gen
        self.pop = np.random.uniform(0, 1, (pop_size, n_vars))
    
    def non_dominated_sort(self, fitness):
        # 非支配排序核心逻辑
        pop_size = fitness.shape[0]
        dominate_count = np.zeros(pop_size)
        dominated_set = [[] for _ in range(pop_size)]
        fronts = [[]]
        
        for p in range(pop_size):
            for q in range(pop_size):
                if p == q:
                    continue
                # 判断p是否支配q
                p_dom_q = all(fitness[p] <= fitness[q]) and any(fitness[p] < fitness[q])
                q_dom_p = all(fitness[q] <= fitness[p]) and any(fitness[q] < fitness[p])
                if p_dom_q:
                    dominated_set[p].append(q)
                elif q_dom_p:
                    dominate_count[p] += 1
            if dominate_count[p] == 0:
                fronts[0].append(p)
        
        i = 0
        while len(fronts[i]) > 0:
            next_front = []
            for p in fronts[i]:
                for q in dominated_set[p]:
                    dominate_count[q] -= 1
                    if dominate_count[q] == 0:
                        next_front.append(q)
            i += 1
            fronts.append(next_front)
        return fronts[:-1]
    
    def crowding_distance(self, front, fitness):
        distance = np.zeros(len(front))
        for obj in range(self.m_obj):
            sorted_indices = sorted(range(len(front)), key=lambda x: fitness[front[x]][obj])
            distance[sorted_indices[0]] = float('inf')
            distance[sorted_indices[-1]] = float('inf')
            obj_min = fitness[front[sorted_indices[0]]][obj]
            obj_max = fitness[front[sorted_indices[-1]]][obj]
            if obj_max - obj_min == 0:
                continue
            for i in range(1, len(sorted_indices)-1):
                distance[sorted_indices[i]] += (fitness[front[sorted_indices[i+1]]][obj] - fitness[front[sorted_indices[i-1]]][obj]) / (obj_max - obj_min)
        return distance
    
    def sbx_crossover(self, parent1, parent2, eta=20):
        child1 = np.copy(parent1)
        child2 = np.copy(parent2)
        for i in range(self.n_vars):
            if np.random.rand() < 0.9:
                if abs(parent1[i] - parent2[i]) > 1e-10:
                    x1 = min(parent1[i], parent2[i])
                    x2 = max(parent1[i], parent2[i])
                    rand = np.random.rand()
                    beta = 1 + (2 * (x1 - 0) / (x2 - x1)) if rand <= 0.5 else 1 + (2 * (1 - x2) / (x2 - x1))
                    beta = beta ** -(eta + 1)
                    c1 = 0.5 * ((x1 + x2) - beta * (x2 - x1))
                    c2 = 0.5 * ((x1 + x2) + beta * (x2 - x1))
                    child1[i] = np.clip(c1, 0, 1)
                    child2[i] = np.clip(c2, 0, 1)
        return child1, child2
    
    def polynomial_mutation(self, individual, eta=20):
        mutated = np.copy(individual)
        for i in range(self.n_vars):
            if np.random.rand() < 1/self.n_vars:
                rand = np.random.rand()
                if rand <= 0.5:
                    delta = (2 * rand) ** (1/(eta+1)) - 1
                else:
                    delta = 1 - (2 * (1 - rand)) ** (1/(eta+1))
                mutated[i] += delta
                mutated[i] = np.clip(mutated[i], 0, 1)
        return mutated
    
    def evolve(self):
        for gen in range(self.max_gen):
            # 计算适应度
            fitness = np.array([dltz1(ind) for ind in self.pop])
            # 非支配排序
            fronts = self.non_dominated_sort(fitness)
            # 生成子代
            offspring = []
            while len(offspring) < self.pop_size:
                p1_idx = np.random.randint(self.pop_size)
                p2_idx = np.random.randint(self.pop_size)
                while p2_idx == p1_idx:
                    p2_idx = np.random.randint(self.pop_size)
                child1, child2 = self.sbx_crossover(self.pop[p1_idx], self.pop[p2_idx])
                child1 = self.polynomial_mutation(child1)
                child2 = self.polynomial_mutation(child2)
                offspring.append(child1)
                offspring.append(child2)
            offspring = np.array(offspring[:self.pop_size])
            # 合并种群
            combined_pop = np.vstack((self.pop, offspring))
            combined_fitness = np.vstack((fitness, np.array([dltz1(ind) for ind in offspring])))
            # 重新排序
            combined_fronts = self.non_dominated_sort(combined_fitness)
            # 选择下一代种群
            new_pop = []
            for front in combined_fronts:
                if len(new_pop) + len(front) <= self.pop_size:
                    new_pop.extend(front)
                else:
                    distances = self.crowding_distance(front, combined_fitness)
                    sorted_front = [x for _,x in sorted(zip(distances, front), key=lambda x: -x[0])]
                    new_pop.extend(sorted_front[:self.pop_size - len(new_pop)])
                    break
            self.pop = combined_pop[new_pop]
        # 返回最终种群和适应度
        final_fitness = np.array([dltz1(ind) for ind in self.pop])
        return self.pop, final_fitness

# 运行示例
if __name__ == "__main__":
    nsga2 = NSGAII(pop_size=100, n_vars=10, m_obj=3, max_gen=200)
    pop, fitness = nsga2.evolve()
    # 可添加可视化代码查看帕累托前沿

内容的提问来源于stack exchange，提问作者Cuchisius