1. 程式人生 > >【尋優演算法】粒子群演算法(PSO) 引數尋優的python實現

【尋優演算法】粒子群演算法(PSO) 引數尋優的python實現

【尋優演算法】粒子群演算法(PSO) 引數尋優的python實現


粒子群優化演算法(Particle swarm optimization,PSO)是模擬鳥群捕食行為的優化演算法。不同於遺傳演算法(Genetic Alogrithm,GA),粒子群演算法是有記憶的,之前迭代過程中的最優位置和最優方向都會保留下來並作用於粒子群的更新【參考資料1】。

一、演算法原理

1、粒子群演算法的名詞解釋

粒子群長度:粒子群長度等於每一個引數取值範圍的大小。
粒子群維度:粒子群維度等於待尋優引數的數量。
粒子群位置:粒子群位置包含引數取值的具體數值。
粒子群方向:粒子群方向表示引數取值的變化方向。
適應度函式:表徵粒子對應的模型評價指標。
pbest:(區域性最優)pbest的長度等於粒子群長度,表示每一個引數取值的變化過程中,到目前為止最優適應度函式值對應的取值。
gbest:(全域性最優)gbest的長度為1,表示到目前為止所有適應度函式值中最優的那個對應的引數取值。
慣性因子 w

w :慣性因子表示粒子保持的運動慣性。
區域性學習因子 c 1 {c_1} :表示每個粒子向該粒子目前為止最優位置運動加速項的權重。
全域性學習因子 c
2 {c_2}
:表示每個粒子向目前為止全域性最優位置運動加速項的權重。

2、粒子更新

粒子方向更新方程為:

v i = w × v i + c 1 × r a n d ( ) × ( p b e s t i x i ) + c 2 × r a n d ( ) × ( g b e s t x i ) {v_i} = w \times {v_i} + {c_1} \times rand() \times (pbes{t_i} - {x_i}) + {c_2} \times rand() \times (gbest - {x_i})

其中 v i {v_i} 表示第i個粒子的方向, p b e s t i pbes{t_i} 表示第i個粒子當前最優引數, g b e s t gbest 表示當前的全域性最優引數, x i x_i 表示第i個粒子的引數值, r a n d ( ) rand() 表示介於 ( 0 , 1 ) (0,1) 之間的隨機數。

粒子位置更新方程為:

x i = x i + v i {x_i} = {x_i} + {v_i}

w w 取值越大說明粒子更新更多受到粒子歷史方向的影響, w = 0 w=0 表示粒子更新只取決於粒子當前位置,區域性最優粒子位置和全域性最優粒子位置,粒子方向本身沒有記憶性。 c 1 {c_1} 表示粒子更新受到區域性最優粒子的影響程度, c 1 = 0 {c_1}=0 表示粒子只有全域性搜尋能力,但缺少區域性搜尋能力。 c 2 {c_2} 表示粒子更新受到全域性最優粒子的影響程度, c 2 = 0 {c_2}=0 表示粒子只有區域性搜尋能力,但缺少全域性搜尋能力。

二、PSO演算法引數尋優的python實現

完整python程式碼和樣本地址:https://github.com/shiluqiang/PSO_python

本博文以非線性SVM為待優化模型,待優化引數為正則化引數 C C 和核引數 σ \sigma ,適應度函式值為3-fold交叉驗證平均值。

## 2. PSO優化演算法
class PSO(object):
    def __init__(self,particle_num,particle_dim,iter_num,c1,c2,w,max_value,min_value):
        '''引數初始化
        particle_num(int):粒子群的粒子數量
        particle_dim(int):粒子維度,對應待尋優引數的個數
        iter_num(int):最大迭代次數
        c1(float):區域性學習因子,表示粒子移動到該粒子歷史最優位置(pbest)的加速項的權重
        c2(float):全域性學習因子,表示粒子移動到所有粒子最優位置(gbest)的加速項的權重
        w(float):慣性因子,表示粒子之前運動方向在本次方向上的慣性
        max_value(float):引數的最大值
        min_value(float):引數的最小值
        '''
        self.particle_num = particle_num
        self.particle_dim = particle_dim
        self.iter_num = iter_num
        self.c1 = c1  ##通常設為2.0
        self.c2 = c2  ##通常設為2.0
        self.w = w    
        self.max_value = max_value
        self.min_value = min_value
        
        
