北京航空航天大学学报

工况	偏心率偏差	轨道倾角偏差/(°)	升交点赤经偏差/(°)	近地点幅角偏差/(°)
下降40%	0.000200	0.000006	0.000015	−0.004397
下降45%	0.000270	0.000005	0.000016	−0.004574
下降50%	0.000373	0.000004	0.000017	−0.004923

工况	偏心率偏差	轨道倾角偏差/(°)	升交点赤经偏差/(°)	近地点幅角偏差/(°)
下降40%	0.000200	0.000006	0.000015	−0.004397
下降45%	0.000270	0.000005	0.000016	−0.004574
下降50%	0.000373	0.000004	0.000017	−0.004923

运载火箭地月转移入轨混合制导方法

PDF下载

郭康 1^,2 , 李惠峰 1^,2^,* , 包为民 3

北京航空航天大学学报 | 2026,52(10): 3592-3601

收起

北京航空航天大学学报 | 2026 , 52 (10) : 3592 -3601

运载火箭地月转移入轨混合制导方法

全屏

郭康1^,2, 李惠峰1^,2^,*, 包为民3

作者信息

¹北京航空航天大学宇航学院，北京 100191

²北京航空航天大学天地往返高效运输技术全国重点实验室，北京 100191

³中国航天科技集团有限公司，北京 100048

通讯作者:

E-mail：lihuifeng@buaa.edu.cn

Hybrid guidance method for launch vehicle trans-lunar injection

Kang GUO1^,2, Huifeng LI1^,2^,*, Weimin BAO3

Affiliations

¹School of Astronautics，Beihang University，Beijing 100191，China

²State Key Laboratory of High-Efficiency Reusable Aerospace Transportation Technology，Beihang University，Beijing 100191，China

³China Aerospace Science and Technology Corporation，Beijing 100048，China

Corresponding Author:

E-mail：lihuifeng@buaa.edu.cn

出版时间: 2026-10-31 doi: 10.13700/j.bh.1001-5965.2025.0353

文章导航

摘要

收起

针对小推重比地月转移入轨任务中长飞行时间和强不确定性对运载火箭闭环制导方法带来的挑战，提出一种混合使用二阶锥规划和邻近最优的制导方法。建立考虑多源不确定性的地月转移入轨问题模型，并提出一种基于实时轨迹规划的混合制导策略。实现一种基于高斯伪谱离散与二阶锥规划的实时轨迹规划方法，综合应用协态映射定理和二阶锥规划对偶输出，能够高精度地估计协态量，为邻近最优方法提供必要基础。提出一种基于邻近最优的轨迹校正方法，在不确定性导致火箭偏离轨迹时，能够为二阶锥规划提供高质量初始参考。开展数值实验，在一个地月转移入轨任务中，与现有制导方法对比，所提制导方法对地月转移入轨任务具有更强的适应性和最优性，混合制导方法的平均计算时间降低30.67%，每制导周期的平均算时间为75.353 ms，为地月转移入轨提供了一条新且可行的技术路线。

关键词

运载火箭制导 / 地月转移入轨 / 实时轨迹规划 / 二阶锥规划 / 邻近最优

Abstract

收起

This paper addresses the critical guidance problem for launch vehicles performing trans-lunar injection missions with low thrust-to-weight ratios, where the extended flight duration and substantial multi-source uncertainties present formidable challenges to conventional closed-loop guidance methods. A novel hybrid guidance approach integrating second-order cone programming with neighboring optimal control is developed to overcome these challenges. The proposed methodology begins with establishing a comprehensive trans-lunar injection problem model that explicitly accounts for multi-source uncertainties. A hybrid guidance strategy based on real-time trajectory planning is formulated to handle these uncertainties effectively. The main innovation is the application of a real-time trajectory planning technique that uses second-order cone programming in conjunction with Gauss pseudospectral discretization. By synergistically applying the costate mapping theorem and exploiting the dual output of second-order cone programming, the method achieves high-precision costate estimation, which provides the essential foundation for the neighboring optimal control component. A trajectory correction method based on neighboring optimal control is developed to address situations where uncertainties cause significant trajectory deviations. This method generates high-quality initial references for the second-order cone programming solver, substantially improving convergence characteristics and computational efficiency. Extensive numerical experiments demonstrate the superiority of the proposed guidance method in trans-lunar injection missions. The computational performance shows remarkable improvement, with the average computation time reduced to 75.353 milliseconds per guidance cycle, representing a 30.67% reduction. In comparison to current methods, this hybrid guidance method’s successful demonstration offers a novel and feasible technical pathway for trans-lunar injection missions.

Key words

launch vehicle guidance / trans-lunar injection / real-time trajectory planning / second-order cone programming / neighboring optimal control

引用本文

郭康, 李惠峰, 包为民. 运载火箭地月转移入轨混合制导方法. 北京航空航天大学学报, 2026 , 52 (10) : 3592 -3601 . DOI: 10.13700/j.bh.1001-5965.2025.0353

Kang GUO, Huifeng LI, Weimin BAO. Hybrid guidance method for launch vehicle trans-lunar injection[J]. Journal of Beijing University of Aeronautics and Astronautics, 2026 , 52 (10) : 3592 -3601 . DOI: 10.13700/j.bh.1001-5965.2025.0353

正文

收起

月球是中国深空探测的第一站，也是未来开展行星探测和小天体探测的前哨^[1]。月球探测是目前世界各国航天工程的重要目标，相关计划包括中国探月工程四期、国际月球科研站、美国 Artemis 计划等。开展月球探测的必要条件是能够可靠、高效地进入地月转移轨道。与近地入轨不同，地月转移入轨需要更大的速度增量，在小推重比条件下飞行时间显著增加。例如，美国空间发射系统（space launch system，SLS）Block 1的地月转移入轨时间可达数百秒，地心转角可达45°^[2]。受地球重力场影响，火箭运动方程将具有更强的非线性，严重违背现有上升段制导方法的基础假设，对制导方法的可靠性和推进剂最优性提出挑战。

火箭上升段制导是一个经典的工程问题，现有方法主要分为3类：迭代制导方法、动力显式制导（powered explicit guidance，PEG）方法、闭环轨迹优化制导方法。迭代制导方法发展自美国国家航空航天局（national aeronautics and space administration，NASA）为土星五号运载火箭研制的迭代制导模式(iterative guidance mode，IGM) ^[3]。

IGM是一种路径自适应的制导方法，在火箭偏离标称轨迹时，基于线性程序角、常值重力加速度等简化假设规划未来飞行轨迹。迭代制导方法对任务和不确定性具有较强适应能力。中国部分长征系列运载火箭亦采用了迭代制导方法，在长征二号F遥八任务中得到了首次验证^[4]。国内外学者针对迭代制导开展了大量的改进工作^[5-8]。例如，文献[5]提出一种向量形式的迭代制导方法，将火箭推力方向作为控制量，与程序角作为控制量对比，其控制形式更加简洁。文献[6]提出一种能够控制入轨姿态的迭代制导方法，将迭代制导中的线性程序角扩展至二次形式，能够保障火箭发动机关机时位置、速度及姿态约束的同时满足。文献[7]提出一种改进的迭代制导方法，在迭代制导中不加简化地使用最优控制形式与横截条件，提高了迭代制导的推进剂最优性，并提出基于高斯-勒让德积分的推力积分数值计算方法，提高了制导算法的积分精度。目前，迭代制导的主要应用场景是大推重比近地轨道入轨任务。在小推重比长飞行时间任务中，迭代制导的线性程序角、常值重力加速度等基础假设不再成立。若使用迭代制导直接入轨地月转移轨道，制导性能可能显著下降^[9]。除近地轨道入轨任务，学者们还将迭代制导方法扩展应用于深空探测等对自主性要求较高的复杂场景^[10]。例如，文献[11]系统分析了月面上升与常规运载火箭上升的特性差异，提出一种专用于月面上升的迭代制导算法，并针对垂直上升段与迭代制导段之间存在的显著姿态角差异问题，设计了一种最优轨迹衔接方法。此外，文献[12]面向火星环境的高不确定性特点，开发了一种具有自适应能力的火星上升迭代制导方法，实现了在复杂火星环境下的高精度入轨控制。

