“地球是人类的摇篮,但是人类不能永远生活在摇篮里”。在2019年5~6月的第十届国际空间轨道设计大赛(Global Trajectory Optimization Competion—GTOC X)上,美国宇航局(NASA)喷气推进实验室(Jet Propulsion Laboratory, JPL)向全世界提出了这样一个问题：假如我们不再满足于做太阳系公民,该如何向银河系移民呢?作为第十届竞赛,其题目恰与赛事发起人Dario Izzo给第一届赛事的总结致辞遥相呼应：”也许有点浪漫的是,我们将GTOC 看作是一种’航海挑战’,我们的银河系是竞赛水域,我们的数学工具是竞赛船只^[1]。”

欧洲航天局(ESA)于2005年发起了国际空间轨道设计大赛,旨在促进深空探测科学与技术发展,探索新型轨道设计理念与方法^[2]。此后,该赛事每隔一到两年举办一次,竞赛期间,全世界最优秀的航天工程师和数学家们将挑战”几乎不可能--的空间轨道设计问题^[3]。历届赛题及其解决方案,成为了科学界一笔宝贵的经验财富。

本届比赛由上届冠军JPL主办,吸引了来自全球的73支队伍报名参赛,其中包括了各国航天机构、高等院校和商业航天公司等。最终,笔者所在的国防科技大学和西安卫星测控中心联队(NUDT & XSCC)以显著优势夺得了此次比赛的冠军,其后依次为清华大学联队、ESA代表队和美国宇航公司代表队(The Aerospace Corporation)等^[4],如图1所示。

本文介绍了比赛的题目,并对竞赛过程中笔者所在团队的解决方案进行了描述,最后总结了比赛的意义及其带来的启发。

本届赛题以银河系移民为背景,要在9千万年的时间里,利用三种星际飞船帮助人类在银河系扩散开来。

题目选择了银河系中的10万颗恒星作为可移民目标,它们分布在直径约21万光年的范围内。移民计划将于银河纪元0年开始,人类需要在1千万年内从太阳系启程,搭乘3艘母船或2艘快船前往银河系各个角落。每艘母船可以携带10个移民舱,它们可以在母船飞越恒星时被释放,移民舱通过速度机动来抵达该恒星,随后它就会被移民。每艘快船只可以和一颗恒星交会,被交会的恒星随即被移民。当一颗恒星被移民2百万年后,就可以发出不超过3艘移民船,每艘移民船飞向一颗恒星,它们随后会被移民。就这样我们从太阳出发飞到初代移民星,又从这里出发飞向更多的恒星,人类在银河系中徐徐扩散,一直持续到银河纪元9千万年时终止^[5]。移民过程如图2所示。

在银河系中,飞船和恒星的运动受到银河系中心的引力作用,但其变化规律不同于近地空间的二体引力模型。题目假设所有的恒星都围绕银心做圆周运动,并根据银河系恒星的观测数据,给出了一个近似的银河系中心引力模型^[5]

其中,$r$为质点与银心的距离,$v_{c}$为该点的圆轨道线速度,$f$为其受到的引力加速度。

在二体模型中,引力与半径为平方反比关系,而本届赛题给出的银河系中心引力与半径为高次多项式反比的关系,因此在求解恒星际转移轨道时不能直接套用二体引力环境下常见的公式与方法。

所有的飞船都使用脉冲推力机动,也就是说,它们的速度变化可以在一瞬间完成。但是,飞船的机动次数和脉冲大小要受限于一些约束条件。一艘母船至多可以进行3次机动,单次机动的脉冲不超过200 km/s,而其总脉冲不能超过500 km/s。每个移民舱只能机动1次,其脉冲不能超过300 km/s。每艘快船可以进行2次机动,单次脉冲不得超过1500 km/s,而总脉冲也不能超过1500km/s。移民船最多可进行5次机动,单次机动不超过175 km/s,总脉冲不能超过400km/s。要注意的是,任意一艘飞船两次机动之间的时间间隔不得小于100万年,移民舱的机动与其母船前一次机动的时间间隔也要满足此约束^[5]。各种飞船的机动如图2所示。

移民结束时,所有的移民序列组成一个完整移民方案,其得分由下式确定

式中,$J$为得分,$B$为时间奖励系数,$N$为被移民恒星数量, $E_r$为径向均匀性指标, $E_{\theta}$为周向均匀性指标, $\Delta V_{\max}$为所有飞船可使用的最大脉冲量, $\Delta V_{\rm used}$为所有飞船实际使用的脉冲量。

由式(2)可以看出,每个移民方案的分数将由三部分决定：时间奖励系数、移民星数量及均匀性、燃料系数。

