火星大气进入段纵向可达区生成的解析同伦法

摘要

本发明公开的火星大气进入段纵向可达区生成的解析同伦法，属于深空探测领域。本发明实现方法包括如下步骤：建立火星大气进入段探测器纵向动力学模型；基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓出最大航程；基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓得到最小航程；基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓得到最大开伞高度；基于解析同伦法，根据延拓得到的最大航程、最小航程和最大开伞高度，通过构造同伦参数，延拓出纵向可达区。本发明能够避免协态初值猜测和内部点可达性验证，进而提高可达区生成效率。

说明书

技术领域

本发明涉及一种火星大气进入段纵向可达区生成的解析同伦法，属于深空探测领域。

背景技术

火星大气进入段指的是探测器从进入火星大气直至降落伞展开的飞行阶段，火星大气进入段可达区分析对于任务设计、着陆点选取及风险评估等均具有重要意义。可达区求解涉及到进入段轨迹优化问题，可利用直接法或间接法求解。对于直接法，已有学者通过网格划分，结合数值优化方法求解大气进入段可达区。该方法需要对每一个网格点单独求取最优化问题，计算量大，求解依赖于初值猜测，求解过程的收敛性及可靠性难以保证。间接法的理论基础是庞特里亚金极小值原理，通过将原优化问题转换为一个等价的两点边值问题求解，由于极小值理论中的协状态缺乏物理意义，且优化问题对协状态初值比较敏感，故间接法存在协态初值猜测困难的问题。特别是对于可以横纵向解耦的小升阻比火星进入探测器，基于其纵向动力学方程的最优控制问题一般均属于bang-bang控制，初值收敛域小且对初值敏感。目前，同伦法在协状态初值猜测方面已有一定的理论基础，在求解小推力星际转移轨道最优问题方面已经有了较为深入的研究。在火星大气进入段轨迹优化方面，已有学者利用同伦法求解了最大开伞高度问题。

发明内容

针对火星大气进入段飞行器的飞行能力分析问题，本发明公开的火星大气进入段纵向可达区生成的解析同伦法要解决的技术问题是：基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓出纵向可达区(开伞点高度-航程剖面)，所述同伦过程能够避免协态初值猜测和内部点可达性验证，进而提高可达区生成效率。

本发明目的是通过下述技术方案实现的。

本发明公开的火星大气进入段纵向可达区生成的解析同伦法，

本发明公开的火星大气进入段纵向可达区生成的解析同伦法，包括如下步骤：

步骤一、建立火星大气进入段探测器纵向动力学模型；

在火星惯性坐标系下，忽略火星自转，取探测器的纵向平面内运动状态为x＝[r,V,γ,s]^T，其中，r探测器质心到火星质心的距离，V为探测器速度大小，γ为飞行路径角，s为航程，则大气进入段无量纲的纵向动力学模型为：

式(1)中，τ为无量纲时间，σ为倾侧角。在无量纲化过程中，长度的量纲单位为火星半径R₀，速度的无量纲单位为其中为火表引力加速度，μ为火星引力常数。时间的无量纲单位为角度的单位为弧度，不需要无量纲化处理。式(1)中，L和D分别为探测器受到的无量纲升力和阻力加速度，分别具有如下形式：

L＝D·L/D(3)

式(1)中，B为探测器的弹道系数，L/D为探测器的升阻比，ρ为火星大气密度，采用如下指数模型：

式(4)中，ρ₀为参考密度，r₀为参考半径，h为探测器的飞行高度，h_s为标高。

步骤二、基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓出最大航程；

取最大航程的辅助优化问题如式(5)、式(6)所示。

问题寻找最优控制使得

满足约束：

其中，σ_max和σ_min分别为倾侧角可以取的最大值及最小值，x₀为探测器在大气进入点的初始状态，V_f为探测器顺利开伞所需满足的末端速度约束。

根据庞特里亚金极小值原理，引入协状态λ＝[λ_r>V>γ>s]^T，问题对应的哈密顿函数为

问题的协状态微分方程为

问题的横截条件如下

由于中不显含时间项，故

由于0≤u_max-u≤u_max-u_min，且终端条件只有速度约束，故将动力学方程按照常值倾侧角σ_min正向积分直到满足终端速度约束为止，即得到问题对应的最优状态和飞行时间。

由式(10)知，将式(9)代入式(10)，得终端速度协状态值

λ_V(τ_f)＝0(11)

根据式(11)和式(9)的协状态终端值、式(8)的协状态微分方程以及正向积分得到的最优状态和飞行时间逆向积分即可得到辅助优化问题对应的最优轨迹协状态值，即辅助优化问题的最优解是解析已知的。

由辅助优化问题延拓得到原最优问题解的关键在于构建同伦参数。将同伦参数ε置于性能指标中，构建如下的优化问题

问题寻找最优控制使得

满足约束：

问题的哈密顿函数为

由于中不显含时间，故成立，则

控制量u以一阶形式出现，故问题的最优解是bang-bang形式，最优解具有如下的形式

问题的协状态微分方程与问题P_s,max相同。横截条件如下

问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z，满足式(13)、式(11)、式(8)和(17)描述的边值条件和微分方程条件，式(15)中的哈密顿函数条件及式(16)中的控制切换条件。

由式(14)知，当ε＝0时，问题等价于问题当ε＝1时，问题等价于最大航程问题。构造递增序列{ε₁,ε₂,...ε_k,ε_k+1,...ε_N}，其中ε₁＝0，ε_N＝1。从ε₁开始，以第k步的解z_k作为第k+1步的问题的初值，依次求解两点边值问题，即能够最终得到原最大航程问题的最优解，即延拓得到最大航程。

步骤三、基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓得到最小航程。

构造辅助约束优化问题如式(18)和式(19)所示。

问题寻找最优控制使得

满足约束：

问题的哈密顿函数如下

由于中不显含时间项，故

由于0≤u-u_min≤u_max-u_min，且终端条件只有速度约束，故将动力学方程按照常值倾侧角σ_max正向积分直到满足终端速度约束为止，即得到问题对应的最优状态和飞行时间，即辅助优化问题的最优解是解析已知的。

问题为拉格朗日问题，性能指标中只包含与状态变量无关的积分项，且中不显含时间项。故问题的横截条件与问题的横截条件相同，即

从式(22)所示的末端条件出发，根据式(8)的微分方程及积分得到的最优状态和飞行时间，逆向积分即得到该辅助优化问题所对应的协状态值。

构建如公式(24)所示的优化问题

问题寻找最优控制使得

满足约束：

问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z，满足式(8)、式(11)、式(22)、和式(23)描述的边值条件和微分方程条件，式(15)中的哈密顿函数条件及式(16)中的控制切换条件。

由式(23)知，当ε＝0时，问题等价于问题当ε＝1时，问题等价于最小航程问题。构造递增序列{ε₁,ε₂,...ε_k,ε_k+1,...ε_N}，其中ε₁＝0，ε_N＝1。从ε₁开始，以第k步的解z_k作为第k+1步的问题的初值，依次求解两点边值问题，即能够最终得到原最小航程问题的最优解，即延拓得到最小航程。

步骤四、基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓得到最大开伞高度；

取性能指标为构建辅助优化问题如公式(25)、(26)所示。

问题寻找最优控制使

满足约束：

由于0≤u²，故问题的最优解为常值σ＝π/2。获取最优状态的方法为从初始状态开始，以常值倾侧角正向积分飞行轨迹。直到满足末端速度天剑为止，得到辅助最优问题的最优状态和飞行时间，即辅助优化问题的最优解是解析已知的。

问题的协态微分方程与问题相同，横截条件与问题相同。根据协状态微分方程逆向积分即得到最优轨迹对应的协状态值。

为从辅助优化问题得到最大开伞高度问题的最优解，取性能指标为构建如公式(27)、(28)所示的优化问题。

问题寻找最优控制使得

满足约束：

问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z，满足式(8)、式(11)、式(22)、和式(28)描述的边值条件和微分方程条件，式(15)中的哈密顿函数条件及式(16)中的控制切换条件。

由式(27)知，当ε＝0时，问题等价于问题当ε＝1时，问题等价于最大开伞高度问题。构造递增序列{ε₁,ε₂,...ε_k,ε_k+1,...ε_N}，其中ε₁＝0，ε_N＝1。从ε₁开始，以第k步的解z_k作为第k+1步的问题的初值，依次求解两点边值问题，即能够最终得到原最大开伞高度问题的最优解，即延拓得到最小航程。

步骤五、基于解析同伦法，根据步骤二延拓得到的最大航程、步骤三延拓得到的最小航程和步骤四延拓得到的最大开伞高度，通过构造同伦参数，延拓出纵向可达区，所述纵向可达区为开伞点高度-航程剖面。

取性能指标为J_sk＝-r(τ_f)，构建如公式(29)、(30)所示的优化模型。

问题RA：寻找使得

满足约束：

问题RA与问题P_h,max哈密顿函数相同，故问题RA的横截条件为

构造两个单调航程序列：{s₀,...s_i,s_i+1,...s_I}与{s₀,...s_j,s_j+1,...s_J}，满足s₀为问题P_h,max对应的航程，s_I为问题P_s,min对应的最小航程，s_J为问题P_s,max对应的最大航程，分别从问题P_h,max的最优解出发，按照序列依次求解问题RA，其中第k次的最优解z_k作为第k+1次问题求解初值，延拓出纵向可达区，所述纵向可达区为开伞点高度-航程剖面。

有益效果：

1、本发明公开的火星大气进入段纵向可达区生成的解析同伦法，在构造最优解已知的辅助优化问题基础上，通过构造合适的同伦参数，能够快速求解最大/小航程问题及最大开伞高度问题，有效避免最优问题协态初值猜测过程，同伦求解过程稳定。

2、本发明公开的火星大气进入段纵向可达区生成的解析同伦法，根据步骤二延拓得到的最大航程、步骤三延拓得到的最小航程和步骤四延拓得到的最大开伞高度，通过构造合适的同伦参数，延拓出纵向可达区，能够避免对可达区内部点可达性的验证，进而提高可达区生成效率。

附图说明

图1为火星大气进入段纵向可达区生成的解析同伦法流程图。

图2为问题中初始协状态值随参数ε的变化。

图3为问题中航程随参数ε的变化关系。

图4为问题中初始协状态值随参数ε的变化。

图5为问题中航程随参数ε的变化关系。

图6为问题中初始协状态值随参数ε的变化。

图7为问题中航程随参数ε的变化关系。

图8为协状态初值随航程的变化关系。

图9为不同初始航迹角下的纵向可达区。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实施实例对发明内容做进一步说明。

本实例为针对火星大气进入段纵向可达区求解问题，利用解析同伦法，在求解最大/小航程和最大开伞高度问题的基础上，进一步利用延拓的方法得到纵向可达区。

火星大气进入段纵向可达区生成的解析同伦法，如图1所示，具体步骤如下：

步骤一、建立火星大气进入段探测器纵向动力学模型；

在火星惯性坐标系下，忽略火星自转，取探测器的纵向平面内运动状态为x＝[r,V,γ,s]^T，其中，r探测器质心到火星质心的距离，V为探测器速度大小，γ为飞行路径角，s为航程，则大气进入段无量纲的纵向动力学模型为：

式(1)中，τ为无量纲时间，σ为倾侧角。在无量纲化过程中，长度的量纲单位为火星半径R₀，速度的无量纲单位为其中为火表引力加速度，μ为火星引力常数。时间的无量纲单位为角度的单位为弧度，不需要无量纲化处理。式(1)中，L和D分别为探测器受到的无量纲升力和阻力加速度，分别具有如下形式：

L＝D·L/D(3)

式(1)中，B为探测器的弹道系数，L/D为探测器的升阻比，ρ为火星大气密度，采用如下指数模型：

式(4)中，ρ₀为参考密度，r₀为参考半径，h为探测器的飞行高度，h_s为标高。

步骤二、基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓出最大航程；

取最大航程的辅助优化问题如式(5)、式(6)所示。

问题寻找最优控制使得

满足约束：其中，σ_max和σ_min分别为倾侧角可以取的最大值及最小值，x₀为探测器在大气进入点的初始状态，V_f为探测器顺利开伞所需满足的末端速度约束。

根据庞特里亚金极小值原理，引入协状态λ＝[λ_r>V>γ>s]^T，问题对应的哈密顿函数为

问题的协状态微分方程为

问题的横截条件如下

由于中不显含时间项，故

由于0≤u_max-u≤u_max-u_min，且终端条件只有速度约束，故将动力学方程按照常值倾侧角σ_min正向积分直到满足终端速度约束为止，即得到问题对应的最优状态和飞行时间。

由式(10)知，将式(9)代入式(10)，得终端速度协状态值

λ_V(τ_f)＝0(11)

根据式(11)和式(9)的协状态终端值、式(8)的协状态微分方程以及正向积分得到的最优状态和飞行时间逆向积分即可得到辅助优化问题对应的最优轨迹协状态值，即辅助优化问题的最优解是解析已知的。

由辅助优化问题延拓得到原最优问题解的关键在于构建同伦参数。将同伦参数ε置于性能指标中，构建如下的优化问题

问题寻找最优控制使得

满足约束：

问题的哈密顿函数为

由于中不显含时间，故成立，则

控制量u以一阶形式出现，故问题的最优解是bang-bang形式，最优解具有如下的形式

问题的协状态微分方程与问题P_s,max相同。横截条件如下

问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z，满足式(13)、式(11)、式(8)和(17)描述的边值条件和微分方程条件，式(15)中的哈密顿函数条件及式(16)中的控制切换条件。

由式(14)知，当ε＝0时，问题等价于问题当ε＝1时，问题等价于最大航程问题。构造递增序列{ε₁,ε₂,...ε_k,ε_k+1,...ε_N}，其中ε₁＝0，ε_N＝1。从ε₁开始，以第k步的解z_k作为第k+1步的问题的初值，依次求解两点边值问题，即能够最终得到原最大航程问题的最优解，即延拓得到最大航程。

步骤三、基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓得到最小航程。

构造辅助约束优化问题如式(18)和式(19)所示。

问题寻找最优控制使得

满足约束：

问题的哈密顿函数如下

由于中不显含时间项，故

由于0≤u-u_min≤u_max-u_min，且终端条件只有速度约束，故将动力学方程按照常值倾侧角σ_max正向积分直到满足终端速度约束为止，即得到问题对应的最优状态和飞行时间，即辅助优化问题的最优解是解析已知的。

问题为拉格朗日问题，性能指标中只包含与状态变量无关的积分项，且中不显含时间项。故问题的横截条件与问题的横截条件相同，即

从式(22)所示的末端条件出发，根据式(8)的微分方程及积分得到的最优状态和飞行时间，逆向积分即得到该辅助优化问题所对应的协状态值。

构建如公式(24)所示的优化问题

问题寻找最优控制使得

满足约束：

问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z，满足式(8)、式(11)、式(22)、和式(23)描述的边值条件和微分方程条件，式(15)中的哈密顿函数条件及式(16)中的控制切换条件。

由式(23)知，当ε＝0时，问题等价于问题当ε＝1时，问题等价于最小航程问题。构造递增序列{ε₁,ε₂,...ε_k,ε_k+1,...ε_N}，其中ε₁＝0，ε_N＝1。从ε₁开始，以第k步的解z_k作为第k+1步的问题的初值，依次求解两点边值问题，即能够最终得到原最小航程问题的最优解，即延拓得到最小航程。

步骤四、基于解析同伦法，从最优解解析已知的辅助问题出发，通过构造同伦参数，延拓得到最大开伞高度；

取性能指标为构建辅助优化问题如公式(25)、(26)所示。

问题寻找最优控制使

满足约束：

由于0≤u²，故问题的最优解为常值σ＝π/2。获取最优状态的方法为从初始状态开始，以常值倾侧角正向积分飞行轨迹。直到满足末端速度天剑为止，得到辅助最优问题的最优状态和飞行时间，即辅助优化问题的最优解是解析已知的。

问题的协态微分方程与问题相同，横截条件与问题相同。根据协状态微分方程逆向积分即得到最优轨迹对应的协状态值。

为从辅助优化问题得到最大开伞高度问题的最优解，取性能指标为构建如公式(27)、(28)所示的优化问题。

问题寻找最优控制使得

满足约束：

问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z，满足式(8)、式(11)、式(22)、和式(28)描述的边值条件和微分方程条件，式(15)中的哈密顿函数条件及式(16)中的控制切换条件。

由式(27)知，当ε＝0时，问题等价于问题当ε＝1时，问题等价于最大开伞高度问题。构造递增序列{ε₁,ε₂,...ε_k,ε_k+1,...ε_N}，其中ε₁＝0，ε_N＝1。从ε₁开始，以第k步的解z_k作为第k+1步的问题的初值，依次求解两点边值问题，即能够最终得到原最大开伞高度问题的最优解，即延拓得到最小航程。

步骤五、基于解析同伦法，根据步骤二延拓得到的最大航程、步骤三延拓得到得最小航程和步骤四延拓得到的最大开伞高度，通过构造同伦参数，延拓出纵向可达区，所述纵向可达区为开伞点高度-航程剖面。

取性能指标为构建如公式(29)、(30)所示的优化模型。

问题RA：寻找使得

满足约束：

问题RA与问题P_h,max哈密顿函数相同，故问题RA的横截条件为

构造两个单调航程序列：{s₀,...s_i,s_i+1,...s_I}与{s₀,...s_j,s_j+1,...s_J}，满足s₀为问题P_h,max对应的航程，s_I为问题P_s,min对应的最小航程，s_J为问题P_s,max对应的最大航程，分别从问题P_h,max的最优解出发，按照序列依次求解问题RA，其中第k次的最优解z_k作为第k+1次问题求解初值，延拓出纵向可达区，所述纵向可达区为开伞点高度-航程剖面。为了保证探测器的安全，约束可达区的最低高度为6km。

步骤六、在MATLAB环境下对上述算法进行仿真分析；

探测器的物理参数如下：质量m为2800kg，参考面积S_ref为15.9m²，升阻比L/D为0.24，阻力系数C_D为1.45，初始条件及开伞条件具体如表1所示。为了给制导控制系统留有控制余量，取最大倾侧角为σ_max为120°，最小倾侧角σ_min为30°。在同伦延拓过程中，以下三个优化子问题的ε均从0开始，其中ε＝0代表辅助优化问题的解。

表1探测器初始条件及开伞条件

参数r₀V₀γ₀s₀V_f数值3.5222×10⁶5800-0.21820450单位mm/sradmm/s

具体仿真结果如下。

对于问题协状态初值随参数ε:0→1的变化如图2所示。由图2可知，协状态随参数ε的增加连续变化，一定程度上代表了求解过程的稳定性，证明了该同伦求解过程的可行性。相比于σ≡σ_min的进入轨迹，航程在ε由0变化到1时增加了约150m，出现这种现象的主要原因在于飞行时间的增加，相比于ε由0.1变化到1的绝大部分过程，最大航程问题对应的末端航程仅增加了3m，如图3所示。故在任务的快速分析设计中，可用σ≡σ_min的进入轨迹近似代替最大航程轨迹来分析最大航程。

对于问题协状态初值随参数ε:0→1的变化分别如图4所示。由图4可知，协状态随参数ε的增加连续变化，一定程度上代表了求解过程的稳定性，证明了该同伦求解过程的可行性。

最小航程轨迹比倾侧角保持为σ≡σ_max的进入轨迹末端航程高135m，如图5所示，这种现象应该是辅助问题求解时采取0.5s步长，积分精度较低引起的。在任务的快速分析设计中，可用σ≡σ_max的进入轨迹近似代替最小航程轨迹。

对于问题协状态初值随参数ε:0→1的变化分别如图6所示。由图6可知，协状态随参数ε的增加连续变化，一定程度上代表了求解过程的稳定性，证明了该同伦求解过程的可行性。由图7可知，在ε由0变化到1的过程中，相比于辅助优化问题，开伞高度由2.65km增至10.23km。

另外选取两个不同的初始航迹角：γ₀＝12°和γ₀＝13°。在分别求取三个子优化问题的基础上分别求解问题RA。其中γ₀＝12°和γ₀＝13°的子问题既可以分别从辅助优化问题出发求解，也可以从γ₀＝12.5°的三个子问题最优解出发，构造{γ_k1,γ_k2,...,γ_kN}的序列，采用同伦法求解，其中γ_k1＝12.5°，γ_kN为12°或13°，也可以是其他任务设计时感兴趣的初始航迹角。其中，γ₀＝12.5°时，求解问题RA得到的初始协状态值随航程的变化关系如图8所示，均为光滑曲线，一定程度上代表了求解过程的稳定性和方法的可行性。三个不同初始航迹角下，纵向可达区分别如图9所示。可以看出，进入角更陡时，航程范围更窄，末端可达最大开伞高度更高。

综上，本发明所提出的火星大气进入段纵向可达区生成的解析同伦法，可从最优解已知的辅助优化问题出发，最终得到火星大气进入段纵向可达区。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。