涵道风扇式无人机是一种特殊的垂直起降无人机，它的风扇置于涵道内部，结构紧凑，安全可靠。这些特性使其逐渐引起各国研究机构的重视^[1]，已有许多典型的机型投入使用，如iSTAR^[2]、HoverEye^[3]等。

在飞行器、船舶等机械系统的运动控制中，操纵面是用于产生作用于系统的力或者力矩的机械装置，如推进器、舵、螺旋桨等。为了满足容错控制和控制重构等需求、提高系统的灵活性和动态响应，通常设计冗余的操纵面^[4]。涵道风扇式无人机的操纵面是底部的气动舵面，气动舵面处于风扇滑流中，通过偏转舵面产生作用于涵道底部的侧向气动力来控制姿态。该类型无人机通常具有冗余的控制舵面，其姿态子系统是一个典型的过驱动(Over-actuated)系统，因此在控制中存在控制分配问题。

操纵面冗余是设计过驱动系统的控制器时要解决的问题之一，一种常见的方法是设计最优控制构成闭环系统，操纵面指令即控制输入已经包含在设计中。另一种设计方法是引入伪控制输入(力或力矩)将调节或跟踪任务和控制分配任务分开，针对调节或跟踪任务设计的控制律作为伪控制输入的期望值，仅指定要产生的力或力矩，然后设计单独的控制分配模块，根据期望的力或力矩求解控制输入，文献[5]讨论了2种设计方法的联系。这种分层设计的优点是可以单独设计控制律，而独立出来的控制分配模块不仅可以协调系统中不同操纵面的作用，还可以处理控制重构和操纵面容错之类的问题，例如文献[6]利用控制分配对过驱动的小型固定翼无人机进行容错控制。本文的讨论是基于这种分层设计方法进行的，并将力矩作为伪控制输入。

目前常用的控制分配方法包括伪逆法^[7-8]、直接分配法^[9]、基于最优化方法的控制分配^[10-11]等，文献[4]对控制分配领域的成果做了详尽的调查。近年来，这些成果在航空^[12]、航天^[13]、航海^[14]、汽车^[15]以及机器人^[16]等领域中得到广泛应用。对于涵道风扇式无人机的控制分配问题，目前主要是用伪逆法求解^[17-19]。由于实际约束的存在，伪逆法无法对任意可达的期望力矩都返回容许控制^[20]，因此并没有在控制分配上达到最优。

在控制分配环节，由于存在约束，操纵面只能产生可达力矩集(Attainable Moment Set, AMS)内的力矩。控制律作为期望力矩，通常设计时已经考虑了AMS的范围，但某些情况下控制律仍可能会超出该范围。例如在飞行包络线的边界附近，操纵面无法产生控制律所要求的力矩，将存在分配误差，通常控制分配算法将返回一个在某种指标下误差最小的次最优解^[4]。传统的控制分配算法返回次最优解时没有考虑控制律的分量的优先级问题，可能会使控制律的某些高优先级分量产生分配误差。

由于涵道风扇式无人机的特殊构型，在飞行中其期望力矩很可能超出AMS的范围。因此应用常规分配方法可能带来分配误差。为了既能对所有可达的期望力矩返回容许控制，又能对不可达的期望力矩做进一步的优化，本文在已有的研究成果的基础上引入一种优先级分配方法。该方法通过先将期望力矩进行矢量分解，并按一定优先级排序，对该优先级序列引入比例因子求解约束最优化问题，使得高优先级分量尽可能无误差分配。当期望力矩可达时，该方法等价于直接分配法，当期望力矩不可达时，该方法引入的比例因子对期望力矩的低优先级分量进行放缩，相当于截断低优先级分量，保证高优先级分量无分配误差。仿真显示所提出的算法对高优先级分量的跟踪效果优于传统算法。最后，本文将所提出的方法应用于涵道风扇式无人机的控制分配中，进行了飞行试验。该试验中，将控制律分解为2部分：一部分与系统抗扰动能力相关，作为高优先级分量；另一部分与瞬态响应相关，作为低优先级分量。试验结果表明，优先级分配法能保证可达的高优先级分量无误差分配，系统保持一定的输出解耦，而传统分配方法则因为执行器饱和而导致输出耦合。

1 控制分配问题 1.1 过驱动系统的控制分配问题

考虑非线性系统

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}) + \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{\delta }})}\\ {\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{h}}(\mathit{\boldsymbol{x}})} \end{array}} \right. $	（1）

式中：x∈Rⁿ为系统状态；δ∈Ω⊂R^p为系统的控制输入，或称为控制向量, Ω为客评控制集；映射f:Rⁿ→Rⁿ，g:Rⁿ×R^p→Rⁿ，h:Rⁿ→R^m。y∈R^m为系统输出。对过驱动系统有m < p。将g(x, δ)做矩阵分解

$ \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{\delta }}) = \mathit{\boldsymbol{GM}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{\delta }}) $	（2）

式中：G∈R^n×m，映射M:Rⁿ×R^p→R^m。引入伪控制输入：