(1)时间系数

本届比赛采用实时计分与在线排名机制,榜单如图1所示。在题目发布一周后,开放提交窗口,持续22天后,关闭提交窗口。时间奖励系数与方案的提交时间有关,其计算方式为

式中,$t_{\rm start}$为提交窗口开放时间,$t_{\rm submission}$为实际提交时间,$t_{\rm end}$为提交窗口关闭时间。

根据式(3),在最开始就提交的方案将得到2倍的分数,到了比赛结束时,该系数将衰减到1.0。另外,这样的时间系数在比赛前期下降得很快,移民方案晚交一天就会产生很大的分数损失。而到了比赛后期,时间系数的变化就很缓慢,一个方案在不同时间提交得到的分数差异不大。时间系数动态变化及实时排行榜机制给竞赛引入了博弈对抗性,增加了比赛的趣味性。

(2)移民星数量及均匀性

得分函数的第二项由被移民恒星的数量和它们的分布均匀性决定。在均匀性相同的情况下,被移民恒星的数量越多越好。均匀性指标对得分的影响至关重要,可以近似认为,均匀性指标每下降1,相当于多移民了1000多颗恒星(均匀性指标越低,此收益还会更高)。但是,要想改善均匀性指标却困难重重,这是因为,均匀性由两个密切耦合的部分构成：径向均匀性和周向均匀性。径向均匀的要求是：将整个银河系由内到外切成30条等宽度的环形区域,每个区域内被移民恒星数要与该区域的面积成正比。而周向均匀的要求是：将整个银河系沿周向切成32块等角度扇形区域,每个区域内被移民恒星的数量要相等。

(3)燃料系数

得分函数的第三项是燃料系数,由所有飞船可用的总燃料和所有飞船实际消耗的总燃料相除得到。这个系数要求在移民过程中使用尽可能少的燃料。该系数对分数也有很大的影响,在比赛初期,一些队伍只移民了5颗恒星便得到了很高的分数,就是最大化利用了燃料系数这一项,从银河系中挑选出5颗最省转移燃料的恒星作为目标,再分别用3艘母船和2艘快船与之交会即可。到了比赛后期,当分数达到了一两千,那么燃料系数每产生0.1的变化,结果就会相差一两百分。

向银河系移民的过程,就像是我们从太阳系发出飞船,在宇宙中播下几粒种子,随后生根发芽,树根长出树干,树干长出树枝,树枝长出树叶,最终形成参天大树的过程。而这几颗”大树--就代表了我们的移民方案,要想得高分就要让大树长得”枝繁叶茂”(数量多),还要能相互配合以”荫庇四方”(分布均匀),另外,在它们生长过程中还要尽量”节水缩肥”(省燃料)。

要得到一个完整的移民方案至少要确定以下几个要素：要移民哪些恒星?恒星被移民的顺序是什么样的?两两之间星际转移的起止时间是多少?要想得到全局最优的移民方案,这些要素都要作为问题的设计变量来考虑,使得这个问题的解空间变得超级庞大。就空间复杂度而言,银河系中有10万颗可移民恒星,每颗恒星有”移民”或”不移民”两种状态,那么就有2¹⁰⁰⁰⁰⁰种可能组合。相对的,一个标准围棋的空间分布约为10¹⁷² 种^[6-7],远低于银河系移民分布数量。如果要设计一个空间复杂度与银河系移民相当的围棋棋盘,其面积至少有半个网球场那么大。对于策略复杂度来说,围棋大约需要考虑250步(决策深度),每一步大约有150种可能落子方式(决策宽度),共有250¹⁵⁰ 种落子序列^[6,8]。在银河系移民问题中,如果要获得全局最优的移民序列,需考虑上千步,每一步至少有近万种可能的选择,其策略复杂度远超围棋。

面对这样大规模的难题,显然不能用简单的暴力搜索来进行全局寻优,那样即便是用尽宇宙诞生以来的全部岁月都无法搜索完成,我们必须将问题进行分解,并以不同的策略来求解各个部分。

移民银河系的第一个问题就是：我们要将这些种子播撒在何方?最开始,我们猜测沿着银河系旋臂分布的一条螺旋线会是比较好的分布,因为旋臂上的恒星数量最多,而且这也能兼顾各个角度和半径的分布。在试算验证的过程中,我们对均匀性指标的认识也逐渐深入,最终通过分析推导得到了理想最优螺旋线的数学表达式^[9]

$ \theta =\frac{2π}{ r_{\max}^{2}-r_{\min}^{2} }r^2+C $

其中,$\theta$为角度,$r$为半径,$r_{\max}$和$r_{\min}$分别为移民范围最大、最小半径,$C$为常数,其形状如图3所示。