PEG方法是航天飞机应用的上升与轨道转移制导方法，基于线性正切推力等简化假设，以半解析方式规划飞行轨迹^[13-15]。与IGM相比，PEG具有更强的任务适应性，兼顾了航天飞机的大推重比上升与小推重比轨道转移任务。然而，PEG在小推重比地月转移入轨中的应用同样受限。例如，文献[16]研究了PEG在SLS Block 1中的应用，针对小推重比入轨地月转移轨道时PEG的性能严重下降甚至迭代发散现象，提出推力方向变化限幅等改进方法，增强了PEG对该工况的适应能力。在实际深空探测任务中，嫦娥五号探测器的月面上升与着陆任务均采用了PEG方法^[17]。该方案通过引入预报策略，能够在PEG算法正式启动前预测飞行器状态并预迭代制导参数，从而实现对姿态角的提前预报，最终保证了飞行器能够及时完成姿态调整以满足轨道匹配的精确要求。

闭环轨迹优化制导方法通过数值方法优化火箭轨迹作为制导指令，根据轨迹优化方法的不同，可进一步分为间接法与直接法。间接法通过推导火箭轨迹优化问题的必要条件，形成两点边值问题并求解的方式规划轨迹，如文献[18]提出的OPGUID (optimal guidance)与Switch制导方法。直接法将火箭轨迹优化问题转化为数学规划问题并求解，是目前广泛使用的火箭任务重构方法^[19-25]。然而，现有文献中的轨迹优化往往只进行一次，随后以开环的方式使用其规划结果，缺乏对不确定性的适应能力。

本文主要创新包括以下3点：

1）提出一种混合使用二阶锥规划和邻近最优的制导策略。在多源不确定性导致火箭持续偏离优化轨迹、二阶锥规划收敛慢时，基于邻近最优的一阶必要条件修正火箭飞行轨迹，作为二阶锥规划优化的初始猜想。与现有方法对比，混合制导方法解决了多源不确定性下闭环轨迹优化收敛慢的问题，提升了计算效率。

2）提出一种基于高斯伪谱法与二阶锥规划的实时轨迹优化方法。与现有方法对比，主要增量在于使用高斯伪谱法的协态映射定理和二阶锥规划的对偶输出，实现了对协态量的高精度估计，为邻近最优方法提供了可靠支撑。

3）在典型地月转移入轨任务中，定量对比了本文制导方法与现有迭代制导方法，说明了本文制导方法在适应性、推进剂最优性等方面的优势，明确了应用场景。

1 地月转移入轨制导问题

收起

1.1 多源不确定性下的火箭运动方程

考虑运载火箭在大气层外的质心运动，在地心惯性坐标系中建立火箭运动方程。状态量包括质心位置

$ \boldsymbol{r} $

、质心速度

$ \boldsymbol{V} $

和质量

$ m $

，控制量包括推力方向

$ \boldsymbol{u} $

，则火箭的质心运动可描述为

$ {\dot{\boldsymbol{r}}}(t){=\boldsymbol{V}}(t) $

$ {\dot{\boldsymbol{V}}}(t)=\boldsymbol{g}[\boldsymbol{r}(t)]+\frac{(1+{\delta }_{T})T}{m(t)}(\boldsymbol{I}+{\boldsymbol{\delta }}_{\boldsymbol{u}})\boldsymbol{u}(t) $

$ \dot{m}(t)=-(1+{\delta }_{R})R $

式中：

$ t $

为时间；

$ \boldsymbol{g} $

为重力加速度；

$ \boldsymbol{I} $

为3维单位矩阵；

$ T $

为推力；

$ R $

为推进剂秒耗量；

$ {\delta }_{T} $

、

$ {\boldsymbol{\delta }}_{\boldsymbol{u}} $

和

$ {\delta }_{R} $

为本文考虑的3种不确定性：推力大小偏差、推力方向偏差、秒耗量偏差。

推力方向

$ \boldsymbol{u} $

受到单位幅值约束：

$ \left|\left|\boldsymbol{u}(t)\right|\right|=1 $

$ {\boldsymbol{\delta }}_{\boldsymbol{u}} $

为三维矩阵，满足以下约束关系：

$ \left|\left|\boldsymbol{I}+{\boldsymbol{\delta }}_{\boldsymbol{u}}\right|\right|=1 $

在火箭飞行过程中，多源不确定性

$ {\delta }_{T} $

、

$ {\boldsymbol{\delta }}_{\boldsymbol{u}} $

和

$ {\delta }_{R} $

将使火箭持续偏离标称轨迹。在缺少不确定性模型、难以准确辨识

$ {\delta }_{T} $

、

$ {\boldsymbol{\delta }}_{\boldsymbol{u}} $

和

$ {\delta }_{R} $

的条件下，实时轨迹规划是提高入轨精度和推进剂最优性的有效方法：通过在每个制导周期中重新规划未来的最优轨迹，能够避免跟踪标称轨迹带来的性能损失。该类方法的关键是如何在有限的计算时间（制导周期）内可靠地规划火箭飞行轨迹。

1.2 地月转移入轨条件

地月转移入轨的任务目标是抵达一个给定的地月转移轨道，由其半长轴

$ {a}_{\mathrm{ref}} $

、偏心率

$ {e}_{\mathrm{ref}} $

、轨道倾角

$ {i}_{\mathrm{ref}} $

、升交点赤经

$ {\mathit{\Omega }}_{\mathrm{ref}} $

、近地点幅角

$ {\omega }_{\mathrm{ref}} $

确定。与近地轨道入轨任务不同的是，地月转移轨道具有更显著的偏心率。为此，本文提出一种基于动量矩

$ \boldsymbol{L} $

与拉普拉斯向量

$ \boldsymbol{A} $

的降阶入轨条件，计算公式分别为

$ \boldsymbol{L}(\boldsymbol{r},\boldsymbol{V})=\boldsymbol{r}\times \boldsymbol{V} $

$ \boldsymbol{A}(\boldsymbol{r},\boldsymbol{V})=\boldsymbol{V}\times \boldsymbol{L}(\boldsymbol{r},\boldsymbol{V})-\mu \boldsymbol{r}/\left|\left|\boldsymbol{r}\right|\right| $

式中：

$ \mu $

为地球引力常数。

本文针对地月转移入轨问题的特点，提出以下形式的入轨条件：

$ \boldsymbol{\beta }[\boldsymbol{r}({t}_{\text{f}}),\boldsymbol{V}({t}_{\text{f}})]=\left[\begin{matrix}\boldsymbol{L}[\boldsymbol{r}({t}_{\text{f}})\quad \boldsymbol{V}({t}_{\text{f}})]-{\boldsymbol{L}}_{\text{ref}}\\\boldsymbol{n}_{1}^{\text{T}}\boldsymbol{A}[\boldsymbol{r}({t}_{\text{f}})\quad \boldsymbol{V}({t}_{\text{f}})]-\boldsymbol{n}_{1}^{\text{T}}{\boldsymbol{A}}_{\text{ref}}\\\boldsymbol{n}_{2}^{\text{T}}\boldsymbol{A}[\boldsymbol{r}({t}_{\text{f}})\quad \boldsymbol{V}({t}_{\text{f}})]-\boldsymbol{n}_{2}^{\text{T}}{\boldsymbol{A}}_{\text{ref}}\end{matrix}\right]={0} $

式中：

$ \boldsymbol{\beta } $

为入轨条件；

$ {\boldsymbol{n}}_{1} $

、

$ {\boldsymbol{n}}_{2} $

分别为入轨条件中第1、第2个方向；

$ {t}_{\text{f}} $

为入轨时刻；

$ {\boldsymbol{L}}_{\text{ref}} $

和

$ {\boldsymbol{A}}_{\text{ref}} $

分别为目标地月转移轨道的动量矩和拉普拉斯向量，计算公式为

$ {\boldsymbol{L}}_{\text{ref}}=\sqrt{\mu {a}_{\text{ref}}(1-e_{\text{ref}}^{2})}\left[\begin{matrix}\sin {i}_{\text{ref}}\sin {\mathit{\Omega }}_{\text{ref}}\\-\sin {i}_{\text{ref}}\cos {\mathit{\Omega }}_{\text{ref}}\\\cos {i}_{\text{ref}}\end{matrix}\right] $

$ {\boldsymbol{A}}_{\text{ref}}=\mu {e}_{\text{ref}} \left[\begin{matrix}\cos {\mathit{\Omega }}_{\text{ref}}\cos {\omega }_{\text{ref}}-\cos {i}_{\text{ref}}\sin {\mathit{\Omega }}_{\text{ref}}\sin {\omega }_{\text{ref}}\\\sin {\mathit{\Omega }}_{\text{ref}}\cos {\omega }_{\text{ref}}+\cos {i}_{\text{ref}}\cos {\mathit{\Omega }}_{\text{ref}}\sin {\omega }_{\text{ref}}\\\sin {{i}}_{\text{ref}}\sin {\omega }_{\text{ref}}\end{matrix}\right]$

与现有文献[25]中使用的入轨条件相比，式(8)只约束了

$ {\boldsymbol{n}}_{1} $

和

$ {\boldsymbol{n}}_{2} $

这2个特定方向的拉普拉斯向量

$ \boldsymbol{A} $

。这一方面降低了约束的数量，同时消除了约束间线性相关性，为邻近最优方法的应用提供了必要条件。考虑到拉普拉斯向量的幅值与偏心率相关，且指向近地点，将

$ {\boldsymbol{n}}_{1} $

和

$ {\boldsymbol{n}}_{2} $

设计为

$ {\boldsymbol{n}}_{1}={\boldsymbol{A}}_{\text{ref}}/\left|\left|{\boldsymbol{A}}_{\text{ref}}\right|\right| $

$ {\boldsymbol{n}}_{2}=\left({\boldsymbol{L}}_{\text{ref}}\times {\boldsymbol{A}}_{\text{ref}}\right)/\left|\left|{\boldsymbol{L}}_{\text{ref}}\times {\boldsymbol{A}}_{\text{ref}}\right|\right| $

对地月转移任务（偏心率

$ {e}_{\mathrm{ref}} $

不为0），

$ {\boldsymbol{n}}_{1} $

和

$ {\boldsymbol{n}}_{2} $

具有良好（非奇异）的定义。式(8)中的约束等价于指定轨道五要素，入轨点的真近点角将优化确定。

1.3 地月转移入轨的最优控制建模

地月转移入轨问题可建模为以下最优控制问题：

问题1　地月转移入轨问题。

$ \begin{cases}\operatorname{minimize} \quad m\left(t_0\right)-m\left(t_{\mathrm{f}}\right) \\\text {subject to} \quad \dot{\boldsymbol{r}}(t)=\boldsymbol{V}(t) \\\begin{cases}\dot{\boldsymbol{V}}(t)=\boldsymbol{g}[\boldsymbol{r}(t)]+\dfrac{T}{m(t)} \boldsymbol{u}(t) \\\dot{m}(t)=-R \\\|\boldsymbol{u}(t)\|=1 \\\boldsymbol{r}\left(t_0\right)=\boldsymbol{r}_0\\\boldsymbol{V}\left(t_0\right)=\boldsymbol{V}_0\\m\left(t_0\right)=m_0 \\\boldsymbol{\beta}\left[\boldsymbol{r}\left(t_{\mathrm{f}}\right), \boldsymbol{V}\left(t_{\mathrm{f}}\right)\right]=\mathbf{0}\end{cases}\end{cases} $

式中：

$ {t}_{0} $

为当前时刻；

$ {\boldsymbol{r}}_{0} $

、

$ {\boldsymbol{V}}_{0} $

和

$ {m}_{0} $

为火箭当前的飞行状态（质心位置、质心速度和质量）。

问题1的性能指标是最小化飞行过程中的推进剂消耗，其解表示一条以推进剂最优方式将火箭从当前状态转移至地月转移轨道的飞行轨迹，是闭环轨迹优化制导方法的基础。

然而，问题1中没有考虑火箭运动中的多源不确定性

$ {\delta }_{T} $

、

$ {\boldsymbol{\delta }}_{\boldsymbol{u}} $

和

$ {\delta }_{R} $

。在使用问题1解中的最优控制量（推力方向）时，火箭将在不确定性作用下持续偏离优化出的标称轨迹，难以精确入轨。

针对不确定性带来的影响，需要闭环轨迹规划方法重复求解问题1进行修正。制导性能的关键在于求解问题的计算效率：求解问题1越快，不确定性在运动方程中的累计作用越小。考虑到在制导周期间，当前时刻

$ {t}_{0} $

和当前飞行状态

$ {\boldsymbol{r}}_{0} $

、

$ {\boldsymbol{V}}_{0} $

和

$ {m}_{0} $

的变化较小，如何利用上一制导周期得到的轨迹初始化当前周期中的轨迹优化是提高计算效率的关键。需要注意的是，在多源不确定性

$ {\delta }_{T} $

、

$ {\boldsymbol{\delta }}_{\boldsymbol{u}} $

和

$ {\delta }_{R} $

的影响下，上一制导周期中的最优轨迹不一定是理想的初始参考。

2 基于实时轨迹优化的混合制导

收起

本文提出一种综合邻近最优和二阶锥规划的混合制导方法：在每个制导周期中（设当前制导周期为

$ k $

），首先，基于当前的飞行状态

$ \boldsymbol{r}_{0}^{k} $

、

$ \boldsymbol{V}_{0}^{k} $

和

$ m_{0}^{k} $

，应用邻近最优方法校正上一制导周期的最优飞行时间

$ t_{\text{f}}^{k-1} $

和飞行轨迹

$ ({\boldsymbol{r}}^{k-1},{\boldsymbol{V}}^{k-1},{m}^{k-1},{\boldsymbol{u}}^{k-1}) $

；然后，基于邻近最优方法得到的初始参考，应用二阶锥规划方法实时求解问题1，得到最优飞行轨迹。该混合制导方法的流程如图1所示。

图1所示制导方法的主要特点是：在二阶锥规划开展较为耗时的迭代前，先使用邻近最优方法根据当前飞行状态更新参考轨迹，降低多源不确定性的影响，使轨迹规划方法能够快速收敛，计算得到抵达目标地月转移轨道的高精度轨迹。

本文提出的混合制导策略同时使用了2种轨迹规划方法：邻近最优方法与基于二阶锥规划的数值优化方法。每个制导周期中的计算过程如图2所示，每个制导周期中将先使用计算效率高的邻近最优方法估计多源不确定性作用下新的最优轨迹，再使用二阶锥规划方法对轨迹进行进一步的校正。

在本文提出的混合制导策略下，邻近最优方法和二阶锥规划方法实现了互补：邻近最优方法计算效率高，但需要准确的协态初始参考才能发挥其优势；二阶锥规划方法虽然计算效率较低，但轨迹精度高，且能够可靠地生成协态参考。混合制导策略充分结合两者优势，应用二阶锥规划方法为邻近最优方法提供高质量初始参考，并应用邻近最优方法快速修正多源干扰带来的偏差，避免二阶锥规划进行过多的耗时迭代，从而实现兼具高轨迹精度和高计算效率的制导方法。

3 基于邻近最优的轨迹校正方法

收起

在多源不确定性导致火箭偏离预期的飞行轨迹时，应用邻近最优方法校正上一制导周期中的最优轨迹，作为本制导周期中的最优参考。具体地，给定上一制导周期中最优飞行时间

$ t_{\text{f}}^{k-1} $

和飞行轨迹

$ ({\boldsymbol{r}}^{k-1},{\boldsymbol{V}}^{k-1},{m}^{k-1},{\boldsymbol{u}}^{k-1}) $

，并给定当前时间

$ t_{0}^{k} $

和飞行状态

$ \boldsymbol{r}_{0}^{k} $

、

$ \boldsymbol{V}_{0}^{k} $

和

$ m_{0}^{k} $

，将初始状态偏差定义为

$ {\text{δ}} \boldsymbol{x}_{0}^{k}=\left[\begin{matrix}\boldsymbol{r}_{0}^{k}\\\boldsymbol{V}_{0}^{k}\\m_{0}^{k}\end{matrix}\right]-\left[\begin{matrix}{\boldsymbol{r}}^{k-1}(t_{0}^{k})\\{\boldsymbol{V}}^{k-1}(t_{0}^{k})\\{m}^{k-1}(t_{0}^{k})\end{matrix}\right] $

邻近最优校正的目标是计算初始状态偏差

$ {\text{δ}} \boldsymbol{x}_{0}^{k} $

下近似最优的飞行时间和飞行轨迹，为二阶锥规划提供高质量的初始参考。

将地月转移入轨问题（问题1）中的入轨时刻

$ {t}_{\text{f}} $

视作优化参数，并重定义控制量，可得到下述等价固定时间、无控制约束的最优控制问题：

问题2　参数化地月转移入轨问题。

$ \begin{cases}\operatorname{minimize}\quad m(0)-m(1) \\\text {subject to}\quad \dot{\boldsymbol{r}}(\tau)=\left(t_{\mathrm{f}}-t_0\right) \boldsymbol{V}(\tau) \\\begin{cases}\dot{\boldsymbol{V}}(\tau)=\left(t_{\mathrm{f}}-t_0\right)\left(\boldsymbol{g}[\boldsymbol{r}(\tau)]+\dfrac{T}{m(\tau)} \boldsymbol{u}^{\prime}(\tau)\right) \\\dot{m}(\tau)=-\left(t_{\mathrm{f}}-t_0\right) R \\\boldsymbol{r}(0)=\boldsymbol{r}_0\\\boldsymbol{V}(0)=\boldsymbol{V}_0\\m(0)=m_0 \\\boldsymbol{\beta}[\boldsymbol{r}(1), \boldsymbol{V}(1)]={\boldsymbol{0}}\end{cases}\end{cases} $

问题2对运动方程的自变量时间

$ t $

进行了变换，将问题1由终止时间

$ {t}_{\text{f}} $

自由的问题转化为固定时间区间的问题。具体地，引入新的时间自变量

$ \tau $

：

$ \tau =\frac{t-{t}_{0}}{{t}_{\text{f}}-{t}_{0}} $

则

$ \tau =0 $

对应

$ t={t}_{0} $

，

$ \tau =1 $

对应

$ t={t}_{\text{f}} $

。此时，入轨时刻

$ {t}_{\text{f}} $

成为待优化的参数。

此外，问题2中将控制量

$ \boldsymbol{u} $

重新定义为

$ \boldsymbol{u}'= ({u}_{a},{u}_{b}) $

，其与原本控制量

$ \boldsymbol{u} $

的关系为

$ \left[\begin{matrix}\cos {u}_{a}\cos {u}_{b}\\\cos {u}_{a}\sin {u}_{b}\\\sin {u}_{a}\end{matrix}\right]=\boldsymbol{u} $

通过上述时间变换和控制量重定义，地月转移入轨问题转化为含参数的固定时间最优控制问题。问题2的必要条件包括：

$ \begin{cases}{\dot{\boldsymbol{x}}}(\tau )=\boldsymbol{f}[\boldsymbol{x}(\tau ),\boldsymbol{u}'(\tau ),{t}_{\text{f}}]\\{\dot{\boldsymbol{\lambda }}}(\tau )=-{\left(\dfrac{\partial {\boldsymbol{H}}}{\partial \boldsymbol{x}}[\boldsymbol{x}(\tau ),\boldsymbol{u}'(\tau ),\boldsymbol{\lambda }(\tau ),{t}_{\text{f}}]\right)}^{\text{T}}\\\dfrac{\partial {\boldsymbol{H}}}{\partial \boldsymbol{u}'}[\boldsymbol{x}(\tau ),\boldsymbol{u}'(\tau ),\boldsymbol{\lambda }(\tau ),{t}_{\text{f}}]={0}\end{cases} $

式中：

$ \boldsymbol{x}=(\boldsymbol{r},\boldsymbol{V},m) $

为状态量；

$ \boldsymbol{f} $

为问题2的运动方程；

$ \boldsymbol{\lambda } $

为协态量；

$ {\boldsymbol{H}}={\boldsymbol{\lambda }}^{{\mathrm{T}}}\boldsymbol{f} $

为最优控制问题的哈密顿函数。

应用邻近最优理论^[21-23]，在最优飞行时间

$ t_{\text{f}}^{k-1} $

和飞行轨迹

$ ({\boldsymbol{r}}^{k-1},{\boldsymbol{V}}^{k-1},{m}^{k-1},{\boldsymbol{u}}^{k-1}) $

附近展开问题2的必要条件，可得到火箭状态修正值

$ {\text{δ}} \boldsymbol{x} $

、协态修正值

$ {\text{δ}} \boldsymbol{\lambda } $

和控制修正值

$ {\text{δ}} \boldsymbol{u}' $

应满足的运动方程。该运动方程可写为紧凑形式：

$ \begin{cases}{\text{δ}} {\dot{\boldsymbol{x}}}(\tau )={\boldsymbol{B}}_{1}(\tau ){\text{δ}} \boldsymbol{x}(\tau )-{\boldsymbol{B}}_{2}(\tau ){\text{δ}} \boldsymbol{\lambda }(\tau )+{\boldsymbol{B}}_{4}(\tau )\text{d}{t}_{\text{f}}\\{\text{δ}} {\dot{\boldsymbol{\lambda }}}(\tau )=-{\boldsymbol{B}}_{3}(\tau ){\text{δ}} \boldsymbol{x}(\tau )-{\boldsymbol{B}}_{1}{}^{\text{T}}(\tau ){\text{δ}} \boldsymbol{\lambda }(\tau )-{\boldsymbol{B}}_{5}(\tau )\text{d}{t}_{\text{f}}\\{\text{δ}} \boldsymbol{u}'(\tau )=-\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}(\tau )\left[{\boldsymbol{H}}_{\boldsymbol{u}'\boldsymbol{x}}^{}(\tau ){\text{δ}} \boldsymbol{x}(\tau )+\right.\\\qquad \left.\boldsymbol{f}_{\boldsymbol{u}'}^{\text{T}}(\tau ){\text{δ}} \boldsymbol{\lambda }(\tau )+{\boldsymbol{H}}_{\boldsymbol{u}'{{t}_{\text{f}}}}(\tau )\text{d}{t}_{\text{f}}\right]\end{cases} $

式中：

$ \text{d}{t}_{\text{f}} $

为飞行时间偏差；带下标的

$ \boldsymbol{f} $

和

$ \boldsymbol{H} $

表示其对各变量的一阶和二阶偏导数，例如：

$ {\boldsymbol{H}}_{\boldsymbol{u}'\boldsymbol{x}}(\tau )=\frac{{\partial }^{2}{\boldsymbol{H}}}{\partial \boldsymbol{u}'\partial \boldsymbol{x}}\left[{\boldsymbol{x}}^{k-1}(\tau ),\boldsymbol{u}{'}^{k-1}(\tau ),{\boldsymbol{\lambda }}^{k-1}(\tau ),t_{\text{f}}^{k-1}\right] $

$ {\boldsymbol{B}}_{i}(i=1,2,\cdots ,5) $

表示以下形式的矩阵：

$ \begin{cases}{\boldsymbol{B}}_{1}(\tau )={\boldsymbol{f}}_{\boldsymbol{x}}(\tau )-{\boldsymbol{f}}_{\boldsymbol{u}'}(\tau )\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}(\tau ){\boldsymbol{H}}_{\boldsymbol{u}'\boldsymbol{x}}(\tau )\\{\boldsymbol{B}}_{2}(\tau )={\boldsymbol{f}}_{\boldsymbol{u}'}(\tau )\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}(\tau )\boldsymbol{f}_{\boldsymbol{u}'}^{\text{T}}(\tau )\\{\boldsymbol{B}}_{3}(\tau )={\boldsymbol{H}}_{\boldsymbol{x}\boldsymbol{x}}(\tau )-{\boldsymbol{H}}_{\boldsymbol{x}\boldsymbol{u}'}(\tau )\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}(\tau ){\boldsymbol{H}}_{\boldsymbol{u}'\boldsymbol{x}}(\tau )\\{\boldsymbol{B}}_{4}(\tau )={\boldsymbol{f}}_{{{t}_{\text{f}}}}(\tau )-{\boldsymbol{f}}_{\boldsymbol{u}'}(\tau )\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}(\tau ){\boldsymbol{H}}_{\boldsymbol{u}{{t}_{\text{f}}}}(\tau )\\{\boldsymbol{B}}_{5}(\tau )={\boldsymbol{H}}_{\boldsymbol{x}{{t}_{\text{f}}}}(\tau )-{\boldsymbol{H}}_{\boldsymbol{x}\boldsymbol{u}'}(\tau )\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}(\tau ){\boldsymbol{H}}_{\boldsymbol{u}'{{t}_{\text{f}}}}(\tau )\end{cases} $

在第

$ k $

个制导周期中，量测得到初始状态偏差

$ {\text{δ}} \boldsymbol{x}_{0}^{k} $

后，若初始协态修正量

$ {\text{δ}} \boldsymbol{\lambda }_{0}^{k} $

和飞行时间修正量

$ \text{d}{t}_{\text{f}} $

已知，则可求解式(19)中的微分方程，从而获得飞行轨迹的修正量

$ {\text{δ}} {\boldsymbol{x}}^{k}(\tau ) $

和

$ {\text{δ}} \boldsymbol{u}{'}^{k}(\tau ) $

。本文应用扫掠法计算初始协态修正量

$ {\text{δ}} \boldsymbol{\lambda }_{0}^{k} $

和飞行时间修正量

$ \text{d}t_{\text{f}}^{k} $

。问题2对应的扫掠变量

$ {\boldsymbol{S}}_{i} $

满足的运动方程为

$ \left[\begin{matrix}{{\dot{\boldsymbol{S}}}}_{1}\\{{\dot{\boldsymbol{S}}}}_{2}\\{{\dot{\boldsymbol{S}}}}_{3}\\{{\dot{\boldsymbol{S}}}}_{4}\\{{\dot{\boldsymbol{S}}}}_{5}\\{{\dot{\boldsymbol{S}}}}_{6}\end{matrix}\right]=\left[\begin{matrix}-{\boldsymbol{B}}_{3}-\boldsymbol{B}_{1}^{\text{T}}{\boldsymbol{S}}_{1}-{\boldsymbol{S}}_{1}{\boldsymbol{B}}_{1}+{\boldsymbol{S}}_{1}{\boldsymbol{B}}_{2}{\boldsymbol{S}}_{1}\\({\boldsymbol{S}}_{1}{\boldsymbol{B}}_{2}-\boldsymbol{S}_{1}^{\text{T}}){\boldsymbol{S}}_{2}\\({\boldsymbol{S}}_{1}{\boldsymbol{B}}_{2}-\boldsymbol{S}_{1}^{\text{T}}){\boldsymbol{S}}_{3}-{\boldsymbol{S}}_{1}{\boldsymbol{B}}_{4}-{\boldsymbol{B}}_{5}\\\boldsymbol{S}_{2}^{\text{T}}{\boldsymbol{B}}_{2}{\boldsymbol{S}}_{2}\\\boldsymbol{S}_{2}^{\text{T}}({\boldsymbol{B}}_{2}{\boldsymbol{S}}_{3}-{\boldsymbol{B}}_{4})\\\boldsymbol{S}_{3}^{\text{T}}({\boldsymbol{B}}_{2}{\boldsymbol{S}}_{3}-{\boldsymbol{B}}_{4})-\boldsymbol{B}_{4}^{\text{T}}{\boldsymbol{S}}_{3}-{\boldsymbol{B}}_{6}\end{matrix}\right] $

式中：

$ {\boldsymbol{B}}_{6}={\boldsymbol{H}}_{{{t}_{\text{f}}}{{t}_{\text{f}}}}-{\boldsymbol{H}}_{{{t}_{\text{f}}}\boldsymbol{u}'}\boldsymbol{H}_{\boldsymbol{u}'\boldsymbol{u}'}^{-1}{\boldsymbol{H}}_{\boldsymbol{u}'{{t}_{\text{f}}}} $

设地月转移入轨条件

$ \boldsymbol{\beta }={0} $

对应的拉格朗日乘子为

$ \boldsymbol{\nu } $

，定义函数：

$ G(\boldsymbol{x},{\boldsymbol{\nu}} )=-m(1)+{\boldsymbol{\nu }}^{\text{T}}\boldsymbol{\beta }\left[\boldsymbol{r}(1),\boldsymbol{V}(1)\right] $

则扫掠变量

$ {\boldsymbol{S}}_{i} $

的终止状态满足以下条件：

$ \begin{cases}{\boldsymbol{S}}_{1}(1)={\boldsymbol{G}}_{\boldsymbol{x}\boldsymbol{x}}\\{\boldsymbol{S}}_{2}(1)=\boldsymbol{\beta }_{\boldsymbol{x}}^{\text{T}}\\{\boldsymbol{S}}_{3}(1)={0}\\{\boldsymbol{S}}_{4}(1)={0}\\{\boldsymbol{S}}_{5}(1)={0}\\{\boldsymbol{S}}_{6}(1)={0}\end{cases} $

应用数值积分方法在时间区间

$ \tau \in [0,1] $

上反向求解式(22)中的状态方程，可获得扫掠变量

$ {\boldsymbol{S}}_{i} $

在

$ \tau =0 $

处的状态。通过求解线性化后的必要条件，初始协态修正量

$ {\text{δ}} \boldsymbol{\lambda }_{0}^{k} $

和飞行时间偏差

$ \text{d}{t}_{\text{f}} $

应满足以下必要条件：

$ {\text{δ}} \boldsymbol{\lambda }_{0}^{k}=\left({\boldsymbol{S}}_{1}-\left[\begin{matrix}{\boldsymbol{S}}_{2} , {\boldsymbol{S}}_{3}\end{matrix}\right]{\left[\begin{matrix}{\boldsymbol{S}}_{4} & {\boldsymbol{S}}_{5}\\\boldsymbol{S}_{5}^{\text{T}} & {\boldsymbol{S}}_{6}\end{matrix}\right]}^{-1}\left[\begin{matrix}\boldsymbol{S}_{2}^{\text{T}}\\\boldsymbol{S}_{3}^{\text{T}}\end{matrix}\right]\right){\text{δ}} \boldsymbol{x}_{0}^{k} $

$ \left[\begin{matrix}\text{d}{\boldsymbol{\nu }}^{k}\\\text{d}t_{\text{f}}^{k}\end{matrix}\right]=-{\left[\begin{matrix}{\boldsymbol{S}}_{4} & {\boldsymbol{S}}_{5}\\\boldsymbol{S}_{5}^{\text{T}} & {\boldsymbol{S}}_{6}\end{matrix}\right]}^{-1}\left[\begin{matrix}\boldsymbol{S}_{2}^{\text{T}}\\\boldsymbol{S}_{3}^{\text{T}}\end{matrix}\right]{\text{δ}} \boldsymbol{x}_{0}^{k} $

式中：为简明起见，

$ {\boldsymbol{S}}_{i} $

均表示

$ {\boldsymbol{S}}_{i}(\tau ) $

在0处的取值；

$ \text{d}{\boldsymbol{\nu }}^{k} $

为拉格朗日乘子

$ \boldsymbol{\nu } $

对应的修正量。

计算得到

$ {\text{δ}} \boldsymbol{\lambda }_{0}^{k} $

和

$ \text{d}{t}_{\text{f}} $

后，可求解(19)中的微分方程，获得状态修正值

$ {\text{δ}} {\boldsymbol{x}}^{k} $

和控制修正值

$ {\text{δ}} \boldsymbol{u}{'}^{k} $

，从而将第

$ k $

个制导周期中的初始猜想设置为

$ \left[\begin{matrix}{\boldsymbol{r}}^{k}(\tau )\\{\boldsymbol{V}}^{k}(\tau )\\{{m}}^{k}(\tau )\end{matrix}\right]=\left[\begin{matrix}{\boldsymbol{r}}^{k-1}(\tau )\\{\boldsymbol{V}}^{k-1}(\tau )\\{m}^{k-1}(\tau )\end{matrix}\right]+{\text{δ}} {\boldsymbol{x}}^{k}(\tau ) $

$ {\boldsymbol{u}}^{k}(\tau )=\left[\begin{matrix}\cos \left[u_{a}^{k-1}(\tau )+{\text{δ}} u_{a}^{k}(\tau )\right]\cos \left[u_{b}^{k-1}(\tau )+{\text{δ}} u_{b}^{k}(\tau )\right]\\\cos \left[u_{a}^{k-1}(\tau )+{\text{δ}} u_{a}^{k}(\tau )\right]\sin \left[u_{b}^{k-1}(\tau )+{\text{δ}} u_{b}^{k}(\tau )\right]\\\sin \left[u_{a}^{k-1}(\tau )+{\text{δ}} u_{a}^{k}(\tau )\right]\end{matrix}\right] $

$ t_{\text{f}}^{k}=t_{\text{f}}^{k-1}+\text{d}t_{\text{f}}^{k} $

式(28)～式(30)中的初始猜想具有以下特点：① 在

$ \tau =0 $

处，初始状态等于火箭当前飞行状态；②满足问题2的线性化一阶必要条件。因此，邻近最优方法提供的初始猜想一方面修正了多源干扰产生初始状态偏差

$ {\text{δ}} \boldsymbol{x}_{0}^{k} $

的影响，另一方面具有近似最优性，可为二阶锥规划提供高质量的初始参考。

4 基于二阶锥规划的轨迹规划方法

收起

考虑问题1中描述的最优控制问题，本节提出一种基于高斯伪谱法和二阶锥规划的轨迹优化方法及一种基于协态映射定理的协态量估计方法。问题1的哈密顿函数可表示为

$ {\boldsymbol{H}}={\boldsymbol{\lambda }}^{\text{T}}\boldsymbol{f}=\boldsymbol{\lambda }_{\boldsymbol{r}}^{\text{T}}{\dot{\boldsymbol{r}}}+\boldsymbol{\lambda }_{\boldsymbol{V}}^{\text{T}}{\dot{\boldsymbol{V}}}+{{\boldsymbol{\lambda}} }_{m}\dot{m} $

式中：

$ \boldsymbol{\lambda }=({\boldsymbol{\lambda }}_{\boldsymbol{r}},{\boldsymbol{\lambda }}_{\boldsymbol{V}},{{\boldsymbol{\lambda}} }_{m}) $

为协态量；

$ {\boldsymbol{\lambda }}_{\boldsymbol{r}} $

、

$ {\boldsymbol{\lambda }}_{\boldsymbol{V}} $

和

$ {{\boldsymbol{\lambda}} }_{m} $

分别为位置协态、速度协态和质量协态。

应用极小值原理，最优控制量

$ {\boldsymbol{u}}^{*} $

的选择应使哈密顿函数

$ H $

最小化，即

$ {\boldsymbol{u}}^{*}(t)=\arg \min {\boldsymbol{H}}=-\frac{{\boldsymbol{\lambda }}_{\boldsymbol{V}}(t)}{\left|\left|{\boldsymbol{\lambda }}_{\boldsymbol{V}}(t)\right|\right|}\qquad {\boldsymbol{\lambda }}_{\boldsymbol{V}}(t)\neq 0 $

注意到上述关系在等式约束

$ \left|\left|\boldsymbol{u}(t)\right|\right|=1 $

或不等式约束

$ \left|\left|\boldsymbol{u}(t)\right|\right|\leqslant 1 $

下均成立。因此，可应用无损凸化方法将问题1中的等式约束松弛为不等式约束，且不改变其最优解的特性，相关理论分析与证明见文献[24]。与问题2相同，同样将时间区间从

$ t\in [{t}_{0},{t}_{\text{f}}] $

变换为

$ \tau \in [0,1] $

，得到如下形式的最优控制问题：

问题3　固定时间区间的地月转移入轨问题。

$ \begin{cases}\operatorname{minimize} \quad m(0)-m(1) \\\text {subject to} \quad \dot{\boldsymbol{x}}(\tau)=\left(t_{\mathrm{f}}-t_0\right) \boldsymbol{f}[\boldsymbol{x}(\tau), \boldsymbol{u}(\tau)] \\\|\boldsymbol{u}(\tau)\| \leqslant 1 \\\boldsymbol{r}(0)=\boldsymbol{r}_0\\\boldsymbol{V}(0)=\boldsymbol{V}_0\\m(0)=m_0 \\\boldsymbol{\beta}\left[\boldsymbol{r}\left(t_{\mathrm{f}}\right), \boldsymbol{V}\left(t_{\mathrm{f}}\right)\right]=0\end{cases} $

在问题3中，初始时间

$ {t}_{0} $

为常量，终止时间

$ {t}_{\text{f}} $

为待优化的问题参数。应用序列线性化方法，可将问题3转化为一系列的二阶锥规划问题，并应用内点法实时求解。

在邻近最优方法提供的初始参考轨迹

$ \left[{\boldsymbol{x}}^{k}(\tau ),\right. \left.{\boldsymbol{u}}^{k}(\tau ),t_{\text{f}}^{k}\right] $

附近将运动方程线性展开，得到以下形式的线性化运动方程：

$ {\dot{\boldsymbol{x}}}(\tau )=\boldsymbol{f}_{\text{linear}}^{k}[\boldsymbol{x}(\tau ),\boldsymbol{u}(\tau ),{t}_{\text{f}}] $

应用高斯伪谱法离散化问题3，离散点数量为

$ N $

，则各离散点

$ {\tau }_{i} $

的位置为

$ {\tau }_{i}=\left\{\begin{matrix}0 & i=1 & \\1 & i=N & \\({P}_{N,i-1}+1)/2 & \text{其他} & \end{matrix}\right. $

式中：

$ {P}_{N,i-1} $

为

$ [-1,1] $

区间上的勒让德-高斯点。

状态量

$ \boldsymbol{x}(\tau ) $

使用前

$ N-1 $

个点组成的拉格朗日插值多项式

$ {l}_{i} $

近似表示：

$ \boldsymbol{x}(\tau )=\sum \limits_{i=1}^{N-1}{l}_{i}(\tau ){\boldsymbol{x}}_{i} $

控制量

$ \boldsymbol{u}(\tau ) $

使用中间

$ N-2 $

个点组成的拉格朗日插值多项式近似表示：

$ \boldsymbol{u}(\tau )=\sum \limits_{i=2}^{N-1}{l}_{i}(\tau ){\boldsymbol{u}}_{i} $

终端状态可使用高斯积分公式获得：

$ {\boldsymbol{x}}_{N}={\boldsymbol{x}}_{0}+\frac{1}{2}\sum \limits_{k=2}^{N-1}{w}_{k}\boldsymbol{f}_{\text{linear}}^{k}({\boldsymbol{x}}_{k},{\boldsymbol{u}}_{k}) $

式中：

$ {w}_{k} $

为高斯伪谱法的积分系数。

令运动方程在中间

$ N-2 $

个点上成立，建立非线性方程组：

$ \sum \limits_{j=1}^{N-1}{\dot{l}}_{j}({\tau }_{i}){\boldsymbol{x}}_{j}=({t}_{\text{f}}-{t}_{0})\boldsymbol{f}_{\text{linear}}^{k}({\boldsymbol{x}}_{i},{\boldsymbol{u}}_{{i}})\qquad i=2,3,\cdots ,N-1 $

将式(20)和式(21)中的非线性方程组写为紧凑形式：

$ {\boldsymbol{X}}_{N}={\boldsymbol{D}}_{N}{\boldsymbol{F}}_{N} $

式中：

$ {\boldsymbol{X}}_{N} $

为状态量组成的矩阵；

$ {\boldsymbol{F}}_{N} $

为状态方程组成的矩阵；

$ {\boldsymbol{D}}_{N} $

为高斯伪谱法的微分方程。

至此，建立一个二阶锥规划问题：

问题4　离散化后的地月转移入轨问题。

$ \begin{cases}\operatorname{minimize} \quad m_1-m_N \\\text { subject to } \quad \boldsymbol{X}_N=\boldsymbol{D}_N \boldsymbol{F}_N \\\left\|\boldsymbol{u}_k\right\| \leqslant 1 \\\boldsymbol{r}_1=\boldsymbol{r}_0\\\boldsymbol{V}_1=\boldsymbol{V}_0\\m_1=m_0 \\\boldsymbol{\beta}_{\text {linear }}\left(\boldsymbol{r}_N, \boldsymbol{V}_N\right)=0\end{cases} $

式中：

$ {\boldsymbol{\beta }}_{\text{linear}} $

为线性化地月转移入轨条件。

求解问题4可以获得新的轨迹

$ \left[{\boldsymbol{x}}^{k}(\tau ),{\boldsymbol{u}}^{k}(\tau ),t_{\text{f}}^{k}\right] $

，通过迭代方式，其收敛后的解将满足问题3的最优必要条件，可作为制导指令使用。

本文应用基于协态映射定理的协态量估计方法，能够基于线性化矩阵方程

$ {\boldsymbol{X}}_{N}={\boldsymbol{D}}_{N}{\boldsymbol{F}}_{N} $

对应的拉格朗日乘子矩阵

$ \boldsymbol{\varLambda } $

估计协态量

$ \boldsymbol{\lambda }(\tau ) $

。具体地，协态量

$ \boldsymbol{\lambda }(\tau ) $

在中间

$ N-2 $

的离散点处的估计值为

$ {\boldsymbol{\lambda }}_{k}={\boldsymbol{\varLambda }}_{k}/{w}_{k}+{\boldsymbol{\varLambda }}_{N}\qquad k=2,3,\cdots ,N-1 $

协态量在最后一个离散点处的估计值为

$ {\boldsymbol{\lambda }}_{N}={\boldsymbol{\varLambda }}_{N} $

随后，可应用拉格朗日插值方法估计协态量的连续形式：

$ \boldsymbol{\lambda }(\tau )=\sum \limits_{k=2}^{N}{l}_{k}(\tau ){\boldsymbol{\lambda }}_{k} $

协态映射定理估计出的协态量

$ \boldsymbol{\lambda }(\tau ) $

为邻近最优方法的应用提供了必要基础。

迭代终止条件设计为迭代间终止时间

$ {t}_{\text{f}} $

和终止速度

$ {\boldsymbol{V}}_{N} $

的相对变化：

$ {\delta }_{\text{opt}}=\min \left\{\frac{\left| {t}_{\text{f}}-t_{\text{f}}^{\text{pre}}\right| }{t_{\text{f}}^{\text{pre}}},\frac{\left|\left|{\boldsymbol{V}}_{N}-\boldsymbol{V}_{N}^{\text{pre}}\right|\right|}{\boldsymbol{V}_{N}^{\text{pre}}}\right\} $

式中：上标pre表示上一次迭代。

当

$ {\delta }_{\text{opt}} $

小于预设的阈值

$ {\epsilon }_{\text{opt}} $

时，二阶锥规划迭代停止，输出制导指令。综上所述，提出的地月转移入轨混合制导方法的算法流程如下：

输入：火箭飞行状态

$ \boldsymbol{r} $

、

$ \boldsymbol{V} $

、

$ m $

，目标轨道参数

$ {\boldsymbol{L}}_{\text{ref}} $

和

$ {\boldsymbol{V}}_{\text{ref}} $

。

输出：推力方向指令

$ \boldsymbol{u} $

。

初始化：

1. 应用线性插值方法，生成状态量和控制量的初始参考

$ ({\boldsymbol{r}}^{0},{\boldsymbol{V}}^{0},{m}^{0},{\boldsymbol{u}}^{0}) $

；

2. 将入轨时间

$ t_{\text{f}}^{0} $

的初始参考设为标称入轨时间。

每个制导周期：

3. 从导航系统获取当前的飞行时间

$ {t}_{0} $

和状态

$ {\boldsymbol{r}}_{0} $

、

$ {\boldsymbol{V}}_{0} $

和

$ {m}_{0} $

；

4. 对第一个制导周期

$ k=1 $

，转步骤10。否则，转步骤5；

5. 应用式(14)计算初始状态偏差

$ {\text{δ}} \boldsymbol{x}_{0}^{k} $

；

6. 求解式(22)和式(25)描述的微分方程，获得

$ {\boldsymbol{S}}_{i}(0) $

；

7. 应用式(26)和式(27)计算初始协态修正量

$ {\text{δ}} \boldsymbol{\lambda }_{0}^{k} $

和飞行时间修正量

$ \text{d}{t}_{\text{f}} $

；

8. 求解式(19)描述的微分方程，获得状态修正值

$ {\text{δ}} {\boldsymbol{x}}^{k} $

和控制修正值

$ {\text{δ}} \boldsymbol{u}{'}^{k} $

；

9. 应用式(28)～式(30)修正初始参考；

10. 基于初始参考线性化问题3获得问题4；

11. 应用二阶锥规划方法求解问题4，并应用式(42)～式(44)计算协态量；

12. 应用式(45)计算判据

$ {\delta }_{\text{opt}} $

；

13. 若

$ {\delta }_{\text{opt}} $

小于预设的阈值

$ {\epsilon }_{\text{opt}} $

，转步骤14；否则，更新初始参考后转步骤10；

14. 输出制导指令

$ \boldsymbol{u} $

，转步骤3。

终止条件：

15. 当待飞时间

$ {t}_{\text{f}}-{t}_{0}\leqslant {t}_{\text{stop}} $

时，制导迭代停止；

当半长轴达到目标时，关闭发动机完成入轨。

5 数值实验与分析

收起

本节应用文献[7]中提供的地月转移工况，验证本文混合制导方法的性能，共进行了2类数值实验：与现有方法的对比及对动力故障的适应能力。

在2类数值实验中，制导方法的参数保持一致，具体包括：制导周期为2 s，制导停止时间

$ {t}_{\text{stop}} $

为60 s；邻近最优中求解微分方程的方法为定步数四阶龙格-库塔方法，步数为50步；二阶锥规划中高斯伪谱法的离散点个数

$ N $

为10，迭代停止阈值

$ {\epsilon }_{\text{opt}} $

为

$ {10}^{-9} $

。混合制导方法使用C语言实现，运行环境为主频为1 GHz的Arm处理器。

5.1 与现有方法的对比

为验证本文混合制导方法对地月转移入轨任务的适应能力，在推力和秒耗量同时下降35%且存在2%未知偏差（

$ \delta_T=\delta_R=0.02 $

）的情况下，对比了本文混合制导方法与文献[2,13]中为SLS Block 1开发的改进PEG方法，结果如图3和图4所示。

在数值仿真中，改进PEG与本文混合制导方法均完成了入轨任务。混合制导方法的入轨时间为651.692 s。在2%的推力和秒耗量偏差下，偏心率偏差为0.000147，近地点幅角偏差为−0.004154°（以半长轴作为入轨条件，故半长轴无偏差）。与离线轨迹规划工具箱优化得到的最优参考对比，燃料消耗相对误差为5.9 kg。与之相对的，改进PEG的燃料消耗相对误差为30.9 kg，较本文混合制导方法消耗更多燃料。上述结果表明，本文混合制导方法具有较好的入轨精度和燃料最优性，能够满足一些工程应用的需求。

为验证本文混合制导方法综合应用2种轨迹规划方法的意义，与纯二阶锥规划方法（即不使用邻近最优修正初始参考）进行了对比。混合制导方法和纯二阶锥规划方法均能完成上述地月转移入轨任务，精度在同一数量级。两者的计算时间如图5所示。

本文混合制导方法在296个制导周期中的平均计算时间为75.353 ms，与纯二阶锥规划方法（108.691 ms）相比，平均计算时间下降30.67%。表明本文混合制导策略有效，应用高计算效率的邻近最优方法降低多源干扰影响后，二阶锥规划方法能够快速迭代收敛。这一方面能够降低箭载计算机的计算负载，另一方面允许更短的制导周期，从而进一步提升对多源干扰的适应能力。

5.2 对动力故障的适应能力

为验证本文混合制导方法对动力故障的适应能力，在推力和秒耗量下降40%、45%、50%的3种工况且存在未知2%偏差的情况下进行了仿真测试，结果如图6、图7和表1所示。

本文混合制导方法在3种不同动力故障下均能够完成制导任务，表1中的各项轨道要素偏差说明入轨精度无明显差异。与离线轨迹规划得到的最优参考对比，推进剂消耗误差小于18.200 kg。改进PEG在这3种动力故障下无法完成制导任务，在推力下降40%和45%的情况下，改进PEG分别于652 s和546 s处迭代发散，单制导周期内迭代20次仍未收敛，在推力50%的情况下，初始迭代不能收敛。上述结果表明，本文混合制导方法对该工况下的动力故障具有较强的适应能力。在现有制导方法无法适应40%及以上推力下降的情况下，本文混合制导方法能够适应50%的推力下降。

6 结　论

收起

1）基于实时轨迹规划的混合轨迹规划方法成功解决了飞行时间在650 s左右的小推重比地月转移入轨问题，且燃料消耗较现有制导方法低。

2）本文提出的邻近最优校正方法有效降低了轨迹优化的迭代次数。使用邻近最优校正后，制导计算时间从108.691 ms降低至75.353 ms，降幅达到30.67%。

3）本文混合制导方法的计算时间满足现有制导周期（2 s）的需求，且存在改善空间。实际应用中可选择更短制导周期，从而进一步提高对多源干扰的适应能力。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

包为民. 航天智能控制技术让运载火箭 “会学习”[J]. 航空学报, 2021, 42(11): 1-10.

Bao W M. Space intelligent control technology enables launch vehicle to “self-learning”[J]. Acta Aeronautica et Astronautica Sinica, 2021, 42(11): 1-10(in Chinese).

Von Der Porten P, Ahmad N, Hawkins M, et al. Powered explicit guidance modifications and enhancements for space launch system Block-1 and Block-1B vehicles[C]//Proceedings of the 41st American Astronautical Society Guidance, Navigation, and Control Conference. Breckenridge: AAS, 2018: 136.

Chandler D C, Smith I E. Development of the iterative guidance mode with its application to various vehicles and missions[J]. Journal of Spacecraft and Rockets, 1967, 4(7): 898-903.

吕新广, 宋征宇. 载人运载火箭迭代制导方法应用研究[J]. 载人航天, 2009, 15(1): 9-14.

Lyu X G, Song Z Y. Study of the iterative guidance engineering application to manned launch vehicle[J]. Manned Spaceflight, 2009, 15(1): 9-14(in Chinese).

茹家欣. 液体运载火箭的一种迭代制导方法[J]. 中国科学(E辑), 2009(4): 696-706.

Ru J X. An iterative guidance method for liquid launch vehicle[J]. Science in China (Series E), 2009(4): 696-706(in Chinese).

施国兴, 吕新广, 巩庆海. 满足多终端约束的二次曲线迭代制导方法研究[J]. 中国空间科学技术, 2018, 38(2): 24-31.

Shi G X, Lyu X G, Gong Q H. Research on quadratic curve IGM for multi-terminal constraints[J]. Chinese Space Science and Technology, 2018, 38(2): 24-31(in Chinese).

马宗占, 许志, 唐硕, 等. 一种改进的运载火箭迭代制导方法[J]. 航空学报, 2021, 42(2): 217-229.

Ma Z Z, Xu Z, Tang S, et al. Improved iterative guidance method for launch vehicles[J]. Acta Aeronautica et Astronautica Sinica, 2021, 42(2): 217-229(in Chinese).

张珂洋, 王鹏, 于琦, 等. 运载火箭助推段大攻角剖面高精度弹道解析解[J]. 北京航空航天大学学报, 2025, 51(10): 3524-3534.

Zhang K Y, Wang P, Yu Q, et al. Analytical solutions for ascent trajectory under high angle-of-attack conditions[J]. Journal of Beijing University of Aeronautics and Astronautics, 2025, 51(10): 3524-3534(in Chinese).

Ahmad N, Hawkins M, Von Der Porten P, et al. Closed loop guidance trade study for space launch system Block-1B vehicle[C]//Proceedings of the AAS/AIAA Astrodynamics Specialist Conference. San Diego: Univelt, 2019.

王大轶, 符方舟, 孟林智, 等. Research of Autonomous Control Technology for Deep Space Probes[J]. 深空探测学报, 2019, 6(4): 317-327.

Wang D Y, Fu F Z, Meng L Z, et al. Research of autonomous control technology for deep space probes[J]. Journal of Deep Space Exploration, 2019, 6(4): 317-327(in Chinese).

巩庆海, 宋征宇, 吕新广. 迭代制导在月面上升段的应用研究[J]. 载人航天, 2015, 21(3): 231-236.

Gong Q H, Song Z Y, Lyu X G. Study on application of iterative guidance in lunar ascent[J]. Manned Spaceflight, 2015, 21(3): 231-236(in Chinese).

宋春, 郭延宁, 郭敏文, 等. 火星上升器入轨段自适应迭代制导方案设计[J]. 深空探测学报(中英文), 2024, 11(1): 31-39.

Song C, Guo Y N, Guo M W, et al. Design of adaptive iterative guidance scheme for Mars ascent vehicle orbiting phase[J]. Journal of Deep Space Exploration, 2024, 11(1): 31-39(in Chinese).

Mchenry R L, Long A L, Cockrell B F, et al. Space shuttle ascent guidance, navigation, and control[J]. Journal of the Astronautical Sciences, 1979, 27: 1-38.

Jaggers R F. Proposed powered explicit guidance thrust integrals derivation/implementation. Mission planning, mission analysis and software formulation[R]. Washington, D. C. : NASA, 1976.

GOODMAN J. Roland jaggers and the development of space shuttle powered explicit guidance (PEG)[C]//Proceedings of the AIAA Scitech 2021 Forum. Reston: AIAA, 2021: 2021.

Hawkins M, Ahmad N. Guidance modifications and enhancements for space launch system Block-1 in support of Artemis I and beyond[C]//Proceedings of the 2020 AAS/AIAA Astrodynamics Specialist Conference. San Diego: Univelt, 2021.

张洪华, 李骥, 于萍, 等. 嫦娥五号月面起飞上升制导导航与控制技术[J]. 中国科学(技术科学), 2021, 51(8): 921-937.

Zhang H H, Li J, Yu P, et al. Guidance navigation and control technology for the lunar ascent vehicle of the Chang’e-5 mission[J]. SCIENTIA SINICA Technologica, 2021, 51(8): 921-937(in Chinese).

Brown K R, Harrold E F, Johnson G W. Rapid optimization of multiple-burn rocket flights: NASA CR-1430[R]. Washington, D. C. : NASA, 1969.

Lu P, Griffin B J, Dukeman G A, et al. Rapid optimal multiburn ascent planning and guidance[J]. Journal of Guidance, Control, and Dynamics, 2008, 31(6): 1656-1664.

陈虹坤, 谭述君, 崔梦真, 等. 大气层外运载火箭上升段轨迹规划的保辛伪谱序列凸优化方法[J]. 北京航空航天大学学报, 2026, 52(4): 1211-1220.

Chen H K, Tan S J, Cui M Z, et al. A symplectic pseudo-spectral successive convex optimization method for trajectory planning of ascent stage of exo-atmosphere launch vehicle[J]. Journal of Beijing University of Aeronautics and Astronautics, 2026, 52(4): 1211-1220(in Chinese).

张志国, 余梦伦, 耿光有, 等. 应用伪谱法的运载火箭在线制导方法研究[J]. 宇航学报, 2017, 38(3): 262-269.

Zhang Z G, Yu M L, Geng G Y, et al. Research on application of pseudo-spectral method in online guidance method for a launch vehicle[J]. Journal of Astronautics, 2017, 38(3): 262-269(in Chinese).

Cheng X M, Li H F, Zhang R. Efficient ascent trajectory optimization using convex models based on the Newton-Kantorovich/Pseudospectral approach[J]. Aerospace Science and Technology, 2017, 66: 140-151.

张永华, 王铮, 王时雨, 等. 基于无迹模型预测静态规划的可回收火箭着陆段鲁棒制导方法[EB/OL]. (2024-09-14)[2025-05-28]. https://doi.org/10.13700/j.bh.1001-5965.2024.0491.

Zhang Y H, Wang Z, Wang S Y, et al. Robust guidance method for reusable rocket landing phase based on UMPSP[EB/OL]. (2024-09-14)[2025-05-28]. https://doi.org/10.13700/j.bh.1001-5965.2024.0491(in Chinese).

Hao Z M, Zhang R. Onboard real-time generation of launch vehicle abort orbits[J]. Journal of Guidance, Control, and Dynamics, 2021, 44(8): 1541-1549.

Wang J W, Zhang R, Li H F. Onboard optimization of multi-arc trajectories with constraints on duration of arcs[J]. Acta Astronautica, 2022, 192: 434-442.

2026年第52卷第10期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.13700/j.bh.1001-5965.2025.0353

接收时间：2025-06-06
首发时间：2026-10-30
出版时间：2026-10-31

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-06-06
录用日期：2025-12-05

基金

作者信息

¹北京航空航天大学宇航学院，北京 100191

²北京航空航天大学天地往返高效运输技术全国重点实验室，北京 100191

³中国航天科技集团有限公司，北京 100048

通讯作者:

E-mail：lihuifeng@buaa.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/bhxbzk/CN/10.13700/j.bh.1001-5965.2025.0353

分享至

全文二维码

扫描看全文

引用本文

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