北京航空航天大学学报

状态约束	速度/ (m·s⁻¹)	航迹角/ (°)	航向角/ (°)	航程/ km	高度/ km	侧向位置/km
初始	1 870	2	−40	0	32	0
终端	1 647	−5	−48	35	31	35

状态约束	速度/ (m·s⁻¹)	航迹角/ (°)	航向角/ (°)	航程/ km	高度/ km	侧向位置/km
初始	1 870	2	−40	0	32	0
终端	1 647	−5	−48	35	31	35

终端状态方法	速度误差/(m·s⁻¹)	航迹角误差/(°)	航向角误差/(°)	航程误差/m	高度误差/m	侧向位置误差/m
DCHU-SC (84)	−2.90×10⁻³	−4.90×10⁻⁵	−2.60×10⁻⁴	0.52	0.14	0.46
MPSC (3 000)	0.10	−5.43×10⁻⁴	1.96×10⁻³	−0.41	−1.34	−0.10
MPSC (84)	−7.11	0.05	0.24	−76.23	−52.18	26.95
GPOPS(84)	−0.02	−4.32×10⁻⁶	−4.18×10⁻⁶	0.54	0.10	0.54
SCP (84)	0.02	0.02	0.03	16.25	2.12	3.26

终端状态方法	速度误差/(m·s⁻¹)	航迹角误差/(°)	航向角误差/(°)	航程误差/m	高度误差/m	侧向位置误差/m
DCHU-SC (84)	−2.90×10⁻³	−4.90×10⁻⁵	−2.60×10⁻⁴	0.52	0.14	0.46
MPSC (3 000)	0.10	−5.43×10⁻⁴	1.96×10⁻³	−0.41	−1.34	−0.10
MPSC (84)	−7.11	0.05	0.24	−76.23	−52.18	26.95
GPOPS(84)	−0.02	−4.32×10⁻⁶	−4.18×10⁻⁶	0.54	0.10	0.54
SCP (84)	0.02	0.02	0.03	16.25	2.12	3.26

基于控制向量参数化的临近空间飞行器轨迹优化

PDF下载

王远卓 , 代洪华 * , 岳晓奎

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

收起

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

基于控制向量参数化的临近空间飞行器轨迹优化

全屏