其中,红色线条表示的是理想螺旋线的空间位置,蓝色线条的长度反映了该弧段内的理想恒星数量。

有了理想分布做指导,接下来就要设计如何将种子播撒到这些地方了。从太阳系出发的5艘飞船中,2艘快船的机动能力很强,可以认为它们能飞到银河系任意地点。不过3艘母船就有点力不从心了,即便飞行9000万年也只能到达太阳附近的1/3个银河系。考虑到母船可以携带很多移民舱,我们决定先安排好母船的飞越位置,再安排快船的交会目标。而母船最多只能机动3次,这就意味着,它进行一次机动需要与多颗恒星”擦身而过”,达到”一石多鸟--的效果,如图2中的蓝色轨迹所示。

银河系中有10万个目标,假使我们只挑选其中4颗来飞越,那也有400亿亿种组合。鉴于此,我们对各母船的搜索区域做了限制,再结合飞越约束条件,筛选掉了大量无效组合,并在国防科技大学天河超级计算机的帮助下找到了一些可用的飞越方案^[9]。这时再依据螺旋线设计原则,挑选出一些比较好的母船飞越目标和快船的交会目标。这样就确定出了所有初代移民星(即种子)的位置。

播撒下初代种子后,接下来就要设计在9000万年内,它们要长出什么样的树干、树枝和树叶。

一开始,我们选择了正向生成树的办法,从种子出发,通过预设的规则让其一层层长出枝叶。经过对生长规则的不断改进,我们用这种方法获得了466分的移民方案。而此时,ESA已经以548分位居榜首,而且其均匀性指标比我们好很多,这就迫使我们反思这种正向生成策略的问题。认真分析后发现,正向生成的过程中我们是不知道最终移民分布的,但是均匀性指标却偏偏是以最终分布为依据计算的,因此在生成期间选定的目标并不一定是全局最优的。这就好比我们没带地图去登山,我们不知道最高的山峰在哪个方向,只能利用眼前的景象先选一个方向走走看看,那这个方向大概率不通向最高峰。

这时,另一种截然相反的想法浮出水面,我们不妨先选好这颗树的叶节点和中间节点,然后反过来生成枝干将它们连接起来,这就是”反向生成树”的策略^[9] (图4)。随后迅速将这个想法付诸实践,按此方法我们很快找到了一个598分的解,其均匀性指标得到了大幅改善。这种反向生成策略的优势就在于,我们能事先设计好移民的最终分布,再来考虑这种分布的实现过程,这样就能直奔主题,以最终得分为导向进行移民。使用反向生成策略,就像是拿着地图去登山,一开始就知道了哪个方向是高峰,不再会被眼前的层峦叠嶂迷惑到了。

在改善了移民分布后,我们开始思考燃料系数是否有改进空间。在移民过程中,被移民恒星的数量几乎是以指数增加的,到了移民后期会有大量飞船可供使用,而这时移民树的主体已接近完成,移民恒星数量的增长并不是很迫切的需求,那我们正好可以利用这些飞船来提高燃料系数,给它们充足的飞行时间来减少燃料消耗。这种措施的效果立竿见影,我们在移民分布几乎不变的情况下,仅靠改善燃料系数就将分数提高到了700分左右。

此时,我们已经与第二名拉开了差距,但团队并没有因此松懈,而是更加投入了。在思维的碰撞中,新的想法、新的创意被接二连三地提了出来。

有队员提出,我们已经生成了一颗树,不妨再对它修剪修剪,然后让它接着生长,对提升分数肯定有帮助。于是,我们很快编写了剪枝算法,删掉一些对得分不利的节点,再利用先前的正向生成算法,让树枝长出一些能增加得分的节点。剪枝与增枝过程如图5所示。

接着,又有队员想到,我们关注最多的是各个节点,其实节点之间的连线也有改善的余地,每条连线代表着一次星际转移,而转移时间直接影响了燃料的消耗。这次,我们以燃料系数为指标,对所有节点之间的转移时间进行优化,使得所有的连线变得更加”节能”。

随后,另一个队员提出,我们在不改变树结构的情况下,还可以让各节点与附近的未移民恒星进行交换。这就像是化学中的置换反应,我们将高分的恒星”拉拢”进来,将低分的恒星”置换”出去。移民树节点置换如图6所示。

经过团队的不懈努力,我们将这些想法逐一实现,结果证明它们都能在一定程度改善我们反向生成的那些”移民树”。我们不断改善初始种子的选择,再对这些方法综合利用,我们的得分一路提升,依次突破了1000分、1500分、2000分......当分数达到2400分后,我们用这些方法却很难再突破2500分了。在这个关头面前,我们已有的几种方法都无能为力,已经到了”按下葫芦起来瓢”的地步,明显可以看到得分函数中的第二项和第三项在此消彼长,难以取得全面突破。

