一种基于零空间的可倾转四旋翼无人机控制分配方法及系统

1.本发明涉及无人机控制技术领域，具体涉及一种基于零空间的可倾转四旋翼无人机控制分配方法及系统。


背景技术：

2.标准四旋翼无法实现位置和姿态的解耦控制，这一劣势限制了其在挑战性任务中的表现，例如，以特定姿态飞行、在悬停状态改变姿态及其他交互任务。为了克服这一缺陷，近年来，一种旋翼可倾转的过驱动四旋翼无人机受到了研究人员的广泛关注。这种无人机的冗余执行器可提供全向推力，从而提高系统的机动性及容错性。然而，冗余配置也导致系统的控制分配存在无数解。控制分配可以根据期望位姿指令将力和转矩合理地分配在执行器上。因此，如何确定控制分配问题的最优解，从而充分利用飞行器的机动能力是一个挑战性的问题。
3.近年来，在过驱动系统的控制分配策略方面已有一定的研究成果。主要的方法分为两类，分别是直接分配方法和数值优化方法。直接分配法利用控制效率矩阵的伪逆确定虚拟控制量的唯一映射，随后利用变量替换法实现该非线性映射的线性化。这种方法的优势在于计算速度快且通常可获得实际控制量的精确解。然而，在飞行器跟踪高机动轨迹时，直接分配方法的解可能超出执行器的速率极限和位置极限，进而使执行器过早进入饱和状态。此外，执行器在奇异点附近变化时，虽然可能满足极限约束，但跳变可能会导致飞行器姿态不稳定。此时则需利用数值优化方法处理多个不同类型的约束条件，然而常规的优化方法在目标函数的近似过程中会引入误差，从而使解算结果不可避免地偏离精确解。


技术实现要素：

