CN114721261A - 一种火箭子级姿态翻转着陆在线制导方法 - Google Patents

一种火箭子级姿态翻转着陆在线制导方法 Download PDF

Info

Publication number
CN114721261A
CN114721261A CN202210247585.1A CN202210247585A CN114721261A CN 114721261 A CN114721261 A CN 114721261A CN 202210247585 A CN202210247585 A CN 202210247585A CN 114721261 A CN114721261 A CN 114721261A
Authority
CN
China
Prior art keywords
landing
convex
attitude
constraint
track
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210247585.1A
Other languages
English (en)
Inventor
任宽
吴胜宝
汪小卫
张烽
焉宁
李扬
刘丙利
胡冬生
张柳
张雪梅
郝宇星
邓思超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Academy of Launch Vehicle Technology CALT
Original Assignee
China Academy of Launch Vehicle Technology CALT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Academy of Launch Vehicle Technology CALT filed Critical China Academy of Launch Vehicle Technology CALT
Priority to CN202210247585.1A priority Critical patent/CN114721261A/zh
Publication of CN114721261A publication Critical patent/CN114721261A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明提供了一种火箭子级姿态翻转着陆在线制导方法，包括：S1，建立火箭末级六自由度着陆动力学模型；S2，构建满足约束条件的火箭末级六自由度动力软着陆轨迹优化模型；S3，将S2轨迹优化模型中的非凸约束通过线性化方式转化为凸约束，得到凸形式的轨迹优化模型；S4，将S3中凸形式的轨迹优化模型进行离散化处理；S5，对S3中的线性化动力学方程添加动力学松弛变量；S6，设计信赖域约束限制参考轨迹的变化范围；S7，确定初始迭代参考轨迹；S8，求解离散凸化模型；S9，重复S8迭代求解使轨迹收敛到最优轨迹，完成一个制导周期采样点的轨迹优化；S10，利用轨迹优化结果更新最优指令，并直接用作制导信号，最终完成火箭子级姿态翻转着陆在线制导。

Description

