2. DBO
基本的蜣螂算法通过模拟蜣螂在自然界中的四种行为（滚动、产卵、觅食和偷窃）来执行种群位置更新。
2.1 滚动蜣螂
在自然界中，蜣螂必须通过太阳导航，使其球滚动的路线尽可能直线。方程(1)用于原始论文中更新滚动蜣螂的位置：
x i ( t + 1 ) = x i ( t ) + α ⋅ k ⋅ x i ( t − 1 ) + b ⋅ Δ x (1) x_i(t + 1) = x_i(t) + \alpha \cdot k \cdot x_i(t - 1) + b \cdot \Delta x \tag{1} xi(t+1)=xi(t)+α⋅k⋅xi(t−1)+b⋅Δx(1)
其中
Δ x = ∣ x i ( t ) − X w ∣ \Delta x = |x_i(t) - X^w| Δx=∣xi(t)−Xw∣
其中 t t t表示当前迭代次数， x i ( t ) x_i(t) xi(t)表示 i i i次迭代中蜣螂的位置。在原始文本中， α \alpha α表示蜣螂是否偏离原始方向， α \alpha α的值以1或-1的概率设置，1表示无偏差，-1表示偏差。 k ∈ ( 0 , 0.2 ] k \in (0, 0.2] k∈(0,0.2]表示缺陷系数。 b b b是一个常数，属于[0, 1]，在代码中赋值为0.3。 X w X^w Xw表示全局最差值。 Δ x \Delta x Δx用于模拟太阳照射，较大的 Δ x \Delta x Δx表示蜣螂离光源更远。在自然界中，当蜣螂遇到障碍物时，它会通过跳舞行为获得新的滚动方向。为了模拟这种情况，原始文章通过概率模拟蜣螂在滚动过程中是否遇到障碍物，当遇到障碍物时，原始文章使用切线函数获得新的滚动方向，以模拟蜣螂的跳舞行为。滚动蜣螂位置的更新变为方程(2)。
x i ( t + 1 ) = x i ( t ) + tan ( θ ) ∣ x i ( t ) − x i ( t − 1 ) ∣ (2) x_i(t + 1) = x_i(t) + \tan(\theta) |x_i(t) - x_i(t - 1)| \tag{2} xi(t+1)=xi(t)+tan(θ)∣xi(t)−xi(t−1)∣(2)
θ ∈ [ 0 , π ] \theta \in [0, \pi] θ∈[0,π]，当 θ = 0 , π / 2 \theta = 0, \pi / 2 θ=0,π/2和 π \pi π时，位置不更新。
2.2 产卵蜣螂
在自然界中，蜣螂选择一个安全区域进行产卵，为了模拟这种行为，原始论文提出了一种边界选择策略来模拟这个区域，定义如下：
L b ∗ = max ( X ∗ × ( 1 − R ) , L b ) (3) Lb^* = \max(X^* \times (1 - R), Lb) \tag{3} Lb∗=max(X∗×(1−R),Lb)(3)
及
U b ∗ = min ( X ∗ × ( 1 + R ) , U b ) Ub^* = \min(X^* \times (1 + R), Ub) Ub∗=min(X∗×(1+R),Ub)
L b ∗ Lb^* Lb∗和 U b ∗ Ub^* Ub∗表示产卵区域的上下界， X ∗ X^* X∗是当前局部最优， R = 1 − t / T max R = 1 - t / T_{\max} R=1−t/Tmax， T max T_{\max} Tmax表示最大迭代次数。一旦产卵蜣螂确定了最佳产卵区域，它就会在这个区域内产卵。在原始文本中，每只产卵蜣螂产卵区域是动态变化的，因此搜索当前最优解所在的区域可以保证，同时防止陷入局部最优。产卵蜣螂的位置更新由方程(4)给出。
x i ( t + 1 ) = X ∗ + b 1 × ( x i ( t ) − L b ∗ ) + b 2 × ( x i ( t ) − U b ∗ ) (4) x_i(t + 1) = X^* + b_1 \times (x_i(t) - Lb^*) + b_2 \times (x_i(t) - Ub^*) \tag{4} xi(t+1)=X∗+b1×(xi(t)−Lb∗)+b2×(xi(t)−Ub∗)(4)
其中 b 1 b_1 b1和 b 2 b_2 b2是大小为1×Dim的随机变量，Dim表示问题的维度。
2.3 觅食蜣螂
在自然界中，蜣螂觅食也会像产卵一样选择一个安全区域。原始文本选择重新定义该区域的具体定义为公式：
L b b = max ( X b × ( 1 − R ) , L b ) (5) Lb^b = \max(X^b \times (1 - R), Lb) \tag{5} Lbb=max(Xb×(1−R),Lb)(5)
及
U b b = min ( X b × ( 1 + R ) , U b ) Ub^b = \min(X^b \times (1 + R), Ub) Ubb=min(Xb×(1+R),Ub)
这里 X b X^b Xb表示最优全局位置， L b b Lb^b Lbb和 U b b Ub^b Ubb表示最佳觅食区域的上下界， L b Lb Lb和 U b Ub Ub表示问题求解的上下界，每次觅食行为的蜣螂对应觅食蜣螂和觅食位置的一个位置更新如下：
x i ( t + 1 ) = x i ( t ) + C 1 × ( x i ( t ) − L b b ) + C 2 × ( x i ( t ) − U b b ) (6) x_i(t + 1) = x_i(t) + C_1 \times (x_i(t) - Lb^b) + C_2 \times (x_i(t) - Ub^b) \tag{6} xi(t+1)=xi(t)+C1×(xi(t)−Lbb)+C2×(xi(t)−Ubb)(6)
C 1 C_1 C1是一个服从正态分布的随机数， C 2 C_2 C2是一个属于[0, 1]的向量，大小为1×Dim。
2.4 偷窃蜣螂
在自然界中，一些蜣螂选择从其他蜣螂那里偷粪球。为了模拟这种行为，原始文章使用最优全局位置 X b X^b Xb作为有争议粪球的位置。偷窃蜣螂的行为定义为偷窃蜣螂的位置更新，偷窃蜣螂的位置更新由以下方程给出：
x i ( t + 1 ) = X b + S ⋅ g ⋅ ( ∣ x i ( t ) − X ∗ ∣ + ∣ x i ( t ) − X b ∣ ) (7) x_i(t + 1) = X^b + S \cdot g \cdot (|x_i(t) - X^*| + |x_i(t) - X^b|) \tag{7} xi(t+1)=Xb+S⋅g⋅(∣xi(t)−X∗∣+∣xi(t)−Xb∣)(7)
S S S表示原始文本中值为0.5的常数， g g g是随机变量的大小，Dim表示问题的维度。在原始文本中，滚动蜣螂、繁殖蜣螂、觅食蜣螂和偷窃蜣螂的数量分别为6、6、7和11。
2.5 DBO算法实现步骤
3. 改进的蜣螂优化算法
3.1 动机
DBO算法在解决优化问题方面优于许多智能优化算法，但对于DBO算法来说，获得某些优化问题的理想最优解仍然是一个挑战，因此提出了一系列改进措施，以提高DBO算法在解决优化问题方面的准确性。
3.2 基于好点集的种群初始化策略
华罗庚和其他著名的中国数学家提出了好点集（Hua Luogeng, 1978）。假设存在这样一个集合
P n ( k ) = { ( { r 1 ( n ) ⋅ k } , { r 2 ( n ) ⋅ k } , . . . , { r s ( n ) ⋅ k } ) , 1 ≤ k ≤ n } P_n(k) = \left\{ \left( \left\{ r_1^{(n)} \cdot k \right\}, \left\{ r_2^{(n)} \cdot k \right\}, ..., \left\{ r_s^{(n)} \cdot k \right\} \right), 1 \leq k \leq n \right\} Pn(k)={({r1(n)⋅k},{r2(n)⋅k},...,{rs(n)⋅k}),1≤k≤n}
在欧几里得空间 H s H_s Hs中，维度为 s s s，其偏差满足 φ ( n ) = C ( r , ε ) n − 1 + ε \varphi(n) = C(r, \varepsilon)n^{-1 + \varepsilon} φ(n)=C(r,ε)n−1+ε，那么我们可以将这个集合称为好点集， r r r是一个好点。好点的值 r r r是
r l = { 2 cos ( 2 m π p ) } , 1 ≤ l ≤ n r^l = \left\{ 2 \cos \left( \frac{2m\pi}{p} \right) \right\}, 1 \leq l \leq n rl={2cos(p2mπ)},1≤l≤n
其中 p p p是满足 p ≥ 2 ⋅ s + 3 p \geq 2 \cdot s + 3 p≥2⋅s+3的最小素数。它可以映射到搜索空间，如方程(11)所示。其中 u p p e r j upper_j upperj是搜索空间的上界， l o w j low_j lowj是下界。
x i ( k ) = ( u p p e r j − l o w j ) [ P n ( k ) ] + l o w j (8) x_i(k) = (upper_j - low_j) \left[ P_n(k) \right] + low_j \tag{8} xi(k)=(upperj−lowj)[Pn(k)]+lowj(8)
图1显示了由好点集生成的500点的初始种群分布，图2显示了由随机策略生成的500点的初始种群分布。为了证明好点集生成的点集比随机生成的点集更均匀，二维坐标轴被划分为小的1 * 1单位方格，以查看每个小方格中的点是否均匀分布。从图3中，我们可以看到好点集生成的点比随机生成的点更均匀。因此，通过引入好点集，可以有效增强种群多样性，这有助于算法摆脱局部最优。
3.3 改进的收敛因子
从图4中，我们可以看到改进的边界收敛因子在早期阶段相对于原始文本中的收敛因子下降得更慢，这可以使最优觅食和产卵边界更大，从而获得更好的全局搜索能力，而在后期阶段改进的边界收敛因子相对于原始文本中的收敛因子下降得更快，从而加速算法的收敛能力。
R = ( cos ( π ∗ ( t T max ) + 1 ) ) ∗ 0.5 (9) R = \left( \cos \left( \pi * \left( \frac{t}{T_{\max}} \right) + 1 \right) \right) * 0.5 \tag{9} R=(cos(π∗(Tmaxt)+1))∗0.5(9)
图4. 收敛因子的比较。
3.4 平衡全局搜索和局部搜索
在原始论文中，觅食蜣螂和产卵蜣螂的数量保持不变。从公式中，我们可以看到产卵蜣螂为算法贡献了局部搜索能力，因为产卵蜣螂的每次更新位置都严格控制在最优产卵边界内，即当前最优位置的邻近性，从公式中，我们可以看到觅食蜣螂为算法贡献了全局搜索能力，因为觅食蜣螂的位置更新不仅限于最优觅食区域， C 1 C_1 C1是一个服从正态分布的随机变量，可能在迭代过程中有较大的步长。我们希望算法在早期阶段获得更好的全局搜索能力，而在后期阶段获得更好的局部搜索能力，这可以减少产卵蜣螂的数量，并在早期阶段增加觅食蜣螂的数量，在后期阶段减少觅食蜣螂的数量。因此，本文提出了两种种群，并设计了一个动态临界值公式方程(10)。
P = 0.2 + 2 ∗ t T max max ∗ 0.6 (10) P = 0.2 + \frac{2 * t}{T_{\max}} \max * 0.6 \tag{10} P=0.2+Tmax2∗tmax∗0.6(10)
在原始论文中，产卵区域的选择存在缺陷，例如当前局部最优位置 X ∗ X^* X∗在每个维度上都有一个分量为0，这将导致最优产卵区域的上界等于下界，即 L b ∗ Lb^* Lb∗和 U b ∗ Ub^* Ub∗在每个维度上都有一个分量为0，这将导致所有产卵蜣螂集中在一个点 X ∗ X^* X∗位置。为了解决这个问题，本文将最优产卵区域定义为以 X ∗ X^* X∗为中心的圆形区域，球的表达式为
( X − X 1 ∗ ) 2 + ( X − X 2 ∗ ) 2 + ⋯ + ( X − X n ∗ ) 2 = R 2 (X - X_1^*)^2 + (X - X_2^*)^2 + \cdots + (X - X_n^*)^2 = R^2 (X−X1∗)2+(X−X2∗)2+⋯+(X−Xn∗)2=R2
其中 R a d i u s = ( U b − L b ) / 2 ∗ R Radius = (Ub - Lb) / 2 * R Radius=(Ub−Lb)/2∗R， U b Ub Ub和 L b Lb Lb分别是优化问题的上下界。在算法的早期迭代中，Radius的值较大，蜣螂位置更新以更充分地探索解空间，以避免陷入局部最优，在算法的后期迭代中，Radius的值较小，蜣螂位置更新以更充分地探索当前局部最优附近的区域，从而加速算法的收敛。产卵蜣螂的位置更新变为：
X 1 ∗ ( t + 1 ) = X 1 ∗ + R a d i u s ∗ K ∗ cos θ 1 X_1^*(t + 1) = X_1^* + Radius * K * \cos \theta_1 X1∗(t+1)=X1∗+Radius∗K∗cosθ1
X 2 ∗ ( t + 1 ) = X 2 ∗ + R a d i u s ∗ K ∗ sin θ 1 cos θ 2 X_2^*(t + 1) = X_2^* + Radius * K * \sin \theta_1 \cos \theta_2 X2∗(t+1)=X2∗+Radius∗K∗sinθ1cosθ2
X 3 ∗ ( t + 1 ) = X 3 ∗ + R a d i u s ∗ K ∗ sin θ 1 sin θ 2 cos θ 3 X_3^*(t + 1) = X_3^* + Radius * K * \sin \theta_1 \sin \theta_2 \cos \theta_3 X3∗(t+1)=X3∗+Radius∗K∗sinθ1sinθ2cosθ3
⋯ \cdots ⋯
X n − 1 ∗ ( t + 1 ) = X n − 1 ∗ + R a d i u s ∗ K ∗ sin θ 1 sin θ 2 ⋯ sin θ n − 2 cos θ n − 1 X_{n-1}^*(t + 1) = X_{n-1}^* + Radius * K * \sin \theta_1 \sin \theta_2 \cdots \sin \theta_{n-2} \cos \theta_{n-1} Xn−1∗(t+1)=Xn−1∗+Radius∗K∗sinθ1sinθ2⋯sinθn−2cosθn−1
X n ∗ ( t + 1 ) = X n ∗ + R a d i u s ∗ K ∗ sin θ 1 sin θ 2 ⋯ sin θ n − 2 sin θ n − 1 X_n^*(t + 1) = X_n^* + Radius * K * \sin \theta_1 \sin \theta_2 \cdots \sin \theta_{n-2} \sin \theta_{n-1} Xn∗(t+1)=Xn∗+Radius∗K∗sinθ1sinθ2⋯sinθn−2sinθn−1
该公式是将n维球坐标转换为欧几里得坐标， K K K是[0, 1]中的随机数， θ 1 . . . θ n − 2 \theta_1 ... \theta_{n-2} θ1...θn−2是[0, 2 π 2\pi 2π]中的随机数。
3.5 基于量子计算的t分布变异
Shor在1994年提出的量子质因数分解算法在当时引起了巨大轰动，其时间复杂度仅为多项式级，因此量子计算越来越受到全世界学者的关注。量子启发式进化算法最初由Han和Kim (2000)开发，他们利用量子比特的特性与经典计算有很大不同，量子比特和传统二进制比特之间的最大区别在于它可以使用量子比特同时表示多个状态，这可以确保搜索空间的良好探索，即使初始点很少。Rui Xu使用量子比特特性优化循环交叉映射以初始化种群，从而获得更好的个体初始化并提高算法解决最优值的准确性。
在本文中，我们提出了一种使用量子计算和变异因子执行各种操作的新型变异策略，以跳出最优局部解 X b X^b Xb，使算法能够增加跳出最优局部解的能力。首先，本文介绍了量子比特，它们是量子计算中包含信息的对象。量子比特可以被分类为单量子比特、双量子比特和多量子比特，本文使用单量子比特。单个量子比特的状态由单位二维列向量 ∣ α ⟩ |\alpha\rangle ∣α⟩给出，即 α 2 + β 2 = 1 \alpha^2 + \beta^2 = 1 α2+β2=1，量子态向量在 [ 0 , 1 ] [0, 1] [0,1]中扮演特殊角色，因为任何量子态向量都可以写成这两个状态的和，即 ∣ ψ ⟩ = α ∣ 0 ⟩ + β ∣ 1 ⟩ |\psi\rangle = \alpha|0\rangle + \beta|1\rangle ∣ψ⟩=α∣0⟩+β∣1⟩，其中 ψ \psi ψ表示量子比特，上述量子比特的特性可以用来构建变异因子，首先建立由多个量子比特组成的量子空间来解决问题，即：
∣ Z R ⟩ = [ α 1 , α 2 , . . . , α dim ] |Z_R\rangle = [\alpha_1, \alpha_2, ..., \alpha_{\text{dim}}] ∣ZR⟩=[α1,α2,...,αdim]
∣ Z R ⟩ = [ β 1 , β 2 , . . . , β dim ] |Z_R\rangle = [\beta_1, \beta_2, ..., \beta_{\text{dim}}] ∣ZR⟩=[β1,β2,...,βdim]
然后最优位置 X b X^b Xb需要映射到量子空间，这在本文中使用以下方法完成：
α j = x j b ∥ x b ∥ (14) \alpha_j = \frac{x_j^b}{\|x^b\|} \tag{14} αj=∥xb∥xjb(14)
β j = P j ⋅ 1 − ( x j b ∥ x b ∥ ) 2 (15) \beta_j = P_j \cdot \sqrt{1 - \left( \frac{x_j^b}{\|x^b\|} \right)^2} \tag{15} βj=Pj⋅1−(∥xb∥xjb)2(15)
其中 x j b x_j^b xjb是最优全局位置 X b X^b Xb在第 j j j维度的值， ∥ x b ∥ \|x^b\| ∥xb∥是模型长度， P P P设置为1或-1，以便将点的多样性最大化，并将量子空间转换为解空间，执行变换后。
X temp = [ α 1 , α 2 , . . . , α dim ] ∗ ∥ X b ∥ X_{\text{temp}} = [\alpha_1, \alpha_2, ..., \alpha_{\text{dim}}] * \|X^b\| Xtemp=[α1,α2,...,αdim]∗∥Xb∥
X temp = [ β 1 , β 2 , . . . , β dim ] ∗ ∥ X b ∥ X_{\text{temp}} = [\beta_1, \beta_2, ..., \beta_{\text{dim}}] * \|X^b\| Xtemp=[β1,β2,...,βdim]∗∥Xb∥
为了更充分地探索解空间，本文使用量子旋转（Xiong et al., 2018）技术旋转 ∣ Z R ⟩ |Z_R\rangle ∣ZR⟩，这可以增加种群的多样性，从而避免陷入局部最优，如下所示：
R ( θ ) = [ cos θ − sin θ sin θ cos θ ] (16) R(\theta) = \begin{bmatrix} \cos \theta - \sin \theta \\ \sin \theta \cos \theta \end{bmatrix} \tag{16} R(θ)=[cosθ−sinθsinθcosθ](16)
R ( θ ) = R ( θ ) [ α 1 , α 2 , . . . , α dim β 1 , β 2 , . . . , β dim ] (17) R(\theta) = R(\theta) \begin{bmatrix} \alpha_1, \alpha_2, ..., \alpha_{\text{dim}} \\ \beta_1, \beta_2, ..., \beta_{\text{dim}} \end{bmatrix} \tag{17} R(θ)=R(θ)[α1,α2,...,αdimβ1,β2,...,βdim](17)
θ ∈ [ 0 , 2 π ] \theta \in [0, 2\pi] θ∈[0,2π]。然后变异通过以下方程(18)进行。
X b = X b + trnd ( t ) ⋅ X temp (18) X^b = X^b + \text{trnd}(t) \cdot X_{\text{temp}} \tag{18} Xb=Xb+trnd(t)⋅Xtemp(18)
X b = X b + trnd ( t ) ⋅ X Rtemp (19) X^b = X^b + \text{trnd}(t) \cdot X_{\text{Rtemp}} \tag{19} Xb=Xb+trnd(t)⋅XRtemp(19)
其中 X Rtemp X_{\text{Rtemp}} XRtemp和 X Rtemp X_{\text{Rtemp}} XRtemp表示旋转后映射到量子比特的解空间中的点， trnd ( Iter ) \text{trnd}(\text{Iter}) trnd(Iter)是服从t分布的方差因子，自由度为Iter（当前迭代次数）；个别位置通过精英选择进行更新。
t分布（Wu et al., 2023），也称为学生分布，是统计学中常见的抽样方法。其分布状态与自由度有关，自由度越小，曲线越平坦。当自由度为1时，t分布转换为柯西分布，当自由度趋于无穷大时，t分布成为标准正态分布。其概率密度函数为方程(19)。图5显示了对应于不同自由度的函数图像。在算法中，当执行变异操作时，如果变异因子为0，则可以使算法更专注于局部探索，反之亦然，专注于全局探索（Jianhua & Zhiheng, 2021）。为了便于观察，定义变异因子服从t分布时的值为[-1, 1]，这使得算法在局部探索阶段更关注局部探索，当变异因子服从t分布的值在[-∞, -1] ∪ [1, +∞]时，使法律更关注全局探索。
trnd ( x , n ) = Γ ( n + 1 2 ) n π × Γ ( n 2 ) ( 1 + x 2 n ) − n + 1 2 , − ∞ < x < + ∞ (19) \text{trnd}(x, n) = \frac{\Gamma \left( \frac{n + 1}{2} \right)}{\sqrt{n \pi} \times \Gamma \left( \frac{n}{2} \right)} \left( 1 + \frac{x^2}{n} \right)^{-\frac{n + 1}{2}}, \quad -\infty < x < +\infty \tag{19} trnd(x,n)=nπ×Γ(2n)Γ(2n+1)(1+nx2)−2n+1,−∞<x<+∞(19)
Γ ( α ) = ∫ 0 ∞ x α − 1 e − x d x (20) \Gamma(\alpha) = \int_0^\infty x^{\alpha - 1} e^{-x} dx \tag{20} Γ(α)=∫0∞xα−1e−xdx(20)
如图5所示，当变异因子服从自由度为0.1的t分布时，算法以0.1的概率导致算法进行局部探索，即，t分布函数在x = -1和x = 1处被t分布函数x = -1和x = 1处的值包围， trnd ( Iter ) \text{trnd}(\text{Iter}) trnd(Iter)是服从自由度为Iter的概率密度函数的t分布。算法进入全局搜索的概率为1 - P，根据图5，随着自由度的逐渐增加，导致算法进入局部探索的概率P也变得更加显著。本文希望算法在第一次迭代中更活跃地进行局部探索，以便在较大的搜索空间中获得更好的搜索空间，而在后期更活跃地进行局部探索，这可以增加算法的迭代速度并保证算法的收敛率，因此本文设计使用迭代次数。自由度的数量也用于确定自由度的值。
图5. 不同自由度的t分布图。
[1] J. K. Xue, B. Shen, Dung beetle optimizer: a new meta‑heuristic algorithm for global optimization, The Journal of Supercomputing. 79(7) (2022) 7305-7336.
[2] F. Zhu, G. S. Li, H. Tang, et al, Dung beetle optimization algorithm based on quantum computing and multi-strategy fusion for solving engineering problems, Expert Systems with Applications. 236 (2024) 121219.