一种火箭子级姿态翻转着陆在线制导方法与流程

1.本发明属于飞行器制导技术领域，特别涉及一种火箭子级姿态翻转着陆在线制导方法。


背景技术：

2.随着航天技术的发展，重复使用已经成为未来航天运输系统的主要发展方向之一。目前，运载火箭一子级重复使用技术发展迅速，但对于更高势能火箭子级(如火箭末级)，回收难度更大，其飞行剖面与运载火箭一子级相比具有较大差别。传统的运载火箭一子级垂直动力着陆阶段可以不需要推力参与飞行器的姿态控制(如falcon9)，主要通过rcs(reaction control system，反应控制系统)和栅格舵来实现姿态的调整，因此大多数学者针对火箭一子级回收着陆轨迹优化模型通常是不考虑姿态的三自由度动力学模型，但由于重复使用火箭末级动力着陆初始时刻以大攻角气动减速状态下降，所以在最后着陆阶段需要通过发动机推力进行姿态翻转机动并完成垂直动力着陆，因此需要设计一种新的着陆制导方法。


技术实现要素：

3.为了克服现有技术中的不足，本发明人进行了锐意研究，提供了一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，在传统三自由度返回着陆制导方法的基础上引入姿态变量，建立了六自由度返回着陆问题，用于解决着陆过程中姿态翻转机动问题，同时，位置、速度、姿态满足终端约束要求，从而完成本发明。
4.本发明提供的技术方案如下：
5.目前，大多数动力着陆轨迹优化方法针对的模型都是类似于火箭一子级的三自由度(3degree of freedom，3dof)质心动力学模型，在着陆过程中并未考虑姿态的变化，这是因为火箭一子级垂直动力着陆阶段姿态变化较小，通常并不需要推力参与飞行器的姿态控制，主要通过rcs和栅格舵来实现姿态的调整。但由于“星舰”式火箭末级着陆段初始姿态处于大攻角气动减速状态，所以它在最后着陆阶段需要通过发动机推力进行姿态翻转机动并完成垂直动力着陆。因此，对于这种飞行器的垂直着陆轨迹优化需要考虑推力对姿态的影响。
6.针对“星舰”式火箭子级动力着陆过程采用推力进行姿态翻转机动的特点，提出了一种基于改进序列凸优化的六自由度轨迹优化方法，该方法引入姿态四元数和角速度变量建立了六自由度着陆轨迹优化问题，通过凸化处理和离散化方式将非凸的轨迹优化问题转化为一系列有限维的二阶锥子问题，设计虚拟控制和动态信赖域更新策略保证二阶锥子问题的可行性和收敛性，然后通过求解一系列二阶锥子问题逐渐收敛至最优轨迹，最后以六自由度动力着陆轨迹优化方法为基础建立火箭动力软着陆滚动时域mpc(模型预测控制，model predictive control)制导框架，通过仿真实验来验证建立的基于轨迹序列凸优化mpc制导方法的有效性。
7.具体方案如下
8.一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，包括如下步骤：
9.s1，引入姿态四元数和转动角速度描述火箭末级的姿态运动，建立火箭末级六自由度着陆动力学模型；
10.s2，以最短着陆时间为目标函数，构建满足约束条件的火箭末级六自由度动力软着陆轨迹优化模型；
11.s3，将步骤s2轨迹优化模型中的非凸约束通过线性化方式转化为凸约束，得到凸形式的轨迹优化模型，所述轨迹优化模型的非凸约束包括动力学方程和推力大小约束；
12.s4，采用一阶保持器的方法将步骤s3中凸形式的轨迹优化模型进行离散化处理，得到离散凸化模型；
13.s5，对步骤s3中线性化动力学方程添加动力学松弛变量，以罚函数的形式将松弛变量添加到目标函数中；
14.s6，设计信赖域约束限制参考轨迹的变化范围；
15.s7，通过初末端状态约束简单线性插值得到初始迭代参考轨迹；
16.s8，步骤s7中的初始迭代参考轨迹作为序列凸化的迭代初值，考虑步骤s5，s6的约束，求解步骤s4中得到的离散凸化模型；
17.s9，重复步骤s8迭代求解使轨迹收敛到最优轨迹，完成一个制导周期采样点的轨迹优化；
18.s10，将步骤s9的一个制导周期采样点的轨迹优化结果作为制导指令生成器，根据当前状态更新执行步骤s2-s9，将轨迹优化结果更新最优指令，并直接用作制导信号，最终完成火箭子级姿态翻转着陆在线制导。
19.根据本发明提供的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，具有以下有益效果：
20.(1)本发明针对火箭末级姿态翻转着陆的特点，在传统三自由度动力着陆模型下引入姿态变量建立六自由度动力着陆轨迹优化问题，提升了任务适应性。
21.(2)本发明对传统固定信赖域的序列凸化方法进行改进，设计了基于线性化误差度量的动态信赖域策略，提升了轨迹优化方法的求解效率。
22.(3)本发明建立了以六自由度动力着陆轨迹优化方法为基础的模型预测制导框架，成功将在线轨迹优化方法运用于制导控制系统中，实现火箭子级姿态翻转着陆全程在线制导。
附图说明
23.图1为本发明的动力着陆示意图；
24.图2为本发明的制导方法框架结构图；
25.图3为本发明的制导方法仿真结果图，其中图3a为6dof三维轨迹图；图3b为6dof纵向平面轨迹投影图；图3c为位置变化曲线图；图3d为速度分量变化曲线图；图3e为推力大小变化曲线图；图3f为质量变化曲线图；图3g为速度大小变化曲线图；图3h为角速度分量变化曲线图。
具体实施方式
26.下面通过对本发明进行详细说明，本发明的特点和优点将随着这些说明而变得更为清楚、明确。
27.在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。
28.本发明提供了一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，包括如下步骤：
29.s1，引入姿态四元数和转动角速度描述火箭末级的姿态运动，建立火箭末级六自由度着陆动力学模型。
30.以目标点地面坐标系f
l
和本体坐标系fb为参考坐标系，描述飞行器与目标点相对位置与姿态关系，如图1所示。考虑火箭末级在最后动力着陆阶段需要通过发动机推力进行一次姿态翻转机动以实现垂直动力着陆，所以在三自由度着陆动力学模型的基础上，引入姿态四元数(q
b/l
)以及转动角速度(ωb)以描述箭体的姿态运动，得到的六自由度着陆动力学模型的具体形式如下：
[0031][0032]
其中，m为飞行器的质量；αm＝1/(i
sp
g0)，i
sp
为比冲，g0为标准地球重力加速度；t为运行时刻，r
l
、v
l
为飞行器在目标点地面坐标系f
l
下的位置、速度矢量；d
l
＝-0.5ρscd||v
l
||v
l
为在目标点地面坐标系f
l
下的气动阻力，ρ大气密度，s为参考面积，cd为阻力系数；tb为在本体坐标系fb下的推力；ωb为在本体坐标系fb下本体坐标系fb相对目标点地面坐标系f
l
的角速度矢量；q
b/l
为目标点地面坐标系f
l
相对本体坐标系fb的单位姿态四元数；g
l
为重力矢量；r
t,b
为推力作用点到质心力矩的力臂；jb为惯性矩阵；c
t/b
(t)为本体坐标系fb到目标点地面坐标系f
l
的旋转矩阵，展开形式为：
[0033][0034]
其中，q0，q1，q2，q3为对应的姿态四元数；
[0035]
ω(ωb)为四元数的乘法展开形式：
[0036]
[0037]
其中，ωb＝[ω
bx ω
by ω
bz
]
[0038]
[ζ
×
]为叉乘展开形式：
[0039][0040]
其中，ζ＝[ζ
x ζ
x ζ
x
]为变量。
[0041]
s2，以最短着陆时间为目标函数，构建满足约束条件的火箭末级六自由度动力软着陆轨迹优化模型。
[0042]
对于动力着陆轨迹优化问题，以最短着陆时间(minimum landing time，mlt)为轨迹优化问题的目标函数，即：
[0043][0044]
tf为着陆时间；
[0045]
同时，在动力软着陆过程中需满足以下约束：
[0046]
a)燃料约束
[0047]mdry
≤m(t)
[0048]
其中，m
dry
为飞行器干重。
[0049]
b)推力约束
[0050]
t
min
≤||tb(t)||2≤t
max
[0051]
cosδ
max
||tb(t)||2≤e1·
tb(t)
[0052]
t
min
和t
max
分别是推力大小的下限和上限，δ
max
为最大推力偏角，e1为单位矩阵。
[0053]
c)斜坡角约束
[0054]
同样，为了防止火箭末级着陆过程中与地面斜坡、凸起等地形发生碰撞，设置了如下斜坡约束：
[0055][0056]h23
＝[e2；e3]
[0057]
其中，ei表示第i个元素为1的单位向量，γ
gs
为最大斜坡角。
[0058]
d)倾斜角约束
[0059]
为了避免优化的轨迹存在过大的倾斜角，需要限制倾斜角大小，
[0060][0061]
其中，θ
max
为最大倾斜角，为了满足姿态翻转机动要求，防止姿态来回大幅摆动，对最大倾斜角θ
max
的约束设置为随时间动态缩小，可以将cosθ
max
设置为与时间t相关的三次多项式，即：
[0062][0063]
其中，k0,k1,k2,k3分别为多项式系数，根据飞行器控制特性选取。
[0064]
e)角速度约束
[0065]
||ωb(t)||2≤ω
max
[0066]
其中，ω
max
为最大角速度；
[0067]
f)边界约束
[0068]
火箭动力软着陆初末端状态约束为：
[0069][0070]
其中，m
wet
是火箭含燃料初始总质量，r
l,i
，v
l,i
分别是初始位置和速度向量，q
b/l,i
，ω
b,i
分别是初始姿态四元数和角速度向量，v
l,f
为终端速度，q
b/l,f
为终端姿态四元数。
[0071]
s3，将步骤s2轨迹优化模型中的非凸约束通过线性化方式转化为凸约束，得到凸形式的轨迹优化模型，所述轨迹优化模型中的非凸约束包括动力学方程和推力大小约束。
[0072]
因为上述动力软着陆轨迹优化问题包含较多非凸约束，主要存在于动力学方程和推力大小约束，该步骤主要针对这两部分的非凸约束转化为凸约束。本发明采用序列凸优化的方法将非凸的轨迹优化问题转化为一个凸形式的二阶锥子问题(socp)，通过多次迭代求解二阶锥子问题最终使优化结果收敛至最优轨迹。
[0073]
a)动力学模型线性化
[0074]
针对非凸形式的动力学方程，采用线性化的方法凸化动力学方程约束。为了方便表示，定义状态变量x(t)∈r
14
和控制变量u(t)∈r3。
[0075][0076]
u(t)＝tb(t)
[0077]
动力学方程可以简写为：
[0078][0079]
将时间变量t标准化为τ∈[0,1]，上式可以表示为：
[0080][0081]
记τ与t的缩比系数以τ为自变量，上式动力学方程可以表示为：
[0082][0083]
为了将整个非线性动力学方程转化为凸约束的形式，可以将动力学方程在参考轨迹(由和构成)进行一阶泰勒展开，得到近似的线性化动力学方程如下：
[0084][0085]
b)控制约束线性化
[0086]
由推力约束形式可知，最小推力约束是一个非凸约束，可以通过线性化将最小推力约束进行凸化处理，定义函数g:r3→
r。
[0087]
g(u(τ))＝t
min-||u(τ)||2≤0
[0088]
同样，将上式在参考轨迹(由和构成)处进行一阶泰勒展开，得到近似的线性化最小推力约束。
[0089][0090]
其中，为参考轨迹控制变量。
[0091]
s4，采用一阶保持器的方法将步骤s3中凸形式的轨迹优化模型进行离散化处理，得到离散凸化模型。
[0092]
为了将连续时间的优化控制问题转化为一个优化参数在有限维度的优化问题，需要将整个问题进行离散化处理，该步骤采用一阶保持器的离散方法，将标准化的时间离散为k个等间距的离散点
[0093][0094]
对控制变量采用一阶保持器的方式对每个步长的控制变量进行插值，在τ∈[τk,τ
k 1
]，使
[0095]
u(τ)＝αk(τ)uk βk(τ)u
k 1
,τ∈[τk,τ
k 1
],k＝0,1,
…
,k-2
[0096][0097][0098]
记φ(τ
k 1
,τk)表示无输入时状态x(τk)到x(τ
k 1
)的状态转移矩阵。状态转移矩阵φ(τ,τk)的微分形式可以表示为：
[0099][0100]
通过上述离散方法，可以将动力学方程离散为如下：
[0101][0102]
ξ为积分变量。
[0103]
s5，对步骤s3中的线性化动力学方程添加动力学松弛变量，以罚函数的形式将松弛变量添加到目标函数中。
[0104]
由于非凸动力学方程线性化过程容易造成子问题人工不可行，所以需要在线性动力学方程上引入动力学松弛变量用于弥补线性化近似误差，即
[0105][0106]
其中，上标i表示迭代序号，k表示离散点位置。为了使动力学线性化误差逐渐收敛，将动力学松弛以罚函数的形式添加到目标函数中，修改目标函数为以下形式：
[0107][0108]
其中，wv为权重系数。
[0109]
s6，设计信赖域约束限制参考轨迹的变化范围。
[0110]
由于线性化展开的局部有效性，即优化变量只有在参考点附近取值才能保证良好的线性近似效果，因此通过设计信赖域约束限制参考轨迹的变化范围。
[0111]
定义参考轨迹变化量为：
[0112][0113][0114][0115]
由于对着陆动力学方程进行了无量纲化，所以可以在子问题加入如下信赖域约束：
[0116][0117]
其中，为信赖域约束，|δxi|1,|δui|1,|δσi|1为状态变量、控制变量和时间变量的变化量。虽然使用固定信赖域约束可以得到较好的优化结果，但固定信赖域约束是一种“刚性硬约束”，当初始参考轨迹质量较低时会使子问题优化计算出现不可行，因此固定信赖域约束通过设置偏大，容易增加序列凸优化迭代计算次数。另一方面，扰动量容许范围在序列凸优化算法迭代过程中是持续变化的，与原非线性函数和对应的线性化函数的近似程度密切相关，因此可以通过线性化误差的度量来调整信赖域变化，提升收敛速度，扩大寻优范围。
[0118]
设计一种简单直接的线性化误差度量方法。定义如下两个函数和和表示在优化变量值取为时，原非线性等式和不等式约束的“满足程度”；而函数表示在参考轨迹取为扰动量取为时，线性化近似等式和不等式约束的“满足程度”。
[0119]
据此，定义针对第k次迭代结果的线性化误差度量系数ρk如下：
[0120][0121]
其中，则表示：针对原问题，第k次迭代结果相对于第k-1次迭代结果在约束满足程度上的改进；表示：针对线性化近似子问题，第k次迭代结果相对于第k-1次迭代结果在约束满足程度上的改进。当系数ρk接近于1时，则认为当前线性化近似精度较高；而当ρk远小于1，则认为当前线性化精度过低。因此，可以根据线性化误差度量系数设计如下动态信赖域策略，将信赖域约束按照以下分段函数进行动态调整：
[0122][0123]
其中，σ1，σ2为信赖域罚参数更新比例系数，满足0＜σ2＜1，σ1＞1。r0，r1，r2为信赖域更新判据参数，满足r0＜r1＜r2＜1。
[0124]
s7，通过初末端状态约束简单线性插值得到初始迭代参考轨迹，以该初始迭代参考轨迹作为序列凸化的迭代初值，考虑步骤s5，s6的约束，求解步骤s4中得到的离散凸化模型；重复迭代求解使轨迹收敛到最优轨迹，完成一个制导周期采样点的轨迹优化。
[0125]
采用上述离散及凸化处理方法，将连续最优控制问题转化为序列凸问题，其内层循环为二阶锥规划问题，外层为迭代循环过程，采用内点法求解器对一系列socp子问题进行求解。设计变量包括最优控制问题的状态变量、控制变量以及着陆经历的时间，约束条件为动力学微分、路径约束以及边界条件。
[0126]
序列凸化算法初始化轨迹可以通过初末端状态约束简单线性插值得到，采用序列凸化算法逐渐迭代收敛到最优轨迹。
[0127]
火箭末级动力软着陆问题的序列凸化算法求解流程如下：
[0128]
[0129][0130]
s8，将一个制导周期采样点的轨迹优化结果作为制导指令生成器，根据当前状态更新执行步骤s2～s7，利用轨迹优化结果更新最优指令，并直接用作制导信号，最终完成火箭子级姿态翻转着陆在线制导。
[0131]
以上文六自由度轨迹优化方法为核心，将优化方法嵌入模型预测控制(model predictive control，mpc)框架。mpc制导是将轨迹优化算法作为制导指令生成器，通过在线滚动时域地执行优化计算不断根据当前状态更新最优指令，并将最新指令直接用作制导信号。mpc制导通过实时采样状态信息执行优化计算更新控制指令，因此对环境参数的不确定性和扰动偏差具有较强的适应性。对于mpc制导方法，其稳定性、精确性和鲁棒性取决于求解相应优化问题的快速性，即优化计算更新频率是mpc算法能力的重要评价指标。
[0132]
考虑如下非线性自治系统模型：
[0133][0134]
其中；x是状态变量，u是控制变量。设在采样时刻ti的飞行器的实际状态为xr(ti)，基于xr(ti)按照前面的轨迹优化方法得到最优控制u
*
(t,xr(ti))，t∈[ti,t
fi
]。按照xr(ti)得到最优控制u
*
(t,xr(ti))作为制导控制指令应用于这个制导周期[ti,t
i 1
]，下一个制导
周期则以采样时刻t
i 1
的实际系统状态为xr(t
i 1
)作为优化算法的初始状态，即以每个制导周期初始时刻系统实际状态xr作为每次优化更新计算的初始状态。可见，mpc框架的最优制导算法可以看作每个制导周期更新初始状态求解轨迹优化问题从而更新制导指令。
[0135]
为方便后面方法论述，可以将mpc框架中每个制导周期的轨迹优化问题抽象为一般数值优化问题(mpc0)：
[0136][0137][0138][0139][0140]
其中，为轨迹优化问题的状态变量，ti为当前问题的更新采样时刻，为当前问题的终端飞行时刻，表示过程约束，u(t)∈u为控制约束，表示初始状态约束；表示终端状态约束。
[0141]
在mpc制导指令优化更新中，可能存在优化计算不可行情况，特别是在接近着陆点时，因为着陆后期通常采用满推力控制，降低了着陆末端的调节能力，若实际状态xr(ti)相对第(i-1)次计算得到x
i-1
(ti)的结果发生了较大的偏差，可能导致轨迹优化问题的多种约束条件无法满足使得当前优化计算的不可行。
[0142]
为了避免上述情况，避免着陆末端的轨迹优化算法出现不可行的情况，可以将轨迹优化问题(mpc0)的终端状态约束进行松弛并在性能指标中加以惩罚，建立一个新的轨迹优化问题rmpc：
[0143][0144][0145][0146][0147]
其中，是松弛的终端状态约束，j0为原性能指标，ω为终端状态变量与期望终端约束偏离度的罚参数。在着陆末端，通过将严苛的终端状态转化为一个软约束，确保优化问题的可行性，同时，将终端状态误差引入目标函数，从而使终端状态误差尽可能小。
[0148]
另一方面，在mpc制导指令的更新计算时，飞行时间都小于上次计算，所以在下一个制导周期可以减少当前使用的离散点数量而不会影响离散精度。因此，为了提高轨迹优化问题的计算效率，设计离散点数量的动态调整策略：
[0149][0150]
其中，floor(
·
)为向下取整函数，ni是第i次执行优化计算的离散点数量，是第i
次优化结果的飞行时间，δt为制导周期时间，是对第(i 1)次优化计算的飞行时间的估计。κ是离散点数量调节因子，其值小于1，这样可以使离散点数量的减小速率小于飞行时间的缩短速率。因此，随着飞行时间的不断减小，离散点数量减少使得问题规模逐渐变小，计算速度逐渐提升。
[0151]
综上，整个mpc制导算法如图2所示，具体流程如下：
[0152][0153][0154]
整个制导方法仿真结果如图3所示。
[0155]
图3a和图3b展示了飞行器在着陆过程中的轨迹和姿态变化，虚线表示着陆轨迹，粗实线表示箭体轴向，细实线表示推力矢量，细实线长度代表推力大小。图3e给出了推力大小的变化曲线，推力大小变化形式近似“bang-bang”控制形式。图3c、图3d、图3f-图3h分别给出了火箭末级着陆过程位置、速度分量、质量、速度大小以及角速度各状态变量的变化，整个mpc计算制导仿真飞行时间为14.707s，仿真结果表示：在设置的初始位置和姿态下，通过推力控制可以着陆到指定地点，且位置、速度和姿态等满足要求。终端水平方向的位置偏差为0.028m；落点竖直方向速度vz为-0.563m/s，水平面方向速度vx为0.171m/s，vy为0.228m/s；终端角速度大小为0.775
°
/s，倾斜角为0.423
°
。火箭末级着陆时最终质量为149293.629kg，燃料消耗为10706.371kg，满足燃料消耗约束。从制导仿真结果可以看出，本发明的制导方法可以通过在线实施更新制导指令达到很高的制导精度，且燃料消耗低，整体制导效果好。
[0156]
以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，
可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。
[0157]
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。