导弹与航天运载技术（中英文）

参数	数值
干质量/t	120
单发动机最大推力/kN	2300
推力调节范围/%	25~100
比冲/s	350
最大发动机摆角/(°)	10
最大摆角速度$/\left({\left({}^{\circ }\right)\cdot {\mathrm{s}}^{-1}}\right)$	20
最大姿态角速度$/\left({\left({}^{\circ }\right)\cdot {\mathrm{s}}^{-1}}\right)$	36

参数	数值
干质量/t	120
单发动机最大推力/kN	2300
推力调节范围/%	25~100
比冲/s	350
最大发动机摆角/(°)	10
最大摆角速度$/\left({\left({}^{\circ }\right)\cdot {\mathrm{s}}^{-1}}\right)$	20
最大姿态角速度$/\left({\left({}^{\circ }\right)\cdot {\mathrm{s}}^{-1}}\right)$	36

类星舰飞行器翻转着陆轨迹优化

PDF下载

鲁鹏 , 李雅轩 , 刘新福

导弹与航天运载技术 | 运载器及导弹总体技术 2024,47(1): 30-38

收起

导弹与航天运载技术 | 运载器及导弹总体技术 2024, 47(1): 30-38

类星舰飞行器翻转着陆轨迹优化

全屏

鲁鹏, 李雅轩, 刘新福

作者信息

北京理工大学宇航学院,北京,100081

鲁鹏(1994—),男,博士,主要研究方向为飞行器轨迹优化、制导与控制技术。

李雅轩(1996—),男,博士研究生,主要研究方向为飞行器轨迹优化、制导与控制技术。

刘新福(1985—),男,研究员,博士生导师,主要研究方向为飞行力学、飞行器轨迹优化、制导与控制技术、模型预测控制与优化算法开发。

Flipping and Landing Trajectory Optimization of Starship-like Vehicles

Peng LU, Yaxuan LI, Xinfu LIU

Affiliations

School of Aerospace Engineering,Beijing Institute of Technology,Beijing,100081

出版时间: 2024-02-25 doi: 10.7654/j.issn.2097-1974.20240107

文章导航

摘要

收起

星舰是美国太空探索技术公司(SpaceX)研制的新一代可重复使用航天运输系统。该类型飞行器的着陆过程与传统航天器截然不同,具有姿态机动大、约束条件多、精度要求高的特点。首先研究了类星舰飞行器翻转着陆段的动力学模型,将发动机摆角响应引入动力学方程来抑制发动机摆角振荡问题。接着,基于凸优化方法设计了可以考虑姿态机动和典型约束条件下类星舰飞行器飞行时间自由的翻转着陆轨迹优化算法。为了提高算法的收敛性,使用虚拟控制、在目标函数中惩罚信赖域并给出了初值选取方法。最后,结合数值优化结果,分析了不同初值选取方法、发动机摆角响应、着陆发动机个数和优化目标选取对翻转着陆过程的影响。

关键词

类星舰飞行器 / 翻转 / 垂直着陆 / 轨迹优化 / 凸优化

Abstract

收起

Starship is a new generation of reusable space transportation system developed by Space Exploration Technologies Corporation (SpaceX). The landing process of this type of spacecraft, which has the characteristics of large attitude maneuver, multiple constraints, and high accuracy requirements, is different from that of traditional spacecraft. The dynamics of the flipping and landing phase of the Starship-like aircraft is studied firstly, and the response of engine gimbal angle is introduced into the dynamics to suppress the oscillation of engine gimbal angle. Then, based on the convex optimization method, an optimization algorithm for the free-final-time flipping and landing trajectory of Starship-like vehicles is designed, which can consider the attitude maneuver and typical constraints. In order to improve the convergence of the algorithm, virtual control and penalized trust region are used, and a method for initial value selection is given. Finally, numerical examples are provided to analyze the influences of initial value selection method, the response of engine gimbal angle, engine number, and the selection of optimization objective functions on the flipping and landing phase.

Key words

starship-like vehicles / flipping / vertical landing / trajectory optimization / convex optimization

引用本文

鲁鹏, 李雅轩, 刘新福. 类星舰飞行器翻转着陆轨迹优化. 导弹与航天运载技术, 2024 , 47 (1) : 30 -38 . DOI: 10.7654/j.issn.2097-1974.20240107

Peng LU, Yaxuan LI, Xinfu LIU. Flipping and Landing Trajectory Optimization of Starship-like Vehicles[J]. Missiles and Space Vehicles, 2024 , 47 (1) : 30 -38 . DOI: 10.7654/j.issn.2097-1974.20240107

正文

收起

0 引言

收起

航天运输系统是开展空间活动的核心支撑, 是发展航天事业、建设航天强国的基本前提,是加速发展空间科学技术、和平开发空间资源和维护国家空间安全的重要保障^[1]。近年来,随着人类空间资源开发需求的快速增长, 可重复使用运载器受到了全球的广泛关注, 世界各主要航天大国积极开展了相关的理论研究和工程实践^[2]。

瞄准星际移民和深空探测的目标, 美国太空探索技术公司(SpaceX)研制了名为星舰(Starship)的新一代可重复使用航天运输系统。2016年9月, SpaceX公司正式对外公布了 “星际运输系统” 方案^[3]。此后,星舰设计方案经历了多次重大改进和演化。自2020年起,通过开展密集的飞行测试, SpaceX公司逐步掌握了${150}\mathrm{\;m}$低空悬停、${10}\mathrm{\;{km}}$高空飞行、翻转机动和垂直定点软着陆等关键技术,项目取得重大进展。值得关注的是, 星舰采取了与传统航天器截然不同的着陆方式: 着陆前,星舰采用${90}^{\circ }$左右大攻角飞行, 在视觉上以腹部平拍的姿态向地面飞行；接近地面时, 星舰发动机开机进行推力矢量控制, 联合首尾四个气动舵, 完成姿态由水平向垂直的快速翻转, 并控制星舰减速和修正位置误差, 实现定点软着陆^[4]。图1展示了该特殊的着陆过程。这种独特的着陆方式对飞行器的轨迹优化设计提出了新的挑战。大范围的姿态机动和复杂的约束条件必须在优化设计翻转着陆飞行轨迹时得到充分的考虑。

近年来, 凸优化方法被广泛地应用于飞行器轨迹优化设计, 特别是在解决星际探测器、运载火箭等常规航天器的着陆轨迹优化问题时发挥积极作用^[5]。美国喷气推进实验室的Açıkmeşe等^[6]通过使用无损松弛和变量替换等方法将带有非凸推力约束的火星着陆轨迹优化问题转化为凸优化问题, 显著提高了轨迹优化算法的可靠性和计算效率。Szmuk等^[7]对六自由度火星动力下降制导问题进行了深入研究, 通过设计序列凸优化算法计算考虑姿态动力学的着陆飞行轨迹。相比于火星着陆, 因为气动力不可忽略, 运载火箭着陆在地球表面的轨迹优化问题求解难度更大。刘新福^[8]将推力和攻角同时作为控制量,把对应的燃料最优着陆问题顺利转化为凸优化问题进行迭代求解, 并理论证明了所采用松弛技术的有效性。杨润秋等^[9]将该研究工作拓展到飞行时间自由的问题, 用高度代替时间作为自变量对优化问题进行重构, 并基于序列凸优化设计了轨迹优化算法, 其中使用的非线性保留技术被证实可以有效提升算法收敛性。

