无人机系统是目前无人系统领域发展较快、水平较高、实际应用较多的一类空中无人系统, 具有成本低、无人员伤亡、操作方便和灵活可靠等特点。与固定翼飞机相比, 四旋翼飞机具有许多优点, 例如, 具备垂直起飞、着陆、悬停，以及室内飞行的能力。随着应用环境的日趋复杂以及任务的日益多样, 单机的能力受到了限制, 多机的编队协同已经成为无人机应用的重要发展趋势。

然而, 四旋翼无人机具有6个自由度, 却只有4个控制输入, 因此, 其动力学不仅是非线性的, 而且是耦合的、欠驱动的, 很难被控制。此外, 无人机在编队飞行的过程中, 会受到阵风、湍流和其他的外界干扰, 且考虑到四旋翼离散动力学的特征值位于单位圆的边界上, 因此, 需要一种先进的控制策略来实现自主编队飞行的稳定, 使编队具有较高的可操作性和对外界干扰的鲁棒性。为此, 人们做了很多努力, 并制定了一些策略来解决这类系统在编队飞行过程中遇到的问题。

近年来, 已经提出了几种策略来应对以上提出的一些问题。Guzey等^[1]开发了一种基于非线性输出反馈神经网络的一致性控制器。Yan等^[2]提出了基于一致性的双闭环滑模编队控制器。魏扬等^[3]考虑闭环系统存在时变干扰的情况下, 设计了用于无人机编队保持的自适应控制器。Zhen等^[4]解决了干扰和不确定性存在时, 编队飞行中无人机的速度和姿态协同控制问题。Ille等^[5]设计了分布式MPC用于跟踪给定的参考轨迹, 并用惩罚因子来避免编队内碰撞的发生。Chevet等^[6]提出基于MPC的编队重构算法, 不仅允许损坏的无人机离开编队, 还允许恢复的无人机重新加入当前编队, 同时避免碰撞。Cai等^[7]研究了一种用于多无人机编队控制的事件触发MPC方案, 可以减少分布式MPC方案的计算负担。Huang等^[8]提出了一种基于卡尔曼滤波和MPC的协同避碰方法, 指导多架无人机协同避碰决策。Wolfe等^[9]设计和测试了用于编队飞行跟踪的分布式多模型MPC控制器。Rosa等^[10]在非线性MPC的框架下, 设计了一种具有避障功能的编队飞行控制器。

MPC是解决约束问题最有效的控制方法之一, 它还可以通过多变量耦合实现多目标最优控制^[11]。但是, 无人机的离散时间模型是临界稳定的, 在MPC的滚动优化过程中, 这种临界稳定性可能会进一步恶化, 这是被许多研究所忽略的。本文针对编队中无人机的TS和RS分别设计多约束预测模型控制器, 以实现多约束和干扰存在的情况下, 编队的最优控制和对于规划路径较好的跟踪; 此外, 通过对预测控制的成本函数进行适当地修改, 不仅可以改善MPC计算中的数值问题, 而且可以确保最优计算过程中闭环系统的稳定性。

1 四旋翼动力学模型及编队算法 1.1 考虑干扰的四旋翼动力学模型

本文考虑的四旋翼结构是完全相同的。其动力学模型可以通过拉格朗日-欧拉的方法获得。这里, 考虑干扰因素的四旋翼非线性数学模型为^[12]

(1)

式中：(φ, θ, ψ)为四旋翼的3个姿态的欧拉角, 分别代表滚转角、俯仰角和偏航角; (x, y, z)为四旋翼的质心在惯性坐标系中的位置坐标; m为四旋翼无人机的总质量; g是重力加速度; U_i(i=1, 2, 3, 4)为对应的控制输入; Ω_r为螺旋桨角速度残差; d_i(i=1, 2, …, 6)表示外部扰动的集合, 且|d_i|≤D_i, D_i(i=1, 2, …, 6)为有界正实数。a_j(j=1, 2, 3, 4, 5), b_k(k=1, 2, 3), u_x, u_y和u_z的表达式如(2)~(3)所示

(2)

(3)

J_r为四旋翼的每个旋翼末端到飞行器重心的距离; (I_x, I_y, I_z)分别为围绕每个轴的转动惯量；l_a为转子的悬臂长度。

1.2 模型的离散化

由(1)式可以看出, 前3个方程描述了四旋翼的旋转动力学, 后面3个方程表示平移动力学, 因此可以将(1)式分成RS和TS。对于第i个无人机的RS, 在同时考虑干扰和不确定性的情况下, 线性时变离散状态空间模型为

