多无人机辅助通信的节能卸载策略

Energy⁃efficient offloading strategy research for multi⁃UAV assisted communication

摘要

针对偏远地区或灾害场景中基础设施不完善、无法提供可靠通信的问题，提出了一种多无人机辅助用户的非正交多址接入中继系统。该系统综合考虑了信息因果、发射功率以及回传时延等约束，联合优化了卸载策略、发射功率、资源分配和轨迹以实现最小化系统能耗这一目标。针对优化问题的非凸性和复杂性，设计了一种两阶段在线资源分配方案进行求解。第一阶段，为消除复杂场景下对不确定信息的依赖，运用李雅普诺夫优化理论将优化问题分解为3个子问题。第二阶段，提出了一种交替迭代优化算法：首先，基于最大权匹配算法获得了用户卸载决策；其次，利用辅助变量法、连续凸逼近和拉格朗日对偶得到了传输功率和频率分配的闭式解；最后，将轨迹规划问题转化为凸问题，并使用凸优化工具进行求解。仿真结果表明：相较于基准算法，本文方案在满足系统长期稳定性约束的前提下，显著降低了系统能耗。

Abstract

To address the problem of inadequate infrastructure and inability to provide reliable communication in remote areas or disaster scenarios， a non orthogonal multiple access relay system with multiple drones assisting users is proposed. The system comprehensively considers constraints such as information causality， transmission power， and return delay， and jointly optimizes the offloading strategy， transmission power， resource allocation， and trajectory to achieve the goal of minimizing system energy consumption. Due to the non convexity and complexity of optimization problems， a two-stage online resource allocation scheme is designed for solving. In the first stage， to eliminate the dependence on uncertain information in complex scenarios， Lyapunov optimization theory is applied to decompose the optimization problem into three sub problems. In the second stage， an alternating iterative optimization algorithm was proposed： firstly， the user offloading decision was obtained based on the maximum weight matching algorithm； secondly， closed form solutions for transmission power and frequency allocation were obtained by utilizing auxiliary variable method， continuous convex approximation， and Lagrangian duality； finally， the trajectory planning problem is transformed into a convex problem and solved using the convex optimization tool. The simulation results show that compared to the benchmark algorithm， the proposed scheme significantly reduces system energy consumption while satisfying the long-term stability constraints of the system.

Graphical abstract

关键词

计算机应用 / 无人机 / 航迹优化 / 资源分配 / 能耗最小化

Key words

computer application / unmanned aerial vehicle / trajectory planning / resource allocation / energy minimization

0 引 言

在偏远山区和受灾地区等地理区域，传统服务模式已经无法满足各种新型互联网应用的需求^［1］。无人机（Unmanned aerial vehicle，UAV）辅助移动边缘计算（Mobile edge computing，MEC）作为一种具备按需操作、灵活部署以及优质链路特性的解决方案，引起了广泛关注^［2］。这一架构通过将任务上传到无人机进行处理，或通过中继传输到基站服务器进行处理，以实现更为即时和高效的数据处理^［3，4］。

现有工作中，文献［5］和文献［6］研究了单无人机充当中继，服务多地面用户的情景，极大提升了服务体验。此外，在文献［7］中无人机作为移动基站，不仅为本地设备提供计算服务，还可为临近小范围覆盖的基站提供计算服务。然而，它们均缺乏对差异化计算需求的充分考虑。在任务较少时，无人机无须考虑中继传输问题，从而可以降低通信能耗。任务增多时，无人机同时充当中继和边缘计算服务器，灵活应对需求，最大限度优化资源配置。

同时，卸载任务的及时处理也至关重要。文献［5-7］未充分考虑任务的执行时间，而文献［8-10］则简化并统一了处理任务的截止日期。在物联网场景中，单一模型并未充分考虑到单个用户卸载任务执行的延迟问题。因此，本文考虑无人机计算能力的限制，引入了计算任务回传时延约束，以确保无人机在通信范围内及时将处理结果返回给移动用户。此外，由于单架无人机的计算能力有限，因此，只有通过部署多架无人机协同服务用户，才能更好地满足每个用户卸载计算任务的需求。

另外，在应急通信场景中频谱资源稀缺问题日益突显。非正交频分多址（Non-orthogonal multiple access，NOMA，）技术通过功率分配和串行干扰消除技术，成功实现了对时间、频率等物理资源的共享和复用^［11］。该技术为克服频谱资源有限的问题提供了解决方案。文献［12］和文献［13］分别提出了联合优化算法，旨在通过优化无人辅助NOMA系统中的资源分配、用户卸载决策和轨迹规划，以尽可能提高服务质量。此外，文献［14］考虑了UAV和基站采用NOMA技术为用户提供服务的场景。在考虑总功率、总带宽和无人机高度等约束下，提出了一个最大最小速率优化问题。然而，NOMA技术的引入通常伴随着信道干扰的增加^［15］，如何降低干扰也成为一项严峻的挑战。

针对上述问题，本文构建了一个多无人机辅助用户的NOMA中继系统，主要贡献如下。

（1）针对偏远地区或应急场景下无法提供可靠通信和服务的问题，建立了多无人机辅助用户的NOMA中继系统。该系统综合考虑了信息因果、发射功率以及回传时延等约束因素，通过联合优化卸载策略、发射功率、资源分配和轨迹规划，以最小化全系统能耗。

（2）该优化问题是一个长期随机优化问题，具有非凸性、非线性和多变量耦合的特点。利用李雅普诺夫（Lyapunov）优化理论，将该问题转化为用户卸载决策、传输功率和频率分配、轨迹规划3个确定性子问题，极大地降低了问题处理难度。

（3）基于3个子问题，提出了一种交替迭代优化算法：首先，基于最大权匹配算法获得了用户卸载决策；其次，利用辅助变量法、连续凸逼近和拉格朗日对偶得到了传输功率和频率分配的闭式解；最后，将轨迹规划问题转化为凸问题，并使用CVX工具进行求解。相较于基准算法，本文方案在满足系统长期稳定性约束下，显著降低了系统能耗。

1 系统模型