一种火箭子级姿态翻转着陆在线制导方法
技术领域
本发明属于飞行器制导技术领域，特别涉及一种火箭子级姿态翻转着陆在线制导方法。
背景技术
随着航天技术的发展，重复使用已经成为未来航天运输系统的主要发展方向之一。目前，运载火箭一子级重复使用技术发展迅速，但对于更高势能火箭子级(如火箭末级)，回收难度更大，其飞行剖面与运载火箭一子级相比具有较大差别。传统的运载火箭一子级垂直动力着陆阶段可以不需要推力参与飞行器的姿态控制(如Falcon9)，主要通过RCS(Reaction Control System，反应控制系统)和栅格舵来实现姿态的调整，因此大多数学者针对火箭一子级回收着陆轨迹优化模型通常是不考虑姿态的三自由度动力学模型，但由于重复使用火箭末级动力着陆初始时刻以大攻角气动减速状态下降，所以在最后着陆阶段需要通过发动机推力进行姿态翻转机动并完成垂直动力着陆，因此需要设计一种新的着陆制导方法。
发明内容
为了克服现有技术中的不足，本发明人进行了锐意研究，提供了一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，在传统三自由度返回着陆制导方法的基础上引入姿态变量，建立了六自由度返回着陆问题，用于解决着陆过程中姿态翻转机动问题，同时，位置、速度、姿态满足终端约束要求，从而完成本发明。
本发明提供的技术方案如下：
目前，大多数动力着陆轨迹优化方法针对的模型都是类似于火箭一子级的三自由度(3Degree of Freedom，3DOF)质心动力学模型，在着陆过程中并未考虑姿态的变化，这是因为火箭一子级垂直动力着陆阶段姿态变化较小，通常并不需要推力参与飞行器的姿态控制，主要通过RCS和栅格舵来实现姿态的调整。但由于“星舰”式火箭末级着陆段初始姿态处于大攻角气动减速状态，所以它在最后着陆阶段需要通过发动机推力进行姿态翻转机动并完成垂直动力着陆。因此，对于这种飞行器的垂直着陆轨迹优化需要考虑推力对姿态的影响。
针对“星舰”式火箭子级动力着陆过程采用推力进行姿态翻转机动的特点，提出了一种基于改进序列凸优化的六自由度轨迹优化方法，该方法引入姿态四元数和角速度变量建立了六自由度着陆轨迹优化问题，通过凸化处理和离散化方式将非凸的轨迹优化问题转化为一系列有限维的二阶锥子问题，设计虚拟控制和动态信赖域更新策略保证二阶锥子问题的可行性和收敛性，然后通过求解一系列二阶锥子问题逐渐收敛至最优轨迹，最后以六自由度动力着陆轨迹优化方法为基础建立火箭动力软着陆滚动时域MPC(模型预测控制，Model Predictive Control)制导框架，通过仿真实验来验证建立的基于轨迹序列凸优化MPC制导方法的有效性。
具体方案如下
一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，包括如下步骤：
S1，引入姿态四元数和转动角速度描述火箭末级的姿态运动，建立火箭末级六自由度着陆动力学模型；
S2，以最短着陆时间为目标函数，构建满足约束条件的火箭末级六自由度动力软着陆轨迹优化模型；
S3，将步骤S2轨迹优化模型中的非凸约束通过线性化方式转化为凸约束，得到凸形式的轨迹优化模型，所述轨迹优化模型的非凸约束包括动力学方程和推力大小约束；
S4，采用一阶保持器的方法将步骤S3中凸形式的轨迹优化模型进行离散化处理，得到离散凸化模型；
S5，对步骤S3中线性化动力学方程添加动力学松弛变量，以罚函数的形式将松弛变量添加到目标函数中；
S6，设计信赖域约束限制参考轨迹的变化范围；
S7，通过初末端状态约束简单线性插值得到初始迭代参考轨迹；
S8，步骤S7中的初始迭代参考轨迹作为序列凸化的迭代初值，考虑步骤S5，S6的约束，求解步骤S4中得到的离散凸化模型；
S9，重复步骤S8迭代求解使轨迹收敛到最优轨迹，完成一个制导周期采样点的轨迹优化；
S10，将步骤S9的一个制导周期采样点的轨迹优化结果作为制导指令生成器，根据当前状态更新执行步骤S2-S9，将轨迹优化结果更新最优指令，并直接用作制导信号，最终完成火箭子级姿态翻转着陆在线制导。
根据本发明提供的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，具有以下有益效果：
(1)本发明针对火箭末级姿态翻转着陆的特点，在传统三自由度动力着陆模型下引入姿态变量建立六自由度动力着陆轨迹优化问题，提升了任务适应性。
(2)本发明对传统固定信赖域的序列凸化方法进行改进，设计了基于线性化误差度量的动态信赖域策略，提升了轨迹优化方法的求解效率。
(3)本发明建立了以六自由度动力着陆轨迹优化方法为基础的模型预测制导框架，成功将在线轨迹优化方法运用于制导控制系统中，实现火箭子级姿态翻转着陆全程在线制导。
附图说明
图1为本发明的动力着陆示意图；
图2为本发明的制导方法框架结构图；
图3为本发明的制导方法仿真结果图，其中图3a为6DOF三维轨迹图；图3b为6DOF纵向平面轨迹投影图；图3c为位置变化曲线图；图3d为速度分量变化曲线图；图3e为推力大小变化曲线图；图3f为质量变化曲线图；图3g为速度大小变化曲线图；图3h为角速度分量变化曲线图。
具体实施方式
下面通过对本发明进行详细说明，本发明的特点和优点将随着这些说明而变得更为清楚、明确。
在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。
本发明提供了一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，包括如下步骤：
S1，引入姿态四元数和转动角速度描述火箭末级的姿态运动，建立火箭末级六自由度着陆动力学模型。
以目标点地面坐标系FL和本体坐标系FB为参考坐标系，描述飞行器与目标点相对位置与姿态关系，如图1所示。考虑火箭末级在最后动力着陆阶段需要通过发动机推力进行一次姿态翻转机动以实现垂直动力着陆，所以在三自由度着陆动力学模型的基础上，引入姿态四元数(qB/L)以及转动角速度(ωB)以描述箭体的姿态运动，得到的六自由度着陆动力学模型的具体形式如下：
Figure BDA0003545658870000041
其中，m为飞行器的质量；αm＝1/(Ispg0)，Isp为比冲，g0为标准地球重力加速度；t为运行时刻，rL、vL为飞行器在目标点地面坐标系FL下的位置、速度矢量；DL＝-0.5ρSCD||vL||vL为在目标点地面坐标系FL下的气动阻力，ρ大气密度，S为参考面积，CD为阻力系数；TB为在本体坐标系FB下的推力；ωB为在本体坐标系FB下本体坐标系FB相对目标点地面坐标系FL的角速度矢量；qB/L为目标点地面坐标系FL相对本体坐标系FB的单位姿态四元数；gL为重力矢量；rT,B为推力作用点到质心力矩的力臂；JB为惯性矩阵；CT/B(t)为本体坐标系FB到目标点地面坐标系FL的旋转矩阵，展开形式为：
Figure BDA0003545658870000051
其中，q0，q1，q2，q3为对应的姿态四元数；
Ω(ωB)为四元数的乘法展开形式：
Figure BDA0003545658870000052
其中，ωB＝[ωBx ωBy ωBz]
[ζ×]为叉乘展开形式：
Figure BDA0003545658870000053
其中，ζ＝[ζx ζx ζx]为变量。
S2，以最短着陆时间为目标函数，构建满足约束条件的火箭末级六自由度动力软着陆轨迹优化模型。
对于动力着陆轨迹优化问题，以最短着陆时间(Minimum Landing Time，MLT)为轨迹优化问题的目标函数，即：
Figure BDA0003545658870000054
tf为着陆时间；
同时，在动力软着陆过程中需满足以下约束：
a)燃料约束
mdry≤m(t)
其中，mdry为飞行器干重。
b)推力约束
Tmin≤||TB(t)||2≤Tmax
cosδmax||TB(t)||2≤e1·TB(t)
Tmin和Tmax分别是推力大小的下限和上限，δmax为最大推力偏角，e1为单位矩阵。
c)斜坡角约束
同样，为了防止火箭末级着陆过程中与地面斜坡、凸起等地形发生碰撞，设置了如下斜坡约束：
Figure BDA0003545658870000061
H23＝[e2；e3]
其中，ei表示第i个元素为1的单位向量，γgs为最大斜坡角。
d)倾斜角约束
为了避免优化的轨迹存在过大的倾斜角，需要限制倾斜角大小，
Figure BDA0003545658870000062
其中，θmax为最大倾斜角，为了满足姿态翻转机动要求，防止姿态来回大幅摆动，对最大倾斜角θmax的约束设置为随时间动态缩小，可以将cosθmax设置为与时间t相关的三次多项式，即：
Figure BDA0003545658870000063
其中，k0,k1,k2,k3分别为多项式系数，根据飞行器控制特性选取。
e)角速度约束
||ωB(t)||2≤ωmax
其中，ωmax为最大角速度；
f)边界约束
火箭动力软着陆初末端状态约束为：
Figure BDA0003545658870000071
其中，mwet是火箭含燃料初始总质量，rL,i，vL,i分别是初始位置和速度向量，qB/L,i，ωB,i分别是初始姿态四元数和角速度向量，vL,f为终端速度，qB/L,f为终端姿态四元数。
S3，将步骤S2轨迹优化模型中的非凸约束通过线性化方式转化为凸约束，得到凸形式的轨迹优化模型，所述轨迹优化模型中的非凸约束包括动力学方程和推力大小约束。
因为上述动力软着陆轨迹优化问题包含较多非凸约束，主要存在于动力学方程和推力大小约束，该步骤主要针对这两部分的非凸约束转化为凸约束。本发明采用序列凸优化的方法将非凸的轨迹优化问题转化为一个凸形式的二阶锥子问题(SOCP)，通过多次迭代求解二阶锥子问题最终使优化结果收敛至最优轨迹。
a)动力学模型线性化
针对非凸形式的动力学方程，采用线性化的方法凸化动力学方程约束。为了方便表示，定义状态变量x(t)∈R14和控制变量u(t)∈R3。
Figure BDA0003545658870000072
u(t)＝TB(t)
动力学方程可以简写为：
Figure BDA0003545658870000073
将时间变量t标准化为τ∈[0,1]，上式可以表示为：
Figure BDA0003545658870000081
记τ与t的缩比系数
Figure BDA0003545658870000082
以τ为自变量，上式动力学方程可以表示为：
Figure BDA0003545658870000083
为了将整个非线性动力学方程转化为凸约束的形式，可以将动力学方程在参考轨迹(由
Figure BDA0003545658870000084
和
Figure BDA0003545658870000085
构成)进行一阶泰勒展开，得到近似的线性化动力学方程如下：
Figure BDA0003545658870000086
b)控制约束线性化
由推力约束形式可知，最小推力约束是一个非凸约束，可以通过线性化将最小推力约束进行凸化处理，定义函数g:R3→R。
g(u(τ))＝Tmin-||u(τ)||2≤0
同样，将上式在参考轨迹(由
Figure BDA0003545658870000087
和
Figure BDA0003545658870000088
构成)处进行一阶泰勒展开，得到近似的线性化最小推力约束。
Figure BDA0003545658870000089
其中，
Figure BDA00035456588700000810
为参考轨迹控制变量。
S4，采用一阶保持器的方法将步骤S3中凸形式的轨迹优化模型进行离散化处理，得到离散凸化模型。
为了将连续时间的优化控制问题转化为一个优化参数在有限维度的优化问题，需要将整个问题进行离散化处理，该步骤采用一阶保持器的离散方法，将标准化的时间离散为K个等间距的离散点
Figure BDA0003545658870000091
对控制变量采用一阶保持器的方式对每个步长的控制变量进行插值，在τ∈[τk,τk+1]，使
u(τ)＝αk(τ)uk+βk(τ)uk+1,τ∈[τk,τk+1],k＝0,1,…,K-2
Figure BDA0003545658870000092
Figure BDA0003545658870000093
记Φ(τk+1,τk)表示无输入时状态x(τk)到x(τk+1)的状态转移矩阵。状态转移矩阵Φ(τ,τk)的微分形式可以表示为：
Figure BDA0003545658870000094
通过上述离散方法，可以将动力学方程离散为如下：
Figure BDA0003545658870000095
ξ为积分变量。
S5，对步骤S3中的线性化动力学方程添加动力学松弛变量，以罚函数的形式将松弛变量添加到目标函数中。
由于非凸动力学方程线性化过程容易造成子问题人工不可行，所以需要在线性动力学方程上引入动力学松弛变量
Figure BDA0003545658870000096
用于弥补线性化近似误差，即
Figure BDA0003545658870000097
其中，上标i表示迭代序号，k表示离散点位置。为了使动力学线性化误差逐渐收敛，将动力学松弛以罚函数的形式添加到目标函数中，修改目标函数为以下形式：
Figure BDA0003545658870000101
其中，wv为权重系数。
S6，设计信赖域约束限制参考轨迹的变化范围。
由于线性化展开的局部有效性，即优化变量只有在参考点附近取值才能保证良好的线性近似效果，因此通过设计信赖域约束限制参考轨迹的变化范围。
定义参考轨迹变化量为：
Figure BDA0003545658870000102
Figure BDA0003545658870000103
Figure BDA0003545658870000104
由于对着陆动力学方程进行了无量纲化，所以可以在子问题加入如下信赖域约束：
Figure BDA0003545658870000105
其中，
Figure BDA0003545658870000106
为信赖域约束，|δxi|1,|δui|1,|δσi|1为状态变量、控制变量和时间变量的变化量。虽然使用固定信赖域约束可以得到较好的优化结果，但固定信赖域约束是一种“刚性硬约束”，当初始参考轨迹质量较低时会使子问题优化计算出现不可行，因此固定信赖域约束通过设置偏大，容易增加序列凸优化迭代计算次数。另一方面，扰动量容许范围在序列凸优化算法迭代过程中是持续变化的，与原非线性函数和对应的线性化函数的近似程度密切相关，因此可以通过线性化误差的度量来调整信赖域变化，提升收敛速度，扩大寻优范围。
设计一种简单直接的线性化误差度量方法。定义如下两个函数
Figure BDA0003545658870000107
和
Figure BDA0003545658870000108
Figure BDA0003545658870000109
表示在优化变量值取为
Figure BDA00035456588700001010
时，原非线性等式和不等式约束的“满足程度”；而函数
Figure BDA00035456588700001011
表示在参考轨迹取为
Figure BDA00035456588700001012
扰动量取为
Figure BDA00035456588700001013
时，线性化近似等式和不等式约束的“满足程度”。
据此，定义针对第k次迭代结果的线性化误差度量系数ρk如下：
Figure BDA00035456588700001014
其中，
Figure BDA00035456588700001015
则表示：针对原问题，第k次迭代结果
Figure BDA00035456588700001016
相对于第k-1次迭代结果
Figure BDA0003545658870000111
在约束满足程度上的改进；
Figure BDA0003545658870000112
表示：针对线性化近似子问题，第k次迭代结果相对于第k-1次迭代结果在约束满足程度上的改进。当系数ρk接近于1时，则认为当前线性化近似精度较高；而当ρk远小于1，则认为当前线性化精度过低。因此，可以根据线性化误差度量系数设计如下动态信赖域策略，将信赖域约束
Figure BDA0003545658870000115
按照以下分段函数进行动态调整：
Figure BDA0003545658870000113
其中，σ1，σ2为信赖域罚参数更新比例系数，满足0＜σ2＜1，σ1＞1。r0，r1，r2为信赖域更新判据参数，满足r0＜r1＜r2＜1。
S7，通过初末端状态约束简单线性插值得到初始迭代参考轨迹，以该初始迭代参考轨迹作为序列凸化的迭代初值，考虑步骤S5，S6的约束，求解步骤S4中得到的离散凸化模型；重复迭代求解使轨迹收敛到最优轨迹，完成一个制导周期采样点的轨迹优化。
采用上述离散及凸化处理方法，将连续最优控制问题转化为序列凸问题，其内层循环为二阶锥规划问题，外层为迭代循环过程，采用内点法求解器对一系列SOCP子问题进行求解。设计变量包括最优控制问题的状态变量、控制变量以及着陆经历的时间，约束条件为动力学微分、路径约束以及边界条件。
序列凸化算法初始化轨迹可以通过初末端状态约束简单线性插值得到，采用序列凸化算法逐渐迭代收敛到最优轨迹。
火箭末级动力软着陆问题的序列凸化算法求解流程如下：
Figure BDA0003545658870000114
Figure BDA0003545658870000121
S8，将一个制导周期采样点的轨迹优化结果作为制导指令生成器，根据当前状态更新执行步骤S2～S7，利用轨迹优化结果更新最优指令，并直接用作制导信号，最终完成火箭子级姿态翻转着陆在线制导。
以上文六自由度轨迹优化方法为核心，将优化方法嵌入模型预测控制(ModelPredictive Control，MPC)框架。MPC制导是将轨迹优化算法作为制导指令生成器，通过在线滚动时域地执行优化计算不断根据当前状态更新最优指令，并将最新指令直接用作制导信号。MPC制导通过实时采样状态信息执行优化计算更新控制指令，因此对环境参数的不确定性和扰动偏差具有较强的适应性。对于MPC制导方法，其稳定性、精确性和鲁棒性取决于求解相应优化问题的快速性，即优化计算更新频率是MPC算法能力的重要评价指标。
考虑如下非线性自治系统模型：
Figure BDA0003545658870000131
其中；x是状态变量，u是控制变量。设在采样时刻ti的飞行器的实际状态为xR(ti)，基于xR(ti)按照前面的轨迹优化方法得到最优控制u*(t,xR(ti))，t∈[ti,tfi]。按照xR(ti)得到最优控制u*(t,xR(ti))作为制导控制指令应用于这个制导周期[ti,ti+1]，下一个制导周期则以采样时刻ti+1的实际系统状态为xR(ti+1)作为优化算法的初始状态，即以每个制导周期初始时刻系统实际状态xR作为每次优化更新计算的初始状态。可见，MPC框架的最优制导算法可以看作每个制导周期更新初始状态求解轨迹优化问题从而更新制导指令。
为方便后面方法论述，可以将MPC框架中每个制导周期的轨迹优化问题抽象为一般数值优化问题(MPC0)：
Figure BDA0003545658870000132
Figure BDA0003545658870000133
Figure BDA0003545658870000134
Figure BDA0003545658870000135
其中，
Figure BDA0003545658870000136
为轨迹优化问题的状态变量，ti为当前问题的更新采样时刻，
Figure BDA0003545658870000137
为当前问题的终端飞行时刻，
Figure BDA0003545658870000138
表示过程约束，u(t)∈U为控制约束，
Figure BDA0003545658870000139
表示初始状态约束；
Figure BDA00035456588700001310
表示终端状态约束。
在MPC制导指令优化更新中，可能存在优化计算不可行情况，特别是在接近着陆点时，因为着陆后期通常采用满推力控制，降低了着陆末端的调节能力，若实际状态xR(ti)相对第(i-1)次计算得到xi-1(ti)的结果发生了较大的偏差，可能导致轨迹优化问题的多种约束条件无法满足使得当前优化计算的不可行。
为了避免上述情况，避免着陆末端的轨迹优化算法出现不可行的情况，可以将轨迹优化问题(MPC0)的终端状态约束进行松弛并在性能指标中加以惩罚，建立一个新的轨迹优化问题RMPC：
Figure BDA0003545658870000141
Figure BDA0003545658870000142
Figure BDA0003545658870000143
Figure BDA0003545658870000144
其中，
Figure BDA0003545658870000145
是松弛的终端状态约束，J0为原性能指标，ω为终端状态变量与期望终端约束偏离度的罚参数。在着陆末端，通过将严苛的终端状态转化为一个软约束，确保优化问题的可行性，同时，将终端状态误差引入目标函数，从而使终端状态误差尽可能小。
另一方面，在MPC制导指令的更新计算时，飞行时间都小于上次计算，所以在下一个制导周期可以减少当前使用的离散点数量而不会影响离散精度。因此，为了提高轨迹优化问题的计算效率，设计离散点数量的动态调整策略：
Figure BDA0003545658870000146
其中，floor(·)为向下取整函数，Ni是第i次执行优化计算的离散点数量，
Figure BDA0003545658870000149
是第i次优化结果的飞行时间，Δt为制导周期时间，
Figure BDA0003545658870000147
是对第(i+1)次优化计算的飞行时间的估计。κ是离散点数量调节因子，其值小于1，这样可以使离散点数量的减小速率小于飞行时间的缩短速率。因此，随着飞行时间的不断减小，离散点数量减少使得问题规模逐渐变小，计算速度逐渐提升。
综上，整个MPC制导算法如图2所示，具体流程如下：
Figure BDA0003545658870000148
Figure BDA0003545658870000151
整个制导方法仿真结果如图3所示。
图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，满足燃料消耗约束。从制导仿真结果可以看出，本发明的制导方法可以通过在线实施更新制导指令达到很高的制导精度，且燃料消耗低，整体制导效果好。
以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (10)