$ \mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{\delta }}) $	（3）

将式(1)化为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}) + \mathit{\boldsymbol{G\tau }}}\\ {\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{h}}(\mathit{\boldsymbol{x}})} \end{array}} \right. $	（4）

式中：伪控制输入τ∈Φ⊂R^m通常为m维广义力或力矩，在本文中统一称τ为力矩。由于操纵面约束，Ω为R^p的子集，称为容许控制集，属于Ω的δ称为容许控制。容许控制集通过M从R^p空间映射到R^m空间得到系统的可达力矩集Φ，若τ∈Φ则称τ可达(Attainable)。

在分层设计方法中，首先由上层控制算法针对式(4)设计的控制律τ_c∈R^m，然后控制分配环节将τ_c作为期望力矩，对给定τ=τ_c，根据式(3)的操纵面模型求解δ∈Ω。注意到，使式(4)渐近稳定的控制律τ_c不一定可达。

实际中更实用的是线性化的操纵面模型。当控制向量δ为零向量时称为零偏转，记为δ₀，将式(3)在某系统状态下在δ₀处做泰勒展开

$ {\mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{B}}(\mathit{\boldsymbol{\delta }} - {\mathit{\boldsymbol{\delta }}_0}) + {\mathit{\boldsymbol{\tau }}_0}} $	（5）

$ {{\mathit{\boldsymbol{\tau }}_0} \buildrel \Delta \over = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{\delta }}_0})} $	（6）

$ \mathit{\boldsymbol{B}} \buildrel \Delta \over = \frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\delta }}}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {\tau _1}}}{{\partial {\delta _1}}}}&{\frac{{\partial {\tau _1}}}{{\partial {\delta _2}}}}& \cdots &{\frac{{\partial {\tau _1}}}{{\partial {\delta _p}}}}\\ {\frac{{\partial {\tau _2}}}{{\partial {\delta _1}}}}&{\frac{{\partial {\tau _2}}}{{\partial {\delta _2}}}}& \cdots &{\frac{{\partial {\tau _2}}}{{\partial {\delta _p}}}}\\ \vdots & \vdots & \ddots & \vdots \\ {\frac{{\partial {\tau _m}}}{{\partial {\delta _1}}}}&{\frac{{\partial {\tau _m}}}{{\partial {\delta _2}}}}& \cdots &{\frac{{\partial {\tau _m}}}{{\partial {\delta _p}}}} \end{array}} \right] $

（7）

式中：B∈R^m×p称为控制效率矩阵，B、τ₀都与系统状态相关。δ₀通常使τ₀为零向量。另外在具体实现中，控制系统以较高的频率运行，将式(3)在上一采样时刻做泰勒展开为

$ \mathit{\boldsymbol{\tau }} - {\mathit{\boldsymbol{\tau }}_{\rm{L}}} = \mathit{\boldsymbol{B}}(\mathit{\boldsymbol{\delta }} - {\mathit{\boldsymbol{\delta }}_{\rm{L}}}) \Rightarrow \Delta \mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{B}}\Delta \mathit{\boldsymbol{\delta }} $	（8）

式中：δ_L和τ_L为上一采样时刻的值。通常将M(x, δ)曲线的离散数据存在计算机中，通过查表可求得曲线斜率进而求得B，即使零偏转处的τ₀不为零向量时同样可以查表获得。因此可以不失一般性地将操纵面模型表示为

$ \mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{B\delta }} $	（9）

并基于式(9)讨论控制分配问题。至此，控制分配问题可描述为：对给定B、Ω及τ_c，求δ，使得τ_c=Bδ且δ∈Ω。

控制分配算法求解得δ后，将作为执行器的指令，执行器使操纵面跟踪该指令。这部分通常集成在一些机电设备中，如无人直升机中所用的力矩舵机。相比于被控对象，通常执行器的动态响应时间非常短，因此认为操纵面状态即为δ。引入控制分配环节后，闭环控制系统如图 1所示。

1.2 涵道风扇式无人机的控制分配问题

本节将1.1节的分层设计方法应用到涵道风扇式无人机中。有关建模问题主要参考文献[17-19]。地面坐标系X_NY_NZ_N、机体坐标系X_BY_BZ_B及按ZYX的顺序定义的欧拉角φ、θ、ψ，如图 2所示，坐标系原点为无人机重心。记绕机体坐标系转动的角速度为ω。在飞行器控制分配中，主要关注和控制输入直接相关的力矩方程^[17]，令式(4)中系统状态x=ω，有

$ \mathit{\boldsymbol{\dot \omega }} = {\mathit{\boldsymbol{I}}^{ - 1}}(\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} - \mathit{\boldsymbol{\omega }} \times \mathit{\boldsymbol{I\omega }}) $	（10）

式中：Γ=Γ_cs+Γ_ext，Γ_cs为操纵面产生的驱动力矩，Γ_ext为其他力矩^[17]；I=diag(I_x, I_y, I_z)为惯性张量。

