四旋翼无人机具有成本较低、设备简单、飞行灵活等特点，近些年来广泛应用于军事和民用领域，如目标侦察、应急救援、农业植保、无人机灯光表演等。随着任务复杂度的增加，单架无人机往往难以满足任务需求，因此无人机集群控制及其应用由此成为当下的一个研究热点，多无人机集群能够提高执行任务的效率、灵活性和容错能力。无人机编队队形控制方法是实现多无人机编队安全飞行的前提和必要手段^[1]，其中集群无人机队形重构是集群无人机所要考虑的一个重要问题，主要是一定数量的无人机根据特定要求变换队形，需要研究无人机路径规划算法和轨迹姿态控制算法等。

对于集群无人机队形重构的路径规划，主要目标是让每架无人机都能从初始位置无碰撞地到达最终位置，保证队形重构过程中代价最小或能耗最优^[2]。这其中既包含目标分配问题，即每架无人机被分配到的目标队形的某一特定位置；又包含轨迹生成问题，保证无碰撞地到达目标位置。许多算法对上述两个问题进行了解耦，即先进行目标分配，分配好目标位置之后再进行起始点到分配好的目标位置之间的轨迹生成，其中目标分配问题多利用匈牙利算法^[3]进行解决。分配好目标位置之后，在轨迹生成部分，常见算法包括图搜索算法A^*[4]及其变体^[5]、自适应随机搜索算法^[6]、模型预测控制^[7]、混合整数规划^[8-9]等。以上算法在单一无人机规划时效果较好，但是在多无人机轨迹规划时普遍存在计算难度加大、规划时间增长、规划效率难以满足实际应用时的实时性要求等缺点，因此，需要探索计算简便、效率高的多无人机路径规划算法。

在无人机集群由同构无人机构成的情况下，每个目标可以任意分配给各架无人机，则称之为无标号问题。若集群中包含异构无人机或某些特定情境下针对某一目标位置只能被分配给特定无人机(如快递运输)，则称之为带标号问题。Kumar等针对无标号问题提出了许多算法，其中CAPT(Concurrent Assignment and Planning of Trajectories)算法^[10]对多架无人机进行并行的目标分配和轨迹规划。考虑到将目标分配和轨迹生成问题解耦计算的问题复杂度，该算法将上述2个问题耦合，提出了包括集中式CAPT(C-CAPT)和分布式CAPT(D-CAPT)2种实现思路，实现了无障碍物情况下的多机无碰撞全局路径规划，极大降低了计算的复杂程度，提高了实时性和实用性。随后Kumar团队的Tang等基于带标号问题，在CAPT算法的基础上相继提出了OMP+CHOP(Optimal Motion Plans+ Circular HOlding Patterns)算法^[11]、HOOP(HOld or take Optimal Plan)算法^[12]，针对无障碍的拥挤环境，提出在无碰撞时采取CAPT算法、将要碰撞时采取圆形保持轨道进行轨迹设计的思路，扩展了CAPT算法的研究情景和研究领域。

对于无人机队形重构的控制问题，主要目标是控制无人机追踪预设轨迹，保证在一定外界干扰的情况下仍不发生偏移。若干扰有界，滑模控制算法在四旋翼无人机控制中应用广泛^[13]，在文献[14]中研究了一种终端滑模控制算法，可以很容易地消除飞行中扰动影响，但是主要缺点是控制律不连续引起的抖振。高阶滑模控制算法如螺旋和超螺旋滑模控制可以产生连续控制信号，抵消连续扰动，文献[15]利用内-外环结构设计了一种超螺旋滑模控制器。文献[16]则提出了利用连续律来控制无人机姿态的新策略，但前提条件必须是扰动边界已知。文献[17]提出了一种新的连续多变量积分滑模控制器。结合连续标称控制器和改进的超螺旋控制器来实现对干扰的有效抑制。基于扰动上界未知问题，引入自适应控制算法，文献[18]设计了一种自适应终端滑模控制策略，实现了有限时间收敛。文献[19]提出的自适应终端滑模有限时间控制算法能有效满足姿态跟踪控制的快速性需求。无人机控制算法大多针对单输入单输出系统进行设计，文献[20]设计多变量有限时间积分滑模控制器，保证了对期望参考输入的高精度快速跟踪控制。采用多变量控制算法可以避免对多输入多输出系统设计时进行解耦^[21]，从而一定程度上提高控制效果如减少滑模控制算法的抖振现象。

本文考虑集群无人机队形重构问题，首先基于四旋翼无人机轨迹姿态数学模型，利用CAPT算法进行集群无人机路径规划，实现对多无人机的并行位置分配和最优轨迹生成；其次，分别考虑干扰影响下的无人机轨迹及姿态控制系统，提出姿态解算及新型连续多变量积分滑模控制器设计算法，实现对无人机轨迹的有限时间精确跟踪；最后，搭建Gazebo-ROS无人机仿真平台，在该平台下实现了12架四旋翼无人机队形重构过程的整体效果仿真，验证所提算法的有效性和实用性。

1 四旋翼无人机轨迹姿态模型

本文所采用的四旋翼无人机模型示意图如图 1所示。x_B、y_B、z_B为无人机机体坐标系，x_W、y_W、z_W为无人机惯性坐标系。F_i(i=1, 2, 3, 4)及M_i(i=1, 2, 3, 4)分别代表 4个电机产生的升力及力矩。充分考虑四旋翼无人机的固有特性，无人机在飞行过程中的动力学因素及受力平衡关系，其轨迹姿态数学模型为^[22]

$ \mathit{\boldsymbol{\dot p}} = \mathit{\boldsymbol{v}} $	（1）

$ \mathit{\boldsymbol{\dot v}} = - g{\mathit{\boldsymbol{e}}_{\text{z}}} + {\tau _{\text{f}}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{e}}_{\text{z}}}/m + {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_1}\left( t \right) $	（2）

$ \mathit{\boldsymbol{ \boldsymbol{\dot \varTheta} }} = \mathit{\boldsymbol{W \boldsymbol{\varOmega} }} $	（3）

$ \mathit{\boldsymbol{I \boldsymbol{\dot \varOmega} }} = - \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \times \mathit{\boldsymbol{I \boldsymbol{\varOmega} }} + \mathit{\boldsymbol{\tau }} + {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_2}\left( t \right) $	（4）

式中：$\boldsymbol{p}=[x, y, z]^{\mathrm{T}} \in \mathbf{R}^{3}, \boldsymbol{v}=\left[v_{x}, v_{y}, v_{z}\right]^{\mathrm{T}} \in $ $ {{\mathbf{R}}^3}, \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = {[\varphi , \theta , \psi ]^{\rm{T}}} \in {{\mathbf{R}}^3}$分别为无人机惯性坐标系下的位置、速度和欧拉角；$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = {\left[ {{\omega _x}, {\omega _y}, {\omega _z}} \right]^{\rm{T}}} \in {{\mathbf{R}}^3}$表示机体坐标系下的姿态角速率；g=9.8 m/s²为地球重力加速度，e_z=[0, 0, 1]^T为常值向量；m为四旋翼无人机质量；τ_f为无人机的总升力；τ为无人机x, y, z这3个方向的控制力矩；$\mathit{\boldsymbol{I}} = {\mathop{\rm diag}\nolimits} \left( {{I_x}, {I_y}, {I_z}} \right) $为无人机惯性常值矩阵；Δ₁(t)、Δ₂(t)为综合干扰。包括模型参数不确定及外界干扰，旋转矩阵R和变换矩阵W定义为