1.一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，包括如下步骤：
S1，引入姿态四元数和转动角速度描述火箭末级的姿态运动，建立火箭末级六自由度着陆动力学模型；
S2，以最短着陆时间为目标函数，构建满足约束条件的火箭末级六自由度动力软着陆轨迹优化模型；
S3，将步骤S2轨迹优化模型中的非凸约束通过线性化方式转化为凸约束，得到凸形式的轨迹优化模型，所述轨迹优化模型中的非凸约束包括动力学方程和推力大小约束；
S4，采用一阶保持器的方法将步骤S3中凸形式的轨迹优化模型进行离散化处理，得到离散凸化模型；
S5，对步骤S3中的线性化动力学方程添加动力学松弛变量，以罚函数的形式将松弛变量添加到目标函数中；
S6，设计信赖域约束限制参考轨迹的变化范围；
S7，通过初末端状态约束简单线性插值得到初始迭代参考轨迹；
S8，步骤S7中的初始迭代参考轨迹作为序列凸化的迭代初值，考虑步骤S5，S6的约束，求解步骤S4中得到的离散凸化模型；
S9，重复步骤S8迭代求解使轨迹收敛到最优轨迹，完成一个制导周期采样点的轨迹优化；
S10，将步骤S9的一个制导周期采样点的轨迹优化结果作为制导指令生成器，根据当前状态更新执行步骤S2-S9，利用轨迹优化结果更新最优指令，并直接用作制导信号，最终完成火箭子级姿态翻转着陆在线制导。
2.根据权利要求1所述的基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S1中，所述六自由度着陆动力学模型表示如下：
Figure FDA0003545658860000021
其中，m为飞行器的质量；
Figure FDA0003545658860000022
Isp为比冲，g0为标准地球重力加速度；t为运行时刻，rL、vL为飞行器在目标点地面坐标系FL下的位置、速度矢量；DL＝-0.5ρSCD||vL||vL为在目标点地面坐标系FL下的气动阻力，ρ大气密度，S为参考面积，CD为阻力系数；TB为在本体坐标系FB下的推力；ωB为在本体坐标系FB下本体坐标系FB相对目标点地面坐标系FL的角速度矢量；qB/L为目标点地面坐标系FL相对本体坐标系FB的单位姿态四元数；gL为重力矢量；rT,B为推力作用点到质心力矩的力臂；JB为惯性矩阵；CT/B(t)为本体坐标系FB到目标点地面坐标系FL的旋转矩阵，展开形式为：
Figure FDA0003545658860000023
其中，q0，q1，q2，q3为对应的姿态四元数；
Ω(ωB)为四元数的乘法展开形式：
Figure FDA0003545658860000024
其中，ωB＝[ωBx ωBy ωBz]
[ζ×]为叉乘展开形式：
Figure FDA0003545658860000025
其中，ζ＝[ζx ζx ζx]为变量。
3.根据权利要求1所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S2中，所述约束条件包括：燃料约束、推力约束、斜坡角约束、倾斜角约束、角速度约束、边界约束；以最短着陆时间为轨迹优化模型的目标函数，即：
Figure FDA0003545658860000031
tf为着陆时间。
4.根据权利要求2所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S3中，将动力学方程约束转化为凸约束的过程如下：定义状态变量x(t)∈R14和控制变量u(t)∈R3；
Figure FDA0003545658860000032
u(t)＝TB(t)
动力学方程简写为：
Figure FDA0003545658860000033
将时间变量t标准化为τ∈[0,1]的动力学方程可以表示为：
Figure FDA0003545658860000034
将动力学方程在参考轨迹进行一阶泰勒展开，得到近似的线性化动力学方程；
将推力大小约束转化为凸约束的过程如下：
非凸的最小推力约束为：g(u(τ))＝Tmin-||u(τ)||2≤0，其中，Tmin是推力大小的下限；将非凸的最小推力约束在参考轨迹进行一阶泰勒展开，得到凸化的最小推力约束。
5.根据权利要求4所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S4中，所述采用一阶保持器的方法将步骤S3中凸形式的轨迹优化模型进行离散化处理，具体通过如下方式实施：采用一阶保持器的离散方法，将标准化的时间离散为K个等间距的离散点；
Figure FDA0003545658860000035
对控制变量采用一阶保持器的方式对每个步长的控制变量进行插值，在τ∈[τk,τk+1]，使
u(τ)＝αk(τ)uk+βk(τ)uk+1,τ∈[τk,τk+1],k＝0,1,…,K-2
Figure FDA0003545658860000041
Figure FDA0003545658860000042
记Φ(τk+1,τk)表示无输入时状态x(τk)到x(τk+1)的状态转移矩阵，状态转移矩阵Φ(τ,τk)的微分形式表示为：
Figure FDA0003545658860000043
通过上述离散方法，将动力学方程离散为：
Figure FDA0003545658860000044
式中，ξ为积分变量。
6.根据权利要求5所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S5中，所述对步骤S3中的线性化动力学方程添加动力学松弛变量，以罚函数的形式将松弛变量添加到目标函数中，具体通过如下方式实施：
在线性动力学方程上引入动力学松弛变量
Figure FDA0003545658860000045
用于弥补线性化近似误差，即：
Figure FDA0003545658860000046
其中，上标i表示迭代序号，k表示离散点位置，为了使动力学线性化误差逐渐收敛，将动力学松弛以罚函数的形式添加到目标函数中，修改后的目标函数为为：
Figure FDA0003545658860000051
其中，wv为权重系数。
7.根据权利要求6所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S6中，通过线性化误差的度量来调整信赖域变化，
信赖域约束如下：
Figure FDA0003545658860000052
其中，
Figure FDA0003545658860000053
为信赖域约束；|δxi|1,|δui|1,|δσi|1为状态变量、控制变量和时间变量的变化量；
根据线性化误差度量系数ρk设计如下动态信赖域策略，将信赖域约束
Figure FDA00035456588600000511
按照以下分段函数进行动态调整：
Figure FDA0003545658860000054
其中，σ1，σ2为信赖域罚参数更新比例系数，满足0＜σ2＜1，σ1＞1；r0，r1，r2为信赖域更新判据参数，满足r0＜r1＜r2＜1。
8.根据权利要求7所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，所述线性化误差度量系数ρk为：
Figure FDA0003545658860000055
其中，
Figure FDA0003545658860000056
表示在优化变量值取为
Figure FDA0003545658860000057
时，原非线性等式和不等式约束的“满足程度”；而函数
Figure FDA0003545658860000058
表示在参考轨迹取为
Figure FDA0003545658860000059
扰动量取为
Figure FDA00035456588600000510
时，线性化近似等式和不等式约束的“满足程度”。
9.根据权利要求8所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S10中，还包括将着陆末端的轨迹优化模型的终端状态约束进行松弛并在性能指标中加以惩罚，建立一个新的轨迹优化模型：
Figure FDA0003545658860000061
Figure FDA0003545658860000062
Figure FDA0003545658860000063
Figure FDA0003545658860000064
其中，
Figure FDA0003545658860000065
为轨迹优化模型的状态变量，ti为当前模型的更新采样时刻，
Figure FDA0003545658860000066
为当前模型的终端飞行时刻，
Figure FDA0003545658860000067
表示过程约束，u(t)∈U为控制约束，
Figure FDA0003545658860000068
是松弛的终端状态约束，J0为原性能指标，ω为终端状态变量与期望终端约束偏离度的罚参数。
10.根据权利要求9所述的一种基于六自由度序列凸优化的火箭子级姿态翻转着陆在线制导方法，其特征在于，步骤S10中，在制导指令的更新时，设计离散点数量的动态调整策略：
Figure FDA0003545658860000069
其中，floor(·)为向下取整函数，Ni是第i次执行优化计算的离散点数量，
Figure FDA00035456588600000610
是第i次优化结果的飞行时间，Δt为制导周期时间，
Figure FDA00035456588600000611
是对第(i+1)次优化计算的飞行时间的估计；κ是离散点数量调节因子，其值小于1。
CN202210247585.1A 2022-03-14 2022-03-14 一种火箭子级姿态翻转着陆在线制导方法 Pending CN114721261A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210247585.1A CN114721261A (zh) 2022-03-14 2022-03-14 一种火箭子级姿态翻转着陆在线制导方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210247585.1A CN114721261A (zh) 2022-03-14 2022-03-14 一种火箭子级姿态翻转着陆在线制导方法

