譚迎龍，喬 兵?，朱嘯宇，靳永強(qiáng)，張慶展

（1.南京航空航天大學(xué)航天工程系，南京210016；2.上海交通大學(xué)航空控制科學(xué)與工程系，上海200240；3.上海宇航系統(tǒng)工程研究所，上海201109）

1 引言

隨著航天技術(shù)的快速發(fā)展，在軌航天器的數(shù)量逐漸增多，導(dǎo)致空間資源日趨緊張。同時(shí)，技術(shù)的進(jìn)步使航天器攜帶設(shè)備發(fā)生故障的概率大幅度下降，燃料的限制已成為制約航天器長(zhǎng)期在軌的主要因素［1］。因此，各航天大國(guó)紛紛提出“航天器在軌加注”的設(shè)想以加強(qiáng)在軌航天器的機(jī)動(dòng)能力和工作壽命［2］，并在此基礎(chǔ)上進(jìn)行了多次“一對(duì)一”航天器在軌加注的驗(yàn)證性實(shí)驗(yàn)，實(shí)驗(yàn)證明在軌加注可行性的同時(shí)，也表明進(jìn)行“一對(duì)一”在軌加注其成本往往比更換一顆新航天器還高［3］。為了獲得更好的效費(fèi)比，一顆服務(wù)航天器對(duì)多顆目標(biāo)航天器進(jìn)行加注的“一對(duì)多”在軌加注模式被提出［3］。研究表明，“一對(duì)多”在軌加注模式優(yōu)于“一對(duì)一”在軌加注模式，因此，“一對(duì)多”模式必將會(huì)成為主流的研究方向［3-4］。隨著目標(biāo)航天器數(shù)量的增加，在有限的時(shí)間內(nèi)，一顆服務(wù)航天器將無(wú)法完成整個(gè)在軌加注任務(wù)，因此根據(jù)任務(wù)的需要，多顆服務(wù)航天器將會(huì)同時(shí)進(jìn)行各自的“一對(duì)多”在軌加注任務(wù)，此時(shí)的在軌加注模式便稱(chēng)為“多對(duì)多”加注模式［3］。

地球靜止軌道（Geostationary Orbit，GEO）運(yùn)行著預(yù)警、通信及中繼衛(wèi)星等大量高價(jià)值航天器，并且由于此類(lèi)航天器的大量發(fā)射導(dǎo)致GEO軌道資源稀缺，且此類(lèi)航天器主要運(yùn)行在GEO帶附近，分布集中。因此，GEO航天器比較適合作為在軌加注目標(biāo)航天器。為了給“一對(duì)多”和“多對(duì)多”在軌加注試驗(yàn)實(shí)施提供理論基礎(chǔ)，已經(jīng)有學(xué)者對(duì)其展開(kāi)了相關(guān)的研究，歐陽(yáng)琪［3］對(duì)地球同步軌道衛(wèi)星群在軌加注任務(wù)規(guī)劃問(wèn)題進(jìn)行了研究，成功地將在軌加注問(wèn)題轉(zhuǎn)化成TSP問(wèn)題求解。Shen等［5-6］解決了基于多圈Lambert問(wèn)題的最優(yōu)雙脈沖問(wèn)題，指出“一對(duì)多”在軌加注調(diào)度的目的是找到總?cè)己淖钚〉淖顑?yōu)加注服務(wù)順序以及最佳時(shí)間分配方案，最優(yōu)服務(wù)順序總是在一些總掃描角最小的序列中。在文獻(xiàn)［6］的基礎(chǔ)上，Shen 等［7］又研究了三種不同的在軌加注模式，研究表明混合P2P加注模式有更少的燃料消耗量。Zhou等［8］對(duì)基于燃料站的“一對(duì)多”加注模式的在軌加注任務(wù)規(guī)劃問(wèn)題進(jìn)行了研究，結(jié)果表明最優(yōu)加注順序是按目標(biāo)航天器真近點(diǎn)角位置依次加注。

