下载:
-
行星探测是人类探索太空最重要的途径之一，月球和火星探测是当前行星探测的热点。当目标行星大气比较稀薄或者没有大气时，探测器将无法使用降落伞降至安全速度，基于反推火箭的动力下降成为实现行星表面软着陆的一种重要手段。由于行星动力下降过程存在诸多不确定性因素，其在线轨迹优化是实现探测器安全、精确着陆所需要解决的重要问题。
近年来，序列凸优化方法成为求解上述在线轨迹优化的有效手段[1-2]，其本质是将无限维空间中的最优控制问题通过离散转化为有限维空间中的非线性规划问题，并利用序列凸优化理论通过凸化处理后迭代求解。Szmuk等[3]通过引入虚拟控制量并且在动力学方程中定义了新的时间状态量，提出了一种终端时间自由的序列凸优化方法计算框架，能够较好地得到最优燃料轨迹；其不足之处是对初始猜想敏感且迭代次数多。Blackmore等[4]通过将燃料最优问题凸化为二阶锥规划（Second-Order-Cone-Programming，SOCP）问题，并考虑了坡度约束，能够较快求解出燃料最优的下降轨迹。上述方法对于终端时间固定的情况能够得到较好的解，但由于其终端时间采用线性搜索方法计算，效率较低。此外，Liu等[5]考虑了气动力作用下的一级火箭回收问题，通过无损凸化技术和动力学方程归一化提高了算法收敛性和快速性。该方法在凸化过程中重定义了控制量，公式推导复杂且引入了额外的约束和计算量，只能处理固定终端时间问题。
有别于上述研究方法，本文探索使用预标记中心差分替代解析表达式计算动力学矩阵，并对动力学方程采用直接线性化替代变量代换后再线性化的凸化方法，最后研究行星动力下降燃料最优控制问题中近似最优终端时间的确定方法，以提高序列凸优化在线轨迹优化方法的工程实用性。
-
以着陆点为原点建立空间直角坐标系，着陆点一般不直接选在地表，而是距离地面稍有一定距离，方便提供安全缓冲余量。x轴方向预先约定，z轴垂直于地表向上，y轴满足右手螺旋关系。行星动力下降可以表述为经典的Bolza型最优控制问题
$$ \left\{ \begin{array}{l} \min {{\theta}} ({{{x}}_{{t_{\rm{f}}}}}) + \int_{{t_0}}^{{{\rm{t}}_f}} {\dfrac{{\left\| {{F}} \right\|}}{{{g_0}{I_{{\rm{sp}}}}}}{\rm{d}}t} \\ {\rm{s}}{\rm{.t}}{\rm{. }}\;\;\;\;{{{x}}_0} = {{x}}({t_0}) \\ \;\;\;\;\;\;\;\;\; \dot {{x}} = {{f}}({{x}},{{u}},t) \\ \;\;\;\;\;\;\;\;\;{F_{\min }} \text{≤} \left\| {{F}} \right\| \text{≤} {F_{\max }} \\ \;\;\;\;\;\;\;\;\;{{{r}}_z} \text{≥} 0 \\ \end{array} \right. $$ (1) 其中终端性能指标为
$$\theta ({{{x}}_{{t_{\rm{f}}}}}) = \left| {{ w^{\rm{T}}}\left( {{ x_{{t_{\rm{f}}}}} - x({t_{\rm{f}}})} \right)} \right|$$ (2) 式（1）中约束分别为初值约束、动力学约束、控制量幅值约束和防撞约束，其中反推发动机为了避免多次开关机带来的可靠性问题，在实际工作中一般不将其完全关闭，故需限定最小推力。行星动力下降动力学方程为
$$\left\{ \begin{array}{l} {\dot{{ r}}} = {{ v}} \\ \dot {{v}} = {{g}} + \dfrac{{{F}}}{m} + {{d}} \\ \dot m = - \dfrac{{\left\| {{F}} \right\|}}{{{g_0}{I_{{\rm{sp}}}}}} \end{array} \right.$$ (3) 其中：g= [0 0 –g（rz）]T是行星的重力加速度矢量。
考虑到动力下降段高度在2~5 km之间，高度相比行星半径尺度较小，重力加速度变化也较小，可以直接通过海拔高度计算。本文以火星引力场为例，其表面重力加速度为gM= 0.377、gE= 3.7 m/s2、RM= 3 397 km，距地表任意高度的重力加速度为
$$g({r_z}) = {g_{\rm{M}}}{\left( {\frac{1}{{1 + {r_z}/{R_{\rm{M}}}}}} \right)^2}$$ (4) 此外，d= [dxdydz]T是3个方向的气动加速度，分为可建模的气动（如飞行器的气动阻力）和无法建模的随机风场扰动。由于火星大气较为复杂，不适合在在线轨迹优化中应用[1]，因此本文采用文献[3~4]的处理方式，在动力下降段中以忽略气动力作用下简化的动力学进行分析，气动扰动可以通过后期设计控制器加以克服[6]。
-
采用直接法求解优化问题式（1）时，需要先将其离散到有限维空间中。当终端时间tf确定时，将飞行时间[0，tf]等距离散为N个时刻，即N–1个子区间，并在这N个时刻对状态量和控制量进行离散，方便施加过程约束。
对于控制力幅值约束
$${F_{\min }} \text{≤} \sqrt {F_x^2 + F_y^2 + F_z^2} \text{≤} {F_{\max }}$$ (5) 因最小推力幅值约束是二阶锥补集，为将非凸约束式（5）凸化，在各个离散点处引入推力松弛因子Fη
$$ \sqrt{{F}_{xi}^{2}+{F}_{yi}^{2}+{F}_{zi}^{2}}\text{≤} {F}_{\eta i},i=0,1,\cdots,N-1$$ (6) 此时，控制幅值约束可以表示为
$$ {F}_{\mathrm{min}}\text{≤} {F}_{\eta i}\text{≤} {F}_{\mathrm{max}},i=0,1,\cdots,N-1$$ (7) 至此非凸约束式（5）被等价转换为凸约束式（6）和式（7），积分型性能指标也可以表达为关于Fη的线性组合
$$\int_{{t_0}}^{{t_f}} {\frac{{\left\| {{ F}} \right\|}}{{{g_0}{I_{{\rm{sp}}}}}}{\rm{d}}t} \approx \frac{{\Delta t}}{{{g_0}{I_{{\rm{sp}}}}}}\sum\limits_{i = 0}^{N - 1} {{F_{\eta i}}} $$ (8) 其中：Δt=（tf–t0）/（N–1），即使用矩形公式将积分化为求和。
-
为将非线性动力学微分方程约束处理为凸约束，需进行凸化处理。文献[5]和文献[7]中处理上述动力学方程通常使用2种凸化技巧：将推力矢量分解为推力大小和方向余弦、对质量取对数重定义新的状态量。前者存在的问题是，方向余弦约束为非凸等式约束，虽然文献[3]中通过极小值原理证明其可无损转化为凸锥约束，但实际数值求解中经常出现不稳定因素，造成收敛速度较慢，且控制变量从原来的3维上升到4维，问题复杂度上升。后者对质量取对数后重定义新的状态量，再对推力约束进行一阶泰勒展开线性化和直接对动力学泰勒展开线性化，本质上都是线性化，并未降低原问题的非线性，仅将非线性从动力学方程转移到过程约束中。取对数后推力不等式约束中存在指数函数，相比原来仅含加、乘法的公式增加了计算量，且人为增加了推导公式和编程的工作量。
本文中不再对动力学使用多余的凸化技巧，直接在第k+1次迭代时于各个离散时间点处，使用第k次优化得到的状态量和控制量的解对非线性动力学方程进行线性化，使其在每个子区间[nΔt，（n+1）Δt]（n= 0，1，···，N–2）内变为线性动力学，即
$\dot {{x}} = {{A}}_n^{(k)}{{x}} + $ ${{B}}_n^{(k)}{{u}} + {{C}}_n^{(k)}$ 。对于动力学矩阵A、B、C的分量使用差分直接计算，第i个微分方程，对第n个离散点处的状态量（控制量）的第j个分量的偏导的差分计算公式为
$$\left\{ \begin{array}{*{20}{l}} {{{A}}_{ij}}({{x}}_n^{(k)},{{u}}_n^{(k)}) = \dfrac{{\partial {f_i}}}{{\partial {{{x}}_j}}} =\\ \dfrac{{{f_i}({{x}}_n^{(k)} + {{{e}}_j} \times {{\delta}},{{u}}_n^{(k)}) - {f_i}({{x}}_n^{(k)} - {{{e}}_j} \times {{\delta}},{{u}}_n^{(k)})}}{{2{{\delta}} }} \\ {{{B}}_{ij}}({{x}}_n^{(k)},{{u}}_n^{(k)}) = \dfrac{{\partial {f_i}}}{{\partial {{{u}}_j}}} =\\ \dfrac{{{f_i}({{x}}_n^{(k)},{{u}}_n^{(k)} + {{{e}}_j} \times {{\delta}} ) - {f_i}({{x}}_n^{(k)},{{u}}_n^{(k)} - {{{e}}_j} \times {{\delta}} )}}{{2{{\delta}} }} \\ {{C}}({{x}}_n^{(k)},{{u}}_n^{(k)}) = {{f}} - {{{{f}}'}_x}{{x}} - {{{{f}}'}_u}{{u}} = {{f}}({{x}}_n^{(k)},\\ {{u}}_n^{(k)}) - {{A}}({{x}}_n^{(k)},{{u}}_n^{(k)}){{x}}_n^{(k)} - {{B}}({{x}}_n^{(k)},{{u}}_n^{(k)}){{u}}_n^{(k)} \end{array} \right.$$ (9) 为了消除不必要的计算量，在设计程序时，提前对每个动力学微分方程显含的状态量和控制量索引进行标注，仅计算这些位置的偏导数，其余元素直接赋值为0。将式（1）中的微分方程约束转换为线性代数方程约束，为表示方便，简记为
$$\left\{ \begin{array}{l}{{{A}}}_{n}^{(k)}={{A}}({{x}}_{n}^{(k)},{{{u}}}_{n}^{(k)}),\\ {{{B}}}_{n}^{(k)}={{B}}({{x}}_{n}^{(k)},{{{u}}}_{n}^{(k)}),\\ {{{C}}}_{n}^{(k)}={{C}}({{x}}_{n}^{(k)},{{{u}}}_{n}^{(k)})\end{array} \right.$$ (10) 为兼顾计算精度和计算效率，使用梯形法施加动力学约束
$$\left\{ \begin{array}{l} {{{k}}_1} = {{A}}_n^{(k)}{{{x}}_n} + {{B}}_n^{(k)}{{{u}}_n} + {{C}}_n^{(k)} \\ {{{k}}_2} = {{A}}_{n + 1}^{(k)}{{{x}}_{n + 1}} + {{B}}_{n + 1}^{(k)}{{{u}}_{n + 1}} + {{C}}_{n + 1}^{(k)} \\ {{{x}}_{n + 1}} = {{{x}}_n} + \dfrac{{\Delta t}}{2}({{{k}}_1} + {{{k}}_2}) \\ \end{array}\!\!\!, \right.\;\;n = 0,2, \cdots,N - 2$$ (11) 启动第一次迭代时，动力学矩阵A，B，C的计算需要提供k= 0时每个离散点n处的状态量x（0）和控制量u（0），即初始猜想。对于状态量，如果同时存在初值约束和终端约束，则状态量的初始猜想直接设计为两个状态的线性插值
$${{x}}_n^{(0)} = \left( {1 - \frac{n}{N}} \right){{{x}}_0} + \frac{n}{N}{{{x}}_{{t_{\rm f}}}},n = 0,1, \cdots,N - 1$$ (12) 如果仅存在初值约束（如质量m），则状态量初始猜想设计为常值向量
$${{x}}_n^{(0)} = {{{x}}_0},n = 0,1, \cdots,N - 1$$ (13) 对于控制量，初始猜测全部初始化为0，即
$${{u}}_n^{(0)} = 0,n = 0,1, \cdots,N - 1$$ (14) 式（12）~（14）简洁的设计过程直接消除了利用数值积分设计的轨迹引入的人为因素和计算量，虽然初始猜想并不符合动力学方程约束，但在每次迭代求解后得到的轨迹将完全满足动力学约束，并在不断迭代中使轨迹收敛到最优解。
-
为验证前文所设计算法的有效性，给定飞行状态和航天器参数如表1所示，求解行星动力下降段燃料最优的轨迹优化问题。
表 1仿真设定参数
Table 1.Simulation parameters
参数 数值 初始坐标r0/m （0，0，10 000） 初始速度v0/（m·s–1） （400，–200，–200） 初始质量m0/kg 50 推力幅值[Fmin，Fmax]/N [100，1 000] 终端坐标rf/m （0，0，0） 终端速度vf/（m·s–1） （0，0，0） 发动机比冲g0Isp/（m·s–1） 2 000 终端时间tf/s 60 离散点个数N 100 设置迭代次数为5次，每步迭代的性能指标J，性能指标绝对偏差e= |Jk+1–Jk|和相对偏差er=e/Jk+1，仿真结果如图1所示。一般地，序列凸优化常以
$ \Vert {{{x}}}_{n}^{(k+1)}-{{{x}}}_{n}^{(k)}\Vert \text{≤} \varepsilon,$ $n=0,1,\cdots,N-1 $ 作为迭代收敛的条件[8]，但在轨迹震荡时易导致无法快速收敛且收敛指标$\varepsilon $ 需根据不同物理量进行独立设计，不具有普适性。分析图1可知，er作为判断收敛指标效果最好，也在一定程度上克服了上述缺点。
图 1性能指标及其绝对偏差和相对偏差随迭代次数变化图
Figure 1.Performance index and its absolute deviation，relative deviation versus iterations
本文采用C++调用SOCP求解器ECOS[9]进行求解。计算环境建立在CPU为i5-8300H（基频2.2 GHz，睿频3.8 GHz）Windows系统的PC上，耗时统计从进入main函数开始到满足序列凸优化迭代条件er＜ 1%退出为止（不含优化结果输出时间），求解500次，均为迭代3步收敛，耗时统计如图2所示，平均耗时0.130 1 s，优化所得控制量与状态量如图3~6所示。
优化得到的推力曲线如图3所示，可以看到推力幅值被精确地约束在[100，1 000]N的区间内，表现为bang-bang控制。速度曲线如图4所示，可看到终端速度被精确约束在0 m/s，图5中质量曲线m（t）的导数是控制力幅值函数||F（t）||2，可以看到m（t）随||F（t）||2也相应分为3段。
图 3发动机3个方向推力曲线
Figure 3.Thrust component in three directions
图 4飞行器速度曲线
Figure 4.Velocity profiles of optimal trajectory
图1证明了用中心差分法建立序列凸优化模型的有效性，良好收敛性和对初始猜测弱敏感性，事实上由图1还可以看出如果以相对偏差作为性能指标，可认为第3步已经收敛，轨迹随迭代次数变化如图7所示，3~5次迭代产生的轨迹基本重合，算例对应的三维轨迹和各离散点处推力矢量如图8所示。
图 5飞行器质量曲线
Figure 5.Mass profile of optimal trajectory.
图 6飞行器位置相对于着陆点变化曲线
Figure 6.Position profiles of optimal trajectory
图 7tf= 60 s轨迹随迭代次数变化
Figure 7.Trajectory versus iterations for case of tf=60 s
图 8tf= 60 s优化轨迹和推力矢量
Figure 8.Optimal trajectory and thrust for case oftf= 60 s
-
如果选取不同的终端时间tf，可以求解出对应的最优控制轨迹如图9所示，并同时可以得到飞行器剩余质量如图10所示，可见存在最优终端时间
$t_{\rm{f}}^*$ 使得在所有最优轨迹中消耗燃料最少，在mf-tf散点图中，可发现右端具有明显的线性递减的特征，左端递增且增速快于右端，极值点必然存在。经过大量拟合测试，发现采用次数分别为–3，0，1的组合多项式${m_{\rm{f}}}({t_{\rm{f}}}) = a{t_{\rm{f}}} + b/t_{\rm{f}}^3 + c$ 进行拟合后，均方误差和确定系数均较理想，且兼顾了极值点计算的简易性。上述问题中拟合参数a= –0.045 04，b= –6.72 × 105，c= 37.1，拟合结果的标准差为0.042 2，确定系数为0.999 1。通过对该函数求极值，可以直接得到其近似最优终端时间
图 9tf∈[50，190]s最优轨迹簇
Figure 9.Optimal trajectories attf∈[50，190]s
图 10不同tf下最优轨迹对应的剩余燃料质量mf
Figure 10.mfof optimal trajectories at differenttf
$$t_{\rm{f}}^{\rm{*}} = \sqrt[\root{10}{{4}}] {{\frac{{3b}}{a}}} = 81.794\rm \;s$$ (15) 由图11可见，降低离散点个数对最优终端时间的确定几乎不产生影响，但可以显著降低计算量，如图12所示。在实际工程应用中，可以通过降低离散点个数的方式或者使用低精度的动力学约束进行最优控制模型的求解，以快速获取（tf，mf）样本点。样本点个数不能少于3个，使用最小二乘逼近算法得到a、b、c3个参数，利用式（15）计算出最优终端时间后，再增大离散点个数进行精规划。
图 11不同离散点个数对应的mf-tf曲线
Figure 11.mf-tfcurves under different number of discrete points
图 12不同离散点个数与动力学约束下的计算耗时
Figure 12.Computing time under different dynamics constraint versus discrete points
为验证上述方法的有效性，用计算机随机产生100组初始状态和终端状态，计算出每组条件下的mf-tf曲线，仿真参数如表2所示。仿真结果如图13所示，在100个飞行条件下，除去2个苛刻的飞行情况无解外，其余98个条件下都找到了最优燃料控制的规划轨迹，且有较高的规划成功率。
表 2随机仿真参数设置
Table 2.Simulation parameters
参数 数值 离散点个数N 30 初始坐标r0/m rx0∈[–300，300]
ry0∈[–300，300]
rz0∈[1500，10000]初始速度v0/（m·s–1） vx0∈[–300，300]
vy0∈[–300，300]
vz0∈[–300，0]终端位置rf/m rf=（0，0，0） 终端速度vf/（m·s–1） vf=（0，0，0） tf步长/s 5 
图 13随机初始状态下最优轨迹簇的mf-tf曲线
Figure 13.mf-tfcurves of optimal trajectories under random initial state
-
记本文中不对动力学改造直接采用中心差分进行线性化建立的序列凸优化算法为P0，另外两种凸化方法分别为P1[10]和P2[11]
$$\left\{ \begin{array}{l} {\rm{P1}}:\;\;\;\min \dfrac{{\Delta t{T_{\max }}}}{{{g_0}{I_{{\rm{sp}}}}}} \displaystyle\sum\limits_{i = 0}^{N - 1} u + \left| {{{{w}}^{\rm{T}}}({{x}}{\rm{(}}{t_f}{\rm{)}} - {{{x}}_{_f}})} \right|\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;\;\;\dfrac{{{T_{\min }}}}{{{T_{\max }}}} \text{≤} u \text{≤} 1,\;{{{r}}_z} \text{≥} 0\\ \;\;\;\;\;\;\;\;\;\;\;\;\;n_x^2 + n_y^2 + n_z^2 \text{≤} 1\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\dot{ r}} = {{v}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\dot{ v}} = {{g}} + \dfrac{{u{T_{\max }}}}{m}{{n}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;m = - \dfrac{{{T_{\max }}}}{{{g_0}{I_{{\rm{sp}}}}}}u \end{array} \right.$$ (16) $$\left\{ \!\!\!\! \begin{array}{l} {\rm{P2}}:\min \;\; - {z_{\rm{f}}} + \left| {{{{w}}^{\rm{T}}}({{x}}({t_{\rm{f}}}) - {{{x}}_{_{\rm{f}}}})} \right|\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\sqrt {a_x^2 + a_y^2 + a_z^2} \text{≤} u,\;{{{r}}_z} \text{≥} 0\\ \;\;\;\;\;\;{\dot{ r}} = {{v}}\\ \;\;\;\;\;\;{\dot{ v}} = {{g}} + {{a}}\\ \;\;\;\;\;\;z = - \dfrac{u}{{{g_0}{I_{{\rm{sp}}}}}}\\ \;\;\;\;\;\;{T_{\min }}{{\rm{e}}^{ - {z^{(k)}}}}\left( {1 - (z - {z^{(k)}})} \right) \text{≤} u \text{≤} {T_{\max }}{{\rm{e}}^{ - {z^{(k)}}}}\left( {1 - (z - {z^{(k)}})} \right) \end{array} (17) \right. $$ P1中将推力分解为推力大小u和推力方向n共计4个控制量进行优化，并使用无损凸化理论将方向余弦等式约束
$n_x^2 + n_y^2 + n_z^2 = 1$ 松弛为$n_x^2 + n_y^2 + n_z^2 \text{≤} 1$ 不等式约束，P2定义了新的状态量为z= lnm，这将原问题动力学的一部分非线性转移到不等式约束中。设置3个模型状态量的初始猜想一致（P2中新状态量z全部初始化为lnm0），控制量都全部初始化为0，使用表1中的任务参数求解P0、P1和P2。3种方法对比结果如图14~16所示。结果显示P0和P1求解结果基本一致，仅推力分量曲线有较小差别，P2收敛得更快，但是并未收敛到全局最优解，P2的最优性依赖于更好的初始猜想。
图 143种凸化方法最优轨迹对比
Figure 14.Comparison of optimal trajectories obtained by three convexification methods
图 153种凸化方法求解得到的推力幅值曲线对比
Figure 15.Comparison of thrust magnitude obtained by three convexification methods
表 33种凸化方法的计算结果
Table 3.Results of three convexification methods
收敛次数 剩余质量mf/kg（剩余质量分数/%） 最优性 P0 3 31.338 6（62.68） 最优 P1 6 31.330 4（62.66） 最优 P2 2 30.446 2（60.89） 次优 为将非线性动力学模型转化为凸问题进行求解，三者都使用了等距离散的分段线性动力学进行近似逼近，线性化误差和离散误差不可避免。使用控制量积分得到的轨迹和规划所得轨迹进行对比，积分步长为0.01 s。发现3种凸化方式规划解和积分解误差随时间变化都在同一数量级内。
图 16P0、P1、P2求得规划轨迹与积分轨迹的误差随时间变化曲线
Figure 16.Error versus time between the integral trajectories and the optimal trajectories obtained by methods P0，P1，P2.
在表4设置的参数下，产生500组随机动力下降任务，分别使用3种方法进行求解，观察终端误差散布。仿真结果如图17~20所示。
表 4终端误差散布仿真参数
Table 4.Simulation parameters
参数 数值 初始位置r0/m rx0∈[–3 000，3 000]
ry0∈[–3 000，3 000]
rz0∈[1 500，10 000]初始速度v0/（m·s–1） vx0∈[–300，300]
vy0∈[–300，300]
vz0∈[–300，0]终端位置rf/m rf=（0，0，0） 终端速度vf/（m·s–1） vf=（0，0，0） 终端时间tf/s [60，200] 离散点N 100 
图 173种方法积分所得终端位置误差散布图
Figure 17.Terminal position error of integral trajectories
图 18P0和P2积分所得终端位置误差散布图
Figure 18.Terminal position error of integral trajectories obtained by convexification methods P0 and P2
图 193种方法积分所得终端速度误差散布图
Figure 19.Terminal velocity error of integral trajectories obtained by three methods
图 20P0和P2积分所得终端位置误差散布图
Figure 20.Terminal velocity error of integral trajectories obtained by convexification methods P0 and P2
根据对比，可以发现P1方法的数值稳定性较差，部分算例达到最大迭代次数20次后仍然没有收敛，导致终端误差较大，P0和P2方法终端误差数量级接近。P0相比于P2的位置终端误差散布更加集中，且速度终端误差在z方向上更小。
另外，离散步长对终端误差也都有一定影响。以P0为例，计算表1任务在不同离散点个数下得到的规划轨迹，并通过积分得到终端误差如图21和图22所示。可以看到，终端误差随离散点个数的增加而不断减少。
图 21积分得到终端位置误差随离散点个数变化图
Figure 21.Terminal position error of integral trajectory versus the number of discrete points
图 22积分得到终端速度误差随离散点个数变化图
Figure 22.Terminal velocity error of integral trajectory versus the number of discrete points
固定离散点个数，针对表5设置的随机参数，采用P0在不同终端时间下各打靶200次，如图23和图24所示，积分所得终端误差散布也处于可接受范围内。
表 5终端误差散布仿真参数
Table 5.Simulation parameters
参数 数值 离散点个数N 100 初始坐标r0/m rx0∈[–3 000，3 000]
ry0∈[–3 000，3 000]
rz0∈[1 500，10 000]初始速度v0/（m·s–1） vx0∈[–300，300]
vy0∈[–300，300]
vz0∈[–300，0]终端位置rf/m rf=（0，0，0） 终端速度vf/（m·s–1） vf=（0，0，0） 
图 23不同固定终端时间下积分所得终端位移误差三维散布图及xy平面图
Figure 23.3D and projection onxyplane of terminal position error dispersion of integral trajectories at different fixed terminal times
图 24不同固定终端时间下积分所得终端速度误差三维散布图及xy平面图
Figure 24.3D and projection on XY plane of terminal velocity error dispersion of integral trajectories at different fixed terminal times
对于收敛性，200个随机初始状态中，P0、P1都求出198个解（第43和181个失败），P2求出199个解（第43个失败），三者收敛性无显著差别。两任务的初始速度大，高度低，规划也失败属于合理现象。虽然第181个任务的解被P2找到，但这条轨迹十分危险，在18.18 s处离地面较近且仍然有较大的水平速度，考虑存在离散误差和环境扰动的情况下可能会导致任务失败，两任务参数如下。
$$ \begin{array}{l} 43:\;{{{r}}_0} = {\left[ {\begin{array}{*{20}{l}} {317.972}&{117.466}&{2\;182.24} \end{array}} \right]^{\rm{T}}},\\ \;\;\;\;\;{{{v}}_0} = {\left[ {\begin{array}{*{20}{c}} {70.836\;5}&{ - 274.877}&{ - 279.446} \end{array}} \right]^{\rm{T}}} \end{array} $$ $$\begin{array}{l} 181:\;{{{r}}_0} = {\left[ {\begin{array}{*{20}{l}} { - 2\;725.7}&{ - 790.948}&{2\;673.82} \end{array}} \right]^{\rm{T}}},\\ \;\;\;\;\;\;\;{{{v}}_0} = {\left[ {\begin{array}{*{20}{c}} {228.953}&{ - 182.388}&{ - 298.709} \end{array}} \right]^{\rm{T}}} \end{array}$$ 
图 25P2解得任务181的规划轨迹
Figure 25.The optimal trajectory of case 181 by using P2 convexification method
使用100个离散点求解，所得的仿真结果如图26~28所示。算法P0平均耗时0.129 2 s，平均3.04步收敛。P2平均耗时0.115 7 s，2步收敛。单步迭代P0耗时更少，总耗时P2略少但优化结果P0更优。若定义P0，P2收敛后剩余燃料质量分别为
${m_{P0{\rm{f}}}}$ 和${m_{P2{\rm{f}}}}$ ，定义指标$J = ({m_{P0{\rm{f}}}} - {m_{P2{\rm{f}}}})/{m_0}$ ，统计打靶结果显示J总为正，即P0总是收敛到相对于P2更好的燃料最优解。
图 26P2解得任务181的速度轨线
Figure 26.Velocity profiles of case 181 by using convexification method P2
图 27200组随机任务打靶耗时统计图
Figure 27.Computing time statistics of 200 random cases
图 28200组随机任务中指标J的统计图
Figure 28.The index J statistics of 200 random cases
以上求解过程都未引入信赖域约束，若增加状态量信赖域约束，在表1仿真参数下采用P0进行求解，每种信赖域各计算100次，图29显示对解最优性并未有明显改善。图30为不同信赖域下求解耗时统计。前3条轨线和无信赖域约束求解结果基本重合，第4条轨迹消耗燃料相比前3条节省了0.032 7 kg，仅占总质量的0.065 4%，但耗时相对较长。在本问题中基本可视信赖域为冗余约束。
图 29不同信赖域约束下优化得到的轨线
Figure 29.Optimal trajectories under different trust regions
图 30不同信赖域下求解耗时统计
Figure 30.Computing time statistics under different trust regions
表 6不同信赖域下的求解结果
Table 6.Results under different trust regions.
编号 信赖域约束${{\delta }}{\rm{ }}({\rm{ }}\left| {{{{x}}^{(k + 1)}} - {{{x}}^{(k)}}} \right| \leqslant {{\delta }}{\rm{ }})$ 收敛次数 平均耗时/s 1 $\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{300}&{300}&{300}&{50} \end{array}} \right]$ 3 0.142 4 2 $\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{200}&{200}&{200}&{20} \end{array}} \right]$ 4 0.192 5 3 $\left[ {\begin{array}{*{20}{c}} {5000}&{3000}&{200}&{200}&{200}&{200}&{20} \end{array}} \right]$ 5 0.255 4 4 $\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{100}&{200}&{200}&{200}&{20} \end{array}} \right]$ 9 0.470 2 5 无 3 0.129 5 -
本文使用中心差分逐次线性化动力学，提高了序列凸优化算法的规范性和工程性，避免因使用推力方向余弦无损凸化和质量方程取对数凸化而引入额外的优化变量，降低了公式推导工作量和编程复杂性。通过对比证明，相对于传统动力学凸化技巧，直接线性化的求解结果并未对算法收敛性、解的最优性和终端误差造成显著影响。同时，采用性能指标相对偏差作为迭代收敛条件，使收敛条件和具体模型解耦，具有更好的普适性。
针对动力下降段最优控制问题，给出最优轨迹簇剩余燃料和终端时间的拟合函数，可解析求出近似最优终端时间，能够有效提高在线轨迹方法的效率。
本文使用的动力学模型较为简单，未考虑更多的过程约束和环境条件，如发动机推力方向约束、落地倾角约束和气动力作用，但可为复杂约束下的轨迹规划算法提供参考初始轨迹，也可为轨迹跟踪控制算法提供标称轨迹。
由于存在线性化误差和离散误差二者的复合作用，主要采用大批次数值仿真计算来评估算法收敛能力，对序列凸优化算法收敛性的严格数学证明尚比较困难，需要进一步进行深入研究。
Central Difference Convexification Method for Soft-Landing Trajectory Optimization in Planetary Powered Descending Phase
-
摘要:针对行星定点软着陆实时在线制导的任务需求，设计了基于序列凸优化的动力下降燃料最优轨迹求解算法。算法采用预标记的中心差分法线性化动力学方程，并提出将性能指标相对偏差作为收敛终止条件，能够快速生成燃料最优轨迹。此外，在分析最优轨迹簇剩余燃料和终端时间关系的基础上给出其拟合函数，作为最优终端时间的近似估算，以减少算法求解未知变量的维数。数值仿真结果表明，与对动力学方程先进行变量代换再线性化的传统凸化方法相比，该算法对初始猜想不敏感，收敛性好且终端误差较小。Abstract:Aiming at the mission of real-time online guidance for fixed-point soft landing on planet, the paper designs an algorithm based on sequential convex optimization that is designed to solve the fuel-optimal trajectory. The pre-labeled central difference algorithm is used to linearize the dynamics equations and the termination condition based on the deviation of index is proposed to judge whether it converges, which can generate a fuel-optimal trajectory quickly. Besides, a fitting function is given to approximately estimate the optimal terminal time by analyzing the relationship between the terminal time of optimal trajectories and their remaining fuel, which can reduce the amount of unknown variables. The simulation results of this algorithm show the weak sensitivity to initial guess, good convergence and small terminal error compared to the traditional convexification method of linearizing the dynamics equations whose variables are substituted at first.Highlights
● The workload of code programming and formula derivation can be reduced and the generality of algorithm can be improved by using pre-labeled central difference method instead of analytically computing the Jacobi matrix. ● m f- t fcurves of optimal trajectories can be fitted by a cluster of functions which can be used to approximately estimate the optimal terminal time. ● Simulation results show that direct linearization of dynamic constraints is better than customized convexification method with variable substitution for planet powered landing. -
图 1性能指标及其绝对偏差和相对偏差随迭代次数变化图
Fig. 1Performance index and its absolute deviation，relative deviation versus iterations
图 12不同离散点个数与动力学约束下的计算耗时
Fig. 12Computing time under different dynamics constraint versus discrete points
图 13随机初始状态下最优轨迹簇的mf-tf曲线
Fig. 13mf-tfcurves of optimal trajectories under random initial state
图 143种凸化方法最优轨迹对比
Fig. 14Comparison of optimal trajectories obtained by three convexification methods
图 153种凸化方法求解得到的推力幅值曲线对比
Fig. 15Comparison of thrust magnitude obtained by three convexification methods
图 16P0、P1、P2求得规划轨迹与积分轨迹的误差随时间变化曲线
Fig. 16Error versus time between the integral trajectories and the optimal trajectories obtained by methods P0，P1，P2.
图 18P0和P2积分所得终端位置误差散布图
Fig. 18Terminal position error of integral trajectories obtained by convexification methods P0 and P2
图 193种方法积分所得终端速度误差散布图
Fig. 19Terminal velocity error of integral trajectories obtained by three methods
图 20P0和P2积分所得终端位置误差散布图
Fig. 20Terminal velocity error of integral trajectories obtained by convexification methods P0 and P2
图 21积分得到终端位置误差随离散点个数变化图
Fig. 21Terminal position error of integral trajectory versus the number of discrete points
图 22积分得到终端速度误差随离散点个数变化图
Fig. 22Terminal velocity error of integral trajectory versus the number of discrete points
图 23不同固定终端时间下积分所得终端位移误差三维散布图及xy平面图
Fig. 233D and projection onxyplane of terminal position error dispersion of integral trajectories at different fixed terminal times
图 24不同固定终端时间下积分所得终端速度误差三维散布图及xy平面图
Fig. 243D and projection on XY plane of terminal velocity error dispersion of integral trajectories at different fixed terminal times
图 25P2解得任务181的规划轨迹
Fig. 25The optimal trajectory of case 181 by using P2 convexification method
图 26P2解得任务181的速度轨线
Fig. 26Velocity profiles of case 181 by using convexification method P2
表 1仿真设定参数
Table 1Simulation parameters
参数 数值 初始坐标r0/m （0，0，10 000） 初始速度v0/（m·s–1） （400，–200，–200） 初始质量m0/kg 50 推力幅值[Fmin，Fmax]/N [100，1 000] 终端坐标rf/m （0，0，0） 终端速度vf/（m·s–1） （0，0，0） 发动机比冲g0Isp/（m·s–1） 2 000 终端时间tf/s 60 离散点个数N 100 
下载:
导出CSV
表 2随机仿真参数设置
Table 2Simulation parameters
参数 数值 离散点个数N 30 初始坐标r0/m rx0∈[–300，300]
ry0∈[–300，300]
rz0∈[1500，10000]初始速度v0/（m·s–1） vx0∈[–300，300]
vy0∈[–300，300]
vz0∈[–300，0]终端位置rf/m rf=（0，0，0） 终端速度vf/（m·s–1） vf=（0，0，0） tf步长/s 5 下载:
导出CSV
表 33种凸化方法的计算结果
Table 3Results of three convexification methods
收敛次数 剩余质量mf/kg（剩余质量分数/%） 最优性 P0 3 31.338 6（62.68） 最优 P1 6 31.330 4（62.66） 最优 P2 2 30.446 2（60.89） 次优 下载:
导出CSV
表 4终端误差散布仿真参数
Table 4Simulation parameters
参数 数值 初始位置r0/m rx0∈[–3 000，3 000]
ry0∈[–3 000，3 000]
rz0∈[1 500，10 000]初始速度v0/（m·s–1） vx0∈[–300，300]
vy0∈[–300，300]
vz0∈[–300，0]终端位置rf/m rf=（0，0，0） 终端速度vf/（m·s–1） vf=（0，0，0） 终端时间tf/s [60，200] 离散点N 100 下载:
导出CSV
表 5终端误差散布仿真参数
Table 5Simulation parameters
参数 数值 离散点个数N 100 初始坐标r0/m rx0∈[–3 000，3 000]
ry0∈[–3 000，3 000]
rz0∈[1 500，10 000]初始速度v0/（m·s–1） vx0∈[–300，300]
vy0∈[–300，300]
vz0∈[–300，0]终端位置rf/m rf=（0，0，0） 终端速度vf/（m·s–1） vf=（0，0，0） 下载:
导出CSV
表 6不同信赖域下的求解结果
Table 6Results under different trust regions.
编号 信赖域约束${{\delta }}{\rm{ }}({\rm{ }}\left| {{{{x}}^{(k + 1)}} - {{{x}}^{(k)}}} \right| \leqslant {{\delta }}{\rm{ }})$ 收敛次数 平均耗时/s 1 $\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{300}&{300}&{300}&{50} \end{array}} \right]$ 3 0.142 4 2 $\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{200}&{200}&{200}&{20} \end{array}} \right]$ 4 0.192 5 3 $\left[ {\begin{array}{*{20}{c}} {5000}&{3000}&{200}&{200}&{200}&{200}&{20} \end{array}} \right]$ 5 0.255 4 4 $\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{100}&{200}&{200}&{200}&{20} \end{array}} \right]$ 9 0.470 2 5 无 3 0.129 5 下载:
导出CSV
-
[1] 崔平远,赵泽端,朱圣英. 火星大气进入段轨迹优化与制导技术研究进展[J]. 宇航学报,2019,40(6):611-620.CUI P Y,ZHAO Z R,ZHU S Y. Research progress of trajectory optimization and guidance techniques for Mars atmospheric entry[J]. Journal of Astronautics,2019,40(6):611-620. [2] 于正湜, 崔平远. 行星着陆自主导航与制导控制研究现状与趋势[J]. 深空探测学报（中英文）, 201, 3（4）: : 345-355.YU Z S, CUI P Y. Research status and developing trend of the autonomous navigation, guidance, and control for planetary landing[J]. Journal of Deep Space Exploration, 2016, 3（4）: 345-355. [3] MICHAEL S, ACIKMESE B. Successive convexification for 6-DoF Mars rocket powered landing with free-final-time[C]//2018 AIAA Guidance, Navigation, and Control Conference. [S. l.]: AIAA, 2018. [4] BLACKMORE L,SCHARF D P. Minimum-landing-error powered-descent guidance for Mars landing using convex optimization[J]. Journal of Guidance Control and Dynamics,2010,33(4):1161-1171.doi:10.2514/1.47202 [5] LIU X. Fuel-optimal rocket landing with aerodynamic controls[J]. Journal of Guidance Control & Dynamics,2018,42(46):1-13. [6] WOLF A A,GRAVES C,POWELL R,et al. Systems for pinpoint landing at Mars[J]. Advances in the Astronautical Sciences,2005,119:2677-2677. [7] 李鑫,欧阳高翔,孙成明,等. 基于凸优化的有限推力远程转移轨迹优化[J]. 航天控制,2016,34(3):19-25.doi:10.3969/j.issn.1006-3242.2016.03.004LI X,OUYANG G X,SUN C M,et al. Optimal remote transfer with finite thrust based on convex optimization[J]. Aerospace Control,2016,34(3):19-25.doi:10.3969/j.issn.1006-3242.2016.03.004 [8] LIU X,SHEN Z,LU P. Entry trajectory optimization by second-order cone programming[J]. Journal of Guidance Control and Dynamics,2015,39(2):1-15. [9] DOMAHIDI A, CHU E, BOYD S. ECOS: an SOCP solver for embedded systems[C]//Proceedings European Control Conference. Zurich: [s. n.], 2013. [10] LU P,LIU X. Autonomous trajectory planning for rendezvous and proximity operations by conic optimization[J]. Journal of Guidance Control & Dynamics,2013,36(2):375-389. [11] 林晓辉,于文进. 基于凸优化理论的含约束月球定点着陆轨道优化[J]. 宇航学报,2013,34(7):901-908.doi:10.3873/j.issn.1000-1328.2013.07.003LIN X H,YU W J. Constrained trajectory optimization for lunar pin-point landing based on convex optimization theory[J]. Journal of Astronautics,2013,34(7):901-908.doi:10.3873/j.issn.1000-1328.2013.07.003 -
点击查看大图
计量
- 文章访问数:840
- HTML全文浏览量:377
- PDF下载量:55
- 被引次数:0