$ \begin{gathered} \mathit{\boldsymbol{R}} = \hfill \\ \left[ {\begin{array}{*{20}{c}} {{\text{c}}\theta {\text{c}}\psi }&{{\text{s}}\theta {\text{c}}\psi {\text{s}}\varphi - {\text{s}}\psi {\text{c}}\varphi }&{{\text{s}}\theta {\text{c}}\psi {\text{c}}\varphi - {\text{s}}\psi {\text{s}}\varphi } \\ {{\text{c}}\theta {\text{s}}\psi }&{{\text{s}}\theta {\text{s}}\psi {\text{s}}\varphi - {\text{c}}\psi {\text{c}}\varphi }&{{\text{s}}\theta {\text{s}}\psi {\text{s}}\varphi - {\text{s}}\psi {\text{s}}\varphi } \\ { - {\text{s}}\theta }&{{\text{c}}\theta {\text{s}}\varphi }&{{\text{c}}\psi {\text{c}}\varphi } \end{array}} \right] \hfill \\ \end{gathered} $

（5）

$ \mathit{\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}} 1&{\sin \varphi \tan \theta }&{\cos \varphi \tan \theta } \\ 0&{\cos \varphi }&{ - \sin \varphi } \\ 0&{\sin \varphi \sec \theta }&{\cos \varphi \sec \theta } \end{array}} \right] $	（6）

式中：$\mathrm{c} \theta=\cos \theta, \mathrm{s} \theta=\sin \theta $，下同。控制量[τ_f，τ₁, τ₂, τ₃]^T与电机转速ω=[ω₁, ω₂, ω₃, ω₄]^T之间的关系为

$ \left[ {\begin{array}{*{20}{c}} {{\tau _{\text{f}}}} \\ {{\tau _1}} \\ {{\tau _2}} \\ {{\tau _3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{k_{\text{F}}}}&{{k_{\text{F}}}}&{{k_{\text{F}}}}&{{k_{\text{F}}}} \\ { - {k_{\text{c}}}}&{{k_{\text{c}}}}&{{k_{\text{c}}}}&{{k_{\text{c}}}} \\ {{k_{\text{c}}}}&{ - {k_{\text{c}}}}&{{k_c}}&{ - {k_c}} \\ {{k_{\text{M}}}}&{{k_{\text{M}}}}&{ - {k_{\text{M}}}}&{ - {k_{\text{M}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\omega _1^2} \\ {\omega _2^2} \\ {\omega _3^2} \\ {\omega _4^2} \end{array}} \right] $

（7）

式中：$k_{\mathrm{c}}=\sqrt{2} / 2 k_{\mathrm{F}} d ; d $为旋翼转动中心到无人机中心的距离；k_F、k_M分别为旋翼的升力系数和扭矩系数。

本文控制目标为：针对无人机轨迹姿态模型式(1)~式(4)，设计集群无人机路径规划算法和轨迹姿态跟踪控制器，实现集群无人机安全、迅速的队形重构。

2 集群无人机路径规划

本节将针对集群无人机路径规划问题，利用CAPT算法及无人机轨迹平滑方法，实现多无人机队形重构过程中的并行目标分配和轨迹生成。

2.1 CAPT算法

CAPT算法^[10]是Kumar团队提出的基于无标号问题的解决方案，可以解决集群由同构个体组成且目标位置可以任意分配的运动规划问题。本文研究的集群无人机队形重构问题，只要求形成特定的队形形状，而不关心每架无人机的具体位置，可以转化成无标号问题，因此可以采用CAPT算法求解。

考虑集群无人机队形重构问题，假设集群无人机规模与无人机目标位置数量相同为N，第i架无人机的初始位置为x_i∈Rⁿ，第j个目标位置g_j∈Rⁿ。解决无人机队形重构问题需要对无人机目标位置进行分配同时产生到达各个目标的安全轨迹。

定义系统状态向量$\mathit{\boldsymbol{X}}(t) = \left[ {{\mathit{\boldsymbol{x}}_1}{{(t)}^{\rm{T}}}, {\mathit{\boldsymbol{x}}_2}{{(t)}^{\rm{T}}}, \cdots } \right., \left. {{\mathit{\boldsymbol{x}}_N}{{(t)}^{\rm{T}}}} \right] $，系统目标状态向量$\mathit{\boldsymbol{G}} = \left[ {\mathit{\boldsymbol{g}}_1^{\rm{T}}, \mathit{\boldsymbol{g}}_2^{\rm{T}}, \cdots , \mathit{\boldsymbol{g}}_N^{\rm{T}}} \right] $。定义目标分配矩阵φ∈R^N×N，为每架无人机分配目标位置：

$ {\varphi _{i,j}} = \left\{ \begin{array}{l} 1\;\;\;\;如果无人机\;i\;被分配到目标\;j\\ 0\;\;\;如果无人机\;i\;没有被分配到目标\;j \end{array} \right. $	（8）

定义增广分配矩阵Φ≡φ$ \otimes $ I_n，$ \otimes $代表求取克罗克内积。

CAPT算法的目的在于规划轨迹p(t)→X(t)(t∈[0，T])，T为所有无人机到达最终位置的终止时间，将在后文进行说明。运动轨迹p(t)满足始端约束条件p(0)=X(0)和终端约束条件Φ^Tp(T)=G。

假设初始位置与目标位置、各无人机之间的距离满足

$ \left\{ \begin{array}{l} \left\| {{\mathit{\boldsymbol{x}}_i}\left( 0 \right) - {\mathit{\boldsymbol{x}}_j}\left( 0 \right)} \right\| > \Delta \;\;\;\forall i \ne j\\ \left\| {{\mathit{\boldsymbol{g}}_i} - {\mathit{\boldsymbol{g}}_j}} \right\| > \Delta \;\;\;\;\;\;\;\;\;\;\;\;\forall i \ne j \end{array} \right. $	（9）

则在满足$\mathit{\Delta}>2 \sqrt{2} R_{\mathrm{d}} $的条件下(其中R_d为圆盘形无人机半径，对于本文中应用的四旋翼模型，需满足假设条件$\mathit{\Delta} > 2\sqrt 2 d $)，各无人机之间满足避撞要求。

CAPT算法将无标号问题看做最优化问题进行求解，确定最优的分配方案和轨迹，最优化目标函数为

$ \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }},\mathit{\boldsymbol{p}}\left( t \right)} \sum\limits_{i = 1}^N {\int\limits_0^T {{{\mathit{\boldsymbol{\dot x}}}_i}{{\left( t \right)}^{\rm{T}}}{{\mathit{\boldsymbol{\dot x}}}_i}\left( t \right){\rm{d}}t} } $	（10）

问题的解由直线轨迹组成，且各无人机走过的距离的平方之和最小。即最优分配方案为

$ \left\{ \begin{array}{l} {\varphi ^ * } = \mathop {{\rm{argmin}}}\limits_\varphi \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {{\varphi _{i,j}}} } {D_{i,j}}\\ {D_{i,j}} = {\left\| {{\mathit{\boldsymbol{x}}_i}\left( 0 \right) - {\mathit{\boldsymbol{g}}_j}} \right\|^2} \end{array} \right. $	（11）

其中分配问题φ^*可以用匈牙利算法进行求解。

所有的机器人都是在一定时间之内到达目标位置，所有无人机到达最终位置的终止时间T为

$ T = \mathop {\max }\limits_i \frac{{{{\left\| {{\mathit{\boldsymbol{x}}_i}\left( 0 \right) - \sum\limits_{j = 1}^M {\varphi _{i,j}^ * {\mathit{\boldsymbol{g}}_j}} } \right\|}_2}}}{{{v_{\max }}}} $	（12）

定义关于时间的多项式函数$ \delta (t) \equiv \sum\limits_{i = 0}^k {{\alpha _i}} {t^i} \in \left[ {0, 1} \right]$，满足δ(0)=0，δ(T)=1。使式(10)中速度平方积分最小的轨迹是速度恒定且满足始终端约束条件的轨迹：

$ {\mathit{\boldsymbol{P}}^ * }\left( t \right) = \left( {1 - \delta \left( t \right)} \right)\mathit{\boldsymbol{X}}\left( 0 \right) + \delta \left( t \right)\mathit{\boldsymbol{ \boldsymbol{\varPhi} G}} $	（13）

2.2 无人机轨迹平滑方法

为了满足四旋翼无人机的动力学方程，需要对轨迹进行平滑，Mellinger和Kumar^[22]证明了最小化位置状态的四阶导可以实现对于轨迹的平滑追踪。代价函数为

$ \mathop {{\rm{min}}}\limits_{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }},\mathit{\boldsymbol{p}}\left( t \right)} \sum\limits_{i = 1}^N {\int\limits_0^T {{{\mathit{\boldsymbol{\ddddot x}}}_i}{{\left( t \right)}^{\text{T}}}{{\mathit{\boldsymbol{\ddddot x}}}_i}\left( t \right){\text{d}}t} } $	（14）

