R 语言手写粒子群算法

粒子群算法简介

简要描述

• 场景：一群鸟在区域内随机搜索食物，区域内只有一块食物。

• 假设条件：所有鸟都不知道食物在哪里，但是知道当前位置离食物有多远。

• 策略：鸟群相互传递信息，让其他鸟知道自己的位置，来判断自己找到的是不是最优解，同时将最优解的信息传递给整个鸟群。

原理

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度 \(V\) 和位置 \(X\)，速度代表移动的快慢，位置代表移动的方向。每个粒子在搜索空间中单独地搜寻最优解，并将其记为当前个体极值 \(pbest\)，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前种群最优解 \(gbest\)，粒子群中的所有粒子根据自己找到的当前个体极值 \(pbest\) 和整个粒子群共享的当前全局最优解 \(gbest\) 来调整自己的速度和位置。

公式

PSO 初始化为一群随机粒子，然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个极值（pbest，gbest）来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。

速度更新公式：

\[ V_{i}^{k+1} = \omega \times V_{i}^{k} + c_{1} \times r_1 \times (pbest_{i}^{k}-X_{i}^{k}) + c_{2} \times r_2 \times (gbest^{k}-X_{i}^{k}) \quad \quad i = 1,\cdots,n \]

注：这只鸟第 \(k+1\) 步的速度 = 上一步自身的速度惯性 + 自我认知部分 + 社会认知部分

位置更新公式：

\[ X_{i}^{k+1}=X_{i}^{k}+V_{i}^{k} \quad \quad i=1,\cdots,n \]

注：这只鸟第 \(k+1\) 步所在的位置 = 第 \(k\) 步所在的位置 + 第 \(k\) 步所在的位置 \(\times\) 运动的时间（每一步运动的时间一般取 1）

其中：

• \(n\) 为粒子的个数；

• \(k\) 为迭代的轮数

• \(\omega\) 为惯性因子，数值范围为 0 到 1 之间；

• \(V_{i}^{k}\) 为第 \(k\) 次迭代时第 \(i\) 个粒子的速度；

• \(X_{i}^{k}\) 为第 \(k\) 次迭代时第 \(i\) 个粒子的位置；

• \(c_1\) 为粒子的个体学习因子，也称为个体加速因子，这个因子越大，鸟儿越倾向于飞往它自己曾去的食物量最多的地方；

• \(c_2\) 为粒子的社会学习因子，也称为社会加速因子，这个因子越大，鸟儿越倾向于飞往其他鸟儿（同伴们）曾去的食物量最多的地方；

• \(r_1\)、\(r_2\) 为 [0, 1] 上的随机数；

• \(f(X)\) 在位置 \(X\) 时的适应度值（一般取目标函数值）；

• \(pbest_{i}^{k}\) 为到第 \(k\) 次迭代为止，第 \(i\) 个粒子经过的最好的位置。最好的位置如何定义？以求解极小值为例，若 \(f(X_{i}^{k}) \leqslant f(pbest_{i}^{k-1})\)，则更新 \(pbest_{i}^{k}=X_{i}^{k}\)。

• \(gbest^{k}\) 为到第 \(k\) 次迭代为止，所有粒子经过的最好的位置，即 \(gbest^{k}=\mathop{\arg\min}\limits_{pbest}f(pbest_{i}^{k})\)；

• \(V_{max}\) 为大于 0 的数，是人工设置的 \(V_{i}\) 的最大值，若 \(V_{i} \geqslant V_{max}\)，则令 \(V_{i}=V_{max}\)，这么做是为了控制每次的步长，以防搜索步长过大错过最优解。

公式解释

\[ V_{i}^{k+1} = {\color{Blue} {\omega \times V_{i}^{k}}} + {\color{Red} {c_{1} \times r_1 \times (pbest_{i}^{k}-X_{i}^{k})}} + {\color{Green} {c_{2} \times r_2 \times (gbest^{k}-X_{i}^{k})}} \quad \quad i=1,\cdots,n \]

