史人赫，寶音賀西，龍 騰，魏 釗

（1. 清華大學(xué)航天航空學(xué)院，北京100084；2. 北京理工大學(xué)宇航學(xué)院，北京100081；3. 飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京100081）

1 引 言

為了保證空間系統(tǒng)在惡劣的空間環(huán)境中能夠長(zhǎng)期穩(wěn)定的在軌運(yùn)行，在軌服務(wù)應(yīng)運(yùn)而生。在軌服務(wù)是指人或者服務(wù)航天器針對(duì)空間系統(tǒng)進(jìn)行的各類空間操作［1］，主要包括在軌裝配、在軌監(jiān)測(cè)、在軌維護(hù)、空間碎片清理等［2］。在軌服務(wù)為空間任務(wù)的實(shí)施提供了更多選擇，具有延長(zhǎng)在軌系統(tǒng)壽命、提升任務(wù)執(zhí)行能力、降低初始成本、增強(qiáng)在軌靈活性和適應(yīng)性等重要意義［3］。目前，以美國(guó)為首的發(fā)達(dá)國(guó)家高度重視在軌服務(wù)技術(shù)的發(fā)展，開展了軌道快車、鳳凰等在軌服務(wù)技術(shù)驗(yàn)證計(jì)劃［4-5］。作為在軌服務(wù)任務(wù)實(shí)施的前端和頂層，在軌服務(wù)任務(wù)規(guī)劃主要通過(guò)求解多空間目標(biāo)交會(huì)問(wèn)題，在滿足各類復(fù)雜約束條件下，確定服務(wù)航天器對(duì)目標(biāo)航天器的軌道轉(zhuǎn)移與交會(huì)方案，使得在軌服務(wù)任務(wù)效益最大化［6］。如何有效求解復(fù)雜在軌服務(wù)任務(wù)規(guī)劃問(wèn)題，成為了國(guó)內(nèi)外研究熱點(diǎn)。

按照在軌服務(wù)任務(wù)的不同，在軌服務(wù)通?？梢苑譃橐粚?duì)一、一對(duì)多、多對(duì)一和多對(duì)多等不同模式。目前，國(guó)內(nèi)外圍繞不同模式下的在軌服務(wù)任務(wù)規(guī)劃技術(shù)開展了廣泛研究。例如，文獻(xiàn)［6］圍繞不同任務(wù)場(chǎng)景下的在軌服務(wù)任務(wù)規(guī)劃模型開展了系統(tǒng)研究；Yu等［7］考慮服務(wù)窗口、末端狀態(tài)、時(shí)間分配等約束，以總速度增量最小為目標(biāo)，建立了離散連續(xù)混合的空間碎片主動(dòng)清除任務(wù)規(guī)劃模型，并通過(guò)定制的粒子群優(yōu)化算法對(duì)任務(wù)規(guī)劃問(wèn)題進(jìn)行了求解；Yang 等［8］建立了考慮最大回報(bào)值的多空間碎片主動(dòng)清理預(yù)先任務(wù)規(guī)劃模型，提出了一種貪婪啟發(fā)式算法，實(shí)現(xiàn)了該問(wèn)題的高效求解；Madakat 等［9］針對(duì)LEO 空間碎片清理問(wèn)題建立了雙目標(biāo)的旅行商問(wèn)題模型，并通過(guò)分支定界方法完成了求解；梁彥剛等［10］采用多目標(biāo)遺傳算法求解在軌服務(wù)航天器任務(wù)指派問(wèn)題，獲取了速度增量-變軌時(shí)間的Pareto前沿；Baranov［11］等人針對(duì)地球靜止軌道空間碎片清理任務(wù)規(guī)劃問(wèn)題，為服務(wù)航天器設(shè)計(jì)了兩階段的最優(yōu)軌道轉(zhuǎn)移方案。

