0 引言

图的Steiner树问题(Steiner tree problem of graph，GSTP)是由最短路径和最小生成树^[1]两个著名的组合优化问题推广而来，广泛应用于超大规模集成电路设计^[2]、通信网络^[3]、组播路由^[4]等诸多工程领域。目前，已有多种求解GSTP的算法。例如，文献[5]为避免编码、解码和全局链路信息等操作，提出一种基于叶交叉机制的遗传算法。文献[6]提出并行2-近似Steiner最小树算法及其基于消息传递接口的分布式实现方法。尽管关于GSTP的求解算法已经取得了显著的研究进展，但现有算法仍有一定的局限性。

哈里斯鹰优化(Harris hawks optimization，HHO)算法具有可调参数少、鲁棒性强等优点^[7-8]，在解决工程问题时展现出较高的应用潜力和灵活性。文献[9]采用混合HHO算法解决移动云计算中的任务调度与分配问题，文献[10]利用多目标HHO算法解决雾计算中的服务部署问题，文献[11]通过并行多目标HHO算法来预测COVID-19患者的死亡风险。

HHO算法具有多样化的搜索策略及自适应特性，在求解GSTP时也展现出较强的适用性。但GSTP属于离散型组合优化问题，而HHO算法在迭代更新过程中生成的个体位置为连续值，这使得其在求解GSTP时需要进行离散化处理。

此外，HHO算法在求解复杂问题时仍面临求解精度低、探索与开发难以平衡、算法易“早熟”等挑战，目前尚未有研究者将该算法应用于GSTP的求解。本文提出一种改进的离散哈里斯鹰优化(improved binary Harris hawks optimization，IBHHO)算法。通过S型函数进行离散化处理，引入Logistic-Sine混合混沌映射增强种群分布的均匀性，采用动态自适应权重策略平衡算法的探索与开发进程，对当前最优个体进行自适应高斯—柯西混合变异扰动，避免算法陷入局部最优。在多个GSTP实例上设计实验，对比分析IBHHO算法、离散遗传算法(binary genetic algorithm，BGA)、离散粒子群优化(binary particle swarm optimization，BPSO)算法和离散哈里鹰优化(binary Harris hawks optimization，BHHO)算法的性能。结果表明，IBHHO算法在精确性和稳定性方面取得了显著增强，证明了所提方法的优越性和有效性。

1 理论知识 1.1 图的Steiner树问题

对于加权无向图G=(V′，E′)，A和S均为顶点集合V′的子集，且A∪S=V′，A∩S=$ \emptyset $。在边集E′上定义权值非负实数函数w: e→R⁺。GSTP旨在图G上寻找一棵跨越所有终端节点的树Tree=(V_T, E_T)，且A⊆V_T⊆V′, E_T⊆E′，使(1)式中的目标函数J最小化, 即

$ J( { Tree })=\sum\limits_{e \in E_{T}} w(e)。$

(1)

将A中的节点称为终端节点，S中的节点称为Steiner点。

1.2 相关性质

相关性质如下。

1) Tree不包含度为0或1的Steiner点。

2) 与度为1的终端节点相邻的点在Tree中。

3) 对于度为2的Steiner节点v，若其相邻的两个节点v₁, v₂间存在边，且w(v, v₁)+w(v, v₂)>w(v₁, v₂)，则Tree不包含v。

4) 对于度为2的终端节点v，若其相邻的两个节点v₁, v₂均为终端节点，则min{w(v, v₁), w(v, v₂)}相应的边在Tree中。

5) 对于度为2的终端节点v₀，若其相邻的两个节点v₁, v₂均为Steiner点，且v₁, v₂间存在边e₁₂，w₀₁+w₁₂<w₀₂，则e₀₁与v₁在Tree中。

2 HHO算法 2.1 全局探索阶段

当q < 0.5时，哈里斯鹰发现猎物；当q≥0.5时，未发现猎物。位置更新公式为