Publications (1)

Publication Number Publication Date
CN114721261A true CN114721261A (zh) 2022-07-08

Family

ID=82238647

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210247585.1A Pending CN114721261A (zh) 2022-03-14 2022-03-14 一种火箭子级姿态翻转着陆在线制导方法

Country Status (1)

Country Link
CN (1) CN114721261A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115421388A (zh) * 2022-09-23 2022-12-02 南京理工大学 一种基于凸优化的远程导弹末级多姿态约束在线轨迹规划方法
CN115586721A (zh) * 2022-10-08 2023-01-10 厦门大学 一种基于模型预测控制的短射程精确制导方法
CN116859981A (zh) * 2023-09-05 2023-10-10 东方空间技术(山东)有限公司 一种运载火箭姿态控制方法、装置及计算设备

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115421388A (zh) * 2022-09-23 2022-12-02 南京理工大学 一种基于凸优化的远程导弹末级多姿态约束在线轨迹规划方法
CN115586721A (zh) * 2022-10-08 2023-01-10 厦门大学 一种基于模型预测控制的短射程精确制导方法
CN116859981A (zh) * 2023-09-05 2023-10-10 东方空间技术(山东)有限公司 一种运载火箭姿态控制方法、装置及计算设备
CN116859981B (zh) * 2023-09-05 2023-12-15 东方空间技术(山东)有限公司 一种运载火箭姿态控制方法、装置及计算设备

