近年来，随着5G商用化的普及，产业与无线通信领域的深度融合，移动互联网、物联网、智能终端应用和智能可穿戴电子设备取得迅猛发展，网络接入请求数目也达到了前所未有的规模。尽管传统地面通信网络能够满足人们大多数的通信需求，但仍然存在着频谱资源短缺，偏远地区缺乏有效覆盖以及基础设施建设成本过高等问题^[1]。充分利用无人机易于部署、低购置和维护成本、高机动性和悬停能力等各种优势辅助通信网络设计，部署可移动的空中通信平台是解决第5代无线通信网络面临问题的有效途径之一^[2]。

在构建无人机通信系统时，最具挑战性的问题是无人机的部署与控制。研究难点主要集中在3个方面：一是如何高效地部署无人机集群；二是如何设计合理的资源分配方案；三是如何根据环境变化灵活地调整无人机的控制策略。为解决上述问题，文献[3]基于K-means聚类方法提出一种针对多空中基站的部署方案以最大化覆盖用户数量。然而，上述研究仅考虑最大化覆盖用户数目，却忽略了服务范围内多用户之间的干扰问题，进而限制了通信系统性能。文献[4]则在多无人机场景下，提出一种能够有效避免干扰的部署策略，以提高通信系统性能。为解决研究场景规模大且时变性高的问题，学者们开始聚焦于以数据驱动为基础的学习类算法。文献[5]针对解决海上无人机群移动边缘计算网络中通信和计算延迟最小化问题，结合深度强化学习理论提出深度Q网络和深度确定性策略梯度算法来优化无人机群的轨迹和虚拟机的配置。

与此同时，未来6G计划旨在实现无线技术从“万物互联”到“万物智联”的演变，对于能够与环境实时交互的自主智能体的实现提出了更多的需求^[6]。设计有效的无人机在线控制策略成为实现6G的关键步骤之一。文献[7]基于在线贪心算法设计一种无人机辅助计算系统。另有学者采用基于强化学习的方法在避免碰撞的前提下，设计出多无人机系统的最优轨迹规划^[8]。文献[9]提出一种联邦深度学习算法，用于分布式处理无人机通信系统中的信息。然而，基于数据驱动的深度学习算法的性能受限于神经网络模型的训练效果。当训练数据集稀缺或病态时，实际控制决策器性能可能不佳。同时，采用“训练-推理”模式的在线学习算法无法针对更新后的数据集实现实时推理。此外，强化学习算法采用反馈机制可以在一定程度上实现与环境的实时交互，但该算法无法充分利用已知模型信息辅助制定实时控制策略，在一定程度上损失了策略的最优性^[10]。综上所述，针对实时性要求高且环境复杂多变的无人机通信场景，学习类算法很难制定出实时最优控制策略，并且对计算能力和硬件支持过高的要求使得该算法在实际工程中难以大规模应用^[6]。

此外，目前大部分研究将固定翼无人机作为辅助通信场景下的主要研究对象^[11]。然而，考虑到地面用户移动的随机性与无人机的能耗等因素，选择机动性更强且具有悬停能力的多旋翼无人机作为空中移动平台更具发展前景。为充分利用无人机的移动性以获得更高的系统总吞吐量收益，人们对于使用3D MIMO技术产生了相当大的兴趣^[12]。文献[13]探讨无人机在大规模MIMO场景下的应用潜力，并针对在无人机和地面基站端同时部署大规模天线阵列的通信场景，建立了极化失配的损失模型。进一步推导出相应的遍历容量下界的闭式表达式，并设计了使遍历容量最大化的最优天线阵列部署策略。

综上所述，本文针对无人机通信网络中多用户间干扰大、动力学模型复杂度高，以及设计实时最优控制策略难等问题展开研究。首先提出一种三维多用户多输入单输出(three-dimensional multi-user multiple-input single-output，3D MU-MISO)场景下与四旋翼无人机姿态角耦合的通信网络模型。其次，在该场景下基于模型预测控制(model predictive control, MPC)框架设计了一种无人机姿态控制和资源分配的在线控制策略。该策略通过滚动式在线求解有限时间窗长的开环控制问题，实现无人机在向终点飞行的过程中最大化所服务系统的平均频谱效率(average spectral efficiency, ASE)，以规避现有在线控制策略中因数据集和训练过程而对性能产生负面影响的问题。最后，通过设计仿真试验验证了所提出的在线控制策略能够在有限降低系统收益的情况下实现对离散控制方案的近似，并且所提出的姿态角耦合复杂通信模型在3种常见的波束成形方式下均能带来系统自由度增益，从而提高无人机通信系统的ASE。

1 模型描述

本文所研究的四旋翼无人机作为移动空中基站，其提供服务的3D MU-MISO通信网络场景如图 1所示。该场景包括一个四旋翼无人机，其装备了一个由M=M_x×M_y根天线组成的均匀平面天线阵列(uniform planar array，UPA)，以及一组单天线可移动用户I ={1, 2, …, I}。为充分利用多天线阵列所带来的通信系统增益，所研究场景考虑在同一时频资源下为地面多用户同时提供通信服务，采用空分多址接入(space division multiple access, SDMA)方式，为每个用户的数据流都分配了一个独立的空间以避免数据流之间的干扰。除考虑四旋翼无人机的三维坐标外，还需要同时考虑无人机姿态角对通信网络的影响。在满足四旋翼无人机的动力学方程的前提下，设计相应的控制算法，以自适应地平衡通信距离和矩形天线阵列的发射功率。为简化分析，假设所有的用户均停留在地面上且位置已知，并且将无人机的控制过程以T为时隙间隔进行离散化的近似处理，若采样时间足够短，无人机在每个时隙内均可看作静止状态^[14-15]。