尽管上述方法能够有效解决常规航天器的着陆轨迹优化问题, 但对类星舰飞行器的翻转着陆轨迹优化并不太适用。原因在于类星舰飞行器在滑翔段和着陆段之间存在独特的翻转过程。大范围的姿态机动导致类星舰飞行器不能被简单地处理为质点模型, 而必须在设计着陆轨迹时考虑姿态动力学和因姿态控制产生的复杂约束条件。但是, 目前针对此类轨迹优化问题的研究非常有限。

本文针对类星舰飞行器的翻转着陆轨迹优化问题进行了深入研究。首先根据类星舰飞行器在翻转着陆阶段的飞行力学特征和任务特点, 建立考虑姿态机动和典型约束条件的翻转着陆轨迹优化问题。然后研究将该优化问题通过凸化方法转化为凸优化问题, 进而设计迭代算法序列求解。通过提供合理的初值、引入虚拟控制和在目标函数中惩罚信赖域半径等方法, 该迭代算法表现出良好的收敛性。针对发动机摆角振荡问题, 本文提出将发动机摆角响应引入动力学, 这种处理被证实能够有效抑制发动机摆角振荡, 并有利于迭代算法的收敛。此外, 依托轨迹优化算法, 本文亦分析了发动机个数和优化目标选取等因素对类星舰飞行器翻转着陆过程的影响, 相关数值结果可以为类星舰飞行器的翻转着陆轨迹优化设计提供参考。

1 翻转着陆轨迹优化问题

收起

1.1 动力学建模

本文将研究类星舰飞行器在二维平面内考虑一个姿态角的翻转着陆轨迹优化问题。再入过程中飞行器控制能力强, 通常可使飞行器再入终端速度矢量和着陆点位于同一垂直平面内, 而着陆过程中的侧向运动通常会消耗额外燃料, 因此在进行类星舰飞行器轨迹优化时可以忽略侧向运动。实际飞行时, 着陆过程中的侧向位置误差可以通过控制进行修正。构造类星舰飞行器翻转着陆轨迹优化问题时会使用到着陆坐标系${Oxy}$和体坐标系${O}_{\mathrm{b}}{x}_{\mathrm{b}}{y}_{\mathrm{b}}$,如图2所示。

图3a展示了飞行器翻转着陆过程受力分析结果, 根据牛顿第二定律建立动力学方程。为了避免在问题求解时出现数值问题, 需要对动力学方程中的物理量进行无量纲化处理。长度的无量纲系数为飞行器初始高度${r}_{y0}$,速度无量纲系数为$\sqrt{{r}_{y0}{g}_{0}}\left({g}_{0}\right.$是地球表面重力加速度),时间无量纲系数为$\sqrt{{r}_{y0}/{g}_{0}}$,质量无量纲系数为初始质量${m}_{0}$。着陆坐标系下,无量纲的飞行器翻转着陆的动力学方程,记为$\dot{\mathbf{x}}\left( t\right)= \mathbf{f}\left({\mathbf{x}\left( t\right),\mathbf{u}\left( t\right)}\right)$,如下:

(1)$\dot{r}= v $

(2)$\dot{\mathbf{v}}= \left({\mathbf{T}+ \mathbf{D}}\right)/m +\mathbf{g}$

(3)$\dot{\theta }= \omega $

(4)$\dot{\omega }= \left({{\mathbf{M}}_{\mathrm{T}}+ {\mathbf{M}}_{\mathrm{D}}}\right)/{J}_{\mathrm{z}}$

(5)$\dot{m}= {\alpha }_{\mathrm{e}}T $

(6)${\dot{\delta }}_{\mathrm{d}}= \left({\delta -{\delta }_{\mathrm{d}}}\right)/{T}_{\mathrm{d}}$

式中$r \in {\mathbb{R}}^{3}$为飞行器位置矢量;$v \in {\mathbb{R}}^{3}$为飞行器速度矢量；$m$为飞行器质量；$\mathbf{T}$为发动机推力矢量；$T$为推力矢量的模,$T =\parallel \mathbf{T}\parallel ;\mathbf{D}$为气动阻力矢量;$\mathbf{g}$为重力加速度矢量；$\theta$为飞行器俯仰角 (定义见图2);$\omega$为俯仰角速度;${M}_{\mathrm{T}}$为推力产生的力矩;${M}_{\mathrm{D}}$为阻力力矩;${J}_{\mathrm{z}}$为飞行器绕质心的转动惯量;${\alpha }_{\mathrm{e}}$为质量消耗率,${\alpha }_{e}= - 1/{I}_{\mathrm{{sp}}}$,${I}_{\mathrm{{sp}}}$为发动机比冲;$\delta$为发动机摆角;${\delta }_{\mathrm{d}}$为$\delta$的响应;${T}_{\mathrm{d}}$为一阶惯性环节的无量纲时间常数。本文将发动机摆角视为控制变量, 它同时控制位置和姿态, 进行轨迹优化时会导致严重的发动机摆角振荡。鉴于此, 本文将发动机摆角响应引入动力学 (见式 (6))。数值仿真结果表明, 这一思路将有效抑制发动机摆角振荡现象,进而提高了2.3节设计的轨迹优化算法的收敛性。$x$为状态变量,$x =$${\left\lbrack \begin{array}{llllll}{\mathbf{r}}^{\mathrm{T}}& {\mathbf{v}}^{\mathrm{T}}& \theta &\omega & m &{\delta }_{\mathrm{d}}\end{array}\right\rbrack }^{\mathrm{T}},\mathbf{u}$为控制变量,$\mathbf{u}= {\left\lbrack \begin{array}{ll} T &\delta \end{array}\right\rbrack }^{\mathrm{T}}$。本文将动力学投影在着陆坐标系中, 因此无量纲的推力$\mathbf{T}$、阻力$\mathbf{D}$以及相应推力力矩${\mathbf{M}}_{\mathrm{T}}$和阻力力矩${\mathbf{M}}_{\mathrm{D}}$的计算公式如下:

(7)$\mathbf{T}= \left\lbrack \begin{matrix}\cos \theta &- \sin \theta \\\sin \theta &\cos \theta \end{matrix}\right\rbrack \cdot T \cdot \left\lbrack \begin{matrix}\cos {\delta }_{\mathrm{d}}\\\sin {\delta }_{\mathrm{d}}\end{matrix}\right\rbrack =$

(8)$\boldsymbol{D}=-C_{\text {drag }}\|\boldsymbol{v}\| \boldsymbol{v}$

