CN110562492B - 探测器火星大气进入轨迹快速生成方法 - Google Patents

Abstract

本发明公开了探测器火星大气进入轨迹快速生成方法，属于深空探测领域。本发明实现包括如下步骤：根据探测器在进入火星大气过程中的控制能力变化特点，将进入飞行过程做分段处理；建立火星大气进入段以速度为自变量的探测器无量纲动力学模型；建立复杂约束下火星大气进入轨迹生成原问题模型；将火星大气进入轨迹生成原问题转化成线性规划问题；构造序列线性规划问题，通过迭代过程求解火星大气进入轨迹生成问题。本发明由于可将约束直接转变为状态及控制量的线性约束，对于进入过程复杂约束具有较强满足能力；由于将进入轨迹生成问题转化成具有良好求解性质的序列线性规划问题，可通过内点法实现进入轨迹高效求解。

Description

探测器火星大气进入轨迹快速生成方法

技术领域

本发明属于深空探测领域，尤其涉及探测器火星大气进入轨迹快速生成方法。

背景技术

从着陆安全和科学探索价值方面考虑，NASA提出了未来火星着陆任务的关键技术是精确着陆。目前为止，已经有八颗探测器在火星表面成功着陆，这些探测器的着陆过程均沿用第一代成功着陆任务(海盗号)的着陆方案，分为大气进入段、伞降段、动力下降段和最终着陆段。大气进入段覆盖了火星EDL(进入、下降、着陆)过程的绝大部分航程，其制导控制精度对于最终的着陆精度有重要影响。目前，火星大气进入段制导控制方法研究以轨迹跟踪方法为主。2012年成功着陆火星表面的火星科学实验室任务在进入段采用的终端控制器即属于跟踪制导控制方法，该任务也是目前为止唯一在进入段采取主动控制的火星表面着陆探测任务。相比于其他任务的百公里量级着陆误差椭圆，火星科学实验室任务的着陆误差椭圆长轴约为20km，最终着陆点距离目标点约2km，进入终端控制器对于最终着陆精度的提高意义重大。

在目前的轨迹跟踪制导方案中，跟踪轨迹在离线情况下设计并存储。由于人类目前对火星了解尚不充足，任务设计过程中会存在先验信息匮乏的问题，实际飞行过程中的不确定性及扰动因素影响会导致探测器真实轨迹与离线设计的标称轨迹之间出现较大的误差。此外，目前的火星大气进入段飞行器均属于小升阻比气动构型，在进入段依靠气动力控制自身运动。考虑到火星大气比较稀薄，进入过程气动力较小，探测器的控制能力较弱，难以有效消除较大跟踪误差，极易出现控制饱和的情形，并进一步增大跟踪误差。据此，研究者提出研究进入轨迹快速生成算法，使探测器能够根据当前状态及目标开伞点状态在线更新参考轨迹，从而减弱累积跟踪误差的影响，合理利用探测器的控制能力，最终减小制导目标偏差

由于大气进入段飞行速度较快，轨迹在线生成方法对求解速度要求严格。大气进入过程约束复杂，除非线性动力学约束外，进入轨迹还需满足路径约束、开伞条件和边界条件等非线性约束。目前，学者研究的轨迹在线生成方法主要有两种：一，通过在阻力加速度-能量(D-e)剖面内对边界进行插值得到满足末端约束的进入轨迹，但插值得到的轨迹难以匹配探测器的跟踪控制能力；二，在倾侧角-能量(σ-e)剖面内生成满足条件的控制曲线并积分得到轨迹，此方法存在难以满足路径约束的困难，不能保证进入过程的安全性。

发明内容

本发明结合火星大气进入段轨迹在线生成的需求，针对现有方法难以同时较好满足复杂路径约束和控制能力约束的问题，提出一种进入轨迹快速生成的序列线性规划方法。该方法通过将进入过程做分段处理，采用速度作为进入过程自变量，将复杂非线性约束转化为线性约束，将纵向轨迹生成问题转化为序列线性规划问题，利用已有凸规划求解器，实现问题的高效求解，得到满足约束的进入轨迹。

本发明具体提供了探测器火星大气进入轨迹快速生成方法，具体包括如下步骤：

步骤1、根据探测器在进入火星大气过程中的控制能力变化特点，将进入飞行过程做分段处理，具体分为初始进入阶段、航程控制段和航向控制段；

步骤2、在火星大气进入过程中的航程控制段和航向控制段，建立以速度为自变量的探测器无量纲动力学模型；

步骤3、将探测器火星大气进入轨迹生成原问题建模成优化问题；

步骤4、将步骤3的优化问题转化成线性规划问题；

步骤5、在步骤4的基础上，构造序列线性规划问题，通过迭代过程求解探测器火星大气进入轨迹生成原问题。

步骤1包括：根据探测器在进入火星大气过程中的控制能力变化特点，将进入飞行过程的控制分为初始进入阶段、航程控制段和航向控制段三个阶段，具体包括如下步骤：