[Author(id=1322941426715210419, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, orderNo=0, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=null, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1322941426799096501, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, authorId=1322941426715210419, language=EN, stringName=Yuanzhuo WANG, firstName=Yuanzhuo, middleName=null, lastName=WANG, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1322941426895565494, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, authorId=1322941426715210419, language=CN, stringName=王远卓, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=西北工业大学航天学院，西安 710072, bio={"content":"

王远卓男，硕士，博士研究生。主要研究方向：临近空间飞行器轨迹优化与制导

"}, bioImg=null, bioContent=

王远卓男，硕士，博士研究生。主要研究方向：临近空间飞行器轨迹优化与制导

, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1322941426551632559, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, xref=null, ext=[AuthorCompanyExt(id=1322941426560021168, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, companyId=1322941426551632559, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China), AuthorCompanyExt(id=1322941426568409777, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, companyId=1322941426551632559, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=西北工业大学航天学院，西安 710072)])]), Author(id=1322941426983645880, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, orderNo=1, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=hhdai@nwpu.edu.cn, emailSecond=null, emailThird=null, correspondingAuthor=1, authorType=1, ext={EN=AuthorExt(id=1322941427067531962, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, authorId=1322941426983645880, language=EN, stringName=Honghua DAI, firstName=Honghua, middleName=null, lastName=DAI, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=*, address=School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1322941427138835131, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, authorId=1322941426983645880, language=CN, stringName=代洪华, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=*, address=西北工业大学航天学院，西安 710072, bio={"content":"

代洪华男，博士，教授。主要研究方向：空天飞行器动力学与制导控制

"}, bioImg=null, bioContent=

代洪华男，博士，教授。主要研究方向：空天飞行器动力学与制导控制

, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1322941426551632559, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, xref=null, ext=[AuthorCompanyExt(id=1322941426560021168, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, companyId=1322941426551632559, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China), AuthorCompanyExt(id=1322941426568409777, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, companyId=1322941426551632559, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=西北工业大学航天学院，西安 710072)])]), Author(id=1322941427201749693, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, orderNo=2, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=null, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1322941427285635775, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, authorId=1322941427201749693, language=EN, stringName=Xiaokui YUE, firstName=Xiaokui, middleName=null, lastName=YUE, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1322941427377910464, tenantId=1146029695717560320, journalId=1322258626896822465, articleId=1322941421258420849, authorId=1322941427201749693, language=CN, stringName=岳晓奎, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=西北工业大学航天学院，西安 710072, bio={"content":"

岳晓奎男，博士，教授。主要研究方向：空天飞行器总体设计、航天动力学与控制

"}, bioImg=null, bioContent=

岳晓奎男，博士，教授。主要研究方向：空天飞行器总体设计、航天动力学与控制

王远卓, 代洪华*, 岳晓奎

作者信息

西北工业大学航天学院，西安 710072

通讯作者:

E-mail：hhdai@nwpu.edu.cn

作者简介:

王远卓男，硕士，博士研究生。主要研究方向：临近空间飞行器轨迹优化与制导

代洪华男，博士，教授。主要研究方向：空天飞行器动力学与制导控制

岳晓奎男，博士，教授。主要研究方向：空天飞行器总体设计、航天动力学与控制

Near-space vehicle trajectory optimization based on control vector parameterization

Yuanzhuo WANG, Honghua DAI*, Xiaokui YUE

Affiliations

School of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China

Corresponding Author:

E-mail: hhdai@nwpu.edu.cn

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

文章导航

摘要

收起

针对计算制导场景下临近空间飞行器的在线轨迹优化问题，提出了一种基于实际动力学方程控制向量参数化的轨迹优化方法。通过翻转Radau伪谱离散非线性动力学方程，建立了含初始和终端状态约束的轨迹优化问题。通过参数化控制向量，将各离散点处控制量的优化问题转化为几个参数向量的优化问题，简化了计算过程，并提升了计算效率。此外，针对参数幅值约束，融合凸规划技术提出了一种基于控制向量参数化的轨迹优化方法。考虑到中制导任务对终端飞行时间不敏感，将优化问题的自变量从时间转换为航程，进一步缩短了计算时间。仿真结果表明，所提方法能够有效满足计算制导中高精度高效率轨迹优化的任务需求。

关键词

临近空间飞行器 / 最优控制 / 轨迹优化 / 凸规划 / 控制向量参数化设计

Abstract

收起

For online trajectory optimization of the near-space vehicle in computational guidance, a trajectory optimization method based on actual dynamic equations is proposed via control vector parameterization. Firstly, a trajectory optimization problem with initial and terminal state constraints is developed by discretizing the nonlinear dynamic equations using the flipped-Radau pseudospectral approach. Secondly, control vector parameterization reduces complexity and boosts computational efficiency by converting the optimization problem for control quantities at each discrete point into an optimization problem for multiple parameter vectors. Thirdly, to deal with parameter amplitude constraints, a trajectory optimization method based on convex programming techniques is proposed via control vector parameterization. Finally, given that the midcourse guidance mission is insensitive to terminal flight time, the optimization problem is transformed from being time-dependent to downrange-dependent, thereby further reducing the computational time. Simulation results demonstrate that the proposed method satisfies the high-efficiency and high-precision requirements for trajectory optimization in computational guidance.

Key words

near-space vehicle / optimal control / trajectory optimization / convex programming / control vector parametric design

引用本文

王远卓, 代洪华, 岳晓奎. 基于控制向量参数化的临近空间飞行器轨迹优化. 北京航空航天大学学报, 2026 , 52 (10) : 3602 -3617 . DOI: 10.13700/j.bh.1001-5965.2025.0247

Yuanzhuo WANG, Honghua DAI, Xiaokui YUE. Near-space vehicle trajectory optimization based on control vector parameterization[J]. Journal of Beijing University of Aeronautics and Astronautics, 2026 , 52 (10) : 3602 -3617 . DOI: 10.13700/j.bh.1001-5965.2025.0247

正文

收起

随着科学技术的迅猛发展，平台体系正加速迈向机械化、信息化与智能化深度融合的崭新阶段。未来，主战平台的战略定位、技术形态以及作战运用模式将迎来革命性的变革^[1]。在进攻型空战场景下，有人/无人异构协同战斗机群需要建立区域制空，保护编队后方的战略/战术轰炸机实现载荷（临近空间飞行器）投放并精确打击高价值目标^[2]。本文聚焦于临近空间滑翔飞行器的制导律研究，致力于为其飞行性能与攻击精度的提升提供坚实的理论支撑与创新的技术方案。

临近空间滑翔飞行器具有飞行速度快、突防能力强、航程远、气动特性复杂等特点^[3]。围绕突防环境下飞行器的制导与轨迹规划，闫循良等^[4]系统梳理了现有轨迹规划方法，并深入分析了不同方法的优缺点。围绕多临近空间飞行器“探测-打击-评估”协同任务需求，刘哲等^[5]实现了时间协同约束下的轨迹规划任务。在助推及滑翔段分别以俯仰角及攻角为控制量，An等^[6]提出了一种面向临近空间飞行器弹道设计与优化的框架，并讨论了设计与优化结果之间的差异。在协同制导框架下，Yu等^[7]提出了一种融合多飞行约束的高超声速滑翔飞行器协同制导策略。

临近空间飞行器的初、末制导律一般分别采用比较成熟的方案制导与比例制导。中制导律的研究方案如下：基于自适应动态规划的中制导律^[8]、基于在线轨迹优化的中制导律^[9]、基于微分对策的中制导律^[10]、数据驱动的中制导律^[11]、基于启发式方法的中制导律^[12]。未来的中制导任务将面临更高的需求：不仅需要精确控制飞行器的终端位置，还需控制终端速度的大小和方向。这要求制导系统必须通过法向与侧向过载2个控制量实现6个状态的精确控制。在机载计算资源受限的约束下，针对该典型欠驱动制导问题^[13]，飞行器一般通过在线解算轨迹优化问题实现各终端状态的精确控制^[14]。因此，本文主要研究中制导的在线轨迹优化方法。

近年来，许多学者对各种轨迹优化方法进行了系统的综述，并深入分析了它们的优缺点^[15-16]。轨迹优化方法主要分为直接法与间接法^[17]。间接法基于最优控制理论中的哈密尔顿函数和一阶必要条件，构建两点边值问题，并通过解算该问题实现轨迹优化的目标^[18]。基于间接打靶轨迹优化架构，刘宇航等^[19]融合人工智能机器学习技术实现了深空探测的高精度高效率轨迹优化。Nurre和Taheri^[20]提出了一种改进的间接优化方法，用于生成固定时间下燃料最优的日心轨迹。间接轨迹优化方法具有解算精度高、满足一阶必要条件等优点，但是存在收敛半径小与计算效率低等缺点^[21]。与间接法不同的是，直接法首先将动力学方程离散化，而后通过非线性规划方法解算^[22]。面向大气层边缘高超声速飞行器补能-跳跃滑翔飞行，雍恩米等^[23]基于Gauss伪谱实现了分段轨迹设计。赵吉松等^[24]融合配点法实现了临近空间飞行器助推-滑翔段的最优轨迹设计，并分析了典型参数的敏感性。尽管上述直接轨迹优化方法能够获得精确的数值解，但通常依赖于离线计算，难以满足轨迹优化问题在线解算的实时性需求。

Lu将轨迹优化的在线解算问题定义为计算制导问题^[14,25]。近20年来，出现了3种在线解算方法：模型预测控制（model predictive control，MPC）、凸规划与模型预测静态规划（model predictive static programming，MPSP）^[26]。MPC广泛应用于航空航天制导控制领域，是一种基于系统模型的先进控制策略，通过在每个控制周期内预测系统的未来行为，优化控制输入保证预测输出满足约束条件并最小化性能指标，从而实现对复杂动力学系统的精确控制。Patel和Subbarao^[27]融合可达性分析与MPC，保证了干扰和不确定条件下高超声速飞行器再入过程的安全性。面向控制能力更强的带襟翼火星登陆器，Engel等^[28]融合MPC提升了终端状态的控制精度。Açıkmeşe及其团队成员对凸规划理论的发展做出了巨大贡献^[29]。该方法借助梯度下降快速收敛到全局最优解，是一种广泛使用的轨迹优化方法。通过凸规划技术，Hayner等^[30]成功实现了视线约束下六自由度无人机轨迹优化问题的快速精确解算。针对小行星着陆任务，Ma等^[31]将深度神经网络与序列凸规划相结合，巧妙消除了雅可比矩阵复杂的计算过程，提出了一种基于机器学习的轨迹优化方法。Benedikter等^[32]基于虚拟控制量与自适应信赖域，依托融合无损凸化和序列凸化技术，构建了一种适用于火箭上升段的高效凸规划轨迹优化架构。针对临近空间飞行器在线制导问题，Wang和Grant^[33]基于现有序列凸规划技术，通过二次规划方法设计最优反馈制导律，实现了参考轨迹的精确跟踪。

Padhi及其团队成员提出了基于偏差动力学方程的MPSP方法^[26]。基于参数化后的控制向量，Dwivedi等^[34]提出了一种模型预测扩展控制（model predictive spread control，MPSC）方法。基于无迹MPSP方法，张永华等^[35]提出了一种风场干扰下用于火箭着陆的鲁棒制导方法。融合谱形式的控制量与小维权重矩阵，谭一廷等^[36]针对防空反导导弹的中制导律设计，提出了一种广义准谱MPSP方法。基于临近空间滑翔飞行器俯仰力矩瞬时平衡假设，Ma等^[37]将这一动力学条件转化为路径约束，进而提出了一种考虑俯仰运动准静态假设的快速MPSP轨迹优化方法。对于临近空间飞行器快速轨迹优化，Wang和Dai^[9]基于信赖域约束提出了一种融合翻转Radau伪谱法（flipped Radau pseudospectral method，FRPM）的改进MPSP方法。2019年，Pan等^[11]证明了上述MPSP方法是一种基于实际动力学方程的牛顿类轨迹优化方法。MPSP方法采用解析递推架构，无需借助非线性规划方法，即可实现轨迹优化问题的快速准确解算。然而，该方法的计算时间受到离散点数量的影响。如果离散点数量不足，可能导致终端状态误差较大，从而影响优化结果的准确性。

受上述MPSP方法启发，笔者提出了基于实际动力学方程控制历史直接更新的准谱轨迹优化方法（quasi-spectral trajectory optimization method via direct control history update，QS-DCHU）^[38]。受上述MPSC与DCHU方法启发，本文进一步提出了一种面向计算制导基于实际动力学方程控制向量参数化的临近空间飞行器轨迹优化方法（direct control history update with spread control，DCHU-SC）。本文的创新点主要归纳为以下2点：

1）通过参数化控制向量，将各离散点处的控制量优化问题转化为参数优化问题，降低了问题的复杂度，并提升了计算效率。此外，与大多数航空航天轨迹优化问题类似，临近空间飞行器的中制导任务对终端飞行时间并不敏感。因此，本文将以时间为自变量的轨迹优化问题转化为以航程为自变量的优化问题，使状态量的个数从6个降低至5个，有效降低了欠驱动程度，进一步提升了优化效率。该方法具有解析递推架构，易转化为C++程序，符合工程实践需求。因此，面向计算制导，本文方法能够同时满足实时性和精确性的需求。

2）针对现有方法未能充分考虑参数幅值约束的局限，本文通过融合凸规划技术提出了一种基于参数化控制向量的轨迹优化方法，进一步改善了方法的适应性。

1 飞行器任务场景与动力学模型

收起

本节首先建立了临近空间飞行器的任务场景，而后建立了以时间为自变量的动力学模型。针对以时间为自变量的轨迹优化问题存在终端飞行时间难以确定、可能无可行解的局限性，建立了以航程为自变量的动力学模型，降低了轨迹优化问题的复杂度，得到更加贴合工程实践的轨迹优化方法。

1.1 临近空间飞行器任务场景

未来高价值目标打击任务中，首先需要有人/无人战斗机群建立进攻型区域制空。在此基础上，编队后方的轰炸机群发射临近空间飞行器实现对目标的精确打击，其任务场景如图1所示。

1.2 以时间/航程为自变量的动力学模型

飞行器以时间为自变量的动力学模型为^[39]

$ \begin{cases}\dot{v}=g\left({n}_{x}-\sin \gamma \right)\\\dot{\gamma }=\dfrac{g}{v}\left({n}_{y}-\cos \gamma \right)\\\dot{\psi }=-\dfrac{g{n}_{{\textit{z}}}}{v\cos \gamma }\\{\dot{r}}_{x}=v\cos \gamma \cos \psi \\{\dot{r}}_{y}=v\sin \gamma \\{\dot{r}}_{{\textit{z}}}=-v\cos \gamma \sin \psi \end{cases} $

式中：6个状态量v、γ、

$\psi $

、r_x、r_y、r_z分别为速度、航迹角、航向角、航程、高度、侧向位置；g为重力加速度；状态向量x=[v，γ，

$\psi $

，r_x，r_y，r_z]^T；控制向量u=[n_y，n_z]^T，n_x、n_y和n_z分别为法向过载和侧向过载。关于动力学模型的参数设定，详见笔者前期的研究工作^[9]。

控制向量选取法向与侧向过载的原因在于：本文后续提出方法源于现有MPSP方法，通过符号递推解确保了快速准确解算。但是，此类方法已被证明属于一种改进的牛顿类方法。因此，该方法存在收敛半径小的问题，并且采用线性化近似模型时非线性程度不宜过高^[11]。为平衡解算效率与现有临近空间飞行器模型中非线性较高的问题，通过提出的轨迹优化方法获取虚拟控制量（法向过载与侧向过载），并反向递推得到实际控制量（攻角与倾侧角）。反向递推公式详见文献[38]，3.1节将验证该思路的可行性及效果。

结合式(1)，航程为自变量的动力学方程为

$ \begin{cases}\dfrac{\text{d}v}{\text{d}{r}_{x}}=\dfrac{g\left({n}_{x}-\sin \gamma \right)}{v\cos \gamma \cos \psi }\\\dfrac{\text{d}\gamma }{\text{d}{r}_{x}}=\dfrac{g\left({{n}}_{y}-\cos \gamma \right)}{{v}^{2}\cos \gamma \cos \psi }\\\dfrac{\text{d}\psi }{\text{d}{r}_{x}}=-\dfrac{g{n}_{{\textit{z}}}}{{v}^{2}{\cos }^{2}\gamma \cos \psi }\\\dfrac{\text{d}{r}_{y}}{\text{d}{r}_{x}}=\dfrac{\tan \gamma }{\cos \psi }\\\dfrac{\text{d}{r}_{{\textit{z}}}}{\text{d}{r}_{x}}=-\tan \psi \end{cases} $

鉴于中段飞行轨迹的高度不一定是单调的，使用航程作为自变量。考虑到滑翔飞行器实际的任务需求，中制导段航程可假设单调递增，因此航程可选取为自变量。

以时间为自变量的轨迹优化问题需要通过法向与侧向过载2个控制量实现6个终端状态的控制；与之不同的是，以航程为自变量的优化问题只需控制5个终端状态。因此，以航程为自变量建立动力学方程降低了原优化问题的维度，更有利于在线解算轨迹优化问题^[40-41]。

1.3 轨迹优化问题

首先，轨迹优化问题的性能指标将在2.2节中详细给出；然后，该问题需要满足初始状态约束和固定的终端状态约束；最后，系统动力学约束需严格遵循式(1)或式(2)。

2 基于控制向量参数化的轨迹优化

收起

基于实际动力学方程，本节提出了一种控制向量参数化的轨迹优化方法。首先，使用翻转Radau伪谱离散动力学方程；然后，通过控制向量参数化实现了轨迹优化问题的精确高效解算；最后，针对参数幅值约束，融合凸规划技术实现了轨迹优化问题的求解。值得注意的是，本节中的推导过程以时间为自变量。若以航程为自变量，只需将所有公式中的时间t换成航程r_x即可。

2.1 翻转Radau伪谱离散

动力学方程式(1)简化表达为

$ {\dot{\boldsymbol{x}}}=\boldsymbol{h}\left(\boldsymbol{x},\boldsymbol{u},t\right) $

式中：h为动力学方程向量；x和u分别为状态和控制向量；t为当前时间。

在参考轨迹附近对式(3)进行一阶泰勒展开^[38]：

$ {\dot{\boldsymbol{x}}}\approx \boldsymbol{H}={\boldsymbol{h}}^{\text{ref}}+\boldsymbol{h}_{{\boldsymbol{x}}}^{\text{ref}}\left(\boldsymbol{x}-{\boldsymbol{x}}^{\text{ref}}\right)+\boldsymbol{h}_{{\boldsymbol{u}}}^{\text{ref}}\left(\boldsymbol{u}-{\boldsymbol{u}}^{\text{ref}}\right) $

式中：H为近似的动力学方程；“ref”表示参考轨迹；

$ {\boldsymbol{h}}_{{\boldsymbol{x}}}^{\text{ref}} $

和

$ {\boldsymbol{h}}_{{\boldsymbol{u}}}^{\text{ref}} $

分别表示动力学方程向量对状态和控制向量的偏导数。

通过翻转Radau伪谱离散式(4)，该过程与笔者先前研究中的方法类似，理论详见文献[9,38]。式(5)～式(7)给出了该方法中的关键公式。

首先，原时间间隔t∈[t₀，t_f]通过归一化被转化至τ∈[−1，1]：

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

式中：t₀和t_f分别为给定的初始和终端时间。

然后，基于加权余量法则动力学约束表达为

$ {\boldsymbol{0}}=\sum \limits_{i=0}^{N-1}{\boldsymbol{D}}_{ki}\text{d}\boldsymbol{X}\left({\tau }_{i}\right)- \frac{{t}_{\text{f}}-{t}_{0}}{2}\boldsymbol{H}\left(\text{d}\boldsymbol{X}\left({\tau }_{k}\right),\text{d}\boldsymbol{U}\left({\tau }_{k}\right),{\tau }_{k};{t}_{0},{t}_{\text{f}}\right)$

式中：X和U分别为离散后的状态和控制向量；τ_i为节点（i=0,1,

$\cdots $

,K）；τ_k为配点（k=1,2,

$\cdots $

,K）；D_ki为微分矩阵中的元素，由参考轨迹计算得到。

最后，基于文献[38]中符号递推的思想，通过归纳构建了各状态间的映射关系：

$ \begin{split}& {\left[\begin{matrix}{\boldsymbol{X}}_{0} & {\boldsymbol{X}}_{1} & \cdots & {\boldsymbol{X}}_{k} & \cdots & {\boldsymbol{X}}_{K}\end{matrix}\right]}^{\text{T}}=\\&\qquad {\left[\begin{matrix}\boldsymbol{I} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}}\\{\boldsymbol{D}}_{10}\boldsymbol{I} & {\boldsymbol{D}}_{11}\boldsymbol{I} - \dfrac{{t}_{\text{f}}-{t}_{0}}{2}\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{1}}}}^{\text{ref}} & \cdots & {\boldsymbol{D}}_{1k}\boldsymbol{I} & \cdots & {\boldsymbol{D}}_{1K}\boldsymbol{I}\\\vdots & \vdots & & \vdots & & \vdots \\{\boldsymbol{D}}_{k0}\boldsymbol{I} & {\boldsymbol{D}}_{k1}\boldsymbol{I} & \cdots & {\boldsymbol{D}}_{kk}\boldsymbol{I} - \dfrac{{t}_{\text{f}}-{t}_{0}}{2}\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{k}}}}^{\text{ref}} & \cdots & {\boldsymbol{D}}_{kK}\boldsymbol{I}\\\vdots & \vdots & & \vdots & & \vdots \\{\boldsymbol{D}}_{K0}\boldsymbol{I} & {\boldsymbol{D}}_{K1}\boldsymbol{I} & \cdots & {\boldsymbol{D}}_{Kk}\boldsymbol{I} & \cdots & {\boldsymbol{D}}_{KK}\boldsymbol{I} - \dfrac{{t}_{\text{f}}-{t}_{0}}{2}\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{K}}}}^{\text{ref}}\end{matrix}\right]}^{-1} \left[\begin{matrix}{\boldsymbol{X}}_{0}\\\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\left[\boldsymbol{h}_{1}^{\text{ref}} - \boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{1}}}}^{\text{ref}}\boldsymbol{X}_{1}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{1}}}}^{\text{ref}}\left({\boldsymbol{U}}_{1}-\boldsymbol{U}_{1}^{\text{ref}}\right)\right]\\\vdots \\\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\left[\boldsymbol{h}_{k}^{\text{ref}} - \boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{k}}}}^{\text{ref}}\boldsymbol{X}_{k}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{k}}}}^{\text{ref}}\left({\boldsymbol{U}}_{k}-\boldsymbol{U}_{k}^{\text{ref}}\right)\right]\\\vdots \\\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\left[\boldsymbol{h}_{K}^{\text{ref}} - \boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{K}}}}^{\text{ref}}\boldsymbol{X}_{K}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{K}}}}^{\text{ref}}\left({\boldsymbol{U}}_{K}-\boldsymbol{U}_{K}^{\text{ref}}\right)\right]\end{matrix}\right]\end{split} $