本文將在傳統(tǒng)“多對(duì)多”在軌加注模式的基礎(chǔ)上引入一個(gè)空間燃料站，構(gòu)成一個(gè)擁有目標(biāo)航天器（Target Spacecraft， TS）、服務(wù)航天器（Servicing Spacecraft， SSc）和空間燃料站（Propellant Depot，PD）的在軌加注網(wǎng)絡(luò)系統(tǒng)，為了敘述方便，該加注網(wǎng)絡(luò)系統(tǒng)所形成的加注模式稱(chēng)為“k-1-n”模式。同時(shí)，將基于“k-1-n”模式，考慮任務(wù)時(shí)間限制、航天器攜帶燃料量限制以及目標(biāo)航天器所需燃料的不同，開(kāi)展對(duì)GEO軌道航天器在軌加注任務(wù)調(diào)度問(wèn)題的研究。

2 數(shù)學(xué)模型

2.1 問(wèn)題描述

基于空間燃料站的“多對(duì)多”在軌加注模式如圖1所示，空間部署著三種航天器。在燃料站軌道（Obrit1）部署著一個(gè)空間燃料站，在加注任務(wù)開(kāi)始前k顆服務(wù)航天器與空間燃料站處于對(duì)接狀態(tài)，其作用是負(fù)責(zé)把空間燃料站的燃料輸送給目標(biāo)航天器。n顆目標(biāo)航天器位于GEO軌道，各目標(biāo)航天器之間有一定的相位差。當(dāng)目標(biāo)航天器提出加注需求，在規(guī)定的時(shí)間內(nèi)，服務(wù)航天器進(jìn)行軌道機(jī)動(dòng)，對(duì)目標(biāo)航天器進(jìn)行加注作業(yè)，當(dāng)服務(wù)航天器的燃料不足以完成加注任務(wù)時(shí)，服務(wù)航天器可在任務(wù)中途返回燃料站獲取燃料后繼續(xù)進(jìn)行加注任務(wù)。服務(wù)航天器完成所有目標(biāo)航天器的加注作業(yè)后需返回空間燃料站，等待下一次的加注任務(wù)，在無(wú)加注任務(wù)的時(shí)間，服務(wù)航天器將從空間燃料站進(jìn)行燃料補(bǔ)給，為下次加注任務(wù)做準(zhǔn)備。

圖1 “k-1-n”在軌加注模式示意圖Fig．1 Diagram of k-1-n on-orbit refueling strategy

2.2 模型分析

本文的共面圓軌道轉(zhuǎn)移將采用多圈Lambert軌道轉(zhuǎn)移，并參考文獻(xiàn)［5］的研究成果優(yōu)化軌道轉(zhuǎn)移的時(shí)間分配方案。當(dāng)r1＝r2時(shí)，即目標(biāo)航天器與服務(wù)航天器同處于GEO軌道，在求得初始真近點(diǎn)角差后，便能求解多圈Lambert軌道轉(zhuǎn)移過(guò)程中最小速度增量和轉(zhuǎn)移時(shí)間的關(guān)系，如圖2虛線(xiàn)所示，初始真近點(diǎn)角差θ＝60o。兩者的關(guān)系并不遞減，但極小值點(diǎn)呈現(xiàn)逐漸減小的趨勢(shì)，文獻(xiàn)［5］利用這一特點(diǎn)，采用軌道轉(zhuǎn)移末端滑行的方法，使最小速度增量與轉(zhuǎn)移時(shí)間呈現(xiàn)出遞減的關(guān)系，近似階躍函數(shù)，如圖2實(shí)線(xiàn)所示，可知在一段時(shí)間內(nèi)燃料消耗量是不變的。因此，可用一個(gè)整數(shù)作為該時(shí)間區(qū)間的索引號(hào)，定義Index＝a表示第a個(gè)時(shí)間區(qū)間，如圖3所示，故可將對(duì)轉(zhuǎn)移時(shí)間的優(yōu)化轉(zhuǎn)化成對(duì)轉(zhuǎn)移時(shí)間區(qū)間索引號(hào)的優(yōu)化。

圖2 最小速度增量和轉(zhuǎn)移時(shí)間的關(guān)系［5］Fig．2 Relationship between the minimum speed increment and transfer time［5］

圖3 最小速度增量和轉(zhuǎn)移時(shí)間關(guān)系（定義索引號(hào)）Fig．3 Relationship between the minimum speed increment and transfer time（with index）

對(duì)于加注任務(wù)而言，完成任務(wù)的時(shí)間越短越好，因此，在給定時(shí)間區(qū)間索引號(hào)后，選擇該區(qū)間最小的轉(zhuǎn)移時(shí)間作為軌道轉(zhuǎn)移時(shí)間的解。