(4)

式中, X_{i, R}(k)=[φ(k), , θ(k), , ψ(k), ]^T为RS的状态向量, U_{i, R}(k)=[U_{i, 2}, U_{i, 3}, U_{i, 4}, Ω_{i, r}]^T为控制输入向量, D_{i, R}(k)=[0, d_{i, R, 1}(k), 0, d_{i, R, 2}(k), 0, d_{i, R, 3}(k)]^T为未知干扰的集合。当RS的采样时间为T_R时, (4)式中的矩阵A_{i, R}, B_{i, R}, B_{D, R}和C_{i, R}分别为

(5)

(6)

(7)

(8)

同理, 对于第i个无人机的TS, 其离散线性化时变模型为

(9)

式中, X_{i, T}(k)=[x(k), , y(k), , z(k), ]^T表示TS的状态向量, U_{i, T}(k)=[U_{i, x}, U_{i, y}, U_{i, z}]=[u_xU_{i, 1}, u_yU_{i, 1}, mg-u_zU_{i, 1}]^T表示TS的控制输入向量, D_{i, T}(k)=[0, d_{i, T, 4}(k), 0, d_{i, T, 5}(k), 0, d_{i, T, 6}(k)]^T为未知干扰的集合。当TS的采样时间为T_T时, (6)式中矩阵A_{i, T}, B_{i, T}, B_{D, T}和C_{i, T}分别为

(10)

(11)

(12)

(13)

1.3 一致性编队算法

考虑n个相同的UAV的编队, 每个UAV由V_i(i∈{1, 2, …, n})表示。它们构成有向图G={V, E}, 其中V={V₁, V₂, …, V_n}是UAV节点的集合, E⊆V×V为图的边集, 其中图的边集由e_ij表示。表示V_i和V_j之间可以相互获得彼此的位置和速度信息。a_ij表示节点V_i和V_j之间连接的权重, a_ij=1表示节点V_i可以收到节点V_j的信息, 否则a_ij=0。V_i的邻点集合由N_i={V_j∈V: e_ij∈E}表示。

以TS的x轴方向为例

(14)

式中, 在只考虑定拓扑结构的系统时, (14)式的一致性控制律为^[12]

(15)

式中, p_ij为第i个无人机和第j个无人机之间的距离。k₁＞0。同理可以得到另外2个轴向的一致性控制律。

2 多约束MPC的设计

本文设计的主要目标是实现无人机编队的稳定和航迹跟踪, 而多无人机编队是一个多约束的复杂系统, 因此, 可以采用MPC的方法进行求解。在本节中, 分别针对每个UAV的TS和RS设计对应的MPC, 以保证UAV的稳定性和路径跟踪能力。如上一节所述, TS通过u_x, u_y和u_z决定了四旋翼的参考方向。因此, 首先对TS进行MPC的设计。

对于编队中第i个UAV, 以TS为例, 考虑干扰情况下的差分模型为

(16)

定义ΔX_{i, T}(k)=X_{i, T}(k)-X_{i, T}(k-1), ΔU_{i, T}(k)=U_{i, T}(k)-U_{i, T}(k-1), 从而, 可以得出一个新的增广系统为

(17)

式中：; , , , , ε(k)=D_{i, T}(k)-D_{i, T}(k-1)。这里, 假设ε(k)为白噪声的集合, 由于在MPC的计算中使用的是状态量的期望值, 可以认为ε(k)=0。用N_{c, T}表示未来控制量个数, N_{p, T}表示预测域长度, 则k时刻的控制量和输出的预测值分别表示为ΔU_{i, T}(k)=[ΔU_{i, T}(k), ΔU_{i, T}(k+1), …, ΔU_{i, T}(k+N_{c, T}-1)], =[Y_{i, T}(k+1|k), Y_{i, T}(k+2|k), …, Y_{i, T}(k+N_{p, T}|k)], 从而, 未来状态的预测方程为

(18)

其中, 矩阵F_{i, ζ}和H_{i, ζ}分别为

(19)

约束模型预测控制的本质是通过求解最优控制轨迹ΔU_{i, T}来最小化给定时间k的成本函数J_{i, T}(k)。当然, 这种最小化是基于J_{i, T}(k)受到一系列线性不等式约束下的实时优化, 通常, 这些约束可以表示为文献[13]中的形式。假如在k时刻, 取新的状态向量矩阵为X_{i, T}(k+i|k)=[ΔX_{i, T}(k+i|k); Y(k+i|k)-(C_{i, T}r_s(k+1|k))], 其中, r_s(k|k)=[x_r(k|k); (k|k); y_r(k|k); (k|k); z_r(k|k); (k|k)]为期望的状态, J_{i, T}(k)和线性不等式约束的定义如下