為了進(jìn)一步提高在軌服務(wù)任務(wù)規(guī)劃求解效率，本文以共面圓軌道衛(wèi)星星群一對(duì)一在軌服務(wù)問(wèn)題為對(duì)象，開展基于網(wǎng)格編碼差分進(jìn)化的在軌服務(wù)任務(wù)規(guī)劃方法研究，在兩層任務(wù)規(guī)劃框架的基礎(chǔ)上，通過(guò)定制一種新型網(wǎng)格編碼生成技術(shù)與差分變異操作機(jī)制，實(shí)現(xiàn)共面圓軌道衛(wèi)星星群一對(duì)一在軌服務(wù)規(guī)劃問(wèn)題的有效求解。

2 在軌服務(wù)任務(wù)規(guī)劃模型

2.1 一對(duì)一在軌服務(wù)任務(wù)規(guī)劃框架

本文研究的在軌服務(wù)任務(wù)規(guī)劃問(wèn)題采用一對(duì)一工作模式，即一顆服務(wù)星對(duì)一顆目標(biāo)星進(jìn)行服務(wù)，且服務(wù)衛(wèi)星星群和目標(biāo)衛(wèi)星星群位于共面圓軌道，各目標(biāo)星的優(yōu)先級(jí)相同。為了實(shí)現(xiàn)在軌服務(wù)任務(wù)規(guī)劃問(wèn)題的高效求解，本文構(gòu)建了如圖1 所示的兩層任務(wù)規(guī)劃框架。其中，任務(wù)分配層通過(guò)定制的網(wǎng)絡(luò)編碼差分進(jìn)化算法，生成滿足一對(duì)一約束的服務(wù)星任務(wù)指派方案；規(guī)劃求解層中，服務(wù)星群按照任務(wù)分配結(jié)果通過(guò)雙脈沖霍曼轉(zhuǎn)移實(shí)現(xiàn)和對(duì)應(yīng)目標(biāo)星的空間交會(huì)，并計(jì)算軌道轉(zhuǎn)移時(shí)間、速度增量等指標(biāo)作為任務(wù)規(guī)劃目標(biāo)函數(shù)；最終根據(jù)優(yōu)化收斂條件，判斷任務(wù)規(guī)劃求解過(guò)程是否終止。

圖1 在軌服務(wù)任務(wù)規(guī)劃框架Fig.1 On-orbit servicing mission planning framework

基于上述在軌服務(wù)任務(wù)規(guī)劃框架，建立如式（1）所示的一對(duì)一在軌服務(wù)任務(wù)規(guī)劃優(yōu)化問(wèn)題的數(shù)學(xué)模型：

其中，X為指派變量矩陣，其元素xij=1 時(shí)表示第i個(gè)服務(wù)星對(duì)第j個(gè)目標(biāo)星進(jìn)行服務(wù)，xij=0 表示不存在服務(wù)關(guān)系；ti和Δvi分別為第i個(gè)服務(wù)星軌道轉(zhuǎn)移時(shí)間與變軌所需速度增量，由空間交會(huì)模型計(jì)算得到；w1和w2分別為在軌服務(wù)任務(wù)時(shí)間和能量指標(biāo)的加權(quán)系數(shù)；t0和v0為歸一化因子，用以消除不同指標(biāo)量級(jí)造成的差別；約束條件g1表示每個(gè)目標(biāo)星必須且只能接受一個(gè)服務(wù)星；約束條件g2表示每個(gè)服務(wù)星最多服務(wù)一個(gè)目標(biāo)星。

2.2 雙脈沖霍曼轉(zhuǎn)移模型

本節(jié)對(duì)在軌服務(wù)任務(wù)規(guī)劃框架采用的雙脈沖霍曼轉(zhuǎn)移模型進(jìn)行簡(jiǎn)要介紹。雙脈沖霍曼轉(zhuǎn)移過(guò)程如圖2 所示，服務(wù)星在軌道半徑為r1的圓軌道C1的任意點(diǎn)P 處產(chǎn)生第一個(gè)速度增量Δv1，進(jìn)入橢圓轉(zhuǎn)移軌道E，并在E 的遠(yuǎn)地點(diǎn)A 產(chǎn)生第二個(gè)速度增量Δv2，進(jìn)入半徑為r2的圓軌道C2。