多无人机辅助用户的NOMA中继系统如图1所示，整个系统由K架搭载服务器的无人机、M个用户以及基站构成。其中，用户可以将任务卸载给无人机，由无人机处理或转发给远端基站处理。定义\forall m\in {\rm M}\triangleq \{\mathrm{1,2},\cdots ,M\}，\forall k\in {\rm K}\triangleq \{\mathrm{1},\mathrm{2},\cdots ,K\}分别表示用户和无人机集合。将任务完成时间T均匀离散成N个时隙，各时隙间隔\tau =T/N足够小，以确保各时隙内无人机和用户的位置基本恒定。考虑一个三维欧几里得空间，定义q_m[n]=(x_m[n],y_m[n])表示n\in {\rm N}\{\mathrm{1,2},\cdots ,N\}时隙用户m的水平坐标，q_{\mathrm{B}}[n]=(x_{\mathrm{B}}[n],y_{\mathrm{B}}[n])为远处基站的水平坐标，q_k[n]=(x_k[n],y_k[n])为无人机k的水平坐标，最大飞行速度为V_{\mathrm{m}\mathrm{a}\mathrm{x}}。此外，无人机之间需维持最小安全飞行距离，以免碰撞。则无人机应满足的移动性约束条件为：

C\mathrm{1}:q_k[\mathrm{1}]=q_{k,I},q_k[N]=q_{k,F},\forall k\in K（1）

C\mathrm{2}:||q_k[n+\mathrm{1}]-q_k[n]||\le V_{\mathrm{m}\mathrm{a}\mathrm{x}}\tau ,\forall k,n（2）

C\mathrm{3}:||q_k[n]-q_{k\mathrm{\text{'}}}[n]||\ge D_{\mathrm{m}\mathrm{i}\mathrm{n}},\forall n,k\ne k\mathrm{\text{'}}（3）

式中：q_{k,I}和q_{k,F}分别为无人机k的初始位置和终止位置；D_{\mathrm{m}\mathrm{i}\mathrm{n}}为最小安全距离。同时，假设所有无人机均以固定高度H（单位：m）离地飞行。

1.1　通信模型

根据文献［16］，对于空-地信道，无人机与用户和基站之间的通信链路均视为视距链路，则无人机k到用户m以及基站之间的信道增益分别为：

h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]=h_{\mathrm{0}}(d_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}{[n])}^{-\mathrm{2}}=\frac{h_{\mathrm{0}}}{||q_k[n]-q_m[n]|{|}^{\mathrm{2}}+H^{\mathrm{2}}}（4）

h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]=h_{\mathrm{0}}(d_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}{[n])}^{-\mathrm{2}}=\frac{h_{\mathrm{0}}}{||q_k[n]-q_{\mathrm{B}}[n]|{|}^{\mathrm{2}}+H^{\mathrm{2}}}（5）

式中：d_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]与d_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]分别为n时隙无人机k与机动用户m和基站间的距离；h_{\mathrm{0}}为参考距离d_{\mathrm{0}}=\mathrm{1}\mathrm{m}处的信道功率增益。

无论用户进行局部计算还是卸载计算，其任务都必须保证完成。用\alpha \in {\{\mathrm{0,1}\}}^{M\times K}表示n时隙的用户接入矩阵，每个元素{\alpha }_{m,k}[n]\in \{\mathrm{0,1}\}表示二进制指标，定义{\alpha }_{m,k}[n]=\mathrm{1}表示无人机k在n时隙为用户m提供服务，否则{\alpha }_{m,k}[n]=\mathrm{0}。受通信距离限制，只有在满足距离约束时，无人机k才会成为用户m的潜在卸载候选对象。假设在任意时隙内，每个用户至多关联一个无人机。因此，调度约束满足：

C\mathrm{4}:\sum _{k\in K}{\mathit{\alpha }}_{m,k}[n]\le \mathrm{1},\forall m,n（6）

C\mathrm{5}:{\mathit{\alpha }}_{m,k}[n]\in \{\mathrm{0,1}\},\forall m,n（7）

为了避免无人机之间的干扰，本文将系统带宽B_{\mathrm{t}\mathrm{o}\mathrm{l}}均匀划分为K个正交的无线子信道，子信道的带宽为W=B_{\mathrm{t}\mathrm{o}\mathrm{l}}/K，各个子信道之间不会产生干扰。假设不同的无人机独占一个子信道，在每个时隙对卸载在同一无人机上的用户按照信道增益做升序排列。定义t时隙内无人机k辅助计算的用户列表为：

\begin{array}{c}\mathrm{l}\mathrm{i}\mathrm{s}{\mathrm{t}}_k[n]=\{l_{\mathrm{1}}^k[n],l_{\mathrm{2}}^k[n],\cdots ,\\ l_u^k[n],\cdots ,l_U^k[n]\},l_u^k(t)\in G_k\end{array}（8）

式中：U为t时隙内利用无人机k完成计算任务的用户总数；l_u^k[n]为与无人机k进行通信的第u个最小信道增益的用户索引；G_k为将计算任务卸载至无人机k的用户联盟。

在每个时隙，无人机首先对多个用户的叠加信号按照信道增益降序进行顺序解码。当无人机k解码至用户u的信号l_u^k[n]时，在NOMA组中强于用户u的信号被一一解码并从叠加信号中移出，因此，这些信号不会对用户u造成干扰，而剩余用户l_{\mathrm{1}}^k[n]~l_{u-\mathrm{1}}^k[n]均视为干扰信息。因此，用户l_u^k[n]向无人机k卸载的最大数据量D_{l_u^k,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]可以表示为^［17］：

D_{l_u^k,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]=W\tau \mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}(\mathrm{1}+\frac{P_{l_u^k,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]h_{l_u^k,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]}{({\vartheta }_{\mathrm{1}}+{\delta }^{\mathrm{2}})})（9）

式中：{\vartheta }_{\mathrm{1}}={\displaystyle \sum _{i=\mathrm{1}}^{u-\mathrm{1}}}{p}_{l_i,k}^{D\mathrm{2}U}[n]h_{l_i,k}^{D\mathrm{2}U}[n]为由剩余用户所造成的干扰；P_{l_u,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]为用户l_u^k[n]向无人机k卸载任务的发射功率；{\delta }^{\mathrm{2}}为加性高斯白噪声的功率谱密度。

同理，用户m通过无人机k中继至基站的最大任务量为：

R_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]={\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+\frac{P_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}{{\delta }^{\mathrm{2}}})（10）

式中：P_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]为用户m通过无人机k向基站卸载任务的传输功率。

1.2　队列积压模型

在第n时隙内，用户m的任务特征可以用二元组表示。其中，I_m[n]为t时隙用户m的任务到达量（服从泊松分布），c_m为执行1 bit计算任务所需CPU数。每个用户都维护一个任务缓冲区，用于存储和处理到达任务。