• 第一部分记忆项是惯性保持部分，粒子沿着当前的速度和方向惯性飞行，不偏移的趋势；

• 第二部分是自我认知部分，粒子受到自身历史最好位置的吸引力，有回到自身历史最好位置的意愿；

• 第三部分是社会认知部分，粒子处在一个种群中，种群中若有有更好的粒子，粒子受到最佳位置的吸引力，有去社会最佳位置的意愿。

参数对粒子群算法的影响

学习因子

\[ V_{i}^{k+1} = \omega V_{i}^{k} + c_{1} r_1 (pbest_{i}^{k}-X_{i}^{k}) + c_{2} r_2 (gbest^{k}-X_{i}^{k}) \]

学习因子都是非负常数，是调整局部最优值和全局最优值权重的参数。

• 若 \(c_{1}=0\)，粒子没有了认知能力，变成只有社会信息的模型，称为全局 PSO 算法。粒子具有拓展探索空间的能力，具有较快的收敛速度，但是缺少局部搜索能力，对于复杂问题的求解更容易落入局部最优；

• 若 \(c_{2}=0\)，粒子间没有社会信息，模型变为只有认知能力的模型，称为局部 PSO 算法。由于个体间没有信息交流，相当于多个粒子独立进行盲目搜索，收敛速度慢，因而得到最优解的可能性较少。

注：在最初提出粒子群算法的论文中指出，个体学习因子和社会学习因子取 2 比较合适。

惯性因子

\(\omega\) 对算法的收敛起到很大的作用，其值越大，粒子飞跃的范围就越广，更容易找到全局最优，但是也会错失局部搜寻的能力。

注：一般来说惯性权重取 0.9 - 1.2 是比较合适的，一般取 0.9 就行

种群数量

粒子群算法的最大特点就是速度快，因此初始种群取 50 - 1000 都是可以的，虽然初始种群越大收敛性会更好，不过太大了也会影响速度。

迭代次数

一般取 100 - 4000，迭代次数太少则解不稳定，太多则浪费时间。对于复杂问题，迭代次数可以相应地提高。

空间维数

粒子搜索的空间维数即为自变量的个数。比如求一元函数的最值就是一维空间，\(D\) 元函数的最值就是 \(D\) 维空间。

注：\(V_{i}\)，\(X_{i}\)，\(pbest_{i}\) 和 \(gbest\) 都是 \(D\) 维向量。

\[ \begin{aligned} V_{i}&=(v_{i1},v_{i2},\cdots,v_{iD})^\top \\ X_{i}&=(x_{i1},x_{i2},\cdots,x_{iD})^\top \\ pbest_{i}&=(p_{i1},p_{i2},\cdots,p_{iD})^\top \\ gbest_{i}&=(g_{1},g_{2},\cdots,g_{D})^\top \end{aligned} \]

位置限制

限制粒子搜索的空间，即自变量的取值范围，对于无约束问题此处可以省略。

速度限制

如果粒子飞行速度过快，则很可能直接飞过最优解位置；但是如果飞行速度过慢，则会使得收敛速度变慢。因此设置合理的速度限制很有必要。

注：参数的确定从来就不是绝对的，如果当前参数得到的结果并不理想的话，可以人为地更改某些参数。

PSO 算法基本流程

PSO 算法的优缺点

优点：

1. 收敛速度较快；

2. 需要调整的参数少，原理简单，容易实现，这是 PSO 的最大优点；

3. 协同搜索，同时利用个体局部信息和群体全局信息进行指导搜索

缺点:

1. 算法局部搜索能力较差，搜索精度不高；

2. 算法不能保证搜索到全局最优解，主要有一下两个方面：

   • 有时粒子群在俯冲过程中会错失全局最优解；

   • 应用 PSO 算法处理高维复杂问题是，算法可能过早收敛。

3. PSO 算法是一种概率算法

---

代码实现

先建立一个目标函数：