$ X(t+1)=\left\{\begin{array}{l} X_{\text {rand }}(t)-r_{1} \mid X_{\text {rand }}(t)- \\ \quad 2 r_{2} X(t) \mid, q \geqslant 0.5, \\ \left(X_{\text {rabbit }}(t)-X_{m}(t)\right)- \\ \quad r_{3}\left(l b+r_{4}(u b-l b)\right), q <0.5, \end{array}\right.$

(2)

式中：X(t)为第t次迭代中鹰的位置；X_rabbit(t)为当前最优个体；ub和lb分别为搜索范围的上界和下界；X_rand(t)为随机选择的个体；r₁, r₂, r₃, r₄和q均为(0, 1)内的随机数；X_m为当前种群的平均位置，可表示为

$ X_{m}(t)=\frac{1}{N} \sum\limits_{i=1}^{N} X_{i}(t), $

(3)

其中：N表示哈里斯鹰种群规模。

2.2 探索与开发转换阶段

当逃逸能量|E|≥1.0时，进入探索阶段；当|E| < 1.0时，进入开发阶段。E的表达式为

$ \begin{gather*} E=2 E_{0}\left(1-\frac{t}{T}\right), \\ E_{0}=2 r-1, \end{gather*}$

(4)

其中：T为最大迭代次数；E₀为初始逃逸能量；逃逸概率r为(0, 1)内的随机数。设N=50，T=500，E的变化如图 1所示。

2.3 开发阶段

当r < 0.5时，猎物逃脱成功；当r≥0.5时，猎物逃脱失败。

以下根据|E|和r来确定狩猎策略。

1) 软包围：当r≥0.5且0.5≤|E|<1.0时，猎物能量充足但无法逃脱，鹰的位置更新公式为

$ \begin{gather*} X(t+1)=\Delta X(t)-E\left|J X_{\mathrm{rabbit}}(t)-X(t)\right|, \\ J=2\left(1-r_{5}\right) , \end{gather*}$

(5)

$ \Delta X(t)=X_{\mathrm{rabbit}}(t)-X(t),$

(6)

其中：ΔX(t)为猎物在第t次迭代中的位置与当前位置之差；J为随机跳跃距离；r₅为(0, 1)内的随机数。

2) 硬包围：当r≥0.5且|E| < 0.5时，猎物体力不足无法逃脱，鹰的位置更新公式为

$ X(t+1)=X_{\text {rabbit }}(t)-E|\Delta X(t)| 。$

(7)

3) 渐进式快速俯冲软包围：当r < 0.5且0.5≤|E| < 1.0时，猎物体力充沛有机会逃脱，鹰的位置更新策略为

$ Y=X_{\text {rabbit }}(t)-E\left|J X_{\text {rabbit }}(t)-X(t)\right| 。$

(8)

若上述策略无效时，引入Levy飞行函数模拟猎物的逃脱模式，此时鹰的更新策略为

$ Z=Y+\boldsymbol{C} \times \operatorname{Levy}(D), $

(9)

式中：D为问题维度，C的大小为1×D，Levy飞行函数的计算公式为

$ \begin{gather*} \operatorname{Levy}(D)=0.01 \times \frac{u \times \sigma}{|v|^{\frac{1}{\beta}}}, \\ \sigma=\left(\frac{\Gamma(1+\beta) \times \sin \left(\frac{{\rm{ \mathsf{ π}}} \beta}{2}\right)}{\Gamma\left(\frac{1+\beta}{2}\right) \times \beta \times 2^{\frac{\beta-1}{2}}}\right)^{\frac{1}{\beta}}, \end{gather*}$

(10)

其中：u, v为(0, 1)内的随机数；β=1.5；Γ为伽马函数。在该阶段鹰的位置更新公式为

$ X(t+1)= \begin{cases}Y, & \text { if } F(Y)<F(X(t)), \\ Z, & \text { if } F(Z)<F(X(t)), \end{cases}$

(11)