用Q_m[n]表示用户m在n时隙上的任务队列积压。在n时刻，定义用户m在本地处理的任务量为D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]，并向无人机卸载的任务量为R_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]。因此，Q_m[n]的演变方程可以表示为：

Q_m[n+\mathrm{1}]=[Q_m[n]-D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}{[n]-\sum _{k\in K}{R}_{m,k}^{D\mathrm{2}U}[n]]}^{+}+I_m[n]（11）

式中：{[·]}^{+}=\mathrm{m}\mathrm{a}\mathrm{x}\{·,\mathrm{0}\}。

为了保证队列的稳定性，要求队列积压在时间平均意义上需满足：

C\mathrm{6}:{\overline{Q}}_m=\underset{N\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}\frac{\mathrm{1}}{N}{\displaystyle \sum _{n=\mathrm{0}}^{N-\mathrm{1}}}E\{Q_m[n]\}\le \infty（12）

1.3　计算模型

任务执行过程的能耗主要包括3个方面：任务计算能耗、通信能耗以及无人机飞行能耗。

1.3.1　计算能耗

假设地面用户和无人机上都采用动态电压和频率放缩技术以自适应控制其计算频率。那么n时隙用户m在本地和在无人机上的计算能耗分别为：

E_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]=\kappa f_{l,m}^{\mathrm{3}}[n]\tau（13）

E_{k,m}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]=\kappa f_{k,m}^{\mathrm{3}}[n]\tau（14）

式中：\kappa为芯片的有效电容系数；f_{l,m}^{}[n]为n时隙用户m的CPU处理频率；f_{k,m}^{}[n]为n时隙无人机k分配给用户m的计算频率。

同时，n时隙用户m在本地执行的计算任务量可表示为：

D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]=\frac{f_{l,m}[n]\tau }{c_m}（15）

然而，无人机上部署的MEC服务器为用户提供并行服务。无人机k上计算资源为f_k[n]，最大计算资源为f_{k,\mathrm{m}\mathrm{a}\mathrm{x}}，可得计算资源满足以下约束：

C\mathrm{7}:\sum _{m\in M}{\alpha }_{m,k}[n]f_{k,m}[n]\le f_{k,\mathrm{m}\mathrm{a}\mathrm{x}}（16）

此外，鉴于基站通常依赖电网供电且具备丰富的计算资源，本文选择忽略基站的能耗和时延^［5］。

1.3.2　通信能耗

用户m卸载任务至无人机k的传输能耗E_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]以及无人机中继至基站的传输能耗E_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]分别表示为：

E_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]={\alpha }_{m,k}[n]P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]\tau（17）

E_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]={\alpha }_{m,k}[n]P_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]\tau（18）

由于无人机通信范围有限，一旦用户超出其覆盖范围，便无法接收任务回传的结果，因此，无人机必须及时处理并回传用户上传的任务。此外，考虑到无人机的计算资源有限，当无法满足回传时延要求时，必须将计算任务卸载到基站上，即用户m卸载的任务量需要与无人机自身计算任务量以及卸载到基站的任务量的总和相同，则用户m计算任务的回传时延约束为：

R_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]=R_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]+\frac{f_{k,m}[n]\tau }{c_m}（19）

因此，根据香农公式以及式（19），可将式（18）中无人机k的上传能耗重写为：

\begin{array}{l}{\tilde{E}}_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]=\frac{\tau {\sigma }^{\mathrm{2}}}{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}({\mathrm{2}}^{\frac{R_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}{W\tau }}-\mathrm{1})=\frac{\tau {\sigma }^{\mathrm{2}}}{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}·\\ \left({\mathrm{2}}^{\frac{{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+\frac{P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]}{({\vartheta }_{\mathrm{1}}+{\delta }^{\mathrm{2}})})-D_{k,m}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]}{W\tau }}-\mathrm{1}\right)\end{array}（20）

此外，由于结果返回的数据量较少，故本文中未考虑计算结果的时延和能耗。

1.3.3　飞行能耗

无人机在执行计算任务期间，自身飞行能耗是整体系统中不可忽视的一项。与文献［18］类似，飞行能耗仅与飞行速度有关。若定义飞行速度为：

v_k[n]=\frac{||q_k[n+\mathrm{1}]-q_k[n]||}{\tau }（21）

则无人机k在n时隙的飞行能耗表示为：

E_k^{\mathrm{f}\mathrm{l}\mathrm{y}}[n]=k_f||v[n]|{|}^{\mathrm{2}},n\in N（22）

式中：k_f=M_g\tau /\mathrm{2}，M_g为UAV的质量。

综上，n时隙用户和无人机总能耗分别表示为：

E_{\mathrm{M}\mathrm{D}\mathrm{s}}[n]=\sum _{k\in K}\sum _{m\in M}(E_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]+E_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n])（23）

\begin{array}{l}{E}_{\mathrm{U}\mathrm{A}\mathrm{V}}[n]=\sum _{k\in K}\{\sum _{m\in M}(E_{k,m}^{comp}[n]+{\tilde{E}}_{k,m}^{U\mathrm{2}B}[n])+\\ \left.\eta E_k^{fly}[n]\right\}\end{array}（24）

式中：\eta是为了缓解飞行能耗差异而引入的惩罚系数。

因此，n时隙系统的总能耗可表示为：

E_{\mathrm{s}}[n]={\omega }_m{E}_{\mathrm{M}\mathrm{D}\mathrm{s}}[n]+{\omega }_k{E}_{\mathrm{U}\mathrm{A}\mathrm{V}}[n]（25）

式中：{\omega }_m\in [\mathrm{0,1}]和{\omega }_k\in [\mathrm{0,1}]分别为用户和无人机能耗的权重因子，且{\omega }_m+{\omega }_k=\mathrm{1}。

1.4　问题描述

在本文中，定义用户卸载策略\mathit{\alpha }、用户上传功率\mathit{P}=\{{\mathit{P}}_{}^{\mathrm{D}\mathrm{2}\mathrm{U}}\}、计算资源分配\mathit{f}=\{{\mathit{f}}_{l,m},{\mathit{f}}_{k,m}\}和轨迹\mathit{q}表示优化的向量参数集合。具体构建优化问题如下：