\[ f(x,y) = \frac{\sin(\sqrt{x^2+y^2})}{\sqrt{x^2+y^2}}+\exp{\{\frac{\cos(2\pi x)+\cos(2\pi y)}{2}}\}-2.71289 \]

可视化目标函数

使用基础函数 persp 画 3D 图：

Code

a <- b <- seq(-2, 2, length=101)

f <- function(a, b) {
  term1 <- sin(sqrt(a^2 + b^2)) / sqrt(a^2 + b^2)
  term2 <- exp(0.5*(cos(2*pi*a) + cos(2*pi*b)))
  term3 <- -2.71289
  
  term1 + term2 + term3
}

z <- outer(a, b, f)

persp(
  a, b, z,
  main="test function",
  col="lightblue",
  theta=30, phi=20,
  r=50, d=0.1,
  expand=0.5,
  # ltheta=90, lphi=120,
  # shade=0.75,
  ticktype="detailed",
  nticks=5
) # produces the 3-D plot

也可以使用 plot3D 包中的 persp3D 函数画 3D 图：

Code

library(plot3D)
persp3D(a, b, z, contour = TRUE, zlim= c(-4, 2), image = FALSE)

求解极大值

确定目标函数

Code

fn1 <- function (x) {
  # 将这个复杂的函数拆解成3部分
  term1 <- sin(sqrt(sum(x^2))) / sqrt(sum(x^2))
  term2 <- exp(0.5*sum(cos((2*pi*x))))
  term3 <- -2.71289
  
  # 将3个部分合并起来得到最终的函数
  term1 + term2 + term3
}
# 这是向量化的目标函数

也将此目标函数作为适应度函数

设置参数

Code

M <- 20 # 粒子个数
D <- 2 # 变量个数（维度）
iters <- 150 # 设置最大迭代次数
vmax <- 0.5 # 最大速度限制
bond <- 2 # 设置目标函数边界
pbeast <- NULL # 个体经过的最佳位置
gbeast <- NULL # 种群经过的最佳位置
gbestDelta <- NULL # 粒子适应度增量
g <- matrix(nrow = iters, ncol = D) # 创建一个空矩阵用以记录每次迭代产生的种群最优解
gDelta <- c() # 用以记录每一次迭代产生的粒子的适应度增量
w <- 0.9 # 设置惯性权重，通常取非负数，一般取0.9-1.2之间
c1 <- c2 <- 2 # 设置加速度常数，一般设为2
alpha <- 0.0001 # 设置最佳适应度值的增量阈值

初始化粒子速度和位置

Code

#--在给定定义域内，随机生成位置矩阵如下
set.seed(1)
x <- matrix(runif(M*D, -bond, bond), nrow = M, ncol = D, byrow = TRUE)

#--在给定的最大速度限制的条件下，随机生成速度矩阵
set.seed(1)
v <- matrix(runif(M*D, -vmax, vmax), nrow = M, ncol = D, byrow = TRUE)

初始化个体最优解和种群最优解

Code

pbest <- x # 个体最优
gbest <- pbest[which.max(fn1(x)), ] # 种群最优

迭代

Code

# 迭代循环
for (t in 1:iters) {
  
  # 每一个粒子的循环
  for (i in 1:M) {
    
    # 更新速度，有速度限制
    # 若更新后速度大于限制，则速度设为最大速度
    v_ = w*v[i, ] + c1*runif(1)*(pbest[i, ] - x[i, ]) + c2*runif(1)*(gbest - x[i, ])
    v[i, ] = v_
    if (sum(v_ > vmax) != 0) {
      v[which(v_ > vmax), ] = vmax
    } else if (sum(v_ < -vmax) != 0) {
      v[which(v_ < -vmax), ] = -vmax
    } else {next}
    
    # 更新位置，有位置限制
    # 若更新位置后超出边界，则不更新位置（不知道这样做对不对）
    x_ = x[i, ] + v[i, ]
    if (sum(abs(x_) < 2)) {
      x[i, ] = x_
    }
    
    # 更新pbest
    # 对于每个粒子，能使适应度函数最大的位置为当前粒子的最优解
    if (fn1(x[i, ]) > fn1(pbest[i, ])) {
      pbest[i, ] = x[i, ]
    }
    
    # 更新gbest
    # 对于所有粒子，能使适应度函数最大的位置为当前种群的最优解
    # 更新粒子适应度增量gbestDelta
    if (fn1(pbest[i, ]) > fn1(gbest)) {
      gbestDelta = abs(fn1(pbest[i, ]) - fn1(gbest)) #粒子群适应度变化
      gbest <- pbest[i, ]    
    }
    
    # 记录每一次迭代的gbest与gbestDelta
    g[t, ] <- gbest
    gDelta[t] <- gbestDelta
    
    # 粒子的停止条件
    # 若对于一个粒子来说，增量gbestDelta存在且小于给定的阈值，则可以认为该粒子到达了最优解
    if (!is.null(gbestDelta)) {
      if (gbestDelta < alpha) break
    }
  }  
}