其中：F(·)为适应度函数。

4) 渐进式快速俯冲硬包围：当r < 0.5且|E| < 0.5时，猎物体力不足，鹰种群先缩短最优位置与平均位置间的距离再进行围攻，其更新策略为

$ Y=X_{\text {rabbit }}(t)-E\left|J X_{\text {rabbit }}(t)-X_{m}(t)\right|, $

(12)

$ X(t+1)= \begin{cases}Y, & \text { if } F(Y)<F(X(t)), \\ Z, & \text { if } F(Z)<F(X(t))。\end{cases} $

(13)

3 改进的HHO算法(IBHHO算法) 3.1 离散化HHO算法

S型函数可以平滑地将连续空间中的值映射为二进制值，从而使得算法能够在离散的解空间内搜索^[12-13]。

在求解GSTP的过程中，哈里斯鹰的位置X_i为连续值，其中X_i^j表示图G中第j个节点的选择状态，S型函数可表示为

$ S_{4}\left(X_{i}^{j}\right)=\frac{1}{1+\mathrm{e}^{-10\left(X_{i}^{j}-0.5\right)}}。$

(14)

通过S型函数，将连续值映射到[0, 1]内的概率值，通过该值判定取1或0。根据S₄(X_i^j)生成二进制值B，其表达式为

$ B= \begin{cases}1, & \text { if rand }<S_{4}\left(X_{i}^{j}\right), \\ 0, & \text { otherwise }, \end{cases}$

(15)

其中：rand为[0, 1]内的随机数。若S₄(X_i^j)>rand，则选择节点j，即将其对应的位置设为1，否则设为0，而终端节点的位置均为1。

3.2 Logistic-Sine混合混沌映射

智能优化算法的初始种群分布不均匀，算法易陷入局部最优解^[14]。为改善种群多样性并增强算法搜索性能，通常引入Logistic或Sine混沌映射^[15]。

Logistic映射计算复杂度较低，适用于大规模计算任务，其表达式为

$ x_{i+1}=\mu x_{i}\left(1-x_{i}\right), $

(16)

其中：x为随机序列；μ为混沌乘数，取值为0.5。

Sine映射有助于增强优化算法的全局搜索能力，其表达式为

$ x_{i+1}=\frac{\mu \sin \left({\rm{ \mathsf{ π}}} x_{i}\right)}{4} 。$

(17)

为开发更优的混沌系统，在种群初始化阶段引入Logistic-Sine混合混沌映射，以增强种群分布的均匀性，表达式为

$ X_{i+1}=\mu X_{i}\left(1-X_{i}\right)+\frac{(4-\mu) \sin \left({\rm{ \mathsf{ π}}} X_{i}\right)}{4} 。$

(18)

3.3 动态自适应逃逸能量

由图 1可知，当迭代次数大于T/2时，E的值恒小于1.0，算法只进行开发操作。此外，E的变化呈线性递减，使得算法难以平衡探索与开发行为。因此，引入动态自适应权重策略，以提高猎物逃逸能量的非线性表达，从而增强算法的寻优能力与搜索精度。通过非线性递减模式设置E₁的值，用余弦函数来描述其值的变化。若E≥E₁，算法进入全局探索阶段，否则进入局部开发阶段。自适应逃逸能量的更新公式为

$E=E_{0} E_{1}, $

(19)

$ \begin{array}{l} E_{1}=E_{0}\left(\frac{w_{\max }-w_{\min }}{2} \cos \left(\frac{{\rm{ \mathsf{ π}}} t}{T}\right)+\right. \\ \left.\frac{w_{\max }+w_{\min }}{2}-w_{\min }\right), \end{array}$

(20)

其中：w_max和w_min的取值分别为2.1和0.1。

3.4 自适应高斯—柯西混合变异扰动