(9)$\boldsymbol{M}_{\mathrm{T}}=-\left(l-l_{\mathrm{cg}}\right) T \sin \delta_{\mathrm{d}}$

(10)$\boldsymbol{M}_{\mathrm{D}}=-\left(l_{\mathrm{cg}}-l_{\mathrm{cp}}\right) C_{\text {drag }}\|\boldsymbol{v}\| \boldsymbol{v}^{\mathrm{T}}\left[\begin{array}{c} -\sin \theta \\ \cos \theta \end{array}\right]$

式中${C}_{\text{drag }}= {0.5}{C}_{\mathrm{D}}\rho {S}_{\text{ref }}/\left({{m}_{0}/{r}_{y0}}\right);\rho$为大气密度;${S}_{\text{ref }}$为飞行器参考面积;${C}_{\mathrm{D}}$为阻力系数; 由于$\rho$和${S}_{\text{ref }}$为有量纲的量,${m}_{0}/{r}_{y0}$这一项用来对$\rho {S}_{\text{ref }}$进行无量纲化处理。

1.2 其它约束条件

飞行器翻转着陆轨迹优化不仅受到动力学方程的约束,还有端点约束和过程约束,接下来详细介绍这些约束。

飞行器初始状态约束可以表述成$\Phi = 0$的形式。本优化问题给定初始位置${\mathbf{r}}_{0}^{\left( 1\right)}$、速度${\mathbf{v}}_{0}^{\left( 1\right)}$、俯仰角${\theta }_{0}^{\left( 1\right)}$、俯仰角速度${\omega }_{0}^{\left( 1\right)}$和初始质量${m}_{0}^{\left( 1\right)}$,本文使用上标 “(1)”标记翻转段对应的量,使用上标“(2)” 标记着陆段对应的量。因此$\Phi$具体表达式如下:

(11)$\Phi \left({{\mathbf{x}}^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right),{\mathbf{u}}^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right)}\right)= \left\lbrack \begin{matrix}{\mathbf{r}}^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right)- {\mathbf{r}}_{0}^{\left( 1\right)} \\{\mathbf{v}}^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right)- {\mathbf{v}}_{0}^{\left( 1\right)} \\{\theta }^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right)- {\theta }_{0}^{\left( 1\right)} \\{\omega }^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right)- {\omega }_{0}^{\left( 1\right)} \\{m}^{\left( 1\right)}\left({t}_{0}^{\left( 1\right)}\right)- {m}_{0}^{\left( 1\right)} \end{matrix}\right\rbrack $

翻转着陆过程中飞行器不能有向上的速度, 即竖直速度${v}_{y}$不能为正,因此得到速度约束:

(12)${v}_{y}^{\left( i\right)}\left( t\right)\leq 0, i = 1,2 $

为了保证翻转过程中飞行器姿态在合理范围, 保证着陆过程飞行器垂直着陆, 翻转着陆过程施加如下关于俯仰角的边界约束:

(13)${\theta }_{\min }^{\left( i\right)} \leq {\theta }^{\left( i\right)}\left( t\right)\leq {\theta }_{\max }^{\left( i\right)}, i = 1,2 $

考虑到飞行器在载人时的舒适性和姿态控制系统的能力,姿态角速度不能超过最大值${\omega }_{\max }$,即有:

(14)$- {\omega }_{\max }\leq {\omega }^{\left( i\right)}\left( t\right)\leq {\omega }_{\max }, i = 1,2 $

发动机的最大摆角记为${\delta }_{\max }$,故摆角需要满足:

(15)$- {\delta }_{\max }\leq {\delta }^{\left( i\right)}\left( t\right)\leq {\delta }_{\max }, i = 1,2 $

发动机摆角速度不能超过最大值${\dot{\delta }}_{\max }$,可以通过约束发动机摆角响应的变化率${\dot{\delta }}_{\mathrm{d}}$实现对发动机摆角速度的限制, 由式 (6) 可知, 需引入如下限制发动机摆角速度的约束:

(16)$- {\dot{\delta }}_{\max }\leq \left({{\delta }^{\left( i\right)}\left( t\right)- {\delta }_{\mathrm{d}}^{\left( i\right)}\left( t\right)}\right)/{T}_{\mathrm{d}}\leq {\dot{\delta }}_{\max }, i = 1,2 $

发动机具有推力调节能力, 单台发动机最大推力记为${T}_{\max }^{\left( i\right)}$,最小推力记为${T}_{\min }^{\left( i\right)}$,因此需要添加如下关于推力的边界约束:

(17)$ n{T}_{\min }^{\left( i\right)} \leq {T}^{\left( i\right)}\left( t\right)\leq n{T}_{\max }^{\left( i\right)}, i = 1,2 $

式中$n$为发动机个数,$n \in {\mathbb{N}}_{+ }$。

为了方便垂直着陆, 每个时刻飞行器和着陆点连线与竖直方向的夹角不能大于给定角度${\gamma }_{\mathrm{{gs}}}$,这就需要添加下滑道约束:

(18)$\begin{Vmatrix}{{\mathbf{r}}^{\left( i\right)}\left( t\right)}\end{Vmatrix}\cos {\gamma }_{\mathrm{{gs}}}\leq {r}_{y}^{\left( i\right)}\left( t\right), i = 1,2 $

式中${r}_{y}$为位置矢量在着陆坐标系$y$轴方向的分量。该下滑道约束是一个二阶锥约束。

过程约束 (12)$\sim$(18) 可以改写为如下不等式约束:

(19)${\mathbf{s}}^{\left( i\right)}\left({{\mathbf{x}}^{\left( i\right)}\left( t\right),{\mathbf{u}}^{\left( i\right)}\left( t\right)}\right)\mathrel{\text{:=}}\left\lbrack \begin{matrix}{v}_{y}^{\left( i\right)}\left( t\right)\\\vdots \\\begin{Vmatrix}{{\mathbf{r}}^{\left( i\right)}\left( t\right)}\end{Vmatrix}\cos {\gamma }_{\mathrm{{gs}}}- {r}_{y}^{\left( i\right)}\left( t\right)\end{matrix}\right\rbrack \leq 0, i = 1,2 $

末端等式约束可以表示为$E = 0$的形式,本优化问题中飞行器将着陆到给定位置${\mathbf{r}}_{f}^{\left( 2\right)}$,与此同时末端速度到达给定值${v}_{f}^{\left( 2\right)}$、俯仰角到达${\theta }_{f}^{\left( 2\right)}$并且俯仰角速度到达${\omega }_{f}^{\left( 2\right)}$。因此,$\mathbf{E}$具有如下形式:

(20)$\mathbf{E}\left({{\mathbf{x}}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right),{\mathbf{u}}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)}\right)= \left\lbrack \begin{matrix}{\mathbf{r}}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)- {\mathbf{r}}_{f}^{\left( 2\right)} \\{\mathbf{v}}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)- {\mathbf{v}}_{f}^{\left( 2\right)} \\{\theta }^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)- {\theta }_{f}^{\left( 2\right)} \\{\omega }^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)- {\omega }_{f}^{\left( 2\right)} \end{matrix}\right\rbrack $