返回结果

Code

gbest # 迭代完毕后的种群最优解

[1] -0.017776262  0.003713483

Code

fn1(gbest) # 目标函数最大值

[1] 0.9965124

绘图

Code

par(mfrow = c(2, 1))
plot(x = 1:iters, y = gDelta, type = "l", main = "适应度增量")
plot(x = 1:iters, y = apply(g, 1, fn1), type = "l", main = "目标函数最大值")

传统 PSO 算法的一些改进

惯性因子动态变化

权重因子递减：

由于较大的权重因子有利于跳出局部最小点，便于全局搜索，而较小的惯性因子则有利于对当前的搜索区域进行精确局部搜索，以利于算法收敛，因此针对 PSO 算法容易早熟以及算法后期易在全局最优解附近产生振荡现象，可以采用递减的权重。

最常用的一种线性权值递减：

\[ \begin{aligned} w(k) &= \frac{(w_{max}-w_{min}) \times (T-k)}{T} + w_{min} \\ &= w_{max}- \frac{k}{T} \times (w_{max}-w_{min}) \end{aligned} \]

其中 \(T\) 为总的迭代次数，\(k\) 为迭代轮数，\(w\) 随着迭代轮数的增加而减少，\(w_{max}\) 和 \(w_{min}\) 为惯性因子 \(w\) 的上限与下限。

收缩因子

常数收缩因子

\[ \begin{aligned} V_{i}^{k+1} &= \varphi \times [\omega V_{i}^{k} + c_{1} r_1 (pbest_{i}^{k}-X_{i}^{k}) + c_{2} r_2 (gbest^{k}-X_{i}^{k})] \\ \varphi &= \frac{2}{\lvert 2 - C- \sqrt {C^2 - 4C} \rvert} \\ C &= c_1 + c_2, \quad C > 4 \end{aligned} \]

注：此函数在 \(C\) 大于 4 之后呈现单调递减的趋势，所以 \(C\) 越大收缩得越厉害。

学习因子 c1 和 c2 决定了微粒本身经验信息和其他微粒的经验信息对微粒运行轨迹的影响，反映了微粒群之间的信息交流。设置较大的 c1 值，会使微粒过多地在局部范围内徘徊，而设置较大的 c2 值，则又会促使微粒过早收敛到局部最小值。为了有效地控制微粒的飞行速度，使算法达到全局探测与局部开采两者间的有效平衡，Clerc 构造了引入收缩因子的 PSO 模型。这种调整方法通过选取合适的参数，确保 PSO 算法的收敛性，同时可取消对速度的边界限制。

设置参数：

Code

