关闭

R语言粒子群优化算法

标签: r语言粒子群算法
455人阅读 评论(0) 收藏 举报
分类:

计算函数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 = 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())
对应参数解释

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:粒子群规模

0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:49255次
    • 积分:1366
    • 等级:
    • 排名:千里之外
    • 原创:89篇
    • 转载:0篇
    • 译文:0篇
    • 评论:3条
    文章存档
    最新评论