2.2 基于符号递推的在线轨迹优化

本节提出了一种面向实际动力学方程具有符号递推解析架构的准确高效轨迹优化方法。定义式(7)动力学方程中的逆矩阵为P。据此，式(7)可进一步整理为

$ \begin{split}&\left[\begin{matrix}{\boldsymbol{X}}_{0}\\{\boldsymbol{X}}_{1}\\\vdots \\{\boldsymbol{X}}_{k}\\\vdots \\{\boldsymbol{X}}_{K}\end{matrix}\right]=\left[\begin{matrix}{\boldsymbol{P}}_{00} & {\boldsymbol{P}}_{01} & \cdots & {\boldsymbol{P}}_{0k} & \cdots & {\boldsymbol{P}}_{0K}\\{\boldsymbol{P}}_{10} & {\boldsymbol{P}}_{11} & \cdots & {\boldsymbol{P}}_{1k} & \cdots & {\boldsymbol{P}}_{1K}\\\vdots & \vdots & & \vdots & & \vdots \\{\boldsymbol{P}}_{k0} & {\boldsymbol{P}}_{k1} & \cdots & {\boldsymbol{P}}_{kk} & \cdots & {\boldsymbol{P}}_{kK}\\\vdots & \vdots & & \vdots & & \vdots \\{\boldsymbol{P}}_{K0} & {\boldsymbol{P}}_{K1} & \cdots & {\boldsymbol{P}}_{Kk} & \cdots & {\boldsymbol{P}}_{KK}\end{matrix}\right]\cdot \\&\qquad \left[\begin{matrix}{\boldsymbol{X}}_{0}\\\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\left(\boldsymbol{h}_{1}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{1}}}}^{\text{ref}}\boldsymbol{X}_{1}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{1}}}}^{\text{ref}}\left({\boldsymbol{U}}_{1}-\boldsymbol{U}_{1}^{\text{ref}}\right)\right)\\\vdots \\\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\left(\boldsymbol{h}_{k}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{k}}}}^{\text{ref}}\boldsymbol{X}_{k}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{k}}}}^{\text{ref}}\left({\boldsymbol{U}}_{k}-\boldsymbol{U}_{k}^{\text{ref}}\right)\right)\\\vdots \\\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\left(\boldsymbol{h}_{K}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{K}}}}^{\text{ref}}\boldsymbol{X}_{K}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{K}}}}^{\text{ref}}\left({\boldsymbol{U}}_{K}-\boldsymbol{U}_{K}^{\text{ref}}\right)\right)\end{matrix}\right]\end{split} $