M <- 20 # 粒子个数
D <- 2 # 变量个数（维度）
iters <- 150 # 最大迭代次数
c1 = 2.6 # 学习因子1
c2 = 1.7 # 学习因子2
c = c1 + c2 # 总加速度
vmax <- 0.5 # 最大速度限制
bond <- 2 # 设置目标函数边界
pbeast <- NULL # 个体经过的最佳位置
gbeast <- NULL # 种群经过的最佳位置
gbestDelta <- NULL # 粒子适应度增量
g <- matrix(0 , nrow = iters, ncol = D) # 用以记录每次迭代产生的种群最优解
gDelta <- c() # 用以记录每一次迭代产生的粒子的适应度增量
w <- 1 # 设置惯性权重，通常取非负数，一般取0.9-1.2之间
maxw = 1.2 # 最大惯性权重
minw = 0.9 # 最小惯性权重
w_ <- c() # 用以记录惯性权重随迭代次数地变化
a <- 0.79 # 收缩因子
alpha <- 0.0001 # 设置最佳适应度值的增量阈值

收缩因子：

Code

ss = 2/(abs(2 - c - sqrt(c^2 - 4*c)))

初始化个体最优解与种群最优解：

Code

#--在给定定义域内，随机生成位置矩阵如下
set.seed(1)
x <- matrix(runif(M*D, -bond, bond), nrow = M, ncol = D, byrow = TRUE)

#--在给定的最大速度限制的条件下，随机生成速度矩阵
set.seed(1)
v <- matrix(runif(M*D, -vmax, vmax), nrow = M, ncol = D, byrow = TRUE)

迭代部分：

Code

# 迭代循环
for (t in 1:iters) {
  
  # 每一个粒子的循环
  for (i in 1:M) {
    
    # 更新速度，有速度限制
    # 使用收缩因子后可以取消对速度地限制
    v[i, ] = ss*w*v[i, ] + c1*runif(1)*(pbest[i, ] - x[i, ]) + c2*runif(1)*(gbest - x[i, ])
     
    # 更新位置，有位置限制
    # 若更新位置后超出边界，则不更新位置（不知道这样做对不对）
    x_ = x[i, ] + v[i, ]
    if (sum(abs(x_) < 2)) {x[i, ] = x_}
    
    # 更新pbest
    # 对于每个粒子，能使适应度函数最大的位置为当前粒子的最优解
    if (fn1(x[i, ]) > fn1(pbest[i, ])) {pbest[i, ] = x[i, ]}
    
    # 更新gbest
    # 对于所有粒子，能使适应度函数最大的位置为当前种群的最优解
    # 更新粒子适应度增量gbestDelta
    if (fn1(pbest[i, ]) > fn1(gbest)) {
      gbestDelta = abs(fn1(pbest[i, ]) - fn1(gbest)) #粒子群适应度变化
      gbest <- pbest[i, ]    
    }
    
    # 记录每一次迭代的gbest与gbestDelta
    g[t, ] <- gbest
    gDelta[t] <- gbestDelta
    
    # 粒子的停止条件
    # 若对于一个粒子来说，增量gbestDelta存在且小于给定的阈值，则可以认为该粒子到达了最优解
    if (!is.null(gbestDelta)) {
      if (gbestDelta < alpha) break
    }
  }
  w = maxw - (t/iters)*(maxw - minw) # 惯性权重线性衰减
  w_[t] = w
}

返回结果：

Code

gbest # 迭代完毕后的种群最优解

[1] 2.216724e-04 3.064636e-05

Code

fn1(gbest) # 目标函数最大值

[1] 1.00539

绘图：

Code

plot(x = 1:iters, y = w_, type = "l", main = "惯性权重衰减")

Code

plot(x = 1:iters, y = gDelta, type = "l", main = "适应度增量")

Code

plot(x = 1:iters, y = apply(g, 1, fn1), type = "l", main = "目标函数最大值")

动态收缩因子

\[ \begin{aligned} V_{i}^{k+1} &= \varphi(k) \times [\omega V_{i}^{k} + c_{1} r_1 (pbest_{i}^{k}-X_{i}^{k}) + c_{2} r_2 (gbest^{k}-X_{i}^{k})] \\ \varphi(k) &= \frac{2 - k/T}{\lvert 2 - C- \sqrt {C^2 - 4C} \rvert} \\ C &= c_1 + c_2, \quad C > 4 \end{aligned} \]

设置参数：

Code

