CN109250153B - 火星大气进入段轨迹最优跟踪制导方法 - Google Patents

Abstract

本发明公开的火星大气进入段轨迹最优跟踪制导方法，属于深空探测领域。本发明具体实现方法包括如下步骤：步骤一、建立火星大气进入段探测器纵向动力学模型；步骤二、通过对路径约束和控制约束凸化处理，选取二次型性能指标，建立轨迹最优跟踪问题模型；步骤三、对步骤二中的轨迹最优跟踪问题模型离散化处理，将其转化成二次规划问题，使所述问题能够用内点法在多项式时间内求解，利用数值方法求解离散化的二次规划问题，提高求解效率，进而实现火星大气进入段实时在线制导。本发明将轨迹跟踪制导问题转化为最优控制问题进行求解，能够提高开伞精度。

Description

火星大气进入段轨迹最优跟踪制导方法

技术领域

本发明涉及一种火星大气进入段轨迹最优跟踪制导方法，属于深空探测领域。

背景技术

从着陆安全和科学探索价值方面考虑，NASA提出了未来火星着陆任务的关键技术是精确着陆。目前为止，已经有七颗探测器在火星表面成功着陆，这些探测器的着陆过程均沿用海盗号的着陆方案，分为大气进入段、伞降段、动力下降段和最终着陆段。大气进入段是着陆过程的重要阶段，其导航、制导精度对最终着陆精度有决定性影响。在以往其次着陆探测任务中，仅有火星科学实验室任务在大气进入段采用了闭环制导控制方案，并将着陆精度从以往的百公里量级提升至十公里量级，但与精确着陆的目标仍相距甚远。火星科学实验室任务在大气进入段采用的是离线设计、存储轨迹及相关控制系数，并在线跟踪轨迹的制导方案。由于火星大气进入段不确定性和扰动多，离线设计的跟踪控制系数难以适应实际飞行情况，易造成较大的跟踪误差以及控制饱和的情况出现，影响最终着陆精度及安全性。

考虑到机载计算机性能的日益提升，计算制导成为新一代进入段制导算法的主流，代表之一是大气进入段预测-修正制导方法。由于预测-修正制导方法需要对动力学反复积分，降低了制导律求解效率，且在过程约束满足方面，预测-修正制导方法仍缺乏完善的技术，特别是对于火星大气进入段的小升阻比探测器。凸优化方法的全局最优性及计算效率使其具备应用于在线制导的潜力。目前，深空领域的凸优化研究包括火星动力下降段精确着陆问题、接近与交会中的轨迹优化问题以及大气进入过程的轨迹优化问题等。

发明内容

为了解决在火星大气进入段存在不确定性和扰动情况下的探测器精确开伞问题。本发明公开的火星大气进入段轨迹最优跟踪制导方法要解决的技术问题是：提供一种基于凸规划的火星大气进入段参考轨迹最优跟踪制导方法，通过将性能指标、路径约束和控制约束处理成凸约束的形式，将轨迹跟踪问题转化为二次规划问题，使所述问题能够用内点法在多项式时间内求解，求解效率高，能够实现实时在线制导。

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

本发明公开的火星大气进入段轨迹最优跟踪制导方法，建立火星大气进入段探测器纵向动力学模型；通过对路径约束和控制约束凸化处理，选取二次型性能指标，建立轨迹最优跟踪问题模型；对轨迹最优跟踪问题模型离散化处理，将其转化成二次规划问题，利用数值方法求解离散化的二次规划问题，提高求解效率，进而实现火星大气进入段实时在线制导。

本发明公开的火星大气进入段轨迹最优跟踪制导方法，包括如下步骤：

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

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

式(1)中，V为探测器速度大小，u＝cosσ为控制量，σ为倾侧角。在无量纲化过程中，长度的量纲单位为火星半径R₀，速度的无量纲单位为

其中

为火表引力加速度，μ为火星引力常数。时间的无量纲单位为

角度的单位为弧度，不需要无量纲化处理。式(1)中，L和D分别为探测器受到的无量纲升力和阻力加速度，分别具有如下形式：

L＝D·L/D (3)

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

其中，ρ₀为参考密度，r₀为参考半径，h_s为标高。

步骤二、通过对路径约束和控制约束凸化处理，选取二次型性能指标，建立轨迹最优跟踪问题模型；

在小扰动下，状态误差的动力学方程如下：

式中，上标“*”代表参考轨迹与参考控制。A(V)为雅克比矩阵，B_u(V)为式(1)中的动力学方程对控制量的偏导数。

大气进入段跟踪制导需满足路径约束，即动压、过载与热流约束。由于动压与过载成比例关系，此处将动压约束与过载约束合并，则路径约束如下：

其中，

为最小动压约束，q为动压，

为热流，n为过载，下标“max”代表上限约束。

除了步骤一中动力学微分方程与路径约束外，大气进入段跟踪制导还应满足边界条件：

δx(V₀)＝δx₀ (9)

其中δx₀为当前制导周期开始时刻，实际状态与标称状态之间的差值。为了保证线性化的有效性，还需满足信赖域条件：

|δx(V)|≤ε_x (10)

|δu(V)|≤ε_u (11)

除此之外，还应满足状态上下限约束及控制能力约束，即

x_min(V)≤x^*(V)+δx(V)≤x_max(V) (12)

u_min(V)≤u^*(V)+δu(V)≤u_max(V) (13)

为了使用最少的燃料消耗实现轨迹跟踪，在每个制导周期内，取性能指标为：

其中，V₀代表当前制导周期的开始时刻，V_f代表开伞时刻。P和Q(V)为对称半正定矩阵，R(V)为一正数。将轨迹跟踪问题建模成最优控制问题Γ0。

问题Γ0：寻找最优控制δu^*(V)，使得

满足约束：

式(5)、式(8)～式(13)

最优控制问题Γ0即为建立的轨迹最优跟踪问题模型。

作为优选，由于大气进入段探测器的控制能力较弱，探测器难以实现全状态跟踪，为了能达到末端航程误差最小，式(14)中的性能指标取为如公式(16)所示形式：

步骤三、对步骤二中的轨迹最优跟踪问题模型离散化处理，将其转化成二次规划问题，利用数值方法求解离散化的二次规划问题，提高求解效率，进而实现实时在线制导。

为了求得问题Γ0的数值解，将速度区间[V₀,V_f]均分为N个区间，每个区间的长度为ΔV＝(V_f-V₀)/N，将式(5)、式(8)～式(13)中所示的约束施加于N+1个节点上。用单调递减序列{V₀,V₁,...,V_N}表示自变量节点，其中V_i＝V₀+iΔV,i＝1,...,N，对应的状态和控制序列分别为{δx₀,δx₁,...,δx_N}和{δu₀,δu₁,...,δu_N}。运用梯形积分公式，动力学约束近似如下：

式中，A_x,i＝A_x(x(V_i),u(V_i))，B_u,i＝B_u(x(V_i),u(V_i))。上式进一步表示为：

式中，

I为单位矩阵。

将优化变量统一用向量δz表示，即：

δz＝[δx₀,δu₀,δx₁,δu₁,...,δx_N,δu_N]^T (19)

问题Γ0的离散形式统一用如公式(20)、(21)所示形式表示：

满足约束：

其中，N_C为总约束数量，C_κ、W_k和l_k的值根据参考轨迹计算。故问题Γ0能够转化成二次规划问题，利用数值方法求解所述离散化的二次规划问题，提高求解效率，进而实现火星大气进入段实时在线制导。

有益效果：

1、本发明公开的火星大气进入段轨迹最优跟踪制导方法，通过对路径约束和控制约束凸化处理，并将轨迹最优跟踪问题模型进行离散化处理，使其适用于二次规划问题。

2、本发明公开的火星大气进入段轨迹最优跟踪制导方法，由于将轨迹跟踪问题建模成二次规划问题，并利用数值方法求解离散化的二次规划问题，提高求解效率，进而实现实时在线制导。

3、本发明公开的火星大气进入段轨迹最优跟踪制导方法，由于将轨迹跟踪制导问题转化为最优控制问题进行求解，提高开伞精度。

附图说明

图1为火星大气进入段轨迹最优跟踪制导方法流程图。

图2为不同案例下倾侧角修正量随速度的变化关系。

图3为末端航程跟踪下开伞点航程误差分布。

图4为开环控制下开伞点航程误差分布。

图5为案例1的制导律求解时间。

具体实施方式

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

本实例为针对火星大气进入段轨迹跟踪问题，基于小扰动假设及凸规划理论，利用SeDuMi，对轨迹跟踪过程中的控制修正量进行优化。

本实施例公开的火星大气进入段轨迹最优跟踪制导方法，如图1所示，具体步骤如下：

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

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

式(1)中，V为探测器速度大小，u＝cosσ为控制量，σ为倾侧角。在无量纲化过程中，长度的量纲单位为火星半径R₀，速度的无量纲单位为

其中

为火表引力加速度，μ为火星引力常数。时间的无量纲单位为

角度的单位为弧度，不需要无量纲化处理。式(1)中，L和D分别为探测器受到的无量纲升力和阻力加速度，分别具有如下形式：

L＝D·L/D (3)

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

其中，ρ₀为参考密度，r₀为参考半径，h_s为标高。

步骤二、通过对路径约束和控制约束凸化处理，选取二次型性能指标，建立轨迹最优跟踪问题模型；

在小扰动下，状态误差的动力学方程如下：

式中，上标“*”代表参考轨迹与参考控制。A(V)为雅克比矩阵，B_u(V)为式(1)中的动力学方程对控制量的偏导数。

大气进入段跟踪制导需满足路径约束，即动压、过载与热流约束。由于动压与过载成比例关系，此处将动压约束与过载约束合并，则路径约束如下：

其中，

为最小动压约束，q为动压，

为热流，n为过载，下标“max”代表上限约束。

除了步骤一中动力学微分方程与路径约束外，大气进入段跟踪制导还应满足边界条件：

δx(V₀)＝δx₀ (9)

其中δx₀为当前制导周期开始时刻，实际状态与标称状态之间的差值。为了保证线性化的有效性，还需满足信赖域条件：

|δx(V)|≤ε_x (10)

|δu(V)|≤ε_u (11)

除此之外，还应满足状态上下限约束及控制能力约束，即

x_min(V)≤x^*(V)+δx(V)≤x_max(V) (12)

u_min(V)≤u^*(V)+δu(V)≤u_max(V) (13)

为了使用最少的燃料消耗实现轨迹跟踪，在每个制导周期内，取性能指标为：

其中，V₀代表当前制导周期的开始时刻，V_f代表开伞时刻。P和Q(V)为对称半正定矩阵，R(V)为一正数。将轨迹跟踪问题建模成最优控制问题Γ0。

问题Γ0：寻找最优控制δu^*(V)，使得

满足约束：

式(5)、式(8)～式(13)

最优控制问题Γ0即为建立的轨迹最优跟踪问题模型。

由于大气进入段探测器的控制能力较弱，探测器难以实现全状态跟踪，为了能达到末端航程误差最小，式(14)中的性能指标取为如公式(16)所示形式：

步骤三、对步骤二中的轨迹最优跟踪问题模型离散化处理，将其转化成二次规划问题，利用数值方法求解离散化的二次规划问题，提高求解效率，进而实现实时在线制导。

为了求得问题Γ0的数值解，将速度区间[V₀,V_f]均分为N个区间，每个区间的长度为ΔV＝(V_f-V₀)/N，将式(5)、式(8)～式(13)中所示的约束施加于N+1个节点上。用单调递减序列{V₀,V₁,...,V_N}表示自变量节点，其中V_i＝V₀+iΔV,i＝1,...,N，对应的状态和控制序列分别为{δx₀,δx₁,...,δx_N}和{δu₀,δu₁,...,δu_N}。运用梯形积分公式，动力学约束近似如下：

式中，A_x,i＝A_x(x(V_i),u(V_i))，B_u,i＝B_u(x(V_i),u(V_i))。上式进一步表示为：

式中，

I为单位矩阵。

将优化变量统一用向量δz表示，即：

δz＝[δx₀,δu₀,δx₁,δu₁,...,δx_N,δu_N]^T (19)

问题Γ0的离散形式统一用如公式(20)、(21)所示形式表示：

满足约束：

其中，N_C为总约束数量，C_κ、W_k和l_k的值根据参考轨迹计算。故问题Γ0能够转化成二次规划问题，利用数值方法求解所述离散化的二次规划问题，提高求解效率，进而实现实时在线制导。

步骤四、利用SeDuMi算法求解轨迹最优跟踪制导问题对应的制导律；

最终，根据步骤一至步骤三所建立的轨迹最优跟踪问题模型，利用SeDuMi算法求解跟踪制导律。由于火星高层大气稀薄，探测器所受的气动力较小，控制力极弱，此时控制的效率较低。本文的研究中，假设初始阶段的探测器飞行是无控的，当探测器所受的阻力加速度大于0.3g_e(g_e为地表引力加速度)时，探测器开始实施制导控制，制导周期取为500m/s，初始节点个数N取为400。探测器的初始状态及误差分布如表1所示，参考轨迹为自标称进入点开伞，以常值倾侧角40°积分动力学方程，直到速度达到450m/s时所对应的轨迹。制导指令采用Sedumi求解。具体仿真参数如表1所示：

表1仿真参数

状态	标称值	分布(3σ)	分布类型
				高度	125km	90m	高斯分布
速度	5500m/s	3m/s	高斯分布
				路径角	-15.2°	0.18°	高斯分布
航程	0km	150m	高斯分布

对于末端航程跟踪问题，取系数分别为：P₁＝1，R₁＝1×10^-6。考虑两个案例：案例1，所有初始误差均取最大正值；案例2，所有初始误差均取最小负值，仿真结果如图2至图5所示。图3和图4为根据表1的初始点分布实施500次蒙特卡洛仿真对应的末端航程误差分布。由于阻力小于0.3g_e时探测器无控，开环制导方案与闭环制导方案没有差别，故下图的结果均为从探测器开始实施控制直到降落伞展开。

由图3可以看出，在跟踪飞行器的末端航程时，相对于图4中开环无控的情况，末端航程误差极大减小，最大值不超过1km，而无控飞行时，末端航程误差高达约18km。

图5给出了案例1下的制导律求解时间，随着飞行时间的推进，因为剩余节点数越来越少，制导律求解所需的时间越来越少，制导律求解时间不超过1.2s。

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

Claims (3)

1.火星大气进入段轨迹最优跟踪制导方法，其特征在于：包括如下步骤，

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

步骤二、通过对路径约束和控制约束凸化处理，选取二次型性能指标，建立轨迹最优跟踪问题模型；

步骤三、对步骤二中的轨迹最优跟踪问题模型离散化处理，将其转化成二次规划问题，利用数值方法求解离散化的二次规划问题，提高求解效率，进而实现实时在线制导；

步骤一具体实现方法为，

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

式(1)中，V为探测器速度大小，u＝cosσ为控制量，σ为倾侧角；在无量纲化过程中，长度的量纲单位为火星半径R₀，速度的无量纲单位为

其中

为火表引力加速度，μ为火星引力常数；时间的无量纲单位为

角度的单位为弧度，不需要无量纲化处理；式(1)中，L和D分别为探测器受到的无量纲升力和阻力加速度，分别具有如下形式：

L＝D·L/D (3)

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

其中，ρ₀为参考密度，r₀为参考半径，h_s为标高；

步骤二具体实现方法为，

在小扰动下，状态误差的动力学方程如下：

式中，上标“*”代表参考轨迹与参考控制；A(V)为雅克比矩阵，B_u(V)为式(1)中的动力学方程对控制量的偏导数；

大气进入段跟踪制导需满足路径约束，即动压、过载与热流约束；由于动压与过载成比例关系，此处将动压约束与过载约束合并，则路径约束如下：

其中，

为最小动压约束，q为动压，

为热流，n为过载，下标“max”代表上限约束；

除了步骤一中动力学微分方程与路径约束外，大气进入段跟踪制导还应满足边界条件：

δx(V₀)＝δx₀ (9)

其中δx₀为当前制导周期开始时刻，实际状态与标称状态之间的差值；为了保证线性化的有效性，还需满足信赖域条件：

|δx(V)|≤ε_x (10)

|δu(V)|≤ε_u (11)

除此之外，还应满足状态上下限约束及控制能力约束，即

x_min(V)≤x^*(V)+δx(V)≤x_max(V) (12)

u_min(V)≤u^*(V)+δu(V)≤u_max(V) (13)

为了使用最少的燃料消耗实现轨迹跟踪，在每个制导周期内，取性能指标为：

其中，V₀代表当前制导周期的开始时刻，V_f代表开伞时刻；P和Q(V)为对称半正定矩阵，R(V)为一正数；将轨迹跟踪问题建模成最优控制问题Γ0；

问题Γ0：寻找最优控制δu^*(V)，使得

满足约束：

式(5)、式(8)～式(13)

最优控制问题Γ0即为建立的轨迹最优跟踪问题模型。

2.如权利要求1所述的火星大气进入段轨迹最优跟踪制导方法，其特征在于：步骤三具体实现方法为，

为了求得问题Γ0的数值解，将速度区间[V₀，V_f]均分为N个区间，每个区间的长度为ΔV＝(V_f-V₀)/N，将式(5)、式(8)～式(13)中所示的约束施加于N+1个节点上；用单调递减序列{V₀，V₁，...，V_N}表示自变量节点，其中V_i＝V₀+iΔV，i＝1，...，N，对应的状态和控制序列分别为{δx₀，δx₁，...，δx_N}和{δu₀，δu₁，...，δu_N}；运用梯形积分公式，动力学约束近似如下：

式中，A_x，i＝A_x(x(V_i)，u(V_i))，B_u，i＝B_u(x(V_i)，u(V_i))；上式进一步表示为：

式中，

I为单位矩阵；

将优化变量统一用向量δz表示，即：

δz＝[δx₀，δu₀，δx₁，δu₁，...，δx_N，δu_N]^T (19)

问题Γ0的离散形式统一用如公式(20)、(21)所示形式表示：

满足约束：

其中，δz^*为离散最优控制量，N_C为总约束数量，C_κ、W_k和l_k的值根据参考轨迹计算；故问题Γ0能够转化成二次规划问题，利用数值方法求解所述离散化的二次规划问题，提高求解效率，进而实现火星大气进入段实时在线制导。

3.如权利要求1或2所述的火星大气进入段轨迹最优跟踪制导方法，其特征在于：为了能达到末端航程误差最小，式(14)中的性能指标取为如公式(16)所示形式：

其中，J₁表示末端航程误差最小条件下的性能指标，P₁为一非负数，R₁(V)为一正数。

Images