2. 中国商用飞机有限责任公司, 上海 200120
自1962年以来，美国、苏联、欧洲等航天大国先后启动了多达54次的火星探测任务。然而，接近2/3的任务以失败告终。一个关键原因是探测器未成功完成捕获制动程序。由于火星捕获制动策略直接决定了探测器携带燃料的多少，因此分析不同的捕获策略对火星探测任务的成功开展至关重要。
火星具有大气，因此捕获制动策略可分为2类：① 大气辅助捕获制动; ② 推力器点火直接制动。闵学龙等[1]提出了直接制动、一次穿越制动和多次穿越制动3种典型捕获策略，各策略均采用脉冲变轨，并且比较了任务耗时、过载和速度增量等指标；李爽等[2]系统性地总结了火星进入、下降、着陆段 (EDL) 中导航、制导与控制技术的历史与现状，分析了各种技术的优缺点；李军锋等[3]分析了基于有限推力模型的火星探测器捕获制动，设计燃料最优捕获、姿态匀速转动和姿态惯性定向3种直接捕获策略，并提出燃料最优捕获时推力近似沿着速度反方向；方宝东等[4]给出了一种火星探测制动捕获段的工程定义，分析了推力沿速度反方向时各误差源对捕获轨道精度的影响。在工程应用中，欧空局的"火星快车号"探测器采用了姿态匀速转动的制动策略，并采用推力沿速度反方向的策略作为故障预案[5]；“金星快车”探测器则利用推力沿速度反方向的制动策略进入大偏心率的捕获轨道[6]。
事实上，直接制动属于轨道机动的一种。Edelbaum[7]和Gobetz[8]早在20世纪60年代就利用极大值原理间接法解决有限推力变轨问题。但是利用间接法推导最优解的过程较为复杂，估计协态变量初值的难度也较大。近来，学者研究如何将启发式优化算法用于轨迹机动问题。蒋方华等[9]采用同伦方法处理小推力轨迹优化问题，借助粒子群算法得到优化变量的值，并进一步作为另一优化算法的初值；为实现燃料最优或时间最优，Conway等[10]将粒子群算法应用到脉冲推进或小推力推进的轨迹优化问题；Bessette等[11]分析发现粒子群算法在收敛速度和可靠性上要优于微分进化算法。
然而，由于理论推导时涉及简化和假设，工程实际中存在导航、执行误差，探测器很难实现燃料最优，甚至有时需启动故障对策预案。轨道机动消耗的燃料和最终的入轨精度之间通常是"矛盾"的，所以有必要讨论捕获制动时二者之间的关系。
本文把燃料消耗和制动精度同时作为优化目标，采用多目标优化的粒子群算法，借助均匀设计生成初始种群，对推力沿着速度反方向、方向固定和姿态匀速转动3种捕获策略进行优化，得到了Pareto意义下的最优解，并拟合建立两目标之间的函数关系。该函数可用于探测器总体设计阶段，便于研究人员根据探测器设计要求选定适宜的控制参数。
1 捕获问题建模 1.1 动力学模型J2000.0火星赤道惯性坐标系下火星、太阳和火卫1、2的位置如图 1所示。r和v分别为探测器的位置和速度矢量, 且r=‖r‖; rs为太阳相对于火星的位置矢量, 且rs=‖rs‖; rps为探测器相对太阳的位置矢量, 且rps=r-rs; r1、r2分别为火卫1和火卫2在火星赤道惯性坐标系的位置矢量, rp1、rp2分别为探测器相对火卫1和火卫2的位置矢量。
|
| 图 1 天体位置示意图 |
考虑火星的非球形引力和上述天体的摄动力, 探测器的轨道动力学方程如下
|
(1) |
式中:F为捕获制动时推力器产生的推力, 且F=|F|, Um为火星的引力势函数, μs、μ1、μ2分别为太阳、火卫1和火卫2的引力常数, Isp为推力器比冲。相较于推力, 太阳光压等产生的摄动力很小, 可忽略不计。
1.2 捕获策略本文讨论推力方向固定、推力方向匀速转动和沿着速度反方向3种不同的直接捕获策略。有限推力变轨时发动机工作时间较长, 不可能瞬间完成变轨。为了实现有限推力变轨, 需要确定点火三要素:点火角、推力方向、推进时间。本文采用近火点预测的方法确定推力方向, 如图 2所示。
|
| 图 2 推力方向示意图 |
近火点处探测器的位置和速度单位矢量分别为ȓp和
|
(2) |
式中
|
(2) |
w为推力旋转角速度, t0为发动机工作的初始时刻, α0为初始推力仰角; 点火角定义为发动机开始工作时探测器的真近点角。
捕获制动过程中偏心率由大于1单调减小到小于1, 所以可以把轨道偏心率作为制动过程结束的判断条件, 偏心率到达预定值后发动机关机, 因此推进时间可不作为优化的变量。有:
1) 对于推力方向固定策略, 推力仰角保持为常数, 即w≡0, 待优化的参数为点火角θ、推力方向α0;
2) 对于推力方向匀速转动策略, 推力仰角线性变化, 即w≡const, 待优化的参数为点火角θ、推力方向w、α0;
3) 对于推力沿着速度反方向策略, 待优化的参数为点火角θ。
2 多目标优化问题及算法 2.1 多目标优化问题的数学模型一般处理多目标的方法是加权[12], 根据各目标的重要性分配每个目标不一样的权重。但是多目标化为单目标只会输出一个给定权值下的最优解, 忽略了其他权值下的可行解。
本文中, 多目标优化问题[13]的数学描述如下
|
(4) |
式中, x为优化控制变量, y为目标矢量, fi(x) 为第i个目标函数, gj(x) 为第j个约束, S为控制变量的解空间。
定义1 设有控制变量u、v, 且u∈S, v∈S, 若fi(u)≤fi(v), i=1, 2, …m, 且至少有其中一个严格不等式成立, 则称u优于v, 即∀i∈{1, 2, …, k}, ui≤vi∧∃i∈{1, 2, …, k}, ui < vi, 记作u≺v;
定义2 若x*∈S, 且在S中不存在比x*更优的解x, 称x*是多目标优化模型的Pareto非劣最优解;
定义3 对于给定的多目标优化问题f(x), 所有Pareto非劣解构成Pareto非劣解集, 记作P*, 即P*={x∈S|¬∃x′∈S, f(x′)≺f(x)};
定义4 对于一个给定的多目标优化问题f(x), 其所有Pareto非劣解所对应的目标向量构成该问题的Pareto前沿, F*={f=(f1(x), …, fm(x))|x∈P*}。
2.2 粒子群优化算法粒子群优化算法是Kennedy和Eberhart[14]受鸟群运动启发提出的智能优化算法, 具有并行搜索、收敛速度快的优点, 被公认为经典的多目标粒子群优化算法。本文即采用了Coello等[15]提出的多目标粒子群优化算法。该算法使用外部的非劣解集储存非劣解并指导粒子的更新, 采用自适应网格丢弃分布密集的非劣解, 利用轮盘赌法选择非劣解, 并引入变异算子增强粒子的探索能力。
优化流程如图 3。
|
| 图 3 多目标粒子群算法流程 |
初始种群的优劣会对粒子群算法的表现产生重大影响。为实现初始种群的多样性和均匀性特点, 本研究采用均匀设计表[16]初始化种群。
该理论提出了n×s的均匀设计表, 用符号Un(qs) 表示, 其中n表示试验数目, s表示试验所取因素即优化参数个数, q为每个因素的水平数, 在这指每个参数的取值个数。由于种群粒子数目较大, 本文利用好格子点法[16]生成均匀设计的近似解。
假设生成100×2的均匀设计表, 把第j个元素的水平k借助变换
|
(5) |
使k取值落在[0,1]之间。图 4所示为均匀设计和蒙特卡洛随机生成的点集, 可见均匀设计生成的点集可以散布整个搜索空间, 具有更好的均匀性。
|
| 图 4 两点集间的比较 |
火星半径为3.398×103 km, 引力常数为4.28×104km3/s2; 探测器质量为4.948×103 kg, 发动机推力3×103N, 比冲为314 s; 探测器进入半径为5×105 km的火星影响球时, 位置矢量为[-0.822-4.724 1.348]×105 km, 速度矢量[0.519 2.756-0.768] km/s, 根据动力学方程积分得到近火点高度为606.30 km, 双曲线的近火点速度为5.45 km/s, 状态转换为轨道六要素后, 偏心率为1.78。目标轨道的近火点高度600 km, 远火点高度8×104 km, 则偏心率为0.908 5。优化指标如下
|
(6) |
式中, rp为实际近火点位置矢量, rp0为目标近火点位置矢量, mend为发动机关机时刻探测器剩余燃料质量, m0为开机时刻探测器的燃料质量。
粒子群算法参数设置 (参考文献[15]):
种群粒子数目: 100
迭代次数: 100
非劣解集粒子数目:250
个体学习系数: 1.0
全局学习系数: 2.0
变异率: 0.1
3.1 推力沿着速度反方向策略捕获制动本质上是降低轨道能量, 根据轨道单位能量的定义
|
(7) |
式中, E为轨道单位能量, a为轨道半长轴, v为轨道线速度, μ为引力常数。可知, 对于给定的引力场, 轨道单位质量的能量只是半长轴的函数。对于双曲线轨道有E>0, 又a < 0, 所以单位能量随着半长轴的减小而减小。若降低轨道单位质量的能量, 需减小轨道半长轴。根据高斯型摄动方程, 有
|
(8) |
fτ为轨道切向摄动力, 则其他条件相同, 最有效率的捕获制动策略应使得摄动力在切向的分力最大, 即推力应沿着速度的反方向。
由图 5可知, 基本上二者可以同时达到最优, Pareto前沿覆盖范围小, 无需拟合。消耗的燃料最少为1 313.35 kg, 此时近火点高度误差最大为61.62 km, 优化参数点火角为48.49°; 燃料消耗最多为1313.41 kg, 高度误差最小61.46 km, 优化参数值为47.70°。图 6表明点火角的非劣解集小, 平均值为48.88°, 标准差为0.23°, 因此一旦不准时开机变会错过最佳点火角。
|
| 图 5 描述近火点高度误差与燃料消耗的Pareto前沿 |
|
| 图 6 推力沿速度反方向时的非劣解集 |
从制动效率来说, 沿速度反方向是最优的。但是优化沿速度反方向策略时参数自由度较低, 当点火角确定, 则速度方向确定, 推力方向也确定。所以此方法实质上处理的是单目标问题, 减小重力损失达到燃料最优的同时牺牲了入轨精度。
3.2 推力方向固定策略待优化变量的搜索空间见表 1。
由图 7可知, 近火点高度误差与燃料消耗之间的互斥性较弱, 近火点高度误差最大为13.52 km, 此时消耗的燃料最少, 为1 337.36 kg, 点火角和推力方向角分别为[49.35°-0.35°]; 消耗燃料最多为1 339.76 kg, 此时的优化参数值为[54.02°-4.08°]。采用推力方向固定策略时, 发动机工作的弧段长, 重力损失较大, 因此燃料消耗较大。
|
| 图 7 描述近火点高度误差与燃料消耗的Pareto前沿 |
燃料消耗和近火点高度误差之间的数值关系近似线性, 拟合得到函数关系y=-0.190x+1 339.8, 除点x=0外, 其余点的残差均很小, 最大不超过0.2 kg, 因此可以用此拟合函数作为最佳估计, 决策时便于在Pareto前沿上选取合适的点。由图 8可知, 非劣解空间紧凑, 说明合适点火角的范围小。
|
| 图 8 推力方向固定时的非劣解集 |
统计非劣解集解得到优化变量的均值和标准差如表 2所示。
待优化的变量的搜索空间见表 3。
由图 9得, 该策略下燃料消耗和近火点高度误差之间的矛盾性强于推力方向固定时, 消耗的燃料质量最少为1 313.11 kg, 此时近火点的误差达到67.16 km, 优化参数点火角、初始推力仰角和旋转角速度分别为[48.30°-21.99° 0.0298°/s]; 近火点高度误差最小为0, 消耗的燃料为1 330.05 kg, 优化参数为[48.27°-3.88° 0.005°/s]。二者之间的数值关系近似二次函数型, 拟合得到函数关系为y=0.003 9x2-0.510 7x+1 330, 除点x=0外, 其余点的残差也较小, 最大也不超过0.2 kg。
|
| 图 9 描述近火点高度误差与燃料消耗的Pareto前沿 |
图 10为Pareto前沿对应的非劣解集, 合适的点火角范围较小。统计位于Pareto前沿的解得到均值和标准差如表 4所示。匀速转动的初始推力仰角为负值, 推力有垂直速度向内的分力, 随着时间增大推力仰角变为正值, 推力有垂直于速度向外的分力, 因此最终抬升了轨道近火点高度。
|
| 图 10 推力方向匀速转动下Pareto前沿对应的非劣解 |
| 优化变量 | 平均值 | 标准差 |
| 点火角/(°) | 48.63 | 0.55 |
| 初始推力仰角/(°) | -10.16 | 4.67 |
| 旋转角速度/(°·s-1) | 0.014 | 0.006 |
绝大部分时间内推力方向与速度的反方向之间的夹角在10°以内，因此匀速转动可看作是沿速度反方向策略下姿态角均匀变化的版本，可以减小重力损失，继而减小燃料消耗。比较图 5和图 9，发现推力方向匀速变化时燃料最省的情况与沿着速度反方向消耗的燃料非常接近，约为1 313 kg，这也说明二者之间的相似性。
对于3种策略, 生成的Pareto前沿均是凸函数，多目标的粒子群算法尤其擅长解决此类问题。统计发现3种策略的点火角的平均值相近，这为一种策略下发动机未按时点火，系统转为另一种策略的可行性提供了依据。
4 结论本文首先介绍了捕获问题的数学模型，明确了捕获制动过程中3种策略各自待优化的参数。然后提出了以近火点高度误差和燃料消耗为双目标的优化问题，采用多目标粒子群算法进行优化并利用均匀设计初始化种群，寻找Pareto意义下的最优解和Pareto前沿，比单目标优化提供更多的决策信息。经仿真，发现Pareto前沿均是凸函数，即两目标之间具有矛盾性。在不考虑姿态控制问题的情况下，比较分析3种策略在制动效果上的表现：
1) 推力方向固定的策略燃料消耗多，近火点高度误差小，二者的Pareto近似线性关系；
2) 沿着速度反方向策略捕获效率高，燃料消耗少，但近火点高度误差较大；
3) 推力方向匀速转动介于两者之间，有抬升近火点高度的趋势，但姿态转动使得推力方向与速度反方向保持较小的角度，效率提高，燃料消耗与近火点高度误差之间的Pareto前沿近似二次函数关系。
以上讨论并没有考虑姿态控制问题，推力方向匀速转动时探测器的姿态不断变化，需启用姿控发动机跟踪预定角度。由于姿态角速度较小且为匀速转动，工程上可以实现。对于推力沿着速度反方向的策略，轨迹的切线方向变化不恒定，实时速度的获取依赖于较高的导航精度，因此系统更为复杂。
推力方向匀速转动的策略可以作为优先选用的策略，当出现故障需紧急制动时可启用推力沿速度反方向策略，控制探测器先被火星捕获。在总体设计阶段，拟合函数可提供重要参考信息，便于决策。
| [1] |
闵学龙, 潘腾, 郭海林.
火星探测器使命轨道捕获策略研究[J]. 航天器工程, 2008, 17 (6): 39–43.
Yan Xuelong, Pan Teng, Guo Hilin. Analysis of Orbit Capture Method for Mars Vehicle[J]. Spacecraft Engineering, 2008, 17(6): 39–43. (in Chinese) |
| [2] |
李爽, 彭玉明, 陆宇平.
火星EDL导航, 制导与控制技术综述与展望[J]. 宇航学报, 2010 (3): 621–627.
Li Shuang, Peng Yuming, Lu Yuping. Review and Project of Mars EDL Navigation Guidance and Control Technologies[J]. Journal of Astronautics, 2010(3): 621–627. (in Chinese) |
| [3] |
李军锋, 龚胜平.
有限推力模型火星探测捕获策略分析[J]. 中国科学:物理学, 力学, 天文学, 2013 (6): 781–786.
Li Junfeng, Gong Shengping. Analysis of Capture Strategies for Mars Explorer with Finite-Thrust[J]. Scientia Sinica (Physica, Mechanica & Astronomica), 2013(6): 781–786. (in Chinese) |
| [4] |
方宝东, 吴美平, 张伟.
火星引力捕获动力学与动态误差分析[J]. 力学学报, 2015, 47 (1): 15–23.
Fang Baodong, Wu Meiping, Zhang Wei. Mars Gravity Capture Dynamic Model and Error Analysis[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(1): 15–23. DOI:10.6052/0459-1879-14-327 (in Chinese) |
| [5] | Hechler M, Yáñez A. Mars Express Orbit Design[J]. Acta Astronautica, 2003, 53(4): 497–507. |
| [6] | Sanchez Perez J M, Canabal R. Review of Venus Express Mission Analysis[C]//18th International Symposium on Space Flight Dynamics, 2004:245 |
| [7] | Edelbaum T N. Optimal Low-Thrust Transfer between Circular and Elliptical Orbits[J]. Journal of ASME, 1962, 84(2): 134–141. |
| [8] | Gobetz F W. Optimal Variable-Thrust Transfer of a Power-Limited Rocket between Neighboring Circular Orbits[J]. AIAA Journal, 1964, 2(2): 339–343. DOI:10.2514/3.2281 |
| [9] | Jiang F, Baoyin H, Li J. Practical Techniques for Low-Thrust Trajectory Optimization with Homotopic Approach[J]. Journal of Guidance, Control, and Dynamics, 2012, 35(1): 245–258. DOI:10.2514/1.52476 |
| [10] | Pontani M, Conway B A. Particle Swarm Optimization Applied to Space Trajectories[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(5): 1429–1441. DOI:10.2514/1.48475 |
| [11] | Bessette C R, Spencer D B. Optimal Space Trajectory Design:a Heuristic-Based Approach[J]. Advances in the Astronautical Sciences, 2006, 124: 1611–1628. |
| [12] | Coello C A C. An Updated Survey of Evolutionary Multiobjective Optimization Techniques:State Of The Art and Future Trends[C]//Proceedings of the Congress on Evolutionary Computation, 1999:3-13 |
| [13] | Zitzler E. Evolutionary Algorithms for Multiobjective Optimization:Methods and Applications[M]. Ithaca, Shaker: , 1999. |
| [14] | Kenndy J, Eberhart R C. Particle Swarm Optimization[C]//Proceedings of IEEE International Conference on Neural Networks, 1995:1942-1948 |
| [15] | Coello C A C, Pulido G T, Lechuga M S. Handling Multiple Objectives with Particle Swarm Optimization[J]. IEEE Trans on Evolutionary Computation, 2004, 8(3): 256–279. DOI:10.1109/TEVC.2004.826067 |
| [16] |
方开泰.
均匀设计[J]. 应用数学学报, 1980, 3 (4): 363–372.
Fang Kaitai. Uniform Design[J]. Acta Mathe-Maticae Applicatae Sinica, 1980, 3(4): 363–372. (in Chinese) |
2. Commercial Aircraft Corporation of China Ltd, Shanghai 200120, China
