如图1所示，假设由N个无人机组成的无人机蜂群分布在三维欧式空间内。$ {\boldsymbol{x}}_{i}=\left[{a}_{i},{b}_{i},{c}_{i}\right] $表示第i个无人机在初始拓扑中的位置(图1中五边形)，$ i=\mathrm{1,2},\cdots ,N $；$\boldsymbol{X} = [{\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2},\cdots ,{\boldsymbol{x}}_{N}{]}^{{\rm{T}}} \in {\mathbb{R}}^{N\times 3}$表示初始拓扑坐标矩阵，${\boldsymbol{z}}_{j} = \left[{a}_{j},{b}_{j},{c}_{j}\right]$表示相同坐标系下目标拓扑中的第j个坐标位置(图1中圆点)，$j=\mathrm{1,2},\cdots ,N $；$ \boldsymbol{Z}=[{\boldsymbol{z}}_{1},{\boldsymbol{z}}_{2},\cdots ,{\boldsymbol{z}}_{N}{]}^{{\rm{T}}}\in {\mathbb{R}}^{N\times 3} $表示目标拓扑坐标矩阵，目标拓扑形状如图1中实线表示。令$ \boldsymbol{M}\left(k\right)=\left[m\right(k{)}_{ij}{]}_{N\times N} $表示从初始拓扑到目标拓扑的映射关系(图1中带箭头虚线)，当初始拓扑中的节点i被指派给目标位置j时，$ m(k{)}_{ij}=1 $，否则$ m(k{)}_{ij}=0 $。假设$ {\boldsymbol{M}}^{{\rm{o}}}\left(k\right)=\left[{m}^{{\rm{o}}}\right(k{)}_{ij}{]}_{N\times N} $表示当前$ {k} $下全局能耗最小的最优映射关系。

图  1  无人机蜂群拓扑构型示意图

假设$ {\boldsymbol{O}}_{x} $和$ {\boldsymbol{O}}_{z} $分别代表初始拓扑和目标拓扑的几何中心(图1中2个正方形)，$ \boldsymbol{k}\in {\mathbb{R}}^{1\times 3} $表示$ \boldsymbol{X} $和$ \boldsymbol{Z} $的相对位置向量，那么：

最优拓扑构型位置$ {{\boldsymbol{O}}^{\rm{o}}}_{z} $可以表示为$ {{\boldsymbol{O}}^{{{\rm{o}}}}}_{z}={\boldsymbol{k}}^{\rm{o}}+{\boldsymbol{O}}_{x} $，其中$ {\boldsymbol{k}}^{\rm{o}} $是待求的相对位置向量。

拓扑构型全局能耗$ {C}_{g} $可以表示为蜂群中所有无人机从初始拓扑位置到目标拓扑位置的欧氏距离之和^[6-8]。假设$ {c}_{ij} $表示从初始拓扑中第i个节点到目标拓扑第j个位置的欧式距离，那么：

假设将目标拓扑平移向量$ {\boldsymbol{k}}'=\left[{a}_{k},{b}_{k},{c}_{k}\right]\in {\mathbb{R}}^{1\times 3} $，那么$ {\boldsymbol{Z}}_{k}=\boldsymbol{Z}+{\boldsymbol{k}}' $，所有的目标拓扑位置$ {\boldsymbol{z}}_{j} $会被平移至$ {\boldsymbol{z}}_{j}+{\boldsymbol{k}}' $。此时平移后第i个节点和第j个目标位置对应的欧氏距离可以表示为：