定义操纵面旋转轴如图 3所示，其中1~4号气动舵面即飞行器操纵面。规定操纵面旋转符合右手系则为正值，否则为负值。零偏转定义在舵面平行于Z_B轴的位置。每个舵面最大转动角度均为±20°。力臂l₁、l₂如图 3所示。驱动力矩可表示为^[17]

$ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{{\rm{cs}}}} = \left[ {\begin{array}{*{20}{l}} {{\varGamma _x}}\\ {{\varGamma _y}}\\ {{\varGamma _z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {l_1}}&0&{{l_1}}&0\\ 0&{ - {l_1}}&0&{{l_1}}\\ {{l_2}}&{{l_2}}&{{l_2}}&{{l_2}} \end{array}} \right]{k_\delta }{V_{\rm{d}}}\mathit{\boldsymbol{\delta }} $	（11）

式中：δ=[δ₁  δ₂  δ₃  δ₄]^T为所有操纵面偏转角组成的控制向量；k_δ为操纵面气动力系数；V_d为涵道底部吹出的气流速度。

取驱动力矩为伪控制输入，即τ=Γ_cs，由式(9)可得

$ \mathit{\boldsymbol{B}} = {k_\delta }{V_{\rm{d}}}\left[ {\begin{array}{*{20}{c}} { - {l_1}}&0&{{l_1}}&0\\ 0&{ - {l_1}}&0&{{l_1}}\\ {{l_2}}&{{l_2}}&{{l_2}}&{{l_2}} \end{array}} \right] $	（12）

式(12)隐含B和系统状态相关。由式(12)易知该姿态子系统是一个过驱动系统，其控制分配问题可基于式(9)、式(12)进行的。讨论涵道风扇式无人机的控制分配问题时，是在如下假设下进行的：

1) 假设已经获得系统的控制律，并且在控制律下式(10)系统渐进稳定。本文的重点是在已知控制律的基础上研究控制分配问题。

2) 假设执行器动态响应时间远远小于被控对象。该无人机所用的力矩舵机带宽远大于被控对象，因此假设是合理的。

3) 假设各个状态下的控制效率矩阵B已经通过仿真或试验数据获得。控制向量和力矩曲线M(x, δ)以离散数据的形式保存在计算机中，系统运行时将通过查表，利用插值法获取B矩阵的元素。

2 控制分配方法

2.1节先介绍常用的2种控制分配方法：伪逆法和直接分配法，将其应用到涵道风扇式飞行器的控制分配问题中, 并对比2种分配方法的可达集。2.2节针对2种常用方法的不足，提出一种优先级分配方法。

2.1 常用控制分配方法

由于操纵面存在约束，容许控制集Ω通常是R^p空间中的凸集。记

$ \mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}}= \left[ {\begin{array}{*{20}{c}} {{\delta _{{\rm{1,min}}}}}\\ {{\delta _{{\rm{2,min}}}}}\\ \vdots \\ {{\delta _{p,{\rm{min}}}}} \end{array}} \right],\mathit{\boldsymbol{\bar \delta }} = \left[ {\begin{array}{*{20}{c}} {{\delta _{{\rm{1,max}}}}}\\ {{\delta _{{\rm{2,max}}}}}\\ \vdots \\ {{\delta _{p,{\rm{max}}}}} \end{array}} \right],\mathit{\boldsymbol{\delta }} = \left[ {\begin{array}{*{20}{c}} {{\delta _1}}\\ {{\delta _2}}\\ \vdots \\ {{\delta _p}} \end{array}} \right] $

（13）

式中：对∀i=1, 2, …, p，若δ_{i, min}≤δ_i≤δ_{i, max}，则表示为δ-≤δ≤δ，则容许控制集为

$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = \{ \mathit{\boldsymbol{\delta }}|\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} \le \mathit{\boldsymbol{\delta }} \le \mathit{\boldsymbol{\bar \delta }}\} $	（14）

式(9)所示操纵面的静态线性模型的重要性在于，可以在一个全局的范围观察一个系统的AMS，AMS是系统的属性，和控制分配算法无关。在具体实现中，在系统的每一个控制周期，控制分配采用的是式(8)，除幅值约束外，还要考虑操纵面的速率约束。容许控制集变为幅值约束和速率约束的交集，相当于在每一个控制周期都有一个局部的幅值约束^[20]。若记上一控制周期的控制向量为δ_L，控制周期为T，操纵面最大执行速率为u_m，有局部约束

$ - T{\mathit{\boldsymbol{u}}_{\rm{m}}} \le \mathit{\boldsymbol{\delta }} - {\mathit{\boldsymbol{\delta }}_{\rm{L}}} \le T{\mathit{\boldsymbol{u}}_{\rm{m}}} $	（15）

考虑式(13)，取交集

$ {{{\mathit{\boldsymbol{\bar \delta }}}^\prime } = {\rm{min}}\{ \mathit{\boldsymbol{\bar \delta }} - {\mathit{\boldsymbol{\delta }}_{\rm{L}}},\quad T{\mathit{\boldsymbol{u}}_{\rm{m}}}\} } $	（16）

