「粒子群算法」Machine Learning学习笔记（七）粒子群优化算法

粒子群算法

粒子群优化 算法( PARTICAL SWARMS OPTIMIZATION)

粒子群优化，是除了蚁群算法，鱼群算法之外的一种群体智能的优化算法，源自对鸟类捕食问题的研究。

鸟类捕食问题

假设区域里就只有一块食物（即通常优化问题中所讲的最优解），鸟群的任务是找到这个食物源。鸟群在整个搜寻的过程中，通过相互传递各自的信息，让其他的鸟知道自己的位置，通过这样的协作，来判断自己找到的是不是最优解，同时也将最优解的信息传递给整个鸟群，最终，整个鸟群都能聚集在食物源周围，即我们所说的找到了最优解，即问题收敛。

PSO算法概述

PSO算法首先在可行解空间中初始化一群粒子，每个粒子都代表极值优化问题的一个潜在最优解，用位置、速度和适应值三项指标表示该粒子特征。

粒子在解空间中运动，通过追踪个体极值Pbest和群体极值Gbest更新个体位置，个体极值是指个体所经历位置中计算得到的适应度值最优位置，群体极值是指种群中的所有粒子搜索到的适应度最优位置。

粒子每更新一次位置，就计算一次适应度值，并且通过比较新粒子的适应度值和个体极值、群体极值的适应度值更新个体极值和群体极值位置。

PSO算法流程

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。粒子群算法的思想相对比较简单，主要分为：1、初始化粒子群；2、评价粒子，即计算适应值；3、寻找个体极值；4、寻找全局最优解；5、修改粒子的速度和位置。下面是程序的流程图：

下面我们具体解释下流程图里面的每一个步骤：

1、初始化

首先，我们需要设置最大的速度区间，防止超出最大的区间。位置信息即为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置。设置群体规模。

2、个体极值与全局最优解

个体极值为每个粒子找到的历史上最优的位置信息，并从这些个体历史最优解中找到一个全局最优解，并与历史最优解比较，选出最佳的作为当前的历史最优解。

3、更新速度和位置的公式

更新公式为：

其中，称为惯性因子，和称为加速常数（学习率），一般取。表示区间上的随机数。表示第个变量的个体极值的第维。表示全局最优解的第维。

4、终止条件

有两种终止条件可以选择，一是最大代数：；二是相邻两代之间的偏差在一个指定的范围内即停止。我们在实验中选择第一种。

R语言实现粒子群优化

计算函数Z = sqrt(x^2+y^2)在x,y属于[-100,100]上的最小值

#绘制XYZ的3d图
library(rgl)  #用RGL包绘制三维交互式图形 
x <- (-100:100)
y <- (-100:100)
#mAPPly调用复合函数，byrow从行到列排
z=matrix(mapply(function(i){mapply(function(v0){return(sqrt(i^2+v0^2))},x)},y),nrow = 201,byrow = T) 
open3d()
surface3d(x,y,z,back = "lines",color = terrain.colors(z^2))
#该函数为锥形，在最低处，z取得最小值。此处，以求解该函数在定义域上的最小值为例，说明粒子群算法的实现过程
#粒子群初始化
#1.由于函数z有x和y两个输入变量，因此针对的是二维空间，在给定定义域的x,y属于【-100。100】上随机生成20个粒子，设置粒
#的最大速度为 30

#初始化粒子群（包括20个粒子）
vmax = 30
pbeast = NULL #历史经过的最合适位置
gbeast = NULL #种群经过的最合适位置
gbest.add = NULL 
w = 1 #设置惯性权重，通常取非负数，用于调节解空间的搜索范围，w = 1 时，算法为基本粒子群算法。w = 0 时，失去粒子对自身速度的记忆。
c1 = c2 =2 #设置加速度常数、
iters = 1000 #设置最大迭代次数
alpha = 0.001 #设置最佳适应度值的增量阈值
#--在给定定义域内，随机生成位置矩阵如下
set.seed(1)
xMat = matrix(c(x=runif(20,-100,100),y=runif(20,-100,100)),byrow = F,ncol = 2,dimnames = list(NULL,c("x","y")))
#--在给定的最大速度限制的条件下，随机生成速度矩阵
set.seed(1)
vMat = matrix(c(x=runif(20,-vmax,vmax),y=runif(20,-vmax,vmax)),byrow = F,ncol = 2,dimnames = list(NULL,c("x","y")))

#每个粒子的适应度
这里由于是求最小值，因此适应函数可以定义为一个增函数，求出对应增函数的最大值
adjusts = apply(xMat,1,function(v){1/sqrt(sum(v^2)+1)})

#更新pbest、gbest
#同时更新所有粒子的位置与速度
#pbest记录每个粒子历史的适应度最高位置
#gbest记录种群历史适应度的最高位置
#更新完成后要计算每个粒子的适应度，以进入循环，当达到迭代次数与最佳适应度小于0.0002时，算法结束
pbest = xMat
pbest = cbind(pbest,adjusts)
gbest = pbest[which.max(pbest[,3]),]
for(k in 1:iters)
{
  #---更新pbest
  #遍历adjusts,如果对应粒子的适应度是历史中最高的，则完成替换
mapply(function(no,adj){if(adj>pbest[no,3])
{pbest[no,] <<- c(xMat[no,],adj)}},1:length #<<-是全局赋值的意思
(adjusts),adjusts)  
 print(pbest)
#更新gbest

if(max(pbest[,3])>gbest[3]){
     gbest.add = max(pbest[,3])-gbest[3]
     gbest = pbest[which.max(pbest[,3]),]  
     print("--更新gbest")
     print(gbest)}       
#画去对应的位置点
plot(xMat[,1],xMat[,2],pch=20,col='blue',xlim=c(-100,100),ylim=c(-100,100))
points(gbest[1],gbest[2],pch=8,col='red')
points(0,0,pch=20,cex=0.5)
points(0,0,pch=21,cex=2)
dev.off()
#更新所有粒子的位置与速度
old.xMat<-xMat
xMat<-xMat+vMat
vMat<-w*vMat+c1*runif(1,0,1)*(pbest[,1:2]-old.xMat)+c2*runif(1,0,1)*
  (matrix(rep(gbest[1:2],20),ncol=2,byrow=T)-old.xMat)
#----如果vMat有值超过了边界值，则设定为边界值
vMat[vMat<(-vmax)]<-(-vmax)
vMat[vMat>vmax]<-vmax
#计算更新后种群中所有粒子的适应度
adjusts<-apply(xMat,1,function(v){1/sqrt(sum(v^2)+1)})
#检查全局适应度的增量，如果小于0.0002，则算法停止
if(!is.null(gbest.add) && gbest.add<0.0002){
  print(paste("k=",k,"算法结束！"))
  break;
} 
}


#粒子群算法的R语言实现

library(pso)
psoobj = psoptim(rep(NA,2),function(x)sqrt(x[1]^2+x[2]^2),lower = c(-100,-100),upper = c(100,100),control = list(s=50))
psoobj
psoptim(par, fn, gr = NULL, ..., lower = -1, upper = 1, control = list())

以Matlab实现粒子群优化

见https://blog.csdn.net/jim_cainiaoxiaolang/article/details/53156878

参考链接

https://blog.csdn.net/jim_cainiaoxiaolang/article/details/53156878

https://blog.csdn.net/dingming001/article/details/76222746

Machine Learning学习笔记（七）粒子群优化算法