由于栖落机动给固定翼飞行器带来的好处，栖落机动的研究得到了越来越多的关注。文献[4]为了验证栖落的可行性，进行了小型滑翔机栖落在电线上的实验。文献[5]研究了变体部件对无人机栖落的影响。文献[6]通过配点法生成了栖落机动的标称轨迹。文献[7]对轨迹的稳定性研究发现，在终端阶段会发生发散，对栖落机动设计了滑模轨迹跟踪控制器。文献[8]对侧滑栖落机动提出了一种自适应增益滑模控制技术。文献[9]采用开放时间最优化方法来研究栖落轨迹优化问题。
目前，栖落机动的研究主要验证了栖落机动的可行性，以及如何设计控制器提高栖落机动的稳定性，没有考虑所设计的栖落机动系统的吸引域问题。吸引域是指非线性系统局部稳定的区域，为了保证飞行器能栖落在目标区域，需要计算栖落机动的吸引域。文献[10]提出了用平方和(Sum-of-Squares，SOS)算法计算动力学系统的吸引域。文献[11]针对栖落机动设计了误差反馈控制律，并将平方和算法用于计算栖落轨迹控制系统的吸引域。受文献[11]启发，本文针对栖落机动设计了更易于求解的分段线性控制律，并对平方和算法进行优化设计，以扩大吸引域。
本文采用广义伪谱法对栖落机动轨迹进行优化并生成标称轨迹，将动力学模型沿标称轨迹线性化，设计了分段线性的轨迹跟踪控制器；采用平方和算法计算系统的吸引域，并且优化吸引域算法以寻找出更大的栖落机动吸引域。
1 飞行器栖落机动的动力学模型 动力学模型的状态变量选取飞行速度V、航迹角μ、迎角α、俯仰角速率q、俯仰角θ、水平位移x和高度h，控制输入选取推力T和升降舵的偏转角度δe。假设推力经过重心且指向机头方向，图 1为飞行器的纵向受力分析图。图中：Xa为速度坐标系的横向位移；Xb为机体坐标系的横向位移；Xg为地面坐标系的横向位移；Za为速度坐标系的纵向位移；Zb为机体坐标系的纵向位移；Zg为地面坐标系的纵向位移；g为重力系数。
 |
| 图 1 飞行器纵向受力分析图 Fig. 1 Longitudinal force analysis diagram of UAV |
| 图选项 |
根据图 1建立速度坐标系下的飞行器纵向动力学方程:
 | (1) |
式中：m为飞行器的质量；M为空气动力和力矩；L为升力；D为阻力；Iy为飞行器俯仰转动惯量。
L、D、M的表达式如下：
 | (2) |
式中：CL、CD和CM分别为升力、阻力和力矩系数；Sa为飞行器的空气动力表面积；ρa为空气密度。
飞行器栖落机动的空气动力学系数由平板模型方法[12-14]得到，可表达为
 | (3) |
式中：le为升降舵空气动力重心到飞行器质心距离；Se为升降舵的空气动力表面积。
栖落机动的动力学模型与常规飞行器的区别主要在于栖落机动具有大迎角。栖落属于一类复杂的机动，是非线性领域控制的一个挑战。
2 栖落机动轨迹跟踪控制律设计 2.1 非线性参考轨迹 将动力学模型(1)记为
 | (4) |