\left\{\begin{array}{l}\underset{\mathit{{\rm X}}[n]}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{N}{\displaystyle \sum _{n=\mathrm{0}}^{N-\mathrm{1}}}{E}_{\mathrm{s}}[n]\\ \mathrm{s}.\mathrm{t}.C\mathrm{1}~C\mathrm{7}\\ C\mathrm{8}:\mathrm{0}\le D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]+\sum _{k\in K}{R}_{m,k}^{D\mathrm{2}U}[n]\le Q_m[n]\\ C\mathrm{9}:P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]\le P_{\mathrm{D}\mathrm{2}\mathrm{U}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\\ C\mathrm{10}:\mathrm{0}\le f_{\mathrm{l},m}[n]\le f_{\mathrm{l}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right.（26）

式中：\mathit{X}[n]=\{\mathit{\alpha },\mathit{P},\mathit{f},\mathit{q}\}。C\mathrm{1}为无人机的初始位置和终止位置约束；C\mathrm{2}为无人机的最大飞行速度约束；C\mathrm{3}为无人机间的防碰撞约束；C\mathrm{4}和C\mathrm{5}为用户设备的调度约束；C\mathrm{6}为系统的长期队列稳定性约束；C\mathrm{7}为无人机资源分配约束；C\mathrm{8}为本地计算任务和卸载任务的总和不能超过每个时隙上用户到达任务队列的积压；C\mathrm{9}为本地资源分配约束；C\mathrm{10}为传输功率约束。

2 优化问题求解

问题式（26）是一个长期随机优化问题，由于时隙和变量之间的耦合使得问题通常难以处理。为解决这一问题，本文提出了一种基于Lyapunov的两阶段在线资源分配方案。首先，运用Lyapunov优化理论构建漂移惩罚函数，将目标优化问题转化为只依赖于当前时隙的问题。其次，将优化问题分解为3个子问题，采用交替迭代优化算法求解。

2.1　队列分析与问题转换

定义系统队列状态为\theta [n]=\{Q_m[n]\}{|}_{m\in M}。则Lyapunov函数为：

U(\theta [n])=\frac{\mathrm{1}}{\mathrm{2}}\sum _{m\in M}(Q_m^{\mathrm{2}}[n])（27）

为了保证系统的稳定性，定义Lyapunov漂移\mathrm{\Delta }U(\theta [n])为：

\mathrm{\Delta }U(\theta [n])=\mathrm{E}[U(\theta [n+\mathrm{1}])-U(\theta [n])|\theta [n]]（28）

进而可以定义权衡系统稳定性和能耗间的drift-plus-penalty函数为：

D[\theta [n]]=\mathrm{\Delta }U[\theta [n]]+V\mathrm{E}[E_s[n]|\theta [n]]（29）

式中：V\ge \mathrm{0}是一个控制参数，用于衡量系统能耗和队列稳定性。

通过最小化漂移加惩罚函数的上限，即可实现保证系统稳定下的能耗最小化。因此，通过合理的放缩可得drift-plus-penalty函数的上界为：

\begin{array}{l}D[\theta [n]]\le C+\mathrm{E}\left\{\sum _{m\in M}{Q}_m[n](I_m[n]-\right.\\ \left.D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]-\sum _{k\in K}{R}_{m,k}^{D\mathrm{2}U}[n])\right\}+V\mathrm{E}[E_{\mathrm{M}\mathrm{D}\mathrm{s}}[n]+\\ V{\omega }_k{E}_{\mathrm{U}\mathrm{A}\mathrm{V}}[n]|\theta [n]]\end{array}（30）

式中：C为存在上界的常数。

忽略上界中的常数项， 将原问题转化为每时隙的优化问题，即：

\left\{\begin{array}{l}\underset{{\rm X}[n]}{\mathrm{m}\mathrm{i}\mathrm{n}}-\sum _{m\in M}{Q}_m[n](D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]+\sum _{k\in K}{R}_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n])+\\ \begin{array}{ccc}& & \end{array}V{\omega }_m{E}_{\mathrm{M}\mathrm{D}\mathrm{s}}[n]+V{\omega }_k{E}_{\mathrm{U}\mathrm{A}\mathrm{V}}[n]\\ \mathrm{s}.\mathrm{t}.C\mathrm{1}~C\mathrm{10}\end{array}\right.（31）

问题式（31）是变量高度耦合的非线性规划问题。为了有效求解该问题，将其分解为为用户卸载策略、传输功率和频率、无人机轨迹3个子问题交替优化，具体求解如下。

2.2　子问题求解方案

2.2.1 卸载决策

在所提系统中，用户与无人机的距离会严重影响系统性能^［19］。当\{\mathit{P},\mathit{f},\mathit{q}\}给定时，用户的卸载决策优化问题可以利用基于图论中的最大权匹配理论进行建模和求解，通过缩短各NOMA组中用户到无人机的距离，改善用户的信道条件。无人机与用户之间的距离关系可以表示为加权二部图G=(V,E)，其中V=({\rm K},{\rm M})；边集为E=\{E_{m,k}(\forall m\in M,k\in K)\}，边的权重为节点间距离的相反数即{\omega }_{m,k}=-||q_m[n]-q_k[n]||。通过假设虚拟节点的方法，将无人机节点按照每时隙的最大服务能力进行拓展，构成新的加权二部图G^{\text{'}}=(V^{\text{'}},E^{\text{'}})，其中V^{\text{'}}=({{\rm K}}^{\text{'}},{\rm M})，K^{\text{'}}=(k_{\mathrm{1}}^{\mathrm{1}},k_{\mathrm{1}}^{\mathrm{2}},\cdots ,k_{\mathrm{1}}^M,\cdots ,k_A^M)，此时用户的卸载决策等价于图论中的最大权匹配问题，可由Kuhn-Munkras算法快速求解。

2.2.2 传输功率以及频率优化

当给定\{\mathit{\alpha },\mathit{q}\}时，用户的传输功率以及用户和无人机的频率优化问题可表示为：

\left\{\begin{array}{l}\underset{\mathbf{P},\mathbf{f}}{\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{k\in K}\sum _{m\in M}\{-Q_m[n]{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+\\ P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]{\vartheta }_{\mathrm{2}})+V{\omega }_m\tau ({\alpha }_{m,k}[n]P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]+\\ \kappa f_{\mathrm{l},m}^{\mathrm{3}}[n])+V{\omega }_k\kappa f_{k,m}^{\mathrm{3}}[n]\tau +\\ {\vartheta }_{\mathrm{3}}\left({\mathrm{2}}^{\frac{{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+P_{m,k}^{D\mathrm{2}U}[n]{\vartheta }_{\mathrm{2}})-\frac{f_{k,m}[n]\tau }{c_m}}{W\tau }}-\mathrm{1}\right)-\\ Q_m(t)\frac{f_{\mathrm{l},m}[n]\tau }{c_m}\}\\ \mathrm{s}.\mathrm{t}.C\mathrm{7},C\mathrm{8},C\mathrm{9},C\mathrm{10}\end{array}\right.（32）

式中：{\vartheta }_{\mathrm{2}}=\frac{h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]}{{\vartheta }_{\mathrm{1}}+{\delta }^{\mathrm{2}}}，{\vartheta }_{\mathrm{3}}=\frac{V{\omega }_k\tau {\sigma }^{\mathrm{2}}}{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}。