(20)

式中：；M_{i, T}的行数和列数分别与约束的个数以及ΔU_{i, T}(k)的元素个数相同；γ_{i, T}为转换后的约束矩阵^[13]。以上最优化存在的前提是满足Kuhn-Tucker条件。为求解上述优化问题, 得到(21)式所示的拉格朗日形式

(21)

式中：R_{s, T}(k|k)=[(C_{i, T}r_s(k+1|k)), (C_{i, T}r_s(k+2|k)), …, (C_{i, T}r_s(k+N_p|k))]为期望的输出矩阵，Y_{i, T}(k)=F_{i, T}X_{i, T}(k|k)+H_{i, T}ΔU_{i, T}(k)。Q_{i, T}和R_{i, T}是由矩阵Q_{i, T}和R_{i, T}组成的对角块阵。这里假设期望轨迹是缓慢变化的, 即Δr_s(k)=r_s(k+i+1|k)-r_s(k+i|k)≈0, (i=1, …, N_p), 将Y_{i, T}(k)代入(21)式中可以得到

(22)

在考虑ΔU_{i, T}和λ_{i, T}的情况下, 求解最优控制和拉格朗日乘子为：

(23)

(24)

在得到最优解后, 只有控制序列ΔU_{i, T}的第一行被用于控制向量中, 即下一刻的控制向量可以表示为: U_{i, T}(k+1)=U_{i, T}(k)+ΔU_{i, T}。假设参考的偏航角为零, 即ψ_{i, r}=0, 则第i个UAV的TS在k时刻的控制量U_{i, 1}可以通过(25)式获得

(25)

假设-π/2＜φ_{i, r}, θ_{i, r}＜π/2, 通过方程(3)可以推出第i个UAV的滚转角和俯仰角的参考值φ_{i, r}和θ_{i, r}分别为

(26)

(27)

最后, 对于第i个UAV, 在各个时刻传递给转动子系统的参考输入如下

(28)

同理, 可以得到RS的最优控制U₂, U₃, U₄和Ω_r, 进而应用到UAV的非线性动态模型中。

3 稳定性分析

(18) 式中存在积分项, 因此, 如果预测范围较大, 则(20)式的最优计算过程中会出现数值问题。从(5)式和(10)式可以看出, 系统矩阵的特征值位于单位圆边界附近, 这会使数值问题变得更糟^[13]。因此, 有必要设计适当的控制器, 在保证系统闭环稳定性的同时, 避免上述的数值问题。本文提出的算法, 是在MPC算法的基础上作适当的修改, 从而实现上述功能。

为了便于分析, 以第i架无人机的TS为例, 对于α_T＞1, 将MPC中的成本函数按以下方式处理

(29)

式中, 由于α_T＞1, 成本函数中的权重系数会随着预测时间的增加而减小, 因此, 在进行最优化的计算过程中, 相较于成本函数(20)式, 修改后的成本函数更加重视当前时间的状态和控制, 从而对应一个预测窗口, 新的状态序列和增量控制矢量选择为

(30)

(31)

矩阵和中的元素可以表示为

(32)

(33)

和增广状态空间模型(17)式类似, (k+m+1|k)和(k+m|k)满足

(34)

因此, 具有多个约束的优化问题可以重新表示为

(35)

式中, M_{T, α}的表达式为

(36)

从(21)式和(29)式可以看出, 这2个成本函数具有相同的形式, 此外, (29)式的最优解与(21)式也一样。因此, 如果选择适当的α_T, 则成本函数(29)式可以用于MPC的计算过程, 同时保证系统矩阵的所有特征值都在单位圆内, 对应状态空间(34)式是稳定的。当N_{P, T}和N_{C, T}足够大时, 优化问题(29)式等效于离散时间线性二次调节器(DLQR)问题^[13]。DLQR可以由代数Riccati方程按相同的权重矩阵Q_{i, T}和R_{i, T}来进行求解

(37)

式中对应的状态反馈增益矩阵K_{T, α}和闭环系统可以表示为

(38)

因此, 可以选择合适的α_T＞1, 来保证(38)式中的特征根。而实际系统的特征值将为, 这意味着(30)式的修改无法保证实际闭环系统的稳定性。但是, 可以通过选择适当的Q_{i, T}和R_{i, T}来解决这个问题。