M <- 20 # 粒子个数
D <- 2 # 变量个数（维度）
iters <- 150 # 最大迭代次数
c1 = 2.6 # 学习因子1
c2 = 1.7 # 学习因子2
c = c1 + c2 # 总加速度
vmax <- 0.5 # 最大速度限制
bond <- 2 # 设置目标函数边界
pbeast <- NULL # 个体经过的最佳位置
gbeast <- NULL # 种群经过的最佳位置
gbestDelta <- NULL # 粒子适应度增量
g <- matrix(0 , nrow = iters, ncol = D) # 用以记录每次迭代产生的种群最优解
gDelta <- c() # 用以记录每一次迭代产生的粒子的适应度增量
w <- 1 # 设置惯性权重，通常取非负数，一般取0.9-1.2之间
maxw = 1.2 # 最大惯性权重
minw = 0.9 # 最小惯性权重
w_ <- c() # 用以记录惯性权重随迭代次数的变化
ss_ <- c() # 用以记录收缩因子随迭代次数的变化
a <- 0.79 # 收缩因子
alpha <- 0.0001 # 设置最佳适应度值的增量阈值

初始化个体最优解与种群最优解：

Code

#--在给定定义域内，随机生成位置矩阵如下
set.seed(1)
x <- matrix(runif(M*D, -bond, bond), nrow = M, ncol = D, byrow = TRUE)

#--在给定的最大速度限制的条件下，随机生成速度矩阵
set.seed(1)
v <- matrix(runif(M*D, -vmax, vmax), nrow = M, ncol = D, byrow = TRUE)

迭代部分：

Code

# 迭代循环
for (t in 1:iters) {

# 动态收缩因子放进循环内
ss = 2*(1 - t/iters)/(abs(2 - c - sqrt(c^2 - 4*c))) 
ss_[t] <- ss

  # 每一个粒子的循环
  for (i in 1:M) {
    
    # 更新速度，有速度限制
    # 使用收缩因子后可以取消对速度地限制
    v[i, ] = ss*w*v[i, ] + c1*runif(1)*(pbest[i, ] - x[i, ]) + c2*runif(1)*(gbest - x[i, ])
     
    # 更新位置，有位置限制
    # 若更新位置后超出边界，则不更新位置（不知道这样做对不对）
    x_ = x[i, ] + v[i, ]
    if (sum(abs(x_) < 2)) {x[i, ] = x_}
    
    # 更新pbest
    # 对于每个粒子，能使适应度函数最大的位置为当前粒子的最优解
    if (fn1(x[i, ]) > fn1(pbest[i, ])) {pbest[i, ] = x[i, ]}
    
    # 更新gbest
    # 对于所有粒子，能使适应度函数最大的位置为当前种群的最优解
    # 更新粒子适应度增量gbestDelta
    if (fn1(pbest[i, ]) > fn1(gbest)) {
      gbestDelta = abs(fn1(pbest[i, ]) - fn1(gbest)) #粒子群适应度变化
      gbest <- pbest[i, ]    
    }
    
    # 记录每一次迭代的gbest与gbestDelta
    g[t, ] <- gbest
    gDelta[t] <- gbestDelta
    
    # 粒子的停止条件
    # 若对于一个粒子来说，增量gbestDelta存在且小于给定的阈值，则可以认为该粒子到达了最优解
    if (!is.null(gbestDelta)) {
      if (gbestDelta < alpha) break
    }
  }

  # 惯性权重线性衰减
  w = maxw - (t/iters)*(maxw - minw) 
  w_[t] = w
}

返回结果：

Code

gbest # 迭代完毕后的种群最优解

[1]  1.164653e-05 -2.965302e-05

Code

fn1(gbest) # 目标函数最大值

[1] 1.005392

绘图：

Code

par(mfrow = c(2, 2))
plot(x = 1:iters, y = w_, type = "l", main = "惯性权重衰减")
plot(x= 1:iters, y = ss_, type = "l", main = "收缩因子衰减")
plot(x = 1:iters, y = gDelta, type = "l", main = "适应度增量")
plot(x = 1:iters, y = apply(g, 1, fn1), type = "l", main = "目标函数最大值")