从式(3)中可以看出，$ {c}_{ij} $与目标拓扑平移向量$ {\boldsymbol{k}}' $相关，而目标拓扑平移$ {\boldsymbol{k}}' $后$ \boldsymbol{X} $和$ \boldsymbol{Z} $新的相对位置向量$ \boldsymbol{k} $会变为$ \boldsymbol{k}+{\boldsymbol{k}}' $，即在原来的相对位置向量上做了修正，因此可认为$ {c}_{ij} $是与$ \boldsymbol{k} $相关的函数，表示为$ c(k{)}_{ij} $，并通过迭代$ \boldsymbol{k} $来求得最优拓扑构型位置。令$ \boldsymbol{C} $表示从初始拓扑到目标拓扑的代价矩阵，即$ \boldsymbol{C}=\left[c\right(k{)}_{ij}{]}_{N\times N} $。根据线性分配算法可知，从初始拓扑到目标拓扑的映射关系和通过线性分配算法求得的全局能耗会随$ c(k{)}_{ij} $变化。因此，全局能耗$ {C}_{g} $也和$ \boldsymbol{k} $相关，表示为：

在拓扑构型场景下，初始拓扑$ \boldsymbol{X} $和目标拓扑初始坐标$ \boldsymbol{Z} $一旦确定，那么就可以得到$ \boldsymbol{k} $和对应的$ c(k{)}_{ij} $，此时最佳映射关系$ {\boldsymbol{M}}^{{\rm{o}}}\left(k\right) $及全局能耗$ C(k{)}_{g} $可以通过线性分配算法求得，从式(4)可以看出，相对位置向量$ \boldsymbol{k} $(拓扑构型位置)和最佳映射关系相互影响，并且都会影响全局能耗。

本文的目标是在最小化全局能耗的前提下求得最优拓扑构型映射关系$ {\boldsymbol{M}}^{{\rm{o}}}\left(k\right) $和最优拓扑构型位置$ {{\boldsymbol{O}}^{\rm{o}}}_{z} $，使式(4)中$ C(k{)}_{g} $在有限的三维空间中最小，即：

基于给出的联合优化模型，本文调用群体智能算法进行求解。

具有简单、灵活、无微分机制和避免局部最优^[10]等优点的群体智能算法特别适合该联合优化求解模型。集群智能算法^[10-13]的主要思想是模仿自然界中动物或者物体在群体中的行为，这类算法的主要特点是使用随机方式来寻找最优解，并在求最优解的过程中共享所有个体的最优解信息。另外这类算法的求解过程与被优化的问题相互独立，适合寻找特定优化问题最优解，而不需要关注问题搜索空间的非线性类型及其约束。

通过调用群体智能算法求解这个联合优化模型，分为初始化和群体智能算法迭代两个阶段，如图2所示。为了公平比较，根据文献[6-8]，确定相对位置向量$ \boldsymbol{k} $后，最佳映射关系$ {\boldsymbol{M}}^{{\rm{o}}}\left(k\right) $及对应的全局最小能耗$ C(k{)}_{g} $可通过线性分配算法求得，通过线性分配算法求得匹配关系和全局最小能耗的具体过程不做详细讨论。

1) 初始化阶段

输入拓扑构型及群体智能算法相关参数，初始化算法的最大迭代次数Max_iter、搜寻节点数N、搜寻区间上限ub、搜寻区间下限lb、初始拓扑坐标矩阵$ \boldsymbol{X} $、目标拓扑坐标矩阵$ \boldsymbol{Z} $，初始化群体智能算法相关参数。

2) 群体智能算法迭代阶段

根据拓扑构型参数及群体智能算法初始化参数或所选群体智能算法位置，更新每个搜索点的相对位置向量$ \boldsymbol{k} $。根据式(3)计算当前$ \boldsymbol{k} $下初始拓扑到目标拓扑的代价矩阵$ \boldsymbol{C} $。根据代价矩阵$ \boldsymbol{C} $及线性分配算法求得当前$ \boldsymbol{k} $对应的映射关系和全局能耗。记录本次迭代中最小全局能耗及对应的映射关系。

重复以上步骤直到满足结束条件。

图  2  群体智能算法拓扑构型流程图

本文提出的优化模型求解框架中，所选群体智能算法中各个搜寻点的位置对应于相对位置向量k。另外，图2中的群体智能算法并不局限于文献[10-13]中的4种算法，现有的群体智能算法均可适配于本框架。所有群体智能算法需要的通用参数包括算法最大迭代次数、搜寻节点数、搜寻区间上限和下限以及用于计算匹配关系和最小全局能耗的初始拓扑位置和最初的目标拓扑位置。最初的目标拓扑位置是指根据初始拓扑坐标矩阵和目标拓扑坐标矩阵选取一个初始的相对位置向量$ \boldsymbol{k} $。

本文选取灰狼优化算法(GWO)(伪代码见算法1)、均衡优化算法(EO)(伪代码见算法2)、穷人富人优化算法(PRO)(伪代码见算法3)来进一步说明基于所提求解框架的具体求解步骤。

算法1　灰狼优化算法(GWO)

输入：目标函数$ C(k{)}_{g} $，最大迭代次数Max_iter

输出：最小能耗及对应的位置$ {X}_{\alpha } $

初始化灰狼狼群大小及初始位置$ {\boldsymbol{X}}_{i}\left(i=\mathrm{1,2},\cdots , N\right) $；

初始化 $ \boldsymbol{a},\boldsymbol{A},\boldsymbol{C} $；

根据$ {\boldsymbol{X}}_{i} $计算目标函数值$ C({\boldsymbol{X}}_{i}{)}_{g} $；

$ {\boldsymbol{X}_{\alpha }}=\alpha $狼(最佳)能耗值对应位置；

$ {\boldsymbol{X}_{\beta }}=\beta $狼(次佳)能耗值对应位置；

$ {\boldsymbol{X}_{\gamma }}=\gamma $狼(次次佳)能耗值对应位置；

while iter < Max_iter do

　 for i=1:N do

　　 通过$ \boldsymbol{X}\left(t+1\right)=\dfrac{{\boldsymbol{X}}_{\alpha }+{\boldsymbol{X}}_{\beta }+{\boldsymbol{X}_{\gamma }}}{3} $更新当前狼位置；

　end

　更新向量$ \boldsymbol{a},\boldsymbol{A},\boldsymbol{C} $；

　计算所有狼所在位置的目标函数值$ C({\boldsymbol{X}}_{i}{)}_{g} $；

　更新$ {\boldsymbol{X}_{\alpha }} $，$ {\boldsymbol{X}_{\beta }} $，$ {\boldsymbol{X}_{\gamma }} $；

　$ {\rm{iter}}={\rm{iter}}+1 $；

end

return $ {\boldsymbol{X}_{\alpha }} $，$ C({\boldsymbol{X}_{\alpha }}{)}_{g} $

算法2　均衡优化算法(EO)

输入：目标函数$ C{\left(k\right)}_{g} $，最大迭代次数Max_iter

输出：最小能耗 $ C({\boldsymbol{C}}_{i}{)}_{g} $及对应位置$ {\boldsymbol{C}}_{i} $

初始化粒子数目N，参数$ {a}_{1}=2 $，$ {a}_{2}=1 $，${\rm{GP}}= 0.5$; 初始化粒子位置$ {C}_{i}^{{\rm{initial}}}={C}_{{\rm{min}}}+{\rm{ran}}{{\rm{d}}}_{i}\left({C}_{{\rm{max}}}-{C}_{{\rm{min}}}\right)， i=\mathrm{1,2},\cdots ,N $；

while Iter < Max_iter do

　 for i=1:N do

　　 计算第i个粒子的目标函数$ C({\boldsymbol{C}}_{i}{)}_{g} $的值；

　　 更新4个最小$ C({\boldsymbol{C}}_{i}{)}_{g} $值及对应的位置$ {\boldsymbol{C}}_{i} $；

　 end

更新平均状态$ {\boldsymbol{C}}_{{\rm{ave}}} = ({\boldsymbol{C}}_{{\rm{eq}}1} + {\boldsymbol{C}}_{{\rm{eq}}2} + {\boldsymbol{C}}_{{\rm{eq}}3} + {\boldsymbol{C}}_{{\rm{eq}}4})/ 4 $，建立位置更新向量可选池${\boldsymbol{C}}_{{\rm{eq,pool}}}=\left\{{\boldsymbol{C}}_{\mathrm{eq}\left(1\right)}, {\boldsymbol{C}}_{\mathrm{eq}\left(2\right)}, {\boldsymbol{C}}_{{\rm{eq}}\left(3\right)}, {\boldsymbol{C}}_{\mathrm{eq}\left(4\right)},{\boldsymbol{C}}_{{\rm{eq}}\left(\mathrm{a}\mathrm{v}\mathrm{e}\right)}\right\}$；

令$ t={(1-\dfrac{{\rm{Iter}}}{{\rm{Max}}\text{\_}{\rm{iter}}})}^{\left({a}_{2}\tfrac{{\rm{Iter}}}{{\rm{Max}}\text{\_}{\rm{iter}}}\right)} $；

　 for i = 1:N do

　　 从均衡状态池随机选1个候选向量${\boldsymbol C _{eq}}$；

　　 生成随机向量$ \boldsymbol{\lambda } $，$ \boldsymbol{r} $；

　　 构建$ \boldsymbol{F}={a}_{1}{\rm{sign}}(\boldsymbol{r}-0.5)[{{\rm{e}}}^{-\boldsymbol{\lambda }t}-1] $；

　　 构建$\bf{GCP} = \left\{ {\begin{array}{*{20}{c}}{0.5{r_1}}&{{r_2} \ge {\rm{GP}}}\\0&{{r_2} < {\rm{GP}}}\end{array}} \right.$；

　　 构建$ {\boldsymbol G _0} = {\bf {GCP}} ({\boldsymbol C _{{\rm{eq}}}} - {\boldsymbol{\lambda}} {\boldsymbol{C}}) $；

　　 构建$ \boldsymbol{G}={\boldsymbol{G}}_{0}\boldsymbol{F} $；

　　 更新$\boldsymbol{C}={\boldsymbol{C}}_{{\rm{eq}}}+(\boldsymbol{C}-{\boldsymbol{C}}_{{\rm{eq}}})\boldsymbol{F}+ \dfrac{\boldsymbol{G}}{\boldsymbol{\lambda }V}\left(1-\boldsymbol{F}\right)$；

　 end

　 $ {\rm{iter}}={\rm{iter}}+1 $；

end

return $ C({\boldsymbol{C}}_{i}{)}_{g} $，$ {\boldsymbol{C}}_{i} $

算法3　穷人富人优化算法(PRO)

输入：目标函数$ C(k{)}_{g} $，最大迭代次数Max_iter，变异概率Pmut

输出：最小能耗$ C(k{)}_{g} $及对应的位置$ {X}_{\alpha } $

初始化穷人富人群体大小及初始位置$ {\bf{XR}}_{i}\left(i= \mathrm{1,2},\cdots ,N\right) $，$ {\bf{XP}}_{i}\left(i=\mathrm{1,2},\cdots ,N\right) $；

对每一个穷人和富人计算目标函数值$ C({\bf{XP}}_{i}{)}_{g} $，$ C({\bf{XR}}_{i}{)}_{g} $，根据目标函数值大小区分穷人群体和富人群体；

while iter < Max_iter do

根据目标函数值区分穷人和富人群体，选择富人最小值$ {\boldsymbol {X}} _{{\rm{rich}},{\rm{best}}}^{{\rm{old}}} $、平均值$ {\boldsymbol{X}} _{{\rm{rich}},{\rm{mean}}}^{{\rm{old}}} $、最大值$ {\boldsymbol{X}} _{{\rm{rich}},{\rm{worst}}}^{{\rm{old}}} $；

for i=1:N do %富人群体

更新位置$ {\boldsymbol{X}} _{{\rm{rich}},{{i}}}^{{\rm{new}}} = {\boldsymbol{X}} _{{\rm{rich}},{{i}}}^{{\rm{old}}} + {{r}}\left[ { {\boldsymbol{X}} _{{\rm{rich}},{\rm{i}}}^{{\rm{old}}} - {\boldsymbol{X}} _{{\rm{poor}},{\rm{best}}}^{{\rm{old}}}} \right]$；

if $ {\rm{rand}} < {P}_{{\rm{mut}}} $ then

　 　 $ {\boldsymbol{X}} _{{\rm{rich}},{{i}}}^{{\rm{new}}} = {\boldsymbol{X}} _{{\rm{rich}},{{i}}}^{{\rm{new}}} + {\rm{rand}}n$；

end

计算新位置的目标函数值；

end

for i=1:N do %穷人群体

计算模式值${\rm{Pattern}} = \dfrac{{ {\boldsymbol{X}} _{{\rm{rich}},{\rm{best}}}^{{\rm{old}}} + {\boldsymbol{X}} _{{\rm{rich}},{\rm{mean}}}^{{\rm{old}}} {\boldsymbol{X}} _{{\rm{rich}},{\rm{worst}}}^{{\rm{old}}}}}{3}$；

更新位置$ {\boldsymbol{X}} _{{\rm{poor}},{{i}}}^{{\rm{new}}} = {\boldsymbol{X}} _{{\rm{poor}},{{i}}}^{{\rm{old}}} + [{{r}}\left( {{\rm{Pattern}}} \right) - {\boldsymbol{X}} _{{\rm{poor}},{{i}}}^{{\rm{old}}}]$；

if $ {\rm{rand}} < {P}_{{\rm{mut}}} $ then

　 $ {\boldsymbol{X}} _{{\rm{poor}},{{i}}}^{{\rm{new}}} = {\boldsymbol{X}} _{{\rm{poor}},{{i}}}^{{\rm{new}}} + {\rm{rand}}n$；

end

　 计算新位置的目标函数值；

　 end

选出最小$ C({\bf{XP}}_{i}{)}_{g} $及对应位置$ {\bf{XP}}_{i} $；

$ {\rm{iter}}={\rm{iter}}+1 $；

end

return $ {\bf{XP}}_{i} $，$ C({\bf{XP}}_{i}{)}_{g} $

在群体智能算法中，需要初始化各个搜寻点的初始位置。如在PSO算法初始化过程中，需要先随机生成所有粒子的初始位置，一般是通过随机函数生成随机的位置。

值得注意的是，初始化过程会影响群体智能算法的收敛速度。在本文的拓扑构型场景下，各个群体智能算法的初始位置可被进一步优化，加速算法收敛速度。从直观的角度来看，目标拓扑和初始拓扑节点重叠的越多，全局能耗就越小。因此在初始化群体智能算法各搜寻点的初始位置时，可以初始拓扑坐标作为各个搜寻点的初始位置，而不是按照原始算法中的随机生成初始位置。另外，在本场景中，目标拓扑位置与初始拓扑位置越接近全局能耗越小，因此可在初始化阶段假设初始拓扑和目标拓扑坐标原点重合，即取$ {k}={{0}} $，并取其中的最小值作为所选群体智能算法的初始最小值。

本文通过MATLAB仿真验证以全局能耗最小化为目标的无人机蜂群拓扑构建联合优化模型，以及本文提出的可以灵活调度不同群体智能算法实现模型优化求解的框架和群体智能算法加速收敛的策略。仿真所运行的MATLAB版本为R2021b 64位，电脑配置为12 th Gen Intel(R) Core(TM) i9 CPU，128 GB RAM，Windows 11 Pro。

在仿真过程中，设定算法最大迭代次数Max_iter=100，每次迭代中搜寻节点数N=32，搜寻区间上限ub=800，搜寻区间下限lb=−800，初始拓扑坐标矩阵$ \boldsymbol{X} $由随机分布在(−400, 400)的三维空间中的随机点组成，目标拓扑坐标矩阵$ \boldsymbol{Z} $为边长为400的立方体，初始拓扑和目标拓扑初始位置的原点重合。PSO算法初始设置为c1=c2=2，w=0.9；EO算法初始值设置为a1=2，a2=1，GP=0.5，V=1；PRO算法初始值设置为Pmut=0.06。

本节选取PSO、EO、PRO、GWO这4种群体智能算法求解联合优化模型。仿真结果如图3所示，其中星号代表节点的初始位置。

图  3  立方体形状目标拓扑构建仿真结果

图3中，实线勾勒出的是GWO算法求得的目标拓扑构型位置上的拓扑外形。虚线表示初始拓扑到GWO算法得到的拓扑构型位置的映射关系。仿真结果表明，通过所提方法可以实现期望的最佳拓扑构型以及最优拓扑构型位置，PRO和GWO最终收敛结果一致，二者得到的全局最小能耗低于PSO算法和EO算法相应值。

为了验证初始位置对收敛速度的影响，根据前面分析，本文进行了2组100次独立仿真，2组仿真分别设置了$ \boldsymbol{k}\ne $0(见图4)和$ \boldsymbol{k}= $0(见图5)以观察各个群体智能算法的收敛速度。仿真结果表明，在$ \boldsymbol{k}\ne $0时，GWO算法、PRO算法和PSO算法平均在35次内收敛，而EO算法平均需要75次才能收敛。在$ \boldsymbol{k}= $0条件下，EO算法、GWO算法和PRO算法均在可以在8次迭代内收敛，而PSO算法依然需要30余次迭代才能收敛。

图  4  $ \boldsymbol{k}\ne {\bf{0}} $时各群体算法收敛速度

图  5  $ \boldsymbol{k}= $0时各群体算法收敛速度

为进一步研究各加速收敛策略的有效性，将GWO算法中各狼的初始位置设置为随机生成(设置1)，各狼的初始位置选取初始拓扑位置且目标拓扑与初始拓扑原点重合，即$ \boldsymbol{k}= $0(设置2)，并进行100次独立仿真以对比各狼初始位置对收敛速度的影响，其余仿真参数如前所述，仿真结果如图6所示。仿真结果表明，目标拓扑的初始位置选取对收敛影响较大，当目标拓扑原点与初始拓扑原点不重合时，算法经过25次迭代后才能收敛，$ \boldsymbol{k}= $0算法可以在8次迭代内收敛，所提的加速收敛策略提高了群体智能算法收敛速度。

图  6  GWO算法不同初值及位置更新收敛对比

为了对比上述4种群体智能算法的性能，本文进一步分析了各算法100次独立仿真的计算时间。图7为上述4种群体智能优化算法100次独立仿真的计算时间箱式图统计。仿真结果表明，EO算法和GWO算法计算时间相对较短，平均计算时间分别是1.6 s/100次迭代和1.7 s/100次迭代，PRO平均计算时间为3.2 s/100次迭代。

图  7  群体智能算法计算时间对比

无人机蜂群已被广泛应用于军用和民用场景，在应用中需要通过拓扑构型来形成特定的拓扑形状以完成指定任务。拓扑构型过程中需要确定最佳映射关系和最优拓扑构型位置。这两个方面相互影响并直接关系到无人机蜂群的全局能量消耗。本文对映射关系和拓扑构型位置进行联合建模，给出了以全局能耗最小化为目标的无人机蜂群拓扑构建联合优化模型，并提出了可以灵活调用不同群体智能算法实现模型优化的求解框架。此外还讨论了群体智能算法加速收敛的策略。