我们对所有的方法重新梳理一遍,分析了各种方法的机理和影响。这时,我们注意到了前面的剪枝增枝、时间优化和置换节点等操作对移民树的各节点产生了很大改动,但是节点之间的连线拓扑结构并没有改变。也就是说,最开始那些长得均匀有致的树枝,经过了多轮局部优化后,可能变得”奇形怪状”。”不破不立！我们干脆将这颗奇形怪状的树给彻底打断,再使用反向生成树的方法将它重新连接起来,这样它的拓扑结构不就得到更新了吗?”队伍里萌生了这样的想法。在大家的帮助下,很快得到印证,重连后的移民树具有良好的结构,我们又能全面地改善各项指标了。移民树拓扑重连的过程如图7所示,重连前后,被选中的移民恒星不变,但是它们之间的连接关系发生了改变。

这种方法的作用,很像是金属工艺里面的”去应力退火”,对工件进行了各种加工操作后,其内部微观组织的体积会发生不均匀变化,由此产生了内应力。去应力退火工艺可以在不改变组织形态的条件下,通过调整内部结构来消除大部分内应力。我们这种打断重连的操作,就相当于消除了移民树中的不合理连接,改善了所有移民序列的拓扑结构。

到了比赛的最后三天,我们已经有两个不同的方案都接近了2800分,但是都很难再取得大幅度的突破了。我们将手里的几个移民方案都画了出来,有队员对其进行了仔细对比后提出,这几种不同方案之间是可以相互借鉴的。比如,A方案在半径方向上的分布较好,而B 方案在角度上的分布较好,要是我们将A方案中的部分移民树提取出来,”移栽”到B方案上,这样不就能取长补短,得到在两个方向上都具有良好分布的方案了吗。移民树嫁接移植的过程如图8所示。

“优秀的方案肯定会有优美的图案。”经过不懈的努力,我们得到这样一幅移民图景——螺旋臂银河系的星图上,繁星密布,红黄蓝绿紫的移民轨迹在其间均匀排布。

最终的冠军方案如图9所示,蓝色背景为移民终止时所有恒星的位置,红色的五角星是太阳系,14个黑点表示初代移民星的位置,左侧的浅蓝、绿、紫3块分别为3艘母船移民的结果,右侧的黄、红2块为2艘快船移民的结果,点与点之间的连线表明有移民船从一颗恒星出发飞到了另一颗恒星(不是真实飞行轨迹),该方案的动画演示见附录。

国际空间轨道设计大赛从发起至今已有15年的时间了^[1],经过多年的积累,特别是由中国力学学会和清华大学发起的国内空间轨道设计竞赛的锻炼^[10-12],终于有中国队伍夺得了该赛事的冠军。表明中国在轨道设计领域已经达到了世界一流水平,国际同行对中国参赛队的成绩和实力给予了高度赞誉。这样的成绩吸引了广泛的社会关注,央视新闻联播、新华社、人民日报和解放军报等媒体对此进行了深度报道。

当然,我们还应该继续向国际同行学习。JPL在本届赛事中提出的题目可谓是”脑洞大开”,当我们还在关注近地空间和太阳系内的深空探测时,JPL已经将视野放到了银河系,并在赛题中提供了详细的银河系动力学模型和丰富的恒星轨道数据,并提出了全新的树形动态搜索问题。我们应该学习这种创新精神,在未来的研究中提出新问题,探索新理论。本届比赛过程也启示我们先进的数学建模思想和手段对解决问题有很大的带动作用,轨道设计进一步发展需要汲取数学新理论。

本届比赛也极大地锻炼了参赛的学生^[13]。笔者所在的冠军团队中,有多名在读博士和硕士研究生,在为期4周的高强度竞赛中,经常需要快速实现一种新的想法和思路,这样的过程迅速激发了他们的创新潜力,提升了动手解决实际问题的能力。竞赛也使学生充分体会到了团队合作攻关的魅力,更培养了学生的执行力,增强了他们攀登世界科技高峰的自信心。

附录

为直观理解银河系移民过程,我们给出了一个演示视频,可扫描二维码观看。

该二维码链接的视频展示了本届比赛的冠军方案,人类从太阳系释放了5 艘(2 绿3 紫飞行曲线) 探测器开始了漫漫移民之路,历经9 000 万年的不懈努力,在银河系10 万颗星星中成功移民了3 798 颗(橙色移民恒星)。