粒子群算法的python實現

　　什么是粒子群算法

　　粒子群算法，也稱粒子群優化算法或鳥群覓食算法(Particle Swarm Optimization，PSO)。由J. Kennedy和R. C. Eberhart等人于1995年提出。其屬于進化算法的一種，也是從隨機解出發，通過迭代尋找最優解，其通過適應度來評價解的品質。

　　這種算法以其實現容易、精度高、收斂快等優點引起了學術界的重視，并且在解決實際問題中展示了其優越性。

　　求解過程

　　PSO通過模擬鳥群的捕食行為完成最優解的求取。

　　假設一群鳥在一個空間捕捉食物。在這個區域里只有一塊食物(對應著最優解)。所有的鳥都不知道食物在那里。但是它們可以判斷自身與食物的大致距離(通過fit值判斷與最優解的距離)。最簡單有效的方法就是搜尋目前離食物最近的鳥的周圍區域。

　　PSO中，問題的每個解都是搜索空間中的一只“鳥”。我們稱之為“粒子”。所有的粒子都有一個由被優化的函數決定的適應值(fitness value)，可以通過種群中適應度最高的粒子的位置和自身與自身最優情況的判斷完成位置的更新。

　　即，PSO 會初始化一群隨機粒子。利用迭代尋找最優解。在每一次迭代中，粒子通過跟蹤兩個"極值"來實現位置的更新。

　　1、粒子本身所找到的最優解，個體極值pbest。

　　2、整個種群目前找到的最優解，全局極值gBest。

　　每個粒子有一個重要的屬性，名為速度，該屬性同樣決定了它們更新的方向和更新的距離。

　　粒子們會追隨當前的最優粒子在解空間中搜索。

　　粒子群算法的求解偽代碼如下：

　　初始化粒子群

　　while 未達到最大迭代次數或最小loss：

　　for each_p in 粒子群:

　　計算適應度

　　如果某個粒子適應度高于歷史上的最佳適應度(pbest)

　　將該值設置為新的pbest

　　選擇所有粒子的最佳適配值的粒子作為gbest

　　for each_p in 粒子群:

　　根據方程(a)計算粒子速度

　　根據方程式(b)更新粒子位置

　　其中方程(a)為：

　　方程(b)為：

　　在方程中，v[i] 是第i個粒子的速度，w是慣性權重(有助于跳出局部最優解)，present[i] 是第i個粒子的當前位置，pbest[i]是第i個粒子的歷史最佳，gbest是全局最佳，rand () 是介于(0， 1)之間的隨機數。c1、c2 是學習因子，通常情況下，c1 = c2 = 2。

　　對于方程(a)而言：

　　代表利用粒子本身所找到的最優解更新自己的速度。

　　代表利用整個種群目前找到的最優解更新自己的速度。

　　實現代碼無錫看婦科的醫院 http://www.ytsgfk120.com/

　　這是一個求取二元一次方程y = -x^2 + 20*x + 10的最大值的例子。

　　import numpy as np

　　class PSO():

　　def __init__(self,pN,dim,max_iter,func):

　　self.w = 0.8 #慣性因子

　　self.c1 = 2 #自身認知因子

　　self.c2 = 2 #社會認知因子

　　self.r1 = 0.6 #自身認知學習率

　　self.r2 = 0.3 #社會認知學習率

　　self.pN = pN #粒子數量

　　self.dim = dim #搜索維度

　　self.max_iter = max_iter #最大迭代次數

　　self.X = np.zeros((self.pN,self.dim)) #初始粒子的位置和速度

　　self.V = np.zeros((self.pN,self.dim))

　　self.pbest = np.zeros((self.pN,self.dim),dtype = float) #粒子歷史最佳位置

　　self.gbest = np.zeros((1,self.dim),dtype = float) #全局最佳位置

　　self.p_bestfit = np.zeros(self.pN) #每個個體的歷史最佳適應值

　　self.fit = -1e15 #全局最佳適應值

　　self.func = func

　　def function(self,x):

　　return self.func(x)

　　def init_pop(self,): #初始化種群

　　for i in range(self.pN):

　　#初始化每一個粒子的位置和速度

　　self.X[i] = np.random.uniform(0,5,[1,self.dim])

　　self.V[i] = np.random.uniform(0,5,[1,self.dim])

　　self.pbest[i] = self.X[i] #初始化歷史最佳位置

　　self.p_bestfit[i] = self.function(self.X[i]) #得到對應的fit值

　　if(self.p_bestfit[i] > self.fit):

　　self.fit = self.p_bestfit[i]

　　self.gbest = self.X[i] #得到全局最佳

　　def update(self):

　　fitness = []

　　for _ in range(self.max_iter):

　　for i in range(self.pN): #更新gbest\pbest

　　temp = self.function(self.X[i]) #獲得當前位置的適應值

　　if( temp > self.p_bestfit[i] ): #更新個體最優

　　self.p_bestfit[i] = temp

　　self.pbest[i] = self.X[i]

　　if(self.p_bestfit[i] > self.fit): #更新全局最優

　　self.gbest = self.X[i]

　　self.fit = self.p_bestfit[i]

　　for i in range(self.pN): #更新權重

　　self.V[i] = self.w*self.V[i] + self.c1*self.r1*(self.pbest[i] - self.X[i]) + \

　　self.c2*self.r2*(self.gbest - self.X[i])

　　self.X[i] = self.X[i] + self.V[i]

　　fitness.append(self.fit)

　　return self.gbest,self.fit

　　def count_func(x):

　　y = -x**2 + 20*x + 10

　　return y

　　pso_example = PSO(pN = 50,dim = 1,max_iter = 300, func = count_func)

　　pso_example.init_pop()

　　x_best,fit_best= pso_example.update()

　　print(x_best,fit_best)