式中：X为非线性模型的状态向量；u为非线性模型的输入量；f(X, u)表示非线性函数。
针对轨迹优化和非线性轨迹，采用广义伪谱法为基础的数值非线性优化程序生成[15-16]。目前的广义伪谱法主要有4种类型，分别为Guass伪谱法、Legendra伪谱法、Rudau伪谱法和Chebyshev伪谱法。本文采用Rudau伪谱法进行轨迹优化，主要原因是：Rudau伪谱法相对其他3种广义伪谱法有着较高的精度，在其他3种广义伪谱法的基础上增加了初始点或者着重点。
动力学模型含有6个状态量，在进行优化计算时，将升降舵的偏转角度δe作为新的状态量，则状态向量可表示为X=[V, α, μ, q, x, h, δe]T，并且设定参考的推力T为常值3.768 N。广义伪谱法是根据整个过程对轨迹状态量的约束和状态方程生成标称轨迹，所以需要对7个状态量根据栖落机动轨迹的要求进行约束。在执行栖落机动的过程中设置参数: 将初始点位置设置为坐标原点，并且根据实际情况设置各个状态的约束范围。升降舵的偏转角度设为状态量，所以将其导数作为输入量u，给定u的约束条件为：τ(u)=(-1, 1)。设定初始时刻上限值和下限值为：Xh(t0)=[9.984 m/s, 0 rad, 0.25 rad, 0 rad/s, 1 m, 0 m, -0.15]T，Xl(t0)=[9.984 m/s, 0 rad, 0.25 rad, 0 rad/s, 0 m, -1 m, -0.15]T。设置状态量过程中的范围如表 1所示。设定终点时刻上限值和下限值为Xh(tf)=[4 m/s, π/2 rad, π/4 rad, 3.5 rad/s, 15 m, 2 m, π/3]T，Xl(tf)=[3 m/s, -π/2 rad, -π/4 rad, -3.5 rad/s, 14 m, 1 m, -π/3]T。栖落终点的位置要根据实际的栖落点确定，如栖落在电线杆上，所以设置栖落终点时给定了一定的约束范围，如表 1所示。
表 1 状态变量过程约束 Table 1 Process constraints of state variables
| 状态变量 | 下限值 | 上限值 |
| V/(m·s-1) | 0 | 25 |
| α/rad | -π/2 | π/2 |
| μ/rad | -π/4 | π/4 |
| q/(rad·s-1) | -3.5 | 3.5 |
| x/m | 0 | 15 |
| h/m | 0 | 10 |
| δe/rad | -π/3 | π/3 |
表选项
选取的二次型优化指标J为
 | (5) |
式中：第1项为积分型，是对整个过程中产生的状态误差和输入误差进行的补偿；第2项为终点型，是对终点产生的误差进行补偿。二次性能指标(5)体现为对栖落机动的综合性能优化。式(5)中：Qf、Q、R分别为对应项所占的权重；x(tf)为栖落终点的目标状态，所以对Qf的选择越大，越要求栖落终点位置越接近指定的位置；u为整个过程的输入，R越大则对输入整体的操控需要的能量越小；x为整个过程的状态变量，Q越大则对整个状态变量要求越平稳。栖落机动状态量变化的幅度较大，所以对状态变量不设限制，整个优化可以理解为以最小的输入达到最优的最终栖落状态。其中，Qf=diag[10, 30, 0, 0, 10, 10, 0], R=90, Q=diag[0, 0, 0, 0, 0, 0, 0]。
飞行器的物理参数来源于文献[14], 具体参数如表 2所示。
表 2 飞行器物理参数 Table 2 Physical parameters of UAV
| 参数 | 数值 |
| 质量m/kg | 0.8 |
| 平均气动弦长c/m | 0.25 |
| 展长b/m | 1 |
| 升力面积Sl/m2 | 0.25 |
| 俯仰转动惯量Iy/(kg·m-2) | 0.1 |
表选项
根据以上约束条件用广义伪谱法即对应MATLAB的GPOPS工具包针对非线性模型生成的标称轨迹如图 2和图 3所示。
 |
| 图 2 状态变量标称曲线 Fig. 2 Nominal curves of state variables |
| 图选项 |
 |
| 图 3 推力和升降舵的偏转角度标称曲线 Fig. 3 Nominal curves of thrust and rudder angle |
| 图选项 |
从图 2中看出，经过轨迹优化设计出栖落机动轨迹的高度和水平位置都可以很好地落在终点约束的目标范围内。
将生成的标称轨迹记为(Xr, ur)。标称轨迹的状态方程为?r=f(Xr, ur)。其中，Xr为GPOPS得到的标称轨迹的状态变量，ur为标称轨迹的输入。将非线性模型(4)沿着标称轨迹(Xr, ur)泰勒展开：
 | (6) |
式中：ΔX和Δu分别为实际和标称轨迹状态量和输入量的差值；?(t)、
 | (7) |
为了将非线性模型分段线性化，定义飞行器栖落机动的时间范围为[t0, tf]，在时间范围上均匀地选取n个时间点{t0, t1, …, tn-1}，且有tn-1=tf。在每个时刻对标称轨迹线性化，所以任意时刻点tq的线性化模型为
 | (8) |
式中：?q和
整个时间内分段线性模型为
 | (9) |
式中：χq(t)相当于随着时间变换的函数，可描述为
 | (10) |
2.2 控制律设计 本节针对栖落机动分段线性模型(9)设计采用分段线性最优轨迹跟踪控制器。针对飞行器栖落机动动力学系统，寻找一个最优控制uc(t)使得X能跟踪上Xr。选取的二次型性能指标J1如下：
 | (11) |
式中：Q≥0, R>0为正定的对角阵常值矩阵；ΔXq为实际轨迹和标称轨迹在q时刻点的状态差值；Δuq为q时刻点的最优控制，需要弥补到实际输入上。
根据线性二次型调节器(Linear Quadratic Regulator，LQR)可得到控制律kq为
 | (12) |
 | (13) |
式中：Sq为满足式(14)的正定对称矩阵，即黎卡提方程：
 | (14) |
所以，系统在q时刻的闭环状态方程为
 | (15) |
式(15)设置了分段线性系统的每个子系统的最优控制律，可以使子系统稳定地收敛至标称轨迹附近。
虽然上述每一个分段线性系统可以保证稳定性，但是不能保证整个栖落系统的稳定性。为了保证整个栖落系统的稳定，需要使得每一个子系统的Lyapunov函数值大于或等于下一个子系统的Lyapunov函数值，即满足：
 | (16) |
式中：VL为Lyapunov函数；Pq和Qq为选取的正定对称矩阵。
下面对式(16)进行证明。
引理1[17]??时变系统?=C(t)x，x(t0)=x0，且t≥t0，C(t)为以时间t的连续和分段连续函数为元的矩阵。则系统一致渐近稳定的充要条件是：存在N>0, c>0，使得对于任意的t0和t≥t0，状态转移矩阵Φ(t, t0)满足：
 | (17) |
证明??构造Lyapunov函数：
 | (18) |
选取时间为(tq-Δt, tq+Δt]的线性子系统。因为Pq为正定对称矩阵，所以将其分解为Pq=UTdiag(λ1, λ2, …, λn)U，λq为Pq的特征值，U为酉矩阵。因此，子系统满足：
 | (19) |
由式(15)闭环子系统可推导出
 | (20) |
再由式(16)可得
 | (21) |
因为
 | (22) |
由式(19)得
 | (23) |
联立式(21)~式(23)得到
 | (24) |
引入非线性系统比较原理，可以证明Ve满足：
 | (25) |
栖落机动的整个过程分为N个子空间且经历的时间为[t0, tf], 所以由引理1可以得到Ve在(tq-Δt, tq+Δt]时间按指数收敛。
又由式(16)可知，在2个子系统切换时刻ts=tq-Δt满足：
 | (26) |
所以，综合可知Ve按照指数形式收敛且在切换时刻后非增，整个系统渐近稳定。由式(13)计算控制律kq后要代到式(16)验算是否满足整个系统的渐近稳定要求，如果不满足，要调整式(11)中Q、R的权值再次验算。这样，满足条件(16)的控制律(13)能保证整个栖落分段线性系统(9)的稳定。????证毕
轨迹跟踪控制系统的结构如图 4所示。
 |
| 图 4 轨迹跟踪控制框图 Fig. 4 Trajectory tracking control block diagram |
| 图选项 |
3 栖落机动轨迹跟踪控制吸引域的计算 吸引域是指非线性系统局部稳定的区域。栖落机动轨迹的吸引域能保证其内的飞行器能够在规定时间内栖落在指定目标区域。2.2节中设计的栖落机动控制律是针对分段线性化的模型设计的。而栖落机动的实际模型是非线性动力学系统(1)。针对非线性系统(1)计算其在控制律(13)下的局部稳定范围即栖落机动轨迹的吸引域是一个很重要的问题。
本节运用平方和算法[18]求出闭环非线性栖落系统标称轨迹的吸引域区域Ω[19]。在求出吸引域后，通过优化平方和算法[20]中求解的变量，对原有吸引域的半径进行了很大程度的扩大，使得能寻找出更大的吸引域。3.1节将介绍在一般求解吸引域方法的基础上求解栖落非线性系统的吸引域。3.2节在3.1节吸引域的基础上通过优化平方和算法的求解变量扩大吸引域。
3.1 轨迹吸引域算法 在非线性时变的情况下，将系统写成如下形式：
 | (27) |
式中：f为关于x和t的多项式函数。设Γ?[t0, tf]×Rn为系统解的集合，给定的目标区域记为g?Rn。当集合Γ终点时刻的状态收敛到g，也就是对于任何的x∈Rn，当(tf, x)∈Γ时，x∈g，即为系统整个运动过程的吸引域。
引理2[10]??若存在一个连续可微的函数: [t0, tf]×Rn→[0，∞)。对于每个t∈[t0, tf]定义：Ωt：={x|V(x, t)≤1}和?Ωt: ={V(x, t)=1}，如果每任意t∈[t0, tf]和x∈?Ωt都满足：
将栖落机动容许的落点范围定义为目标区域g，定义栖落机动标称轨迹xr的候选吸引域区域为
 | (28) |
式中：X为实际和标称轨迹的差量；V(X, t)为正定的；ρ(t)为时变大于零的标量；Ω(ρ，t)为吸引域区域。选择二次型V(X, t)=XTSX，S为2.2节用LQR求出的常值矩阵。
根据引理2，要使式(28)为吸引域，要满足ρ(t)≤1和式(29)：
 | (29) |
式(29)可由引理2将V(x, t)=V(X, t)/ρ(t)推导可得。由式(29)可知，当t=tq-1即当V(X, tq-1)=ρ(tq-1)时，在下一个时刻t=tq时满足：
 | (30) |
即当t=tq时状态量仍在吸引域内。
选择适当的ρ(t)使得Ω(ρ，t)满足：
 | (31) |
将上述求吸引域的过程整理成最优求解的形式如下：
 | (32) |
将式(32)表示的条件转化为
 | (33) |
最终将条件转化为
 | (34) |
引理3[21]??给定多项式集{f1…fs}，{g1…gl}，{h1…hu}∈Rn，下面的2个条件是等价的:
1) 集合
2) 存在多项式f∈p{f1…fs}, g∈M{g1…gl}, h∈I{h1…hu}满足f+g2+h=0，Rn表示所有实数域上n元多项式的集合。
再根据引理3，将条件(34)转化为
 | (35) |
式中：s1，s2，s3∈Σn，为了简化使用平方和算法的计算，取s3=1并将s1移项。平方和算法要求未知量之间是线性关系，所以将s2转化s2(X，t)。将ρ(t)按照一定的步长增加，通过求出s2(X，t)有解找出最大的ρ(t)，所以式(35)转化为
 | (36) |
s2(X, t)设定格式如下：
 |
式中：O为6×6的实对称矩阵。
将栖落的时间[t0, tf]分为N段，计算栖落机动标称轨迹的吸引域。根据栖落终点的误差范围给定终点吸引域半径ρf的值。假设在[tf-1+Δt, tf]这段时间内的吸引半径都是ρf，计算tf-1时刻的吸引域半径ρf-1，这时用ρf-1按照线性的模式去估计时间[tf-1-Δt, tf-1+Δt]内的半径。
在计算整个标称轨迹吸引域时，按照t=tq的时间取样。因为ρf是根据栖落终点范围约束给定的，所以从tf时刻开始倒着往前计算。由于式(36)左边是对ρ(t)和s2(X，t)的双线性函数，通过不断递增当前时刻ρ(tq)，直到s2(X，t)无解，即是最大的吸引域半径。
3.2 扩大吸引域半径的优化算法 3.1节推导中V(X, t)=XTSX，S为2.2节求解黎卡提方程(14)得到的正定常值矩阵。之所以选择LQR里的S矩阵为Lyapunov函数的系数矩阵，是因为这样能保证系统是渐近稳定的，即
将式(12)代入模型系统(6)，得到的系统闭环状态方程为
 | (37) |
 | (38) |
 | (39) |
 | (40) |