在HHO算法迭代后期，个体快速同化，种群搜索停滞，导致算法陷入局部最优解，为此通常引入柯西变异和高斯变异等操作算子。柯西变异虽然具有较大的搜索范围，但过大的步长可能导致个体偏离最优值。高斯变异以较大的概率生成较小的变异值，使得其只在较小的范围内具有良好的搜索能力。因此，在求解GSTP的过程中，根据上述两种变异操作算子的优缺点，提出了自适应高斯—柯西混合变异扰动策略。为降低算法的时间复杂度，该操作仅对当前最优个体进行扰动，其计算公式为

$ \begin{align*} & H(t)=X_{\text {rabbit }}(t)\left(1+\mu_{1} \operatorname{Gauss}(\sigma)+\right. \\ & \left.\mu_{2} \operatorname{Cauchy}(\sigma)\right), \end{align*}$

(21)

$ \mu_{1}=\frac{t}{T}, \mu_{2}=1-\frac{t}{T}, $

(22)

其中：H(t)为扰动后的位置；μ₁, μ₂为权重系数；Gauss(σ)和Cauchy(σ)分别为高斯变异算子和柯西变异算子。

3.5 IBHHO算法流程

IBHHO算法流程如图 2所示。

4 实验设计 4.1 参数设置

选用BHHO、BGA和BPSO^[16]算法与本文IBHHO算法进行对比实验，设四种算法中N=50，T=500。BGA算法的突变概率为0.2；BPSO算法的自我学习因子和全局学习因子均为1.5，最大速度和最小速度分别为1和-1；IBHHO算法的ub和lb分别为1和0。

4.2 算法性能分析

选用OR-Library^[17]中关于GSTP的标准测试实例，每个实例独立运行20次。算法性能评估指标采用最优值(Best)和误差率(RPD)，其中误差率计算公式为

$ R P D=\frac{B e s t-O P T}{O P T} \times 100 \%, $

(23)

其中：OPT为已知最优解。

实验1  根据边集数选取OR-Library中的部分i080实例，四种算法的实验结果对比见表 1。

表 1中|V′|为节点总数，|E′|为边集总数，|A|为终端节点总数，OPT为已知最优解，Best为求得最优Steiner树时的最小权重。

由表 1可知，BGA、BPSO和BHHO算法分别在0、5、13个i080实例上达到已知最优解，其RPD较大。IBHHO算法在所有实例上均能达到已知最优解，RPD均为0，表明IBHHO算法在解决GSTP时具有显著优势。

图 3为四种算法在不同终端节点数实例(i080-031、i080-144、i080-243、i080-334)上的适应度曲线。可以看出，相比BGA、BPSO和BHHO算法，IBHHO算法在大多数实例上收敛速度更快，最优解精度较高且更加稳定。

实验2  根据边集数选取OR-Library中的部分i160实例，四种算法的实验结果对比见表 2。

从表 2可以看出，BGA、BPSO和BHHO算法达到已知最优解的占比分别为0、10%和15%，其RPD相对较大。IBHHO算法达到已知最优解的占比为70%，且在绝大多数实例上RPD为0，表明IBHHO算法性能较好。

图 4为四种算法在不同终端节点数实例(i160-141和i160-315)上的适应度曲线。可以看出，IBHHO算法的求解精度更高，并且IBHHO算法曲线在取得最优值时的最小迭代次数更小，说明该算法的收敛速度更快，证明了改进策略对算法寻优能力的提升。

5 结语

本文围绕GSTP展开研究，提出一种融合多策略的HHO算法。通过S型函数离散化解的搜索空间，利用Logistic-Sine混合混沌映射优化种群初始化过程，引入动态自适应权重策略增强猎物逃逸能量的非线性表达，对当前最优个体采用自适应高斯—柯西混合变异扰动，避免算法陷入局部最优。最后，在OR-Library标准测试集上设计对比实验，通过分析实验结果可知，本文提出的IBHHO算法在求解GSTP时具有较高的求解精度和较快的收敛速度。在未来研究中，可以将该算法应用于其他组合优化问题，使算法具有更强的适用性。