本文考虑一个实际的应用场景，位于森林、海上，或是受灾区域等缺少基站覆盖的通信场景，无人机作为可移动的空中基站为地面多用户提供通信服务。请注意，本文研究的重点是装备有UPA的四旋翼无人机在3D MU-MISO场景下的自主在线控制策略设计问题。因此，仅围绕单四旋翼无人机展开且不涉及传统的地面基站。

1.1 无人机动力学模型与能量消耗模型 1.1.1 无人机控制系统模型

四旋翼无人机由4个螺旋桨通过刚性十字架相连组成，其中可分为以相反方向旋转的2组螺旋桨。通过增加或减小总推力，四旋翼系统可实现向不同方向旋转和移动^[16]。具体关于四旋翼无人机的动力学方程请参考文献[17-18]。

根据现代控制理论^[19]，无人机动态系统被建模为一组微分方程，并可以用状态空间模型的形式表示。定义四旋翼无人机在惯性坐标系下在时刻t的位置为$\boldsymbol{q}(t)=[x(t), y(t), z(t)]^{\mathrm{T}} \in \mathbb{R}^3$ 和速度为 $\dot{\boldsymbol{q}}(t)=[\dot{x}(t), \dot{y}(t), \dot{z}(t)]^{\mathrm{T}} \in \mathbb{R}^3$。同样地，定义在时刻 t 机体坐标系下的横滚角、俯仰角和偏航角分别为$\boldsymbol{\varPhi}(t)=[\boldsymbol{\phi}(t), \theta(t), \psi(t)]^{\mathrm{T}}$，对应的角速度为$\dot{\boldsymbol{\varPhi}}(t)=[\dot{\boldsymbol{\phi}}(t), \dot{\boldsymbol{\theta}}(t), \dot{\psi}(t)]^{\mathrm{T}}$。无人机控制系统中的状态向量可以定义为

$ \begin{gathered} \boldsymbol{x}(t)=[x(t), y(t), z(t), \dot{x}(t), \dot{y}(t), \dot{z}(t), \phi(t), \\ \theta(t), \psi(t), \dot{\phi}(t), \dot{\theta}(t), \dot{\psi}(t)]^{\mathrm{T}}, \end{gathered} $

(1)

其中控制向量定义为u(t)=[u₁(t), u₂(t), u₃(t), u₄(t)]^T。

为方便表述，可以得到常微分方程模型的紧凑形式，该模型引入了一个时不变的非线性函数f

$ \dot{\boldsymbol{x}}(t)=f(\boldsymbol{x}(t), \boldsymbol{u}(t)), $

(2)

其中ẋ(t)表示x(t)相对于t的导数。

从MPC的角度来看，使用欧拉法从上述时间不变非线性函数中构建近似的离散时间动态模型，以Δ_T为时间单位的离散动态系统可以表示为

$ \boldsymbol{x}[t+1]=f(\boldsymbol{x}[t], \boldsymbol{u}[t]), $

(3)

其中: x[t]表示在时隙t时离散化动态模型的状态向量，u[t]表示离散时间控制向量。上述变量在一个时隙中被视为常数。本文的后续研究将在离散时间尺度上展开，且表示为[·]。

1.1.2 能量消耗模型

参照文献[20]中所述，考虑无人机装备4个电池供电的无刷马达，对于每一个电机马达n ∈ {1, 2, 3, 4}位于时隙t时的能量损耗可以描述为

$ \begin{aligned} & P_n^{\mathrm{con}}[t]=U_n[t] I_n[t] \\ & =c_4 \omega_n[t]^4+c_3 \omega_n[t]^3+c_2 \omega_n[t]^2+c_1 \omega_n[t]+c_0, \end{aligned} $

(4)

其中：I_n[t]与U_n[t]分别表示电机马达n∈{1, 2, 3, 4}位于时隙t的电流和电压。并且有

$ \begin{gathered} c_0=I_0^2 R_0, \quad c_1=30 K_E I_0 / \pi, \quad c_2=2 C_m R_0 I_0 / K_T, \\ c_3=30 C_m K_E /\left(\pi K_T\right), \quad c_4=C_m^2 R_0 / K_T^2, \end{gathered} $

其中：$K_E \triangleq\left(U_0-I_0 R_0 / K_v U_0\right)$为反电动势常数，$K_T \triangleq 9.55 K_E$为转矩常数，R₀为电机电阻。

因此，本文考虑的四旋翼无人机在时隙t的总能量损耗可以表示为

$ P^{\mathrm{con}}[t]=\sum\limits_{n=1}^4 P_n^{\mathrm{con}}[t] . $

(5)

1.2 无人机通信模型 1.2.1 无线信道模型