上述问题是一个非凸优化问题。为了解决这一问题，首先处理目标函数中的非凸指数项问题。为了便于处理，引入辅助变量{\psi }_{k,m}[n]：

{\mathrm{2}}^{\frac{{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]{\vartheta }_{\mathrm{2}})-D_{k,m}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]}{W\tau }}\le {\psi }_{k,m}[n]（33）

由于式（33）的指数项存在P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]和f_{k,m}[n]的耦合关系，难以求解，于是将式（33）进一步化简得：

\begin{array}{l}{\alpha }_{m,k}[n]\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]{\vartheta }_{\mathrm{2}})\le \\ \mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}{\psi }_{k,m}[n]+\frac{f_{k,m}[n]}{W{c}_m}\end{array}（34）

显然，式（34）的左侧中\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]{\vartheta }_{\mathrm{2}})函数关于P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]是非凸的，定义其下界为{\mathrm{ℒ}}_{\mathrm{1}}(P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}})。对于第j次迭代产生的局部点P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},j}[n]，根据一阶泰勒展开式f(x)\ge f(x_{\mathrm{0}})+\frac{\partial f(x)}{\partial x_{\mathrm{0}}}(x-x_{\mathrm{0}})可得其下界为：

\begin{array}{l}{\mathrm{ℒ}}_{\mathrm{1}}(P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}})=\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}(\mathrm{1}+P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},j}[n]{\vartheta }_{\mathrm{2}})+\\ \frac{h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]}{({\vartheta }_{\mathrm{3}}+P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},j}[n]h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n])\mathrm{l}\mathrm{n}\mathrm{2}}·\\ (P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]-P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},j}[n])\end{array}（35）

式中：{\vartheta }_{\mathrm{3}}={\vartheta }_{\mathrm{1}}+{\delta }^{\mathrm{2}}。

将式（35）代入式（34），可将式（33）转化为：

C\mathrm{11}:{\alpha }_{m,k}[n]{\mathrm{ℒ}}_{\mathrm{1}}(P_{\mathrm{m},\mathrm{k}}^{\mathrm{D}\mathrm{2}\mathrm{U}})\le \mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}{\psi }_{k,m}[n]+\frac{f_{k,m}[n]}{W{c}_m}（36）

此外，引入辅助变量后，根据约束C\mathrm{9}中P_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]的范围约束可得：

P_{k,m}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]=\frac{{\sigma }^{\mathrm{2}}}{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}({\psi }_{k,m}[n]-\mathrm{1})\le P_{\mathrm{U}\mathrm{2}\mathrm{B}}^{\mathrm{m}\mathrm{a}\mathrm{x}}（37）

由此可知，辅助变量{\psi }_{k,m}[n]的取值范围为：

C\mathrm{12}:\mathrm{1}\le {\psi }_{k,m}[n]\le \frac{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]P_{\mathrm{U}\mathrm{2}\mathrm{B}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}{{\sigma }^{\mathrm{2}}}+\mathrm{1}（38）

此时，问题式（32）的目标函数已经转化为凸函数，并且非凸约束C\mathrm{8}可整理为：

C\mathrm{13}:\mathrm{0}\le D_{m,\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]+\sum _{k\in K}{\mathrm{ℒ}}_{\mathrm{1}}(P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}})·{\alpha }_{m,k}[n]W\tau \le Q_m[n]（39）

优化问题式（32）已被转化为如下优化问题：

\left\{\begin{array}{l}\underset{\mathbf{P},\mathbf{f}}{\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{k\in K}\sum _{m\in M}\{-Q_m[n]{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+\\ P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]{\vartheta }_{\mathrm{2}})+V{\omega }_m\tau ({\alpha }_{m,k}[n]P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]+\\ \kappa f_{\mathrm{l},m}^{\mathrm{3}}[n])+V{\omega }_k\kappa f_{k,m}^{\mathrm{3}}[n]\tau +{\vartheta }_{\mathrm{3}}({\psi }_{k,m}[n]-\\ \mathrm{1})-Q_m(t)\frac{f_{\mathrm{l},m}[n]\tau }{c_m}\}\\ \mathrm{s}.\mathrm{t}.C\mathrm{7},C\mathrm{9},C\mathrm{11},C\mathrm{12},C\mathrm{13}\end{array}\right.(\mathrm{40})

记最优辅助变量、无人机最优计算频率、本地最优计算频率以及最优传输功率分别表示为{\psi }_m^{\mathrm{*}}[n]\mathrm{、}{f}_{k,m}^{\mathrm{*}}[n]\mathrm{、}{f}_{\mathrm{l},m}^{\mathrm{*}}[n]\mathrm{、}{P}_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},\mathrm{*}}[n]。定义{\lambda }_{m,k}^{\mathrm{1}}\mathrm{、}{\lambda }_{m,k}^{\mathrm{2}}\mathrm{、}{\lambda }_{m,k}^{\mathrm{3}}分别表示与约束条件C\mathrm{7}\mathrm{、}C\mathrm{11}和C\mathrm{13}相关的非负拉格朗日乘子，得到问题最优解如下：

{\psi }_{k,m}^{\mathrm{*}}[n]={[\frac{{\lambda }_{m,k}^{\mathrm{2}}}{\mathrm{l}\mathrm{n}\mathrm{2}{\vartheta }_{\mathrm{3}}}]}_{\mathrm{1}}^{\frac{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]P_{\mathrm{U}\mathrm{2}\mathrm{B}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}{{\sigma }^{\mathrm{2}}}+\mathrm{1}}（41）

f_{k,m}^{\mathrm{*}}[n]=\sqrt[]{\frac{{\lambda }_{m,k}^{\mathrm{2}}-{\lambda }_{m,k}^{\mathrm{1}}{\alpha }_{m,k}[n]c_{m}W}{\mathrm{3}V{\omega }_k\kappa \tau c_{m}W}}（42）

