iT邦幫忙

0

[PSO文獻2]Adaptive Particle Swarm Optimization

2020.06.24 更新錯誤的Schwefel,現在已經收斂到跟原作者一樣好了


0. 前言:

本篇主要是紀錄Adaptive Particle Swarm Optimization的閱讀心得。

1. 動機:

最近腦子動到想用PSO來訓練神經網路,所以蒐集了一些相關資料,預計將分享5篇文章

2. 正文

原文的粒子定義如下,在一個D維度解空間中,粒子i的速度和位置定義如下:
https://ithelp.ithome.com.tw/upload/images/20200604/20124766n3GSsfTwGR.png
https://ithelp.ithome.com.tw/upload/images/20200604/20124766kaf55wZFWQ.png
粒子i在任意維度d的速度和位置更新公式:
https://ithelp.ithome.com.tw/upload/images/20200604/20124766qLmTd5c5LK.png
https://ithelp.ithome.com.tw/upload/images/20200604/20124766MjptRxiopA.png
原文提到PSO有兩個版本,分別是global-version PSO(GPSO)和local-version PSO(LPSO)。在Advances of Computational Intelligence in Industrial Systems的第9頁和Manufacturing Optimization through Intelligent Techniques的第41頁有提到GPSO和LPSO的差別以及優缺點:

  • GPSO的C1和C2不為0,速度快但容易陷入區域最佳解
  • LPSO的C2為0(不考慮其他粒子的最佳解),速度慢但不容易陷入區域最佳解

3. APSO的計算流程與算式

相較於傳統PSO,APSO在計算過程中多了演化狀態評估(Evolutionary State Estimation, ESE)和菁英學習策略(Elitist Learning Strategy, ELS)兩個步驟。

  • ESE可以根據粒子在每次迭代的分布情形f來調整加速度常數C、慣性權重W,進而加速粒子收斂
  • 對於gBest提供專屬的學習策略,進而提升所有粒子的求解效率

3.1 整體流程

不囉嗦,直接貼圖。如果刪去流程圖的第三步驟,那就是傳統PSO
https://ithelp.ithome.com.tw/upload/images/20200620/20124766A3124XJeJg.png

作者定義的加速度限制為解空間的20%,也就是
https://ithelp.ithome.com.tw/upload/images/20200620/20124766MeyUxMnzfh.png

由於作者沒有特別說明,因此初始加速度我令為
https://ithelp.ithome.com.tw/upload/images/20200620/20124766vjc286Z4jG.png

3.2 演化狀態評估ESE

根據每次迭代粒子的分布情形來調整w、c1、c2
https://ithelp.ithome.com.tw/upload/images/20200620/20124766XoI9HW4OaY.png

3.2.1 計算粒子彼此間的距離

https://ithelp.ithome.com.tw/upload/images/20200620/20124766x4RPu6SCNc.png
很明顯這是在求L2範數(歐式距離),然後為了編程方便,我把式子簡化如下
https://ithelp.ithome.com.tw/upload/images/20200620/20124766XMuTAaXf06.png

3.2.2 計算演化因子(evolutionary factor, f)

https://ithelp.ithome.com.tw/upload/images/20200620/20124766dg92HA4Zyd.png
https://ithelp.ithome.com.tw/upload/images/20200620/20124766jZucD5EKJB.png
https://ithelp.ithome.com.tw/upload/images/20200620/20124766NHA5ubMMn6.png

  • dg是當前適應值最好的粒子的di
  • 如果分母為0則直接定義為收斂狀態S3

3.2.3 模糊推論

  • 作者對粒子群分布狀況作了4個類別:探勘(exploration, S1)、開發(exploitation, S2)、收斂(Convergence, S3)、跳出(Jumping-out, S4)
  • 對f作模糊推論,得知當前粒子分布是屬於哪一類別
    https://ithelp.ithome.com.tw/upload/images/20200620/20124766JAfP0sa4p5.png
    https://ithelp.ithome.com.tw/upload/images/20200620/20124766mIzfwaoVlY.png