當(dāng)r1≠r2時(shí)，多圈Lambert軌道轉(zhuǎn)移過(guò)程中，最小速度增量與轉(zhuǎn)移時(shí)間關(guān)系的極小值并未呈現(xiàn)出逐漸減小的趨勢(shì)，無(wú)法采用滑行的手段將兩者變換成遞減的關(guān)系，故不能采用處理r1＝r2的方法處理r1≠r2的情況。因此，該條件下軌道轉(zhuǎn)移所需速度增量與轉(zhuǎn)移時(shí)間取兩者關(guān)系的第一個(gè)極小值點(diǎn)。

2.3 優(yōu)化模型

航天器在軌加注服務(wù)系統(tǒng)中服務(wù)航天器的數(shù)量為 k， 序號(hào)集合為 S ＝ ｛s1，s2，…，sk｝，si的最大燃料攜帶量為 Qi，其中 i＝ 1，2，…，k；目標(biāo)航天器的數(shù)量為n，序號(hào)集合為 X ＝ ｛x1，x2，…，xn｝，xj的燃料需求量為 fj，j ＝ 1，2，…，n ；燃料站的數(shù)量為1，其序號(hào)為x0。服務(wù)航天器需要在時(shí)間內(nèi)完成對(duì)所有目標(biāo)航天器的加注作業(yè)。

“k-1-n”在軌加注系統(tǒng)擁有k顆服務(wù)航天器，因此可將整個(gè)在軌加注任務(wù)分成k個(gè)子任務(wù)討論，即第i個(gè)子任務(wù)為si所執(zhí)行的加注任務(wù)。因此，需對(duì)子任務(wù)進(jìn)行優(yōu)化，使在軌加注總?cè)蝿?wù)用于軌道轉(zhuǎn)移的燃耗最小。設(shè)計(jì)優(yōu)化變量Ji、Ai，如式（1）：

式中：Ji為si的加注目標(biāo)和加注順序，如J1＝［1，2，3］，表示s1需要加注的目標(biāo)航天器序號(hào)為 x1、x2、x3， 加注順序?yàn)椋?x1→ x2→ x3；Ai為si執(zhí)行對(duì)應(yīng)目標(biāo)航天器加注作業(yè)所分配的軌道轉(zhuǎn)移時(shí)間區(qū)間索引號(hào)。

服務(wù)航天器si在對(duì)目標(biāo)航天器執(zhí)行加注作業(yè)的過(guò)程中可能出現(xiàn)兩種情況：第一種為si在執(zhí)行對(duì)1的加注作業(yè)后，由于燃料的限制，無(wú)法完成對(duì)的加注作業(yè)，需返回燃料站補(bǔ)充燃料后再機(jī)動(dòng)至為其加注；第二種為si在執(zhí)行對(duì)的加注作業(yè)后，有充足的燃料完成對(duì)的加注作業(yè)無(wú)需返回燃料站。因此，si在完成對(duì)的加注作業(yè)后，需要根據(jù)服務(wù)航天器的燃料剩余量與完成對(duì)執(zhí)行加注作業(yè)燃料需求量的關(guān)系對(duì)的加注作業(yè)進(jìn)行評(píng)估。評(píng)估策略如下：當(dāng)滿(mǎn)足式（2）時(shí)， si在完成的加注作業(yè)后無(wú)需返回燃料站進(jìn)行燃料補(bǔ)充，該情況下設(shè)決策變量＝ 0，si為執(zhí)行加注作業(yè)是從燃料站出發(fā)，故＝0。

即si用于軌道轉(zhuǎn)移的時(shí)間ti＝ ［，…，，tlii＋1］， 其中，tlii＋1為si完成任務(wù)后返回燃料站的軌道轉(zhuǎn)移時(shí)間。完成任務(wù)的總用時(shí)Ti和軌道機(jī)動(dòng)的燃料消耗量分別如式（10）、式（11）所示：

以軌道轉(zhuǎn)移燃耗為優(yōu)化目標(biāo)的航天器在軌加注任務(wù)調(diào)度問(wèn)題考慮三類(lèi)約束條件，第一類(lèi)為加注關(guān)系約束，在加注任務(wù)中一顆目標(biāo)航天器只能被服務(wù)一次，且每顆目標(biāo)航天器都要被服務(wù)。因此，加注關(guān)系約束如式（12）所示：