$ {{\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}}^\prime } = {\rm{max}}\{ \mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} - {\mathit{\boldsymbol{\delta }}_{\rm{L}}}, - T{\mathit{\boldsymbol{u}}_{\rm{m}}}\} } $	（17）

记Δδ=δ-δ_L，则在每一控制周期进行控制分配时，容许控制集为

$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = \{ \mathit{\boldsymbol{\delta }}|{\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}}^\prime } \le \mathit{\boldsymbol{\delta }} \le {\mathit{\boldsymbol{\bar \delta }}^\prime }\} $	（18）

结合式(8)，控制分配问题和仅存在幅值约束的情形类似，区别是前者分配的起点是上一控制周期的控制向量δ_L，而后者分配的起点是R^p空间的原点。考虑速率约束是系统的具体实现问题，仅考虑幅值约束得到的结论可应用到考虑速率约束的情形，因此在本文的讨论中不提及操纵面的速率约束，而在具体实现中默认考虑了速率约束。

1) 伪逆法

若不考虑任何操纵面约束，对给定τ_c，取δ=B⁺τ_c，其中B⁺=B^T(BB^T)^-1为B的伪逆(Pseudo-inverse)，这种控制分配称为伪逆法。文献[17-19]在涵道风扇式无人机上所用的控制分配方法都等价于伪逆法。由于操纵面约束的存在，对任意给定τ_c∈Φ，δ可能不在容许控制集内，即δ∉Ω，实际执行的控制输入是δ经过限幅后的值。

2) 直接分配法

为了对任意给定τ_c∈Φ，都返回一个容许控制，设τ_max∈Φ是AMS中矢量τ_c方向上幅值最大的矢量，定义比例因子α=‖τ_max‖²/‖τ_c‖²。求解一组α、δ满足ατ_c=Bδ且δ∈Ω。用最优化问题的形式来表达，可表示为^[21]

$ \begin{array}{l} \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{\delta }},\alpha } {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \alpha \\ {\rm{s}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\begin{array}{*{20}{l}} {\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} \le \mathit{\boldsymbol{\delta }} \le \mathit{\boldsymbol{\bar \delta }}}\\ {\alpha {\mathit{\boldsymbol{\tau }}_{\rm{c}}} = \mathit{\boldsymbol{B\delta }}}\\ {0 \le \alpha \le 1} \end{array}} \right. \end{array} $

（19）

求解式(19)即为直接分配法。具体算法实现中，α的取值可以放宽为大于等于零的实数，将式(19)化为线性规划问题的标准形式，用单纯形法求解该最优化问题^[21]。该方法可对任意可达的τ_c无误差分配，即返回一个容许控制，使τ_c=Bδ。对不可达的τ_c，仅能满足ατ_c=Bδ，0≤α < 1，分配误差为τ_c-Bδ。比例因子的作用是减小τ_c的幅值，因此该方法具有方向保持的特性。

在对比2种分配算法之前先引入所需的一些符号。用∂(·)表示包围某点集的凸包，包围容许控制集、AMS的凸包分别记为∂(Ω)、∂(Φ)。对所有可达的τ_c，某一分配算法返回的所有容许控制的集合记为Θ。记Π为Θ通过控制效率矩阵B从R^p空间映射到R^m空间得到的集合，称为某分配算法的可达集。对m=3的情形，使用Π和Φ的体积之比的百分数表示某分配算法的品质因数γ，作为衡量算法优劣的指标^[20]。

对本文所研究的涵道风扇式无人机，根据飞行试验测得在悬停状态下的B矩阵为

$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} { - 0.539{\kern 1pt} {\kern 1pt} {\kern 1pt} 3}&0&{0.539{\kern 1pt} {\kern 1pt} {\kern 1pt} 3}&0\\ 0&{ - 0.539{\kern 1pt} {\kern 1pt} {\kern 1pt} 3}&0&{0.539{\kern 1pt} {\kern 1pt} {\kern 1pt} 3}\\ {0.209{\kern 1pt} {\kern 1pt} {\kern 1pt} 9}&{0.209{\kern 1pt} {\kern 1pt} {\kern 1pt} 9}&{0.209{\kern 1pt} {\kern 1pt} {\kern 1pt} 9}&{0.209{\kern 1pt} {\kern 1pt} {\kern 1pt} 9} \end{array}} \right] $

（20）

其操纵面约束为

$ \mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} = \left[ {\begin{array}{*{20}{l}} { - {{20}^\circ }}\\ { - {{20}^\circ }}\\ { - {{20}^\circ }}\\ { - {{20}^\circ }} \end{array}} \right],\mathit{\boldsymbol{\bar \delta }} = \left[ {\begin{array}{*{20}{c}} {{{20}^\circ }}\\ {{{20}^\circ }}\\ {{{20}^\circ }}\\ {{{20}^\circ }} \end{array}} \right] $

（21）

分别绘制出伪逆法和直接分配法的可达集Π₁、Π₂，如图 4所示。由算法流程易知直接分配法的可达集Π₂=Φ，分别计算Π₁、Π₂的品质因数为γ₁=71.99%，γ₂=100%。对Φ中的所有力矩，任何无法返回容许控制的控制分配方法都没有充分利用飞机的控制性能^[20]，因此在涵道风扇式无人机的控制分配中广泛应用的伪逆法，实际上牺牲了部分控制能力。

2.2 优先级控制分配方法

由2.1节的分析，直接分配法虽然品质因数比伪逆法更大，但对不可达的τ_c，将被保方向地限制在∂(Φ)上，其分配误差为(1-α)τ_c，即每一个分量都有分配误差。然而，在某些情况下为了满足控制需求，有些分量应尽可能无分配误差。为此提出一种优先级分配方法，先将期望力矩进行优先级分解，若τ_c不可达，将其与∂(Φ)的交点进一步优化，使得τ_c的高优先级分量尽可能无误差分配，仅低优先级分量产生分配误差。

2.2.1 控制律的优先级分解

期望力矩由式(4)所示系统的控制律τ_c给出，τ_c∈R^m，可以在R^m空间按坐标轴或不同控制效果做矢量分解，分解后按控制目标划分优先级。

本节以无人机中应用广泛的主动扰动抑制控制(ADRC)^[22]为例，该方法的核心是利用扩张状态观测器(ESO)对“总和扰动”进行估计并补偿，使对象模型变成“积分器串联型”线性系统。将式(10)改写为

$ \left\{ {\begin{array}{*{20}{l}} {\ddot \varphi = {f_1}(\varphi ,\theta ,\psi ,\dot \varphi ,\dot \theta ,\dot \psi ) + {W_1} + I_x^{ - 1}{\tau _x}}\\ {\ddot \theta = {f_2}(\varphi ,\theta ,\psi ,\dot \varphi ,\dot \theta ,\dot \psi ) + {W_2} + I_y^{ - 1}{\tau _y}}\\ {\ddot \psi = {f_3}(\varphi ,\theta ,\psi ,\dot \varphi ,\dot \theta ,\dot \psi ) + {W_3} + I_z^{ - 1}{\tau _z}}\\ {\mathit{\boldsymbol{y}} = {{\left[ {\begin{array}{*{20}{l}} \varphi &\theta &\psi \end{array}} \right]}^{\rm{T}}}} \end{array}} \right. $

（22）

式中：${\left[ {{\tau _x}\;\;{\tau _y}\;\;{\tau _z}} \right]^{\rm{T}}} \buildrel \Delta \over = \mathit{\boldsymbol{\tau }}, {\left[ {{W_1}\;\;{W_2}\;\;{W_3}} \right]^{\rm{T}}} \buildrel \Delta \over = \mathit{\boldsymbol{W}}{\rm{ }} $为外部扰动，$ {f_i}\left( {\varphi , {\rm{ }}\theta , {\rm{ }}\psi , {\rm{ }}\dot \varphi , {\rm{ }}\dot \theta , {\rm{ }}\dot \psi } \right)\left( {i = 1, 2, 3} \right)$为由旋转产生的非线性项，简记为f_i。下面以滚转通道为例，给出ADRC控制律。ADRC控制器包含ESO、跟踪微分器(TD)、非线性误差反馈(NLSEF)，其离散形式可分别设计为

$ \left\{ {\begin{array}{*{20}{l}} {e = [\varphi (k) - {z_{{\rm{r1}}}}(k)]/{\varepsilon ^2}}\\ {{z_{{\rm{r1}}}}(k + 1) = {z_{{\rm{r1}}}}(k) + T({z_{{\rm{r2}}}}(k) + \varepsilon {g_1}(e))}\\ {{z_{{\rm{r2}}}}(k + 1) = {z_{{\rm{r2}}}}(k) + T({z_{{\rm{r3}}}}(k) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {g_2}(e) + I_x^{ - 1}{\tau _x}(k))}\\ {{z_{{\rm{r3}}}}(k + 1) = {z_{{\rm{r3}}}}(k) + T{g_3}(e)/\varepsilon } \end{array}} \right. $

（23）

$ \left\{ \begin{array}{l} {v_{{\rm{r1}}}}(k + 1) = {v_{{\rm{r1}}}}(k) + T{v_{{\rm{r2}}}}(k)\\ {v_{{\rm{r2}}}}(k + 1) = {v_{{\rm{r2}}}}(k) + T{\rm{fst}}({v_{{\rm{r1}}}}(k) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varphi _{\rm{c}}},{v_{{\rm{r2}}}}(k),r,h) \end{array} \right. $

（24）

$ \left\{ {\begin{array}{*{20}{l}} {{e_1} = {v_{{\rm{r1}}}}(k) - {z_{{\rm{r1}}}}(k)}\\ {{e_2} = {v_{{\rm{r2}}}}(k) - {z_{{\rm{r2}}}}(k)}\\ {{\tau _{x{\rm{c}}}} = {k_1}{\rm{fal}}({e_1},{\sigma _1},\xi ) + {k_2}{\rm{fal}}({e_2},{\sigma _2},\xi ) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {I_x}{z_{{\rm{r3}}}}(k)} \end{array}} \right. $

（25）

式中：T为控制周期；φ_c为给定滚转角，且函数

$ \left\{ {\begin{array}{*{20}{l}} {{g_i}(e) = {{\rm{ \mathsf{ β} }}_i}{{[e]}^{i{\zeta _i} - (i - 1)}}}&{i = 1,2,3}\\ {{{[e]}^\chi } \buildrel \Delta \over = {\rm{sign}} (e)|e{|^\chi }} \end{array}} \right. $	（26）

$ \left\{ \begin{array}{l} {\rm{fst}} ({v_{{\rm{r1}}}},{v_{{\rm{r2}}}},r,h) = - \left\{ {\begin{array}{*{20}{l}} {ra/d}&{|a| \le d}\\ {r {\rm{sign}} (a)}&{|a| > d} \end{array}} \right.\\ a \buildrel \Delta \over = \left\{ {\begin{array}{*{20}{l}} {{v_{{\rm{r2}}}} + {\rm{sign}} (n)({a_0} - d)/2}&{|n| > {d_0}}\\ {{v_{{\rm{r2}}}} + n/h}&{|n| \le {d_0}} \end{array}} \right.\\ \begin{array}{*{20}{l}} {d \buildrel \Delta \over = rh,{d_0} \buildrel \Delta \over = dh,n \buildrel \Delta \over = {v_{{\rm{r1}}}} + h{v_{{\rm{r2}}}}}\\ {{a_0} \buildrel \Delta \over = {{({d^2} + 8r|n|)}^{1/2}}} \end{array} \end{array} \right. $

（27）

$ {\rm{fal}}(e,\sigma ,\xi ) = \left\{ {\begin{array}{*{20}{l}} {|e{|^\sigma } {\rm{sign}} (e)}&{|e| > \xi }\\ {e/{\xi ^{1 - \sigma }}}&{|e| \le \xi ,} \end{array}\xi > 0} \right. $	（28）

TD以及NLSEF的参数 r、h、σ₁、σ₂、ξ、k₁、k₂的选取原则可参考文献[22]。ESO的输出z_r1、z_r2、z_r3将在有限时间内收敛于φ、${\dot \varphi } $、f₁+W₁，参数ε、$\zeta i $、β_i的选取原则以及其收敛性证明参考文献[23]。对滚转通道的控制律τ_xc做优先级分解τ_xc=τ_x1+τ_x2，其中

$ {{\tau _{x1}} \buildrel \Delta \over = - {I_x}{z_{{\rm{r3}}}}(k)} $	（29）

$ {{\tau _{x2}} \buildrel \Delta \over = {k_1}{\rm{fal}}({e_1},{\sigma _1},\xi ) + {k_2}{\rm{fal}}({e_2},{\sigma _2},\xi )} $	（30）

同理，对俯仰、偏航通道设计ADRC控制器，得控制律τ_yc、τ_zc，同样做类似的分解，得总的控制律

$ {\mathit{\boldsymbol{\tau }}_{\rm{c}}} = {\mathit{\boldsymbol{\tau }}_1} + {\mathit{\boldsymbol{\tau }}_2} $	（31）

$ {\mathit{\boldsymbol{\tau }}_{\rm{c}}} \buildrel \Delta \over = \left[ {\begin{array}{*{20}{l}} {{\tau _{x{\rm{c}}}}}\\ {{\tau _{y{\rm{c}}}}}\\ {{\tau _{z{\rm{c}}}}} \end{array}} \right],{\tau _1} \buildrel \Delta \over = \left[ {\begin{array}{*{20}{l}} {{\tau _{x1}}}\\ {{\tau _{y1}}}\\ {{\tau _{z1}}} \end{array}} \right],{\tau _2} \buildrel \Delta \over = \left[ {\begin{array}{*{20}{l}} {{\tau _{x2}}}\\ {{\tau _{y2}}}\\ {{\tau _{z2}}} \end{array}} \right] $

（32）

其中：τ₁包含对系统总和扰动的补偿，τ₂通常是应用状态反馈设计的，因此和系统动态响应相关。在该控制律下，系统能按照预期运行的首要前提是τ₁可以抵消系统总和扰动，将非线性系统化为线性系统，而后才可以设计误差反馈τ₂，使系统稳定(或跟踪输入)。因此将控制律分解后，可以认为τ₁分量优先级较高，而τ₂分量优先级较低。

应当指出，这些分解是人为规定的，还可以根据不同需求做其他分解。基于反馈线性化方法设计的控制律中，都可以做和上述方法类似的分解，因为其本质都是先将非线性系统化为线性系统，再对该线性系统设计误差反馈律，例如非线性动态逆(NDI)及增量非线性动态逆(INDI)^[24]，其思路和ADRC是相似的。文献[25-26]讨论了NDI控制律的分解，以及文献[27]讨论了串级PID控制律的分解。

2.2.2 优先级控制分配

在控制分配环节控制律作为期望力矩，假设期望力矩已经分解为τ_c=τ₁+τ₂+…+τ_k，其中从τ₁~τ_k优先级依次降低。为保证高优先级分量不产生分配误差，考虑仅对优先级最低的τ_k进行缩放，引入比例因子α，求解式(33):

$ \begin{array}{l} \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{\delta }},\alpha } \alpha \\ {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\begin{array}{*{20}{l}} {\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} \le \mathit{\boldsymbol{\delta }} \le \mathit{\boldsymbol{\bar \delta }}}\\ {\mathit{\boldsymbol{B\delta }} = {\mathit{\boldsymbol{\tau }}_1} + {\mathit{\boldsymbol{\tau }}_2} + \cdots + \alpha {\mathit{\boldsymbol{\tau }}_k}}\\ {0 \le \alpha \le 1} \end{array}} \right. \end{array} $

（33）

1) 若式(33)有解

则控制分配结束，返回该最优解δ。此时若α=1，表示τ_c可达，求解式(33)等价于对τ_c用直接分配法求解容许控制。若0≤α < 1，表示τ_c不可达但τ₁+τ₂+…+ατ_k在∂(Φ)上，仅改变τ_k的幅值即可返回容许控制，因此仅在τ_k方向有分配误差。

2) 若式(33)无解

则在τ_c中把低优先级分量τ_k去掉，此时期望力矩变为τ′_c=τ₁+τ₂+…+τ_k-1，此种情况表明τ′_c不可达。同理，仅对优先级最低的τ_k-1进行缩放，求解式(34):

$ \begin{array}{l} \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{\delta }},\alpha } \alpha \\ {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\begin{array}{*{20}{l}} {\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} \le \mathit{\boldsymbol{\delta }} \le \mathit{\boldsymbol{\bar \delta }}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{B\delta }} = {\mathit{\boldsymbol{\tau }}_1} + {\mathit{\boldsymbol{\tau }}_2} + \cdots + \alpha {\mathit{\boldsymbol{\tau }}_{k - 1}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0 \le \alpha \le 1} \end{array}} \right. \end{array} $

（34）

若有解则分配结束，若无解，则去掉τ_k-1重复上述过程。注意到若重复了k-1次该过程之后仍然无解，说明τ₁不可达，第k次求解时等价于对τ₁用直接分配法求解控制向量。因此该问题在有限次循环后一定有解。将整个过程表示为从s=k开始求解式(35):

$ \begin{array}{l} \mathop {{\rm{max}}}\limits_{\mathit{\boldsymbol{\delta }},{\alpha _s}} {\alpha _s}\\ {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {\begin{array}{*{20}{l}} {\mathrm{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\mathit{\boldsymbol{{ }\!\!\delta\!\!\text{ } }}}} \le \mathit{\boldsymbol{\delta }} \le \mathit{\boldsymbol{\bar \delta }}}\\ \begin{array}{l} \mathit{\boldsymbol{B\delta }} = {\alpha _1}{\mathit{\boldsymbol{\tau }}_1} + {\alpha _2}{\mathit{\boldsymbol{\tau }}_2} + \cdots + {\alpha _k}{\mathit{\boldsymbol{\tau }}_k}\\ {\alpha _j} = \left\{ {\begin{array}{*{20}{l}} 0&{j < s}\\ 1&{j > s} \end{array},j = 1,2, \cdots ,k} \right. \end{array}\\ {{\kern 1pt} 0 \le {\alpha _s} \le 1} \end{array}} \right. \end{array} $

（35）

式中:s=k, k-1, …, 1，若无可行解，则令s=s-1，循环求解式(35)直到有可行解时停止，即为优先级分配方法。当τ_c可达时，求解式(35)等价于直接分配法；当τ_c不可达时，若在i次求解式(35)后返回最优解，则τ₁, τ₂, …, τ_k-i都是无误差分配，其中0≤i≤k。当控制律中没有做优先级分解时，等价于τ_c=τ₁，优先级分配法退化为直接分配法。从式(35)可知该方法的可达集和直接分配法一样，品质因数为100%。

在具体实现中，为了减小计算量，求解式(35)前，可先用伪逆法计算一个参考值，若该参考值已经使操纵面饱和，再求解式(35)，否则输出该参考值即可。应该指出，这一步参考值的计算不是必须的。控制分配流程图如图 5所示，虚线方框内是可选的步骤，蓝色方框内是优先级分配方法流程。

图 6给出了一种2维力矩空间的例子，若τ_c=τ₁+τ₂，且m=2，τ_c不可达。分别用直接分配法、优先级分配法求解控制向量，再通过B映射到R²空间得τ_d、τ_p。将τ_d和τ_p沿τ₁、τ₂方向做矢量分解，如图 6所示。τ_d在每个分量方向上都有分配误差。τ_p则仅在低优先级分量τ₂方向上有分配误差，而在高优先级分量τ₁方向没有分配误差。

3 仿真验证

给定期望力矩，使用不同控制分配算法进行仿真。设m=3，p=4，控制效率矩阵B和操纵面约束分别如式(20)、式(21)所示。并设已按优先级分解的期望力矩为τ_c=τ₁+τ₂，其中高优先级分量τ₁=[0   0  0.25]^T、低优先级分量τ₂=[0.34cos(πt)  0.34sin(πt)  0]^T，该期望力矩随时间变化，涵盖了可达和不可达的情形。利用Simulink的Nonlinear Second-Order Actuator模块作为操纵面执行机构进行仿真。求解控制分配时采用式(8)，速率约束为400 (°)/s，操纵面起始状态为R^p原点。分别用伪逆法、直接分配法、优先级分配法求解控制向量，并通过B映射到力矩空间，得到实际执行的伪控制输入，结果如图 7所示。和上述分析一致，相对于伪逆法和直接分配法，优先级分配法虽然对低优先级分量τ₂的分配误差最大，但对高优先级分量τ₁无分配误差。

本文所研究的涵道风扇式无人机的飞行控制的MATLAB & Simulink仿真程序已开源，地址：https://github.com/mengchaoheng/Plan-D，在该仿真中可以看到更多控制分配算法的对比及其与不同控制器结合的效果。

4 飞行试验

在图 8所示的涵道风扇式无人机上做不同分配算法的对照试验，无人机参数如表 1所示，试验时应用ADRC进行姿态控制，控制律为如式(31)所示，闭环系统如图 9所示，其中扰动观测器ESO的输出z₁、z₂、z₃分别表示估计的欧拉角、欧拉角的导数、系统总和扰动。Iz₃即估计的扰动力矩，扰动补偿为τ₁=-Iz₃。每个通道的NLSEF类似地由式(30)给出，总的NLSEF的输出为τ₂。TD的输出v₁、v₂表示给定欧拉角[φ_c  θ_c  ψ_c]^T的过渡过程及其微分信号。试验中除控制分配算法不同外，控制律、参考指令、飞行环境均相同。系统控制周期T=0.01 s，在定高且锁尾模式下飞行。

将τ₁设为最高优先级，τ₂设为低优先级。操纵面物理约束如式(21)所示，为方便后续试验中添加9°的操纵面扰动，试验时在控制程序中统一设定操纵面最大偏转角为11°。如前文所述，忽略执行器动力学，舵面偏转角用其执行器的指令δ近似。姿态角参考输入为滚转、俯仰通道20°的阶跃信号，航向角保持起飞时锁定的角度，统一为正北方向(零度偏航角)。

若不添加扰动，分别使用伪逆法、直接分配法、优先级法进行控制分配，姿态角的阶跃响应曲线如图 10所示，相应的操纵面偏转角响应曲线如图 11所示。由图 10可知3种方法下的阶跃响应性能指标比较接近。比较3种方法的操纵面曲线，可知品质因数较小的伪逆法更容易饱和，而直接分配法和优先级分配法由于其可达集是整个AMS边界所包围的力矩空间，因而操纵面饱和的时间相对较短。

若添加恒定扰动δ_disturb=[9°   9°  9°  9°]^T到操纵面，用于产生偏航通道的扰动力矩，ESO的输出τ₁将主要用于补偿该力矩。待系统进入稳态后开始计时，姿态角参考输入同上文。伪逆法和直接分配法的姿态角响应曲线如图 12所示，优先级分配法的姿态角响应曲线如图 13所示，3种方法相应的操纵面偏转角响应曲线如图 14所示。

由图 12可知伪逆法和直接分配法都在输入参考指令后航向角由于扰动的影响持续减小，最终使系统发散，3 s后已经由回收设备介入。而由图 13可知优先级分配法使系统仍稳定跟踪参考输入，和图 10的结果相比仅响应速度变慢。由图 14可知，系统在阶跃信号到来之前，操纵面偏转角稳定在δ_disturb附近，表明τ₁是可达的。在第1 s~第1.5 s的瞬态响应期间，操纵面饱和即表明参考指令使得误差反馈输出的τ₂叠加τ₁的结果τ₁+τ₂不可达。在优先级分配法中，τ₁被设置为高优先级，使操纵面仅在补偿扰动所需的偏转角附近动作，保证高优先级分量尽可能的得到执行，因此始终能抵抗扰动的影响。而伪逆法和直接分配法都由于误差反馈的影响，操纵面在其约束范围内动作，高优先级分量的执行得不到保证，最终导致输出耦合，表现为输入阶跃参考指令后系统发散。

5 结论

在涵道风扇式无人机的控制分配中，对本文所讨论的几种控制分配方法有如下结论：

1) 优先级分配法、直接分配法的品质因数都为100%，而常规的伪逆法仅为71.99%。

2) 优先级分配法可以保证控制律中高优先级分量尽可能无误差分配，而仅在低优先级分量上产生分配误差，具体表现为使系统可一定程度上防止因操纵面约束引起的系统输出耦合，尽可能保持输出解耦，这是直接分配法和伪逆法所没有的特性。