原文沒有完整說明規則庫內容,我只能透過內文的範例進行推論,結果如下
https://ithelp.ithome.com.tw/upload/images/20200620/20124766MGyuiXT5eU.png

  • 必須遵守S1=>S2=>S3=>S4=>S1=>...這個順序
  • 為了求解穩定,能不改變狀態就不改變
  • 即使要改變狀態,每次也只向右1步或往左1步
  • 初始狀態我是令S1

3.2.4 慣性權重

https://ithelp.ithome.com.tw/upload/images/20200620/20124766q08HyC5cd4.png
原文的初始慣性權重為0.9

3.2.5 加速度常數

加速度常數的範圍需落在,初始值為2.0
https://ithelp.ithome.com.tw/upload/images/20200620/20124766y6M09x5gml.png
https://ithelp.ithome.com.tw/upload/images/20200620/201247662v9AYeE99A.png
如果c1+c2大於指定範圍,則用下式調整
https://ithelp.ithome.com.tw/upload/images/20200620/20124766SOoyyd9aXk.png

根據不同的狀態,以加速度率(acceleration rate)調整加速度常數

  • 加速度率不會隨著迭代而改變
    https://ithelp.ithome.com.tw/upload/images/20200620/201247661ofgj7UvlX.png
    https://ithelp.ithome.com.tw/upload/images/20200620/201247660x0tCe4Kah.png
    https://ithelp.ithome.com.tw/upload/images/20200620/20124766gelknoXVSO.png

原文有明確定義所謂的輕微是指0.5倍的加速度率,因此可以得到下面結論
https://ithelp.ithome.com.tw/upload/images/20200620/20124766xEWjehHqE7.png

3.3 菁英學習策略ELS

想像一下,某個粒子要更新自己的加速度,但是他發現自己就是pBest,同時他又是gBest,所以加速度更新公式的三個子項中就有2項是0。由於這個粒子在該次代沒有任何可以學習的對象,因此他本次移動是毫無根據的,這會間接導致粒子子求解速度被拖慢,因此制定一個擾動策略(類似基因算法的突變)來嘗試得到更好的gBest。當然,不是每次擾動都會得到更好的gBest,但我們可以取代最差的粒子,使得整體表現變得更均勻
https://ithelp.ithome.com.tw/upload/images/20200620/20124766TdXKzO2ihy.png

  • 只有在狀態為S3時才執行ELS!!!
  • 1)令P=gBest,並且隨機取得某一個維度位置d
  • 2)改變位置d的數值。其中平均值u=0
    https://ithelp.ithome.com.tw/upload/images/20200620/20124766k4DQ6mnyD6.png
    https://ithelp.ithome.com.tw/upload/images/20200620/20124766SMVT4GrIby.png
    https://ithelp.ithome.com.tw/upload/images/20200620/20124766Lh8V2f8XmV.png
  • 3)確保位置d仍在解空間內
  • 4)計算P的適應值v
  • 5-1)如果適應值v比先前的gBest還要好,則P取代gBest
  • 5-2)如果適應值v沒有先前的gBest好,則P取代當前適應值最差的粒子

我自己對5-2作了修改,如下

  • 新5-2)如果適應值v沒有先前的gBest好,但確定勝過當前適應值最差的粒子才取代,同時如果勝過最差的pBest則也進行替換
  • 原因是在沒有確認誰強誰弱就直接取代,怎麼想都覺得有點奇怪

4. 實現

4.1 APSO代碼

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)