### 2.1 粒子群初始化
    def swarm_origin(self):
        '''粒子群初始化
        input:self(object):PSO類
        output:particle_loc(list):粒子群位置列表
               particle_dir(list):粒子群方向列表
        '''
        particle_loc = []
        particle_dir = []
        for i in range(self.particle_num):
            tmp1 = []
            tmp2 = []
            for j in range(self.particle_dim):
                a = random.random()
                b = random.random()
                tmp1.append(a * (self.max_value - self.min_value) + self.min_value)
                tmp2.append(b)
            particle_loc.append(tmp1)
            particle_dir.append(tmp2)
        
        return particle_loc,particle_dir

## 2.2 計算適應度函式數值列表;初始化pbest_parameters和gbest_parameter   
    def fitness(self,particle_loc):
        '''計算適應度函式值
        input:self(object):PSOparticle_loc(list):粒子群位置列表
        output:fitness_value(list):適應度函式值列表
        '''
        fitness_value = []
        ### 1.適應度函式為RBF_SVM3_fold交叉校驗平均值
        for i in range(self.particle_num):
            rbf_svm = svm.SVC(kernel = 'rbf', C = particle_loc[i][0], gamma = particle_loc[i][1])
            cv_scores = cross_validation.cross_val_score(rbf_svm,trainX,trainY,cv =3,scoring = 'accuracy')
            fitness_value.append(cv_scores.mean())
        ### 2. 當前粒子群最優適應度函式值和對應的引數
        current_fitness = 0.0
        current_parameter = []
        for i in range(self.particle_num):
            if current_fitness < fitness_value[i]:
                current_fitness = fitness_value[i]
                current_parameter = particle_loc[i]

        return fitness_value,current_fitness,current_parameter 
        

## 2.3  粒子位置更新 
    def updata(self,particle_loc,particle_dir,gbest_parameter,pbest_parameters):
        '''粒子群位置更新
        input:self(object):PSOparticle_loc(list):粒子群位置列表
              particle_dir(list):粒子群方向列表
              gbest_parameter(list):全域性最優引數
              pbest_parameters(list):每個粒子的歷史最優值
        output:particle_loc(list):新的粒子群位置列表
               particle_dir(list):新的粒子群方向列表
        '''
        ## 1.計算新的量子群方向和粒子群位置
        for i in range(self.particle_num): 
            a1 = [x * self.w for x in particle_dir[i]]
            a2 = [y * self.c1 * random.random() for y in list(np.array(pbest_parameters[i]) - np.array(particle_loc[i]))]
            a3 = [z * self.c2 * random.random() for z in list(np.array(gbest_parameter) - np.array(particle_dir[i]))]
            particle_dir[i] = list(np.array(a1) + np.array(a2) + np.array(a3))
#            particle_dir[i] = self.w * particle_dir[i] + self.c1 * random.random() * (pbest_parameters[i] - particle_loc[i]) + self.c2 * random.random() * (gbest_parameter - particle_dir[i])
            particle_loc[i] = list(np.array(particle_loc[i]) + np.array(particle_dir[i]))
            
        ## 2.將更新後的量子位置引數固定在[min_value,max_value]內 
        ### 2.1 每個引數的取值列表
        parameter_list = []
        for i in range(self.particle_dim):
            tmp1 = []
            for j in range(self.particle_num):
                tmp1.append(particle_loc[j][i])
            parameter_list.append(tmp1)
        ### 2.2 每個引數取值的最大值、最小值、平均值   
        value = []
        for i in range(self.particle_dim):
            tmp2 = []
            tmp2.append(max(parameter_list[i]))
            tmp2.append(min(parameter_list[i]))
            value.append(tmp2)
        
        for i in range(self.particle_num):
            for j in range(self.particle_dim):
                particle_loc[i][j] = (particle_loc[i][j] - value[j][1])/(value[j][0] - value[j][1]) * (self.max_value - self.min_value) + self.min_value
                
        return particle_loc,particle_dir

## 2.4 畫出適應度函式值變化圖
    def plot(self,results):
        '''畫圖
        '''
        X = []
        Y = []
        for i in range(self.iter_num):
            X.append(i + 1)
            Y.append(results[i])
        plt.plot(X,Y)
        plt.xlabel('Number of iteration',size = 15)
        plt.ylabel('Value of CV',size = 15)
        plt.title('PSO_RBF_SVM parameter optimization')
        plt.show() 
        
## 2.5 主函式        
    def main(self):
        '''主函式
        '''
        results = []
        best_fitness = 0.0 
        ## 1、粒子群初始化
        particle_loc,particle_dir = self.swarm_origin()
        ## 2、初始化gbest_parameter、pbest_parameters、fitness_value列表
        ### 2.1 gbest_parameter
        gbest_parameter = []
        for