MOPSO 多目标粒子群python实现

参考：

https://blog.csdn.net/m0_38097087/article/details/79818348

http://yarpiz.com/59/ypea121-mopso

Coello C A C , Pulido G T , Lechuga M S . Handling multiple objectives with particle swarm optimization[J]. IEEE Transactions on Evolutionary Computation, 2004, 8(3):256-279.。

是按照上面的第三个的那篇论文实习实现的

论文里面的速度更新公式没有用到惯性c，好像是因为他有个变异算子(还没实现)。

算法流程

初始化群体粒子群的位置和速度，计算适应值
评价粒子的适应度和Pareto支配关系
将非劣解保存到Archive中去
计算Archive集中的拥挤度，为粒子选择gbest
更新粒子的速度、位置、适应度、pbest
更新Archive
满足结束条件，则结束；否则，转到第4步继续循环。

import numpy as np
class Partical:
    #每次都是对一个粒子进行操作
    def __init__(self,x_min,x_max,max_v,min_v,fitness):
        self.dim=len(x_min) #获得变量数
        self.max_v=max_v
        self.min_v=min_v
        self.x_min=x_min
        self.x_max=x_max
        
        self.dominated=False #表示是否被支配
        
        self.pos=np.zeros(self.dim)
        self.pbest=np.zeros(self.dim)
        self.initPos(x_min,x_max)#初始化粒子坐标&pbest
        
        self._v=np.zeros(self.dim)

        self.fitness=fitness
        self.cur_fitness=fitness(self.pos)#当前的适应度
        self.bestFitness=fitness(self.pos)#初始化pbest
        
    def _updateFit(self):
        if isDominates(np.array(self.cur_fitness),np.array(self.bestFitness)):
            self.bestFitness=self.cur_fitness
            self.pbest=self.pos
        elif isDominates(np.array(self.bestFitness),np.array(self.cur_fitness)):
            pass
        else:
            #互不支配随机选择一个
            if np.random.random()<0.5:
                self.bestFitness=self.cur_fitness
                self.pbest=self.pos
    
    def _updatePos(self):
        self.pos=self.pos+self._v
        for i in range(self.dim):
            self.pos[i]=min(self.pos[i],self.x_max[i])
            self.pos[i]=max(self.pos[i],self.x_min[i])
            
    def _updateV(self,w,c1,c2,gbest):
        '''这里进行的都是数组的运算'''
        self._v=w*self._v+c1*np.random.random()*(self.pbest-self.pos)+c2*np.random.random()*(gbest-self.pos)
        for i in range(self.dim):
            self._v[i]=min(self._v[i],self.max_v[i])
            self._v[i]=max(self._v[i],self.min_v[i])
            
    def initPos(self,x_min,x_max):
        for i in range(self.dim):
            self.pos[i]=np.random.random()*(x_max[i]-x_min[i])+ x_min[i]#unifom是左闭右开取不到最大值
            self.pbest[i]=self.pos[i]
            
    def getPbest(self):
        return self.pbest
    
    def getBestFit(self):
        return self.bestFitness
    
    def setIsDominated(self,r):
        self.dominated=r
        
    def getIsDominated(self):
        return self.dominated
    
    def update(self,w,c1,c2,gbest):
        self._updateV(w,c1,c2,gbest)
        self._updatePos()
        self.cur_fitness=self.fitness(self.pos)
        
        #先不加变异的部分，加的话应该加在这里
        self._updateFit()

判断是否被支配

链接一的实现感觉不大对

def isDominates(x,y):
    #x是否支配y
    return (x<=y).all() and (x<y).any()
def determinDomination(pop):#确定是否支配,会存在重复比较(已经确定了被支配的劣解实际可以不参与之后的比较)
    for i in range(len(pop)):
        pop[i].setIsDominated(False)
    for i in range(0,len(pop)-1):
        for j in range(i+1,len(pop)):
            if isDominates(np.array(pop[i].cur_fitness),np.array(pop[j].cur_fitness)):
                pop[j].setIsDominated(True)#j被i支配
            if isDominates(np.array(pop[j].cur_fitness),np.array(pop[i].cur_fitness)):
                pop[i].setIsDominated(True)

处理存档的方法

如果外部存档为空，那么非支配解直接进入
如果新的解被存档中的解支配，就丢弃这个解
如果新的解不被存档中的任意一个解支配就进入存档中
如果存档中的解被这个新的解支配就删去存档中被支配的解
如果存档满了调用自适应网格程序

全局最优的粒子的选择

论文中写得方法是用10除以每个网格中的粒子数作为每个网格的适应度，然后使用轮盘赌的方法选出一个网格，如果网格中的粒子数目不止一个就随机选择一个。