class APSO():
    def __init__(self, fit_func, num_dim, num_particle=20, max_iter=500,
                 x_max=1, x_min=-1, w_max=0.9, w_min=0.4, c1=2.0, c2=2.0, k=0.2):
        self.fit_func = fit_func        
        self.num_dim = num_dim
        self.num_particle = num_particle
        self.max_iter = max_iter
        
        self.x_max = x_max
        self.x_min = x_min
        self.w = w_max        
        self.w_max = w_max
        self.w_min = w_min
        self.c1 = c1
        self.c2 = c2
        self.k = k
        self._iter = 1
        self.bestX_idx = -1
        self.worstX_idx = -1
        # case1. 僅初始化
        self.delta = np.random.uniform(high=0.1, low=0.05, size=1)
        self.state = 1
        self.rulebase = [[1, 1, 4, 1], 
                         [2, 2, 2, 1], 
                         [2, 3, 3, 3], 
                         [4, 3, 4, 4]]
        self.gBest_curve = np.zeros(self.max_iter)
 
        self.X = np.random.uniform(size=[self.num_particle, self.num_dim])*(self.x_max-self.x_min) + self.x_min
        self.V = np.random.uniform(size=[self.num_particle, self.num_dim])*(self.x_max-self.x_min) + self.x_min
        self.v_max = self.k*(self.x_max-self.x_min)
        
        self.pBest_X = self.X.copy()
        self.pBest_score = self.fit_func(self.X).copy()
        self.X_score = self.pBest_score.copy()
        self.gBest_X = self.pBest_X[self.pBest_score.argmin()].copy()
        self.gBest_score = self.pBest_score.min().copy()
        self.gBest_curve[0] = self.gBest_score.copy()
        self.bestX_idx = self.pBest_score.argmin()
        self.worstX_idx = self.pBest_score.argmax()

    def ESE(self):
        # step 1.
        d_set = np.zeros(self.num_particle)
                      
        for i in range(self.num_particle):
            d_set[i] = np.sum(np.linalg.norm((self.X[i, :].reshape(1, self.num_dim)-self.X), ord=2, axis=1))\
                /(self.num_particle-1)
        
        # step 2.        
        dmax = np.max(d_set)
        dmin = np.min(d_set)
        dg = d_set[self.bestX_idx].copy()        
        # 我自己額外加的,防止崩潰
        if dmax==dmin==dg==0:
            f = 0
        else:
            f = (dg-dmin)/(dmax-dmin)

        if not(0<=f<=1):
            print('f='+str(np.round(f, 3))+' 超出範圍[0, 1]!!!')            
            if f>1:
                f = 1
            elif f<0:
                f = 0 
            
        # step 3.
        singleton = np.zeros(4)
        if 0<=f<=0.4:
            singleton[0] = 0
        elif 0.4<f<=0.6:
            singleton[0] = 5*f - 2
        elif 0.6<f<=0.7:
            singleton[0] = 1
        elif 0.7<f<=0.8:
            singleton[0] = -10*f + 8
        elif 0.8<f<=1:
            singleton[0] = 0

        if 0<=f<=0.2:
            singleton[1] = 0
        elif 0.2<f<=0.3:
            singleton[1] = 10*f - 2
        elif 0.3<f<=0.4:
            singleton[1] = 1
        elif 0.4<f<=0.6:
            singleton[1] = -5*f + 3
        elif 0.6<f<=1:
            singleton[1] = 0

        if 0<=f<=0.1:
            singleton[2] = 1
        elif 0.1<f<=0.3:
            singleton[2] = -5*f + 1.5
        elif 0.3<f<=1:
            singleton[2] = 0
            
        if 0<=f<=0.7:
            singleton[3] = 0
        elif 0.7<f<=0.9:
            singleton[3] = 5*f - 3.5
        elif 0.9<f<=1:
            singleton[3] = 1
        
        if np.max(singleton)==np.min(singleton)==0:
            print('因為f超出範圍[0, 1],所以模糊推論異常!!!')
        
        t = np.argmax(singleton)
        t_1 = self.state - 1
        self.state = self.rulebase[int(t)][int(t_1)]
    
        # (10)
        self.w = 1/(1+1.5*np.exp(-2.6*f))
        if not(self.w_min<=self.w<=self.w_max):
            print('w='+str(np.round(self.w, 3))+' 超出範圍[0.4, 0.9]!!!')
            if self.w<self.w_min:
                self.w = self.w_min
            elif self.w>self.w_max:
                self.w = self.w_max
                   
        # (11)
        if self.state==1:
            self.c1 = self.c1 + self.delta
            self.c2 = self.c2 - self.delta
        if self.state==2:
            self.c1 = self.c1 + 0.5*self.delta
            self.c2 = self.c2 - 0.5*self.delta
        if self.state==3:
            self.c1 = self.c1 + 0.5*self.delta
            self.c2 = self.c2 + 0.5*self.delta
            self.ELS()
        if self.state==4:
            self.c1 = self.c1 - self.delta
            self.c2 = self.c2 + self.delta
        if not(1.5<=self.c1<=2.5):
            # print('c1='+str(np.round(self.c1, 3))+' 超出範圍[1.5, 2.5]!!!') 
            if self.c1<1.5:
                self.c1 = 1.5
            elif self.c1>2.5:
                self.c1 = 2.5
        if not(1.5<=self.c2<=2.5):
            # print('c2='+str(np.round(self.c2, 3))+' 超出範圍[1.5, 2.5]!!!')  
            if self.c2<1.5:
                self.c2 = 1.5
            elif self.c2>2.5:
                self.c2 = 2.5
            
        # (12)
        if not(3<=self.c1+self.c2<=4):
            # print('c1='+str(np.round(self.c1, 3))+' + c2='+str(np.round(self.c2, 3))+' 超出範圍[3, 4]!!!')
            self.c1, self.c2 = 4.0*self.c1/(self.c1+self.c2), 4.0*self.c2/(self.c1+self.c2)            
    def ELS(self):
        rho = 1 - ((1-0.1)*self._iter/self.max_iter)
        d = np.random.randint(low=0, high=self.num_dim, size=1)
        P = self.gBest_X.copy()
        P[d] = P[d] + (self.x_max[d]-self.x_min[d])*np.random.normal(loc=0.0, scale=rho**2, size=1)
        
        if P[d]>self.x_max[d]: P[d] = self.x_max[d]
        if P[d]<self.x_min[d]: P[d] = self.x_min[d]
        
        score = self.fit_func(P.reshape(1, -1))
        
        if score<self.gBest_score:
            self.gBest_X = P.copy()
            self.gBest_score = score.copy()
        else: 
            # # case1. 原文版本
            # self.X[self.worstX_idx, :] = P.copy()
            
            # case2. 我的版本
            if score<self.X_score[self.worstX_idx]:
                self.X[self.worstX_idx, :] = P.copy()
            if score<np.max(self.pBest_score):
                idx = np.argmax(self.pBest_score)
                self.pBest_score[idx] = score.copy()
                self.pBest_X[idx, :] = P.copy()  

    def opt(self):
        while(self._iter<self.max_iter):   
            # case2. 每代更新
            # self.delta = np.random.uniform(high=0.1, low=0.05, size=1)
            self.bestX_idx = self.X_score.argmin()
            self.worstX_idx = self.X_score.argmax()
            
            self.ESE()
            R1 = np.random.uniform(size=(self.num_particle, self.num_dim))
            R2 = np.random.uniform(size=(self.num_particle, self.num_dim))
            for i in range(self.num_particle):                              
                self.V[i, :] = self.w * self.V[i, :] \
                        + self.c1 * (self.pBest_X[i, :] - self.X[i, :]) * R1[i, :] \
                        + self.c2 * (self.gBest_X - self.X[i, :]) * R2[i, :]               
                self.V[i, self.v_max < self.V[i, :]] = self.v_max[self.v_max < self.V[i, :]]
                self.V[i, -self.v_max > self.V[i, :]] = -self.v_max[-self.v_max > self.V[i, :]]
                
                self.X[i, :] = self.X[i, :] + self.V[i, :]
                self.X[i, self.x_max < self.X[i, :]] = self.x_max[self.x_max < self.X[i, :]]
                self.X[i, self.x_min > self.X[i, :]] = self.x_min[self.x_min > self.X[i, :]]
                
                score = self.fit_func(self.X[i, :])
                self.X_score[i] = score.copy()
                if score<self.pBest_score[i]:
                    self.pBest_X[i, :] = self.X[i, :].copy()
                    self.pBest_score[i] = score.copy()
                    if score<self.gBest_score:                        
                        self.gBest_X = self.X[i, :].copy()
                        self.gBest_score = score.copy()
                                       
            self.gBest_curve[i] = self.gBest_score.copy()         
            self._iter = self._iter + 1            


    def plot_curve(self):
        plt.figure()
        plt.title('loss curve ['+str(round(self.gBest_curve[-1], 3))+']')
        plt.plot(self.gBest_curve, label='loss')
        plt.grid()
        plt.legend()
        plt.show()      

    