根据文献[21]的研究结果可知，在无人机飞行高度为35 m以上且其对地覆盖范围区域的直径小于100 m的情况下，无人机与地面用户节点间通信链路主要由直射信号分量所构成，可合理假设无人机通信模型为莱斯衰落信道。因此，第i个用户的无人机空-地信道增益在每个离散时隙t可以建模为

$ \boldsymbol{h}_i[t]=\sqrt{\beta_i[t]} \boldsymbol{g}_i[t] . $

(6)

式中：β_i[t]表示在t时刻的大尺度信道增益，包括信号衰减、路径损耗和阴影效应；而g_i[t]表示在t时刻小尺度衰落系数。假设每时每刻用户的具体位置可以通过GPS完美获取。g_i[t]遵循莱斯分布，且具有归一化功率$\mathbb{E}\left\|\boldsymbol{g}_i[t]\right\|^2=1$，可由下式得到

$ g_i[t]=\sqrt{\frac{K_a}{K_a+1}} \bar{g}_i[t]+\sqrt{\frac{1}{K_a+1}} \tilde{g}_i[t], $

(7)

其中: ḡ_i[t]代表直射路径信道分量，$\tilde{\boldsymbol{g}}_i[t]$表示非视距信道分量，其元素是独立同分布的随机变量，且服从于$\mathcal{C N}$(0, 1)分布。

1.2.2 姿态角耦合信道模型

基于上文所描述的四旋翼无人机动态系统，本文针对3D MU-MISO场景下与四旋翼无人机姿态角耦合的问题进行建模。该模型考虑无人机控制系统中横滚角，俯仰角和偏航角3个角度状态与通信信道之间的耦合关系。

如图 2所示，由于发射天线阵列的任意旋转会改变无人机与地面用户之间直射链路向量的方位角和俯仰角的大小，因此它们必须相对于旋转后的坐标轴进行计算。使用上标×表示发射天线阵列的旋转坐标系的分量。

根据旋转矩阵的定义^[22]，发射天线阵列经过旋转之后的新的坐标系的单位方向向量可表示为

$ \left(\begin{array}{lll} \hat{x}^{\times} & \hat{y}^{\times} & \hat{z}^{\times} \end{array}\right)=\left(\begin{array}{lll} \hat{x} & \hat{y} & \hat{z} \end{array}\right) \boldsymbol{R}(\phi[t], \theta[t], \psi[t]), $

(8)

其中：$\boldsymbol{R}(\phi[t], \theta[t], \psi[t])$是3×3的旋转矩阵，$\left(\begin{array}{lll} \hat{x} & \hat{y} & \hat{z} \end{array}\right)$为旋转之前坐标系下的单位坐标，$\left(\begin{array}{lll} \hat{x}^x & \hat{y}^x & \hat{z}^x \end{array}\right)$为旋转之后新坐标系下的单位坐标。在上述旋转过后的新坐标系下，接收器(即第i个用户)的位置坐标更新为

$ \begin{gathered} \left(\begin{array}{l} x_i^{\times}[t] \\ y_i^{\times}[t] \\ z_i^{\times}[t] \end{array}\right)=(\boldsymbol{R}(\phi[t], \theta[t], \psi[t]))^{\mathrm{T}}\left(\begin{array}{l} x_i[t] \\ y_i[t] \\ z_i[t] \end{array}\right) \\ =\left(\begin{array}{l} \boldsymbol{R}_{11} x_i[t]+\boldsymbol{R}_{21} y_i[t]+\boldsymbol{R}_{31} z_i[t] \\ \boldsymbol{R}_{12} x_i[t]+\boldsymbol{R}_{22} y_i[t]+\boldsymbol{R}_{32} z_i[t] \\ \boldsymbol{R}_{13} x_i[t]+\boldsymbol{R}_{23} y_i[t]+\boldsymbol{R}_{33} z_i[t] \end{array}\right) . \end{gathered} $

(9)

在无人机控制系统中，相对于发射天线阵列旋转后的坐标轴的方位角和俯仰角分别可以计算为

$ \begin{aligned} \zeta_i^{\times}[t] & =\cos ^{-1}\left(\hat{y}_i^{\times}[t] \cdot \overrightarrow{\boldsymbol{q q}_i}[t]\right) \\ & =\cos ^{-1}\left(\frac{y_i^{\times}[t]}{\sqrt{x_i[t]^2+y_i[t]^2}}\right), \end{aligned} $

(10)

$ \begin{aligned} \iota_i^{\times}[t] & =\cos ^{-1}\left(-\hat{z}_i^{\times}[t] \cdot \overrightarrow{\boldsymbol{q} \boldsymbol{q}_i}[t]\right) \\ & =\cos ^{-1}\left(-\frac{z_i^{\times}[t]}{\sqrt{x_i[t]^2+y_i[t]^2+z_i[t]^2}}\right), \end{aligned} $

(11)

其中: $\overrightarrow{\boldsymbol{q q}_i}[t]$表示在当前时刻t无人机与用户i所连成的向量，ζ_i^×[t]∈[0, 2π]表示无人机与用户i之间在t时刻更新后的方位角(无人机与用户之间的方向向量在无人机自身坐标系x-y面上的投影与无人机自身坐标系下y轴正半轴之间的夹角)，ι_i^×[t]∈[0, π]表示更新后的俯仰角(无人机与用户之间的方向向量与无人机自身坐标系下z轴负半轴之间的夹角)。