由式(14)可知
 | (41) |
 | (42) |
因为Q≥0, S>0, R>0，所以
 | (43) |
 | (44) |
在式(36)中，V(X, t)中的S和吸引域半径ρ(tq)必须固定一个才能求解。3.1节求解吸引域半径时是固定了这个时刻的S矩阵寻找最大吸引域半径ρ(tq)。
为了寻找出更大的吸引域半径，先固定吸引域的初始半径，按照3.1节的方法求解下一时刻的吸引域半径ρ(tq)，求出最大的ρ(tq)且保存下此时刻的s2(X，t)。再固定求出的最大的ρ(tq)和s2(X，t)，通过按设定的方法不断改变S矩阵以扩大吸引域。但是，改变S矩阵时一定要保证
 | (45) |
直到式(45)左边不属于平方和时，选取无解时前一时刻有解的S矩阵。
改变S矩阵的方法是给S矩阵乘上一个常数k(0 < k < 1), k按照一定步长从1往下递减直到kS矩阵使得式(45)不成立。这样选取改变S矩阵的原因有：
1) 乘上一个0k < 1的常数k，没有改变S特征值的符号即可以保证V(X, t)正定和
证明如下：
 | (46) |
因为只是选取的S矩阵发生了变化，而由LQR控制器得到的状态方程没有变化，所以
 | (47) |
 | (48) |
 | (49) |
 | (50) |
因此，0 < k < 1时，
2) 选择一定步长从1往下递减的原因是：在最优化S矩阵的过程中，搜索出的最优S矩阵可以变相地扩大吸引域半径。
推导过程如下：
 | (51) |
当S矩阵乘上k
 | (52) |
两边同时除以k
 | (53) |
又因为0 < k < 1，所以ρ(t)/k比原来的ρ(t)大，所以减小S矩阵可以扩大吸引域半径。
4 仿真结果与分析 4.1 栖落机动轨迹跟踪控制仿真 固定翼飞行器的几何参数如表 2所示，仿真采用栖落机动非线性模型(1)，非线性模型的标称输入为由GPOPS得到的标称轨迹的推力和升降舵的偏转角度。气动参数由式(2)、式(3)计算。
LQR控制器采用分段线性化的形式设计，整个栖落的时间设定为2 s，每隔0.05 s选取出参考点进行控制律设计。R=diag[9, 45]，Q=diag[18, 100, 16, 100, 20, 50]; 联立式(16)求得控制律。将求得的控制律加入到非线性模型的输入上，即可得到经LQR控制器的跟踪轨迹。飞行器初始理想状态X*(t0)=[9.984 m/s, 0 rad, 0.25 rad, 0 rad/s, 0 m, 0 m], 第2节的标称轨迹的初始状态也是按照理想状态设定，得到了理想的标称轨迹。为了检验跟踪控制器的效果，在设计实际的初始状态上加上了偏差。加上偏差的初始状态为[11 m/s，0.1 rad，0.4 rad，0 rad/s，-1 m，-0.7 m]。仿真结果如图 5和图 6所示。
 |
| 图 5 控制输入曲线 Fig. 5 Control input curves |
| 图选项 |
 |