由于燃料有限,终端质量需要满足${m}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)\geq {m}_{\mathrm{{dry}}}$,${m}_{\mathrm{{dry}}}$为燃料耗尽时飞行器的质量。

飞行器在翻转段和着陆段交接点处的飞行状态应当连续,连接约束可以表示为$\Psi = 0$,其中:

(21)$\Psi =\left\lbrack \begin{matrix}{t}_{f}^{\left( 1\right)} -{t}_{0}^{\left( 2\right)} \\{\mathbf{r}}^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {\mathbf{r}}^{\left( 2\right)}\left({t}_{0}^{\left( 2\right)}\right)\\{\mathbf{v}}^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {\mathbf{v}}^{\left( 2\right)}\left({t}_{0}^{\left( 2\right)}\right)\\{\theta }^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {\theta }^{\left( 2\right)}\left({t}_{0}^{\left( 2\right)}\right)\\{\omega }^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {\omega }^{\left( 2\right)}\left({t}_{0}^{\left( 2\right)}\right)\\{m}^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {m}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)\end{matrix}\right\rbrack $

翻转段耗时、着陆段耗时和翻转着陆段总飞行时间满足约束:

(22)${t}_{f}^{\left( i\right)} -{t}_{0}^{\left( i\right)} \geq 0, i = 1,2,{t}_{\min }\leq {t}_{f}^{\left( 2\right)} -{t}_{0}^{\left( 1\right)} \leq {t}_{\max },{t}_{0}^{\left( 1\right)} = 0 $

式中${t}_{\min }$为总飞行时间下界;${t}_{\max }$为总飞行时间上界。

1.3 翻转着陆轨迹优化问题

在优化类星舰飞行器翻转着陆的轨迹时, 为节约燃料, 同时保证翻转结束时飞行器高度不能太低, 设计了如下目标函数:

(23)$ J =- \mu {r}_{y}^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {m}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)$

式中$\mu$为惩罚系数,调节$\mu$可以改变翻转段结束时飞行器的离地高度。当$\mu = 0$时,该优化目标将转化为翻转着陆过程燃料最优。本文3.4节将详细探讨该惩罚系数的取值对翻转着陆轨迹的影响。

综上所述,可以构造翻转着陆轨迹优化问题${P0}$:

${P0}: \;\mathop{\min }\limits_{{{\mathbf{u}}^{\left( i\right)},{t}_{f}^{\left( i\right)}}}- \mu {r}_{y}^{\left( 1\right)}\left({t}_{f}^{\left( 1\right)}\right)- {m}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)\\\text{ s.t.}\left\{\begin{array}{l}{\dot{\mathbf{x}}}^{\left( i\right)}\left( t\right)= {\mathbf{f}}^{\left( i\right)}\left({{\mathbf{x}}^{\left( i\right)}\left( t\right),{\mathbf{u}}^{\left( i\right)}\left( t\right)}\right), i = 1,2 \\\Phi = 0,{\mathbf{s}}^{\left( i\right)}\left( t\right)\leq 0, E = 0,\Psi = 0, i = 1,2 \\{m}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)\geq {m}_{\text{dry }},{t}_{f}^{\left( i\right)} -{t}_{0}^{\left( i\right)} \geq 0, i = 1,2 \\{t}_{\text{min }}\leq {t}_{f}^{\left( 2\right)} -{t}_{0}^{\left( 1\right)} \leq {t}_{\max },{t}_{0}^{\left( 1\right)} = 0 \end{array}\right.$

问题${P0}$的本质是多阶段飞行时间自由的最优控制问题。因为动力学约束非线性很强,${P0}$求解难度较大。

2 凸优化求解

收起

本节将介绍使用改进的序列凸优化(Sequential Convex Programming, SCP)算法求解飞行时间自由的类星舰飞行器翻转着陆轨迹优化问题的方法。虽然序列凸优化算法无法保证收敛, 但是可以使用一些方法提高其收敛性。为了提高收敛性,求解问题${P0}$时引入虚拟控制,惩罚信赖域半径^[10],并给出了初值选取方法。为了公式的简洁性, 本节不再添加区分翻转段和着陆段的上标 “(i)”。

为了将飞行时间自由的最优控制问题${P0}$转化为等价的飞行时间固定的最优控制问题,引入变量$\tau =(t -$$\left.{t}_{0}\right)/{\Delta t}$,其中${\Delta t}= {t}_{f}- {t}_{0}$表示该阶段耗时,${t}_{0}$是该阶段的初始时间,${t}_{f}$是该阶段的终端时间。将动力学方程转换为以$\tau \in \left\lbrack {0,1}\right\rbrack$为自变量的形式:

(24)${\mathbf{x}}^{\prime }\left( t\right)= \frac{\mathrm{d}\mathbf{x}\left( t\right)}{\mathrm{d}\tau }= \frac{\mathrm{d}t}{\mathrm{\;d}\tau }\frac{\mathrm{d}\mathbf{x}\left( t\right)}{\mathrm{d}t}= {\Delta t}\cdot \mathbf{f}\left({\mathbf{x}\left( t\right),\mathbf{u}\left( t\right)}\right)$

式中上标 “$"$表示变量对$\tau$求导。状态变量和控制变量是$t$的函数,同样也是$\tau$的函数,故有:

(25)${\mathbf{x}}^{\prime }\left(\tau \right)= \mathbf{F}\left({\mathbf{x}\left(\tau \right),\mathbf{u}\left(\tau \right),{\Delta t}}\right)= {\Delta t}\cdot \mathbf{f}\left({\mathbf{x}\left(\tau \right),\mathbf{u}\left(\tau \right)}\right)$

当以$\tau$作为自变量时,初始状态约束是$\tau = 0$时的状态约束,终端状态约束是$\tau = 1$时的状态约束。

2.1 线性化

线性化可以将上述非线性动力学方程(25)转化为线性动力学方程。将方程(25)在上一次迭代得到的解附近泰勒展开并保留常数项和一阶项, 用展开后的线性方程近似原本的非线性动力学方程。这里将上一次迭代得到的解记为${\mathbf{z}}^{\left\lbrack k\right\rbrack }\left(\tau \right)= {\left\lbrack {\left({\mathbf{x}}^{\left\lbrack k\right\rbrack }\right)}^{\mathrm{T}}\;{\left({\mathbf{u}}^{\left\lbrack k\right\rbrack }\right)}^{\mathrm{T}}\right\rbrack }^{\mathrm{T}}$,时间记为$\Delta {t}^{\left\lbrack k\right\rbrack }$,那么在${z}^{\left\lbrack k\right\rbrack }\left(\tau \right)$和$\Delta {t}^{\left\lbrack k\right\rbrack }$附近线性化方程 (25), 可以得到如下线性动力学方程:

(26)${\mathbf{x}}^{\prime }\left(\tau \right)= \mathbf{A}\left(\tau \right)\mathbf{x}\left(\tau \right)+ \mathbf{B}\left(\tau \right)\mathbf{u}\left(\tau \right)+ \mathbf{C}\left(\tau \right){\Delta t}+ \mathbf{d}\left(\tau \right)$

式中

$\mathbf{A}\left(\tau \right)= {\left.\frac{\partial \mathbf{F}}{\partial \mathbf{x}}\right|}_{{\mathbf{z}}^{\left\lbrack k\right\rbrack }\left(\tau \right),\Delta {t}^{\left\lbrack k\right\rbrack }},\mathbf{B}\left(\tau \right)= {\left.\frac{\partial \mathbf{F}}{\partial \mathbf{u}}\right|}_{{\mathbf{z}}^{\left\lbrack k\right\rbrack }\left(\tau \right),\Delta {t}^{\left\lbrack k\right\rbrack }}\\\mathbf{C}\left(\tau \right)= {\left.\frac{\partial \mathbf{F}}{\partial {\Delta t}}\right|}_{{\mathbf{z}}^{\left\lbrack k\right\rbrack }\left(\tau \right)},\mathbf{d}\left(\tau \right)= -\mathbf{A}\left(\tau \right){\mathbf{x}}^{\left\lbrack k\right\rbrack }\left(\tau \right)- \mathbf{B}\left(\tau \right){\mathbf{u}}^{\left\lbrack k\right\rbrack }\left(\tau \right)$

为了构造凸优化问题, 除了动力学方程外其它非凸约束也需要进行线性化处理,例如: 假设${s}_{n}\left( z\right)$是关于$z$的非凸函数,则如下约束是非凸约束:

(27)${s}_{n}\left({z\left(\tau \right)}\right)\leq 0 $

使用一阶泰勒展开,略去高阶项,${s}_{n}\left({z\left(\tau \right)}\right)\leq 0$可以由如下线性不等式约束近似:

(28)${s}_{n}\left({{z}^{\left\lbrack k\right\rbrack }\left(\tau \right)}\right)+ {\frac{\partial {s}_{n}}{\partial z}}_{{z}^{\left\lbrack k1\right\rbrack }\left(\tau \right)}\left({z\left(\tau \right)- {z}^{\left\lbrack k\right\rbrack }\left(\tau \right)}\right)\leq 0 $

而问题${P0}$中除了动力学约束外没有其它非凸约束。

2.2 离散

为了将连续时间最优控制问题转化为非线性规划问题,需要对连续时间最优控制问题进行离散。将$\tau$所在的区间$\left\lbrack {0,1}\right\rbrack$离散成$N$份,$0 ={\tau }_{0}< {\tau }_{1}< \cdots <{\tau }_{N}=$1,式中${\tau }_{0},{\tau }_{1},\cdots ,{\tau }_{N}$是离散节点,与自变量对应的状态离散点记为${x}_{0},{x}_{1},\cdots ,{x}_{N}$,控制离散点记为${u}_{0},{u}_{1},\cdots ,{u}_{N}$。离散动力学的数值方法有很多^[11],包括欧拉法、梯形法、零阶保持器和一阶保持器等,本文使用梯形法。线性化后的动力学方程(26)经梯形法离散得到线性等式:

(29)$\left\lbrack {2\mathbf{I}/{h}_{j}- \mathbf{A}\left({\tau }_{j + 1}\right)}\right\rbrack {\mathbf{x}}_{j + 1}= \left\lbrack {2\mathbf{I}/{h}_{j}+ \mathbf{A}\left({\tau }_{j}\right)}\right\rbrack {\mathbf{x}}_{j}+ \mathbf{B}\left({\tau }_{j}\right){\mathbf{u}}_{j}+ \\\mathbf{B}\left({\tau }_{j + 1}\right){\mathbf{u}}_{j + 1}+ \mathbf{C}\left({\tau }_{j}\right){\Delta t}+ \mathbf{C}\left({\tau }_{j + 1}\right){\Delta t}+ \\\mathbf{d}\left({\tau }_{j}\right)+ \mathbf{d}\left({\tau }_{j + 1}\right), j = 0,1,\cdots ,\left({N - 1}\right)$

式中${h}_{j}= {\tau }_{j + 1}- {\tau }_{j}$。

动力学以外的其它约束,例如${s}_{n}\left({z\left(\tau \right)}\right)\leq 0$,离散形式是每个离散节点上都满足相应约束:

(30)${s}_{n}\left({z\left({\tau }_{j}\right)}\right)\leq 0, j = 0,1,\cdots , N $

2.3 虚拟控制和惩罚信赖域

非线性动力学和其它非线性约束经过线性化后, 可能使原本可行的问题变成不可行的问题。为了解决由线性化导致的问题不可行, 在动力学方程 (29) 上添加虚拟控制。离散的动力学方程(29)添加虚拟控制${v}_{\mathrm{c}}\left(\tau \right)\in {\mathbb{R}}^{8}$后得到等式约束 (31),记为${\mathbf{f}}_{\mathrm{d}}\left({{\mathbf{z}}_{j},{\mathbf{z}}_{j + 1},{\Delta t},{\mathbf{v}}_{\mathrm{c}j}}\right)= 0$,式中${\mathbf{z}}_{j}= {\left\lbrack \begin{array}{ll}{\mathbf{x}}_{j}^{\mathrm{T}}& {\mathbf{u}}_{j}^{\mathrm{T}}\end{array}\right\rbrack }^{\mathrm{T}}$。

(31)$\left\lbrack {2\mathbf{I}/{h}_{j}- \mathbf{A}\left({\tau }_{j + 1}\right)}\right\rbrack {\mathbf{x}}_{j + 1}= \left\lbrack {2\mathbf{I}/{h}_{j}+ \mathbf{A}\left({\tau }_{j}\right)}\right\rbrack {\mathbf{x}}_{j}+ \mathbf{B}\left({\tau }_{j}\right){\mathbf{u}}_{j}+ \mathbf{B}\left({\tau }_{j + 1}\right)+ \\{\mathbf{u}}_{j + 1}+ \mathbf{C}\left({\tau }_{j}\right){\Delta t}+ \mathbf{C}\left({\tau }_{j + 1}\right){\Delta t}+ \mathbf{d}\left({\tau }_{j}\right)+ \\\mathbf{d}\left({\tau }_{j + 1}\right)+ {\mathbf{v}}_{\mathrm{c}j}, j = 0,1,\cdots ,\left({N - 1}\right)$