随后根据更新后的方位角和俯仰角，得到无人机经过姿态调整之后的转向矢量矩阵为

$ \begin{aligned} & \bar{g}_i^{\times}[t]= \\ & \frac{1}{\sqrt{M}}\left[1, \cdots, \mathrm{e}^{\mathrm{j} \frac{2 \pi}{w} d\left(\left(m_x-1\right) \sin \zeta_i^{\times}[t] \sin u_i^{\times}[t]+\left(m_y-1\right) \sin u_i^{\times}[t] \cos \zeta_i^{\times}[t]\right)}\right. \\ & \left.\cdots, \mathrm{e}^{\mathrm{j} \frac{2 \pi}{w} d\left(\left(M_x-1\right) \sin \zeta_i^{\times}[t] \sin x_i^{\times}[t]+\left(M_y-1\right) \sin u_i^{\times}[t] \cos \zeta_i^{\times}[t]\right)}\right]^{\mathrm{T}} . \end{aligned} $

(12)

从而得到无人机姿态与通信信道质量之间的耦合关系表达式，并代入式(7)中得到在上述模型下的信道增益h_i^×[t]，以及在时隙t内的系统ASER_i^×[t]用于组成后续内容中待控制问题的目标函数

$ \begin{aligned} & R_i^{\times}[t]= \\ & \log _2\left(1+\frac{\eta_i[t]\left|\boldsymbol{h}_i^{\times}[t]^{\mathrm{H}} \boldsymbol{w}_i^{\times}[t]\right|^2}{\sum\limits_{s \in \hat{N}|i|} \eta_s[t]\left|\boldsymbol{h}_i^{\times}[t]^{\mathrm{H}} \boldsymbol{w}_s^{\times}[t]\right|^2+\sigma_i^2}\right), \end{aligned} $

(13)

其中：$n_i[t] \sim \mathcal{C} \mathcal{N}\left(0, \sigma_i^2\right) \in \mathbb{C}$表示用户i在时隙t中的加性复高斯白噪声，σ²代表噪声功率，变量η_i[t]是无人机在时隙t向用户i的发射功率。

1.3 研究问题描述

现有的无人机轨迹策略设计通常采用离散路径规划方法^{[17, 23]}。具体而言，将无人机在运行时间T期间的轨迹离散化为N_T个航路点，并将时间T划分为N_T个等长的时隙，其中每个时隙的持续时间为Δ_T，且使得T=N_TΔ_T。基于四旋翼无人机的动力学方程组，式(13)所描述的耦合信道模型下的ASE，以及式(5)所得到的四旋翼无人机能量消耗模型，在当前迭代时隙t₀下$\mathcal{H}$ ={t₀, t₀+1, …, N_T}的视野范围内考虑各种约束情况下构建出需要解决的最优姿态控制决策问题(optimal control problem, OCP)，如下所示

$ \begin{aligned} & \max _{X, U, \eta} \sum\limits_{t=t_0}^{N_{T^{-1}}}\left[Q_1 \sum\limits_{i=1}^I R_i^{\times}[t]+Q_2\left\|q[t]-q_{\mathrm{des}}\right\|^2+\right. \\ & \left.Q_3 P^{\mathrm{con}}[t]\right]+Q_4\left\|\left(q\left[N_T\right]-q_{\mathrm{des}}\right)\right\|^2, \end{aligned} $

(14)

$ \begin{array}{ll} \text { s. t. } & C_0: \boldsymbol{x}[t+1]=f(\boldsymbol{x}[t], \boldsymbol{u}[t]), \quad t \in \mathcal{H}, \\ & C_1: \boldsymbol{x}\left[t_0\right]=\boldsymbol{x}_0, \\ & C_2: 0 \leqslant u_1[t] \leqslant u_{1 \max }, t \in \mathcal{H}, \\ & C_3:\left|u_n[t]\right| \leqslant u_{n \max }, n=\{2, 3, 4\}, t \in \mathcal{H}, \\ & C_4: h_{\min } \leqslant z[t], |\phi[t]| \leqslant \phi_{\max }, |\theta[t]| \leqslant \theta_{\max }, t \in \mathcal{H}, \\ & C_5:\|\Delta(\dot{x}[t], \dot{y}[t], \dot{z}[t])\| \leqslant \delta a_{\max }, t \in \mathcal{H}, \\ & C_6:\|(x[t], y[t])\| \leqslant V_{\operatorname{hmax}}, \|z[t]\| \leqslant V_{\max }, t \in \mathcal{H}, \\ & C_7: R_i^{\times}[t] \geqslant Q_{\min }, t \in \mathcal{H}, i \in I, \\ & C_8: \sum\limits_{i=1}^l \eta_i[t] \leqslant P_{\max }, t \in \mathcal{H}, \\ & C_9: \boldsymbol{q}\left[N_T\right]=\boldsymbol{q}_{\mathrm{des}}, \end{array} $

