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源碼