式中：P_kk与式(7)中逆矩阵的对应元素保持一致。

基于式(8)，可进一步归纳递推关系：

$ \begin{split}&{\boldsymbol{X}}_{k}={\boldsymbol{P}}_{k0}{\boldsymbol{X}}_{0}+\frac{{t}_{\text{f}}-{t}_{0}}{2}\cdot \\&\qquad \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{ki}\left[\boldsymbol{h}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left({\boldsymbol{U}}_{i}-\boldsymbol{U}_{i}^{\text{ref}}\right)\right]\end{split} $

通过式(9)，终端时刻的状态转移方程为

$ \begin{split}&{\boldsymbol{X}}_{K}={\boldsymbol{P}}_{K0}{\boldsymbol{X}}_{0}+\frac{{t}_{\text{f}}-{t}_{0}}{2}\cdot \\&\qquad \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\left[\boldsymbol{h}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left({\boldsymbol{U}}_{i}-\boldsymbol{U}_{i}^{\text{ref}}\right)\right]\end{split} $

面向计算制导，为在线高精度解算轨迹优化问题，将各离散点处控制向量的优化问题转化为几个参数向量的优化问题，降低了维度与复杂性，提升了解算效率。参数化后的控制向量表达为

$ {\boldsymbol{U}}_{i}={\boldsymbol{p}}_{N+1}\tau _{i}^{N}+{\boldsymbol{p}}_{N}\tau _{i}^{N-1}+\cdots +{\boldsymbol{p}}_{n+1}\tau _{i}^{n}+\cdots +{\boldsymbol{p}}_{1} $

式中：p为参数矩阵（p=[p_N+1p_N

$\cdots $

p₁]），由N+1个待优化的参数列向量组成。

控制历史曲线设计为式(11)的原因为：多项式的形式不仅确保了曲线的连续性和光滑性，还为后续的内环姿态跟踪控制系统提供了良好的基础。此外，光滑的控制曲线决定了光滑的轨迹曲线，从而减小了阻力、增加了航程、降低了结构应力，进一步提升了飞行器的整体性能和可靠性。

为了简化推导，令N=2：

$ {\boldsymbol{U}}_{i}=\boldsymbol{a}\tau _{i}^{2}+\boldsymbol{b}{\tau }_{i}+\boldsymbol{c} $

式中：a、b、c分别对应离散轨迹优化问题中3个待优化的参数向量。因此，轨迹优化问题的优化变量由各离散点处的控制向量转化为3个参数向量。式(12)呈现出典型的抛物线特征，a决定了抛物线的开口大小及方向；a与b共同决定了抛物线对称轴的位置；c决定了初始0时刻控制量的大小。

为保证控制曲线的光滑性，二次型离散性能指标J设计为

$ J=\frac{{t}_{\text{f}}-{t}_{0}}{4}\left({\boldsymbol{a}}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{a}}}\boldsymbol{a}+{\boldsymbol{b}}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{b}}}\boldsymbol{b}+{\boldsymbol{c}}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{c}}}\boldsymbol{c}\right) $

式中：K_a、K_b与K_c分别为与3个参数向量对应的惩罚矩阵。

为证明通过式(13)进行迭代的收敛性，可将其改写为

$ \begin{split}& J=\frac{{t}_{\text{f}}-{t}_{0}}{4}\left[{\left({\boldsymbol{a}}_{\text{ref}}-\text{d}\boldsymbol{a}\right)}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{a}}}\left({\boldsymbol{a}}_{\text{ref}}-\text{d}\boldsymbol{a}\right)+{\left({\boldsymbol{b}}_{\text{ref}}-\text{d}\boldsymbol{b}\right)}^{\text{T}}\cdot \right.\\& \qquad \left. {\boldsymbol{K}}_{{\boldsymbol{b}}}\left({\boldsymbol{b}}_{\text{ref}}-\text{d}\boldsymbol{b}\right)+{\left({\boldsymbol{c}}_{\text{ref}}-\text{d}\boldsymbol{c}\right)}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{c}}}\left({\boldsymbol{c}}_{\text{ref}}-\text{d}\boldsymbol{c}\right)\right]\end{split} $

式中：“ref”表示参考值；设j表示a、b、c向量中的某一个分量，dj表示偏差值，值得注意的是，需要j与dj异号：

$ \frac{\partial J}{\partial j} \gt 0,\frac{\partial J}{\partial \text{d}j} \lt 0\rightarrow \begin{cases}J\searrow ,j\searrow ,\text{d}j\nearrow \\J\nearrow ,j\nearrow ,\text{d}j\searrow \end{cases} $

通过式(15)即可证明轨迹收敛^[9]。

结合式(10)中的终端状态转移方程与式(13)中的性能指标，增广性能指标

$ \overline{J} $

表达为

$ \begin{split}& \overline{J}=\frac{{t}_{\text{f}}-{t}_{0}}{4}\left({\boldsymbol{a}}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{a}}}\boldsymbol{a}+{\boldsymbol{b}}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{b}}}\boldsymbol{b}+{\boldsymbol{c}}^{\text{T}}{\boldsymbol{K}}_{{\boldsymbol{c}}}\boldsymbol{c}\right)+{\boldsymbol{\lambda }}^{\text{T}}\cdot \\&\quad \left\{{\boldsymbol{X}}_{K}-{\boldsymbol{P}}_{K0}{\boldsymbol{X}}_{0}-\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\cdot \right.\\&\quad\left. \displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\left[\boldsymbol{h}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left(\boldsymbol{a}\tau _{i}^{2}+\boldsymbol{b}{\tau }_{i}+\boldsymbol{c}-\boldsymbol{U}_{i}^{\text{ref}}\right)\right]\right\}\end{split} $

式中：

${\boldsymbol{\lambda }} $

为拉格朗日乘子向量。

不考虑参数幅值约束，3个参数向量的一阶最优必要条件为

$ \frac{\partial \overline{J}}{\partial \boldsymbol{a}}=\frac{{t}_{\text{f}}-{t}_{0}}{2}{\boldsymbol{K}}_{{\boldsymbol{a}}}\boldsymbol{a}-\frac{{t}_{\text{f}}-{t}_{0}}{2}{\left(\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2}\right)}^{\text{T}}\qquad \boldsymbol{\lambda }={\boldsymbol{0}} $

$ \frac{\partial \overline{J}}{\partial \boldsymbol{b}}=\frac{{t}_{\text{f}}-{t}_{0}}{2}{\boldsymbol{K}}_{{\boldsymbol{b}}}\boldsymbol{b}-\frac{{t}_{\text{f}}-{t}_{0}}{2}{\left(\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i}\right)}^{\text{T}}\qquad \boldsymbol{\lambda }={\boldsymbol{0}} $

$ \frac{\partial \overline{J}}{\partial \boldsymbol{c}}=\frac{{t}_{\text{f}}-{t}_{0}}{2}{\boldsymbol{K}}_{{\boldsymbol{c}}}\boldsymbol{c}-\frac{{t}_{\text{f}}-{t}_{0}}{2}{\left(\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\right)}^{\text{T}}\qquad \boldsymbol{\lambda }={\boldsymbol{0}} $

通过式(17)～式(19)，可以得到

$ \boldsymbol{a}=\boldsymbol{\lambda }\boldsymbol{K}_{{\boldsymbol{a}}}^{-1}{\left(\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2}\right)}^{\text{T}} $

$ \boldsymbol{b}=\boldsymbol{\lambda }\boldsymbol{K}_{{\boldsymbol{b}}}^{-1}{\left(\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i}\right)}^{\text{T}} $

$ \boldsymbol{c}=\boldsymbol{\lambda }\boldsymbol{K}_{{\boldsymbol{c}}}^{-1}{\left(\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\right)}^{\text{T}} $

结合式(12)与式(20)～式(22)：

借助式(23)，式(10)进一步转化为

$ \begin{split}&{\boldsymbol{X}}_{K}={\boldsymbol{P}}_{K0}{\boldsymbol{X}}_{0}+\frac{{t}_{\text{f}}-{t}_{0}}{2}\cdot \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\\&\qquad\left(\boldsymbol{h}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}}+\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left\{\left[\tau _{i}^{2}\boldsymbol{K}_{{\boldsymbol{a}}}^{-1}{\left(\displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2}\right)}^{\text{T}}+\right.\right.\right.\\&\qquad \left.\left.\left.{\tau }_{i}\boldsymbol{K}_{{\boldsymbol{b}}}^{-1}{\left(\displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i}\right)}^{\text{T}}+\right.\right.\right.\\&\qquad \left.\left.\left.\boldsymbol{K}_{{\boldsymbol{c}}}^{-1}{\left(\displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\right)}^{\text{T}}\right]\left(\boldsymbol{\lambda }-\boldsymbol{U}_{i}^{\text{ref}}\right)\right\}\right)\end{split} $

而后，通过式(24)解算符号形式的向量

$\boldsymbol{\lambda } $

：

$ \begin{split}& \boldsymbol{\lambda }=\left\{\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left[\tau _{i}^{2}\boldsymbol{K}_{{\boldsymbol{a}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2}\right)}^{\text{T}}+\right.\right.\\&\quad\left.\left.{\tau }_{i}\boldsymbol{K}_{{\boldsymbol{b}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i}\right)}^{\text{T}} + \boldsymbol{K}_{{\boldsymbol{c}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\right)}^{\text{T}}\right]\right\}^{-1}\cdot\\&\quad \left\{{\boldsymbol{X}}_{K} - {\boldsymbol{P}}_{K0}{\boldsymbol{X}}_{0} - \dfrac{{t}_{\text{f}} - {t}_{0}}{2}\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\left(\boldsymbol{h}_{i}^{\text{ref}} - \boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}} - \boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{U}_{i}^{\text{ref}}\right)\right\}\end{split} $

根据式(23)与式(25)，各离散点处的控制历史直接更新为

$ \begin{split}&{\boldsymbol{U}}_{i} = \left[\tau _{i}^{2}\boldsymbol{K}_{{\boldsymbol{a}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2}\right)}^{\text{T}}+{\tau }_{i}\boldsymbol{K}_{{\boldsymbol{b}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i}\right)}^{\text{T}}+\right.\\&\qquad \left.\boldsymbol{K}_{{\boldsymbol{c}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\right)}^{\text{T}}\right]\cdot \\&\qquad \left\{\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left[\tau _{i}^{2}\boldsymbol{K}_{{\boldsymbol{a}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2}\right)}^{\text{T}}+\right.\right.\\&\qquad \left.\left.{\tau }_{i}\boldsymbol{K}_{{\boldsymbol{b}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i}\right)}^{\text{T}} + \boldsymbol{K}_{{\boldsymbol{c}}}^{-1}{\left(\displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\right)}^{\text{T}}\right] \right\}^{-1} \cdot\\&\qquad \bigg\{{\boldsymbol{X}}_{K}-{\boldsymbol{P}}_{K0}{\boldsymbol{X}}_{0}-\dfrac{{t}_{\text{f}}-{t}_{0}}{2}\cdot \\&\qquad \displaystyle \sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\left[\boldsymbol{h}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{U}_{i}^{\text{ref}}\right]\bigg\}\end{split} $

根据上述推导，本文基于实际动力学方程，通过参数化控制向量，构建了一种面向计算制导的高效精确轨迹优化方法。值得注意的是，式(4)通过一阶泰勒展开得到线性形式的近似动力学方程与实际非线性动力学方程存在一定的误差。因此，本文方法借鉴了序列凸规划或MPSP中序列迭代的思想，实现了各终端状态的精确控制。

最后，上述轨迹优化问题的截止条件为达到最大迭代次数或满足绝对误差容限：

$ s=\text{Norm}\left(\left[\begin{matrix}{\boldsymbol{x}}_{0}-\boldsymbol{x}_{0}^{\text{ref}} & \cdots & {\boldsymbol{x}}_{i}-\boldsymbol{x}_{i}^{\text{ref}} & \cdots & {\boldsymbol{x}}_{K}-\boldsymbol{x}_{K}^{\text{ref}}\end{matrix}\right]\right) \lt {s}_{\text{abs}} $

式中：Norm()表示2-范数；s为相邻2次迭代轨迹间的距离；s_abs为绝对误差容限。

2.3 融合凸规划考虑参数幅值约束的轨迹优化

2.2节方法中3个一阶必要条件成立的前提是未考虑参数幅值约束。为有效处理该约束，给出了下列凸规划（二次规划）问题。

参数幅值约束表达为

$p_{\min }^{l}\leqslant p_{n}^{l}\leqslant p_{\max }^{l}$

式中：l表示第n个参数向量的第l个分量；

$p_{\min }^{l} $

和

$p_{\text{max}}^{l} $

分别为该分量的最小值和最大值。值得注意的是，参数幅值约束和控制输入约束密切相关。因此，在轨迹优化问题中引入参数幅值约束，能够更精准地契合实际的中制导场景。

通过式(10)与式(12)，建立动力学等式约束：

$\begin{split}&{\boldsymbol{X}}_{K}-{\boldsymbol{P}}_{K0}{\boldsymbol{X}}_{0}-\frac{{t}_{\text{f}}-{t}_{0}}{2}\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\left[\boldsymbol{h}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{x}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{X}_{i}^{\text{ref}}-\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\boldsymbol{U}_{i}^{\text{ref}}\right]=\\&\qquad \frac{{t}_{\text{f}}-{t}_{0}}{2}\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\left(\boldsymbol{a}\tau _{i}^{2}+\boldsymbol{b}{\tau }_{i}+\boldsymbol{c}\right)=\\&\qquad \frac{{t}_{\text{f}}-{t}_{0}}{2}\left[\begin{matrix}\displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\tau _{i}^{2} & \displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}{\tau }_{i} & \displaystyle\sum \limits_{i=1}^{K}{\boldsymbol{P}}_{Ki}\boldsymbol{h}_{{{\boldsymbol{u}}}_{{\boldsymbol{i}}}}^{\text{ref}}\end{matrix}\right] \left[\begin{matrix}\boldsymbol{a}\\\boldsymbol{b}\\\boldsymbol{c}\end{matrix}\right]\end{split} $

针对式(13)中的性能指标、式(28)中的参数幅值约束、式(29)中的等式约束及初始和终端状态约束，通过凸规划（二次规划）理论实现优化问题的精确高效解算。

需要注意的是，2.2节方法通过符号递推解直接更新控制历史，因此具备较高的计算效率，适用于在线解算计算制导下的轨迹优化问题。尽管该方法未能有效处理参数幅值约束（控制输入约束），但通过设计二次型性能指标，能够有效保证控制曲线的光滑性。在一般情况下，只要轨迹优化问题存在可行解，通常不会触碰输入约束。3.3节证明了该情况。

与之不同的是，2.3节方法因为考虑了额外的参数幅值约束，不具备符号递推解，仅能通过凸规划（二次规划）技术解算轨迹优化问题，牺牲了一定的解算效率^[42]。因此，在实际应用中，需结合战场态势动态评估并综合考虑，从2.2节或本节方法中选取合适的轨迹优化方法。

3 仿真验证与方法比较

收起

在一致的条件下，本节首先验证了基于实际动力学方程控制向量参数化的临近空间飞行器轨迹优化方法，并与现有通用伪谱最优控制软件（general pseudospectral optimal control software，GPOPS）轨迹优化方法、序列凸规划方法、基于偏差动力学方程的MPSC方法进行了对比。然后，针对大部分航空航天轨迹优化问题对终端飞行时间不敏感而对计算效率敏感的特点，将本文方法的自变量由时间改为航程，降低了优化问题的复杂度，并提升了优化效率。最后，面向参数幅值约束，对2.2节与2.3节中提出的2种方法进行了对比与分析。需要注意的是，2.2节方法未考虑参数幅值约束可构建符号递推解，因此不需要使用求解器。与之不同，2.3节方法因为处理幅值约束而不满足式(17)～式(19)中的一阶必要条件，因此可采用MATLAB环境下的CVXGEN求解器解算凸规划问题。

对于临近空间飞行器受限的机载计算资源，最大迭代次数设为20；绝对误差容限设为10；翻转Radau伪谱中离散点个数设为84；K_a、K_b、K_c这3个与参数向量对应的惩罚矩阵均设为单位阵；初始控制历史猜想（初值）采用比例制导生成；在仿真结果中，第7次迭代曲线表示最终优化得到的轨迹（见图2～图9）。初始和终端状态约束设定如表1所示。

在完成轨迹优化后，需要通过开环积分验证本文方法的终端控制精度。开环积分采用欧拉积分，步长设为0.001 s。

3.1 以时间为自变量控制向量参数化的轨迹优化

不考虑参数幅值约束，本节验证了以时间为自变量控制向量参数化的轨迹优化方法。首先，对2.2节方法进行了验证与分析。然后，将本文方法与现有MPSC方法^[34]、基于Radau伪谱的非线性轨迹优化方法（，GPOPS）^[43]及序列凸规划方法进行对比，说明了本文方法的有效性及意义。仿真中，终端飞行时间设为28.28 s，初始和终端状态约束如表1所示。图2～图4为轨迹曲线。

图2～图4中的第1次迭代曲线表示初始控制历史猜想积分得到的轨迹。可以明显看出，初始猜想与最终得到的最优解存在较大误差。表明即使初值较差，本文方法凭借较强的鲁棒性仍能保证终端状态的收敛。此外，终端速度误差为−2.90×10⁻³ m/s；航迹角与航向角误差分别为−4.90×10⁻⁵（°）与−2.60×10⁻⁴（°）；三轴位置误差分别为0.52 m、0.14 m与0.46 m。因此，以上数据充分说明了本文方法具有较高的终端状态控制精度。

图5～图7为航迹角、航向角与性能指标随时间的变化曲线。

图5表明，航迹角曲线在第3次迭代时收敛；图6和图7显示，航向角与性能指标曲线的收敛速度稍慢，直至第5次迭代才完成收敛。最终的性能指标为30.28。多次迭代的目的是在确保终端状态高精度控制的前提下降低性能指标。虽然该过程增加了计算时间，但是本文方法在MATLAB环境中单次优化的时间约为200 ms，若转化成C++程序，推测耗时可进一步降低至50 ms左右^[44]。即使迭代次数增加，轨迹优化的计算时间也不会过多增加，满足临近空间飞行器机载计算资源受限下的计算制导问题。

通过7次迭代，3个参数向量表示为