式中${\mathbf{v}}_{\mathrm{c}j}= {\mathbf{v}}_{\mathrm{c}}\left({\tau }_{j}\right)$,虚拟控制${\mathbf{v}}_{\mathrm{c}j}$可以保证任意下一个节点处的状态${\mathbf{x}}_{j + 1}$都是可以达到的。需要注意,如果序列凸优化算法收敛时虚拟控制的值不为零, 则得到的结果就不满足线性化前的动力学约束。为了使程序收敛时虚拟控制为零, 需要在目标函数中增加关于虚拟控制的惩罚项${J}_{\mathrm{{vc}}}$,并且需要为虚拟控制设置一个较大的惩罚系数${\lambda }_{\mathrm{{vc}}}$,本文使用的${J}_{\mathrm{{vc}}}$如下:

(32)${J}_{\mathrm{{vc}}}= {\lambda }_{\mathrm{{vc}}}\mathop{\sum }\limits_{{j = 0}}^{{N - 1}}{\begin{Vmatrix}{\mathbf{v}}_{\mathrm{c}j}\end{Vmatrix}}_{\infty }$

为了保证线性化的合理性, 需要添加信赖域约束。本文使用可变信赖域, 信赖域作为惩罚项放在目标函数中, 即在目标函数中添加关于信赖域的惩罚项${J}_{\mathrm{{tr}}}$:

(33)${J}_{\mathrm{{tr}}}= {\lambda }_{\mathrm{{tr}}}\mathop{\sum }\limits_{{j = 0}}^{N}\left({{\begin{Vmatrix}{\mathbf{x}}_{j}- {\mathbf{x}}_{j}^{\left\lbrack k\right\rbrack }\end{Vmatrix}}_{\infty }+ {\begin{Vmatrix}{\mathbf{u}}_{j}- {\mathbf{u}}_{j}^{\left\lbrack k\right\rbrack }\end{Vmatrix}}_{\infty }}\right)+ {\lambda }_{\mathrm{{tr}}}\left({{\Delta t}- \Delta {t}^{\left\lbrack k\right\rbrack }}\right)$

式中${\lambda }_{\mathrm{{tr}}}$为惩罚系数。这种以添加惩罚项的方式处理信赖域半径的方法通常可以使问题收敛得更快。${\lambda }_{\mathrm{{tr}}}$取值一般较小,因为当该值较大时, 前后两次迭代的结果相差较小, 会导致序列凸优算法收敛速度变慢, 并且很容易收敛到初值附近的可行解。

至此, 可以得到时间自由的两阶段轨迹优化问题${P0}$被凸化后的优化问题${P1}$。本文在翻转段将$\tau$离散为${N}_{1}$等份,在着陆段将$\tau$离散为${N}_{2}$等份。