| 图 6 状态变量跟踪控制曲线 Fig. 6 State variable tracking control curves |
| 图选项 |
由图 5和图 6可知，在LQR控制器的作用下，由非线性模型生成闭环曲线的状态量均能跟上参考轨迹。所以，设计的LQR控制器是有效的。控制输入范围分别在T∈[2.7, 4.4]N, δe∈[-0.7, 0.3]rad，是合理的输入量范围。栖落的最终位置精度要求很高。在飞行器栖落机动的跟踪曲线中，由GPOPS生成的参考轨迹的栖落位置是(14.80，1.162) m，在控制器下的飞行器实际栖落位置是(14.79，1.160) m。由此可见，即使初始位置偏差较大，经过LQR控制器的作用，飞行器能准确地降落在参考轨迹栖落位置附近。
4.2 栖落机动轨迹吸引域和扩大吸引域仿真 基于非线性模型(1)，在栖落的时间2 s内每隔0.05 s选取一个时间点联立式(36)，计算每个时间段的吸引域。栖落机动吸引域是从栖落终点时刻开始倒推计算的。因为飞行器栖落机动要求降落在参考轨迹终点的附近，需要设定栖落位置允许的误差。栖落终点的容许位置误差设为Δx < 0.15 m, Δh < 0.15 m。根据栖落终点位置的容许误差范围，设定终点2 s的栖落收敛半径为0.15 m。已知2 s的半径，即可根据式(36)推出1.95 s时刻的吸引域半径，如此倒推出整个过程的吸引域半径。飞行器栖落机动的整个过程共计算了41个吸引域半径。
为了形象地展示出吸引域的图形，将吸引域半径投影到xh平面上，即根据吸引域的表达式V(X, t)≤ρ(t)，将V(X, t)=XTSX中的X除了水平位置和高度都设置为0，右边是吸引域半径，即可得到xh平面的投影。
根据式(36)得到的吸引域比较小。根据3.2节扩大吸引域的算法，通过减小S矩阵来扩大吸引域。图 7为根据计算的41个吸引域半径在xh平面的投影，其中绿色是没有扩大前的吸引域，黑色是扩大后的吸引域，红色是由GPOPS根据约束条件生成的参考轨迹水平位置和高度的曲线。
 |
| 图 7 吸引域扩大对比 Fig. 7 Comparison of expanded attraction domain |
| 图选项 |
由图 7可知，通过优化计算吸引域的方法，可以找出飞行器栖落机动更大的吸引域，这样使得飞行器能够在更大的初始状态偏差下保证栖落在预定目标范围。
为了验证吸引域的可靠性，选取了4个在0 s吸引域的不同初始状态，仿真结果如图 8所示，栖落轨迹全部在所求吸引域内，并且终点范围收敛到目标区域验证了吸引域的有效性。
 |