问题(14)的最优轨迹与式(10)所求直线轨迹相同，δ(t)为高阶时间多项式。

本文中为了实现对轨迹的精确平滑追踪，同时对位置、速度和加速度进行控制，由文献[23]中的最小化三次可微运动的方法，将δ(t)表示为5阶时间多项式函数。利用文献[23]中给出的不同终端约束条件下的最优状态轨迹参数表达式，在固定位置、固定速度和固定加速度的情况下，将位置状态求取p(t)→X(t)(t∈[0, T])进一步扩展成求取包括位置、速度、加速度在内的状态轨迹$ \mathit{\boldsymbol{s}}(k) = \left[ {{p_k}, {v_k}, {a_k}} \right]$，其中：p_k、v_k、α_k分别为X、Y、Z坐标系中单架无人机在某一时刻某一坐标维度的位置、速度和加速度。注意此时只针对X-Y-Z坐标系中的单架无人机轨迹进行求取，即求得的是单架无人机在某一时刻某一坐标维度的位置、速度和加速度状态信息。

状态轨迹参数为

$ \left[ {\begin{array}{*{20}{l}} {{\alpha _1}}\\ {{\alpha _2}}\\ {{\alpha _3}} \end{array}} \right] = \frac{1}{{{T^5}}}\left[ {\begin{array}{*{20}{c}} {720}&{ - 360T}&{60{T^2}}\\ { - 360T}&{168{T^2}}&{ - 24{T^3}}\\ {60{T^2}}&{ - 24{T^3}}&{3{T^4}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\Delta p}\\ {\Delta v}\\ {\Delta a} \end{array}} \right] $	（15）

定初始状态向量$\mathit{\boldsymbol{s}}(0) = \left[ {{p_0}, {v_0}, {a_0}} \right] $。终端状态向量$ \mathit{\boldsymbol{s}}(T) = \left[ {{p_{\rm{f}}}, {v_{\rm{f}}}, {a_{\rm{f}}}} \right]$，则

$ \left[ {\begin{array}{*{20}{c}} {\Delta p}\\ {\Delta v}\\ {\Delta a} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{p_f} - {p_0} - {v_0}T - \frac{1}{2}{a_0}{T^2}}\\ {{v_{\rm{f}}} - {v_0} - {a_0}T}\\ {{a_{\rm{f}}} - {a_0}} \end{array}} \right] $	（16）

最优状态轨迹为

$ \begin{array}{l} {\mathit{\boldsymbol{s}}^ * }\left( t \right) = \left[ {\begin{array}{*{20}{c}} {p\left( t \right)}\\ {v\left( t \right)}\\ {a\left( t \right)} \end{array}} \right] = \\ \;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {\frac{{{\alpha _1}}}{{120}}{t^5} + \frac{{{\alpha _2}}}{{24}}{t^4} + \frac{{{\alpha _3}}}{6}{t^3} + \frac{{{a_0}}}{2}{t^2} + {v_0}t + {p_0}}\\ {\frac{{{\alpha _1}}}{{24}}{t^4} + \frac{{{\alpha _2}}}{6}{t^3} + \frac{{{\alpha _3}}}{2}{t^2} + {a_0}t + {v_0}}\\ {\frac{{{\alpha _1}}}{6}{t^3} + \frac{{{\alpha _2}}}{2}{t^2} + {\alpha _3}t + {a_0}} \end{array}} \right] \end{array} $

（17）

在本文研究的集群无人机队形重构问题中，采用“悬停-运动-悬停”的运动模式，给定s(0)=[p₀, 0, 0]和s(T)=[p_f, 0, 0]，代入式(15)~式(17)中，得到最优状态轨迹：

$ {\mathit{\boldsymbol{s}}^ * }\left( t \right) = \left[ {\begin{array}{*{20}{c}} {p\left( t \right)}\\ {v\left( t \right)}\\ {a\left( t \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{{{\alpha _1}}}{{120}}{t^5} + \frac{{{\alpha _2}}}{{24}}{t^4} + \frac{{{\alpha _3}}}{6}{t^3} + {p_0}}\\ {\frac{{{\alpha _1}}}{{24}}{t^4} + \frac{{{\alpha _2}}}{6}{t^3} + \frac{{{\alpha _3}}}{2}{t^2}}\\ {\frac{{{\alpha _1}}}{6}{t^3} + \frac{{{\alpha _2}}}{2}{t^2} + {\alpha _3}t} \end{array}} \right] $

（18）

$ \left[ {\begin{array}{*{20}{l}} {{\alpha _1}}\\ {{\alpha _2}}\\ {{\alpha _3}} \end{array}} \right] = \frac{1}{{{T^5}}}\left[ {\begin{array}{*{20}{c}} {720}&{ - 360T}&{60{T^2}}\\ { - 360T}&{168{T^2}}&{ - 24{T^3}}\\ {60{T^2}}&{ - 24{T^3}}&{3{T^4}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{p_{\rm{f}}} - {p_0}}\\ 0\\ 0 \end{array}} \right] $	（19）

3 四旋翼无人机轨迹姿态控制

基于第2节路径规划所设计的无人机轨迹，本节将分别针对无人机位置环及姿态环，设计有限时间多变量滑模控制及姿态解算算法，实现干扰影响下的无人机位置-姿态跟踪控制。

无人机位置-姿态控制原理如图 2所示，设定在飞行过程中，期望位置p_ref=[x_ref, y_ref, z_ref]^T∈R³与期望偏航姿态角ψ_ref已知。首先考虑干扰影响下的无人机位置控制问题，设计位置环有限时间多变量积分滑模控制，实现外界干扰及模型不确定等综合干扰下的无人机位置跟踪控制；其次，通过姿态解算，得出期望俯仰姿态角θ_ref和滚转姿态角φ_ref，获得姿态内环的期望跟踪指令Θ_ref=[φ_ref，θ_ref，ψ_ref]^T；最后，考虑干扰影响下的无人机姿态控制问题，设计有限时间姿态控制器，实现综合干扰影响下的无人机姿态跟踪控制。图 2中p, v, Θ, Ω为无人机位置、速度、姿态角、姿态角速度实际值，acc_ref为期望加速度。

3.1 位置环控制器设计

定义位置误差e₁=p－p_ref，基于式(1)和式(2)，可得

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot e}}}_1} = {\mathit{\boldsymbol{e}}_2}\\ {{\mathit{\boldsymbol{\dot e}}}_2} = \underbrace {{\tau _{\rm{f}}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{e}}_{\rm{z}}}/m}_\mathit{\boldsymbol{u}} + \underbrace {\left( { - g{\mathit{\boldsymbol{e}}_{\rm{z}}} - {{\mathit{\boldsymbol{\ddot p}}}_{{\rm{ref}}}}} \right)}_{\mathit{\boldsymbol{f}}\left( t \right)} + \underbrace {{\Delta _1}\left( t \right)}_{\mathit{\boldsymbol{d}}\left( t \right)} = \\ \;\;\;\;\;\;\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{f}}\left( t \right) + \mathit{\boldsymbol{d}}\left( t \right) \end{array} \right. $

（20）

式中：$ \boldsymbol{e}_{2}=\boldsymbol{v}-\dot{\boldsymbol{p}}_{\mathrm{ref}} ; \boldsymbol{f}(t)$为已知项；d(t)为综合干扰。

假设1   综合干扰d(t)是导数有界且上界已知，即存在正实数L，使得$\|\dot{\boldsymbol{d}}(t)\| \leqslant L $。

本文所设计的控制器包括2部分：标称控制器u_n及补偿控制器u_d，即

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{u}}_{\rm{n}}} + {\mathit{\boldsymbol{u}}_{\rm{d}}} $	（21）

其中，标称控制器u_n用于保证无综合干扰(外加干扰及模型参数不确定综合)影响下系统的稳定跟踪控制；而补偿控制器u_d则用于保证对综合干扰的抑制。故位置环实际控制输入为

$ {\tau _{\rm{f}}} = m\mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{e}}_{\rm{z}}}} \right) $	（22）

下面将分别介绍标称控制器u_n和补偿控制器u_d的设计过程。

步骤1   标称控制器设计

针对标称系统：

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_1} = {\mathit{\boldsymbol{e}}_2}}\\ {{{\mathit{\boldsymbol{\dot e}}}_2} = {\mathit{\boldsymbol{u}}_{\rm{n}}} + \mathit{\boldsymbol{f}}\left( t \right)} \end{array}} \right. $	（23）

引理1^[16]   针对积分链系统：

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot x}}}_1} = {\mathit{\boldsymbol{x}}_2}}\\ {{{\mathit{\boldsymbol{\dot x}}}_2} = {\mathit{\boldsymbol{x}}_3}}\\ {\; \vdots }\\ {{{\mathit{\boldsymbol{\dot x}}}_n} = \mathit{\boldsymbol{\tau }}} \end{array}} \right. $	（24）

定义k₁, k₂, …, k_n>0使$ {s^n} + {k_n}{s^{n - 1}} + \cdots + {k_2}s + {k_1}$是Hurwitz的，且存在ε∈(0, 1)使得对于每个r∈(1－ε, 1)，如果采用如下反馈控制器：

$ \mathit{\boldsymbol{\tau }} = - {k_1}{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{{r_1}}}\frac{{{\mathit{\boldsymbol{x}}_1}}}{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}} - \cdots - {k_n}{\left\| {{\mathit{\boldsymbol{x}}_n}} \right\|^{{r_n}}}\frac{{{\mathit{\boldsymbol{x}}_n}}}{{\left\| {{\mathit{\boldsymbol{x}}_n}} \right\|}} $	（25）

式中:r₁, r₂, …, r_n满足

$ {r_{i - 1}} = \frac{{{r_i}{r_{i + 1}}}}{{2{r_{i + 1}} - {r_i}}}\;\;\;\;i = 2,3, \cdots ,n $	（26）

且r_n=r, r_n+1=1，那么系统式(24)是有限时间收敛到平衡点的。

引理2^[24]   针对系统$ \dot{x}=f(t, x)$，如果存在连续可微的正定函数V(x)：D→Rⁿ，正常数λ及0 < τ < 1，使得$ \dot{V}(x)+\lambda V^{\tau}(x) \leqslant 0$成立，则系统是有限时间稳定的，收敛时间为

$ {T_{{\rm{reach}}}} \le \frac{{{V^{1 - \tau }}\left( {{x_0}} \right)}}{{\lambda \left( {1 - \tau } \right)}} $	（27）

式中:V(x₀)为V(x)的初值；T_reach为收敛时间。

考虑标称系统式(23)，设计标称控制器为

$ {\mathit{\boldsymbol{u}}_{\rm{n}}} = - {k_1}{\left\| {{\mathit{\boldsymbol{e}}_1}} \right\|^{{r_1}}}\frac{{{\mathit{\boldsymbol{e}}_1}}}{{\left\| {{\mathit{\boldsymbol{e}}_1}} \right\|}} - {k_2}{\left\| {{\mathit{\boldsymbol{e}}_2}} \right\|^{{r_2}}}\frac{{{\mathit{\boldsymbol{e}}_2}}}{{\left\| {{\mathit{\boldsymbol{e}}_2}} \right\|}} - \mathit{\boldsymbol{f}}\left( t \right) $	（28）

式中：k₁, k₂>0。

将式(28)代入标称系统(23)可得

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot e}}}_1} = {\mathit{\boldsymbol{e}}_2}\\ {{\mathit{\boldsymbol{\dot e}}}_2} = - {k_1}{\left\| {{\mathit{\boldsymbol{e}}_1}} \right\|^{{r_1}}}\frac{{{\mathit{\boldsymbol{e}}_1}}}{{\left\| {{\mathit{\boldsymbol{e}}_1}} \right\|}} - {k_2}{\left\| {{\mathit{\boldsymbol{e}}_2}} \right\|^{{r_2}}}\frac{{{\mathit{\boldsymbol{e}}_2}}}{{\left\| {{\mathit{\boldsymbol{e}}_2}} \right\|}} \end{array} \right. $

（29）

故由引理1及式(29)可知，标称控制器式(28)能够保证系统式(23)的有限时间收敛。然而系统式(23)未考虑外界干扰和系统内部不确定等综合干扰的影响，因此，下面将给出连续补偿控制器设计过程，用于对系统的综合干扰进行处理。

步骤2   补偿控制器设计

设计积分滑模面：

$ \mathit{\boldsymbol{s}} = {\mathit{\boldsymbol{e}}_2}\left( t \right) - {\mathit{\boldsymbol{e}}_2}\left( 0 \right) - \int_0^t {\left( {{\mathit{\boldsymbol{u}}_{\rm{n}}} + \mathit{\boldsymbol{f}}\left( t \right)} \right){\rm{d}}t} $	（30）

式中：e₂(0)为e₂(t)的初值。由式(30)可以看出，当t=0时，滑模面s=0，即系统初始状态就位于滑模面上，避免了趋近滑模的过程，能够保证整个过程的鲁棒性。

对积分滑模面(30)求导，可得

$ \mathit{\boldsymbol{\dot s}} = \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{d}}\left( t \right) - {\mathit{\boldsymbol{u}}_{\rm{n}}} $	（31）

由式(31)可知，当$ \mathit{\boldsymbol{\dot s}} = {\mathbf{0}}$时，可得等效控制为

$ {\mathit{\boldsymbol{u}}_{{\rm{eq}}}} = {\mathit{\boldsymbol{u}}_{\rm{n}}} - \mathit{\boldsymbol{d}}\left( t \right) $	（32）

代入系统(20)可得

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot e}}}_1} = {\mathit{\boldsymbol{e}}_2}}\\ {{{\mathit{\boldsymbol{\dot e}}}_2} = {\mathit{\boldsymbol{u}}_{\rm{n}}} + \mathit{\boldsymbol{f}}\left( t \right)} \end{array}} \right. $	（33）

通过比较式(33)及式(23)可得，当系统状态到达滑模面s后，综合干扰影响下的系统即转化为标称控制系统，基于所设计的标称控制器(28)即可实现系统的有限时间稳定。因此，此时的控制目标转变为：设计补偿控制器u_d，保证系统状态有限时间内到达滑模面s。

对式(31)再次求导，可得

$ \mathit{\boldsymbol{\ddot s}} = \mathit{\boldsymbol{\tilde \omega }} + \mathit{\boldsymbol{D}} $	（34）

式中：$\tilde{\boldsymbol{\omega}}=\dot{\boldsymbol{u}}_{\mathrm{d}} ; \boldsymbol{D}=\dot{\boldsymbol{d}}(t) $。

根据式(31)和式(34)，设计虚拟控制$ \mathit{\boldsymbol{\widetilde \omega }}$为

$ \mathit{\boldsymbol{\tilde \omega }} = - \alpha \left( {\frac{\mathit{\boldsymbol{s}}}{{\left\| \mathit{\boldsymbol{s}} \right\|}} + 0.5\frac{{\mathit{\boldsymbol{\dot s}}}}{{\left\| {\mathit{\boldsymbol{\dot s}}} \right\|}}} \right) $	（35）

式中：α>2L为正常数。

而真实控制器为

$ {\mathit{\boldsymbol{u}}_{\rm{d}}} = \int {\mathit{\boldsymbol{\tilde \omega }}{\rm{d}}\tau } $	（36）

由式(35)和式(36)可以看出，不连续量出现在虚拟控制输入$ \tilde{\boldsymbol{\omega}}=\dot{\boldsymbol{u}}_{\mathrm{d}} ; \boldsymbol{D}=\dot{\boldsymbol{d}}(t)$中，而真正的补偿控制输入u_d是连续的，因此能够有效减少抖振，提高控制精度。

定理1  考虑综合干扰影响下的跟踪控制系统式(20)，在假设1成立的条件下，如果控制器设计为式(21)，其中标称控制器u_n设计如式(28)所示，补偿控制器u_d设计如式(36)所示，那么存在一系列常数0 < r₁ < 1，0 < r₂ < 1，k₁>0, k₂>0，α>2L，使得系统状态跟踪误差是有限时间收敛到零的。

证明   为了后文书写方便，定义变量x₁=s, ${\mathit{\boldsymbol{x}}_{\rm{2}}}\mathit{\boldsymbol{ = \dot s}} $，故式(34)转化为

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot x}}}_1} = {\mathit{\boldsymbol{x}}_2}}\\ {{{\mathit{\boldsymbol{\dot x}}}_2} = \mathit{\boldsymbol{\tilde \omega }} + \mathit{\boldsymbol{D}}} \end{array}} \right. $	（37）

定义以下Lyapunov函数：

$ \begin{array}{*{20}{c}} {V = {\alpha ^2}\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{x}}_1} + \gamma \mathit{\boldsymbol{x}}_1^{\rm{T}}{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}}{\mathit{\boldsymbol{x}}_2} + }\\ {\alpha \left\| {{\mathit{\boldsymbol{x}}_1}} \right\|\mathit{\boldsymbol{x}}_2^{\rm{T}}{\mathit{\boldsymbol{x}}_2} + \frac{1}{4}{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^4}} \end{array} $

（38）

式中：γ为待设计正常数。式(38)可重新书写为

$ V = \left\| {{\mathit{\boldsymbol{x}}_1}} \right\|{\mathit{\boldsymbol{z}}^{\rm{T}}}\mathit{\boldsymbol{Az}} + \frac{1}{4}{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^4} $	（39）

式中：$\mathit{\boldsymbol{z}} = {\left[ {\frac{{{\mathit{\boldsymbol{x}}_1}}}{{\parallel {\mathit{\boldsymbol{x}}_1}{\parallel ^{1/2}}}}, {\mathit{\boldsymbol{x}}_2}} \right]^{\rm{T}}} $，矩阵A计算为

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {{\alpha ^2}}&{\gamma /2}\\ {\gamma /2}&\alpha \end{array}} \right] $	（40）

为保证矩阵A的正定，参数γ需设计为$ |\gamma| <2 \alpha^{3 / 2}$。考虑到$\lambda_{\min }(\boldsymbol{A})\|\mathit{\boldsymbol{z}}\|^{2} \leqslant \mathit{\boldsymbol{z}}^{\mathrm{T}} \boldsymbol{A} \mathit{\boldsymbol{z}} \leqslant \lambda_{\max }(\boldsymbol{A})\|\boldsymbol{z}\|^{2} $，考虑到$ \left\|\boldsymbol{x}_{2}\right\|^{4}$，由式(39)可得