${P1}: \;\mathop{\min }\limits_{{{\mathbf{u}}^{\left( i\right)},\Delta {\mathbf{t}}^{\left( i\right)}}}- \mu {r}_{y}^{\left( 1\right)}\left({\tau }_{{N}_{1}}^{\left( 1\right)}\right)- {m}^{\left( 2\right)}\left({\tau }_{N2}^{\left( 2\right)}\right)+ {J}_{\mathrm{{vc}}}+ {J}_{\mathrm{{tr}}}\\\text{ s.t.}\left\{\begin{array}{l}{\mathbf{f}}_{\mathrm{d}}^{\left( i\right)}\left({{\mathbf{z}}_{j}^{\left( i\right)},{\mathbf{z}}_{j + 1}^{\left( i\right)},\Delta {t}^{\left( i\right)},{\mathbf{v}}_{\mathrm{c}j}^{\left( i\right)}}\right)= 0, j = 0,1,\cdots , N - 1 \\\Phi = 0,{\mathbf{s}}^{\left( i\right)}\left({\tau }_{j}\right)\leq 0, j = 0,1,\cdots , N \\ E = 0,\Psi = 0,\Delta {t}^{\left( i\right)} \geq 0, i = 1,2 \\{m}^{\left( 2\right)}\left({t}_{f}^{\left( 2\right)}\right)\geq {m}_{\text{dry }},{t}_{\text{min }}\leq \Delta {t}^{\left( 1\right)} +\Delta {t}^{\left( 2\right)} \leq {t}_{\text{max }}\end{array}\right.$

序列求解凸优化问题${P1}$,当${J}_{\mathrm{{vc}}}< {\epsilon }_{\mathrm{{vc}}}$并且${J}_{\mathrm{{tr}}}< {\epsilon }_{\mathrm{{tr}}}$时程序收敛,即可得到问题${P0}$的解,${\epsilon }_{\mathrm{{vc}}}$和${\epsilon }_{\mathrm{{tr}}}$是给定的收敛判据, 都是很小的正值。本文算法流程如图4所示。

2.4 初值选取

对于动力学和其它约束较为简单的优化问题, 在使用序列凸优化算法求解时, 通常可以直接使用线性初值。而对于本文构造的动力学非线性很强的优化问题, 由于线性初值表现不佳, 很容易造成算法1收敛慢或无法收敛的情况。为了提高算法1的收敛性, 在这里给出一种选取更好初值的方法。

在求翻转段初值时, 先忽略翻转过程中气动阻力产生的力矩, 然后给定发动机推力和摆角, 积分到俯仰角等于给定的${\theta }_{\mathrm{s}}$且垂直速度等于${v}_{\mathrm{{ys}}}$时结束,积分结果作为翻转段初值。${\theta }_{\mathrm{s}}$取${90}^{\circ }$附近的值便于后续垂直着陆。速度${v}_{\mathrm{{ys}}}$不能太大,需要保证着陆时速度能减到0。在计算着陆段初值时,求解在问题${P0}$着陆段基础上简化得到的凸优化问题, 以简化问题的解作为着陆段初值。具体来说就是给定简化问题的着陆时间, 忽略气动力,关于速度的动力学方程不考虑质量变化, 不考虑姿态动力学和发动机摆角,并满足部分过程和终端约束,即求解如下凸优化问题:

$\mathop{\min }\limits_{\mathbf{T}}- m\left({t}_{f}\right)+ {\int }_{{t}_{0}}^{{t}_{f}}\varepsilon \left( t\right)\mathrm{d}t \\\text{ s.t.}\left\{\begin{array}{l}\dot{\mathbf{r}}\left( t\right)= \mathbf{v}\left( t\right),\dot{\mathbf{v}}\left( t\right)= \mathbf{T}\left( t\right)/{m}_{0}+ \mathbf{g},\dot{m}\left( t\right)= -{\alpha }_{\mathrm{e}}T\left( t\right)\\{T}_{\min }\leq T\left( t\right)\leq {T}_{\max },\begin{Vmatrix}{\mathbf{T}\left( t\right)}\end{Vmatrix}\cos {\gamma }_{\mathrm{T}}\leq {T}_{y}\left( t\right)+ \varepsilon \left( t\right)\\\parallel \mathbf{T}\left( t\right)\parallel \leq T\left( t\right),{r}_{y}\left( t\right)\geq 0,\mathbf{x}\left({t}_{0}\right)= {\mathbf{x}}_{0},\mathbf{r}\left({t}_{f}\right)= {\mathbf{r}}_{f},\mathbf{v}\left({t}_{f}\right)= {\mathbf{v}}_{f}\end{array}\right.$

式中${\mathbf{x}}_{0}$取翻转段末端状态;${\gamma }_{\mathrm{T}}$为推力方向和$y$轴间允许的最大夹角;$\varepsilon$为松弛变量; 不使用$\varepsilon$计算初值的优化问题可能无解。上述优化问题中不包括$\theta ,\omega ,\delta ,{\delta }_{\mathrm{d}}$。俯仰角$\theta$根据推力方向计算,引入推力方向约束就是为了避免通过推力方向计算的俯仰角过大,$\omega$通过对$\theta$进行差商得到,$\delta \left( t\right)= {\delta }_{\mathrm{d}}\left( t\right)= 0$。求解上述凸优化问题可以得到一个较好的初值, 且计算效率高。仿真部分将比较本初值选取方法和线性初值选取方法对算法1的影响。

3 数值仿真

收起

本节将通过数值仿真来分析类星舰飞行器翻转着陆问题的特性, 希望能为相关任务设计提供一定的参考。本节数值仿真主要参考星舰公开的参数^[13],部分重要参数见表1。本文在实现算法1时使用了凸优化求解器ECOS^[14]。

飞行器高$l ={50}\mathrm{\;m},{l}_{\mathrm{{cg}}}= {0.6l},{l}_{\mathrm{{cp}}}= {0.55l}$。飞行器初始位置${\mathbf{r}}_{0}= {\left\lbrack \begin{array}{ll}{360}& {1200}\end{array}\right\rbrack }^{\mathrm{T}}\mathrm{m}$,初始速度${\mathbf{v}}_{0}= \lbrack -{18.82}$${\left.-{106.73}\right\rbrack }^{\mathrm{T}}\mathrm{m}/\mathrm{s},\;{\theta }_{0}= {170}^{\circ },\;{\omega }_{0}= 0\;\left({}^{\circ }\right)/\mathrm{s},\;{m}_{0}= {132}\mathrm{t}$。用于计算的初值${\theta }_{\mathrm{s}}= {80}^{\circ },{v}_{y\mathrm{\;s}}= -{20}\mathrm{\;m}/\mathrm{s}$。翻转过程使用3台发动机,着陆使用2台发动机。末端位置${\mathbf{r}}_{f}=$${\left\lbrack \begin{array}{ll} 0 & 0 \end{array}\right\rbrack }^{\mathrm{T}}\mathrm{m}$,末端速度${v}_{f}= {\left\lbrack \begin{array}{ll} 0 &- {0.1}\end{array}\right\rbrack }^{\mathrm{T}}\mathrm{m}/\mathrm{s},{\theta }_{f}= {90}^{\circ },{\omega }_{f}=$$0\left({}^{\circ }\right)/\mathrm{s}$。边界约束中${\theta }_{\max }^{\left( 1\right)} ={170}^{\circ },{\theta }_{\min }^{\left( 1\right)} ={45}^{\circ },{\theta }_{\max }^{\left( 2\right)} =$${105}^{\circ },{\theta }_{\min }^{\left( 2\right)} ={75}^{\circ }$。算法1参数:${N}_{1}= {40},{N}_{2}= {50}$,${\lambda }_{\mathrm{{vc}}}= {10}^{3},{\lambda }_{\mathrm{{tr}}}= 5 \times {10}^{-5}$,收敛判据${\epsilon }_{\mathrm{{vc}}}= 1 \times {10}^{-8},{\epsilon }_{\mathrm{{tr}}}=$$1 \times {10}^{-8}$,目标函数中$\mu ={0.3}$。

使用图4算法求解该算例时, 迭代7次收敛, 翻转段耗时${2.43}\mathrm{\;s}$,着陆段耗时${13.61}\mathrm{\;s}$,整个过程消耗燃料${9.773}\mathrm{t}$。图5展示了优化结果:图5a展示了翻转着陆过程的路径点和推力方向。图5b展示了俯仰角和俯仰角速度的变化, 最大的俯仰角速度小于给定上界36 。优化的控制量绘制在图5$\mathrm{c}$和图5$\mathrm{\;d}$中,由于目标函数考虑了让翻转过程高度差较小, 所以翻转时先以最大推力和最大摆角加速翻转, 接着使用小推力翻转, 同时摆角开始减小为接下来减小俯仰角速度做准备, 然后在着陆段开始时以最大推力和最小摆角为翻转过程减速, 即降低俯仰角速度。整个翻转着陆过程推力几乎是砰一砰(bang-bang)形式,发动机摆角不大于10°。

接下来改变部分仿真条件, 深入研究类星舰飞行器翻转着陆轨迹优化问题的性质, 并将对应结果与图5中的结果进行比较。

3.1 不同初值选取方法的影响

对于复杂的轨迹优化问题, 初值选取对于序列凸优化算法的收敛性有较大影响。这一小节将比较线性初值和本文初值选取方法对算法1的影响。

使用线性初值, 算法1迭代9次后收敛, 比使用本文初值选取方法多了2次迭代, 从图6a可以看出使用本文初值收敛更快。翻转段耗时${2.38}\mathrm{\;s}$,着陆段耗时${14.02}\mathrm{\;s}$,整个过程消耗燃料${9.947}\mathrm{t}$,比使用本文初值多消耗${174}\mathrm{\;{kg}}$。图6c展示了推力剖面的对比,使用线性初值时, 推力不是砰一砰形式。此外, 使用线性初值时, 发动机摆角剖面也不如使用本文初值时那么平缓, 如图6d所示。大量仿真实例表明算法1使用本文初值选取方法比线性初值收敛性更好, 在有些使用线性初值无法收敛的情况下,使用本文初值也能收敛。本文初值选取方法不仅提高了算法的收敛性, 而且优化的推力和发动机摆角剖面更适合工程使用。

3.2 发动机摆角响应相关约束的作用

在3.1节中说明了引入发动机摆角响应和相关约束的作用——避免求解的发动机摆角振荡。这一小节将给出不使用发动机摆角响应${\delta }_{\mathrm{d}}$而其它部分都相同时的计算结果。

如前所述,使用${\delta }_{\mathrm{d}}$的算例迭代7次收敛,而不使用${\delta }_{\mathrm{d}}$的算例迭代20次未收敛。图7给出了使用和不使用${\delta }_{\mathrm{d}}$计算结果的对比(不使用摆角响应的算例20次迭代未收敛,这里画了第20次迭代中的解),图7a和图7$\mathrm{c}$显示两个算例计算的推力和速度剖面几乎相同,但是发动机摆角差别较大,不使用${\delta }_{\mathrm{d}}$的算例得到的发动机摆角 (见图7b中红色实线) 着陆前出现了剧烈的振荡。从图7$\mathrm{\;d}$可以看出不使用${\delta }_{\mathrm{d}}$时,${J}_{\mathrm{{tr}}}$还未小于${\epsilon }_{\mathrm{{tr}}}$就不再减小,因此不使用${\delta }_{\mathrm{d}}$的算例没有收敛。发动机摆角剧烈振荡是此算例不好收敛的主要原因。

3.3 着陆使用发动机台数的影响

仿真使用的飞行器的发动机推力较大, 足以实现单台发动机着陆。这一小节给出单台发动机和两台发动机着陆时的对比结果。翻转段仍然使用三台发动机。

使用单台和两台发动机着陆的算例都是迭代7次收敛。图8展示了使用单台和两台发动机着陆时优化的控制量, 两个算例推力变化趋势相近, 两台发动机着陆时大推力持续时间更短。单发动机着陆消耗燃料${10.727}\mathrm{t}$,两台发动机着陆消耗燃料${9.773}\mathrm{t}$,少消耗${954}\mathrm{\;{kg}}$燃料。由此可见飞行器的推力设计对其着陆燃料消耗也是有一定影响的。

3.4 先翻转再着陆和边翻转边着陆的比较

上文的仿真算例中,目标函数中的惩罚系数$\mu$都取0.3,当目标函数中$\mu$取0时目标函数就变成了最省燃料, 因为不考虑对翻转过程高度差的优化, 可能出现翻转刚结束飞行器就接近地面的情况, 这里把这种优化目标对应的优化问题${P0}$称为边翻转边着陆问题;$\mu ={0.3}$时目标函数表示先翻转再着陆的同时节省燃料,这种优化目标对应的问题${P0}$就是本文主要讨论的先翻转再着陆问题。本节将通过仿真比较两种着陆方式的结果。

先翻转再着陆和边翻转边着陆两个算例都迭代7次收敛,边翻转边着陆时燃料消耗${8.844}\mathrm{t}$,比先翻转再着陆少消耗${929}\mathrm{\;{kg}}$燃料。图9展示了两个算例的仿真结果对比。图9a是飞行器俯仰角和俯仰角速度变化曲线,可以看出边翻转边着陆$\left({\mu = 0}\right)$时翻转过程变长,翻转过程变慢。图9b显示了质量变化,边翻转边着陆更省燃料,因为其目标函数就是最省燃料。图9$\mathrm{c}$和图9$\mathrm{\;d}$是两个算例推力和发动机摆角的对比。边翻转边着陆的前期使用最小推力, 因为不需要快速翻转；着陆前先最小推力再最大推力。这个边翻转边着陆的算例中翻转段耗时较短, 有些算例中边翻转边着陆问题中翻转过程可以占整个飞行时间的一半以上。

4 结论

收起

本文对类星舰飞行器翻转着陆轨迹优化问题进行了建模, 设计了用于求解该问题的序列凸优化算法。为了提高算法的收敛性, 该算法引入虚拟控制, 惩罚信赖域, 并给出了较好的初值计算方法。通过数值仿真可以发现, 在求解考虑姿态动力学的类星舰飞行器翻转着陆轨迹优化问题这样一个强非线性轨迹优化问题时, 与一般的序列凸优化算法相比本文设计的算法收敛性有很大提高。本文给出的翻转着陆过程的轨迹优化结果可以为相关飞行器的研发提供一定的参考。

基金

收起

国家自然科学基金(11972076)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

吴燕生. 中国航天运输系统的发展与未来[J]. 导弹与航天运载技术, 2007(5): 1-4.

Yansheng

. Development and future of space transportation system of China[J]. Missiles and Space Vehicles, 2007(5): 1-4.

[2]

王小军. 中国航天运输系统未来发展展望[J]. 导弹与航天运载技术, 2021(1): 1-6.

WANG

Xiaojun

. Future development of space transportation system of china[J]. Missiles and Space Vehicles, 2021(1): 1-6.

[3]

齐环环, 韩虹. SpaceX“超重-星舰”运输系统研制进展[J]. 国际太空, 2022(6): 36-39.

Huanhuan

, HAN

Hong

. Development progress of SpaceX's super-heavy starship transportation system[J]. Space International, 2022(6): 36-39.

[4]

李志文, 张磊, 李亮, 等. 星舰气动布局性能特点分析[J]. 空气动力学学报, 2022,40(5): 1-14.

Zhiwen

, ZHANG

Lei

, LI

Liang

, et al. Performance and chara-cteristics analysis on the Starship aerodynamic configuration[J]. Acta Aerodynamica Sinica, 2022,40(5): 1-14.

[5]

宋征宇, 王聪. 运载火箭返回着陆在线轨迹规划技术发展[J]. 宇航总体技术, 2019,3(6): 1-12.

SONG

Zhengyu

, WANG

Cong

. Development of online trajectory planning technology for launch vehicle return and landing[J]. Astronautical Systems and Engineering Technology, 2019,3(6): 1-12.

[6]

AÇIKMEȘE

, PLOEN

S R

. Convex programming approach to powered descent guidance for mars landing[J]. Journal of Guidance, Control, and Dynamics, 2007,30(5): 1353-1366.

[7]

SZMUK

, EREN

, AÇIKMEŞE

. Successive convexification for mars 6-dof powered descent landing guidance[C]// AIAA Guidance, Navigation, and Control Conference. Texas: 2017.

[8]

LIU

. Fuel-optimal rocket landing with aerodynamic controls[J]. Journal of Guidance, Control, and Dynamics, 2019,42(1): 65-77.

[9]

YANG

, LIU

. Fuel-optimal powered descent guidance with free final-time and path constraints[J]. Acta Astronautica, 2020,172: 70-81.

[10]

SZMUK

, REYNOLDS

, AÇIKMEŞE

. Successive convexification for real-time six-degree-of-freedom powered descent guidance with state-triggered constraints[J]. Journal of Guidance, Control, and Dynamics, 2020,43(8): 1399-1413.

[11]

BETTS

J T

. Practical methods for optimal control and estimation using nonlinear programming[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2010.

[12]

MALYUTA

, REYNOLDS

, SZMUK

, et al. Discretization performance and accuracy analysis for the rocket powered descent guidance problem[C]// AIAA Scitech 2019 Forum. San Diego: 2019.

[13]

Starship users guide[EB/OL].（2020-04-01）[2023-02-24]. https://www.spacex.com/media/starship_users_guide_v1.pdf. https://www.spacex.com/media/starship_users_guide_v1.pdf

[14]

DOMAHIDI

, CHU

, BOYD

. ECOS: An SOCP solver for embedded[C]// European Control Conference(ECC). Switzerland: 2013.

2024年第47卷第1期

PDF下载

455

178

引用本文

BibTeX

文章信息

doi: 10.7654/j.issn.2097-1974.20240107

接收时间：2023-03-17
首发时间：2025-07-04
出版时间：2024-02-25

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-03-17
修回日期：2023-07-03

基金

国家自然科学基金(11972076)

作者信息

北京理工大学宇航学院,北京,100081

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/ddyht/CN/10.7654/j.issn.2097-1974.20240107

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT