CN112149225B - 基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法 - Google Patents

Abstract

本发明公开了一种基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法，建立火星探测器定点着陆轨迹优化问题并将其转化为凸优化问题，引入Hermite‑Simpson离散格式对所述凸优化问题进行离散化，推导并建立了基于Hermite‑Simpson离散格式的火星探测器定点着陆轨迹凸优化问题，并采用凸优化算法求解离散的火星探测器定点着陆轨迹凸优化问题。本发明能够在几乎不降低求解效率的情况下，显著提高凸优化方法求解着陆轨迹优化时的精度。

Description

基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法

技术领域

本发明属于飞行动力学领域，特别涉及了一种飞行器定点着陆轨迹凸优化方法。

背景技术

轨迹在线优化技术是航天器高精度着陆的一项重要关键技术。凸优化由于不需要优化初值、计算量小、优化速度快、可靠性高(经过有限次数的迭代即可达到全局最优解，传统优化方法可能需要很多次迭代才能接近最优解)等优点，近年来成为轨迹优化乃至最优控制领域的一个研究热点。传统的轨迹凸优化方法基于梯形离散格式，存在精度偏低的弊端，精度低带来的直接问题是离散最优轨迹与实际积分得到的轨迹差异明显，导致实际轨迹的终端约束(着陆位置、着陆速度等)和路径约束不能严格满足要求。近期提出的高阶离散格式(比如伪谱法离散格式)虽然精度较高，但是计算量较大，降低了凸优化技术的实时性。

发明内容

为了解决上述背景技术提到的技术问题，本发明提出了基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法，提高凸优化方法求解着陆轨迹优化时的精度。

为了实现上述技术目的，本发明的技术方案为：

基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法，建立火星探测器定点着陆轨迹优化问题并将其转化为凸优化问题，引入Hermite-Simpson离散格式对所述凸优化问题进行离散化，推导并建立基于Hermite-Simpson离散格式的火星探测器定点着陆轨迹凸优化问题，并采用凸优化算法求解离散的火星探测器定点着陆轨迹凸优化问题。

进一步地，所述火星探测器定点着陆轨迹优化问题的目标函数如下：

min J＝m(t_f)

着陆器的状态方程如下：

其中，J为目标函数，r为着陆器的位置矢量，v为速度矢量，字母上方的点表示一阶导数，g为火星重力加速度矢量，m为着陆器的质量，t_f为着陆时间，m(t_f)表示t_f时刻的m取值；参数I_sp为燃料比冲，g_e为地球海平面处的重力加速度，φ为推力器的安装角；T_c为推力器的推力，T_c＝nT_maxT_rcosφ·e，T_max为单个推力器的最大推力，n为推力器的个数，e为沿总推力方向的单位矢量，T_r为每个推力器的相对推力；

火星探测器在着陆过程中需要满足的视线角约束如下：

其中，r_x(t),r_y(t),r_z(t)为着陆器的位置向量r的3轴分量，t为时间，参数为着陆器的视线角下限。

进一步地，将火星探测器定点着陆轨迹优化问题转化为凸优化问题，定义以下3个变量：

z＝log(m)

其中，Γ为T_c的模；

目标函数如下：

其中，括号中的t表示对应参数在t时刻的取值，t₀为初始时间；

状态方程如下：

三个路径约束如下：

||u_c(t)||₂≤σ(t)

其中，ρ_l和ρ_u分别为Γ的下边界和上边界，z_l和z_u分别为z的下边界和上边界；上述第三个路径约束称为视线角约束。

进一步地，引入Hermite-Simpson离散格式对所述凸优化问题进行离散化的过程如下：

对于如下线性状态方程：

其中，x为状态变量，u为控制变量，c为常数向量，A为系数矩阵，B为系数矩阵；定义如下变量：

其中，x_i为第i个离散点的状态变量，u_i为第i个离散点的控制变量，为区间中点处的u_i取值，N为离散点数；

其中，I_n为单位矩阵，h_i＝t_i+1-t_i，t_i为第i个离散点对应的时刻，Δt＝t_f-t₀；定义如下向量：

其中，

则状态方程的离散形式如下：

A_dynx_dyn＝b_dyn

对于如下线性形式的目标函数：

其中，x(t_f)为t_f时刻的状态变量，c_M为目标函数端点项系数，c_x和c_x分别为目标函数积分项的系数，上标T表示转置；该目标函数的离散形式如下：

其中，c_M为目标函数端点项系数，上标T表示转置；D₁＝[1,…,1]_1×N；

其中：

将第一个路径约束进行离散并写为二阶锥约束的标准形式：

其中，u_x,u_y,u_z为u的3轴分量；在区间中点处，该约束的离散形式如下：

其中，为u_x,u_y,u_z,σ在区间中点处的值；

第二个路径约束的左端约束的离散形式如下：

其中，b_2E＝[b₂；-1],w_i＝[x_i；σ_i]；

其中，z_l,i为z_l在第i个离散点的值，σ_i为第i个离散点的σ值；在区间中点处，该约束的离散形式如下：

其中，

其中，为z_l,i在区间中点处的值；

第二个路径约束的右端约束的离散形式如下：

其中，

c′_2E＝[c′₂；-1]

d₂＝ρ_uexp(-z_u,i)·(1+z_u,i)

c′₂＝[0 0 0 0 0 0 -ρ_uexp(-z_u,i)]^T

其中，z_u,i为z_u在第i个离散点处的值；在区间中点处，该约束的离散形式如下：

其中，

其中，为z_u,i在区间中点处的值；

第三个路径约束(视线角约束)的离散形式的如下：

其中，r_x,r_y,r_z为位置向量r的3轴离散分量；在区间中点处，该约束的离散形式如下：

其中，

采用上述技术方案带来的有益效果：

本发明为轨迹凸优化技术引入一种基于Hermite-Simpson格式的高阶离散格式，推导并建立基于Hermite-Simpson格式的着陆轨迹凸优化方法，能够显著提高凸优化方法求解飞行器着陆轨迹优化问题的精度(相对传统的凸优化离散格式，在计算量几乎不增加的情况下，能够将轨迹优化误差降低约5个数量级)。

附图说明

图1为实施例中火星探测器最优推力加速度随时间变化曲线图；

图2为实施例中火星探测器最优推力大小随时间变化曲线图；

图3为实施例中火星探测器最优着陆轨迹的位置分量随时间变化曲线图；

图4为实施例中火星探测器最优着陆轨迹的速度分量随时间变化曲线图；

图5为实施例中火星探测器最优着陆轨迹的质量随时间变化曲线图；

图6为实施例中火星探测器最优着陆轨迹的高度-航程曲线图；

图7为实施例中优化的目标函数与Radau伪谱凸优化方法的对比图；

图8为实施例中计算耗时与Radau伪谱凸优化方法的对比图；

图9为实施例中最大位置误差与Radau伪谱凸优化方法的对比图；

图10为实施例中最大速度误差与Radau伪谱凸优化方法的对比图；

图11为实施例中初始参数不确定情况下火星探测器最优推力加速度随时间变化曲线图；

图12为实施例中初始参数不确定情况下火星探测器最优推力大小随时间变化曲线图；

图13为实施例中初始参数不确定情况下火星探测器最优着陆轨迹图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

本发明设计了基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法，建立火星探测器定点着陆轨迹优化问题并将其转化为凸优化问题，引入Hermite-Simpson离散格式对所述凸优化问题进行离散化，建立基于Hermite-Simpson离散格式的火星探测器定点着陆轨迹凸优化问题，并采用凸优化算法求解离散的火星探测器定点着陆轨迹凸优化问题。

1.火星探测器着陆轨迹凸优化问题描述

火星探测器着陆段轨迹优化问题的目标函数为

min J＝m(t_f) (1)

状态方程为

其中，r∈R³为着陆器的位置矢量，v∈R³为速度矢量，二者均在火星地面固连坐标系中描述。火星重力加速度矢量为g＝[00-3.7114]^Tm/s²。变量m为着陆器的质量，m＝1905kg。着陆时间t_f为给定值，t_f＝81s。参数α定义如下

参数I_sp＝225s为燃料比冲，g_e＝9.807m/s2为地球海平面处的重力加速度，φ＝27deg为推力器的安装角。推力器的推力表达式如下

T_c＝nT_maxT_rcosφ·e (4)

其中，T_max＝3.1kN为单个推力器的最大推力，n＝6为推力器个数，e为沿推力方向的单位矢量。每个推力器的相对推力T_r＝||T_c||T_max需要满足如下约束

T_l≤T_r≤T_u (5)

其中参数T_l＝0.3，T_u＝0.8。着陆器的位置和速度边界条件为

其中初始时间t₀＝0。探测器在着陆过程中还需要满足另一个路径约束(即视线角约束)

其中变量r_x(t),r_y(t),r_z(t)为着陆器的位置向量r的3个分量，t为时间，参数为着陆器下降过程中允许的最小高低角，/>

以二阶锥凸优化为例，目标函数为

约束条件为

其中：x∈Rⁿ为优化变量，c₀∈Rⁿ为目标函数系数，A₀∈R^m×n和b₀∈R^m×1为描述等式约束的常数矩阵，c_i∈R^n×1和d_i∈R为描述锥约束的参数。

为了能够采用二阶锥凸优化的方法求解火星探测器着陆轨迹优化问题，首先需要将着陆轨迹优化问题转化为凸优化问题。为此，定义如下变量

其中：Γ为为推力矢量的模，即Γ＝||T_c||，满足如下约束

ρ_l≤Γ≤ρ_u (11)

其中：参数ρ_l＝0.3T_max,ρ_u＝0.8T_max。利用上述变量，目标函数可写为

路径约束为

||u_c(t)||₂≤σ(t) (14)

ρ_le^-z(t)≤σ(t)≤ρ_ue^-z(t) (15)

方程(15)描述的约束不是凸约束。为了将其转化为凸约束，将方程(15)的左端在z_l处展开，保留二阶项，右端在z_u处展开，保留一阶项，得到

其中：z_l和z_u分别为变量z的下边界和上边界，其离散值为

z_l,i＝log(m₀-αρ_lt_i),i＝0,…,N

z_u,i＝log(m₀-αρ_ut_i),i＝0,…,N

其中：m₀为着陆器的初始质量，t_i为时间的离散值，N表示离散节点个数。

视线角约束(7)经过变形可以写成如下凸优化形式

由方程(12)～(14)、(16)和(17)描述的问题为着陆轨迹凸优化问题。

2.传统离散格式与高精度离散格式

采用凸优化求解轨迹优化问题时，需要将连续轨迹优化问题离散化，在离散点上进行优化。轨迹优化问题本质上属于最优控制问题。不失一般性，这里以Bolza型最优控制问题为例，介绍轨迹优化的离散格式。Bolza型最优控制问题描述为：求解最优控制u(t)＝u^*(t)，使得如下Bolza型目标函数最小化

并且满足状态方程

端点约束

Φ[x(t₀),t₀,x(t_f),t_f]＝0, (20)

以及路径约束

C[x(t),u(t),t]≤0,t∈[t₀,t_f], (21)

其中，M表示目标函数中的Mayer函数，L表示目标函数中的Lagrange函数，Φ表示端点约束函数，C表示路径约束函数。

目前，传统的轨迹凸优化方法采用的离散格式是梯形格式(二阶格式)。采用梯形格式将方程(18)-(21)所描述的轨迹优化问题在离散节点上进行离散，得到的优化变量为(x₀,x₁,···,x_N；u₀,u₁,···,u_N；t₀,t_f)，目标函数为

约束条件为

C_i＝C(x_i,u_i,t_i)≤0,(i＝0,1,···,N) (24)

Φ(x₀,t₀,x_f,t_f)＝0 (25)

其中，Δt＝t_f-t₀,h_i＝t_i+1-t_i，x_i＝x(t_i)，L_i＝L(x_i,u_i,t_i)，f_i＝f(x_i,u_i,t_i)。

梯形离散格式虽然形式简单，但是精度较低，需要较多节点才能提高精度。为了克服该缺点，本专利引入Hemite-Simpson格式(四阶格式)。Hemite-Simpson格式需要用到区间中点的变量和函数值，为此需要将区间中点的控制变量作为优化变量，并且通常在区间中点处添加路径约束，区间中点的集合为

采用Hemite-Simpson格式得到的NLP的优化变量为(x₀,x₁,···,x_N；u₀,u₁,···,u_N；t₀,t_f),目标函数为

约束条件为

/>

C_i＝C(x_i,u_i,τ_i；t₀,t_f)≤0,(i＝0,1,···,N) (29)

Φ(x₀,t₀,x_f,t_f)＝0 (31)

其中

3.基于Hemite-Simpson离散格式的着陆轨迹凸优化方法

为了应用凸优化算法求解离散后的轨迹优化问题，需要将目标函数和约束全部改写为凸优化的标准形式。以状态方程离散格式为例，考虑到凸优化问题的状态方程为线性形式的状态方程，对于如下一般形式的线性状态方程

其中，x∈R^n×1为状态变量，u∈R^m×1为控制变量，c∈R^n×1为常数向量，A∈R^n×n为系数矩阵，B∈R^n×m为系数矩阵。将方程(32)代入到方程(28)得到

其中，I_n为n×n的单位矩阵。为了简化方程(33)的形式，定义如下变量

/>

其中，对于i＝0,1,···,N-1，上式中的各项的定义如下

定义如下向量

其中

那么，方程(33)可以写出如下标准形式

A_dynx_dyn＝b_dyn (37)

将目标函数写成一般线性形式

其中，c_M∈R^n×1为目标函数端点项系数，c_x∈R^n×1和c_x∈R^n×1分别为目标函数积分项的系数，u∈R^m×1为控制向量。将方程(38)代入方程(27)可得到

为了简化方程(39)的形式，定义如下矩阵

其中，对于i＝0,1,···,N-1，各项定义如下

应用上述定义，目标函数可以写为

其中

D₁＝[1,…,1]_1×N

对于将路径约束(14)，在节点上进行离散然后写为二阶锥约束的标准形式

其中u_x,u_y,u_z为控制变量u_c的3个分量。在区间中点处，该约束为

其中表示控制变量u_x,u_y,u_z,σ在区间中点中间节点/>处的取值，为优化变量。对于路径约束(16)的左端约束，将其离散化得到

将其改写成另一种更普遍的形式得到

其中

对方程(45)进行扩展得到

其中

方程(46)可写成如下锥约束的标准形式

在区间点中间处，该约束的离散形式为

其中

根据Hemite-Simpson离散格式可知

整理得到

其中

将方程(50)代入方程(48)得到

采用如下记法对系数矩阵进行扩展

那么，方程(51)可写成如下形式

方程(52)可进一步写成如下二阶锥约束的标准形式

对于路径约束(16)的右端项，其离散格式为

σ_i≤ρ_uexp(-z_u,i)[1-(z_i-z_u,i)],i＝0,···,N (54)

整理得到

σ_i≤c′₂ ^Tx_i+d₂,i＝0,···,N (55)

式中

c′₂＝[0 0 0 0 0 0 -ρ_uexp(-z_u,i)]^T

d₂＝ρ_uexp(-z_u,i)·(1+z_u,i)

方程(55)可进一步改写为

式中

c′_2E＝[c′₂；-1]

在区间点中间处，该约束的离散形式为

其中和/>分别表示变量z_i+1和z_u在节点/>处的取值。整理变形得到

式中

将代入方程(58)得到

记

那么上式可改写为

对于视线角约束的离散形式可写为如下二阶锥约束的标准形式

在区间点中间处，该约束为

其中：表示/>在中间节点/>处取值。

将其改写为向量形式，得到

式中

将代入方程(63)得到如下二阶锥约束的标准形式

由方程(37)、(39)、(42)、(43)、(47)、(53)、(56)、(60)、(61)和(64)描述的凸优化问题为采用Hemite-Simpson格式离散得到的火星探测器定点着陆轨迹凸优化问题，可采用一些高效的凸优化方法(比如内点法)快速求解。

本实施例采用基于内点法的求解器ECOS求解该问题，仿真结果如图1-6所示。图1和图2给出优化出的最优控制变量(加速度分量和推力大小)随时间变化曲线。其中，“○”为节点处的离散最优解，“*”为各区间中点处的离散最优解，细实线为根据离散最优解插值得到的连续最优控制变量。图3为最优轨迹的位置分量随时间变化曲线，图4为最优轨迹的速度分量随时间变化曲线，图5为最优轨迹对应的着陆质量随时间变化曲线。其中，“○”为节点处的离散最优解，细实线为根据图1和图2所示的最优控制变量采用四阶龙格-库塔数值积分方法对状态方程积分得到的结果。可见，数值积分结果与离散最优解非常一致。图6为最优轨迹曲线。可见，着陆器着陆过程中准确满足视线角要求。

表1给出本发明方法的离散最优解和数值积分结果在轨迹终端的差异情况以及计算效率情况。为了对比，表中还给出了传统离散格式(即梯形离散格式)的结果。其中传统离散格式的结果来源于航天器动力学与控制领域的顶级期刊Journal of GuidanceControl&Dynamics(Journal of Guidance Control&Dynamics,2018,41(2):1-15)。由表1可知，本发明方法的离散最优解与数值积分结果非常一致，二者在轨迹终端的误差非常小。与传统离散格式相比，本发明方法将终端位置误差和速度误差降低了大约5个数量级，优化的最优燃料消耗量减少0.87kg，优化耗时有轻微增加。可见，本发明方法具有很高的精度和效率。

表1本发明方法与传统凸优化方法的优化结果对比(N＝40)

图7给出在采用不同离散节点数目的情况下，本发明方法优化的最优燃耗与Radau伪谱凸优化方法的结果对比。可见，本发明方法优化的最优燃耗与Radau伪谱法的结果比较接近，但是比Radau伪谱法优化的燃耗略小。图8给出对于不同离散节点数，本发明方法的计算耗时与Radau伪谱法的结果对比。可见，本发明方法的计算耗时比Radau伪谱法低大约1个数量级。图9和图10分别给出对于不同离散节点数，本发明方法的最大位置误差和最大速度误差与Radau伪谱法的结果对比。可见，本发明方法将误差减小了4～5个数量级。

为了测试本专利方法的鲁棒性，假设着陆轨迹的初始状态参数和着陆时间的不确定性范围及其分布规律如表2所示。采用蒙特卡洛方法对着陆轨迹进行优化仿真，图11～13给出25条着陆轨迹的仿真结果。其中图11为最优推力加速度曲线，图12为最优推力大小曲线，图13为最优着陆轨迹曲线。可见，存在扰动时，本发明方法仍然能够可靠低优化出最优轨迹，反应了方法的鲁棒性。

表2着陆轨迹参数不确定性

参数	不确定范围	分布规律	单位
				||r(t0)||2	±25	均匀分布	m
||v(t0)||2	±10	均匀分布	m/s
				m(t0)	±1	均匀分布	kg
tf	±3	均匀分布	s

实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (3)

1.基于高精度离散格式的火星探测器定点着陆轨迹凸优化方法，其特征在于，建立火星探测器定点着陆轨迹优化问题并将其转化为凸优化问题，引入Hermite-Simpson离散格式对所述凸优化问题进行离散化，推导并建立基于Hermite-Simpson离散格式的火星探测器定点着陆轨迹凸优化问题，并采用凸优化算法求解离散的火星探测器定点着陆轨迹凸优化问题；

所述火星探测器定点着陆轨迹优化问题的目标函数如下：

min J＝m(t_f)

着陆器的状态方程如下：

其中，J为目标函数，r为着陆器的位置矢量，v为速度矢量，字母上方的点表示一阶导数，g为火星重力加速度矢量，m为着陆器的质量，t_f为着陆时间，m(t_f)表示t_f时刻的m取值；参数I_sp为燃料比冲，g_e为地球海平面处的重力加速度，φ为推力器的安装角；T_c为推力器的推力，T_c＝nT_maxT_rcosφ·e，T_max为单个推力器的最大推力，n为推力器的个数，e为沿总推力方向的单位矢量，T_r为每个推力器的相对推力；

火星探测器在着陆过程中需要满足的视线角约束如下：

其中，r_x(t),r_y(t),r_z(t)为着陆器的位置向量r的3轴分量，t为时间，参数为着陆器的视线角下限。

2.根据权利要求1所述火星探测器定点着陆轨迹凸优化方法，其特征在于，将火星探测器定点着陆轨迹优化问题转化为凸优化问题，定义以下3个变量：

z＝log(m)

其中，Γ为T_c的模；

目标函数如下：

其中，括号中的t表示对应参数在t时刻的取值，t₀为初始时间；

状态方程如下：

第一个路径约束如下：

||u_c(t)||₂≤σ(t)

第二个路径约束如下：

第三个路径约束如下：

其中，ρ_l和ρ_u分别为Γ的下边界和上边界，z_l和z_u分别为z的下边界和上边界；上述第三个路径约束称为视线角约束。

3.根据权利要求2所述火星探测器定点着陆轨迹凸优化方法，其特征在于，引入Hermite-Simpson离散格式对所述凸优化问题进行离散化的过程如下：

对于如下线性状态方程：

其中，x为状态变量，u为控制变量，c为常数向量，A为系数矩阵，B为系数矩阵；定义如下变量：

其中，x_i为第i个离散点的状态变量，u_i为第i个离散点的控制变量，为区间中点处的u_i取值，N为离散点数；

其中，I_n为单位矩阵，h_i＝t_i+1-t_i，t_i为第i个离散点对应的时刻，Δt＝t_f-t₀；定义如下向量：

其中，

则状态方程的离散形式如下：

A_dynx_dyn＝b_dyn

对于如下线性形式的目标函数：

其中，x(t_f)为t_f时刻的状态变量，c_M为目标函数端点项系数，c_x和c_x分别为目标函数积分项的系数，上标T表示转置；该目标函数的离散格式如下：

其中，D₁＝[1,…,1]_1×N；

其中：

将第一个路径约束进行离散并写为二阶锥约束的标准形式：

其中，u_x,u_y,u_z为u的3轴分量；在区间中点处，该约束的离散形式如下：

其中，为u_x,u_y,u_z,σ在区间中点处的值；

第二个路径约束的左端约束的离散形式如下：

其中，b_2E＝[b₂；-1],w_i＝[x_i；σ_i]；

其中，z_l,i为z_l在第i个离散点的值，σ_i为第i个离散点的σ值；在区间中点处，该约束的离散形式如下：

其中，

其中，为z_l,i在区间中点处的值；

第二个路径约束的右端约束的离散形式如下：

其中，

c′_2E＝[c′₂；-1]

d₂＝ρ_uexp(-z_u,i)·(1+z_u,i)

c′₂＝[0 0 0 0 0 0 -ρ_uexp(-z_u,i)]^T

其中，z_u,i为z_u在第i个离散点处的值；在区间中点处，该约束的离散形式如下：

其中，

其中，为z_u,i在区间中点处的值；

第三个路径约束的离散形式的如下：

其中，r_x,r_y,r_z为位置向量r的3轴离散分量；在区间中点处，该约束的离散形式如下：

其中，