智能优化算法（源码）-生物地理学算法（Biogeography-Based Optimization，BBO)

获取更多资讯，赶快关注上面的公众号吧！

文章目录

介绍
生物地理学数学模型
基于生物地理学的优化
- 迁移
- 变异
- BBO定义和算法

生物地理学是研究生物有机体地理分布的科学，20世界60年代控制生物分布的数学方程首次被提出并得以发展。工程师的想法是向大自然学习，这咋一定程度上推动了生物地理学在优化问题上的应用。正如生物遗传学中的数学启发了遗传算法(GAs)的发展，生物神经元中的数学启发了人工神经网络的发展一样，Dan Simon认为生物地理学中数学也是发展一个新领域的基础，于是在2008年提出了基于生物地理学的优化(BBO)。扫码关注公众号，后台回复“生物地理”或“BBO”可以获取Matlab代码。

介绍

生物地理学的数学模型描述了物种如何从一个岛屿迁移到另一个岛屿，新物种如何出现，以及物种如何灭绝的过程。“岛”这个词在这里是描述性的，不能从字面上直接理解为“岛屿”，也就是说，岛屿是在地理上与其他栖息地隔离的任何栖息地。

非常适合生物物种居住的地理区域被认为具有高的栖息地适宜指数(HSI)。与HSI相关的特征包括降雨量、植被多样性、地形特征多样性、土地面积和温度等因素。表征可居住性的变量称为适宜性指数变量(SIVs)。SIVs可以看作是栖息地的自变量，HSI可以看作是因变量。HSI高的生境物种数量多，而HSI低的生境物种数量少。

简单地讲，HSI越高，该栖息地物种越多，整个环境也越趋于饱和，物种更多的迁出而不是迁入，反之，HSI越低，该栖息地物种越少，物种更多的迁入而不是迁出，物种的迁入就会提升HSI，而如果HSI持续保持在低水平，那么居住在那里的物种将趋向于灭绝，这将进一步为更多的迁入打开道路。因此，低HSI生境的物种分布比高HSI生境的物种分布更动态。这就类似于一般的求解过程，较优解就是HSI高的栖息地，较劣解就是HSI低的栖息地，高HSI的解更加稳定，相比于低HSI的解不会变化那么大，而是分享给低HSI更多特征，低HSI的解接受来自较优解的大量新的特征，这些特征的加入将提升解的质量，这就是整个BBO算法的核心原理。

生物地理学数学模型

图1为单一栖息地的物种丰富度模型，迁入率λ\lambdaλ和μ\muμ都是栖息地物种数量的函数。

考虑一下迁入曲线。当栖息地中没有物种时，可能的最大迁入率为III。随着物种数量的增加，栖息地变得更加拥挤，能够成功生存到栖息地的物种越来越少，迁入率下降。栖息地能支持的最大物种数量是SmaxS_{max}Smax，在该点迁入率变为零。

现在考虑迁出曲线。如果栖息地没有物种，那么迁出率当前为零。随着物种数量的增加，栖息地变得更加拥挤，更多的物种能够离开栖息地去探索其他可能的居住地，迁出率也随之增加。当物种达到栖息地所能承载的最大数量时，最大迁移率是EEE。

物种的平衡数是S0S_0S0，在该点，迁入率和迁出率相等，当然有一些扰动事件来打破这种平衡，比如自然灾害、疾病等。在图1中迁入率和迁出率都采用了直线，但是通常情况下它们都是更加复杂的曲线，采用直线的目的主要是为了更加清晰地表达迁入和迁出过程。

现在，考虑栖息地包括SSS种物种的概率为PsP_sPs,PsP_sPs从时间ttt到(t+Δt)(t+\Delta t)(t+Δt)的变化过程如下：