R程序包 pso

参数

R 程序包 pso 中的 psoptim 可以实现粒子群优化算法，输入 ?psoptim 以查看具体用法。

该函数默认求解目标函数的最小值，若求最大值则需要一些处理，具体如下。

• par：是由优化问题的维度所定义长度的向量，不需要在 par 向量中提供真实的值，用 NA 就可以。设置这个参数主要是为了保持与 optim 函数和兼容性，但是如果将 par 设置成 lower 和 upper 之间的值，那么第一个粒子的位置将会由 par 决定。

• fn：用于最小化（或者最大化）的目标函数，参数向量做为它的第一个参数，基于它实现了最小化。它应该返回一个标量结果。

• gr：使用 BFGS 局部搜索所返回的梯度的函数。 如果为 NULL，将使用有限差分近似。

• ...：需要进一步传递给 fn 和 gr 的参数。

• lower：所有变量的下界。

• upper：所有变量的上界。

• control：控制参数的列表。

默认情况下，该函数使用粒子群算法处理最小化问题，如果设置参数 control$fnscale 为负则该函数便可以处理最大化问题。 默认的 control 参数遵循了 2007 年提出的标准 PSO 标准，但是代码也对 2011 年提出的标准的 PSO 提供了支持，限制了最大速度，当所有粒子收敛到一个区域内时算法会重启，同时，使用 BFGS 来处理局部搜索问题。

control参数是一个列表（list），可以提供以下任何组件：

• trace：非负整数，如果是正的，就会跟踪优化过程中产生的信息，默认为 0

• fnscale：为目标函数的尺度乘数，设置之后优化的目标函数变为 fn/fnscale。由于 psoptim 默认求目标函数最小值，所以可以将 fnscale 设置为负数（通常取 -1）即可将最小化问题转变为最大化问题。要注意的是返回的结果为 fn/fnscale 的最小值，还需要乘以 fnscale 才能得到原目标函数的最值。

• maxit：迭代的最大数量，默认为 1000。

• maxf：函数求值（没有考虑在数值梯度计算中执行的任何函数求值）的最大数量，默认为 Inf。

• abstol：绝对收敛容忍度，该方法收敛一次所获得的最佳适应值小于或者等于abstol，默认为 -Inf。

• reltol：重新启动的容忍度。 一旦最佳粒子与所有其他粒子之间的最大距离小于 reltol*d，算法就会重新启动。 默认为 0，即禁用重新启动检查。

• s：粒子群规模（粒子个数）。默认为 floor(10+2*sqrt(length(par))) 除非类型为”SPSO2011”，在这种情况下默认为 40。

• 其他

使用pso求解

Code

library(pso)
solu1 <- psoptim(
  rep(NA, 2), 
  fn1, 
  lower = c(-2, -2), 
  upper = c(2, 2), 
  control = list(
    fnscale = -1, 
    p = 1, 
    s = 20
  ) 
)

Code

solu1$par # 最大值点

[1] 3.143326e-09 1.461398e-10

Code

fn1(solu1$par) # 最大值

[1] 1.005392

Code

(-1)*solu1$value # 或者使用这种方法返回最大值

[1] 1.005392

参考

[1] https://blog.csdn.net/flyfish866/article/details/110207881

[2] https://zhuanlan.zhihu.com/p/398856271

[3] https://blog.csdn.net/qq_40886952/article/details/117387270?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1-117387270-blog-104594060.pc_relevant_default&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1-117387270-blog-104594060.pc_relevant_default&utm_relevant_index=2

[4] https://blog.csdn.net/dingming001/article/details/76222746

[5] https://blog.csdn.net/daaikuaichuan/article/details/81382794?spm=1001.2101.3001.6650.6&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7Edefault-6-81382794-blog-110207881.pc_relevant_default&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7Edefault-6-81382794-blog-110207881.pc_relevant_default&utm_relevant_index=9