Similar Documents

Publication Publication Date Title
CN114721261A (zh) 一种火箭子级姿态翻转着陆在线制导方法
CN102073755B (zh) 近空间高超声速飞行器运动控制仿真方法
Liu et al. Immersion and invariance-based output feedback control of air-breathing hypersonic vehicles
CN109062241B (zh) 基于线性伪谱模型预测控制的自主全射向再入制导方法
CN111580535B (zh) 一种基于凸优化的再入轨迹三维剖面规划方法及系统
Wang et al. Nonlinear aeroelastic control of very flexible aircraft using model updating
CN110015446B (zh) 一种半解析的火星进入制导方法
CN110750850A (zh) 强约束复杂任务条件下的三维剖面优化设计方法、系统及介质
Saranathan et al. Relaxed autonomously switched hybrid system approach to indirect multiphase aerospace trajectory optimization
Benedikter et al. Convex optimization of launch vehicle ascent trajectory with heat-flux and splash-down constraints
CN112395689A (zh) 基于凸优化的火箭故障后在线重构方法
Jin et al. Development and validation of linear covariance analysis tool for atmospheric entry
CN115437406A (zh) 基于强化学习算法的飞行器再入跟踪制导方法
CN113602535A (zh) 一种微纳卫星在轨自主交会控制的方法及计算机设备
CN115755598A (zh) 一种智能航天器集群分布式模型预测路径规划方法
Benedikter et al. Convex optimization of launch vehicle ascent trajectory with heat-flux and splash-down constraints
CN113504723B (zh) 一种基于逆强化学习的运载火箭减载控制方法
Wang et al. A new spacecraft attitude stabilization mechanism using deep reinforcement learning method
Cheng et al. Hover-to-cruise transition control for high-speed level flight of ducted fan UAV
Karelahti et al. Automated generation of realistic near-optimal aircraft trajectories
CN116088549B (zh) 一种尾座式垂直起降无人机姿态控制方法
Zhang et al. Prediction-based guidance algorithm for high-lift reentry vehicles
Zhao et al. Optimal guidance for lunar soft landing with dynamic low-resolution image sequences
CN115848648A (zh) 一种动能撞击小行星的伪谱轨道优化方法和装置
de Paiva et al. A robust pitch attitude controller for AURORA's semi-autonomous robotic airship

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination