蜜獾算法（HBA，Honey Badger Algorithm）

2021年由Hashim等人提出（論文：Honey Badger Algorithm: A New Metaheuristic Algorithm for Solving Optimization Problems）。模擬蜜獾在自然界中的智能捕食行為，屬于群體智能優化算法（與粒子群PSO、遺傳算法GA同屬一類），適用于連續優化問題。

一行為模型

蜜獾的兩種主要行為驅動算法設計：

（1）挖掘模式：

針對地下獵物：通過嗅覺定位 + 動態挖掘路徑。

（2）采蜜模式：

跟隨導蜜鳥定位蜂巢 → 直接奔襲蜜源。

二??數學模型與迭代步驟

2.1關鍵變量定義：

蜜獾位置：候選解向量 $x_i$

獵物強度（I）：

$I_i = \frac{?{F(x_i) - \min(F)}}{\max(F) - \min(F)} + \epsilon$

$F(x_i)$ ：個體 $x_i$ 的目標函數值（如適應度）

$F_{\min},F_{\max}$ ?：當前種群中最小和最大函數值，兩者做差可以將 $I$ 進行歸一化處理。

$\epsilon$ ：極小常數（防止分母為零）

上述公式代標強度的范圍最大是1+ $\epsilon$ ，最小是? $\epsilon$ 。

函數值 $F(x_i)$ 反映解的質量。優化問題中， $F(x_i)$ 越小（化最小問題）表示解越優，更優解代表更接近全局最優解，代表與“獵物”的隱含距離更小，因此隱含距離越小，強度越小，接近程度越大。

個體	函數值 F(x)F(x)F(x)	與最優解距離	計算強度 IiI_iIi?	強度 vs 距離關系
A（最優）	$F_{\min}=5$	最近	$I_A=\frac{5-5}{20-5}+0.01=0.01$	強度最小
B	$F=15$	中等	$I_B=\frac{15-5}{20-5}+0.01\approx 0.67$	強度中等
C（最差）	$F_{\max}=20$	最遠	$I_C=\frac{20-5}{20-5}+0.01=1.01$	強度最大

強度越小越好。

距離大???F(xi)差???(F(xi)?Fmin?)?↑??? $I_{i}$ ↑???算法反饋：加大搜索步長

這種設計以極低計算代價實現了蜜獾“嗅到遠距離獵物時采取更大動作”的生物智能模擬。

為什么不用幾何距離？

避免高維距離計算、函數值差距能推廣到離散/非幾何空間，泛用性高、算法運行中已計算 $F(x_i)$ ，復用數據無需額外開銷，效率高。

2.2探索階段（Exploration）- 氣味擴散

基于獵物強度的全局搜索：

$\vec{x}_{new} = \vec{x}_{prey} + F \cdot \vec{r}_1 \cdot \alpha \cdot I \cdot | \cos(2\pi r_2) \cdot [1 - \cos(2\pi r_3)] |$

$\vec{x}_{prey}$ ：當前最佳獵物位置

$F$ ：控制搜索方向（±1隨機切換）

$r_1, r_2, r_3$ ?：隨機向量

$\alpha$ ：密度因子（隨迭代遞減）

強度 $I$ 的作用：

$I_{i}$ 大（解質量差、距離遠）→ 乘以 $I_{i}$ ? 后步長增大 → 加強全局探索

$I_{i}$ 小（解質量好、距離近）→ 乘以 $I_{i}$ 后步長減小 → 精細局部開發

2.3開發階段（Exploitation）- 精確捕食

局部挖掘（模擬蜜獾洞穴內捕獵）：

$\vec{x}_{new} = \vec{x}_{prey} + F \cdot \vec{r}_4 \cdot \mathit{\alpha} \cdot | \vec{x}_{prey} - \vec{x}_{i} |$

$\vec{r}_4$ ?：隨機向量

$| \vec{x}_{prey} - \vec{x}_{i} |$ ：個體到獵物的距離

通過密度因子 $\alpha$ 動態平衡探索與開發：

$\alpha = C \cdot \exp\left(-\frac{t}{T_{max}}\right)$

（ $C$ 為常數， $t$ 為當前迭代， $T_{max}$ ? 為最大迭代）

決策規則：若 $rand < \delta$ （ $\delta$ 為切換閾值），使用挖掘模式；否則使用采蜜模式。

三實現方法

3.1過程示例

初始化種群位置 x_i, i=1,2,...,N
計算適應度值 F(x_i)
while t < T_max:更新密度因子 α計算獵物強度 I_ifor each 蜜獾個體:if rand < δ:執行探索階段（氣味擴散）else:執行開發階段（挖掘捕食）更新當前位置 x_i計算新適應度 F(x_i)更新全局最優解 x_preyt = t + 1
返回最優解 x_prey

3.2參數設置

參數	意義	推薦值
$N$	種群規模	30~50
$\delta$	模式切換閾值	0.8
$C$	密度衰減常數	2.0
$T_{max}$	最大迭代次數	100~500

3.3優勢和局限性

優勢	局限性
全局探索能力強（氣味機制）	對高維問題敏感
局部開發高效（動態密度因子）	參數需經驗調整（如δ）
收斂速度快于傳統算法（PSO等）	易陷局部最優的改進變體
代碼實現簡單

四應用場景?

工程優化：結構設計、PID控制器調參

人工智能：神經網絡超參數優化

能源管理：光伏陣列最大功率點跟蹤（MPPT）

研究方向：離散化改進（如HBA-TSP）、多目標版本（MO-HBA）

五?在CEC2017測試函數集上使用HBA的python代碼

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocatorclass HBA:def __init__(self, obj_func, dim=30, pop_size=100, max_iter=1000, lb=-100, ub=100, C=2, delta=0.8):"""蜜獾算法 (Honey Badger Algorithm)參數:obj_func -- 目標函數dim -- 問題維度pop_size -- 種群大小max_iter -- 最大迭代次數lb -- 變量下界ub -- 變量上界C -- 密度因子常數delta -- 模式切換閾值"""self.obj_func = obj_funcself.dim = dimself.pop_size = pop_sizeself.max_iter = max_iterself.lb = lbself.ub = ubself.C = Cself.delta = delta# 初始化種群self.positions = np.random.uniform(lb, ub, (pop_size, dim))self.fitness = np.zeros(pop_size)self.best_position = np.zeros(dim)self.best_fitness = np.infself.convergence_curve = np.zeros(max_iter)def evaluate(self):"""評估所有個體適應度"""for i in range(self.pop_size):self.fitness[i] = self.obj_func(self.positions[i])if self.fitness[i] < self.best_fitness:self.best_fitness = self.fitness[i]self.best_position = self.positions[i].copy()def optimize(self):"""執行優化過程"""self.evaluate()  # 初始評估for t in range(self.max_iter):# 1. 計算獵物強度f_min = np.min(self.fitness)f_max = np.max(self.fitness)I = (self.fitness - f_min) / (f_max - f_min + 1e-25)  # 避免除零I = 1 + I  # 確保強度至少為1# 2. 更新密度因子alpha = self.C * np.exp(-t / self.max_iter)# 3. 更新每個個體位置for i in range(self.pop_size):# 為每個個體生成隨機方向因子F = 1 if np.random.rand() < 0.5 else -1  # 隨機方向（在循環內部重新定義）if np.random.rand() < self.delta:  # 挖掘模式(探索)r1, r2, r3 = np.random.rand(3)  # 只需要3個隨機數# 氣味導向的位置更新new_pos = self.best_position + F * r1 * alpha * I[i] * \np.abs(np.cos(2 * np.pi * r2) * (1 - np.cos(2 * np.pi * r3)))else:  # 采蜜模式(開發)r4 = np.random.rand()# 動態挖掘的位置更新new_pos = self.best_position + F * r4 * alpha * np.abs(self.best_position - self.positions[i])# 越界處理new_pos = np.clip(new_pos, self.lb, self.ub)# 計算新位置適應度new_fitness = self.obj_func(new_pos)# 更新個體位置和適應度if new_fitness < self.fitness[i]:self.positions[i] = new_posself.fitness[i] = new_fitnessif new_fitness < self.best_fitness:self.best_position = new_pos.copy()self.best_fitness = new_fitness# 記錄當前最優解self.convergence_curve[t] = self.best_fitnessprint(f"Iteration {t + 1}/{self.max_iter} - Best Fitness: {self.best_fitness:.6e}")return self.best_position, self.best_fitness, self.convergence_curve# ===================== CEC2017測試函數集實現 ========================
def cec17_f1(x):"""Shifted and Rotated Bent Cigar Function (Function 1)"""d = len(x)z = x - 100  # Shift to new optimumrotated_z = z @ rot_matrix[d] if d in rot_matrix else z  # Rotationreturn rotated_z[0] ** 2 + 1e6 * np.sum(rotated_z[1:] ** 2) + 100def cec17_f3(x):"""Shifted and Rotated Rosenbrock's Function (Function 3)"""d = len(x)z = 0.5 * (x - 40)  # Shift and scalerotated_z = z @ rot_matrix[d] if d in rot_matrix else z  # Rotationsum_val = 0for i in range(d - 1):sum_val += 100 * (rotated_z[i] ** 2 - rotated_z[i + 1]) ** 2 + (rotated_z[i] - 1) ** 2return sum_val + 300def cec17_f5(x):"""Shifted and Rotated Ackley's Function (Function 5)"""d = len(x)z = x + 50  # Shift to new optimumrotated_z = (z @ rot_matrix[d]) if d in rot_matrix else z  # Rotationsum1 = np.sum(rotated_z ** 2)sum2 = np.sum(np.cos(2 * np.pi * rotated_z))return -20 * np.exp(-0.2 * np.sqrt(sum1 / d)) - np.exp(sum2 / d) + 20 + np.e + 500# ===================== 輔助函數 ========================
def generate_rotation_matrix(dim):"""生成隨機旋轉矩陣（模擬CEC2017的旋轉特性）"""H = np.random.randn(dim, dim)Q, R = np.linalg.qr(H)return Qdef plot_results(convergence, func_name, dim):"""繪制收斂曲線"""plt.figure(figsize=(10, 6))plt.semilogy(convergence, 'b-', linewidth=2)plt.title(f'HBA on {func_name} - Dimension {dim}\nFinal Fitness: {convergence[-1]:.4e}', fontsize=12)plt.xlabel('Iteration', fontsize=12)plt.ylabel('Fitness (log scale)', fontsize=12)plt.grid(True, which='both', linestyle='--')plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))plt.tight_layout()plt.show()# ===================== 主程序 ========================
if __name__ == "__main__":# 實驗參數dim = 30  # 問題維度max_iter = 1000  # 最大迭代次數pop_size = 100  # 種群大小# 預生成旋轉矩陣 (模擬CEC2017特性)np.random.seed(42)rot_matrix = {dim: generate_rotation_matrix(dim)}# 定義要測試的函數test_functions = {"CEC2017 F1 (Bent Cigar)": cec17_f1,"CEC2017 F3 (Rosenbrock)": cec17_f3,"CEC2017 F5 (Ackley)": cec17_f5}# 對每個測試函數運行HBAresults = {}for name, func in test_functions.items():print(f"\n{'=' * 80}")print(f"Optimizing Function: {name}")print(f"{'=' * 80}")# 初始化并運行HBAoptimizer = HBA(obj_func=func,dim=dim,pop_size=pop_size,max_iter=max_iter,lb=-100,ub=100)best_sol, best_fit, convergence = optimizer.optimize()results[name] = convergence# 打印最終結果print(f"\n{'*' * 80}")print(f"{name} Optimization Result:")print(f"Best Solution: {best_sol[:5]}...")  # 只打印前5維print(f"Best Fitness: {best_fit:.6e}")print(f"Optimal Value Found: {func(best_sol):.6e}")print(f"{'*' * 80}\n")# 繪制收斂曲線plot_results(convergence, name, dim)# 比較所有函數的收斂曲線plt.figure(figsize=(12, 8))for name, conv in results.items():plt.semilogy(conv, label=name, linewidth=2)plt.title(f'HBA Performance on CEC2017 Functions (Dim={dim})', fontsize=16)plt.xlabel('Iteration', fontsize=12)plt.ylabel('Fitness (log scale)', fontsize=12)plt.legend(fontsize=12)plt.grid(True, which='both', linestyle='--')plt.tight_layout()plt.savefig('hba_cec2017_convergence.png', dpi=300)plt.show()

輸出值為：