其中：$\boldsymbol{X}=\left[\boldsymbol{x}\left[t_0\right]^{\mathrm{T}}, \boldsymbol{x}\left[t_0+1\right]^{\mathrm{T}}, \cdots, \boldsymbol{x}\left[N_T\right]^{\mathrm{T}}\right]^{\mathrm{T}} \in$ $\mathbb{R}^{\left(N_T-t_0+1\right) n_x}$ 和 $\boldsymbol{U}=\left[\boldsymbol{u}\left[t_0\right]^{\mathrm{T}}, \boldsymbol{u}\left[t_0+1\right]^{\mathrm{T}}, \cdots, \boldsymbol{u}\left[N_T-\right.\right.$ $\left.1]^{\mathrm{T}}\right]^{\mathrm{T}} \in \mathbb{R}^{\left(N_T-t_0\right) n_u}$表示当前时隙下OCP所对应的优化变量，包括在视野范围$\mathcal{H}$内的所有状态和控制向量。Q₁, Q₂, Q₃和Q₄为目标函数的权衡因子，通过恰当设定这些值以平衡多个目标。第1项最大化所有航路点N_T内由无人机为用户提供通信服务的总ASE之和；通过控制Q₂，第2项使得无人机在时隙t位置坐标与目标终点之间的差距最小；第3项则限制了无人机导航中的能量消耗；最后1项使得在最后时隙的状态与目标终点的距离最小。

对于约束条件而言，C₀代表四旋翼无人机状态空间模型的离散动态方程，该式为高度非线性等式约束，使得上述控制问题的求解变得十分棘手。C₁给出了系统状态向量的初始值。C₂和C₃用于限制电机旋转的最大角速度。为安全起见，限制C₄被强制执行，以限制无人机的飞行高度、翻转角度和俯仰角度，其中h_min是最低飞行的安全高度，ϕ_max和θ_max分别是ϕ[t]和θ[t]的最大允许值。C₅限制无人机在采样时间内的最大加速度a_max。C₆要求无人机在安全水平速度V_hmax和垂直速度V_vmax以下飞行。C₇保证地面用户的通信服务质量，其中Q_min是每个用户的目标速率阈值。C₈确保无人机在视野范围内的各时隙t中的传输功率不超过最大值P_max。C₉确保无人机准确抵达终点。

2 在线控制策略 2.1 基于MPC的在线控制策略

图 3为基于MPC的无人机控制系统的整体框架。假设在每次迭代开始时，可以通过安置在系统内的传感器获取系统的实时数据，随后使用扩展卡尔曼滤波器(extended Kalman filter, EKF)得到当前系统状态的完美估计值 $\hat{x}_0$，并将其输入至MPC在线控制器计算得到当前时刻系统的控制输入量 $\boldsymbol{u}_0^*\left(\hat{\boldsymbol{x}}_0\right)$，重复上述步骤直至系统达到稳态或者实现控制目的。

2.2 MPC在线控制器设计

基于问题(14)，设定N₀为MPC控制器的预测时间窗长度。区别于现有的在线控制策略，基于MPC的在线控制策略要求在当前MPC迭代时隙t₀求解$\mathcal{H} \triangleq\left\{t_0, \cdots, t_0+N_0\right\}$视野范围内所构成的开环控制问题

$ \begin{aligned} &\max _{\boldsymbol{X}, U, \boldsymbol{\eta}}^{t_0+N_0-1} \sum\limits_{t=t_0} l(\boldsymbol{x}[t], \boldsymbol{u}[t], \boldsymbol{\eta}[t])+E\left(\boldsymbol{x}\left[t_0+N_0\right]\right), \\ &\begin{array}{ll} \text { s.t. } & \forall t \in\{0, \cdots, n-1\}, \\ & \boldsymbol{x}[t+1]=f(\boldsymbol{x}[t], \boldsymbol{u}[t]), \\ & g(\boldsymbol{x}[t], \boldsymbol{u}[t])=0, \\ & h(\boldsymbol{x}[t], \boldsymbol{u}[t]) \leqslant 0, \\ & (\boldsymbol{x}[t], \boldsymbol{u}[t]) \in \mathbb{X} \times \mathbb{U}, \end{array} \end{aligned} $

(15)

其中: 目标函数中l(·)与E(·)分别为代价函数(对应式(14)目标函数中前3项)与终端约束函数(对应于式(14)目标函数中最后1项)，对应于公式目标函数中前3项与最后项。而对于上述非线性规划问题中的约束条件，g(·)描述了无人机控制系统的动态方程组和h(·)描述了所包含的所有不等式约束，并有

$ \begin{aligned} & \mathbb{X}=\left\{\boldsymbol{x}[t] \in \mathbb{R}^{n_x} \mid \boldsymbol{x}_{\min } \leqslant \boldsymbol{x}[t] \leqslant \boldsymbol{x}_{\max }\right\}, \\ & \mathbb{U}=\left\{\boldsymbol{u}[t] \in \mathbb{R}^{n_u} \mid \boldsymbol{u}_{\min } \leqslant \boldsymbol{u}[t] \leqslant \boldsymbol{u}_{\max }\right\} . \end{aligned} $

(16)

随后，根据内点法，引入平滑常数τ>0和松弛变量s_j, j={1, …, n_h}用于改善牛顿法求解时对初值的敏感程度，写出非线性规划问题的KKT条件为

$ \begin{gathered} \nabla F\left(\boldsymbol{w}^*\right)+\nabla G\left(\boldsymbol{w}^*\right) \boldsymbol{\lambda}^*+\nabla H\left(\boldsymbol{w}^*\right) \boldsymbol{\mu}^*=0, \\ G\left(\boldsymbol{w}^*\right)=0, \\ \boldsymbol{\mu}_j^* s_j^*-\tau=0, \quad j \in\left\{1, \cdots, n_h\right\}, \\ h_j\left(\boldsymbol{w}^*\right)+s_j^*=0, \quad j \in\left\{1, \cdots, n_h\right\}, \end{gathered} $

(17)

其中: 优化变量向量w^T=(x[t₀]^T, u[t₀]^T, x[t₀+1]^T, u[t₀+1]^T, …, u[t₀+N₀]^T, x[t₀+N₀]^T)，以及对应的等式对偶乘子为λ^T=(λ₁^T, λ₂^T, …, λ_N0+1^T)和对应的不等式对偶乘子为μ^T=(μ₁^T, μ₂^T, …, μ_nh^T)，n_h为不等式约束的个数。为简化描述，定义如下变量为

$ z:=\left[\begin{array}{c} w \\ \boldsymbol{\lambda} \\ \boldsymbol{\mu} \\ s \end{array}\right], F(z):=\left[\begin{array}{c} \nabla_w \mathcal{L}(w, \boldsymbol{\lambda}, \boldsymbol{\mu}) \\ G(w) \\ \boldsymbol{\mu}_j s_j-\tau \\ h_j(w)+s_j \end{array}\right], $

(18)

其中: $z \in \mathbb{R}^{N_0\left(n_x+n_u+1\right)+2 n_h}, F: \mathbb{R}^{N_0\left(n_x+n_u+1\right)+2 n_h} \rightarrow \mathbb{R}^{N_0\left(n_x+n_u+1\right)+2 n_h}$。随后总结上述非线性求零点问题为

$ F(z)=0 . $

(19)

从初始值z₀开始，牛顿法通过在线性化当前迭代处的非线性方程，生成一系列迭代{z_k}_k=0^∞。

上述是基于内点法求解MPC非线性规划问题的主要推导过程，同时也是为了调整原问题从而使其更加适合于非线性求解器的求解形式。为确保无人机能够在准确抵达终点的过程中尽可能地调整其姿态，为地面用户提供最大的通信服务，需要将式(15)所对应的控制问题分为2个阶段。阶段1为最大化通信服务阶段，即当无人机从起点出发之后且与终点存在一定距离差时，恰当地设定目标函数中各项的权重参数(Q₁, Q₂, Q₃, Q₄)，使无人机能够实现在抵达终点的过程中最大程度地为地面用户提供通信服务，其中权重参数如文献[19]所述，以及通过多次重复实验所确定。阶段2为准确抵达终点阶段，即当无人机即将抵达终点时，调整其在线控制策略，随着控制策略迭代步数的增加，逐步减少MPC所考虑的视野范围N₀，并且在目标函数中增加关于视野末端抵达终点的约束条件，在有限降低系统收益的情况下指导无人机准确抵达终点。基于MPC的无人机轨迹优化算法的完整步骤如算法1所示。

算法1   基于MPC的无人机在线控制策略算法
初始化:
初始化系统待优化变量X, U, η, 预测时间窗长N0, 无人机初始状态x0, 终止阈值∈, 终点坐标, 令MPC最大迭代步数MAXMPCiter, 当前迭代次数MPCiter=0与时钟t=0。
在线:
1:使用EKF获得系统当前MPC迭代时刻t0状态变量的估计值, 即x[t0]。
2:if无人机当前坐标与xdes相差较大时then
3:   在t0时刻, 得非线性规划问题(17)。
4:   采用牛顿法求解问题(19), 若达到最大迭代数或对偶乘子满足‖μ‖ < ∈, 终止迭代。
5:   将求解的ut0输入至实际系统。
6:   将所求得的在t∈{t0+1, t0+2, …, t0+N0}时刻的优化变量解值作为初始值赋值给下一步在t∈{t0, t0+1, …, t0+N0}时刻的系统优化变量。
7:else
8:   在t0时刻, 得非线性规划问题(17)。并且添加等式约束x[t0+N0]=xdes。
9:   采用牛顿法求解问题(19), 若达到最大迭代步数或对偶乘子满足‖μ‖ < ∈, 终止迭代。
10:  令N0=N0-1。
11:  将求解的ut0输入至实际系统。
12:  将所求得的在t∈{t0+1, t0+2, …, t0+N0}时刻的优化变量解值作为初始值赋值给下一步在t∈{t0, t0+1, …, t0+N0}时刻的系统优化变量。
13:end if
14:令MPCiter = MPCiter+1, t=t+1。
15:检查MPCiter < MAXMPCiter或N0>0, 算法终止; 否则返回步骤1继续迭代。

3 仿真结果与分析

本节展示所提出的3D MU-MISO场景下的姿态角耦合通信模型，以及基于MPC的无人机自主姿态控制及资源分配的在线控制策略，通过仿真实验评估该在线控制策略在耦合通信模型下的性能，并对实验结果进行分析。

实验的硬件环境为：系统类型为64位Windows 10操作系统，CPU为AMD Ryzen 5 PRO 4650G 3.700 GHz，内存为16.0 GB，软件为R2020b版本Matlab，使用IPOPT^[24]和CasADi^[25]工具进行实验仿真。

选取100 m×100 m的区域作为仿真场景，该区域内分布有3个固定地面设备与可移动无人机基站设备，假设用户状态信息和无人机状态信息都可以通过传感器准确获得。其用户坐标分别是：(x₁, y₁)=(25, 25)，(x₂, y₂)=(60, 70)和(x₃, y₃)=(70, 50)。四旋翼无人机上装配有M_x=M_y=4的矩阵天线阵列可同时为多个地面用户提供通信服务，同时满足用户的最小通信需求和无人机的最大功率损耗约束。其中无人机从起点q_A=(50, 0, 100)飞至终点q_B=(50, 100, 100)。通信信道假设服从前文所描述的莱斯信道，并且设定天线总发射功率为P_max=100 W。MPC每次迭代所涉及的视野范围N₀=25，本章所阐述的结果都是对莱斯信道的小尺度衰落做50次实现求其平均值。

为简化描述，关于四旋翼无人机各参数设置具体参考文献[18, 20]，而无人机通信系统的仿真参数设置则参考文献[26]，全部仿真参数如表 1所示。

3.1 无人机姿态变化对通信系统性能影响的仿真分析

首先，在固定四旋翼无人机与其服务地面用户的三维坐标的场景下，设置100组无人机姿态角变化组，其中每一组中的横滚角、俯仰角与偏航角均在各自状态变化范围内等间隔采样。将状态数值代入前文提出的通信场景下的信道质量与不同波束成形方式下的ASE计算公式(13)中，从而获得地面各个用户信道质量与ASE的变化图，用于验证姿态角与信道质量之间的非线性耦合关系，说明针对姿态角进行优化的必要性。

图 4展示了无人机姿态改变时，其与地面用户之间无线信道质量的变化以及采用3种常见的波束成形方式所得到的系统ASE变化情况。当固定四旋翼无人机的三维坐标，改变其姿态角时，可以明显观察到无人机与所服务地面用户之间信道增益的显著变化，如图 4(a)所示。同时，采用3种常见的波束成形方式所获得的系统ASE也随之发生显著的改变，如图 4(b)、4(c)和4(d)所示。在空间上，用户2和用户3的相对靠近导致信道之间的相互干扰较大。基于MRT波束成形方式无法有效降低相邻用户之间的干扰信号，从而导致采用该方式的系统ASE为三者之中最低。另外2种波束成形方式则在一定程度上牺牲了算法的简易性，以换取更高的系统ASE。由无人机姿态角改变引起的信道质量和系统ASE的显著变化，强调了建立姿态角与系统通信质量之间耦合关系的重要性。通过优化姿态角可以提高系统ASE，从而获得更好的通信性能。

3.2 无人机在线控制策略仿真分析

其次，通过仿真实验评估基于MPC的在线控制策略的性能，将其与2种基准控制方案相比较，即基准方案1：采用固定轨迹方案，四旋翼无人机沿起点到终点的直线飞行且在飞行过程中优化其姿态角以及天线的发射功率，以实现最大化该方案下所获得的通信系统的ASE；基准方案2：参考文献[18]中所描述的离散路径规划方案，将整条无人机轨迹离散成为N_T个航路点进行整体优化。然而，考虑到基于该方案所设计的在线控制策略需要在当前时刻t₀来临时解决包含高维待优化变量的问题(14)，会消耗过多的计算资源且缺乏实际意义。因此，此处考虑一种离线控制策略(non-MPC)，假设在无人机控制策略设计之初可获得整个飞行持续时间内的信道状态信息，并通过内点法计算得到该理想状态下关于整个飞行过程的离线控制策略，即在线控制策略的理论性能上界，作为后续实验的基准性能指标。

图 5展示了采用了基准1，基准2和MPC方案的3种不同控制策略分别在不同波束成形方式下的四旋翼无人机的飞行轨迹，该飞行轨迹选自多组相互独立信道实现中的一组实验结果，其余实验结果均有相同的特性，不再展示。从图中可以看出，采用non-MPC会使得无人机尽可能靠近每个用户后抵达终点，提高用户与无人机之间的信道质量。其中，由于ZF non-MPC和RZF non-MPC策略具有足够长的视野范围且用户分布信息已知，在制定飞行策略时会选择在平稳飞行后主动靠近用户较多的位置，最终抵达终点。而MRT non-MPC则是主动靠近用户1后飞向终点。对于MPC方案，由于视野长度受限，会尽可能地优化其有限视野范围内的最优路径，因此会采取先靠近用户1后飞至用户数较多的位置，最后飞向终点。可以说，MPC方案是对离线控制策略的一种在线近似，能够根据当前时隙内系统状态和用户位置，自适应地调整系统输入以获得更高的系统收益。