步骤1-1，初始进入阶段：进入过程开始阶段升力加速度基本为零，从控制效率和节约燃料的角度综合考虑，在火星大气进入段初期，为初始进入阶段，探测器采用常值控制量进入，所述常值控制量由地面站确定并在进入大气前发送给探测器，所述火星大气进入段初期的满足条件为：D＜0.2g_e，其中D为探测器所受阻力加速度大小，g_e为地表引力加速度；

步骤1-2，航程控制段：当D≥0.2g_e且探测器速度大于1100m/s时，进入航程控制段，探测器已经具备了一定的控制能力，根据探测器的当前状态及目标末端状态，快速生成并更新参考轨迹，以达到减小末端航程误差的目的，所述末端航程误差包括纵程误差和横程误差，在大气阻力的作用下，探测器的速度单调减小；

步骤1-3，航向控制段：当探测器的速度小于1100m/s时，航程跟踪控制的效果大大降低，进入航向控制段，此后转向航向角调整阶段，根据航向角控制律调整倾侧角的大小和方向，以达到减小横程误差的目的。

步骤1-1至步骤1-3为典型的火星大气进入制导控制过程，本发明所述的方法应用于航程控制段和航向控制段，故对于航程控制段和航向控制段，执行步骤2～步骤5。

步骤2包括如下步骤：

步骤2-1，建立如下火星大气进入段以速度为自变量的探测器无量纲动力学模型：

其中，r为火星质心到探测器质心的距离，γ为探测器的路径角，θ为探测器的火星经度，φ为探测器的火星纬度，ψ为探测器的方位角，V为探测器相对于火星表面的速度大小，L为探测器的升力加速度大小，D为探测器所受阻力加速度大小，σ为倾侧角；

各变量的无量纲化处理方法为：长度的无量纲单位为火星半径R₀＝3397200m，速度的无量纲单位为

其中

代表火星表面引力加速度，μ_M＝4.2828×10¹³m³/s²为火星引力常数，加速度大小的无量纲单位为g₀，各角度的单位均为弧度，无需做进一步处理；

步骤2-2，火星大气密度ρ采取如下指数密度模型：

其中，密度的无量纲单位为

m_scale＝7.84×10¹⁵kg为质量的无量纲单位，ρ₀＝(2.0×10^-4kg·m^-3)/ρ_scale为参考大气密度，r₀＝(3.4372×10⁶m)/R₀为参考高度，h_s＝(7.5×10³m)/R₀为标高；

步骤2-3，无量纲的阻力和升力加速度计算公式分别如下：

L＝D·(L/D)，(3)

其中，B表示弹道系数，弹道系数的无量纲单位为

L/D为升阻比大小；

步骤2-4，用

代表航程，代替公式(1)中的经度和纬度来描述探测器水平位置，其中，V₀为探测器初始速度大小，V_f为探测器在开伞点的速度大小。得到如下以速度为自变量的纵向动力学方程：

步骤2-5，将公式(4)简化写成如下紧凑形式：

其中，x＝[r γ s]^T为纵向运动状态，

代表状态x相对于速度的一阶导数，

为系统矩阵，

为控制矩阵，u＝cosσ为控制变量。

步骤3包括：

步骤3-1，确定约束项：探测器进入火星大气过程需满足的约束包括热流、过载及动压等路径约束，初始点及开伞点状态等边界约束，以及降落伞展开条件约束。此外，为了减少控制饱和情况对飞行过程的影响，设计轨迹时还需考虑控制余量约束；将探测器进入火星大气过程需满足的约束统一表达如下：

x∈[x_min,x_max] (7)

u∈[u_min,u_max] (8)

x(V₀)＝x₀,x(V_f)＝x_f (9)

q_f∈[q_f,min,q_f,max],M_f∈[M_f,min,M_f,max] (10)

其中，公式(6)为动力学约束，公式(7)为状态空间约束，公式(8)代表控制约束，公式(9)为边界条件约束，公式(10)代表开伞条件约束，M代表马赫数。公式(11)至公式(13)表示路径约束，x_min为状态下限，x_max为状态上限，u_min为控制下限，u_max为控制上限，x(V₀)为V₀时刻的状态，x₀为初始状态，x(V_f)为V_f时刻的状态，x_f为开伞点状态，q_f为开伞点动压大小，q_f,min为开伞点动压下限，q_f,max为开伞点动压上限，M_f为开伞点马赫数，M_f,min为开伞点马赫数下限，M_f,max为开伞点马赫数上限，q代表动压，q_max代表动压上限，

代表热流，

代表动压上限，n代表过载，n_max代表过载上限，

为热流系数，r_n代表探测器防热大底的半径；

公式(9)中的终端状态约束只包含末端航程约束；

步骤3-2，确定性能指标：为在开伞条件允许范围内提高开伞点高度，轨迹生成时以最大化末端高度为目标，将生成火星大气进入轨迹原问题建模成如下优化问题P0：

寻找最优控制u^*，使得：

满足约束：公式(6)至公式(13)；其中，J₀＝-r(V_f)为问题P₀对应的性能指标。

问题P0为一复杂约束下的非线性优化问题，为实现该问题的快速求解，本发明接下来将问题P0转化为一个具有良好求解性质的凸问题。

步骤4包括如下步骤：

步骤4-1，将公式(6)在参考轨迹x^ref处作如下的近似线性化处理：

公式(15)中，参考轨迹x^ref初始值由设计者给定，后续可取迭代过程中产生的轨迹为参考轨迹，具体可参考步骤5，状态矩阵A(x^ref)计算公式为：

其中，a₁₁表示矩阵A(x^ref)第一行第一列的元素；为了保证近似线性化处理过程的合理性，要求状态需满足如下附加的置信区间约束：

|x-x^ref|≤δ (23)

其中，参数

即要求每一个状态变量分别满足对应的置信区间范围约束；

公式(7)、公式(8)和公式(9)本身即为凸约束的形式，不再做凸化处理；公式(10)～公式(13)通过转换，转化为凸约束的形式，具体过程如下：

在公式(2)中的指数密度模型下，开伞条件约束中的开伞动压约束转换为对开伞高度的约束：

其中，r_f,ref为开伞时刻探测器质心距离火星质心的参考距离；

开伞马赫数约束近似转换为对末端速度的约束；

路径约束的转换如下：

等价于r≥r₁，

其中，参数r₁的计算公式如下：

q≤q_max等价于r≥r₂，其中，参数r₂的计算公式如下：

n≤n_max等价于r≥r₃，其中，参数r₃的计算公式如下：

取参数r_L＝max(r₁,r₂,r₃)，当r≥r_L时，公式(11)～公式(13)同时满足；

在纵向动力学背景下，公式(9)中的末端状态约束条件指航程约束，即：

s(V_f)＝s_f (28)

其中，s(V_f)为开伞点的航程大小，s_f为目标末端航程；

考虑到进入段探测器控制能力较弱的特点，为了保证序列求解过程中各子问题解的存在性，将航程约束进行松弛，引入新的松弛变量Θ_s、松弛系数k_s及松弛精度指标ε_s，并将性能指标调整如下：

J₁＝-r(V_f)+k_sΘ_s (29)

并满足如下约束：

|s(V_f)-s_f|≤Θ_s (30)

0≤Θ_s≤ε_s (31)

经过步骤4-1的处理，将问题P0转换成问题P1，其中，J₁为问题P1对应的性能指标，r(V_f)为开伞时刻对应的探测器质心与火星质心之间的距离大小。

由于终端速度约束已经在自变量中作为区间终点，故松弛处理过的问题P0中边界约束将只包含初始状态约束；

步骤4-2，经过步骤4-1的处理，将问题P0转换成如下问题P1：

寻找最优控制u^*，以及最优控制对应的最优状态x^*，使得：

满足约束：

x∈[x_min,x_max] (34)

u∈[u_min,u_max] (35)

x(V₀)＝x₀ (36)

r_L≤r (38)

|x-x^ref|≤δ (39)

|s(V_f)-s_f|≤Θ_s (40)

Θ_s≤ε_s (41)

问题P1的约束及性能指标为状态及控制量的线性函数形式，因此问题P1为线性规划问题。由于动力学的线性化过程是在参考轨迹附近展开的，需要迭代过程得到最优解。

步骤4-1中，公式(16)中的各元素分别为：

步骤5包括如下步骤：

问题P1与问题P0最大的区别在于动力学的线性化处理，只有当x^ref与真实轨迹x足够接近时，线性化后的动力学与真实动力学之间的误差足够小，问题P1的解与问题P0的解才能更加接近。以下通过构造序列子优化问题并依次求解的方法，来逼近问题P0的解。具体求解步骤为：

步骤5-1，k为迭代次数，k＝0时，根据公式(36)中的初始状态及常值控制量u₀，按照公式(4)中的纵向动力学方程，在速度区间[V₀,V_f]内积分得到初始轨迹x⁽⁰⁾；

步骤5-2，k＞0时，求解如下问题P2，得到{x^(k),u^(k)}：

问题P2：寻找最优控制u^*，使得

满足约束：

步骤5-3，检查以下收敛条件是否满足：

max|x^(k)-x^(k-1)|≤ε,k＞1 (44)

其中，x^(k)为第k次迭代所得的状态，x^(k-1)为第k-1次迭代所得的状态，ε为给定的精度，当公式(44)中的条件满足时，执行步骤5-4；否则，令k＝k+1并返回步骤5-2；

步骤5-4，得到问题P0的近似最优解x^*＝x^(k)，u^*＝u^(k)，u^(k)为第k次迭代所得的控制量大小，对问题P2离散化，得到最终结果。

步骤5-4包括：

步骤5-4-1，将速度区间[V₀,V_f]均分为N个区间，步长为ΔV＝(V_f-V₀)/N，离散的节点为{V₀,V₁,...V_N}，相应的状态和控制量分别离散为序列{x₀,x₁,...,x_N}和序列{u₀,u₁,...,u_N}，其中，V_N为第N个节点处的速度大小，x_N为第N个节点处的状态变量，u_N为第N个节点处的控制量大小，运用梯形积分公式，动力学约束表示如下：

其中，x_i为第i个节点处的状态，x_i-1为第i-1个节点处的状态，

为第k-1次迭代中第i个节点处的状态向量，

为第k-1次迭代中第i个节点处的系统矩阵，

为第k-1次迭代中第i个节点处的状态矩阵，

为第k-1次迭代中第i个节点处的控制矩阵；

步骤5-4-2，将公式(45)做如下处理：

其中，

为第k-1次迭代中第i个节点处的矩阵系数，

为第k-1次迭代中第i-1个节点处的矩阵系数，I为单位矩阵；

步骤5-4-3，需要优化设计的变量是{x₀,x₁,...,x_N}、{u₀,u₁,...,u_N}以及松弛变量Θ_s，即待优化变量y＝{x₀,x₁,...,x_N,u₀,u₁,...,u_N,Θ_s}，问题P2转化为如下线性规划问题P3：

寻找最优变量y^*，使得：

满足约束：

E_iy＝0,i＝1,2,...,N_c1 (48)

其中，J₃为问题P3对应的性能指标，当(E_i,M_j,P_j)取不同的值时，公式(48)和公式(49)即能够转化为公式(43)中的具体约束形式；

公式(47)中，

为性能指标系数；公式(48)中，

为第i个等式约束系数，n_i为矩阵行数；公式(49)中，

为第i个不等式约束系数，

为约束矩阵，n_j为矩阵行数，N_c1为线性等式约束的总个数，N_c2为线性不等式约束的总个数。区间N的个数越大，离散问题与原问题越接近，问题P3属于线性规划问题，可运用内点法在多项式时间内求解。

有益效果

1、本发明的一种星大气进入轨迹快速生成的序列线性规划方法，由于可将约束直接转变为状态及控制量的线性约束，对于进入过程复杂约束具有较强满足能力；

2、本发明的一星大气进入轨迹快速生成的序列线性规划方法，由于将进入轨迹生成问题转化成具有良好求解性质的序列线性规划问题，可通过内点法实现进入轨迹高效求解。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1为本发明方法流程图；

图2为控制曲线随进入速度的变化关系；

图3为开伞高度随迭代次数的变化关系；

图4为不同迭代次数下探测器高度随进入速度的变化关系图；

图5为不同迭代次数下探测器过载随进入速度的变化关系图；

图6为不同迭代次数下探测器动压随进入速度的变化关系图；

图7为不同迭代次数下探测器热流随进入速度的变化关系图。

具体实施方式

本实例为针对火星大气进入段轨迹跟踪制导方案中的参考轨迹快速生成问题，在性能指标及约束分析的基础上建立轨迹生成问题原模型，通过将轨迹生成问题转化成序列线性规划问，利用内点法，实现复杂约束下进入轨迹快速生成。

本发明具体提供了探测器火星大气进入轨迹快速生成方法，如图1所示，具体步骤如下：

步骤1、根据探测器在进入火星大气过程中的控制能力变化特点，将进入飞行过程做分段处理；

根据探测器在进入火星大气过程中的控制能力变化特点，将进入过程的控制分为初始进入阶段、航程控制段和航向控制段三个阶段。具体分段方法如下：

1)初始进入段：进入过程开始阶段升力加速度基本为零，从控制效率和节约燃料的角度综合考虑，在火星大气进入段初期(满足条件：D＜0.2g_e，其中D为探测器所受阻力加速度大小，g_e为地表引力加速度)，探测器采用常值倾侧角45°进入。

2)航程控制段：当D≥0.2g_e且速度大于1100m/s时，探测器已经具备了一定的控制能力，制导目标为减小末端航程误差，包括纵程误差和横程误差。在大气阻力的作用下，探测器的速度单调减小。

3)航向控制段：当探测器的速度小于1100m/s时，航程跟踪控制的效果大大降低，此后转向航向角调整阶段，制导控制目标转变为减小横程误差。

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

以速度为自变量构建大气进入段三自由度动力学方程及相应的纵向运动方程。以速度为自变量的探测器三自由度无量纲动力学方程如下：

其中，r为火星质心到探测器质心的距离，γ为探测器的路径角，θ为探测器的火星经度，φ为探测器的火星纬度，ψ为探测器的方位角，V为探测器相对于火星表面的速度大小，L为探测器的升力加速度大小，D为探测器所受阻力加速度大小，σ为倾侧角。各变量的无量纲化处理方法为：长度的无量纲单位为火星半径R₀＝3397200m，速度的无量纲单位为

其中

代表火星表面引力加速度，μ_M＝4.2828×10¹³m³/s²为火星引力常数，加速度大小的无量纲单位为g₀，各角度的单位均为弧度，无需做进一步处理。

大气密度采取如下指数密度模型：

其中，密度的无量纲单位为

m_scale＝7.84×10¹⁵kg为质量的无量纲单位，ρ₀＝(2.0×10^-4kg·m^-3)/ρ_scale为参考大气密度，r₀＝(3.4372×10⁶m)/R₀为参考高度，h_s＝(7.5×10³m)/R₀为标高。

无量纲的阻力和升力加速度分别如下：

其中，弹道系数B＝(0.00685m²·kg^-1)/B_scale，弹道系数的无量纲单位为

L/D＝0.24为升阻比大小。

用

代表航程，代替式(1)中的经度和纬度来描述探测器水平位置，得到以速度为自变量的纵向动力学方程如下：

为简化问题描述，将上式写成如下紧凑形式：

其中，x＝[r γ s]^T为纵向运动状态，

u＝cosσ为控制变量。

步骤3、建立火星大气进入轨迹生成原问题模型；

将火星大气进入轨迹生成问题建模为优化问题，下面将分别确定该优化问题的约束项及性能指标。

1)确定约束项。进入过程需满足的约束主要为热流、过载及动压等路径约束，初始点及开伞点状态等边界约束，以及降落伞展开条件约束。此外，为了减少控制饱和情况对飞行过程的影响，设计轨迹时还需考虑控制余量约束。各约束可统一表达如下：

x∈[x_min,x_max] (7)

u∈[u_min,u_max] (8)

x(V₀)＝x₀,x(V_f)＝x_f (9)

q_f∈[q_f,min,q_f,max],M_f∈[M_f,min,M_f,max] (10)

其中，式(6)为动力学约束，式(7)为状态空间约束，式(8)代表控制约束，式(9)为边界条件约束，式(10)代表开伞条件约束，M代表马赫数。式(11)至式(13)表示路径约束。式(9)中的终端状态约束只包含末端航程约束。下标“0”代表火星大气进入点，下标“f”代表开伞点，下标“min”和“max”分别代表约束的下限和上限。q代表动压，M代表马赫数，

代表热流，n代表过载，

为热流系数，r_n代表探测器防热大底的半径。边界约束值具体如表1所示，开伞条件、控制量及路径约束如表2所示。

表1

状态	大气进入点值	开伞点值	单位
				高度	125	自由，待最大化	km
速度	5500	450	m/s
				路径角	-15.2	自由	度
经度	90	101.5361	度
				纬度	-28	-28.4038	度
方位角	95	自由	度

表2

状态	最大值	最小值	单位
				开伞马赫数	2.5	1.1	km
开伞动压	180	800	Pa
				过载	12	-	g<sub>e</sub>
动压	23	-	kPa
				热流	90	-	W/cm<sup>2</sup>
倾侧角大小	120	30	度

2)确定性能指标。为在开伞条件允许范围内提高开伞点高度，轨迹生成时以最大化末端高度为目标，取性能指标如式(14)所示。

综上，可将轨迹生成问题建模成问题P0：

寻找最优控制u^*，使得：

满足约束：式(6)至式(13)

问题P0为一复杂约束下的非线性优化问题，为实现该问题的快速求解，本发明接下来将问题P0转化为一个具有良好求解性质的凸问题。

步骤4、将火星大气进入轨迹生成原问题转化成线性规划问题；

将式(6)在参考轨迹x^ref处作如下的近似线性化处理：

其中，上标“ref”代表参考轨迹。式(15)中，

式(16)中的各系数分别为：

为了保证上述线性化过程的合理性，要求状态需满足如下附加的置信区间约束：

|x-x^ref|≤δ (23)

式中，

即要求每一个状态变量分别满足对应的置信区间范围约束。本实施案例中，取δ＝[10000/R₀ 20π/180 0.2s_togo,0]^T。

式(7)、式(8)和式(9)本身即为凸约束的形式，不必再做凸化处理。式(10)～式(13)可通过简单的变换，转化为凸约束的形式，具体过程如下：

在式(2)中的指数密度模型下，开伞动压约束可转换为对开伞高度的约束：

开伞马赫数约束可近似转换为对末端速度的约束。由于本发明中的末端速度是给定的，且落在马赫数要求的范围内，故该条件自动满足。路径约束的转换如下：

等价于r≥r₁，其中：

q≤q_max等价于r≥r₂，其中：

n≤n_max等价于r≥r₃，其中：

取r_L＝max(r₁,r₂,r₃)，当r≥r_L时，式(11)～式(13)同时满足。

在纵向动力学背景下，式(9)中的末端状态约束条件主要指航程约束，即：

s(V_f)＝s_f (28)

考虑到进入段探测器控制能力较弱的特点，为了保证序列求解过程中各子问题解的存在性，此处将末端航程约束进行松弛，引入新的松弛变量Θ_s、松弛系数k_s及松弛精度指标ε_s，并将性能指标调整如下：

J₁＝-r(V_f)+k_sΘ_s (29)

其中，

|s(V_f)-s_f|≤Θ_s (30)

0≤Θ_s≤ε_s (31)

由于终端速度约束已经在自变量中作为区间终点，故松弛处理过的问题P0中边界约束将只包含初始状态约束。

综上，可将问题P0转换成问题P1：

寻找最优控制u^*，使得：

满足约束：

x∈[x_min,x_max] (34)

u∈[u_min,u_max] (35)

x(V₀)＝x₀ (36)

r_L≤r (38)

|x-x^ref|≤δ (39)

|s(V_f)-s_f|≤Θ_s (40)

Θ_s≤ε_s (41)

可以看出，问题P1的约束及性能指标为状态及控制量的线性函数形式，故该问题为线性规划问题。由于动力学的线性化过程是在参考轨迹附近展开的，需要迭代过程得到最优解。

步骤5、构造序列线性规划问题，通过迭代过程求解火星大气进入轨迹生成问题。

问题P1与问题P0最大的区别在于动力学的线性化处理，只有当x^ref与真实轨迹x足够接近时，线性化后的动力学与真实动力学之间的误差足够小，问题P1的解与问题P0的解才能更加接近。以下通过构造序列子优化问题并依次求解的方法，来逼近问题P0的解。具体求解步骤为：

1)k＝0时，根据式(36)中的初始状态及常值控制量u₀＝cos(50°)，按照式(4)中的动力学，在速度区间[V₀,V_f]内积分得到初始轨迹x⁽⁰⁾。

2)k＞0时，求解如下问题P2，得到{x^(k),u^(k)}。

问题P2：寻找最优控制u^*，使得

满足约束：

3)检查以下收敛条件是否满足：

max|x^(k)-x^(k-1)|≤ε,k＞1 (44)

式中，ε＝[(100m)/R₀0.05π/180(100m)/R₀]^T为给定的精度，k为迭代次数，上标“k-1”、“k”分别代表第k-1、第k次迭代。当式(44)中的条件满足时，执行“4)”；否则，令k＝k+1并返回“2)”。

4)得到问题P0的近似最优解x^*＝x^(k)，u^*＝u^(k)。

为应用数值求解器求解问题P2，需先对其做离散化。具体离散化过程为：将速度区间[V₀,V_f]均分为N＝400个区间，步长为ΔV＝(V_f-V₀)/N，离散的节点为{V₀,V₁,...V_N}。相应的状态和控制量分别离散为序列{x₀,x₁,...,x_N}和序列{u₀,u₁,...,u_N}。运用梯形积分公式，动力学约束可表示如下：

其中，下标中的i和i+1表示第i和i+1个节点，上标k-1表示第k-1次的迭代结果。

为方便数值处理，将式(45)做如下处理：

其中，

需要优化设计的变量是{x₀,x₁,...,x_N}、{u₀,u₁,...,u_N}以及松弛变量Θ_s，即y＝{x₀,x₁,...,x_N,u₀,u₁,...,u_N,Θ_s}。问题P2可转化为如下线性规划问题P3：

寻找最优变量y^*，使得：

满足约束：

E_iy＝0,i＝1,2,...,N_c1 (48)

当(E_i,M_j,P_j)取不同的值时，式(48)和式(49)即可转化为式(43)中的具体约束形式。式(47)中，

式(48)中，

式(49)中，

N_c1为线性等式约束的总个数，N_c2为线性不等式约束的总个数。区间N的个数越大，离散问题与原问题越接近。问题P3属于线性规划问题，本实施案例中求解方法采用CVX中的ECOS求解器，程序的运行环境为Interl(R)Core(TM)i7-4790，频率为3.6GHz。

图2至图7为仿真结果。其中，图2为控制曲线随进入速度的变化关系，其中50°常值倾侧角对应初始轨迹。由局部放大图可以看出轨迹的具体收敛过程。第一次迭代得到的控制结果与最优结果虽然接近，但仍不是最优解。而第二次和第三次的得到的最优解已十分接近，可认为该迭代求解过程已经收敛；图3为开伞高度随迭代次数的变化关系；可以看出最优轨迹对应的末端高度为13.89km；图4为不同迭代次数下探测器高度随进入速度的变化关系图，可以看出，轨迹从初始参考轨迹开始进行第一次迭代时，高度-速度曲线已与最优的高度-速度曲线十分接近。在精度要求更低的情况下，迭代收敛的步数会更少，轨迹生成的速度更快；图5至图7为初始轨迹与迭代轨迹的路径约束满足情况，相对于初始轨迹，迭代得到的最优轨迹可以满足任务要求的过载、热流以及动压约束。而初始轨迹的峰值过载超过了过载上限约束，威胁探测器的安全。

表3所示为迭代过程中相邻轨迹之间的各节点状态误差最大值及求解时间。

表3

迭代次数	求解时间/s	|Δr|,km	|Δγ|,(°)	|Δs|,km
					1	0.6240	1.5067	2.5588	4.9374
2	0.2650	0.0687	0.4375	0.4336
					3	0.2660	0.0024	0.0221	0.0081

综上，本发明所提出的一种探测器火星大气进入轨迹快速生成方法可以在短时间内生成复杂约束下的进入轨迹，具备在线实时运用的潜力。

本发明提供了探测器火星大气进入轨迹快速生成方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (6)

1.探测器火星大气进入轨迹快速生成方法，其特征在于：包括如下步骤，

步骤1、根据探测器在进入火星大气过程中的控制能力变化特点，将进入飞行过程做分段处理，具体分为初始进入阶段、航程控制段和航向控制段；

步骤2、在火星大气进入过程中的航程控制段和航向控制段，建立以速度为自变量的探测器无量纲动力学模型；

步骤3、将探测器火星大气进入轨迹生成原问题建模成优化问题；

步骤4、将步骤3的优化问题转化成线性规划问题；

步骤5、在步骤4的基础上，构造序列线性规划问题，通过迭代过程求解探测器火星大气进入轨迹生成原问题；

步骤1包括：根据探测器在进入火星大气过程中的控制能力变化特点，将进入飞行过程的控制分为初始进入阶段、航程控制段和航向控制段三个阶段，具体包括如下步骤：

步骤1-1，初始进入阶段：在火星大气进入段初期，为初始进入阶段，探测器采用常值控制量进入，所述常值控制量由地面站确定并在进入大气前发送给探测器，所述火星大气进入段初期的满足条件为：D＜0.2g_e，其中D为探测器所受阻力加速度大小，g_e为地表引力加速度；

步骤1-2，航程控制段：当D≥0.2g_e且探测器速度大于1100m/s时，进入航程控制段；

步骤1-3，航向控制段：当探测器的速度小于1100m/s时，进入航向控制段；

步骤2包括如下步骤：

步骤2-1，建立如下火星大气进入段以速度为自变量的探测器无量纲动力学模型：

其中，r为火星质心到探测器质心的距离，γ为探测器的路径角，θ为探测器的火星经度，φ为探测器的火星纬度，ψ为探测器的方位角，V为探测器相对于火星表面的速度大小，L为探测器的升力加速度大小，D为探测器所受阻力加速度大小，σ为倾侧角；

各变量的无量纲化处理方法为：长度的无量纲单位为火星半径R₀＝3397200m，速度的无量纲单位为

其中

代表火星表面引力加速度，μ_M＝4.2828×10¹³m³/s²为火星引力常数，加速度大小的无量纲单位为g₀，各角度的单位均为弧度，无需做进一步处理；

步骤2-2，火星大气密度ρ采取如下指数密度模型：

其中，密度的无量纲单位为

m_scale＝7.84×10¹⁵kg为质量的无量纲单位，ρ₀＝(2.0×10^-4kg·m^-3)/ρ_scale为参考大气密度，r₀＝(3.4372×10⁶m)/R₀为参考高度，h_s＝(7.5×10³m)/R₀为标高；

步骤2-3，无量纲的阻力和升力加速度计算公式分别如下：

L＝D·(L/D)， (3)

其中，B表示弹道系数，弹道系数的无量纲单位为

L/D为升阻比大小；

步骤2-4，用

代表航程，代替公式(1)中的经度和纬度来描述探测器水平位置，其中，V₀为探测器初始速度大小，V_f为探测器在开伞点的速度大小；得到如下以速度为自变量的纵向动力学方程：

步骤2-5，将公式(4)简化写成如下紧凑形式：

(5)其中，x＝[r γ s]^T为纵向运动状态，

代表状态x相对于速度的一阶导数，

为系统矩阵，

为控制矩阵，u＝cosσ为控制变量。

2.如权利要求1所述的方法，其特征在于，步骤3包括：

步骤3-1，确定约束项：将探测器进入火星大气过程需满足的约束统一表达如下：

x∈[x_min,x_max] (7)

u∈[u_min,u_max] (8)

x(V₀)＝x₀,x(V_f)＝x_f (9)

q_f∈[q_f,min,q_f,max],M_f∈[M_f,min,M_f,max] (10)

其中，公式(6)为动力学约束，公式(7)为状态空间约束，公式(8)代表控制约束，公式(9)为边界条件约束，公式(10)代表开伞条件约束，M代表马赫数；公式(11)至公式(13)表示路径约束，x_min为状态下限，x_max为状态上限，u_min为控制下限，u_max为控制上限，x(V₀)为V₀时刻的状态，x₀为初始状态，x(V_f)为V_f时刻的状态，x_f为开伞点状态，q_f为开伞点动压大小，q_f,min为开伞点动压下限，q_f,max为开伞点动压上限，M_f为开伞点马赫数，M_f,min为开伞点马赫数下限，M_f,max为开伞点马赫数上限，q代表动压，q_max代表动压上限，

代表热流，

代表动压上限，n代表过载，n_max代表过载上限，

为热流系数，r_n代表探测器防热大底的半径；

公式(9)中的终端状态约束只包含末端航程约束；

步骤3-2，确定性能指标：为在开伞条件允许范围内提高开伞点高度，轨迹生成时以最大化末端高度为目标，将生成火星大气进入轨迹原问题建模成如下优化问题P0：

寻找最优控制u^*，使得：

满足约束：公式(6)至公式(13)，其中，J₀＝-r(V_f)为问题P₀对应的性能指标。

3.如权利要求2所述的方法，其特征在于，步骤4包括如下步骤：

步骤4-1，将公式(6)在参考轨迹x^ref处作如下的近似线性化处理：

公式(15)中，状态矩阵A(x^ref)计算公式为：

其中，a₁₁表示矩阵A(x^ref)第一行第一列的元素；为了保证近似线性化处理过程的合理性，要求状态需满足如下附加的置信区间约束：

|x-x^ref|≤δ (23)

其中，参数

即要求每一个状态变量分别满足对应的置信区间范围约束；

公式(7)、公式(8)和公式(9)本身即为凸约束的形式，不再做凸化处理；公式(10)～公式(13)通过转换，转化为凸约束的形式，具体过程如下：

在公式(2)中的指数密度模型下，开伞条件约束中的开伞动压约束转换为对开伞高度的约束：

其中，r_f,ref为开伞时刻探测器质心距离火星质心的参考距离；

开伞马赫数约束近似转换为对末端速度的约束；

路径约束的转换如下：

等价于r≥r₁，其中，参数r₁的计算公式如下：

q≤q_max等价于r≥r₂，其中，参数r₂的计算公式如下：

n≤n_max等价于r≥r₃，其中，参数r₃的计算公式如下：

取参数r_L＝max(r₁,r₂,r₃)，当r≥r_L时，公式(11)～公式(13)同时满足；

在纵向动力学背景下，公式(9)中的末端状态约束条件指航程约束，即：

s(V_f)＝s_f (28)

其中，s(V_f)为开伞点的航程大小，s_f为目标末端航程；

将航程约束进行松弛，引入新的松弛变量Θ_s、松弛系数k_s及松弛精度指标ε_s，并将性能指标调整如下：

J₁＝-r(V_f)+k_sΘ_s (29)

并满足如下约束：

|s(V_f)-s_f|≤Θ_s (30)

0≤Θ_s≤ε_s (31)

经过步骤4-1的处理，将问题P0转换成问题P1，其中，J₁为问题P1对应的性能指标，r(V_f)为开伞时刻对应的探测器质心与火星质心之间的距离大小；

由于终端速度约束已经在自变量中作为区间终点，则松弛处理过的问题P0中边界约束将只包含初始状态约束；

步骤4-2，经过步骤4-1的处理，问题P0转换成如下问题P1：

寻找最优控制u^*，以及最优控制对应的最优状态x^*，使得：

满足约束：

x∈[x_min,x_max] (34)

u∈[u_min,u_max] (35)

x(V₀)＝x₀ (36)

r_L≤r (38)

|x-x^ref|≤δ (39)

|s(V_f)-s_f|≤Θ_s (40)

Θ_s≤ε_s (41)

问题P1的约束及性能指标为状态及控制量的线性函数形式，因此问题P1为线性规划问题。

4.如权利要求3所述的方法，其特征在于，步骤4-1中，公式(16)中的各元素分别为：

5.如权利要求4所述的方法，其特征在于，步骤5包括如下步骤：

步骤5-1，k为迭代次数，k＝0时，根据公式(36)中的初始状态及常值控制量u₀，按照公式(4)中的纵向动力学方程，在速度区间[V₀,V_f]内积分得到初始轨迹x⁽⁰⁾；

步骤5-2，k＞0时，求解如下问题P2，得到{x^(k),u^(k)}：

问题P2：寻找最优控制u^*，使得

满足约束：

步骤5-3，检查以下收敛条件是否满足：

max|x^(k)-x^(k-1)|≤ε,k＞1 (44)

其中，x^(k)为第k次迭代所得的状态，x^(k-1)为第k-1次迭代所得的状态，ε为给定的精度，当公式(44)中的条件满足时，执行步骤5-4；否则，令k＝k+1并返回步骤5-2；

步骤5-4，得到问题P0的近似最优解x^*＝x^(k)，u^*＝u^(k)，u^(k)为第k次迭代所得的控制量大小，对问题P2离散化，得到最终结果。

6.如权利要求5所述的方法，其特征在于，步骤5-4包括：

步骤5-4-1，将速度区间[V₀,V_f]均分为N个区间，步长为ΔV＝(V_f-V₀)/N，离散的节点为{V₀,V₁,...V_N}，相应的状态和控制量分别离散为序列{x₀,x₁,...,x_N}和序列{u₀,u₁,...,u_N}，其中，V_N为第N个节点处的速度大小，x_N为第N个节点处的状态变量，u_N为第N个节点处的控制量大小，运用梯形积分公式，动力学约束表示如下：

其中，x_i为第i个节点处的状态，x_i-1为第i-1个节点处的状态，

为第k-1次迭代中第i个节点处的状态向量，

为第k-1次迭代中第i个节点处的系统矩阵，

为第k-1次迭代中第i个节点处的状态矩阵，

为第k-1次迭代中第i个节点处的控制矩阵；

步骤5-4-2，将公式(45)做如下处理：

其中，

为第k-1次迭代中第i个节点处的矩阵系数，

为第k-1次迭代中第i-1个节点处的矩阵系数，I为单位矩阵；

步骤5-4-3，需要优化设计的变量是{x₀,x₁,...,x_N}、{u₀,u₁,...,u_N}以及松弛变量Θ_s，即待优化变量y＝{x₀,x₁,...,x_N,u₀,u₁,...,u_N,Θ_s}，问题P2转化为如下线性规划问题P3：

寻找最优变量y^*，使得：

满足约束：

E_iy＝0,i＝1,2,...,N_c1 (48)

其中，J₃为问题P3对应的性能指标，当(E_i,M_j,P_j)取不同的值时，公式(48)和公式(49)即能够转化为公式(43)中的具体约束形式；

公式(47)中，

为性能指标系数；公式(48)中，

为第i个等式约束系数，n_i为矩阵行数；公式(49)中，

为第i个不等式约束系数，

为约束矩阵，n_j为矩阵行数，N_c1为线性等式约束的总个数，N_c2为线性不等式约束的总个数。

Images