R语言粒子群优化算法

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

绘制x，y，z的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 = 30pbeast = 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 = xMatpbest = 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)#更新gbestif(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<-xMatxMat<-xMat+vMatvMat<-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))psoobjpsoptim(par, fn, gr = NULL, ..., lower = -1, upper = 1, control = list())

对应参数解释

par 由优化问题的维度定义长度的向量，不需要在par向量中提供真实的值，用NA就够了。包含这个参数主要是为了保持。

与optim函数和兼容性，但是如果将par设置成lower和upper之间的值，那么第一个粒子的位置将会由par提供的位置决定
用于最小化（或者最大化函数），参数向量做为它的第一个参数，基于它实现了最小化。它应该返回一个标量结果、
…需要进一步传递给fn和gr参数
lower 所有变量的下界
upper 所有变量的上界
control 控制参数的列表，默认情况下，该函数使用粒子群算法处理最小化问题，如果设置参数control$fnscale为负
该函数便可以处理最大化问题。默认的control参数遵循了2007年提出的标准PSO标准，但是代码也对2011年提出的标准的PSO提供了支持，限制了最大速度，当所有粒子收敛到一个区域内时算法会重启，同时，使用BFGS来处理局部搜索问题该参数包含如下部分参数，通过list向主函数提供
trace:非负整数，如果是正的，就会跟踪优化过程中产生的信息，默认为0
fnscale:在优化的过程中，对fn和gr产生值的全面扩展。如果是负值，就将转化成最大值问题
maxit:迭代的最大数量，默认为1000
maxf:函数求值（没有考虑在数值梯度计算中执行的任何函数求值）的最大数量，默认为Inf
abstol:绝对收敛容忍度，该方法收敛一次所获得的最佳适应值小于或者等于abstol,默认为-Inf
s:粒子群规模