4.2 測試代碼

from APSO import APSO
import numpy as np
import time
import pandas as pd
np.random.seed(42)

def Sphere(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    return np.sum(x**2, axis=1)

def Schwefel_P222(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    return np.sum(np.abs(x), axis=1) + np.prod(np.abs(x), axis=1)

def Quadric(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    
    outer = 0
    for i in range(x.shape[1]):
        inner = np.sum(x[:, :i+1], axis=1)**2
        outer = outer + inner
    
    return outer

def Rosenbrock(x):
    if x.ndim==1:
        x = x.reshape(1, -1) 
    
    left = x[:, :-1].copy()
    right = x[:, 1:].copy()
    
    return np.sum(100*(right - left**2)**2 + (left-1)**2, axis=1)

def Step(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    return np.sum(np.round((x+0.5), 0)**2, axis=1)

def Quadric_Noise(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
        
    matrix = np.arange(x.shape[1])+1
     
    return np.sum((x**4)*matrix, axis=1)+np.random.rand(x.shape[0])

def Schwefel(x):
    if x.ndim==1:
        x = x.reshape(1, -1)        
     
    return -1*np.sum(x*np.sin(np.abs(x)**.5), axis=1)

def Rastrigin(x):
    if x.ndim==1:
        x = x.reshape(1, -1) 
    
    return np.sum(x**2 - 10*np.cos(2*np.pi*x) + 10, axis=1)

def Noncontinuous_Rastrigin(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    
    outlier = np.abs(x)>=0.5
    x[outlier] = np.round(2*x[outlier])/2
    
    return np.sum(x**2 - 10*np.cos(2*np.pi*x) + 10, axis=1)

def Ackley(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    
    left = 20*np.exp(-0.2*(np.sum(x**2, axis=1)/x.shape[1])**.5)
    right = np.exp(np.sum(np.cos(2*np.pi*x), axis=1)/x.shape[1])
    
    return -left - right + 20 + np.e

def Griewank(x):
    if x.ndim==1:
        x = x.reshape(1, -1)
    left = np.sum(x**2, axis=1)/4000
    right = np.prod( np.cos(x/((np.arange(x.shape[1])+1)**.5)), axis=1)
    return left - right + 1

d = 30
g = 1000
p = 20
times = 30
table = np.zeros((2, 11))
for i in range(times):
    x_max = 100*np.ones(d)
    x_min = -100*np.ones(d)
    optimizer = APSO(fit_func=Sphere, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 0] += optimizer.gBest_score
    table[1, 0] += end - start
    
    x_max = 10*np.ones(d)
    x_min = -10*np.ones(d)
    optimizer = APSO(fit_func=Schwefel_P222, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 1] += optimizer.gBest_score
    table[1, 1] += end - start
    
    x_max = 100*np.ones(d)
    x_min = -100*np.ones(d)
    optimizer = APSO(fit_func=Quadric, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 2] += optimizer.gBest_score
    table[1, 2] += end - start
    
    x_max = 10*np.ones(d)
    x_min = -10*np.ones(d)
    optimizer = APSO(fit_func=Rosenbrock, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 3] += optimizer.gBest_score
    table[1, 3] += end - start
    
    x_max = 100*np.ones(d)
    x_min = -100*np.ones(d)
    optimizer = APSO(fit_func=Step, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 4] += optimizer.gBest_score
    table[1, 4] += end - start
    
    x_max = 1.28*np.ones(d)
    x_min = -1.28*np.ones(d)
    optimizer = APSO(fit_func=Quadric_Noise, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 5] += optimizer.gBest_score
    table[1, 5] += end - start
    
    x_max = 500*np.ones(d)
    x_min = -500*np.ones(d)
    optimizer = APSO(fit_func=Schwefel, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 6] += optimizer.gBest_score
    table[1, 6] += end - start
    
    x_max = 5.12*np.ones(d)
    x_min = -5.12*np.ones(d)
    optimizer = APSO(fit_func=Rastrigin, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 7] += optimizer.gBest_score
    table[1, 7] += end - start
    
    x_max = 5.12*np.ones(d)
    x_min = -5.12*np.ones(d)
    optimizer = APSO(fit_func=Noncontinuous_Rastrigin, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 8] += optimizer.gBest_score
    table[1, 8] += end - start
    
    x_max = 32*np.ones(d)
    x_min = -32*np.ones(d)
    optimizer = APSO(fit_func=Ackley, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 9] += optimizer.gBest_score
    table[1, 9] += end - start
    
    x_max = 600*np.ones(d)
    x_min = -600*np.ones(d)
    optimizer = APSO(fit_func=Griewank, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)
    start = time.time()
    optimizer.opt()
    end = time.time()
    table[0, 10] += optimizer.gBest_score
    table[1, 10] += end - start
    
table = table / times
table = pd.DataFrame(table)
table.columns=['Sphere', 'Schwefel_P222', 'Quadric', 'Rosenbrock', 'Step', 'Quadric_Noise', 'Schwefel', 
                'Rastrigin', 'Noncontinuous_Rastrigin', 'Ackley', 'Griewank']
table.index = ['mean score', 'mean time']

4.3 討論

這邊主要是探討原文的ELS和我的ELS間對計算結果的差異,以及delta是要在每代都更新一次或者在初始化階段完成定義後就不再更動,因此會有四個組合的比較。

因為不想等太久,所以

  • 共同環境條件為:解空間維度d=30; 最大迭代g=1000; 粒子數p=20; 執行次數times=30; 每跑完一條舊更新gBest
  • 這個條件一定跑不到理想
    https://ithelp.ithome.com.tw/upload/images/20200724/201247666hHzRHOcmZ.png
    https://ithelp.ithome.com.tw/upload/images/20200724/20124766AuiEsAMbUK.png
    可以看出我自己定義的ELS效果比原文好,而delta對於最終結果感覺沒有什麼關係

然後是採我的ELS+delta部每代更新的方式,將最大迭代g改成10000後其他條件不變的結果
https://ithelp.ithome.com.tw/upload/images/20200724/20124766dPOjTeJlWS.png

5. 結論

  • 魔改後的效果比原作者好,雖然我知道這樣改有點奇怪就是了哈哈
  • 改成10000代後跑得有夠慢,建議用cupy或者GPU版的pytorch改寫後放在GPU上面跑會比較好

圖片
  直播研討會
圖片
{{ item.channelVendor }} {{ item.webinarstarted }} |
{{ formatDate(item.duration) }}
直播中

尚未有邦友留言

立即登入留言