MOEA/D 的python代码实现

MOEA/D 的python代码实现
主要是参考了：
《MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition》
和《Multiobjective Optimization Problems With
Complicated Pareto Sets, MOEA/D and NSGA-II》两篇文章
主要是把第一篇文章用的SBX交叉算子换成了第二篇文章的差分进化和高斯变异，（但是也没有完全实现第二篇文章提到的做法）。
测试函数用的ZDT1测试函数 ，如果想测试别的还得自己按照我给的格式慢慢弄了。
代码只实现了二维的权重向量初始化方式，三维的还没有想好怎么等…
别的就是基本上参考了作者给出的.m代码实现啦，还给出了相应的解释。

import numpy as np
import matplotlib.pyplot as plt






class MOEAD():
    def __init__(self,objDim,parDim):

        self.objDim = objDim                # 问题的维数是多少
        self.parDim = parDim                # 输入变量的维数是多少
        self.popSize = 100                  # 种群的大小
        self.niche = 30                     # 邻居的个数
        self.iteration = 100                # 迭代次数
        self.F = 0.5                        # 变异概率
        self.CR = 0.5                       # 交叉概率


        self.idealPoint = np.full((self.objDim),np.inf)       # 初始化理想点
        self.weightVct = np.zeros((self.popSize,self.objDim))   # 权重向量
        self.nicheIndex = np.zeros((self.popSize,self.niche),dtype=int)   # 每个权重向量对应的邻居
        self.population = np.zeros((self.popSize,self.parDim))  # 初始种群
        self.fit = np.full((self.popSize,self.objDim),np.inf)   # 种群的适应度


    # 对每个个体计算适应度
    def cal_fitness(self,fitness_function,flag = 1,chrome = None):
        # 初始化的时候 直接计算
        if(flag == 1):
            for i in range(0,self.popSize):
                self.fit[i] = fitness_function(self.population[i].copy())
        # 需要单独对某个个体计算时
        else:
            return fitness_function(chrome)

    
    # 可能可以改进的地方
    def Init_Pop(self,lb,ub,fitness_function):
        delta = ub - lb
        self.population = np.random.rand(self.popSize, self.parDim) * delta + lb
        self.cal_fitness(fitness_function=fitness_function)


    # index:对应的要计算权重向量的下标
    def init_neighbors(self,index):
        distances = np.linalg.norm(self.weightVct - self.weightVct[index], axis=1)
        return np.argsort(distances)[:self.niche]


    # 自己手动写
    # 这个论文里每一个weight是(i / popsize,(popsize - i)/popsize)
    # 很显然三维就不行了..
    def init_weight_vector(self):
        if(self.objDim == 2):
            # 计算二维权重向量
            for i in range(0,self.popSize):
                self.weightVct[i][0] = i / self.popSize
                self.weightVct[i][1] = (self.popSize - i) / self.popSize
            # 找到对应的邻居
            for i in range(0,self.popSize):
                self.nicheIndex[i] = self.init_neighbors(index=i)
        else:
            raise "还没考虑二维以上问题写法"
        
    # 生成参考点 z1,z2,z3...
    # 在初始化之后找到每一维问题的最小值
    def init_reference_point(self):
        # 找到 arr1 中每一列的最小值
        min_values_arr1 = np.min(self.fit, axis=0)
        # 将 arr2 与 min_values_arr1 比较，保留更小的值
        self.idealPoint = np.minimum(min_values_arr1, self.idealPoint)


    # 计算切比雪夫距离
    # weight指的是：对应的权重向量
    # ref_point:对应的参考点
    def chebyshev_aggregation(self,x, weight, ref_point):
        return np.max(weight * np.abs(x - ref_point))
    

    # 更新index的子问题对应的邻居的解
    def update_neighbors_solution(self,index,offspring,offspring_eval):
        for j in self.nicheIndex[index]:
            if self.chebyshev_aggregation(offspring_eval, self.weightVct[j], self.idealPoint) < \
                self.chebyshev_aggregation(self.fit[j], self.weightVct[j], self.idealPoint):
                    print("找到了")
                    self.population[j] = offspring.copy()
                    self.fit[j] = offspring_eval.copy()
    # 加权求和的方法:论文中没有用到 后续可以用
    # def ws():pass

    # 交叉 - 变异

    # index:表示现在在优化的子问题序号
    # F_Scale:缩放因子（DE算法独有的参数）
    def crossover_DE(self,index,lb,ub,F_Scale=0.5):
        # 1. 首先找到3个 index的邻居:
        # 原文要求:不要包括index本身 且3个数不重复
        # 第0个数肯定是本身下标，因为欧氏距离是0
        select_indices =  np.random.choice(range(1,self.niche), size=3, replace=False)
        # 选中的染色体的下标
        select_chrome_index = self.nicheIndex[index][select_indices]
        # 2.生成新的染色体
        new_chrome = self.population[select_chrome_index[0]] + \
                    F_Scale * (self.population[select_chrome_index[1]] - self.population[select_chrome_index[2]])
        
        # 3.对染色体的每个维度进行变异
        random_sequence = np.random.rand(self.parDim)
        # 根据随机序列生成新数组
        random_sequence = (random_sequence <= self.F).astype(int)
        
        # 如果新数组全为0，则随机将一个位置置为1
        if np.all(random_sequence == 0):
            random_index = np.random.randint(0, self.parDim)
            random_sequence[random_index] = 1
        
        # 对数组进行限制
        clipped_array = np.clip(self.population[index] * random_sequence + \
                                new_chrome * (1 - random_sequence),lb, ub)
        return clipped_array

    # 高斯变异
    def mutation_Gauniss(self,chrome, lb, ub, prob):
        
        mutated_chromosome = chrome.copy()  # 复制染色体
        # 生成与染色体相同大小的随机矩阵，其中的元素服从均匀分布
        random_matrix = np.random.rand(*chrome.shape)
        # 根据变异概率和随机矩阵来确定哪些元素发生变异
        mutation_indices = random_matrix < prob
        # 对需要变异的元素应用高斯变异
        # 对需要变异的元素应用高斯变异，均值为当前染色体在该位置的值
        mutation_means = mutated_chromosome[mutation_indices]
        mutation_values = np.random.normal(mutation_means, (ub - lb) / 20)
        mutated_chromosome[mutation_indices] = mutation_values
        
        # 确保变异后的值在定义域范围内
        mutated_chromosome = np.clip(mutated_chromosome, lb, ub)
        
        return mutated_chromosome

    # 计算每一个DE得到的子问题的结果
    def evalute(self,fitness_function,index,lb,ub,F_Scale=0.5,mut_prob = 0):
        # 交叉
        cross_chrome = self.crossover_DE(index = index,lb=lb,ub=ub,F_Scale=F_Scale)
        #变异
        new_chrome = self.mutation_Gauniss(chrome=cross_chrome,lb=lb,ub=ub,prob=1/self.parDim)
        fit = self.cal_fitness(fitness_function=fitness_function,flag=0,chrome=new_chrome)
        return new_chrome,fit
    
    # 支配的计算
    def is_dominated(self, a, b):
        return np.all(a >= b) and np.any(a > b)

    # 找到当前的前沿
    def find_pareto_front(self):
        pareto_front = []
        for i in range(self.popSize):
            dominated = False
            for j in range(self.popSize):
                if self.is_dominated(self.fit[i], self.fit[j]):
                    dominated = True
                    break
            if not dominated:
                pareto_front.append(self.fit[i])
        return np.array(pareto_front)

    def plot_pareto_front(self, gen, pareto_front):
        plt.clf()
        plt.scatter(self.fit[:, 0], self.fit[:, 1], label='Population')
        plt.scatter(pareto_front[:, 0], pareto_front[:, 1], color='red', label='Pareto Front')
        plt.xlabel('Objective 1')
        plt.ylabel('Objective 2')
        plt.title(f'Pareto Front at Generation {gen}')
        # Plot ideal Pareto front
        f1_values = np.linspace(0, 1, 100)
        f2_values = 1 - np.sqrt(f1_values)
        plt.plot(f1_values, f2_values, 'g--', label='Ideal Pareto Front')
        plt.legend()
        plt.pause(0.1)  # Pause to update the plot

    def run(self,fitness_function,lb,ub):

        self.Init_Pop(fitness_function=fitness_function,lb=lb,ub=ub)    # 随机初始化种群并计算适应度
        self.init_weight_vector()       # 初始化权重向量和对应的邻居
        self.init_reference_point()     # 初始化参考点

        # 进行迭代
        plt.ion()  # Turn on interactive mode for dynamic plotting
        for i in range(self.iteration):
            for j in range(self.popSize):
                new_chrom,new_fit = self.evalute(fitness_function=fitness_function,index=j,
                                                    lb=lb,ub=ub)
                # 是否可以更新参考点
                self.idealPoint = np.minimum(self.idealPoint,new_fit)
                # 更新他的邻居
                self.update_neighbors_solution(index=j,offspring=new_chrom,offspring_eval=new_fit)
            pareto_front = self.find_pareto_front()
            self.plot_pareto_front(i, pareto_front)
        
        plt.ioff()  # Turn off interactive mode
        plt.show()

# objDIM:2
# parDim: >= 2
def ZDT1(chromosome):
    n = len(chromosome)
    f1 = chromosome[0]
    g = 1 + 9 * np.sum(chromosome[1:]) / (n - 1)
    h = 1 - np.sqrt(f1 / g)
    f2 = g * h
    return np.array([f1, f2])

if __name__ == "__main__":
    MOEA_ZDT = MOEAD(objDim = 2,parDim = 30)
    MOEA_ZDT.run(fitness_function=ZDT1,lb=0,ub=1)