第二類(lèi)約束為時(shí)間約束，根據(jù)加注任務(wù)的要求，服務(wù)航天器必須在規(guī)定的時(shí)間內(nèi)完成加注任務(wù)，且為了考慮工程實(shí)際，需要對(duì)每次軌道轉(zhuǎn)移時(shí)間進(jìn)行最大值約束［8］，時(shí)間約束如式（13）、式（14）所示：

式（14）中，tmax為允許的最大轉(zhuǎn)移時(shí)間。

第三類(lèi)約束為燃料約束，服務(wù)航天器的燃料需求量需小于服務(wù)航天器的最大燃料攜帶量，燃料約束如式（15）所示：

對(duì)于整個(gè)在軌加注任務(wù)而言，用于軌道機(jī)動(dòng)所消耗的燃料越少越好。所以，本文將通過(guò)優(yōu)化Ji和Ai，在滿(mǎn)足約束條件的前提下，使整個(gè)加注任務(wù)軌道機(jī)動(dòng)的燃料消耗量最小化。在軌加注任務(wù)的數(shù)學(xué)模型如式（16）所示：

3 問(wèn)題求解

基于“k-1-n”模式的在軌加注調(diào)度問(wèn)題屬于整數(shù)規(guī)劃問(wèn)題。該調(diào)度問(wèn)題與車(chē)輛路徑問(wèn)題（Vehicle Routing Problem，VRP）和作業(yè)車(chē)間調(diào)度問(wèn)題（Job-shop Scheduling Problem，JSP）等具有相似性，遺傳算法對(duì)于此類(lèi)問(wèn)題的求解具有較強(qiáng)的優(yōu)勢(shì)。因此，本文將采用遺傳算法求解在軌加注調(diào)度問(wèn)題。

3.1 遺傳算子設(shè)計(jì)

1）編碼

染色體基因分為兩段，順序段基因J和時(shí)間段基因 A， 其中 J ＝ ［0，J1，0，J2，0，…，0，Jk，0］，AA ＝ ［A1，A2，…，Ak］ 。

2）選擇

本文采用隨機(jī)遍歷抽樣法（Stochastic Universal Sampling，SUS）進(jìn)行選擇。對(duì)于染色體適應(yīng)度的計(jì)算取燃料消耗量的倒數(shù)。

3）交叉

染色體包含順序段和時(shí)間段，由于兩段染色體存在差異，所以進(jìn)行交叉操作時(shí)分為順序段和時(shí)間段。

順序段基因的交叉方式為隨機(jī)選取兩個(gè)交叉位，保留雜交位外的基因，再搜索交叉位中間的基因是否為0，如果為0，則保留相應(yīng)基因位上的基因，如果不為0，則交換相應(yīng)基因位的基因。在交叉后對(duì)順序段染色體進(jìn)行基因去重操作。

時(shí)間段基因由于沒(méi)有基因值不一樣的約束，因此交叉方式為隨機(jī)選取兩個(gè)基因位，將基因位中間的基因串交換。

4）變異

順序段基因的變異方式為隨機(jī)選取兩個(gè)需要變異的基因位，將兩個(gè)基因位上的基因互換。

時(shí)間段基因的變異方式為隨機(jī)選取一個(gè)基因位，將該基因位為上的基因換成其它合法值。

3.2 算法流程

采用遺傳算法求解在軌加注任務(wù)調(diào)度問(wèn)題的算法流程如圖4所示。首先確定每個(gè)子任務(wù)的分配方案，即確定每顆服務(wù)航天器的加注目標(biāo)、加注順序以及每次軌道轉(zhuǎn)移所分配的時(shí)間，初始化種群。然后服務(wù)航天器按任務(wù)分配方案執(zhí)行加注作業(yè)，在每次加注作業(yè)的初始時(shí)刻需對(duì)此次加注作業(yè)進(jìn)行評(píng)估，當(dāng)服務(wù)航天器燃料不足時(shí)需返回空間燃料站補(bǔ)充燃料，直至所有目標(biāo)航天器都被加注。最后更新子任務(wù)分配方案，在找到最優(yōu)分配方案后，輸出結(jié)果。

圖4 k-1-n模式在軌加注任務(wù)調(diào)度算法流程圖Fig．4 Algorithm flowchart for k-1-n on-orbit refueling strategy