定理1  令σ_T=α_T^-1β_T, (0＜β_T＜1), Q_σ=σ²Q_{i, T}+(1-σ_T²)S_∞, R_σ=σ_T²R_{i, T}, 通过选择一个合适的β_T, (39)式的最优解可以确保闭环系统(, )的稳定性。

(39)

证明  与(35)式类似, (39)式的代数Riccati方程为

(40)

将Q_σ和R_σ代入(40)式中可以得到

(41)

(41) 式两边同时乘以α_T²β_T², 可以得到

(42)

式中对应的状态反馈增益矩阵K_{T, σ}和闭环系统为

(43)

从(43)式闭环系统方程可以看出, 通过选择合适β_T, 能确保特征值, 也就意味着保证了闭环系统是稳定的。

综上所述, 通过(29)式和(39)式对原始成本函数(21)式进行连续变换, 等效于在最优化的求解过程中, 得到了(43)式所示的最终闭环函数。如前所述, 这样做的目的不仅保证最优计算过程中所使用模型的稳定性, 而且还能够保证闭环系统的稳定性。

4 仿真结果

以3架四旋翼无人机构成的编队为例, 这里提出一种围绕半径为r的圆构成的正三角形编队, 其几何结构如图 1所示。

其中, UAV₁能够接收所有的指令, 其主要任务是跟踪规划好的轨迹。无人机的主要参数见表 1。

编队中3架无人机初始状态分别为: X_{1, T}(0)=[3, 0, , 0, 0.1, 0];X_{2, T}(0)=[0, 0, 2, 0, 0.1, 0];X_{3, T}(0)=[0, 0, 0, 0, 0, 0], 单位: m; X_{1, R}(0)=X_{2, R}(0)=X_{3, R}(0)=[0, 0, 0, 0, 0, 0], 单位: rad。由于物理约束的存在, 编队中无人机所受到的约束有: 0≤U₁≤11.23, |U₂|≤5.61, |U₃|≤5.61, |U₄|≤0.16, 。仿真中, 分别在10, 15, 20, 25~30 s, 30 s及以后的时间施加干扰, 30 s之前的干扰值相同, 分别为d₁=d₂=d₃=3 rad/s², d₄=d₅=d₆=2.7 m/s²; 30 s及以后的干扰大小为d₁=d₂=d₃=0.67 rad/s², d₄=d₅=d₆=0.67 m/s²。假设x, y, z的期望轨迹为R₁(t), R₂(t), R₃(t), 分别为:

无人机编队的主要参数为: k_x=k_y=k_z=0.01。MPC在滚动优化过程中的主要参数为: , R_{i, T}=R_{i, R}=0.01I₆, I₆表示单位阵, N_{c, T}=N_{c, R}=1, N_{p, T}=10, N_{p, R}=4, α_T=1.6, β_T=0.9, α_R=1.2, β_R=0.92, r=2。得到仿真结果如下所示。

由图 2中2个子系统的闭环极点分布可以看出, 本文MMPC不仅能够保证所有闭环极点位于单位圆内, 相较于RMPC, 闭环极点能够更靠近原点, 也就保证了更多的稳定度。

通过图 3~7可以看出, 虽然干扰的存在对姿态角和三轴轨迹跟踪造成了一定程度的影响, 特别是对姿态角影响较为明显, 但由于控制器对于干扰存在鲁棒性, 使得编队的队形影响较小, 且能在很短的时间内恢复到正常状态。

图 8进一步证明了, 本文设计的控制器, 能够在很好处理约束的前提下, 实现编队稳定的轨迹跟踪和对干扰较好的鲁棒能力。

表 2则可以看出, MMPC的条件数比RMPC的要小很多, 证明了MMPC能够在实现上述控制效果的同时, 计算的复杂度能够大大降低。

5 结论

本文针对四旋翼的编队问题, 提出了一种改进的多约束模型预测控制器, 该控制器是基于误差的增广模型进行设计。主要结论如下:

1) 将四旋翼的三维空间模型分为TS和RS, 并分别针对2个子系统设计相应的多约束MPC。

2) 通过对成本函数进行合理的修改, 不仅保证了滚动优化过程中所使用模型的渐进稳定, 还保证了闭环系统的稳定性。

3) 利用增广模型的特点能够消除外部干扰对于无人机路径跟踪的影响。

4) 相较于RMPC, MMPC的计算复杂度能够大大降低, 并且能够保证在多约束存在的情况下, 实现编队较好的路径跟踪。仿真结果也证明了所述方法的有效性。

未来的工作会考虑不确定性、避障和编队重构等问题。