進化演算法之粒子群（PSO）
概述：
正如其他的那些進化演算法一樣，粒子群演算法的靈感同樣來自於大自然。它由Eberhart和Kennedy共同提出的，其基本思想來自於他們早期對許多鳥類的群體行為進行建模模擬研究結果的啟發。
假設在一塊廣袤無垠的棲息地上有一群自由自在的鳥兒和一堆豐盛的食物，可是鳥兒卻不知道食物在哪，現在我們的目標是讓所有的鳥兒找到這堆美食，那該怎麼辦呢？Kennedy認為鳥之間存在著相互交換資訊，於是他們在模擬中添加了一些內容：每個個體（鳥兒）能夠通過一定規則估計自身位子的適應值，並能記住自己當前所找到的最好位置，記作：“區域性最優pBest”。此外還應記住整個群裡中所有鳥兒找到的最優位置，記作：“全域性最優gBest”。
基本數學概念：
設
在每次迭代中，粒子根據以下公式進行更新速度和位置：
其中i = 1,2,3…m；d = 1,2,3…D；k是迭代次數；r1和r2為[ 0 ，1]之間的隨機數，為了保持群體的多樣性；c1和c2為學習因子，也稱加速因子，其使粒子具有自我總結和向群體中優秀個體學習的能力，從而向自己的歷史最優點以及群體內最優點靠近，這兩個引數對粒子群演算法的收斂起的作用不是很大，但適當調整這兩個引數可以減少區域性最優的困擾。
粒子群演算法的速度公式右邊包含了3個部分：第一部分是粒子之前的速度向量，第二部分是粒子向自身最優位置靠近的速度向量，第三部分為粒子向群體最優位置靠近的速度向量。如果沒了後面兩部分，粒子將會保持相同的速度向量朝一個方向飛行，這樣粒子有極大的可能找不到最優解，如果沒了第一部分時粒子的飛行速度將僅僅由它們當前位置和歷史最優位子決定，則速度自身是無記憶的，且整個群體沒有向空間拓展的功能，除非全域性最優解在初始群體的空間內，否則無法搜尋到全域性最優。
演算法流程：
1）開始
2）設定最大迭代次數（這裡迭代次數人為設定，也可以自行根據條件來跳出迴圈）
3）設定群體個數（大群體搜尋快計算慢，小群體搜尋慢計算快）
4）初始化每個粒子的初始位子z
5）初始化每個粒子的初始速度v
6）計算每個粒子的初始位置的適應值fun（z），通過比較粒子間的適應值來選擇初始群體最優位置gbest
7）將每個粒子各自的初始位置賦值給它們的pbest
8）k=0
9）k=k+1
10）根據速度更新公式更新每個粒子的新速度v
11）根據位置更新公式更新每個粒子的新位置z
12）每個粒子通過對比新舊位置更新pbest
13）更新粒子群體的最優位置gbest
14）從（9）開始迭代，直到達到迭代次數後跳出迴圈
演算法效果展示：
輸入函式f(x) = -1/2*((x1-50)(x1-50)+(x2-25)(x2-25)) , 我們很容易知道其最值為0，分別為x1=50和x2=25的時候：
輸入函式f(x) = -1/2*(x1*x1+x2*x2+x3*x3) , 我們很容易知道其最值為0，分別為x1=0,x2=0,x3=0的時候：
R程式碼
PSO<-function(circle = 1000,biomass = 1000,c1,c2,position_range,vector_range,FUN,purpose = 'max'){
dimension<-nrow(position_range);#計算維度
position <- matrix(nrow = biomass,ncol = dimension)#初始化位置
for (i in 1:biomass) {
position[i,]<-apply(position_range, 1, FUN = function(x){
return( runif(n = 1,min = x[1],max= x[2]))
})
}
vector <- matrix(nrow = biomass,ncol = dimension)##初始化速度
for (i in 1:biomass) {
vector[i,]<-apply(vector_range, 1, FUN = function(x){
return( runif(n = 1,min = x[1],max= x[2]))
})
}
pbest<-position;#初始化每個粒子最優點
fun<-match.fun(FUN)
adopter<-apply(position, 1, FUN = fun)#計算全域性最優點
gbest<-matrix(data = rep(x = position[which(adopter == match.fun(purpose)(adopter,na.rm = T))[1],],biomass),nrow = biomass,byrow = T)
for (i in 1:circle) {
vector<-vector+c1*runif(n = 1,min = 0,max = 1)*(pbest-position)+c2*runif(n = 1,min = 0,max = 1)*(gbest-position)#更新速度
temp<-position+vector#重新計算位置
for (i in 1:dimension) {#控制粒子的位置和速度在域內
temp[which(temp[,i]<position_range[i,1]),]<-position_range[i,1]
temp[which(temp[,i]>position_range[i,2]),]<-position_range[i,2]
vector[which(vector[,i]<vector_range[i,1]),]<-vector_range[i,1]
vector[which(vector[,i]>vector_range[i,2]),]<-vector_range[i,2]
}
if (purpose == 'max')#更新每個粒子的最優點
c<-which(apply(temp, 1, FUN =fun)>apply(position, 1, FUN =fun))
else if (purpose == 'min')
c<-which(apply(temp, 1, FUN =fun)<apply(position, 1, FUN =fun))
pbest[c,]<-temp[c,]
position<-temp#更新位置
adopter<-apply(position, 1, FUN = fun)#更新全域性最優點
gbest<-matrix(data = rep(x = position[which(adopter == match.fun(purpose)(adopter,na.rm = T))[1],],biomass),nrow = biomass,byrow = T)
}
return(list(best = gbest[1,],value = fun(gbest[1,])))
}