4 算例

4.1 問(wèn)題配置

選取14顆GEO衛(wèi)星作為目標(biāo)航天器，衛(wèi)星的軌道根數(shù)如表1所示（數(shù)據(jù)參考文獻(xiàn)［5］和文獻(xiàn)［6］）。目標(biāo)航天器所需要的燃料分別為［260，240，230，270，230，200，240，220，270，280，210，230，240，230］，單位為 kg。 燃料站軌道半徑為39 164 km，初始真近點(diǎn)角為0°；服務(wù)航天器的數(shù)量為2，兩顆服務(wù)航天器的干重均為500 kg，最大燃料攜帶量為1500 kg，初始時(shí)刻燃料攜帶量為1500 kg，推進(jìn)系統(tǒng)比沖300 s。設(shè)服務(wù)航天器為單個(gè)目標(biāo)航天器進(jìn)行加注時(shí)加注操作時(shí)間為1 d，從燃料站獲取燃料的操作時(shí)間為1 d，規(guī)定的任務(wù)總時(shí)間約束為40 d，設(shè)置遺傳算法的種群規(guī)模為100，迭代代數(shù)為200，交叉概率為 0.9，變異概率為0.1，代溝為0.9。

4.2 近似最優(yōu)解

根據(jù)給出的條件，遺傳算法經(jīng)過(guò)200次的迭代后（如圖5所示），獲得近似最優(yōu)解。服務(wù)航天器s1用于軌道轉(zhuǎn)移的燃耗為156.1400 kg，完成任務(wù)的時(shí)間為39.1694 d，加注目標(biāo)及加注順序?yàn)椋踴9，x11，x13，x1，x2，x4，x14，x6，x8］， 每次軌道轉(zhuǎn)移燃耗和轉(zhuǎn)移時(shí)間如表2所示。服務(wù)航天器s2用于軌道轉(zhuǎn)移的燃耗為108.9544 kg，完成任務(wù)的時(shí)間為22.5938 d，加注目標(biāo)及加注順序?yàn)?［x3，x10，x12，x5，x7］，每次軌道轉(zhuǎn)移燃耗和轉(zhuǎn)移時(shí)間如表3所示。

表1 GEO航天器的軌道根數(shù)［5-6］Table 1 Orbit elements of GEO spacecraft［5-6］

圖5 遺傳算法優(yōu)化過(guò)程Fig．5 Optimization process of genetic algorithm

表2 服務(wù)航天器s1的軌道轉(zhuǎn)移燃耗和時(shí)間Table 2 Consumption of fuel and time for servicing spacecraft s1

表3 服務(wù)航天器s2的軌道轉(zhuǎn)移燃耗和時(shí)間Table 3 Consumption of fuel and time for servicing spacecraft s2

對(duì)問(wèn)題求解獲得決策變量 α1＝ ［0，0，0，0，0，1，0，0，0］，α2＝ ［0，0，0，0，0］。 所以服務(wù)航天器s1在對(duì)目標(biāo)航天器x4執(zhí)行加注作業(yè)的過(guò)程中需返回空間燃料站補(bǔ)充燃料。服務(wù)航天器在軌加注順序示意圖如圖6所示。

圖6 服務(wù)衛(wèi)星加注示意圖Fig．6 Diagram of servicing spacecraft refueling

4.3 結(jié)果分析