Ps(t+Δt)=Ps(t)(1−λsΔt−μsΔt)+Ps−1λs−1Δt+Ps+1μs+1Δt(1)\begin{aligned} P_{s}(t+\Delta t)=P_{s}(t)(1-&\left.\lambda_{s} \Delta t-\mu_{s} \Delta t\right) \\ &+P_{s-1} \lambda_{s-1} \Delta t+P_{s+1} \mu_{s+1} \Delta t \end{aligned}\tag{1} Ps(t+Δt)=Ps(t)(1−λsΔt−μsΔt)+Ps−1λs−1Δt+Ps+1μs+1Δt(1)
其中λs\lambda_{s}λs和μs\mu_{s}μs当栖息地中有SSS种物种时的迁入和迁出率。该等式成立的原因是为了保证在时间(t+Δt)(t+\Delta t)(t+Δt)有SSS种物种，以下条件必须满足其中之一：

在时间ttt有SSS种物种，且在ttt和(t+Δt)(t+\Delta t)(t+Δt)之间没有迁入或迁出发生；
在时间ttt有(S−1)(S-1)(S−1)种物种，且有一种物种迁入；
在时间ttt有(S+1)(S+1)(S+1)种物种，且有一种物种迁出。

假设Δt\Delta tΔt足够小，从而可以忽略多个迁入或迁出的概率。取Δt→0\Delta t \rightarrow 0Δt→0时式（1）的极限就得到了式（2）：
P˙s={−(λs+μs)Ps+μs+1Ps+1,S=0−(λs+μs)Ps+λs−1Ps−1+μs+1Ps+1,1≤S≤Smax⁡−1−(λs+μs)Ps+λs−1Ps−1S=Smax⁡(2)\dot{P}_{s}= \begin{cases}-\left(\lambda_{s}+\mu_{s}\right) P_{s}+\mu_{s+1} P_{s+1}, & S=0 \\ -\left(\lambda_{s}+\mu_{s}\right) P_{s}+\lambda_{s-1} P_{s-1}+\mu_{s+1} P_{s+1}, & 1 \leq S \leq S_{\max }-1 \\ -\left(\lambda_{s}+\mu_{s}\right) P_{s}+\lambda_{s-1} P_{s-1} & S=S_{\max }\end{cases}\tag{2} P˙s=⎩⎪⎨⎪⎧−(λs+μs)Ps+μs+1Ps+1,−(λs+μs)Ps+λs−1Ps−1+μs+1Ps+1,−(λs+μs)Ps+λs−1Ps−1S=01≤S≤Smax−1S=Smax(2)
定义n=Smax⁡n=S_{\max }n=Smax，P=[P0…Pn]TP=\left[\begin{array}{lll}P_{0} & \ldots & P_{n}\end{array}\right]^{T}P=[P0…Pn]T,可以将P˙s\dot{P}_{s}P˙s等式(对于 S=0,...,nS=0,...,nS=0,...,n)重组成矩阵形式：

P˙=AP(3)\dot{P}=A P\tag{3} P˙=AP(3)

其中矩阵AAA如式(4)所示。
A=[−(λ0+μ0)μ10…0λ0−(λ1+μ1)μ2⋱⋮⋮⋱⋱⋱⋮⋮⋱λn−2−(λn−1+μn−1)μn0…0λn−1−(λn+μn)](4)A=\left[\begin{array}{ccccc} -\left(\lambda_{0}+\mu_{0}\right) & \mu_{1} & 0 & \ldots & 0 \\ \lambda_{0} & -\left(\lambda_{1}+\mu_{1}\right) & \mu_{2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \lambda_{n-2} & -\left(\lambda_{n-1}+\mu_{n-1}\right) & \mu_{n} \\ 0 & \ldots & 0 & \lambda_{n-1} & -\left(\lambda_{n}+\mu_{n}\right) \end{array}\right]\tag{4} A=⎣⎢⎢⎢⎢⎢⎢⎢⎡−(λ0+μ0)λ0⋮⋮0μ1−(λ1+μ1)⋱⋱…0μ2⋱λn−20…⋱⋱−(λn−1+μn−1)λn−10⋮⋮μn−(λn+μn)⎦⎥⎥⎥⎥⎥⎥⎥⎤(4)
对于图1中的直线，有：
μk=Eknλk=I(1−kn).(5)\begin{aligned} &\mu_{k}=\frac{E k}{n} \\ &\lambda_{k}=I\left(1-\frac{k}{n}\right) . \end{aligned}\tag{5} μk=nEkλk=I(1−nk).(5)

考虑一种特殊的情况：E=IE=IE=I，这种情况下，有：
λk+μk=E(6)\lambda_{k}+\mu_{k}=E\tag{6} λk+μk=E(6)
且矩阵AAA变为：
A=E[−11n0…0nn−12n⋱⋮⋮⋱⋱⋱⋮⋮⋱2n−1nn0…01n−1]=EA′(7)\begin{aligned} A &=E\left[\begin{array}{ccccc} -1 & \frac{1}{n} & 0 & \ldots & 0 \\ \frac{n}{n} & -1 & \frac{2}{n} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \frac{2}{n} & -1 & \frac{n}{n} \\ 0 & \ldots & 0 & \frac{1}{n} & -1 \end{array}\right] \\ &=E A^{\prime} \end{aligned}\tag{7} A=E⎣⎢⎢⎢⎢⎢⎢⎢⎡−1nn⋮⋮0n1−1⋱⋱…0n2⋱n20…⋱⋱−1n10⋮⋮nn−1⎦⎥⎥⎥⎥⎥⎥⎥⎤=EA′(7)

1)观察1：000是A′A^{\prime}A′的特征值，特征向量为：

v=[v1…vn+1]Tvi={n!(n−1−i)!(i−1)!(i=1,…,i′)vn+2−i(i=i′+1,…,n+1)(8)\begin{aligned} v &=\left[\begin{array}{lll} v_{1} & \ldots & v_{n+1} \end{array}\right]^{T} \\ v_{i} &= \begin{cases}\frac{n !}{(n-1-i) !(i-1) !} & \left(i=1, \ldots, i^{\prime}\right) \\ v_{n+2-i} & \left(i=i^{\prime}+1, \ldots, n+1\right)\end{cases} \end{aligned}\tag{8} vvi=[v1…vn+1]T={(n−1−i)!(i−1)!n!vn+2−i(i=1,…,i′)(i=i′+1,…,n+1)(8)

其中i′i^{\prime}i′是大于等于(n+1)/2(n+1)/2(n+1)/2的最小整数，即i′=ceil⁡((n+1)/2)i^{\prime}=\operatorname{ceil}((n+1) / 2)i′=ceil((n+1)/2)。

这个观察可以通过特征值方程A′v=kvA^{\prime} v=k vA′v=kv的解来验证。例如，n=4n=4n=4时，有

v=[14641]T(9)v=\left[\begin{array}{lllll} 1 & 4 & 6 & 4 & 1 \end{array}\right]^{T}\tag{9} v=[14641]T(9)

n=5n=5n=5时，有：

v=[16101051]T(10)v=\left[\begin{array}{lllll} 1 & 6 & 10 & 10 & 5 & 1 \end{array}\right]^{T}\tag{10} v=[16101051]T(10)

2）推论1：A′A^{\prime}A′的特征值为：
k={0,−2n,−4n,…−2}(11)k=\left\{0, \quad \frac{-2}{n}, \quad \frac{-4}{n}, \quad \ldots \quad-2\right\}\tag{11} k={0,n−2,n−4,…−2}(11)

然该推论还没被证明，但是目前看对于已经测试过的所有nnn的值都是正确的。

定理1：各物种数目的概率的稳态值由下式给出：

P(∞)=v∑i=1n+1vi(12)P(\infty)=\frac{v}{\sum_{i=1}^{n+1} v_{i}}\tag{12} P(∞)=∑i=1n+1viv(12)

基于生物地理学的优化

本节将介绍如何将上述的生物地理学理论用于离散优化问题。其实上面的这一堆公式都可以不用关心其推导过程，只要能够理解什么是HSI、什么是SIV、什么是迁出率和什么是迁入率就够了。

迁移

解向量中的每个值可以认为是一个SIV，较优解对应高HSI的栖息地，较劣解对应低HSI的栖息地，所以HSI就是适应度值。假设每个解（栖息地）具有相同的物种曲线（简单起见，E=IE=IE=I），但是解所代表的SSS值取决于其HSI。图2中的S1S_1S1表示一个较劣解，栖息地具有较少的物种，迁入率也比较高，而S2S_2S2表示一个较优解，栖息地具有较多的物种，迁出率也比较高。

使用每个解的迁出率和迁入率来在栖息地之间概率地共享信息。以概率PmodP_{mod}Pmod，基于其他解更新每个解。如果一个给定的解被选择进行更新，那么就使用其迁入率λ\lambdaλ来概率性地决定是否更新解中每个适宜性指数变量SIV，即每个维度是不是要更新。如果选择了给定解SiS_iSi的某个SIVSIVSIV，那么就使用其他解的迁出率μ\muμ来概率性地决定哪个解迁移一个随机选择的SIV到解SiS_iSi中。

BBO的迁移策略类似于遗传算法和进化策略中的全局重组，但也有不同。在进化策略中，全局重组用于创造新的解，而BBO迁移却是用于更改现有解。

变异

灾难性事件可以彻底改变一个自然栖息地的HSI。它们还可能导致物种数量与平衡值不同。因此，栖息地的HSI可能会因为明显的随机事件而突然改变。我们在BBO中将其建模为SIV变异，并使用物种计数概率来确变异变率。

每种物种计数的概率由式(2)中给出的微分方程决定。通过观察图2物种曲线上的平衡点，我们可以看到低物种计数和高物种计数的概率都相对较低。这也可以从定理1推导出来。中等物种数有很高的概率，因为它们接近平衡点。

举一个例子，考虑Smax=10S_{max}=10Smax=10的情况，式（2）的稳态解和初始条件P(0)P(0)P(0)无关，可以通过数值计算或定理1计算得到，如式(13)。
P(∞)≈[0.0010.0010.0440.1170.2050.2460.2050.1170.0440.0010.001]T(13)P(\infty) \approx\left[\begin{array}{llllllllll} 0.001 & 0.001 & 0.044 & 0.117 & 0.205 & 0.246 & 0.205 & 0.117 & 0.044 & 0.001 & 0.001 \end{array}\right]^{T}\tag{13} P(∞)≈[0.0010.0010.0440.1170.2050.2460.2050.1170.0440.0010.001]T(13)

每个种群成员都有一个相关的概率，表明其作为给定问题解的可能性。非常高或非常低HSI的解都不太可能出现，中间值却比较可能出现。如果给定的解就有较低的概率PsP_sPs，那么其就可能变异成其他解。相反，一个具有较高概率的解就不太可能变异成其他解。这可以通过与解概率成反比的变异率mmm来实现：
m(S)=mmax⁡(1−PsPmax⁡)(14)m(S)=m_{\max }\left(\frac{1-P_{s}}{P_{\max }}\right)\tag{14} m(S)=mmax(Pmax1−Ps)(14)
其中mmaxm_{max}mmax是用户自定义参数，这种变异机制倾向于增加种群的多样性，同时使得低HSI的解进行质量提升，也使得高HSI的解优上加优。请注意，优化过程将使用精英主义方法来保存在BBO过程中具有最优解的栖息地的特征，所以即使变异破坏了其HSI，也可以在需要时进行恢复。

BBO定义和算法

定义1：栖息地H∈SIVmH \in \mathrm{SIV}^{m}H∈SIVm是一个mmm维的整数向量，表达给定问题的一个合理解。

定义2：适宜性指数变量SIV∈CS I V \in CSIV∈C是一个整数。

定义3：栖息地适宜度指数HSI: H→RH \rightarrow RH→R是解的优劣度量，即适应度函数。

定义4：生态系统HnH^nHn包括nnn个栖息地的组。

定义5：迁入率λ(HSI):R→R\lambda(\mathrm{HSI}): R \rightarrow Rλ(HSI):R→R是HSIHSIHSI的单调非增函数，λi\lambda_iλi正比于邻近栖息地的SIVs迁移至栖息地HiH_iHi的可能性。

定义6：迁出率μ(HSI):R→R\mu(\mathrm{HSI}): R \rightarrow Rμ(HSI):R→R是HSIHSIHSI的单调非减函数，μi\mu_iμi正比于栖息地HiH_iHi的SIVs迁移至邻近栖息地的可能性。

定义7：栖息地更新Ω(λ,μ):Hn→H\Omega(\lambda, \mu): H^{n} \rightarrow HΩ(λ,μ):Hn→H是基于生态系统HnH^nHn调整栖息地HHH的概率算子。HHH被改变的概率正比于其迁入率λ\lambdaλ，修改源来自于HjH_jHj的概率正比于迁出率μj\mu jμj。

栖息地改进可大致描述如下：

定义8：变异M(λ,μ):H→HM(\lambda, \mu): H \rightarrow HM(λ,μ):H→H是基于栖息地的存在先验概率而进行随机修正栖息地SIVs的概率算子。变异过程如下：

定义9：生态系统转移函数Ψ=\Psi=Ψ= (m,n,λ,μ,Ω,M):Hn→Hn(m, n, \lambda, \mu, \Omega, M): H^{n} \rightarrow H^{n}(m,n,λ,μ,Ω,M):Hn→Hn是一个6元数组，来迭代优化生态系统。生态系统转移函数可以写为：
Ψ=λn∘μn∘Ωn∘HSIn∘Mn∘HSIn(15)\Psi=\lambda^{n} \circ \mu^{n} \circ \Omega^{n} \circ \mathrm{HSI}^{n} \circ M^{n} \circ \mathrm{HSI}^{n}\tag{15} Ψ=λn∘μn∘Ωn∘HSIn∘Mn∘HSIn(15)

换句话说，生态系统转移函数从计算每个栖息地的迁入和迁出率开始。然后，对每个栖息地进行修改，然后重新计算HSI。最后，进行变异，然后再次计算每个栖息地的HSI。

定义10：BBO算法BBO=(I,Ψ,T)\mathrm{BBO}=(\mathcal{I}, \Psi, T)BBO=(I,Ψ,T)是一个3元组。I:∅→{Hn,HSIn}\mathcal{I}: \emptyset \rightarrow\left\{H^{n}, \mathrm{HSI}^{n}\right\}I:∅→{Hn,HSIn}是一个用于创建栖息地初始生态系统和计算对应的HSI的函数。Ψ\PsiΨ是生态系统转移函数，T:Hn→{T:H^{n} \rightarrow\{T:Hn→{ true, false }\}}是终止准则。

BBO算法如下：

初始化BBO参数。初始化最大物种计数SmaxS_{max}Smax，最大迁移率EEE和III，最大变异率mmaxm_{max}mmax，以及精英参数。
随机初始化一组栖息地，每个栖息地对应于一个可能解。
对每个栖息地，将HSIHSIHSI映射为物种数量SSS、迁入率λ\lambdaλ和迁出率μ\muμ。
概率地使用迁入和迁出率来改进每个非精英解，然后重新计算HSI。
对于每个栖息地，使用式(2)更新其物种计数的概率，然后基于概率对每个非精英栖息地进行变异，重新计算每个HSI。
跳转到3进行下次迭代，此循环可以在预定义的代数之后终止，或者在找到一个可接受的问题解决方案之后终止。