圖2 霍曼轉(zhuǎn)移過(guò)程示意圖Fig.2 Illustration of Hohmann transfer

霍曼轉(zhuǎn)移過(guò)程兩次脈沖所需的速度增量如式（2）所示：

其中，vC1和vC2分別為半徑r1和r2圓軌道的速度。對(duì)應(yīng)的霍曼轉(zhuǎn)移變軌時(shí)間為：

其中，μe= 398600.5 km3/s2為地球常數(shù)。

此外，為了實(shí)現(xiàn)服務(wù)星與目標(biāo)星的交會(huì)，在軌道轉(zhuǎn)移前目標(biāo)星需要超前服務(wù)星一個(gè)提前角θH，因此服務(wù)星需要根據(jù)和目標(biāo)星的相對(duì)相位在停泊軌道上等待一段時(shí)間。因此，服務(wù)星的實(shí)際軌道轉(zhuǎn)移時(shí)間為等待時(shí)間和變軌時(shí)間之和。文獻(xiàn)［12］中給出了雙脈沖霍曼轉(zhuǎn)移時(shí)間和速度增量的詳細(xì)計(jì)算方法。

3 網(wǎng)格編碼差分進(jìn)化算法

3.1 總體流程

為了高效求解式（1）中的在軌服務(wù)任務(wù)規(guī)劃問(wèn)題，本研究在標(biāo)準(zhǔn)差分進(jìn)化算法的基礎(chǔ)上［13］，提出了一種網(wǎng)格編碼差分進(jìn)化算法（Grid Coding based Differential Evolution，GCDE）。該方法通過(guò)一種隨機(jī)序列逐次枚舉的指派矩陣生成算法，保證每一代個(gè)體自動(dòng)滿足規(guī)劃約束條件，并通過(guò)定制網(wǎng)格編碼差分變異操作，實(shí)現(xiàn)設(shè)計(jì)空間的高效探索。GCDE的具體步驟敘述如下。

步驟1：初始化。利用NP個(gè)維度為n×m的0-1任務(wù)指派矩陣作為每一代的種群Pk，G，種群中的每個(gè)個(gè)體（即任務(wù)指派矩陣）表示為：

其中，k為個(gè)體編號(hào)，G為進(jìn)化的代數(shù)，Xk，G分別為第k個(gè)個(gè)體對(duì)應(yīng)的任務(wù)指派矩陣。3.2 節(jié)中給出了Xk，G的具體生成方式。

步驟2：網(wǎng)格編碼差分變異。為了產(chǎn)生變異向量，算法從當(dāng)前種群中隨機(jī)選擇三個(gè)個(gè)體Xr1，G，Xr2，G，Xr3，G，在此基礎(chǔ)上，采用式（5）所示的差分變異方法產(chǎn)生變異向量。

式（5）中，vk，G為受體向量，⊕和?為網(wǎng)格編碼差分變異操作算子。其中，?表示廣義差分運(yùn)算，用于計(jì)算Xr2，G和Xr3，G的差異性；⊕為廣義變異運(yùn)算，用于根據(jù)差分結(jié)果實(shí)現(xiàn)Xr1，G的擾動(dòng)。3.3 節(jié)中給出了網(wǎng)格編碼差分變異操作的詳細(xì)介紹。

步驟3：選擇。按照式（6）進(jìn)行選擇，確定vk，G是否作為下一代的個(gè)體成員。若vk，G的任務(wù)規(guī)劃目標(biāo)函數(shù)值不劣于當(dāng)前個(gè)體Xk，G，則將其替代Xk，G作為新的種群個(gè)體參與下一輪迭代；反之則不會(huì)保留vk，G。

步驟4：對(duì)當(dāng)前種群中的每一個(gè)個(gè)體重復(fù)上述優(yōu)化過(guò)程，直至達(dá)到模型最大調(diào)用次數(shù)或最大進(jìn)化代數(shù)，最終輸出最優(yōu)任務(wù)規(guī)劃結(jié)果。

3.2 基于序列逐次枚舉的網(wǎng)格編碼方法

為了減少計(jì)算成本，保證任務(wù)規(guī)劃結(jié)果的可行性，本研究引入隨機(jī)序列逐次枚舉思想［14］，實(shí)現(xiàn)Xk，G的快速可靠生成。圖3中給出了Xk，G生成過(guò)程的示意圖，具體步驟包括：

步驟1：生成n×m的二維平面網(wǎng)格，令網(wǎng)格列編號(hào)q=1，可用行編號(hào)r=[1，2，...，n]。

步驟2：在第q列網(wǎng)格中，隨機(jī)從r中選取一個(gè)元素r(t)，為坐標(biāo)為[r(t)，q]的網(wǎng)格賦值為1，為當(dāng)前列其他網(wǎng)格賦值為0，并從r中刪除r(t)。

步驟3：若q≠m，則q=q+1 返回步驟1.1；否則算法終止，輸出指派矩陣個(gè)體Xk，G。

可以看出，按照上述方法生成的任務(wù)指派矩陣自動(dòng)滿足在軌服務(wù)任務(wù)規(guī)劃問(wèn)題中的一對(duì)一約束條件，且交換任務(wù)指派矩陣網(wǎng)格編碼中任意兩列或兩行，不改變?nèi)蝿?wù)規(guī)劃結(jié)果的可行性，從而有效提升了在軌服務(wù)任務(wù)規(guī)劃約束優(yōu)化問(wèn)題的求解效率和魯棒性。

圖3 網(wǎng)格編碼生成過(guò)程示意圖Fig.3 Illustration of grid coding generation

3.3 網(wǎng)格編碼差分變異操作

作為差分進(jìn)化算法的核心，差分變異通過(guò)尋找兩個(gè)隨機(jī)個(gè)體之間的差別對(duì)目標(biāo)個(gè)體進(jìn)行擾動(dòng)，從而實(shí)現(xiàn)設(shè)計(jì)空間探索［15］。與傳統(tǒng)連續(xù)優(yōu)化問(wèn)題或者二進(jìn)制離散編碼形式不同，GCDE 采用平面二維網(wǎng)格編碼方式構(gòu)建子代個(gè)體，因此無(wú)法直接通過(guò)傳統(tǒng)差分操作產(chǎn)生變異個(gè)體向量。因此，本研究根據(jù)網(wǎng)格編碼交換任意兩列不破壞一對(duì)一任務(wù)約束的特點(diǎn)，定制了一種新型網(wǎng)格編碼差分變異操作，其示意圖如圖4所示，具體步驟敘述如下。

圖4 網(wǎng)格編碼差分變異操作示意圖Fig.4 Illustration of grid coding differential mutation operation

步驟1：對(duì)于隨機(jī)選擇的個(gè)體Xr2，G和Xr3，G，分別比較其網(wǎng)格編碼，找到任務(wù)分配結(jié)果不同的列并記錄為集合：

其中，dz為Xr2，G和Xr3，G中第z個(gè)不同編碼列在網(wǎng)格編碼矩陣中的索引序號(hào)。

步驟2：生成0～1 之間的隨機(jī)數(shù)rand[0，1]，并與差分變異因子F進(jìn)行比較。若rand[0，1]≤F，則將Qd中的索引序號(hào)順序隨機(jī)排列，重新生成一組索引F越大則表示當(dāng)前個(gè)體具有更大的概率實(shí)現(xiàn)差分變異。

步驟3：在Xr1，G中，根據(jù)原始索引Qd找到對(duì)應(yīng)的需要變異的列，并根據(jù)新的索引向量Q?d交換列的順序，生成變異后的受體向量vi，G。

4 工程案例