從實(shí)驗(yàn)結(jié)果可知，服務(wù)航天器在返回空間燃料站補(bǔ)充燃料前的最優(yōu)加注順序是按一種真近點(diǎn)角位置依次加注的順序，當(dāng)服務(wù)航天器返回空間燃料站后的加注順序不再完全符合這種規(guī)律。這主要是所有目標(biāo)航天器處于同一圓軌道，服務(wù)航天器按照真近點(diǎn)角位置依次加注，軌道轉(zhuǎn)移的燃耗最小，這一結(jié)果與 Zhou［8］等的結(jié)論一致。但是，空間燃料站與目標(biāo)航天器運(yùn)行在不同高度的共面圓軌道上，服務(wù)航天器在執(zhí)行加注任務(wù)的過(guò)程中需往返燃料站補(bǔ)充燃料，導(dǎo)致空間燃料站與目標(biāo)衛(wèi)星的真近點(diǎn)角差會(huì)發(fā)生變化，改變服務(wù)航天器的下一個(gè)最優(yōu)加注目標(biāo)。并且根據(jù)仿真計(jì)算數(shù)據(jù)可知，不同半徑圓軌道間的軌道轉(zhuǎn)移所需的燃料大于同軌間軌道轉(zhuǎn)移所需的燃料，為了減小服務(wù)航天器往返空間燃料站的燃耗，需要考慮空間燃料站與目標(biāo)航天器的真近點(diǎn)角差，而非只考慮目標(biāo)航天器之間的真近點(diǎn)角差，致使服務(wù)航天器的最優(yōu)加注順序不再符合最優(yōu)加注順序規(guī)律。同時(shí)，軌道轉(zhuǎn)移的燃耗與目標(biāo)航天器的質(zhì)量成正比，當(dāng)某顆目標(biāo)航天器的燃料需求量相對(duì)于其它目標(biāo)足夠大時(shí)，其將會(huì)有較高的優(yōu)先級(jí)。

5 結(jié)論

1）在軌加注調(diào)度問(wèn)題與VRP、JSP等調(diào)度問(wèn)題有著較高的相似度，因此可以借鑒VRP、JSP的求解方法解決在軌加注調(diào)度問(wèn)題。

2）以14顆GEO衛(wèi)星作為加注對(duì)象，運(yùn)用遺傳算法求解結(jié)果證明，利用求解VRP與JSP的方法求解在軌加注作業(yè)調(diào)度問(wèn)題具備可行性。

3）“k-1-n”在軌加注任務(wù)調(diào)度問(wèn)題模型的建立以及問(wèn)題的求解可為后續(xù)解決“k-m-n”在軌加注問(wèn)題提供基礎(chǔ)，甚至可在此基礎(chǔ)上進(jìn)一步分析帶時(shí)間窗口等約束條件的在軌加注問(wèn)題。

參考文獻(xiàn)（References）

［1］ 林巧.GEO航天器在軌燃料補(bǔ)給任務(wù)調(diào)度及優(yōu)化［D］.南京：南京航空航天大學(xué)，2016.Lin Q.Optimal Scheduling of On-orbit Refuelling for GEO Spacecraft［D］.Nanjing：Nanjing University of Aeronautics and Astronautics， 2016.（in Chinese）

［2］ 朱嘯宇，喬兵，張慶展，等.一種基于空間燃料站的GEO航天器在軌加注任務(wù)調(diào)度算法［J］.工程科學(xué)與技術(shù)，2017，49（增刊1）：186-194.Zhu X Y，Qiao B，Zhang Q Z，et al.A reusable on-orbit refueling mode and mission scheduling algorithm for GEO spacecraft through a space fuel station ［J］.Advanced Engineering Sciences， 2017，49（Supp.1）： 186-194.（in Chinese）

［3］ 歐陽(yáng)琦，姚文，陳小前.地球同步軌道衛(wèi)星群在軌加注任務(wù)規(guī)劃［J］. 宇航學(xué)報(bào)，2010，31（12）：2629-2634.Ou Y Q，Yao W，Chen X Q.Mission Programming of on-orbit refueling for geosynchronous satellites ［J］.Journal of Astronautics， 2010， 31（12）： 2629-2634.（in Chinese）

［4］ Tsiotras P， Nailly A.Comparison between peer-to-peer and single spacecraft refueling strategies for spacecraft in circular orbits［C］／／Proceedings of the 2005 Infotech ＠ Aerospace Conference.Arlington：American Institute of Aeronautics and Astronautics， 2005： 26-29.

［5］ Shen H， Tsiotras P.Optimal two-impulse rendezvous using multiple-revolution Lambert solutions ［J］.Journal of Guidance， Control， and Dynamics， 2003，26（1）：50-61.

［6］ Shen H.Optimal Scheduling For Satellite Refuelling In Circular Orbits ［D］.Georgia Institute of Technology，2003.

［7］ Shen H，Tsiotras P.Peer-to-Peer refueling for circular satellite constellation ［J］.Journal of Guidance， Control， and Dynamics， 2005， 28（6）： 1220-1230.

［8］ Zhou H，Yan Y，Huang X，et al.Multi-objective planning of a multiple geostationary spacecraft refuelling mission ［J］.Engineering Optimization， 2017， 49（3）： 531-548.