但在参考链接二的实现中是用P=exp(-beta*N); P=P/sum(P）;作为每个网格的适应度（概率）然后再用轮盘赌。其中beta表示选择压力，可能beta越小选择压力越大？？因为他们的概率差别越小。原因应该就是为了轮盘赌的通用性，存档中删去粒子也是要用轮盘赌的。

参考链接一用的是10，但除了立方没有很明白，[为了增加选择全局最佳的压力，我试了一下效果会好特别多]。需要再去看下论文。然后他的轮盘赌是所有粒子进行的。我实现的是网格进行轮盘赌。

外部存档和网格的实现

#感觉这个就应该用单例模式，不会写。。
class Archive:
    def __init__(self,partical,size,mesh_div=10):
        self.mesh_div=mesh_div #网格的数量
        self.particals=partical
        self.size=size#存档的大小
        determinDomination(self.particals)#判断初始化种群之间的支配关系
        tempArchive=[x for x in self.particals if x.dominated!=True]#将初始化种群中的非支配解放入外部存档中
        self.archive=copy.deepcopy(tempArchive)
        self.mesh_id=np.zeros(len(self.archive))#初始化存档中粒子的网格编号
        
        self.inflation=0.1
        self.fitness=[par.cur_fitness for par in self.archive] #每个粒子的适应度值不止一个,list里面只存数值
        self.f_max=np.max(np.array(self.fitness),axis=0)#求网格的取值范围
        self.f_min=np.min(np.array(self.fitness),axis=0)          
        
    def createGrid(self):
        '''存档中要是只有一个粒子就返回吧，也可以给他创建个网格'''
        if len(self.archive)==1:
            return
        self.f_max=np.max(np.array(self.fitness),axis=0)#求网格的取值范围
        self.f_min=np.min(np.array(self.fitness),axis=0)  
        dc=self.f_max-self.f_min
        f_max=self.f_max+self.inflation*dc
        f_min=self.f_min-self.inflation*dc #比范围稍微大一点,都是数组运算
        self.mesh_id=np.zeros(len(self.archive))        
        for i in range(len(self.archive)):
                   self.mesh_id[i]=self._cal_mesh_id(self.fitness[i])
        
    def _cal_mesh_id(self,fit):
        id_=0
        for i in range(len(fit)):#分别对每个目标函数进行计算
            '''首先，将每个维度按照等分因子进行等分离散化，获取粒子在各维度上的编号。按照10进制将每一个维度编号等比相加
            （如过用户自定义了mesh_div_num的值，则按照自定义），计算出值,第一次就只有一个非支配'''
            id_dim=int((fit[i]-self.f_min[i])*self.mesh_div/(self.f_max[i]-self.f_min[i]))#为什么乘了粒子数，我没理解
            id_ = id_ + id_dim*(self.mesh_div**i)
        return id_
    
    def RouletteWheelSelection(self,prob):#轮盘赌
        p=np.random.random()
       # c=np.cumsum(p) #反正我都要遍历，遍历的时候累加好了。。。
        cunsum=0
        for i in range(len(prob)):
            cunsum+=prob[i]
            if p<=cunsum:
                return i
        
    def selectLeader(self):
        mesh_id=self.mesh_id
        unique_elements,counts_elements = np.unique(mesh_id,return_counts=True)#所占的网格的id
        print(unique_elements,counts_elements)
        counts_elements=10/(counts_elements**3) #按照论文的说法
        prob=counts_elements/np.sum(counts_elements)
        idx=self.RouletteWheelSelection(prob)#通过轮盘赌选择一个网格
        smi=np.where(self.mesh_id==unique_elements[idx])#获得该网格的粒子索引值
        smi=list(smi)[0]
        select=np.random.randint(0,len(smi))
        return self.archive[smi[select]]
    
    def deletePartical(self):
        mesh_id=self.mesh_id
        unique_elements,counts_elements = np.unique(mesh_id,return_counts=True)#所占的网格的id
        p=np.exp(counts_elements)
        p=p/np.sum(p)
        idx=self.RouletteWheelSelection(p)       
        smi=np.where(self.mesh_id==unique_elements[idx])#获得该网格的粒子索引值，返回值是tuple
        smi=list(smi)[0]
        select=np.random.randint(0,len(smi))       
        del self.archive[smi[select]]
        self.mesh_id=np.delete(self.mesh_id,smi[select])
        del self.fitness[smi[select]]
    
    def update(self,particals):
        self.particals=particals
        determinDomination(self.particals)        
        #将本次迭代的非支配解加到存档中
        curr_temp=[x for x in self.particals if x.dominated!=True]
        curr=copy.deepcopy(curr_temp)
        self.archive+=curr
        
        determinDomination(self.archive)
        curr_temp=[x for x in self.archive if x.dominated!=True]
        self.archive=copy.deepcopy(curr_temp)
        self.fitness=[par.cur_fitness for par in self.archive] 
        self.createGrid() #如果粒子的范围超出了原来的范围更新网格（我索引计划全重新算，这里就直接create吧）
        while len(self.archive)>self.size:
            #按照密度利用轮盘赌删掉一些粒子
            self.deletePartical()

主流程

class Mopso:
    def __init__(self,pop,fitnessFunction,w,c1,c2,max_,min_,thresh,mesh_div=10):
        self.w,self.c1,self.c2=w,c1,c2
        self.mesh_div=mesh_div
        self.pop=pop#种群大小
        self.thresh=thresh
        self.max_=max_
        self.min_=min_
        self.max_v=(max_-min_)*0.05 
        self.min_v=(max_-min_)*0.05*(-1)
        self.fitnessFunction=fitnessFunction#适应度函数这里是个数组
        
        #初始化种群,坐标，速度,pbest...
        self.particals=[Partical(self.min_,self.max_,self.max_v,self.min_v,self.fitnessFunction) for i in range(self.pop)]
        #初始化外部存档，
        self.archive=[]
       
    def evaluation_fitness(self):
        fitness_curr=[]  #
        for i in range((self.in_).shape[0]):
            fitness_curr.append(fitness_(self.in_[i]))
        self.fitness_=np.array(fitness_curr)
        
    def update_(self):
        #更新粒子坐标，速度，适应值，个体最优，全局最优，外部存档
        for part in self.particals:
            gbestPart=self.archive.selectLeader()#选择全局最优
            gbest=gbestPart.pos
            part.update(self.w,self.c1,self.c2,gbest)        
        #对存档进行更新
        self.archive.update(self.particals)            
            
    def done(self,cycle_):
        print("in")
        self.archive=Archive(self.particals,self.thresh)
        self.archive.createGrid()
        for i in range(cycle_):
            print("="*10,i,"="*10)
            self.update_()
        return self.archive

我的参数设置

w=0.4#惯性，论文里用的是0.4
    c1=1#局部速度
    c2=2#全局速度
    particals=200
    gengeration=200#迭代次数
    mesh_div=10#网格等分数量
    thresh=100#外部存档数
    var_num=30#决策变量数目
    min_=np.zeros(30)
    max_=np.ones(30)
    mopso_ = Mopso(particals,fitness_,w,c1,c2,max_,min_,thresh,mesh_div) #粒子群实例化
    pareto_front = mopso_.done(gengeration) #经过cycle_轮迭代后，外部存档里面的结果

结果

试了一下ZDT1函数

c1=1;c2=2(同样的参数又跑了几次，拟合效果一点也不稳定。。下图这次是唯一比较好的一次)

MOPSO 多目标粒子群python实现

红色的是正确答案，蓝色是拟合的帕累托前沿

应该还是有问题，为什么这么短,参数的影响非常大

当我把c2改成0.5时，线是边长了，但拟合效果emmmm什么垃圾

MOPSO 多目标粒子群python实现

c1=0.5;c2=1;粒子数80，迭代300次

MOPSO 多目标粒子群python实现

论文结尾说他们发现算法对惯性权重十分敏感所以有了变异算子。。emmm看来有必要实现一下

MOPSO 多目标粒子群python实现

算法流程

外部存档和网格的实现

主流程

结果

推荐参数

继续阅读

一种以超体积测量（S度量）为选择准则的EMO算法-SMS-EMOA

多目标量子粒子群优化的经济排放调度问题（Matlab代码实现）

学习笔记--基于范围选择的进化多目标优化（PESA-II）

多目标粒子群_深度学习经典算法 | 粒子群算法详解

【学习笔记】多目标优化问题分解成若干简单多目标子问题--MOEA/D-M2M算法简介主要思想

多目标优化算法之SPEA2

多目标优化问题中常见分解方法的理解

单目标&多目标优化算法的测试函数与解

HyperVolume多目标评价指标

多目标进化算法的性能评价指标总结（一）

利用NSGA-II算法求解多元多目标函数优化问题

多目标优化之帕累托最优

基于C#的多目标进化算法平台MOEAPlat实现

【学习笔记】基于参考点的非支配排序--NSGA-III简介问题详述关于多目标问题的两个想法NSGA-III算法

【智能优化算法-差分进化算法】基于多种交叉策略和变异策略的差分进化算法求解单目标优化问题含Matlab源码