本節(jié)通過(guò)某衛(wèi)星星群一對(duì)一在軌服務(wù)任務(wù)規(guī)劃案例，對(duì)本文提出的GCDE方法的有效性和實(shí)用性進(jìn)行驗(yàn)證。本案例中，12顆服務(wù)星均勻分布在半長(zhǎng)軸為8378.14 km、軌道傾角為60°的圓形軌道上，8顆目標(biāo)星均勻分布在同一平面內(nèi)半長(zhǎng)軸為8678.14 km的圓形軌道上，其初始軌道根數(shù)如表1 所示。采用GCDE 方法求解本案例中的在軌服務(wù)任務(wù)規(guī)劃問(wèn)題。其中，時(shí)間和能量指標(biāo)的加權(quán)系數(shù)w1和w2均取0.5，歸一化因子取t0=2592000 s 和v0=1000 m/s，GCDE 種群規(guī)模為96，最大進(jìn)化代數(shù)為10，差分變異因子F=0.8，得到優(yōu)化后的任務(wù)規(guī)劃結(jié)果如表2所示。

從結(jié)果中可以看出，服務(wù)星群總體上優(yōu)先選擇軌道運(yùn)行方向上相位差最小的臨近目標(biāo)星進(jìn)行服務(wù)，從而同時(shí)降低在軌服務(wù)任務(wù)完成時(shí)間和燃料消耗。此外，優(yōu)化后的任務(wù)規(guī)劃方案滿足一對(duì)一在軌服務(wù)模式下的任務(wù)指派約束，且所需的最大速度增量為122.0 m/s、最長(zhǎng)變軌時(shí)間為40225.0 s，符合工程實(shí)際要求。為了進(jìn)一步說(shuō)明GCDE 的優(yōu)化性能，本研究還采用MATLAB 提供的整數(shù)編碼遺傳算法（GA）工具箱求解上述在軌服務(wù)任務(wù)規(guī)劃問(wèn)題（GA種群規(guī)模為100，最大進(jìn)化代數(shù)為100，其它參數(shù)采用MATLAB GA 工具箱默認(rèn)設(shè)置），得到GCDE 與GA 優(yōu)化結(jié)果對(duì)比如表3所示。表3中的結(jié)果表明，本文提出的GCDE 方法在求解衛(wèi)星星群一對(duì)一在軌服務(wù)任務(wù)規(guī)劃問(wèn)題時(shí)具有更好的全局收斂性，且其計(jì)算成本相比于傳統(tǒng)GA 減少了約90%，從而驗(yàn)證了GCDE的有效性和實(shí)用性。

表1 服務(wù)星和目標(biāo)星初始軌道根數(shù)Table 1 Initial orbital elements of the service satellites and target satellites

表2 任務(wù)規(guī)劃結(jié)果Table 2 Results of mission planning

表3 GCDE和GA優(yōu)化結(jié)果對(duì)比Table 3 Comparison of GCDE and GA optimization results

5 結(jié) 論

本文提出了一種基于網(wǎng)格編碼差分進(jìn)化的在軌服務(wù)任務(wù)規(guī)劃方法，能夠?qū)崿F(xiàn)共面圓軌道衛(wèi)星星群一對(duì)一任務(wù)規(guī)劃問(wèn)題的有效求解。本研究采用任務(wù)分配-霍曼轉(zhuǎn)移兩層結(jié)構(gòu)建立了規(guī)劃框架，并提出了一種新型網(wǎng)格編碼差分進(jìn)化算法，通過(guò)定制網(wǎng)格編碼機(jī)制與差分變異操作，實(shí)現(xiàn)最優(yōu)可行任務(wù)規(guī)劃方案的高效求解。仿真案例表明，本文提出的網(wǎng)格編碼差分進(jìn)化方法能夠?qū)崿F(xiàn)共面圓軌道衛(wèi)星星群一對(duì)一在軌服務(wù)任務(wù)規(guī)劃問(wèn)題的有效求解，且計(jì)算效率較傳統(tǒng)遺傳算法具有顯著提升，從而驗(yàn)證了本文研究工作的有效性與工程實(shí)用性。