表 2展示了分别采用了基准1、基准2和MPC方案的3种不同控制策略在不同波束成形方式下所获得的系统ASE。表中数值均为对50组相互独立随机信道实现下的结果取平均得到的。其中，基准1采用固定的飞行轨迹，从而损失了无人机控制系统的自由度，导致通信系统性能下降；而采用离线控制策略的基准2假设整个飞行路径中信道的状态信息已知，且视野范围足够长，因此获得了最大的系统ASE。相比之下，基于MPC的控制方案作为一种在线近似控制策略，获得了接近离线控制策略的系统ASE。

这表明基于MPC的控制方案能够在有限降低系统收益的情况下实现对四旋翼无人机的高效闭环控制，适用于更实际、多变的通信系统，并有效避免了模型系统不匹配的问题。因此，本文后续实验均基于MPC在线控制方案进行。

3.3 无人机通信系统模型仿真分析

后续将在3D MU-MISO场景下，对基于四旋翼无人机动态系统所提出的姿态角耦合通信模型与非耦合的原始通信模型进行性能对比，相关仿真结果均为对50组相互独立随机信道实现下的结果取平均得到的。其中，非耦合的原始通信模型实验忽略了姿态角对通信系统ASE的影响，仅针对与无人机三维坐标相关的状态进行优化，并将优化后的状态值代入实际通信模型中进行结果分析。图 6展示了天线的总发射功率与采用3种常见的波束成形方式在2种模型场景下所获得的系统ASE之间的关系。其中，设定天线阵列满足M_x=M_y=4，MPC的视野范围为N₀=25。实验结果表明，无论采取何种天线总发射功率大小和波束成形方式，本文所提出的通信模型均能够有效提升无人机通信系统的ASE。这是因为，该模型能够有效利用由四旋翼无人机的姿态角所带来的系统自由度增益，通过调整无人机的姿态使得其与地面用户之间所构成的信道的俯仰角和方位角最优，即使在天线总发射功率较小时也能获取更高的系统ASE。同时，随着总发射功率的增长，复杂系统所带来的ASE提升能够继续保持。此外，实验还表明采用ZF和RZF 2种波束成形方式获得的系统增益远高于MRT，这是因为前2种算法能够有效抑制用户之间的干扰，从而带来更高的系统ASE。而MRT算法复杂度较低，只能注重于最大化当前用户的ASE，受到其余用户的干扰较大。

图 7展示了天线数目与采用3种常见的波束成形方式在2种模型场景下所获得的系统ASE之间的关系。其中，设定天线总发射功率为P_max=100 W，MPC的视野范围为N₀=25。

与图 6类似，本文所提出的耦合复杂通信网络模型相较于原始的通信模型，无论采取何种天线阵列数目和波束成形方式，都能够有效提升无人机通信系统的ASE。这是因为，耦合复杂通信网络模型具有更高的系统自由度，随着发射天线数目的增加，复杂系统所带来的ASE提升能够继续保持。而与发射功率图不同的是，MRT方式则在增加天线数目时能够持续提升ASE，并且其增长速率逐渐减缓。其余2种波束成形方式也表现出类似的增长特性。

图 8展示了在不同波束成形方式下，采用2种不同模型的系统目标函数归一化后值的变化趋势。实线表示在姿态角耦合通信模型下，应用3种波束成形方式的目标函数归一化后值的收敛趋势。

可以清楚地看出，所提出的MPC方案在该场景下能够确保在MPC每一步迭代时目标函数值下降，并最后收敛至稳定状态。相比之下，原始通信模型下(MRT without，ZF without和RZF without)的性能表现较差。在阶段1时，MPC每一步迭代能够保证其目标函数值的下降，并逐步收敛至稳态。而在阶段2中，由于增加了关于抵达终点的相关约束，优化变量搜索空间被缩减，因此系统最优性受到了损失。该通信质量的损失无法通过调整无人机姿态得到弥补。值得一提的是，在阶段1到阶段2的转变过程中，6条曲线都存在1个跳变。这是由于在阶段2中引入终点等式约束，导致最优性受到损失，使得阶段2的第一次迭代所得到的状态变量与阶段1的最后一次迭代所得到的状态变量不同，前者的目标函数值相对较差。

4 总结

本文基于MPC控制算法，针对多天线四旋翼无人机同时为多用户服务的姿态控制及资源分配问题，提出具体的在线控制方案。该方案将整个控制过程分为2个阶段，能够根据当前时刻系统状态变量和控制变量的更新，通过联合优化无人机的三维坐标和姿态角，在当前时刻制定出视野范围以内的最优控制策略，使得其在满足四旋翼无人机动力学方程的条件下，为整个系统提供更高的ASE，重复迭代控制过程直至完成任务。仿真结果表明，所提出的基于MPC的控制策略能够指导无人机准确抵达终点。同时，验证了该在线控制策略能够在有限降低系统收益的情况下实现对离散控制方案的近似，从而适用于更加实际且多变的通信场景中。所提出的四旋翼无人机姿态角耦合通信模型能够带来系统自由度收益，并提高系统整体的ASE。后续可将本研究扩展至多无人机的复杂合作场景，并考虑控制系统的随机性，以制定无人机的鲁棒在线控制策略。