$ \begin{array}{l} V \le {\lambda _{\max }}\left( \mathit{\boldsymbol{A}} \right)\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^2} + \left\| {{\mathit{\boldsymbol{x}}_1}} \right\|{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right) + \frac{1}{4}{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^4} \le \\ \;\;\;{\lambda _{\max }}\left( \mathit{\boldsymbol{A}} \right)\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^2} + \frac{{{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^4}}}{2}} \right) + \frac{1}{4}{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^4} = \\ \;\;\;\frac{3}{2}{\lambda _{\max }}\left( \mathit{\boldsymbol{A}} \right){\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^2} + \frac{1}{2}\left( {{\lambda _{\max }}\left( \mathit{\boldsymbol{A}} \right) + \frac{1}{2}} \right){\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^4} = \\ \;\;{\mathit{\boldsymbol{\theta }}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\theta }} \end{array} $

（41）

式中：$\mathit{\boldsymbol{\theta = }}{\left[ {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|, {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right]^{\rm{T}}} $，矩阵P₁为

$ {\mathit{\boldsymbol{P}}_1} = \left[ {\begin{array}{*{20}{c}} {\frac{3}{2}{\lambda _{\max }}\left( \mathit{\boldsymbol{{\rm A}}} \right)}&0\\ 0&{\frac{1}{2}\left( {{\lambda _{\max }}\left( \mathit{\boldsymbol{{\rm A}}} \right) + \frac{1}{2}} \right)} \end{array}} \right] $	（42）

由式(42)可知，如果矩阵A>0那么矩阵P₁>0。由于

$ {\lambda _{min}}\left( {{\mathit{\boldsymbol{P}}_1}} \right){\left\| \mathit{\boldsymbol{\theta }} \right\|^2} \le {\mathit{\boldsymbol{\theta }}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\theta }} \le {\lambda _{max}}\left( {{\mathit{\boldsymbol{P}}_1}} \right){\left\| \mathit{\boldsymbol{\theta }} \right\|^2} $	（43）

因此，式(41)可以写为

$ \begin{array}{l} {\lambda _{min}}\left( \mathit{\boldsymbol{A}} \right)\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^2} + \left\| {{\mathit{\boldsymbol{x}}_1}} \right\|{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right) + \frac{1}{4}{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^4} \le \\ \;\;\;\;\;\;\;\;V \le {\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_1}} \right){\left\| \mathit{\boldsymbol{\theta }} \right\|^2} \le \\ \;\;\;\;\;\;\;\;{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_1}} \right){\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}} + \left\| {{\mathit{\boldsymbol{x}}_2}} \right\|} \right)^4} \end{array} $

（44）

由此可见，Lyapunov函数V是正定的。对V求导，可得

$ \begin{array}{l} \dot V = 2{\alpha ^2}\mathit{\boldsymbol{x}}_1^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_1} + \frac{3}{2}\gamma {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{1/2}}\mathit{\boldsymbol{x}}_2^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_1} + \gamma {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{1/2}}\mathit{\boldsymbol{x}}_1^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_2} + \\ \;\;\;\;\;\;\mathit{\boldsymbol{x}}_2^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_2}{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^2} + \alpha \left( {\frac{{\mathit{\boldsymbol{x}}_1^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_1}}}{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}}\mathit{\boldsymbol{x}}_2^{\rm{T}}{\mathit{\boldsymbol{x}}_2} + 2\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|\mathit{\boldsymbol{x}}_2^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_2}} \right) = \\ \;\;\;\;\;\;\left( {2{\alpha ^2}\mathit{\boldsymbol{x}}_1^{\rm{T}} + \frac{3}{2}\gamma {{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}}\mathit{\boldsymbol{x}}_2^T + \alpha \frac{{\mathit{\boldsymbol{x}}_1^{\rm{T}}{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}}}{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}}} \right){{\mathit{\boldsymbol{\dot x}}}_1} + \\ \;\;\;\;\;\;\left( {\gamma {{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}}\mathit{\boldsymbol{x}}_1^{\rm{T}} + 2\alpha \left\| {{\mathit{\boldsymbol{x}}_1}} \right\|\mathit{\boldsymbol{x}}_2^{\rm{T}} + \mathit{\boldsymbol{x}}_2^{\rm{T}}{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right){{\mathit{\boldsymbol{\dot x}}}_2} = \\ \;\;\;\;\;\; - \gamma \left( {\alpha + 0.5\alpha \frac{{\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{x}}_2}}}{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}} - \frac{{\mathit{\boldsymbol{x}}_1^{\rm{T}}\mathit{\boldsymbol{D}}}}{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}}} \right){\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^{3/2}} - \\ \;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|\left[ {\alpha \left( {\alpha - \frac{{\mathit{\boldsymbol{x}}_2^{\rm{T}}\mathit{\boldsymbol{D}}}}{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}}} \right)\left\| {{\mathit{\boldsymbol{x}}_1}} \right\| - \frac{3}{2}\gamma {{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}}\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|} \right] + \\ \;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|\left[ {\left( {\frac{1}{2}\alpha - \frac{{\mathit{\boldsymbol{x}}_2^{\rm{T}}\mathit{\boldsymbol{D}}}}{{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}}} \right){{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right] \end{array} $

（45）

由于α>2L，且||D||≤L，式(45)可以写为

$ \begin{array}{l} \dot V \le - \left\| {{\mathit{\boldsymbol{x}}_2}} \right\|\left[ {\alpha \left( {\alpha - 2L} \right)\left\| {{\mathit{\boldsymbol{x}}_1}} \right\| - \frac{3}{2}\gamma {{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}}\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|} \right] - \\ \;\;\;\;\;\left( {\frac{1}{2}\alpha - L} \right){\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^3} - \gamma \left( {\frac{1}{2}\alpha - L} \right){\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{3/2}} = \\ \;\;\;\;\; - \gamma \left( {1/2\alpha - L} \right){\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{3/2}} - \left\| {{\mathit{\boldsymbol{x}}_2}} \right\|{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PB}} \end{array} $

（46）

式中：$\boldsymbol{B}=\left[\begin{array}{l} \left.\left\|\boldsymbol{x}_{1}\right\|^{1 / 2}, \left\|\boldsymbol{x}_{2}\right\|\right]^{\mathrm{T}} \end{array}\right. $，矩阵P为

$ \mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} {2\alpha \left( {\frac{1}{2}\alpha - L} \right)}&{ - \frac{3}{4}\gamma }\\ { - \frac{3}{4}\gamma }&{\frac{1}{2}\alpha - L} \end{array}} \right] $	（47）

如果参数α, γ满足下列条件：

$ \left\{ \begin{array}{l} \alpha > 2L\\ 0 < \gamma < \frac{{4\sqrt 2 }}{3}\sqrt \alpha \left( {\frac{1}{2}\alpha - L} \right) \end{array} \right. $	（48）

那么，矩阵P是正定的，且由式(46)可以看出$ {\dot V}$是负定的。由于矩阵P是正定的，那么以下不等式成立

$ \begin{array}{l} {\lambda _{\min }}\left( \mathit{\boldsymbol{P}} \right)\left( {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\| + {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right) \le \\ \;\;\;\;\;\;\;{\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{BP}} \le {\lambda _{\max }}\left( \mathit{\boldsymbol{P}} \right)\left( {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\| + {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right) \end{array} $

（49）

基于式(44)和式(49)，式(46)可以写为

$ \begin{array}{l} \dot V \le - \gamma \left( {\frac{1}{2}\alpha - L} \right){\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{3/2}} - \\ \;\;\;\;\;\;\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|{\lambda _{\min }}\left( \mathit{\boldsymbol{P}} \right)\left( {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\| + {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^2}} \right) \le \\ \;\;\;\;\;\; - \gamma {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^{3/2}}\left( {\frac{1}{2}\alpha - L} \right) - {\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^3}{\lambda _{\min }}\left( \mathit{\boldsymbol{P}} \right) \le \\ \;\;\;\;\;\; - K\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{3/2}} + {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^3}} \right) \le \\ \;\;\;\;\;\; - \frac{K}{{{2^{2/3}}}}\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^{1/2}} + {{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|}^3}} \right) \le \\ \;\;\;\;\;\; - \frac{K}{{{2^{2/3}}{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_1}} \right)}}{V^{3/4}} \end{array} $

（50）

式中：$ K=\min \left\{\lambda_{\min }(\boldsymbol{P}), \gamma\left(\frac{1}{2} \alpha-L\right)\right\}$。基于式(50)及引理2可得，系统状态x₁、x₂会在有限时间收敛到零，时间为$ \frac{4 V^{1 / 4}(0)}{\varepsilon}$，其中$\varepsilon = \frac{K}{{{2^{2/3}}\lambda _{\max }^{1/4}\left( {{\mathit{\boldsymbol{P}}_1}} \right)}}, V(0) $为V的初始值。基于定义${\mathit{\boldsymbol{x}}_1} = \mathit{\boldsymbol{s}}, {\mathit{\boldsymbol{x}}_2}\mathit{\boldsymbol{ = \dot s}} $，即可得有限时间内s=0。基于式(31)~式(33)的描述，当到达滑模面s后，系统控制等同于标称系统控制。故基于引理1可得，系统跟踪误差是有限时间收敛到零的。定理1得证。

3.2 姿态解算算法

四旋翼无人机无法通过控制力矩直接对飞行轨迹进行控制，需要对式(21)中位置控制器得到的控制指令进行姿态解算，转化为期望的飞行姿态。本文中采用悬停控制思想，追踪期望位置p_ref和期望偏航角ψ_ref=ψ₀，姿态解算算法^[25]为

$ \left\{ \begin{array}{l} {\varphi _{{\rm{ref}}}} = arcsin\frac{{{\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_{\rm{1}}}}}sin{\psi _{{\rm{ref}}}} - {\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_{\rm{2}}}}}cos{\psi _{{\rm{ref}}}}}}{{{\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_1}}} + {\rm{ac}}{{\rm{c}}_{{\rm{ref}}_2^2}} + {{\left( {{\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_3}}} + g} \right)}^2}}}\\ {\theta _{{\rm{ref}}}} = \arctan \frac{{{\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_{\rm{1}}}}}sin{\psi _{{\rm{ref}}}} - {\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_{\rm{2}}}}}\sin {\psi _{{\rm{ref}}}}}}{{{\rm{ac}}{{\rm{c}}_{{\rm{re}}{{\rm{f}}_3}}} + g}} \end{array} \right. $

（51）

式中：φ_ref、θ_ref、ψ_ref分别为期望的滚转角、俯仰角和偏航角。且期望加速度acc_ref=－ge_z+u+Δ₁(t)，acc_refj(j=1, 2, 3)表示向量acc_ref的第j个元素。

3.3 姿态控制器设计

基于式(3)和式(4)，定义姿态误差$ {\mathit{\boldsymbol{e}}_3} = \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} - {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_{{\rm{ref}}}}, {\mathit{\boldsymbol{e}}_4} = \mathit{\boldsymbol{W \boldsymbol{\varOmega} }} - {{\mathit{\boldsymbol{\dot \Theta}}}_{{\rm{ref}}}}$，则可得

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot e}}}_3} = {\mathit{\boldsymbol{e}}_4}\\ {{\mathit{\boldsymbol{\dot e}}}_4} = \underbrace {\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{I}}^{ - 1}}\mathit{\boldsymbol{\tau }}}_{{\mathit{\boldsymbol{u}}_{\rm{a}}}} + \underbrace { - \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{I}}^{ - 1}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \times \mathit{\boldsymbol{I \boldsymbol{\varOmega} }} - {{\mathit{\boldsymbol{ \boldsymbol{\ddot \varTheta} }}}_{{\rm{ref}}}}}_{{\mathit{\boldsymbol{f}}_{\rm{a}}}\left( t \right)} + \\ \;\;\;\;\;\;\underbrace {\mathit{\boldsymbol{\dot W \boldsymbol{\varOmega} }} \times \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{I}}^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_2}}_{{\mathit{\boldsymbol{d}}_{\rm{a}}}\left( t \right)} = \\ \;\;\;\;\;\;{\mathit{\boldsymbol{u}}_{\rm{a}}} + {\mathit{\boldsymbol{f}}_{\rm{a}}}\left( t \right) + {\mathit{\boldsymbol{d}}_{\rm{a}}}\left( t \right) \end{array} \right. $

（52）

等同于位置环控制器设计，可进行姿态跟踪控制器u_a的设计。姿态控制器包括2部分：

$ {\mathit{\boldsymbol{u}}_{\rm{a}}} = {\mathit{\boldsymbol{u}}_{{\rm{na}}}} + {\mathit{\boldsymbol{u}}_{{\rm{da}}}} $

式中:标称控制器u_na用于保证无综合干扰(外加干扰及模型参数不确定综合)影响下系统的稳定跟踪控制；而补偿控制器u_da则用于保证对综合干扰的抑制。标称控制器及补偿控制器分别设计为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_{{\rm{na}}}} = - {k_3}{{\left\| {{\mathit{\boldsymbol{e}}_3}} \right\|}^{{r_3}}}\frac{{{\mathit{\boldsymbol{e}}_3}}}{{\left\| {{\mathit{\boldsymbol{e}}_3}} \right\|}} - {k_4}{{\left\| {{\mathit{\boldsymbol{e}}_4}} \right\|}^{{r_4}}}\frac{{{e_4}}}{{\left\| {{\mathit{\boldsymbol{e}}_4}} \right\|}} - {\mathit{\boldsymbol{f}}_{\rm{a}}}\left( t \right)}\\ {{\mathit{\boldsymbol{u}}_{{\rm{da}}}} = \int {{{\mathit{\boldsymbol{\tilde \omega }}}_{\rm{a}}}{\rm{d}}\tau } } \end{array}} \right. $

（53）

式中：虚拟控制量$ {\mathit{\boldsymbol{\tilde{\omega}}}}_{\mathrm{a}}=-\beta\left(\frac{ {\mathit{\boldsymbol{s}}}_{\mathrm{a}}}{\left\| {\mathit{\boldsymbol{s}}}_{\mathrm{a}}\right\|}+0.5 \frac{{\mathit{\boldsymbol{\dot{s}}}}_{\mathrm{a}}}{\left\|{\mathit{\boldsymbol{\dot{s}}}}_{\mathrm{a}}\right\|}\right)$，${\mathit{\boldsymbol{s}}_{\rm{a}}} = {\mathit{\boldsymbol{e}}_4} - {\mathit{\boldsymbol{e}}_4}(0) - \int_0^t {\left( {{\mathit{\boldsymbol{u}}_{{\rm{na}}}} + {\mathit{\boldsymbol{f}}_{\rm{a}}}(t)} \right)} {\rm{d}}t $。证明过程与位置环控制相同，此处不再赘述。

最终得到姿态环控制力矩为

$ \mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{W}}^{ - 1}}{\mathit{\boldsymbol{u}}_{\rm{a}}} $	（54）

3.4 仿真校验

为验证设计的轨迹姿态控制算法的有效性，在MATLAB环境下对单架无人机进行仿真，仿真采样频率设置为0.001，仿真过程中各参数选取如下：

1) 四旋翼无人机物理参数：质量m=0.5 kg，惯性参数I_x=2.3×10^－3 kg·m²，I_y=2.4×10^－3 kg·m²，I_z=2.6×10^－3 kg·m²，旋翼转动中心到无人机中心的距离d=0.225 m，旋翼的升力系数k_F=5.55×10^－8 N·m· s²/rad²，旋翼的扭矩系数k_M=2.5×10^－9 N·s²/rad²。

2) 期望轨迹p_ref=[sin t, cos t, sin t]，期望偏航角ψ_ref=0 °，所受干扰为Δ₁=Δ₂=[0.1sin t，0.1cos t，0.1(sin t+cos t)]，初始值p₀=[2, －3, 2]^T, v₀=[0, 0, 0]^T，Θ₀=[0.2, －0.3, 0.2]^T，Ω₀=[0, 0, 0]^T。

3) 控制器参数：位置控制器k₁=5, k₂=4, r₁=1/2, r₂=5/7, α=5。姿态控制器k₃=12, k₄=5, r₃=1/2, r₄=5/7, β=4。

仿真结果如图 3~图 5所示。由图 3和图 4可知，在综合干扰影响的情况下，位置跟踪误差及姿态跟踪误差均可以在有限时间内收敛到零。

总控制量[τ_f，τ₁，τ₂，τ₃]^T的仿真结果如图 5所示，可以看到控制输入连续且稳定，系统在一定外界干扰的情况下仍能保持稳定，且取得良好的轨迹姿态跟踪效果。由图 3~图 5的仿真结果可以有效验证所提控制算法的有效性。

4 仿真平台搭建与效果验证 4.1 仿真平台基本结构

集群无人机队形重构即一定数量的无人机根据特定任务需求，从一个队形变换形成另一个新的队形形状。需要考虑的问题涉及目标位置分配、路径规划、避撞要求、姿态轨迹跟踪控制等多个方面。要对上述问题分别考虑控制策略，并以一种直观可靠的方式对控制算法进行验证和效果对比，就需要进行综合虚拟仿真实验。

本文采用的Gazebo仿真软件，是基于ROS操作系统的一款三维物理仿真软件，与ROS操作系统配套使用，可实现完整的“建模-仿真-可视化”的仿真流程，与MATLAB仿真相比，Gazebo仿真更贴近真实物理情景，仿真的实时性和真实性更强，同时能够克服MATLAB仿真曲线表示缺乏直观性的缺点，实现仿真结果的三维立体视景演示。Gazebo软件中可以搭建近似于现实世界的仿真场景，构建多种传感器仿真模型，搭建包括无人机在内的各种机器人模型，包含多个物理仿真引擎，由此可以为机器人模型添加现实世界的物理性质，使仿真对象满足贴近真实的运动学和动力学物理要求。

仿真平台的基本框架如图 6所示。其中ROS作为机器人开源的元操作系统，主要承担对于控制器算法的实现，通过编写路径规划算法、无人机控制算法等C++代码，实现对无人机的位置、姿态等状态的控制；Gazebo仿真软件作为一个可视化的图形用户界面，需要在其中搭建仿真场景、传感器模型、仿真无人机等，通过Gazebo-ROS bridge与ROS系统建立连接，实现无人机模型状态的获取与发送，可视化地演示集群无人机队形重构的整体流程效果。平台各个部分描述见4.2节和4.3节。

4.2 搭建仿真场景

Gazebo仿真软件是一个功能完善的3D机器人仿真软件，通过Gazebo搭建无人机仿真的整体场景，包括搭建仿真环境和无人机模型，添加传感器、执行器等各种插件，实现ROS与Gazebo的连接等环节，主要操作流程如图 7所示。其中URDF文件是Gazebo-ROS中的一种用于描述机器人组件的XML文件，该文件的主要功能就是描述机器人各个部分(link)与关节(joint)的属性，创建机器人模型。.world文件是对仿真环境的描述，通过.launch文件启动并加载仿真场景，由此实现Gazebo中仿真模型的整体搭建。详细操作步骤参考Gazebo软件官方说明^[26]，此处不再赘述。

4.3 集群无人机飞行控制的实现

本文中基于无人机灯光秀场景，设计实现对集群无人机进行队形重构的整体流程演示，每架无人机都由相同的控制系统框架进行控制，总体仿真演示方案如图 8所示，其核心包括轨迹规划模块和轨迹跟踪模块，其中轨迹规划模块包括队形设计和路径规划2个子模块，队形设计子模块针对无人机表演的队形形状进行队形设计，路径规划子模块为每架无人机从起始点到目标队形最终点的轨迹进行规划设计；轨迹跟踪模块，主要利用本文第3节中的轨迹姿态控制器，通过内环的姿态控制和外环的位置控制，确保每架无人机都能够安全地、无碰撞地追踪给定轨迹朝着各自目标位置飞行。

4.4 仿真验证

以12架四旋翼无人机为例，对上述提出的路径规划算法及轨迹姿态控制器的有效性进行数值仿真验证。12架四旋翼无人机集群飞行的飞行示意图如图 9所示，两队形之间的详细队形重构过程如图 10所示。图中标注实时追踪2架无人机(无人机6和无人机7)的实时位置变化。由图 9可以看出仿真实现了连续的队形重构过程，显示了各队形重构后无人机的位置变化情况，其中由队形1到队形2的变换过程中无人机6的位置几乎不变，体现了队形位置最优分配保证效率最高的思想。由图 10可以看出在两队形重构过程中，无人机间均保持一定的飞行距离，确保相互之间不会发生碰撞。总结可得基于本文的控制算法，利用CAPT算法进行集群无人机路径规划，可以实现多无人机的连续、安全、无碰撞的队形重构过程。

5 结论

1) 研究了路径规划CAPT算法和无人机轨迹平滑方法，针对无人机集群队形重构的具体场景需求，实现多无人机安全迅速的并行目标分配和轨迹生成。

2) 针对无人机轨迹姿态跟踪问题，提出了一种有限时间多变量积分滑模控制算法，实现了对无人机期望轨迹的有限时间精确跟踪。

3) 基于Gazebo-ROS搭建了无人机仿真平台，基于实际物理引擎及多种动力学插件实现了集群无人机队形重构“建模-仿真-可视化”的一体化虚拟仿真演示，对上述控制算法及路径规划算法进行了有效性验证。