第 3 章 混合算法

在第2章中说明过，在基本BAS算法及变体的基础上，做出了大的改动的算法在本章说明。

3.1 BSO

关于BSO，主要的参考资料为王糖糖同学的论文Beetle Swarm Optimization Algorithm:Theory and Application。作为一个捣鼓了两天算法的能源人，论文对我而言很是吃力。因为我此前只听说过粒子群算法，却从未了解过其原理。在囫囵看过几篇入门粒子群算法的博文后，再阅读上述的论文，就看到了很多有意思的地方。

首先，假设大家都阅读过第2章，对BAS有了一定的了解。那么，接下来就能看看，BSO是如何把BAS和PSO结合起来的。此处，只把我的认知写在此处，如有问题，还请大家指出。

3.1.1 粒子群算法流程

我们先说一下标准的粒子群算法的几个核心更新关系。

粒子在第 \(k+1\) 回合中，其速度与第 \(k\) 回合的速度关系如式(3.1):

\[\begin{equation} v_{i,k+1} = \omega v_{i,k} + c_1*rand_1*(pbest_i-x_{i,k})+c_2*rand_2*(gbest_k-x_{i,k}) \tag{3.1} \end{equation}\]

那么第 \(k+1\) 回合中粒子位置如式(3.2)更新。

\[\begin{equation} x_{i,k} = x_{i,k} + v_{i,k} \tag{3.2} \end{equation}\]

其中，\(i=1,2,\cdots,S\)，\(S\) 为设定的粒子数目。而\(k=1,2,\cdots,n\)，\(n\) 为设定的最大循环数。\(\omega\) 表示粒子速度的惯性，也就是本回合速度对下回合速度更新的影响。 \(c1,c2\)则是学习因子，\(pbest_i-x_{i,k}\) 表示粒子\(i\) 当前的位置和其历史最优的位置的差距，而 \(gbest_k-x_{i,k}\)则表示，粒子\(i\) 当前的位置和历史全局最优的粒子的位置（所有粒子在历史时间跨度中最好的位置）的差距。

那么，这样更新了粒子的位置之后，自然要计算这一回合的适应度，也就是目标函数的值。如果存在适应度比自身历史或全局历史更优，那么就相应地更新\(pbest\)和\(gbest\)。然后，就重复上述过程，直到设定的终止条件满足。

\(\omega\)动态变化，能取得好的效果，一般，它如式(3.3)所示变化。

\[\begin{equation} \omega = \omega_{max} - (\omega_{max}-\omega_{min})\frac{k}{n} \tag{3.3} \end{equation}\]

当然，如果 \(k\) 取0，也就是初始回合，那么\(\omega = \omega_{max}\)。如果 \(k\) 取\(n\) （也就是最终的回合）,那么\(\omega = \omega_{min}\)。所以，这也是个权重线性递减的过程。

3.1.2 与BAS的结合

我们假想一下，要把PSO和BAS结合，从什么方向入手呢？

作为门外汉，我可能会想，那就多放置几只天牛，让它们像粒子一样，信息共享，但是本身的更新寻优还按照BAS的来（后面会讲，这类似于BSAS，但不是BSO）。

不过呢，BSO给了一种很丰满和交错的方式，让BAS和PSO不是那么的泾渭分明。

首先，位置更新时，按照PSO和BAS的方式各自更新后加权得到新的位置。也就是，利用了BAS的触须方向和步长，同时也利用了粒子（在BSO中就是天牛）的速度。如式(3.4):

\[\begin{equation} X_{s}^{k+1} = X_{s}^{k} + \lambda V_{s}^k + (1 - \lambda)\xi_{s}^k \tag{3.4} \end{equation}\]

第一项 \(X_{s}^{k}\) 自然是表明粒子 \(s\) 在 \(k\) 回合的位置了，第二项 \(\lambda V_{s}^k\) 也就是PSO中的天牛速度乘以权重 \(\lambda\)。那么第三项中的 \(\xi_{s}^k\) 应该是BAS中的步长乘以方向向量才对。但是，方向并不是BAS中随机生成的方向。

我们看看 \(\xi_s^k\) 怎么得到。在文中，给出了式(3.5)。

\[\begin{equation} \xi_s^{k+1} = \delta^kV_s^ksign(f(X_{rs}^k) - f(X_{ls}^k)) \tag{3.5} \end{equation}\]

这就是我们熟悉的BAS更新的方式，不过，把随机生成的天牛方向，换为了粒子速度 \(V_s^k\)。在BSO中，并没有BAS中成方向的操作。有同学可能会问，既然没有随机方向，那么天牛的左右须坐标怎么来的呢？ 和式(3.5)中使用速度替代方向的操作相同，计算左右须坐标时，也是用的速度，来代替方向。即，如式(3.6)。

\[\begin{equation} \begin{split} X_{rs}^k&=X_{s}^k+V_s^k\frac{d}{2} \\ X_{ls}^k&=X_{s}^k-V_s^k\frac{d}{2} \end{split} \tag{3.6} \end{equation}\]

这样，就完成了BSO算法的更新过程。里面的速度，和粒子群算法的更新方式相同。

注明，原文中使用的上一回合的左右须坐标， 更新这回合的左右须坐标。但是根据代码来看，是文章存在了脚标的小错误，应该是用天牛坐标，更新须坐标。

王糖糖同学还给出了一个动态的（BAS）步长更新系数，如式(3.7):

\[\begin{equation} \eta = \omega_{\delta_1}(\frac{\omega_{\delta_1}}{\omega_{\delta_0}})^{1/(1+10k/n)} \tag{3.7} \end{equation}\]

默认参数为\(\omega_{\delta_1} = 0.4\)， \(\omega_{\delta_0} = 0.9\)。即，步长在回合初始，衰减最快，回合临终时，衰减最慢。

3.1.3 总结

此前，在写BSAS算法说明的时候，我把BSO加入到了待完成列表内。看了双方的英文名，是在是太巧了。BSAS是Beetle Swarm Antennae Search，而BSO是Beetle Swarm Optimization。彼时，我还不知晓PSO是啥含义。现在看来，这样的名称所包含的意义是千差万别的。后来，在复现BSO的过程中（此句为马后炮，毕竟是有了王糖糖同学提供的代码），也深刻地体会到两者的不同，下面我们来做一下罗列。

• 区别
  • BSAS：是BAS的衍生，移动和搜索的方式也很BAS。一言蔽之，BSAS就像天牛有了k对须而不是一对，在每回合，它必须要沿着这些须的方向都探索一遍。所谓的Swarm，是因为每回合有k只（虚拟的）天牛在探索这k对须的方向。
  • BSO：在粒子群算法的大框架下，由速度的更新（PSO）和步长的更新（BAS）共同决定下一步的天牛位置。这其中的Swarm，可以视作是真正独立的群体，它们间共享信息。
• 联系
  • BSO和BSAS保留着BAS的搜寻和移动方式
  • BSO和BSAS都是意图通过每回合搜索更多的位置（也可以视为方向），达到更优的效果。

3.1.4 算法调用说明

BSOoptim(fn, init = NULL, constr = NULL, 
         lower = c(-50, -50), upper = c(50, 50), n = 300, 
         s = floor(10 + 2 *sqrt(length(lower))), 
         w = c(0.9, 0.4), 
         w_vs = 0.4, 
         step = 10,
         step_w = c(0.9, 0.4), 
         c = 8, 
         v = c(-5.12, 5.12), 
         trace = T,
         seed = NULL, 
         pen = 1e+06)

在本节写函数调用，是因为函数内有部分参数和本节的公式变量名称有出入。上面的代码是函数的默认调用形式，其中\(s\)表示设定的天牛或者粒子数目，\(w\)也就是式(3.3)中提到的\(\omega_{max}\)和\(\omega_{min}\)。w_vs 是\(\lambda\)，也就是式(3.4)中天牛速度和移动步长之间的权重，直观地理解为weight between velocity and step-size。step和c仍然是表示步长和步长与须距离之比。而\(step_w\)为一个向量，表示的是式(3.7)中的\(\omega_{\delta_1} = 0.4\)， \(\omega_{\delta_0} = 0.9\)。

一个简单的例子如下：

library(rBAS)
mich <- function(x){
y1 <- -sin(x[1])*(sin((x[1]^2)/pi))^20
y2 <- -sin(x[2])*(sin((2*x[2]^2)/pi))^20
return(y1+y2)
}
result <-
 BSOoptim(fn = mich,
           init = NULL,
           lower = c(-6,0),
           upper = c(-1,2),
           n = 100,
           step = 5,
           s = 10,seed = 1, trace = F)
result$par; result$value

## [1] -4.965998  1.570796

## [1] -1.967851

3.2 bBAS

bBAS算法由阮月同学提出，主要针对于0-1规划问题。当然，对于实数空间上的其他优化问题也可以使用，但是会牺牲一部分的精度（后面会谈到分辨率问题）。

3.2.1 算法流程

step1: 根据计算（或者初始随机生成）的天牛坐标，计算天牛质心和两触须的目标函数值，取最小的作为下一回合的位置。其中，搜索距离 \(d\) 的更新方式与BAS算法相同。可参考式(2.2) 与式(2.4)。然后得到了当前最佳的天牛位置 \(x_{best}\)。

step2: 两个速度向量的计算（或者初始随机生成）；如式(3.8)所示，\(\overrightarrow{V_i^0}\) 表示天牛状态（在对应的二进制位上）由1转0的速度（或者称之为可能性，概率），\(\overrightarrow{V_i^1}\) 表示天牛由0转1的速度。

\[\begin{equation} V_i^c = \begin{cases} \overrightarrow{V_i^0}, & if\quad x_i = 1 \\ \overrightarrow{V_i^1}, & if\quad x_i = 0 \end{cases} \tag{3.8} \end{equation}\]

那么这两个速度分量该如何计算呢？ 假设按照step1得到的最优位置（二进制位表示）向量长度为 \(m+n\) ，有 \(m\) 个元素为1， \(n\) 个元素为0。那么，在 \(x_{best}=1\) 的这 \(m\) 个元素上，速度分量的更新方式如式(3.9)。

\[\begin{equation} \begin{split} \text{if}\quad & x_{best}^i = 1\quad \text{then}\\ &V_i^1 = wV_i^1+cr\\ &V_i^0 = wV_i^0-cr \end{split} \tag{3.9} \end{equation}\]

在 \(x_{best}=0\) 的这 \(n\) 个元素上，速度分量的更新方式如式(3.10)所示。

\[\begin{equation} \begin{split} \text{if}\quad & x_{best}^i = 0\quad \text{then}\\ &V_i^1 = wV_i^1-cr\\ &V_i^0 = wV_i^0+cr \end{split} \tag{3.10} \end{equation}\]

其中， \(w\) 是惯性项，为前一时刻速度分量的影响权重， \(c\) 为一个固定的常数， \(r\) 为 \((0,1)\) 区间的随机数。

step3: 有了速度分量后，如何更新 \(x_{best}\) 呢？ 和实数空间上的优化问题不同，如果直接把速度与对应的位置向量相加，这势必不符合会使得原有的位置向量脱离0-1的范畴。实际上，更新方式如式(3.11)所示。

\[\begin{equation} x_i(t+1)=\begin{cases} \overline{x_i}(t) & if\quad rand() \leq S(V_i)\\ x_i(t) & if\quad rand() > S(V_i) \end{cases} \tag{3.11} \end{equation}\]

当生成的服从均匀分布 \(U(0,1)\) 的随机数向量（长度与速度相同），某些元素小于等于状态转换概率 \(S(V_i)\) 时，那么下一时刻的位置等于在这些元素对应位置上的位置取非。反之，保持不变。例如，位置向量 \(x_t = (1,1,1,0,0,1,0)\) ，如果其 2,3,5元素上满足随机数小于等于转换概率条件，那么下一时刻的位置为 \(x_{t+1}=(1,0,0,0,1,1,0)\) 。

速度参与更新的方式，表现在计算状态转换概率 \(S(V_i)\) 上，类似于Sigmoid函数，该概率计算方式如式(3.12)所示。

\[\begin{equation} S(V_i) = \frac{1}{1+e^{-V_i}} \tag{3.12} \end{equation}\]

3.2.2 总结

bBAS算法与原始的BAS算法差异很大，具体表现在，已经不使用左右须作为类似梯度的指导下一步移动的方向的指标了，而是比较三者的大小，择优者作为最佳位置。

此外，不管是0-1问题，还是实数空间的问题，都会（通过扩大维度）转换为0-1问题来求解。因此，算法本身没有步长，速度也仅仅是用于求取0-1之间相互转换的概率。

最后，由于实数空间的优化问题会被乘以一个分辨率向量后，转为二进制向量。比如，区间(-2.048,2.048)内的优化问题，会被转换为(0,4096)的空间，然后再把正整数进行二进制的转换。这样会导致，bBAS算法的搜索空间粒度较大（粗糙）。

“(对实数空间的优化问题)bbas的效果是比bpso好点，但是肯定没有实数空间的算法好”

— 阮月

因此，如果优化参数不包含整型变量，或者简单地可以转为非整型变量（如压力容器例子约束的构造，见4.2.3.2节），则优先考虑其他的算法。

更细致与形象地讨论分辨率问题，在4.5.1节。