4.鉴于以上问题，本发明提出一种基于零空间的可倾转四旋翼无人机控制分配方法及系统，用以解决现有的可倾转四旋翼无人机控制分配方法不能得到准确的执行器控制指令的问题。
5.根据本发明的一方面，提供一种基于零空间的可倾转四旋翼无人机控制分配方法，该方法包括以下步骤：
6.步骤一、根据牛顿-欧拉方程建立可倾转四旋翼无人机的动力学模型；所述动力学模型表示为：
[0007][0008]
其中，m表示无人机飞行器的质量；j表示无人机飞行器的惯性矩阵；表示无人机飞行器在机体坐标系下的质心线性加速度；表示无人机飞行器在机体坐标系下的角加速度；g表示重力加速度；向量e3＝[0；0；1]；ω表示无人机飞行器在机体坐标系下的角速度；f
drag
和m
drag
分别表示飞行器质心处受到的空气阻力及力矩；表示从机体坐标系到世界坐
标系的旋转矩阵，i3表示3阶单位矩阵；uf(αi,ni)和um(αi,ni)分别表示飞行器质心处受到的合外推力及力矩，αi表示第i号旋翼绕倾转轴的倾转角，ni表示第i号旋翼的转速，i＝1、2、3、4；
[0009]
步骤二、对于给定的无人机期望飞行轨迹，根据所述动力学模型获得的期望控制输入ud(αi,ni)为：
[0010][0011]
将期望控制输入ud(αi,ni)分解为：含时变参数的控制效率矩阵γ(αi)左乘旋翼推力向量t(ni)，即：
[0012]
ud(αi,ni)＝γ(αi)t(ni)
[0013]
步骤三、将控制效率矩阵γ(αi)中的时变参数αi转移到旋翼推力向量t(ni)中，以获得虚拟推力向量n(αi,ni)，并建立虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系，所述等式关系为：
[0014][0015]
式中，矩阵b表示不含时变参数的控制效率矩阵，表示矩阵b的右伪逆矩阵；v表示矩阵b的零空间的一组正交基；k表示调整因子且在飞行初始时刻的值为0；
[0016]
步骤四、对虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系进行求解，获得四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni，从而解算出无人机的期望控制输入ud(αi,ni)。
[0017]
进一步地，步骤一中飞行器质心处受到的空气阻力表示为：
[0018][0019]
其中，cf≥0表示空气阻力系数；ρ
air
表示空气密度；s
uav
表示飞行器迎风面积；表示飞行器在机体坐标系下的质心线速度；
[0020]
飞行器质心处受到的空气阻力矩m
drag
表示为：
[0021][0022]
其中，cm≥0表示空气阻力矩系数；
[0023]
飞行器质心处受到的合外推力uf(αi,ni)表示为：
[0024][0025]
其中，kf≥0表示旋翼推力系数；表示机体坐标系下螺旋桨推力的列矢量，
[0026][0027]
λi表示机体坐标系的x轴沿顺时针方向旋转到第i号旋翼的夹角，且：
[0028]
飞行器质心处受到的合外推力矩表示为：
[0029][0030]
其中，km≥0表示旋翼反扭矩系数；表示机体坐标系下螺旋桨几何中心的坐标，
[0031][0032]
其中，l表示机体坐标系原点到旋翼坐标系原点的距离。
[0033]
进一步地，步骤四中获得的四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni即共八个元素中，若超过无人机执行器物理极限的元素个数为一个或两个，则采用非线性重分配方法找到所述等式关系的一组精确解：
[0034][0035]
式中n即为虚拟推力向量n(αi,ni)。
[0036]
进一步地，采用非线性重分配方法找到所述等式关系的一组精确解的具体步骤包括：将虚拟推力向量n(αi,ni)分块为：其中，n
d1
对应推力的超限部分，n
d2
对应剩余部分；同时，对正交基矩阵v也进行分块，分块规则与虚拟推力向量n(αi,ni)的分块规则相对应；则虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系变为：
[0037][0038]
其中，v1、v2分别代表对应于n
d1
、n
d2
的正交基分块矩阵；进而得到调整因子k的计算式为：k＝v
1-(n
d1-n1)，将调整因子k重新代入上述矩阵等式，获得重分配的虚拟推力向量：
[0039][0040]
最后经非线性变换得到执行器控制指令的精确解：
[0041][0042]
进一步地，步骤四中获得的四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni即共八个元素中，若超过无人机执行器物理极限的元素个数为两个以上，则采用二次规划法找到所述等式关系的一组接近于精确解的最优解：
[0043][0044]
式中，α
i0
,n
i0
分别表示上一时刻的各个旋翼绕各自倾转轴的倾转角和旋翼转速；
△
αi、
△
(n
i2
)分别表示各个旋翼绕各自倾转轴的倾转角变化率和旋翼转速变化率。
[0045]
进一步地，步骤四中采用二次规划法找到所述等式关系的一组接近于精确解的最优解的具体步骤包括：
[0046]
构造二次规划法所需的等式约束：
[0047][0048]
式中，x＝[
△
αi,
△
(n
i2
)]表示控制指令的变化率向量；d为高阶非线性项；
[0049]
构造二次规划法所需的目标函数：
[0050]
g(x,d,k)＝x
t
px+d
t
qd+k
t
rk
[0051]
式中，p,q,r分别是规模为8维、8维、2维的参数对角矩阵；
[0052]
构造二次规划法所需的极值约束：
[0053][0054]
式中，
△
α
max
表示倾转角变化率的极限值，n
max
表示转速的极限值；
[0055]
将上述等式约束、目标函数、极值约束输入二次规划求解器，获得当前时刻的
△
αi,
△
ni，以获得最优解αi和ni。
[0056]
根据本发明的另一方面，提供一种基于零空间的可倾转四旋翼无人机控制分配系统，该系统包括：
[0057]
模型建立模块，其配置成根据牛顿-欧拉方程建立可倾转四旋翼无人机的动力学模型；所述动力学模型表示为：
[0058][0059]
其中，m表示无人机飞行器的质量；j表示无人机飞行器的惯性矩阵；表示无人机飞行器在机体坐标系下的质心线性加速度；表示无人机飞行器在机体坐标系下的角加速度；g表示重力加速度；向量e3＝[0；0；1]；ω表示无人机飞行器在机体坐标系下的角速度；f
drag
和m
drag
分别表示飞行器质心处受到的空气阻力及力矩；表示从机体坐标系到世界坐标系的旋转矩阵，i3表示3阶单位矩阵；uf(αi,ni)和um(αi,ni)分别表示飞行器质心处受到的合外推力及力矩，αi表示第i号旋翼绕倾转轴的倾转角，ni表示第i号旋翼的转速，i＝1、2、3、4；其中，
[0060]
飞行器质心处受到的空气阻力表示为：
[0061][0062]
其中，cf≥0表示空气阻力系数；ρ
air
表示空气密度；s
uav
表示飞行器迎风面积；表示飞行器在机体坐标系下的质心线速度；
[0063]
飞行器质心处受到的空气阻力矩m
drag
表示为：
[0064][0065]
其中，cm≥0表示空气阻力矩系数；
[0066]
飞行器质心处受到的合外推力uf(αi,ni)表示为：
[0067][0068]
其中，kf≥0表示旋翼推力系数；表示机体坐标系下螺旋桨推力的列矢量，
[0069][0070]
λi表示机体坐标系的x轴沿顺时针方向旋转到第i号旋翼的夹角，且：
[0071]
飞行器质心处受到的合外推力矩表示为：
[0072][0073]
其中，km≥0表示旋翼反扭矩系数；表示机体坐标系下螺旋桨几何中心的坐标，
[0074][0075]
其中，l表示机体坐标系原点到旋翼坐标系原点的距离；
[0076]
期望控制输入变换模块，其配置成对于给定的无人机期望飞行轨迹，根据所述动力学模型获得的期望控制输入ud(αi,ni)为：
[0077][0078]
将期望控制输入ud(αi,ni)分解为：含时变参数的控制效率矩阵γ(αi)左乘旋翼推力向量t(ni)，即：
[0079]
ud(αi,ni)＝γ(αi)t(ni)；
[0080]
将控制效率矩阵γ(αi)中的时变参数αi转移到旋翼推力向量t(ni)中，以获得虚拟推力向量n(αi,ni)，并建立虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系，所述等式关系为：
[0081][0082]
式中，矩阵b表示不含时变参数的控制效率矩阵，表示矩阵b的右伪逆矩阵；v表示矩阵b的零空间的一组正交基；k表示调整因子且在飞行初始时刻的值为0；
[0083]
期望控制输入解算模块，其配置成对虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系进行求解，获得四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni，从而解算出无人机的期望控制输入ud(αi,ni)。
[0084]
进一步地，所述期望控制输入解算模块中获得的四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni即共八个元素中，若超过无人机执行器物理极限的元素个数为一个或两个，则采用非线性重分配方法找到所述等式关系的一组精确解；具体包括：将虚拟推力向量n(αi,ni)分块为：其中，n
d1
对应推力的超限部分，n
d2
对应剩余部
分；同时，对正交基矩阵v也进行分块，分块规则与虚拟推力向量n(αi,ni)的分块规则相对应；则虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系变为：
[0085][0086]
其中，v1、v2分别代表对应于n
d1
、n
d2
的正交基分块矩阵；进而得到调整因子k的计算式为：k＝v
1-(n
d1-n1)，将调整因子k重新代入上述矩阵等式，获得重分配的虚拟推力向量：
[0087][0088]
最后经非线性变换得到执行器控制指令的精确解：
[0089][0090]
进一步地，所述期望控制输入解算模块中获得的四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni即共八个元素中，若超过无人机执行器物理极限的元素个数为两个以上，则采用二次规划法找到所述等式关系的一组接近于精确解的最优解：
[0091][0092]
式中，α
i0
,n
i0
分别表示上一时刻的各个旋翼绕各自倾转轴的倾转角和旋翼转速；
△
αi、
△
(n
i2
)分别表示各个旋翼绕各自倾转轴的倾转角变化率和旋翼转速变化率。
[0093]
进一步地，所述期望控制输入解算模块中采用二次规划法找到所述等式关系的一组接近于精确解的最优解的具体步骤包括：
[0094]
构造二次规划法所需的等式约束：
[0095][0096]
式中，x＝[
△
αi,
△
(n
i2
)]表示控制指令的变化率向量；d为高阶非线性项；
[0097]
构造二次规划法所需的目标函数：
[0098]
g(x,d,k)＝x
t
px+d
t
qd+k
t
rk
[0099]
式中，p,q,r分别是规模为8维、8维、2维的参数对角矩阵；
[0100]
构造二次规划法所需的极值约束：
[0101][0102]
式中，
△
α
max
表示倾转角变化率的极限值，n
max
表示转速的极限值；
[0103]
将上述等式约束、目标函数、极值约束输入二次规划求解器，获得当前时刻的
△
αi,
△
ni，以获得最优解αi和ni。
[0104]
本发明的有益技术效果是：
[0105]
本发明对于可倾转四旋翼无人机的动力学模型进行模型变换，得到含时变参数的控制效率矩阵，再利用变量代换将控制效率矩阵中的时变参数转移到虚拟推力向量中；利用矩阵伪逆与零空间的性质得到虚拟推力向量的显式表达；在虚拟推力向量中的期望指令
超限个数不大于2个时，采用实时性佳的非线性重分配方法得到执行器控制指令的精确解；在虚拟推力向量中的期望指令超限个数大于2个时，基于上一时刻控制指令的精确解进行二次规划，进而得到指令的变化量及当前时刻指令的最优解。本发明方法的实时性好，且所得解为精确解或最优解；适用于可倾转四旋翼飞行的任意阶段，且解算有效率在99％以上；本发明方法在无人机跟踪高机动及奇异轨迹时亦适用。
附图说明
[0106]
图1是本发明实施例中可倾转四旋翼无人机的结构示意图；
[0107]
图2是本发明实施例中可倾转四旋翼无人机的坐标系建立示意图；
[0108]
图3是本发明实施例中对比方法da方法的真实跟踪轨迹图；
[0109]
图4是本发明实施例中对比方法da方法的旋翼转速示意图；
[0110]
图5是本发明实施例中对比方法da方法的倾转角度示意图；
[0111]
图6是本发明实施例方法的真实跟踪轨迹图；
[0112]
图7是本发明实施例方法的旋翼转速示意图；
[0113]
图8是本发明实施例方法的倾转角度示意图；
[0114]
图9是本发明实施例一种基于零空间的可倾转四旋翼无人机控制分配系统的结构示意图。
具体实施方式
[0115]
为了使本技术领域的人员更好地理解本发明方案，在下文中将结合附图对本发明的示范性实施方式或实施例进行描述。显然，所描述的实施方式或实施例仅仅是本发明一部分的实施方式或实施例，而不是全部的。基于本发明中的实施方式或实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式或实施例，都应当属于本发明保护的范围。
[0116]
本发明首先利用坐标变换和动力学分析推导出可倾转四旋翼无人机的运动学模型和动力学模型；将模型变换得到含时变参数的控制效率矩阵，再利用变量代换将控制效率矩阵中的时变参数转移到虚拟推力向量中；利用矩阵伪逆与零空间的性质得到虚拟推力向量的显式表达；在虚拟推力向量中的期望指令超限个数不大于2个时，将超限部分与未超限部分所对应的矩阵分块，进而得到调整因子的计算式，再将调整因子重新代入分块矩阵等式，得到重分配的虚拟推力向量，最后经非线性变换得到执行器控制指令的精确解；在虚拟推力向量中的期望指令超限个数大于2个时，基于上一时刻控制指令的精确解进行二次规划，进而得到指令的变化量及当前时刻指令的最优解。
[0117]
本发明实施例提供一种基于零空间的可倾转四旋翼无人机控制分配方法，该方法包括以下步骤：
[0118]
步骤一、建立世界坐标系机体坐标系和旋翼坐标系其中，旋翼坐标系的原点固定于i号电机质心；i号旋翼绕倾转轴的倾转角为αi，相邻序号的旋翼倾转方向相同；4个旋翼的初始位置垂直于旋转平面向上，1、3号螺旋桨沿逆时针方向旋转，2、4号螺旋桨沿顺时针方向旋转；
[0119]
根据本发明实施例，如图1所示，可倾转四旋翼无人机包括机身1、飞行控制器2、
gps定位模块3、四个无刷电机4-1、4-2、4-3、4-4和四个倾转舵机5-1、5-2、5-3、5-4，倾转舵机设置在四旋翼无人机的机臂上，当收到电机转速指令和舵机偏转指令后，无刷电机的转子以指定速度转动，倾转舵机产生不同角度的偏转角，进而带动无刷电机围绕机臂整体转动。如图2所示，建立地面坐标系机体坐标系和旋翼坐标系旋翼坐标系的原点固定于i号电机质心，轴过i号电机所在机臂的轴心且正方向朝外，轴正方向垂直向上，轴符合右手螺旋定则。
[0120]
步骤二、根据牛顿-欧拉方程建立可倾转四旋翼无人机的动力学模型为：
[0121][0122]
其中，m表示无人机飞行器的质量；j＝diag(j
xx
,j
yy
,j
zz
)表示无人机飞行器的惯性矩阵(j
xx
,j
yy
,j
zz
分别表示机体在绕主轴xx、yy、zz旋转时的转动惯量)，表示无人机飞行器在机体坐标系下的质心线性加速度(a
x
；ay；az分别表示飞行器质心在机体系的x轴、y轴、z轴方向上的线性加速度)；表示无人机飞行器在机体坐标系下的角加速度(a
roll
；a
pitch
；a
yaw
分别表示飞行器在机体系下进行横滚、俯仰、偏航的角加速度)；g表示重力加速度；向量e3＝[0；0；1]；ω＝[p；q；r]表示无人机飞行器在机体坐标系下的角速度(p；q；r分别表示飞行器在机体系下进行横滚、俯仰、偏航的角速度)；f
drag
和m
drag
分别表示飞行器质心处受到的空气阻力及力矩；表示从机体坐标系到世界坐标系的旋转矩阵，i3表示3阶单位矩阵；uf(αi,ni)和um(αi,ni)分别表示飞行器质心处受到的合外推力及力矩，αi表示第i号旋翼绕倾转轴的倾转角，ni表示第i号旋翼的转速，i＝1、2、3、4。
[0123]
步骤三、进一步展开动力学模型中机体所受空气阻力f
drag
、空气阻力矩m
drag
、合外推力uf、合外推力矩um；
[0124]
根据本发明实施例，机体所受空气阻力f
drag
、空气阻力矩m
drag
与空气密度ρ
air
、飞行器迎风面积s
uav
、飞行线速度飞行角加速度有关
[0125][0126][0127]
其中，cf,cm≥0分别代表空气阻力及阻力矩系数。
[0128]
机体所受合外推力uf由4个螺旋桨产生的总推力组成，即
[0129][0130]
其中，kf≥0是旋翼推力系数，ni是i号旋翼的转速，λi是xb沿顺时针方向旋转到i号旋翼的夹角，且表示机体坐标系下螺旋桨推力的列矢
量(分别用s(
·
),c(
·
),t(
·
)表示sin(
·
),cos(
·
),tan(
·
))：
[0131][0132]
其中，r
x
(
·
),ry(
·
),rz(
·
)是右手系下的轴向旋转矩阵：
[0133][0134]
机体所受合外推力矩um由4个螺旋桨产生的总推力矩及其反扭力矩组成，即
[0135][0136]
其中，km≥0是旋翼反扭矩系数，表示机体坐标系下螺旋桨几何中心的坐标：
[0137][0138]
其中，l为ob到的距离。
[0139]
步骤四、将期望控制输入分解为含时变参数的控制效率矩阵左乘推力向量的形式；
[0140]
根据本发明实施例，在给定期望飞行轨迹的条件下，根据动力学模型可以获得期望控制输入ud(αi,ni)。引入控制效率矩阵γ(αi,ni)，用于将控制输入ud(αi,ni)映射到旋翼推力向量t(ni)：
[0141][0142]
其中，γ(αi)的行与空间自由度x,y,z,roll,pitch,yaw对应，列与螺旋桨编号i＝1,2,3,4对应。
[0143]
根据机体所受合外推力uf、合外推力矩um得到控制效率矩阵的显式表达式：
[0144][0145]
步骤五、将控制效率矩阵中的时变参数转移到推力向量中，此时推力向量变为虚拟推力向量，进一步建立虚拟推力向量与零空间中自由向量的等式关系；
[0146]
根据本发明实施例，利用变量替换法实现重映射，将旋翼推力t(ni)∈r4×1转换为虚拟推力n(αi,ni)∈r8×1，而变量形式的控制效率矩阵γ(αi)则变为标量形式的控制效率矩阵b∈r6×8。
[0147]
ud(αi,ni)＝bn(αi,ni)
[0148][0149]
由于b具有右伪逆故而有
[0150][0151]
又因为需利用零空间的性质将上式展开：
[0152][0153]
其中，i∈r8×8,v∈r8×2是b的零空间的一组正交基，k＝[k
1 k2]
t
是一组调整因子且在飞行的初始时刻k＝[0 0]。vk是位于b的零空间中的一组向量，且其在被b映射后为零向量，即b(vk)＝0。进而，虚拟推力向量与零空间中自由向量的等式关系为：
[0154][0155]
步骤六、执行器接收的控制指令αi,ni(i＝1,2,3,4)分别是4个舵机的倾转角和4个电机的转速。在一般情况下，该指令不能超过执行器的物理极限，否则容易导致执行器故障并使无人机失稳。然而，由于可倾转四旋翼无人机具备2个额外的自由度，所以其标量控制分配矩阵的零空间中具备2个正交基。进一步的，在虚拟推力向量中超限的元素个数μ不超过2个时，就可以采用实时性佳的非线性重分配方法找到一组精确解；
[0156]
根据本发明实施例，定义期望虚拟推力并将系统分块为其中n
d1
对应期望推力的超限部分，n
d2
对应剩余部分；同时，对正交基矩阵v也进行分块，且分块规则与虚拟推力矩阵的分块规则相对应；则原系统由变为：
[0157][0158]
其中，代表实际虚拟推力且v＝[v
1 v2]
t
，v1、v2分别代表对应于n
d1
、n
d2
的分块矩阵。由于虚拟推力向量的特殊形式，其内的超限元素总是成对出现，例如在n1》n
max
时，n(αi,ni)中的第1个元素和第二个元素均超出极限，此时v1∈r2×2，v2∈r6×2，且v1具备逆矩阵v
1-，v2具备左伪逆
[0159]
由上式可得k的计算式k＝v
1-(n
d1-n1)，进而实际虚拟推力可利用下式进行重分配：
[0160][0161]
由上述变换可知，当n
d1
超出约束时，调整因子k总能使超限部分重新满足约束，且重分配的结果仍为精确解。在得到n中所有元素的值之后，利用下式即可解算出当前时刻所需的电机转速αi与舵机倾转角ni。
[0162][0163]
式中n即为虚拟推力向量n(αi,ni)。
[0164]
步骤七、在虚拟推力向量中超限的元素个数μ超过2个时，由于正交基的个数不足以将指令重新分配至约束空间内部，所以采用优化方法中实时性较好的二次规划法找到一组接近于精确解的最优解。
[0165]
根据本发明实施例，采用优化方法中实时性较好的二次规划法找到一组接近于精确解的最优解的具体过程包括：
[0166]
步骤七一、构造二次规划法所需的等式约束。由于无人机初始动作一般是平稳起飞指令，不会导致执行器控制指令超限，因此初始阶段的控制指令一般是由步骤五及步骤六计算得到。在二次规划法被触发后，定义上一时刻的执行器控制指令为α
i0
,n
i0
(已在上一时刻被计算出)。进一步，利用雅可比矩阵对非线性等式进行线性展开，得到该二次规划问题的等式约束如下式所示：
[0167][0168]
其中，x＝[
△
αi,
△
(n
i2
)]表示控制指令的变化率向量，d为高阶非线性项。
[0169]
步骤七二、构造二次规划法所需的目标函数。在本发明中，目标函数的规划对象为执行器指令状态、高阶非线性项和零空间参数向量：
[0170]
g(x,d,k)＝x
t
px+d
t
qd+k
t
rk
[0171]
其中，p,q,r分别是规模为8维、8维、2维的参数对角矩阵。
[0172]
步骤七三、构造二次规划法所需的极值约束。根据执行器指令的约束空间π＝{
△
αi≤
△
α
max
,ni∈[0,n
max
]}确定执行器指令状态x＝[
△
αi,
△
(n
i2
)]的极限值，其中
△
α
max
表示舵机倾转角变化率的极限值，n
max
表示电机转速的极限值：
[0173][0174]
步骤七四、将等式约束、目标函数、极值约束输入二次规划求解器，得到x当前时刻的取值
△
αi,
△
ni，最后经下式得到当前时刻所需的电机转速αi与舵机倾转角ni：
[0175][0176]
本发明所述的基于零空间的控制分配算法伪代码如下所示。
[0177][0178]
进一步通过实验验证本发明的技术效果。
[0179]
飞行控制器的处理器为arm系列stm32h7型号单片机，处理频率达到500hz，满足控制分配解算的要求。将本发明方法与直接分配法(da)进行了对比。在假设无外部扰动的情况下，进行轨迹跟踪与姿态变换实验，验证了本发明方法的效果。预设参数如表1所示。
[0180]
表1 仿真实例相关参数
[0181][0182][0183]
所跟踪的期望轨迹如下：
[0184]
xd＝0,0《t≤30
[0185]
yd＝0,0《t≤30
[0186][0187]
[0188][0189]
ψd＝0,0《t≤30
[0190]
da方法的真实跟踪轨迹如图3所示，旋翼转速及倾转角度分别如图4、图5所示；本发明方法的真实跟踪轨迹如图6所示，旋翼转速及倾转角度如图7、8所示。从图中可以看出，da方法在跟踪期望轨迹时出现了发散，而本发明方法使无人机平稳地完成了轨迹跟踪且跟踪误差在2％误差带以内。
[0191]
本发明另一实施例提供一种基于零空间的可倾转四旋翼无人机控制分配系统，如图9所示，该系统包括：
[0192]
模型建立模块10，其配置成根据牛顿-欧拉方程建立可倾转四旋翼无人机的动力学模型；所述动力学模型表示为：
[0193][0194]
其中，m表示无人机飞行器的质量；j表示无人机飞行器的惯性矩阵；表示无人机飞行器在机体坐标系下的质心线性加速度；表示无人机飞行器在机体坐标系下的角加速度；g表示重力加速度；向量e3＝[0；0；1]；ω表示无人机飞行器在机体坐标系下的角速度；f
drag
和m
drag
分别表示飞行器质心处受到的空气阻力及力矩；表示从机体坐标系到世界坐标系的旋转矩阵，i3表示3阶单位矩阵；uf(αi,ni)和um(αi,ni)分别表示飞行器质心处受到的合外推力及力矩，αi表示第i号旋翼绕倾转轴的倾转角，ni表示第i号旋翼的转速，i＝1、2、3、4；其中，
[0195]
飞行器质心处受到的空气阻力表示为：
[0196][0197]
其中，cf≥0表示空气阻力系数；ρ
air
表示空气密度；s
uav
表示飞行器迎风面积；表示飞行器在机体坐标系下的质心线速度；
[0198]
飞行器质心处受到的空气阻力矩m
drag
表示为：
[0199][0200]
其中，cm≥0表示空气阻力矩系数；
[0201]
飞行器质心处受到的合外推力uf(αi,ni)表示为：
[0202][0203]
其中，kf≥0表示旋翼推力系数；表示机体坐标系下螺旋桨推力的列矢量，
[0204][0205]
λi表示机体坐标系的x轴沿顺时针方向旋转到第i号旋翼的夹角，且：
[0206]
飞行器质心处受到的合外推力矩表示为：
[0207][0208]
其中，km≥0表示旋翼反扭矩系数；表示机体坐标系下螺旋桨几何中心的坐标，
[0209][0210]
其中，l表示机体坐标系原点到旋翼坐标系原点的距离；
[0211]
期望控制输入变换模块20，其配置成对于给定的无人机期望飞行轨迹，根据所述动力学模型获得的期望控制输入ud(αi,ni)为：
[0212][0213]
将期望控制输入ud(αi,ni)分解为：含时变参数的控制效率矩阵γ(αi)左乘旋翼推力向量t(ni)，即：
[0214]
ud(αi,ni)＝γ(αi)t(ni)；
[0215]
将控制效率矩阵γ(αi)中的时变参数αi转移到旋翼推力向量t(ni)中，以获得虚拟推力向量n(αi,ni)，并建立虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系，所述等式关系为：
[0216][0217]
式中，矩阵b表示不含时变参数的控制效率矩阵，表示矩阵b的右伪逆矩阵；v表示矩阵b的零空间的一组正交基；k表示调整因子且在飞行初始时刻的值为0；
[0218]
期望控制输入解算模块30，其配置成对虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系进行求解，获得四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni，从而解算出无人机的期望控制输入ud(αi,ni)。
[0219]
本实施例中，可选地，所述期望控制输入解算模块中获得的四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni即共八个元素中，若超过无人机执行器物理极限的元素个数为一个或两个，则采用非线性重分配方法找到所述等式关系的一组精确解；具体包括：将虚拟推力向量n(αi,ni)分块为：其中，n
d1
对应推力的超限部分，n
d2
对应剩余部分；同时，对正交基矩阵v也进行分块，分块规则与虚拟推力向量n(αi,ni)的分块规则相对应；则虚拟推力向量n(αi,ni)与零空间中自由向量的等式关系变为：
[0220][0221]
其中，v1、v2分别代表对应于n
d1
、n
d2
的正交基分块矩阵；进而得到调整因子k的计算式为：k＝v
1-(n
d1-n1)，将调整因子k重新代入上述矩阵等式，获得重分配的虚拟推力向量：
[0222][0223]
最后经非线性变换得到执行器控制指令的精确解：
[0224][0225]
本实施例中，可选地，所述期望控制输入解算模块中获得的四个旋翼绕各自倾转轴的倾转角αi和旋翼的转速ni即共八个元素中，若超过无人机执行器物理极限的元素个数为两个以上，则采用二次规划法找到所述等式关系的一组接近于精确解的最优解：
[0226][0227]
式中，α
i0
,n
i0
分别表示上一时刻的各个旋翼绕各自倾转轴的倾转角和旋翼转速；
△
αi、
△
(n
i2
)分别表示各个旋翼绕各自倾转轴的倾转角变化率和旋翼转速变化率。
[0228]
本实施例中，可选地，所述期望控制输入解算模块中采用二次规划法找到所述等式关系的一组接近于精确解的最优解的具体步骤包括：
[0229]
构造二次规划法所需的等式约束：
[0230][0231]
式中，x＝[
△
αi,
△
(n
i2
)]表示控制指令的变化率向量；d为高阶非线性项；
[0232]
构造二次规划法所需的目标函数：
[0233]
g(x,d,k)＝x
t
px+d
t
qd+k
t
rk
[0234]
式中，p,q,r分别是规模为8维、8维、2维的参数对角矩阵；
[0235]
构造二次规划法所需的极值约束：
[0236][0237]
式中，
△
α
max
表示倾转角变化率的极限值，n
max
表示转速的极限值；
[0238]
将上述等式约束、目标函数、极值约束输入二次规划求解器，获得当前时刻的
△
αi,
△
ni，以获得最优解αi和ni。
[0239]
本实施例所述一种基于零空间的可倾转四旋翼无人机控制分配系统的功能可以由前述一种基于零空间的可倾转四旋翼无人机控制分配方法说明，因此本实施例未详述部分，可参见以上方法实施例，在此不再赘述。
[0240]
尽管根据有限数量的实施例描述了本发明，但是受益于上面的描述，本技术领域内的技术人员明白，在由此描述的本发明的范围内，可以设想其它实施例。对于本发明的范围，对本发明所做的公开是说明性的，而非限制性的，本发明的范围由所附权利要求书限定。