| 图 8 栖落轨迹及吸引域 Fig. 8 Perching trajectory and attraction domain |
| 图选项 |
5 结论 1) 广义伪谱法可以根据对飞行状态的约束生成栖落机动标称轨迹。
2) 栖落轨迹跟踪控制仿真的结果表明，采用LQR设计的轨迹跟踪控制器可以很好地跟踪上由GPOPS生成的标称轨迹，验证了LQR控制器的有效性。
3) 吸引域的仿真结果表明，通过优化求解吸引域的过程可以寻找出更大的栖落机动轨迹的吸引域。对扩大的吸引域进行的栖落机动控制仿真表明，起始于初始吸引域内的轨迹均保留在吸引域内，验证了其有效性。
本文是针对一条栖落轨迹进行控制器及吸引域的设计。后续工作中将设计多条轨迹形成栖落轨迹库，以此为基础研究吸引域计算，扩大飞行器实现定点栖落的容许初始范围，并逐步展开栖落机动实验验证的研究。后续科研工作中，也将进一步考虑参数不确定性和干扰条件对控制器的影响，以及如何提高控制器的鲁棒性。
参考文献
| [1] | POPE M, TKIMES C W, JIANG H, et al. A multimodal robot for perching and climbing on vertical outdoor surfaces[J]. IEEE Transactions on Robotics, 2017, 33(1): 38-48. DOI:10.1109/TRO.2016.2623346 |
| [2] | FEROSKHAN M, TIAUW H. Control strategy of sideslip perching maneuver under dynamic stall influence[J]. Aerospace Science and Technology, 2018, 72: 150-163. DOI:10.1016/j.ast.2017.11.002 |
| [3] | COLIN G, ANTONY W, THOMS R. Perched landing manoeuvres with a variable sweep wing UAV[J]. Aerospace Science and Technology, 2017, 71: 510-520. DOI:10.1016/j.ast.2017.09.034 |
| [4] | CORY R, TEDRAKE R.Experiments in fixedwing UAV perching[C]//Proceedings of the 33th Chinese Control Conference, 2014: 2047-2056. |
| [5] | 袁亮, 何真, 王月. 变体无人机栖落机动建模与轨迹优化[J]. 南京航空航天大学学报, 2018, 50(2): 267-275. YUAN L, HE Z, WANG Y. Modeling and trajectory optimization of perching maneuvers for morphing UAV[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2018, 50(2): 267-275. (in Chinese) |
| [6] | ROBEST J W, CORY R, TEDRAKE R.On the controlability of fixedwing perching[C]//American Control Conference, 2009: 2018-2023. |
| [7] | VENKAESWARA R, TIAUW H. Optimization stability analysis and trajectory tracking of perching maneuvers[J]. Journal of Guidance, Control, and Dynamics, 2014, 37(7): 879-888. |
| [8] | MIR F, TIAUW H G. Performance of sideslip perching maneuver with an adaptive gain control feature[J]. Aerospace Science and Technology, 2018, 78: 648-660. DOI:10.1016/j.ast.2018.05.026 |
| [9] | MIHIR V, SHAWN K, GONZALO G, et al.Fixed wing aircraft perching[C]//AIAA Guidance Navigation and Control Conference.Reston: AIAA, 2017: 14-34. |
| [10] | MARK M, IAN R, RUSS T.Invariant funnels around trajectories using sum-of-squares programming[C]//The International Federation of Automatic Control, 2014: 9218-9223. |
| [11] | TEDRAKE R, MANCHESTE I R, TOBENKIN M, et al. LQR trees: Feedback motion planning via sums-of-square verification[J]. International Journal of Robotics Research, 2010, 29(8): 1025-1038. |
| [12] | PUPPOLO M, REYNOLDS R, JACOB J.Comparison of three aerodynamic model used in simulation of high angle of attack UAV perching maneuver[C]//Proceedings of the AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.Reston: AIAA, 2006: 1291-1302. |
| [13] | TAHK M J, HAN S, LEE B Y, et al.Trajectory optimization and control algorithm of longitudinal perch landing assisted by thruster[C]//Proceedings of the IEEE Control Conference.Piscataway: IEEE Press, 2017: 2247-2252. |
| [14] | 李达. 飞行器栖落机动飞行轨迹优化与控制[D]. 南京: 南京航空航天大学, 2017. LI D.Trajectory optimization and control of perching maneuvers for aircraft[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2017(in Chinese). |
| [15] | KAUSHIK H, MOHAN R, PRAKASH K A. Utilization of wind shear for powering unmanned aerial vehicles in surveillance application: A numerical optimization study[J]. Energy Procedia, 2016, 90: 349-359. DOI:10.1016/j.egypro.2016.11.202 |
| [16] | RAO A V, BENSON D A, DQRBY C, et al. Algorithm 902:GPOPS, a MATLAB software for solving multiple phase optimal control problems using the Gauss pseudospectral method[J]. ACM Transactions on Mathematical Software, 2010, 37(2): 163-172. |
| [17] | BARAK B, KENLNER J, STEURER D.Dictionary learning and tensor decomposition via the sum-of-squares method[C]//Forty-Seventh ACM Symposium on Theory of Computing, 2015: 143-151. |
| [18] | ZHOU Y, KUANG C, SUN Y.Nonlinear H-infinity control for spacecraft attitude maneuver based on SOS optimization[C]//13th International Conference on Computer Science & Education, 2018: 46-51. |
| [19] | 隋吉超, 罗飞. 基于平方和规划法的一种估计系统吸引域的改进算法[J]. 科学技术与工程, 2012, 12(5): 978-980. SUI J C, LUO F. An improved algorithm for estimating the domain of attraction of a system based on the square sum programming[J]. Science Technology and Engineering, 2012, 12(5): 978-980. DOI:10.3969/j.issn.1671-1815.2012.05.002 (in Chinese) |
| [20] | RADIAN F, YINGJEN C, MOTOYASU T. An SOS based control Lyapunov function design for polynomial fuzzy control of nonlinear systems[J]. IEEE Transactions on Fuzzy Systems, 2017, 25(4): 775-787. DOI:10.1109/TFUZZ.2016.2578339 |
| [21] | 李瑞莲. 一类非线性切换系统的吸引域估计及控制器设计[D]. 沈阳: 东北大学, 2009. LI R L.Domain of attraction estimation and control for a class of nonlinear switched systems[D].Shenyang: Northeastern University, 2009(in Chinese). |