$ \begin{cases}\boldsymbol{a}=\left[\begin{matrix}0.287\;3\\-0.031\;5\end{matrix}\right]\\\boldsymbol{b}=\left[\begin{matrix}-0.134\;9\\-0.774\;0\end{matrix}\right]\\\boldsymbol{c}=\left[\begin{matrix}0.135\;0\\0.907\;5\end{matrix}\right]\end{cases} $

法向/侧向过载曲线如图8和图9所示。

通过式(5)与式(12)，可得控制历史表示为

$ \begin{split}& {n}_{y}\left(t\right)=0.287\;3 {\left(\frac{t}{14.14}-1\right)}^{2}-\\&\qquad 0.134\;9 \left(\frac{t}{14.14}-1\right)+0.135\;0\end{split} $

$ \begin{split}& {n}_{{\textit{z}}}\left(t\right)=-0.031\;5 {\left(\frac{t}{14.14}-1\right)}^{2}-\\& \qquad 0.774\;0 \left(\frac{t}{14.14}-1\right)+0.907\;5\end{split} $

通过式(31)与式(32)可以看出，控制曲线较为光滑，未出现斜率及幅值过大的情况，这一特性为后续的内环姿态跟踪控制系统提供了有利条件。

根据参考文献[38]，实际控制量攻角/倾侧角曲线反向解算结果如图10和图11所示。

从图10和图11可以看出，攻角曲线的幅值较小，不会导致失速的现象；鉴于倾侧角的范围为[−90°,90°]，虽然曲线存在细微的抖动，但仍能满足控制需求。攻角与倾侧角曲线较为光滑的原因在于式(31)和式(32)中选取的控制历史曲线本身是光滑的。因此，通过虚拟控制量（法向/侧向过载）反向解算实际控制量（攻角/倾侧角）的方案是可行的。

然后，基于MPSC、序列凸规划（sequential convex programming，SCP）与GPOPS这3种对比方法，通过仿真说明本文方法（DCHU-SC）能够在计算资源受限的约束下实现终端状态的精确控制。因为协态映射定理的存在，融合Radau伪谱与非线性规划的GPOPS轨迹优化方法生成的解满足一阶最优必要条件，可视为最优解。但是，该方法使用了大规模的非线性规划，一般离线计算，不适合在线解算轨迹优化问题。

MPSC方法的理论详见文献[28]。该方法中，3个待优化的参数向量更新方程为

$ {\left[\begin{matrix}\boldsymbol{a} & \boldsymbol{b} & \boldsymbol{c}\end{matrix}\right]}^{\text{T}}={\left[\begin{matrix}{\boldsymbol{C}}_{{\textit{z}}} & {\boldsymbol{D}}_{{\textit{z}}} & {\boldsymbol{E}}_{{\textit{z}}}\end{matrix}\right]}^{-1}\left[{\boldsymbol{B}}_{\lambda }-\text{d}{\boldsymbol{X}}_{K}\right] $

式中：C_z、D_z、E_z分别为与3个参数向量相关的矩阵；B_λ表示通过参考控制历史与敏感度矩阵计算得到的参数矩阵；dX_K表示终端状态向量与参考值的偏差，

$ \text{d}{\boldsymbol{X}}_{K}={\boldsymbol{X}}_{K}-\boldsymbol{X}_{K}^{\text{ref}} $

。

通过dX_K的公式可以看出，该方法基于偏差动力学方程进行推导。与之不同，本文方法基于实际动力学方程进行推导。图12～图14为4种不同方法的轨迹优化结果。其中，括号内的数字表示离散点个数，该参数影响SCP、DCHU-SC及MPSC这3种方法的计算效率与精度。

从图12～图14可以明显看出，相比于其他4条轨迹，采用84个离散点的MPSC方法存在较大误差。终端状态误差如表2所示。

从表2可以看出，仿真“MPSC (84)”的终端状态误差较大，位置误差超过了96 m，原因是终端速度未能精确控制。但是，仿真“MPSC (3 000)”不存在该问题，保证了较高的控制精度。产生该现象的原因为：当离散点个数从84个增加至3 000个，优化精度提升。但是，由于离散点个数直接决定轨迹优化所需的计算时间，随着个数的增加，解算时间也随之增加。在MATLAB环境下，仿真“MPSC (3 000)”的单轮轨迹优化时间约为5.5 s，推测在C++环境中约为1.8 s^[44]。即使在初始控制历史猜想较好的情况下，考虑到每次轨迹优化大约需要7轮迭代，其在C++环境下的优化时间为12～13 s。临近空间飞行器的平均速度约为1700 m/s，优化过程中大约飞行21 km，这意味着优化结果在实际飞行中无法直接使用。因此，现有基于偏差动力学方程的MPSC方法并不适用于机载计算资源受限下的飞行器在线轨迹优化任务。此外，与仿真“MPSC (84)”相比，仿真“SCP (84)”虽然降低了终端状态误差，但是位置误差仍大于16 m；同时，通过凸规划技术解算轨迹优化问题需进行非线性规划，牺牲了一定的解算效率^[42]。

与现有采用3 000个离散点的MPSC方法相比，本文方法（DCHU-SC）仅采用84个离散点，缩短了优化时间，同时进一步提升了控制精度。此外，本文方法生成的轨迹与GPOPS最优解基本重合，从仿真角度验证了轨迹的最优性，与2.2节中的理论推导高度契合。

上述分析表明，本文方法仅需优化3个参数向量即可获得最优解，降低了优化问题的复杂性，并提升了解算效率。因此，该方法在确保优化效率的同时，能够实现准确的终端状态控制，兼具高效性与精确性，为面向计算制导的轨迹优化领域提供了一种创新且实用的解决方案。

3.2 以航程为自变量控制向量参数化的轨迹优化

由于大多数航空航天轨迹优化问题对终端飞行时间不敏感，本节将验证以航程为自变量的轨迹优化效果。以航程为自变量而不是以时间为自变量的优点是：以时间为自变量的轨迹优化问题需通过法向与侧向过载2个控制量实现终端速度、航迹角、航向角、航程、高度与侧向位置6个状态的精确控制；与之不同，以航程为自变量的优化问题需通过2个控制量实现终端速度、航迹角、航向角、高度与侧向位置5个状态的控制。因此，相比于以时间为自变量，以航程为自变量的优化问题复杂度低、收敛半径大。即使初始控制历史猜想较差，轨迹也更容易收敛。

仿真初始设定如表1所示。终端速度误差为−0.02 m/s；航迹角与航向角误差分别为−6.50×10⁻⁴(°)与−5.45×10⁻⁴(°)；三轴位置误差分别为3.04 m、−0.29 m与3.45 m。从终端状态误差的统计结果来看，各状态均实现了高精度的控制。

然后，说明以航程为自变量方法的优越性。由于控制量直接作用于式(2)中的航迹/航向角微分方程，以上2个状态与控制量之间存在高度关联性。2个角度的变化曲线如图15和图16所示。控制曲线如图17和图18所示。

与图5、图6、图8与图9中以时间为自变量的航迹角、航向角、过载曲线相比，图15～图18表明以航程为自变量的轨迹优化方法具有更快的收敛速度。此外，本文方法生成的轨迹幅值较小、更为光滑，表明该方法不依赖于较大的收敛半径，进一步凸显了其性能的优越性。

图19表示性能指标的变化曲线。图7与图19的对比结果显示，误差容限相同时，以时间为自变量的曲线在第5次迭代时趋于收敛，最终迭代次数为7次；相比之下，以航程为自变量的曲线在第3次迭代时已实现收敛，且最终的迭代次数为5次。需要注意，尽管图中显示轨迹已收敛，但后续的几次迭代仍然是必要的。这些额外的迭代能够进一步提升终端状态的控制精度，从而确保最终结果的准确性和可靠性。

为验证本文方法对初始状态摄动的鲁棒性，开展了2 000次蒙特卡罗模拟测试。其中，在表1初始状态的基础上，初始高度摄动范围为[−300，300] m，侧向位置摄动范围为[−100，100] m，速度摄动范围为[−2，2] m/s，航迹/航向角摄动范围为[−1°，1°]，期望的终端状态均保持不变。需要注意的是，未施加初始航程摄动的原因是航程为自变量。蒙特卡罗仿真结果如图20～图24所示。

仿真结果表明，最大终端位置误差不超过75 m，最大航迹角与航向角误差不超过0.02°。此外，终端速度误差不超过0.3 m/s。因此，以上数据表明各终端状态实现了高精度收敛，满足中制导结束时的精度要求。因此，在对终端飞行时间敏感度相对较低的航空航天轨迹优化场景中，采用航程作为自变量取代传统的时间自变量，能够增强迭代过程的鲁棒性、提升收敛速度与优化效率。在面向计算制导的快速轨迹优化问题中，实时性是至关重要的核心指标。本文方法能够缩短计算时间，这一改进具有重要意义。

针对不确定建模下的鲁棒性进行蒙特卡罗模拟。部分论文通过非线性无迹变换的思想实现了鲁棒轨迹优化的效果，但导致了问题规模的增大，延长了计算时间^[45]。需要注意的是，轨迹优化方法普遍为开环方法，上述仿真已经证明本文方法具有快速性，因此可通过在线多次计算实现终端状态的精确控制^[38]。基于标准模型单次轨迹优化的输出，仅验证不确定模型开环积分的对比结果。首先，依据表1中的初始和终端状态约束构建参考轨迹。然后，考虑到模型不确定性可能导致的建模误差，对飞行器的阻力系数C_D、参考面积S及质量m分别施加±5%的摄动。在此基础上，利用上述随机摄动模型，通过2 000次蒙特卡罗模拟进行开环积分验证，以评估本文方法对不确定模型的适应性和鲁棒性，仿真结果如图25所示。

图25表明，尽管存在模型摄动，但本文提出的轨迹优化方法终端状态误差仍较小，表现出良好的适应性及鲁棒性。由于该方法生成的控制曲线是最优且光滑的，结合临近空间滑翔飞行器中制导段简单、稳定的气动条件，其受摄动的影响较小，因此仍能保持轨迹的稳定性。尽管在开环积分条件下，终端状态误差受到摄动的影响有所增加，但该方法具备快速计算的优势，能够通过多次迭代优化，实现高精度的终端状态控制^[38]。

3.3 融合凸规划考虑参数幅值约束的轨迹优化

本节将对比2.2节与2.3节中2种不同轨迹优化方法的结果。其中，2.2节提出未考虑参数幅值约束，但是式(13)中设计的二次型性能指标保证了控制曲线的光滑性。在轨迹优化问题存在可行解时，该方法通常不会违反参数幅值约束。与之不同，2.3节方法基于凸规划理论，成功实现了参数幅值约束的处理。

仿真以航程为自变量，初始和终端状态约束如表1所示。3个参数的幅值约束均设置为[−5，5]。其中，“SC”代表2.2节提出基于扩展控制的轨迹优化方法，其终端状态误差统计结果同3.2节。“CP”代表2.3节中基于凸规划的轨迹优化方法，其终端状态误差基本与仿真“SC”一致。

2种方法最终生成的控制历史基本接近：

$ {n}_{y}\left({r}_{x}\right) = 0.27 {\left(\frac{{r}_{x}}{17\;500} - 1\right)}^{2} - 0.15 \left(\frac{{r}_{x}}{17\;500} - 1\right)+0.15 $

$ {n}_{{\textit{z}}}\left({r}_{x}\right) = 0.24 {\left(\frac{{r}_{x}}{17\;500} - 1\right)}^{2} - 0.79 \left(\frac{{r}_{x}}{17\;500} - 1\right)+0.85 $

尽管在优化过程中存在参数幅值约束，但2.2节方法最终生成的控制曲线并未触发该约束。原因在于二次型性能指标确保了控制曲线的光滑性。因此，在优化问题存在可行解时，该方法通常不会触发参数幅值约束。

在构建航程-法向/侧向过载曲线后，进一步以时间为自变量，对式(1)中的实际动力学方程数值积分。法向/侧向过载曲线如图26和图27所示。与前述结论相似，2种方法生成的轨迹基本一致。

2.3节方法主要用于处理参数饱和的问题。然而，在上述仿真中，所有参数均未达到其幅值约束，因此未能充分展示该方法在处理饱和问题时的有效性。为了更全面地验证该方法的实际应用价值，后续仿真将针对参数饱和的情况展开研究。3个参数的幅值约束分别设置为：K_a∈[−3，3]，K_b∈[−0.7，0.7]，K_c∈[−2，2]。通过凸规划方法解算含有参数幅值约束的优化问题，控制历史表达为

$ {n}_{y}\left({r}_{x}\right) = 0.30 {\left(\frac{{r}_{x}}{17\;500} - 1\right)}^{2} - 0.15 \left(\frac{{r}_{x}}{17\;500} - 1\right)+0.14 $

$ {n}_{{\textit{z}}}\left({r}_{x}\right) = -1.13 {\left(\frac{{r}_{x}}{17\;500} - 1\right)}^{2} - 0.70 \left(\frac{{r}_{x}}{17\;500} - 1\right) + 1.31 $

从式(36)和式(37)可以看出，参数向量K_b触碰了参数幅值约束，证明了本文方法具备幅值约束的处理能力。数值仿真结果表明，三轴位置误差分别为−8.81 m、1.27 m与−9.86 m；速度误差为0.06 m/s；航迹角与航向角误差分别为1.96×10⁻³(°)与−1.41×10⁻²(°)。需要特别指出的是，本文方法因为采用符号递推解析架构，从而实现了较高的计算效率。相比之下，凸规划类方法需要将优化问题转化为工具箱要求的标准形式并通过非线性规划方法解算，这一过程在一定程度上牺牲了求解效率^[42]。因此，在实际应用中，需综合权衡约束条件与计算精度及效率之间的关系，以选择最合适的轨迹优化策略。

与上述仿真类似，以时间为自变量对动力学方程数值积分，法向/侧向过载曲线如图28和图29所示。

图28和图29的对比结果表明，侧向过载曲线在采用凸规划的改进方法后发生了显著变化。这说明该方法通过调整控制曲线，有效处理了参数幅值约束问题，进而实现了快速准确的轨迹优化。

4 结　论

收起

本文提出了一种处理初始和“硬”终端状态约束的在线轨迹优化方法，并得出以下结论：

1）基于实际动力学方程控制向量参数化的临近空间飞行器轨迹优化方法实现了终端状态的精确控制。终端速度误差为−2.90×10⁻³ m/s；航迹角与航向角误差分别为−4.90×10⁻⁵(°)与−2.60×10⁻⁴(°)；三轴位置误差分别为0.52 m、0.14 m与0.46 m，生成的轨迹与GPOPS所得的最优解高度一致，从仿真角度验证了其最优性，并与2.2节中的理论相契合。与现有MPSC方法相比，本文方法在离散点数量相同的条件下，能够实现终端状态的高精度控制。同时，该方法在C++环境下完成优化的时间仅为300～400 ms，满足计算制导下轨迹优化问题在线解算的实时性需求。

2）鉴于临近空间飞行器的中制导任务对终端飞行时间的敏感度较低，本文将以时间为自变量的轨迹优化问题转化为以航程为自变量的优化问题，降低了复杂度。该转换不仅减少了迭代次数，还提升了收敛速度。因此，本文方法适用于面向计算制导的轨迹优化任务。

3）针对参数幅值约束，提出了一种融合凸规划控制向量参数化的轨迹优化方法。由于额外约束的引入，导致问题的复杂度增加，从而产生了一定的终端状态误差。因此在实际应用中，需要平衡约束条件与计算精度之间的关系，以选择合适的轨迹优化策略。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

孙智孝, 杨晟琦, 朴海音, 等. 未来智能空战发展综述[J]. 航空学报, 2021, 42(8): 28-42.

Sun Z X, Yang S Q, Piao H Y, et al. A survey of air combat artificial intelligence[J]. Acta Aeronautica et Astronautica Sinica, 2021, 42(8): 28-42(in Chinese).

傅莉, 谢福怀, 孟光磊, 等. 基于滚动时域的无人机空战决策专家系统[J]. 北京航空航天大学学报, 2015, 41(11): 1994-1999.

Fu L, Xie F H, Meng G L, et al. An UAV air-combat decision expert system based on receding horizon control[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(11): 1994-1999(in Chinese).

Hu Y D, Liu Z C, Jiang C W, et al. Maneuver mode parametric modeling based on trajectory curve evolution laws for hypersonic glide vehicles[J]. Aerospace Science and Technology, 2025, 157: 109856.

闫循良, 王培臣, 郭杨. 再入滑翔机动突防轨迹规划与制导方法研究综述[J]. 航空学报, 2025, 46(17): 239-266.

Yan X L, Wang P C, Guo Y. Review of research on trajectory planning and guidance methods for entry glide maneuvering penetration[J]. Acta Aeronautica et Astronautica Sinica, 2025, 46(17): 239-266(in Chinese).

刘哲, 陆浩然, 郑伟, 等. 多滑翔飞行器时间协同轨迹快速规划[J]. 航空学报, 2021, 42(11): 310-324.

Liu Z, Lu H R, Zheng W, et al. Rapid time-coordination trajectory planning method for multi-glide vehicles[J]. Acta Aeronautica et Astronautica Sinica, 2021, 42(11): 310-324(in Chinese).

An K, Guo Z Y, Xu X P, et al. A framework of trajectory design and optimization for the hypersonic gliding vehicle[J]. Aerospace Science and Technology, 2020, 106: 106110.

Yu J L, Dong X W, Li Q D, et al. Cooperative guidance strategy for multiple hypersonic gliding vehicles system[J]. Chinese Journal of Aeronautics, 2020, 33(3): 990-1005.

孙景亮, 刘春生. 基于自适应动态规划的导弹制导律研究综述[J]. 自动化学报, 2017, 43(7): 1101-1113.

Sun J L, Liu C S. An overview on the adaptive dynamic programming based missile guidance law[J]. Acta Automatica Sinica, 2017, 43(7): 1101-1113(in Chinese).

Wang Y Z, Dai H H. Secure model predictive static programming with initial value generator for online computational guidance of near-space vehicles[J]. Aerospace Science and Technology, 2025, 156: 109768.

王雨琪, 宁国栋, 王晓峰, 等. 基于微分对策的临近空间飞行器机动突防策略[J]. 航空学报, 2020, 41(增刊2): 69-78.

Wang Y Q, Ning G D, Wang X F, et al. Maneuver penetration strategy of near space vehicle based on differential game[J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(Sup 2): 69-78(in Chinese).

Pan B F, Ma Y Y, Yan R. Newton-type methods in computational guidance[J]. Journal of Guidance, Control, and Dynamics, 2019, 42(2): 377-383.

卢青, 周军, 呼卫军. 基于改进A*算法的滑翔飞行器轨迹规划[J]. 系统工程与电子技术, 2016, 38(12): 2758-2763.

Lu Q, Zhou J, Hu W J. Trajectory planning for gliding aircraft based on improved A* algorithm[J]. Systems Engineering and Electronics, 2016, 38(12): 2758-2763(in Chinese).

Wang Z, Bao W M, Li H F. Second-order dynamic sliding-mode control for nonminimum phase underactuated hypersonic vehicles[J]. IEEE Transactions on Industrial Electronics, 2017, 64(4): 3105-3112.

Lu P. Introducing computational guidance and control[J]. Journal of Guidance, Control, and Dynamics, 2017, 40(2): 193.

Malyuta D, Yu Y, Elango P, et al. Advances in trajectory optimization for space vehicle control[J]. Annual Reviews in Control, 2021, 52: 282-315.

Shirazi A, Ceberio J, Lozano J A. Spacecraft trajectory optimization: a review of models, objectives, approaches and solutions[J]. Progress in Aerospace Sciences, 2018, 102: 76-98.

Kazim M, Hong J, Kim M G, et al. Recent advances in path integral control for trajectory optimization: an overview in theoretical and algorithmic perspectives[J]. Annual Reviews in Control, 2024, 57: 100931.

Wu D, Cheng L, Jiang F H, et al. Analytical costate estimation by a reference trajectory-based least-squares method[J]. Journal of Guidance, Control, and Dynamics, 2022, 45(8): 1529-1537.

刘宇航, 杨洪伟, 李爽. 小推力最优轨迹协态估计的高效机器学习方法[J]. 宇航学报, 2022, 43(5): 593-602.

Liu Y H, Yang H W, Li S. Efficient machine learning method for co-state estimation of low-thrust optimal trajectories[J]. Journal of Astronautics, 2022, 43(5): 593-602(in Chinese).

Nurre N P, Taheri E. Comparison of indirect and direct methods for operationally constrained low-thrust trajectories[J]. Journal of Guidance, Control, and Dynamics, 2025, 48(4): 885-902.

武迪, 郭祥, 蒋方华, 等. 连续小推力轨迹优化间接法研究进展综述[J]. 宇航学报, 2024, 45(7): 995-1007.

Wu D, Guo X, Jiang F H, et al. Progress of indirect method for continuous low-thrust trajectory optimization[J]. Journal of Astronautics, 2024, 45(7): 995-1007(in Chinese).

任鹏飞, 王洪波, 周国峰. 基于自适应伪谱法的高超声速飞行器再入轨迹优化[J]. 北京航空航天大学学报, 2019, 45(11): 2257-2265.

Ren P F, Wang H B, Zhou G F. Reentry trajectory optimization for hypersonic vehicle based on adaptive pseudospectral method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(11): 2257-2265(in Chinese).

雍恩米, 魏桐, 刘深深, 等. 一种基于伪谱法的助推-补能跳跃滑翔轨迹优化方法[J]. 宇航学报, 2024, 45(6): 893-902.

Yong E M, Wei T, Liu S S, et al. A trajectory optimization method for boost-energy supplementing skip gliding motion base on pseudospectral method[J]. Journal of Astronautics, 2024, 45(6): 893-902(in Chinese).

赵吉松, 李畅, 王晨曦, 等. 吸气式高超声速飞行器助推-巡航轨迹优化[J]. 宇航学报, 2021, 42(7): 820-829.

Zhao J S, Li C, Wang C X, et al. Boost-cruise trajectory optimization of air breathing hypersonic vehicles[J]. Journal of Astronautics, 2021, 42(7): 820-829(in Chinese).

Lu P. What is guidance?[J]. Journal of Guidance, Control, and Dynamics, 2021, 44(7): 1237-1238.

Padhi R, Banerjee A, Mathavaraj S, et al. Computational guidance using model predictive static programming for challenging space missions: an introductory tutorial with example scenarios[J]. IEEE Control Systems Magazine, 2024, 44(2): 55-80.

Patel J, Subbarao K. Reachable-set constrained model predictive control for safe hypersonic reentry[J]. Journal of Guidance, Control, and Dynamics, 2025, 48(6): 1253-1265.

Engel D L, Putnam Z R, Dutta S. Assessment of model predictive control for Mars entry vehicles with flaps[J]. Journal of Spacecraft and Rockets, 2025, 62(5): 1528-1545.

Malyuta D, Reynolds T P, Szmuk M, et al. Convex optimization for trajectory generation: a tutorial on generating dynamically feasible trajectories reliably and efficiently[J]. IEEE Control Systems Magazine, 2022, 42(5): 40-113.

Hayner C R, Carson III J M, Açıkmeşe B, et al. Continuous-time line-of-sight constrained trajectory planning for 6-degree of freedom systems[J]. IEEE Robotics and Automation Letters, 2025, 10(5): 4332-4339.

Ma Y Y, Pan B F, Tan Q D. Neural convex optimization for real-time trajectory generation of asteroid landings[J]. Acta Astronautica, 2025, 229: 606-615.

Benedikter B, Zavoli A, Colasurdo G, et al. Convex approach to three-dimensional launch vehicle ascent trajectory optimization[J]. Journal of Guidance, Control, and Dynamics, 2021, 44(6): 1116-1131.

Wang Z B, Grant M J. Autonomous entry guidance for hypersonic vehicles by convex optimization[J]. Journal of Spacecraft and Rockets, 2018, 55(4): 993-1006.

Dwivedi P N, Bhattacharya A, Padhi R. Suboptimal midcourse guidance of interceptors for high-speed targets with alignment angle constraint[J]. Journal of Guidance, Control, and Dynamics, 2011, 34(3): 860-877.

张永华, 王铮, 王时雨, 等. 基于无迹模型预测静态规划的可回收火箭着陆段鲁棒制导方法[J/OL]. 北京航空航天大学学报, 2024 (2024-09-14)[2025-04-18]. 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[J/OL]. Journal of Beijing University of Aeronautics and Astronautics, 2024 (2024-09-14) [2025-04-18]. https://doi.org/10.13700/j.bh.1001-5965.2024.0491(in Chinese).

谭一廷, 荆武兴, 高长生. 考虑零控交班视窗角约束的拦截中制导设计[J]. 宇航学报, 2021, 42(10): 1257-1270.

Tan Y T, Jing W X, Gao C S. Design of interception midcourse guidance with zero effort handover visual angle constraint[J]. Journal of Astronautics, 2021, 42(10): 1257-1270(in Chinese).

Ma S D, Yang Y X, Yang H, et al. Trajectory optimization of hypersonic vehicle considering the quasi-static assumption of pitch motion[J]. Aerospace Science and Technology, 2024, 146: 108969.

王远卓, 代洪华, 岳晓奎. 基于实际动力学方程的飞行器准谱轨迹优化[J]. 航空学报, 2025, 46(22): 184-203.

Wang Y Z, Dai H H, Yue X K. Quasi-spectral trajectory optimization for vehicle based on actual dynamic equations[J]. Acta Aeronautica et Astronautica Sinica, 2025, 46(22): 184-203(in Chinese).

Lu P, Doman D B, Schierman J D. Adaptive terminal guidance for hypervelocity impact in specified direction[J]. Journal of Guidance, Control, and Dynamics, 2006, 29(2): 269-278.

王远卓. 高超声速飞行器邻近最优制导律研究[D]. 北京: 北京航空航天大学, 2024: 13-19.

Wang Y Z. Research on neighboring optimal guidance law for hypersonic vehicle[D]. Beijing: Beihang University, 2024: 13-19(in Chinese).

王远卓, 张冉, 李惠峰. 基于协态估计的火箭动力下降邻近最优制导[J]. 宇航学报, 2024, 45(5): 741-752.

Wang Y Z, Zhang R, Li H F. Neighboring optimal guidance for rocket powered descent based on costate estimation[J]. Journal of Astronautics, 2024, 45(5): 741-752(in Chinese).

Liu X, Li S, Xin M. Pseudospectral convex optimization based model predictive static programming for constrained guidance[J]. IEEE Transactions on Aerospace and Electronic Systems, 2023, 59(3): 2232-2244.

Patterson M A, Rao A V. GPOPS-II: a MATLAB software for solving multiple-phase optimal control problems using hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming[J]. ACM Transactions on Mathematical Software, 2014, 41(1): 1-37.

Agamawi Y M, Rao A V. CGPOPS: a C++ software for solving multiple-phase optimal control problems using adaptive Gaussian quadrature collocation and sparse nonlinear programming[J]. ACM Transactions on Mathematical Software, 2020, 46(3): 1-38.

Mathavaraj S, Padhi R. Quasi-spectral unscented MPSP guidance for robust soft-landing on asteroid[J]. Journal of Optimization Theory and Applications, 2021, 191(2): 823-845.

2026年第52卷第10期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.13700/j.bh.1001-5965.2025.0247

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

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-04-29
录用日期：2025-07-25

基金

作者信息

西北工业大学航天学院，西安 710072

通讯作者:

E-mail：hhdai@nwpu.edu.cn

参考文献

分享链接

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

分享至

全文二维码

扫描看全文

引用本文

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