f_{l,m}^{\mathrm{*}}[n]={[\sqrt[]{\frac{Q_m[n]-{\lambda }_{m,k}^{\mathrm{3}}}{\mathrm{3}V{\omega }_k\kappa c_m}}]}_{\mathrm{0}}^{f_l^{\mathrm{m}\mathrm{a}\mathrm{x}}}（43）

\begin{array}{c}{P}_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},\mathrm{*}}[n]=\\ {[\frac{Q_m[n]W\tau {\theta }_{\mathrm{4}}}{V{\omega }_m\tau \mathrm{l}\mathrm{n}\mathrm{2}{\theta }_{\mathrm{4}}+({\lambda }_{m,k}^{\mathrm{2}}+{\lambda }_{m,k}^{\mathrm{3}}W\tau )h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]}-\frac{\mathrm{1}}{{\theta }_{\mathrm{2}}}]}_{\mathrm{0}}^{P_{\mathrm{D}\mathrm{2}\mathrm{U}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}\end{array}（44）

式中：{[x]}_a^b=\mathrm{m}\mathrm{i}\mathrm{n}\{b,\mathrm{m}\mathrm{a}\mathrm{x}\{a,\mathrm{0}\}\}；{\theta }_{\mathrm{4}}={\theta }_{\mathrm{3}}+P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U},j}[n]h_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]。

在获得局部计算资源分配以及传输功率分配后，使用次梯度法更新对偶变量。

2.2.3 无人机轨迹优化

当给定\{\mathit{\alpha },\mathit{P},\mathit{f}\}时，无人机轨迹优化问题可表示为：

\left\{\begin{array}{l}\underset{\mathbf{q}}{\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{k\in K}\sum _{m\in M}\{-Q_m[n]{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+\\ {\gamma }_{\mathrm{1}}{h}_{m,k}^{D\mathrm{2}U}[n])+\frac{V{\omega }_k\tau {\sigma }^{\mathrm{2}}}{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}·\\ ({\mathrm{2}}^{\frac{{\alpha }_{m,k}[n]W\tau \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+{\gamma }_{\mathrm{1}}{h}_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n])-D_{k,m}^{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}[n]}{W\tau }}-\\ \mathrm{1})+\eta V{\omega }_k{k}_f||v[n]|{|}^{\mathrm{2}}\}\\ \mathrm{s}.\mathrm{t}.C\mathrm{1}~C\mathrm{3},C\mathrm{8},C\mathrm{9}\end{array}\right.（45）

式中：{\gamma }_{\mathrm{1}}=P_{m,k}^{\mathrm{D}\mathrm{2}\mathrm{U}}[n]/({\vartheta }_{\mathrm{1}}+{\delta }^{\mathrm{2}})。

问题式（45）与问题式（32）的求解方式类似。首先通过引入辅助变量{\phi }_{k,m}[n]和\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{1}+{\gamma }_{\mathrm{1}}{h}_{m,k}^{D\mathrm{2}U}[n])\le {\xi }_{k,m}[n]简化指数函数问题为：

C\mathrm{14}:{\xi }_{k,m}[n]\le \mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{2}}{\phi }_{k,m}[n]+f_{k,m}[n]/c_{m}W（46）

接下来处理非凸约束。首先处理由引入辅助变量产生的非凸约束问题。同理，通过引入辅助变量{\phi }_{k,m}[n]，约束C\mathrm{9}可以转化为：

C\mathrm{15}:\mathrm{1}\le {\phi }_{k,m}[n]\le \frac{{\gamma }_{\mathrm{3}}}{||q_k[n]-q_{\mathrm{B}}|{|}^{\mathrm{2}}+H^{\mathrm{2}}}+\mathrm{1}（47）

式中：{\gamma }_{\mathrm{3}}=\frac{P_{\mathrm{U}\mathrm{2}\mathrm{B}}^{\mathrm{m}\mathrm{a}\mathrm{x}}{h}_{\mathrm{0}}}{{\sigma }^{\mathrm{2}}\tau }。

通过移项和整理可将约束C\mathrm{15}转化为：

\begin{array}{l}(||q_k[n]-q_{\mathrm{B}}|{|}^{\mathrm{2}}+{\phi }_{k,m}{[n])}^{\mathrm{2}}+(||q_k[n]-q_B|{|}^{\mathrm{2}}-\\ {\phi }_{k,m}{[n])}^{\mathrm{2}}\le \mathrm{4}[{\gamma }_{\mathrm{3}}+||q_k[n]-q_{\mathrm{B}}|{|}^{\mathrm{2}}-\\ H^{\mathrm{2}}{\phi }_{k,m}[n]+H^{\mathrm{2}}]\end{array}（48）

接下来处理非凸约束C\mathrm{3}，对C\mathrm{3}的左侧进行一阶泰勒展开可得其凸约束为：

\begin{array}{l}C\mathrm{16}:D_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}\le -||q_k^r[n]-q_{k\mathrm{\text{'}}}^{\mathrm{r}}[n]|{|}^{\mathrm{2}}+\\ \mathrm{2}(q_k^r[n]-q_{k\mathrm{\text{'}}}^r{[n])}^T(q_k^r[n]-q_{k\mathrm{\text{'}}}^{}[n])\end{array}（49）

因此，优化问题（47）最终被重写为：

\left\{\begin{array}{l}\underset{\mathbf{q}}{\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{k\in K}\sum _{m\in M}\{-Q_m[n]{\alpha }_{m,k}[n]W\tau {\xi }_{k,m}[n]+\\ \frac{V{\omega }_k\tau {\sigma }^{\mathrm{2}}}{h_{k,\mathrm{o}\mathrm{f}\mathrm{f}}^{\mathrm{U}\mathrm{2}\mathrm{B}}[n]}({\phi }_{k,m}[n]-\mathrm{1})+\eta V{\omega }_k{k}_f||v[n]|{|}^{\mathrm{2}}\}\\ \mathrm{s}.\mathrm{t}.C\mathrm{1},C\mathrm{2},C\mathrm{8},C\mathrm{9},C\mathrm{14},C\mathrm{15},C\mathrm{16}\end{array}\right.（50）

显然，式（50）是有关轨迹规划的凸优化问题，可利用Matlab中的CVX求解器求得最优解。3个子问题可联合迭代求解，具体步骤见算法1。

算法1 基于Lyapunov的在线资源分配算法

初始化变量集合\mathrm{X}[\mathrm{0}]=\{\mathrm{f}[\mathrm{0}],\mathrm{P}[\mathrm{0}],\mathrm{q}[\mathrm{0}],\mathrm{R}[\mathrm{0}]\}，{\mathrm{Q}}_{\mathrm{m}}[\mathrm{0}]，{\mathrm{L}}_{\mathrm{m}}[\mathrm{0}]；设置迭代次数\mathrm{n}=\mathrm{0}

1.Repeat

2.固定其他优化变量，求解子问题1，得到D2U最优卸载决策

3.固定{\mathrm{q}}_{\mathrm{k}}[\mathrm{n}]，求解式（42）~（44），得到最优频率和功率分配{\mathrm{f}}_{\mathrm{k},\mathrm{m}}^{\mathrm{*}}[\mathrm{n}],{\mathrm{P}}_{\mathrm{m},\mathrm{k}}^{\mathrm{D}\mathrm{2}\mathrm{U},\mathrm{*}}[\mathrm{n}],{\mathrm{R}}_{\mathrm{n}}^{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},\mathrm{*}}[\mathrm{n}]

4.固定其他优化变量，求解式（50），得到最优无人机飞行轨迹{\mathrm{q}}_{\mathrm{k}}^{\mathrm{*}}[\mathrm{n}]

5.由式（9），（11）和（15）更新任务队列

6. 更新时隙\mathrm{n}=\mathrm{n}+\mathrm{1}

7. until 飞行任务完成时间\mathrm{n}=\mathrm{N}

8.获得最优解{\mathrm{f}}_{\mathrm{k},\mathrm{m}}^{\mathrm{*}}[\mathrm{n}],{\mathrm{P}}_{\mathrm{m},\mathrm{k}}^{\mathrm{D}\mathrm{2}\mathrm{U},\mathrm{*}}[\mathrm{n}],{\mathrm{R}}_{\mathrm{n}}^{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},\mathrm{*}}[\mathrm{n}],{\mathrm{f}}_{\mathrm{l},\mathrm{m}}^{\mathrm{*}}[\mathrm{n}],{\mathrm{q}}_{\mathrm{k}}^{\mathrm{*}}[\mathrm{n}]

3 仿真分析

本文通过Matlab仿真平台验证所提方案的有效性。模拟一个多无人机多用户的计算卸载系统，假设地面机动用户随机分布在\mathrm{100}\mathrm{m}\times \mathrm{100}\mathrm{m}的矩形区域内，且在分布区域内按一定的轨迹运动。规定场景中存在K=\mathrm{3}架无人机，无人机的位置坐标依次表示为q_{\mathrm{1}}[\mathrm{1}]=(\mathrm{0,0})，q_{\mathrm{2}}[\mathrm{1}]=(\mathrm{100,0})，q_{\mathrm{3}}[\mathrm{1}]=(\mathrm{0,100})。无人机任务完成总时间为T=\mathrm{2}\mathrm{s}，平均为N=\mathrm{50}个时隙。系统其余相应参数设置如表1所示^［20］。

为了体现本文算法的优越性，首先分析算法的参数对最终效果的影响，之后将本文算法与以下卸载方案相比较。

（1）本地计算：所有的用户总是在本地执行计算任务。

（2）随机卸载：移动设备在无人机覆盖范围内时，用户随机选择在本地计算还是卸载计算。

（3）固定无人机飞行轨迹：无人机从初始位置出发，沿顺时针方向飞行，飞行轨迹固定为一个半圆。

图2为本文方案中李雅普诺夫控制参数、权重因子以及系统平均队列长度的关系。定义平均系统队列长度为\sum _{m\in M}\sum _{n\in N}(Q_m[n]+L_m[n])/(MN)。在系统优化的过程中，V值扮演着重要的角色，其目的在于平衡队列的稳定性与系统效用之间的复杂关系。从图2中可以观察到，当值较小时，其对V队列积压的影响相对较小，然而，随着V值的增大，队列积压显著增长。这一现象主要归因于在卸载任务和计算资源有限的情况下，为实现能耗的平衡，任务队列需要积压更多任务。同时，从图2可以明确地观察到权重因子{\omega }_k对系统队列的影响，随着{\omega }_k的增加，加权后的无人机计算能耗在整体系统能耗中所占比重逐渐增大，致使更多的用户优先选择本地计算方式，导致任务堆积。

图3展示了本文方案中V值和{\omega }_k对系统效用的影响。从图3中可以明显观察出，随着V值的增大，系统的能耗呈现逐渐降低并趋于稳定的趋势。这是由于较小的V值导致卸载任务的执行开销相对较低，从而使得更多任务选择卸载执行，忽略了系统效用和队列稳定性之间的权衡。然而，随着V值逐渐增大，卸载任务的执行能耗显著增加，导致系统减少了卸载任务的数量，转而更多地选择本地执行计算任务，这在一定程度上忽略了对队列稳定性的考虑。

综合图2和图3可得，通过调整V值设置可以实现优化能耗和队列囤积的最优效果，调整{\omega }_k可以均衡无人机与用户端之间的能耗，本文方案对队列稳定性和长期系统平均能耗的联合优化，可以实现整体最优的结果。

图4呈现了不同方案下系统平均能耗随用户任务到达量变化的趋势。观察图4可知，系统能耗随任务到达量的增加逐渐上升。这归因于通信资源和计算资源在系统中均受到限制，当任务量达到饱和状态时，产生了网络拥塞现象，从而推高了系统能耗。相较之下，其他方案随着用户计算能力的提高，相应卸载到无人机执行的任务量会减少，从而降低了传输成本，进而降低了能耗，但降低幅度较本地计算小。综合4种对比方案的结果，本文方案在优化功耗性能方面表现出显著的优势。这意味着在实际应用中，该方案能够更有效地降低系统的整体能耗水平。

图5展示了在不同控制参数V和到达任务量I_m[n]下，随着用户个数的增加系统平均能耗的变化情况。从图5中可以看出，随着用户数的增加系统能耗呈现指数型增长趋势，这是由于用户数的增加导致了更加激烈的系统资源竞争，使得处理单位数据量所需的能耗也快速增加。此外，还可以发现控制参数越大系统能耗越低，到达数据量越大则系统能耗越高，这一结论与图4中所得结论一致。进一步观察不同参数下能耗的变化情况可以发现，无论控制参数的取值如何，数据量较大的能耗均高于数据量较小的情况。这表明能耗受数据量的影响要大于控制参数，这是由于在任务密集型场景中，过多的任务势必会对系统队列产生强烈的冲击，进而提升能耗，这是无法通过合理调整控制参数而弥补的。但也同时表明，本文算法在多种场景下都具有良好的优化效果，验证了算法的有效性和鲁棒性。

图6为无人机的飞行轨迹，3架无人机在不同的起始位置出发，并最终回到各自的起始位置。观察图6可知，无人机在飞行过程中向相应的用户集群靠近，以优先为邻近用户提供服务。然而，由于计算任务的排队堆积，导致优化轨迹呈现出不平滑、弯曲和曲线较多的特征。当大量计算任务在某一时间段分配给无人机时，会引发无人机位置的跃变。此外，观察到3架无人机的轨迹趋向于遵循防碰撞约束，以避免彼此碰撞。这是本文所考虑的关键因素之一。由于用户在不同时间段产生随机的计算任务，因此，无人机1和无人机2的轨迹图像在某些时间段可能会发生重叠。最后可以看出，为了在有限的任务完成时间内执行飞行任务，无人机会在飞行一段时间后选择高速且直线的轨迹返回到最终位置。

4 结束语

本文研究了具有硬性截止期限的多无人机协同陆地设施辅助的非正交多址接入中继系统卸载决策、资源分配以及轨迹规划问题，通过考虑机动用户任务的硬性截止期限问题，以及用户任务量不同时无人机的两种接入模式，当任务量较少时，无人机仅作为边缘服务器，反之，无人机作为具备中继功能的移动基站，更加切合实际场景的应用要求。之后，通过Lyapunov优化理论、SCA技术和常规凸优化方法将最小化系统能耗的目标优化问题分解为3个子问题逐步求解。仿真结果表明：该算法在能耗方面明显优于基准方案，并保证了算法的收敛性能。

参考文献

[1] Cisco. Cisco visual networking index: forecast and trends, 2017-2022[R/OL]. [2018-11-27].

[2] Zeng Y, Zhang R, Lim T J. Wireless communications with unmanned aerial vehicles: opportunities and challenges[J]. IEEE Communications Magazine, 2016, 54(5): 36-42.

[3] 朱思峰, 赵明阳, 柴争义. 边缘计算场景中基于粒子群优化算法的计算卸载[J]. 吉林大学学报: 工学版, 2022, 52(11): 2698-2705.

[4] Zhu Si-feng, Zhao Ming-yang, Chai Zheng-yi. Computing offloading scheme based on particle swarm optimization algorithm in edge computing scene[J]. Journal of Jilin University (Engineering and Technology Edition), 2022, 52(11): 2698-2705.

[5] Gao B, Zhou Z, Liu F, et al. Winning at the starting line: joint network selection and service placement for mobile edge computing[C]∥IEEE Conference on Computer Communications, Paris, France, 2019: 1459-1467.

[6] 吴启晖, 吴伟. 无人机辅助边缘计算的能量效率最大化算法设计[J]. 通信学报, 2020, 41(10): 15-24.

[7] Wu Qi-hui, Wu Wei. Algorithm design on energy efficiency maximization for UAV-assisted edge computing[J]. Journal on Communications, 2020, 41(10): 15-24.

[8] Fang R, Cui J, Jin S, et al. Optimal node placement and resource allocation for UAV relaying network[J]. IEEE Communications Letters, 2018, 22(4): 808-811.

[9] Liu Y, Xie S, Zhang Y. Cooperative offloading and resource management for UAV-enabled mobile edge computing in power IoT system[J]. IEEE Transactions on Vehicular Technology, 2020, 69(10): 12229-12239.

[10] Zhang T, Xu Y, Loo J, et al. Joint computation and communication design for UAV-assisted mobile edge computing in IoT[J]. IEEE Transactions on Industrial Informatics, 2020, 16(8): 5505-5516.

[11] Li M, Cheng N, Gao J, et al. Energy-efficient UAV-assisted mobile edge computing: resource allocation and trajectory optimization[J]. IEEE Transactions on Vehicular Technology, 2020, 69(3): 3424-3438.

[12] Liu Y, Xiong K, Ni Q, et al. UAV-assisted wireless powered cooperative mobile edge computing: joint offloading, CPU control, and trajectory optimization[J]. IEEE Internet of Things Journal, 2020, 7(4): 2777-2790.

[13] Liu X, Wang J, Zhao N, et al. Placement and power allocation for NOMA-UAV networks[J]. IEEE Wireless Communication Letters, 2019, 8(3): 965-968.

[14] Selim M, Rihan M, Yang Y, et al. On the outage probability and power control of D2D underlaying NOMA UAV-assisted networks[J]. IEEE Access, 2019, 7: 16525-16536.

[15] Zhang X C, Zhang J, Xiong J, et al. Energy-efficient multi-UAV-enabled multiaccess edge computing incorporating NOMA[J]. IEEE Internet of Things Journal, 2020, 7(6): 5613-5627.

[16] 张海君, 张资政, 隆克平. 基于移动边缘计算的 NOMA 异构网络资源分配[J]. 通信学报, 2020, 41(4): 27-33.

[17] Zhang Hai-jun, Zhang Zi-zheng, Long Ke-ping. Resource allocation in NOMA heterogeneous network based on MEC[J]. Journal on Communications, 2020, 41(4): 27-33.

[18] Xue Z, Wu Q, Feng Z.Sum rate maximization in UAV-enabled mobile relay networks[C]∥International Conference on Wireless Communications and Signal Processing, Hangzhou, China, 2018: 1-6.

[19] Zeng Y, Zhang R.Energy-efficient UAV communication with trajectory optimization[J]. IEEE Transactions on Wireless Communication, 2017, 16(6): 3747-3760.

[20] 李国权, 林金朝, 徐勇正, 等.无人机辅助的 NOMA 网络用户分组与分组分配算法[J]. 通信学报, 2020, 41(9): 22-28.

[21] Li Guo-quan, Lin Jin-zhao, Xu Yong-zheng, et al. User grouping and power allocation algorithm for UAV-aided NOMA network[J]. Journal on Communications, 2020, 41(9): 22-28.

[22] Zeng Y, Xu J, Zhang R. Energy minimization for wireless communication with rotary-wing UAV[J]. IEEE Transactions on Wireless Communications, 2019, 18(4): 2329-2345.

[23] Do N, Da Costa D B, Duong T Q,et al.A BNBF user selection scheme for NOMA-based cooperative relaying systems with SWIPT[J]. IEEE Communications Letters, 2017, 21(3): 664-667.

[24] Ren Y, Xie Z, Ding Z, et al. Computation offloading game in multiple UAV-enabled mobile edge computing networks[J]. IET Communications, 2021, 15(10): 1392-1401.

基金资助

陕西省重点研发计划项目(2024NC-YBXM-206)