李志亮, 李小將, 張東來
(1. 航天工程大學研究生院, 北京 101416; 2. 航天工程大學宇航科學與技術系, 北京 101416; 3. 酒泉衛星發射中心, 甘肅 酒泉 732750)
敏捷成像衛星與傳統成像衛星相比具有滾動、俯仰和偏航姿態機動性能,對地觀測能力更強,任務沖突的解決方式更多,已成為新型衛星發展的重要方向[1]。有效的規劃和調度是敏捷成像衛星發揮其優勢的重要支撐,也是國內外學者研究的熱點問題。目前,多數研究[1-4]在確定條件下展開,即假設與調度相關的信息都是確定的。然而不確定因素在敏捷成像衛星任務調度中時刻存在,并且影響調度方案的執行。文獻[5]總結了4種典型的不確定因素,包括新任務的到達、目標觀測和數傳機會的變化、衛星資源狀態的不確定性、任務屬性的變化。為應對不確定因素的影響,許多學者研究了前攝式調度方法,前攝式調度[6]是指在調度之初預先考慮不確定因素,將其作為先驗知識引進調度中,生成適應性或魯棒性調度方案。有效地表示不確定因素并構建魯棒性指標是建立前攝式調度魯棒模型的關鍵。文獻[7]基于隨機云層預測模型研究了EO-1衛星的魯棒性規劃問題。文獻[8]針對云層覆蓋的不確定因素,提出一種基于鄰域的魯棒性指標,建立了約束滿足模型并采用基于分級優化策略的隨機變鄰域禁忌搜索算法進行求解。文獻[9]在文獻[8]基礎上,增加了總任務執行時間指標,建立了多目標魯棒調度模型,采用NSGA-Ⅱ算法求解。文獻[10]將總時間重疊度和總任務執行時間作為魯棒性指標,建立了多目標魯棒性模型,同樣采用NSGA-Ⅱ算法求解。文獻[11]考慮將任務調度到多個可見窗口以增加觀測成功的概率,構建了魯棒模型。文獻[12]考慮任務間衛星姿態轉換約束和重要時間窗預留規則,提出一種基于資源預留的魯棒性任務規劃方法。
上述研究中,文獻[7-12]均沒有考慮衛星資源失效對調度的影響,不能滿足衛星資源狀態不確定條件下調度系統具有較強魯棒性的需求。此外,觀測任務的選擇、排序以及觀測時間的確定增加敏捷成像衛星前攝式調度求解的復雜性,傳統求解算法[9-10]難以適用。
針對以往研究的不足,本文考慮衛星資源失效、應急任務加入兩類不確定因素,構建了期望收益和松弛時間兩個魯棒性指標,將這兩類指標作為優化目標構建了前攝式調度魯棒模型,針對該模型的多目標優化特性,設計了一種多目標離散差分進化求解算法。并通過仿真計算驗證了算法的優越性。
敏捷成像衛星觀測任務的執行伴隨著兩類主要的不確定因素:衛星資源失效和應急任務加入。前攝式調度是應對不確定因素的重要方法,敏捷成像衛星前攝式調度是指將衛星資源失效概率、任務執行主從窗口等信息作為先驗知識,構建兼顧收益和魯棒性的優化目標,同時結合衛星資源和任務執行約束條件,將有限的衛星資源分配給多個觀測任務,并確定任務執行時間,生成魯棒性調度方案。本文首先構建期望收益、松弛時間兩個魯棒性指標,進而建立前攝式調度魯棒模型。
本文考慮的魯棒性指標包括期望收益和松弛時間。
1.1.1 期望收益
衛星資源失效具有很強的隨機性,其發生難以預測。文獻[13]假設某一時刻至多一顆衛星失效,且失效事件相互獨立,但是沒有給出失效概率計算方法。文獻[14]假設機器故障概率與工作時間呈正相關,對機器故障概率進行了描述。借鑒該描述,本文增加空閑時間失效概率,假設衛星資源失效滿足指數分布,采用下式計算衛星失效概率:
(1)
式中,pn(t)表示t時刻衛星sn失效的概率;rth表示衛星sn第h次失效后維修結束時間,且rth≥0,調度初始時刻,記rth=0;FTl表示衛星sn在[rth,t]內的第l個空閑時間段;tsFTl和teFTl表示該時間段的起止時間;NFT表示[rth,t]內空閑時間段的數量;AWTn表示衛星sn的平均正常工作時間;ARTn表示衛星sn的平均維修時間。

(2)


考慮了衛星資源失效的概率,以及任務執行的潛在可用從窗口,本文采用下式表示任務ti的期望收益:
(3)

1.1.2 松弛時間
研究表明調度方案中預留松弛時間能夠吸收任務變化帶來的擾動,增強調度的魯棒性[15]。借鑒文獻[15]對松弛時間的定義,本文將敏捷成像衛星前攝式調度中的松弛時間定義為:設任務ti、tj為相鄰任務,松弛時間是指從任務ti到tj的姿態轉換結束到任務tj開始觀測的時間段,描述如下式:
slackij=otj-oti-di-trij-tst
(4)
式中,oti、otj分別表示ti、tj的開始觀測時間;di表示ti的觀測時長;trij表示衛星從ti到tj的姿態機動時間;tst表示衛星姿態機動穩定時間。
本文在敏捷成像衛星調度領域經典文獻[1]的基礎上進行合理擴展,考慮多顆衛星多個軌道圈次、任務間姿態轉換、以及星上存儲和能量等現實條件,從以下4個方面描述模型。
1.2.1 參數定義
本文模型相關參數定義如表1所示。

表1 相關參數定義
1.2.2 決策變量
本文采用的決策變量包括兩個:任務ti是否被觀測,任務ti和tj是否為相鄰任務,描述如下式:
(6)
1.2.3 目標函數
根據文中提出的期望收益和松弛時間指標,敏捷成像衛星前攝式調度魯棒模型的優化目標為最大化調度期望收益和總松弛時間,其表達式為
(7)
(8)
式中,F1表示最大化調度期望收益;F2表示最大化總松弛時間。
1.2.4 約束條件
考慮任務間姿態轉換約束、星上存儲和能量約束等,本文模型包括以下5方面的約束條件:
(1) 任務執行的唯一性約束

(9)

(10)
(2) 任務執行的時間窗口約束

(11)
(3) 相鄰任務間姿態轉換時間約束
?ti,tj∈T,i≠j:
(12)
(4) 衛星單個軌道圈次的存儲約束

(13)
(5) 衛星單個軌道圈次的能量約束
?sn∈S,?k∈OBn,?ti,tj∈T,i≠j:
(14)
上述模型將調度期望收益和總松弛時間作為優化目標,同時考慮了衛星資源和觀測任務的多重復雜約束,是較為復雜的多目標優化模型,本文采用多目標離散差分進化算法進行求解。
敏捷成像衛星前攝式調度魯棒模型的求解是一個多目標優化問題。NSGA-Ⅱ算法在多目標優化上取得了不錯的效果[9,16],同時,NSGA-Ⅱ算法在尋優中容易陷入局部最優,而且基于擁擠距離的多樣性評價準則沒有考慮目標函數值階數的差別,容易出現多個解聚集的情況。本文考慮差分進化的全局探索優勢[17]和擾動操作的局部搜索能力,從離散差分進化、外部存檔更新策略、Pareto解集的評價方面設計MDDE求解算法,并給出了算法的實現步驟。
離散差分進化基本思路是首先采用離散形式對解進行編碼,然后生成算法的初始種群,進而采用變異、交叉、選擇算子得到下一代個體,通過競爭獲得解空間內更好的解。
解的編碼和初始種群生成是離散差分進化的基礎。對于解的編碼,本文根據任務執行順序采用自然數編碼方式,建立衛星與任務的對應關系,例如,x=[3 5 1 0 2 4 9 0 6 8 7],其中數字0用以區分不同的衛星。對于初始種群的生成,本文在貪婪算法[1]的基礎上,增加兩個隨機因素:隨機選擇衛星軌道圈次、輪盤賭法選擇候選任務,采用隨機貪婪算法(random greedy algorithm,RGA)生成規模為NP的初始種群。
2.1.1 變異算子

(15)

對于變異個體的產生,DesCons[18]、InsertMut[19]等操作取得了不錯的尋優效果,但是這些操作所針對的調度問題沒有任務執行的時間窗口約束,不能被本文直接采用。因此,本文提出適用于敏捷成像衛星調度問題的DC操作和F操作。
DC操作的基本思路是首先從每顆衛星的任務序列中隨機刪除比例為dps的任務,dps∈(0,1),然后從未安排的任務集合中隨機選擇比例為cps的任務,將選擇的任務按照調度約束條件插入到現有序列中,得到變異個體。
F操作的基本思路是判斷相鄰任務間的松弛時間是否大于閾值MinD·slc,MinD表示所有任務的最小觀測時長,slc表示擾動強度,slc∈[1,10],如果小于MinD·slc,則刪除相鄰任務中wi/di較小的任務,得到變異個體。
2.1.2 交叉算子
根據離散差分進化的機理,試驗個體可采用下式獲得[19]:
(16)

基于塊結構的整體移動策略[19]是一種有效的試驗個體產生操作。在敏捷成像衛星調度問題中,由于時間窗口約束,任務序列同樣具有塊結構特征,但是塊結構的整體移動容易產生不可行解。為此,本文提出一種考慮局部修復的CR操作,基本思路是對于每顆衛星的任務序列,隨機選擇一個切點,計算切點處任務的觀測結束時間,將滿足姿態轉換的切點右側任務序列進行交換,然后刪除交換后任務序列中重復的任務,得到試驗個體。
2.1.3 選擇算子

(17)

為保留算法進化中優質解,本文采用外部存檔保存搜索到的非支配解,設外部存檔為EA,EA={ea1,ea2,…,eaNE},NE表示EA的規模,每進化一代,EA根據當前新解y進行更新。同時,為避免EA中非支配解過于集中,引進個體間的歸一化鄰域距離,當EA的規模大于NP時,根據歸一化鄰域距離刪除擁擠的解,以保持解的多樣性。
對于兩個目標函數的多目標優化問題,設兩個目標函數分別為f1(·)和f2(·),根據Pareto解集的定義,EA中任意兩個解互相不構成支配關系,即對于?eak,eak+1∈EA,若f1(eak)

圖1 新解y與外部存檔EA的相對關系Fig.1 Relative relationship between new solution y and external archive EA
由圖1可知,EA的更新包括4種情況:①當y與EA的每個解均互不支配時,將y加入EA;②當y支配EA的部分解時,將被支配的解從EA刪除,并將y加入EA;③當y支配EA的全部解時,將EA的所有解刪除,并將y加入EA;④當y被EA的任一解支配時,則不更新EA。
為保持EA中非支配解的多樣性,需量化評價解之間的擁擠程度,文獻[16]提出一種擁擠距離計算公式,該公式具有簡單高效的計算特點,但是沒有考慮目標函數值階數的差別,容易出現多個解聚集的情況。為此,本文在計算相鄰解的距離時對目標函數值進行歸一化處理,稱之為歸一化鄰域距離,其計算公式為
(18)

當外部存檔超過一定規模時,需根據歸一化鄰域距離計算出外部存檔中最擁擠的解,將EA中所有解的歸一化鄰域距離記為集合ND,ND={nd1,nd2,…,ndNE-1},設ND中最小的距離為ndk,比較ndk左側和右側的鄰域,若ndk-1
MDDE算法的求解過程就是不斷構造Pareto解集,并向解空間的最優前沿靠近的過程。因此,對多目標差分進化算法性能評價的重要依據是Pareto解集的優劣。不同于單目標優化,多目標優化很難用單一的指標來衡量,通常會考慮3方面的性能指標[20]:收斂性、均勻性和多樣性。針對這3方面的性能指標,本文結合敏捷成像衛星前攝式調度的特點,對文獻[20]提出的指標進行了改進,采用以下3個量化指標對Pareto解集進行評價,這里記待評價的Pareto解集為SP,其規模為NS,初始種群的非劣解集為InitS。
2.3.1 世代距離
世代距離[20](generational distance,GD)用于評價求解所得Pareto解集與真實Pareto最優解集的接近程度,該指標需要參考Pareto最優解集才能確定,而現實中調度問題解空間的最優解集難以獲得。本文將初始種群的非劣解集作為參考集合以評價Pareto解集的收斂性,GD計算公式為
(19)
(20)
式中,rdk為非支配解xk(xk∈SP)到InitS的最小歐式距離;GD表示SP相對InitS的進化距離,與文獻[20]不同,本文中GD越大,SP越接近Pareto最優前沿;當GD=0時,SP即為InitS。
2.3.2 分散間距
分散間距(spacing,SP)[20]用于評價Pareto解集中各個解分布的均勻性和多樣性,與文獻[20]不同的是,本文考慮目標函數值階數的差別,采用歸一化鄰域距離計算相鄰解之間的距離,SP計算公式為
(21)

2.3.3 最大分布廣度
最大分布廣度(maximum spread,MS)[20]用于衡量Pareto解集的散布范圍,本文在MS的計算時將初始解集的非劣解集作為參考集合,計算公式為
(22)

綜合上述的分析和設計,MDDE算法實現的具體步驟如下:
步驟1初始化算法的相關參數:種群規模NP、最大迭代次數Itermax、變異概率Pm、交叉概率Pc;采用隨機貪婪算法生成初始種群。
步驟2對初始種群進行基于Pareto的非支配排序,將排序后初始種群的最優前沿作為外部存檔EA。


步驟5判斷是否達到最大迭代次數Itermax,如果達到,則輸出外部存檔EA作為Pareto非劣解集,并對解集的收斂性、均勻性、多樣性進行評價,如果沒有達到,則重復步驟步驟3和步驟4。
首先設置敏捷成像衛星、成像任務和MDDE算法的相關參數,然后采用Matlab R2010b實現MDDE算法,并對比分析MDDE與NSGA-Ⅱ[9,16]算法的求解結果。
仿真算例中敏捷成像衛星的主要參數參考Pleiades衛星[1],成像任務采取隨機生成的方式獲得,考慮到點目標是成像的基本單元,條帶目標可以看作具有一定觀測時長的點目標,區域目標通過分解可由點目標和條帶目標組成,本文選擇在一定范圍內隨機生成不同規模的點目標任務集[21]。衛星和任務的主要參數如表2所示。

表2 衛星和任務的主要參數
表2中,rand([x,y])表示x和y之間服從均勻分布的隨機數。通過組合產生5組算例(衛星數~任務數):2~50、3~100、4~200、5~300、6~600。
算法參數設置如下:種群規模NP=100,進化代數Itermax=100,變異概率Pm=0.7,交叉概率Pc=0.5,解構的擾動強度dps=0.4,重構的擾動強度cps=0.6,插入松弛時間的擾動強度slc=2。
采用Matlab實現MDDE算法,計算機配置為Intel Core i7-3770 CPU @3.4GHz、8核4GB內存,操作系統為Windows 7、編程環境為Matlab R2010b。為驗證MDDE的優越性,將MDDE與NSGA-Ⅱ算法[9,16]在5組算例上進行對比,以世代距離、分散間距、最大分布廣度作為算法性能對比的指標,需要說明的是,這3個指標均無量綱。為避免隨機因素的影響,每種算法在每組算例上獨立運行20次并取平均值,對比如表3所示。

表3 MDDE與NSGA-Ⅱ算法的性能對比
由表3可知,在GD指標上,MDDE比NSGA-Ⅱ高14.32%左右,說明MDDE算法求解得到的Pareto解集更接近真實最優前沿;在SP指標上,MDDE比NSGA-Ⅱ低5.17%左右,說明MDDE算法求解得到的Pareto解集中各個解的分布更均勻;在MS指標上,MDDE比NSGA-Ⅱ高11.79%左右,說明MDDE算法求解得到的Pareto解集的最大分布范圍更廣。
以算例4~200為例,MDDE與NSGA-Ⅱ算法的Pareto求解結果對比如圖2所示。圖2中,RGA表示采用隨機貪婪算法求解得到的初始種群非劣解集,NSGA-Ⅱ-100表示采用NSGA-Ⅱ算法迭代100次求解得到的Pareto解集,MDDE-1、MDDE-50、MDDE-100分別表示采用MDDE算法迭代1、50、100次求解得到的Pareto解集。

圖2 MDDE與NSGA-Ⅱ算法Pareto求解結果對比Fig.2 Pareto solution comparison of MDDE with NSGA-Ⅱ algorithm
由圖2可知,MDDE算法從初始種群開始,通過多次迭代,不斷向解空間最優前沿靠近,且MDDE算法經過100次迭代求解所得解集優于NSGA-Ⅱ算法經過100次迭代求解所得解集。
在運行時間方面,MDDE與NSGA-Ⅱ算法在5個算例上的結果對比如圖3所示。由圖3可知,MDDE的運行時間在5個算例上比NSGA-Ⅱ平均降低了9.72%左右,說明MDDE比NSGA-Ⅱ的求解效率更高。同時可以看出,從算例2~50到6~600,隨著算例規模的變大,算法的運行時間大幅增長,這是由于敏捷成像衛星前攝式調度問題具有NP-Hard特性,解空間隨問題規模呈指數級增大,算法運行時間也大幅增長。

圖3 MDDE與NSGA-Ⅱ算法運行時間對比Fig.3 Average CPU time comparison of MDDE and NSGA-Ⅱ algorithm
綜上結果分析,MDDE算法在求解結果和求解效率上均優于NSGA-Ⅱ算法,驗證了MDDE算法的優越性。
本文針對敏捷成像衛星前攝式調度中存在的衛星資源失效、應急任務加入兩類不確定因素,構建了以期望收益和松弛時間為優化目標的調度魯棒模型,針對該模型的多目標優化特性,提出了一種MDDE求解算法,仿真結果表明3個Pareto解集評價指標上MDDE比NSGA-Ⅱ提高了10.42%左右,在運行時間上降低了9.72%左右,驗證了MDDE算法的優越性。同時,本文研究沒有考慮擾動事件發生后如何對方案進行調整,在后續工作中,可進一步研究反應式調度模型和算法,以期采用較小幅度的方案調整應對多種擾動事件的影響。
[2] HABET D, VASQUEZ M, VIMONT Y. Bounding the optimum for the problem of scheduling the photographs of an agile earth observing satellite[J]. Computational Optimization and Applications, 2010, 47(2): 307-333.
[3] 姜維, 龐秀麗, 郝會成. 成像衛星協同任務規劃模型與算法[J]. 系統工程與電技術, 2013, 35(10): 2093-2101.
JIANG W, PANG X L, HAO H C. Collaborative scheduling model and algorithm for imaging satellite network[J]. Systems Engineering and Electronics, 2013, 35(10): 2093-2101.
[4] XU R, CHEN H, LIANG X, et al. Priority-based constructive algorithms for scheduling agile earth observation satellites with total priority maximization[J]. Expert Systems with Applications, 2016, 51(C): 195-206.
[5] PEMBERTON J C, GREENWALD L G. On the need for dynamic scheduling of imaging satellites[J]. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences, 2002, 34(1): 165-171.
[6] LAMBRECHTS O, DEMEULEMEESTER E, HERROELEN W. Proactive and reactive strategies for resource-constrained project scheduling with uncertain resource availabilities[J]. Journal of Scheduling, 2008, 11(2): 121-136.
[7] ABRAMSON M, CARR F, CARR D, et al. Robust planning for the Earth Observing-1 (EO-1) mission[C]∥Proc.of the AIAA Aerospace Conference, 2013: 33-36.
[8] 王軍民,李菊芳,譚躍進.不確定條件下衛星魯棒性調度問題[J].系統工程,2007,25(12): 94-99.
WANG J M, LI J F, TAN Y J. The robust satellite scheduling problem under uncertainty[J].Systems Engineering,2007,25(12): 94-99.
[9] NIU X, TANG H, WU L, et al. Imaging-duration embedded dynamic scheduling of earth observation satellites for emergent events[J].Mathematical Problems in Engineering,2015,4:1-31.
[10] 鄧潤, 唐宏, 單越, 等. 面向應急任務衛星魯棒性規劃模型及算法[J]. 遙感信息, 2014, 29(5): 25-31.
DENG R, TANG H, SHAN Y, et al. Emergency task oriented satellite robust mission planning model and algorithm[J]. Remote Sensing Information, 2014, 29(5): 25-31.
[11] WANG J, DEMEULEMEESTER E, QIU D. A pure proactive scheduling algorithm for multiple earth observation satellites under uncertainties of clouds[J]. Computers & Operations Research, 2016,74(C):1-13.
[12] 張忠山, 譚躍進, 義余江, 等. 基于資源預留的成像衛星魯棒性任務規劃方法[J]. 系統工程理論與實踐, 2016, 36(6): 1544-1554.
ZHANG Z S, TAN Y J, YI Y J, et al. A robustness task planning method for imaging satellite based on resources reservation[J]. Systems Engineering Theory & Practice, 2016, 36(6): 1544-1554.
[13] ZHU X, WANG J, QIN X, et al. Fault-tolerant scheduling for real-time tasks on multiple earth observation satellites[J]. IEEE Trans.on Parallel & Distributed Systems, 2014, 26(1): 1-15.
[14] 何偉. 機器故障下柔性Job Shop調度研究[D]. 重慶: 重慶大學, 2012.
HE W. Study on flexible job shop scheduling with machine breakdown[D]. Chongqing: Chongqing University, 2012.
[15] DAVENPORT A J, BECK J C. Slack-based techniques for robust schedules[C]∥Proc.of the 6th European Conference on Planning, 2001: 43-49.
[16] DEB K, PRATAP A, AGARWAL S, et al. A fast and elitist multi-objective genetic algorithm: NSGA-Ⅱ[J]. IEEE Trans.on Evolutionary Computation, 2002, 6(2): 182-197.
[17] STORN R, PRICE K. Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of Global Optimization, 1997, 1(4): 341-359.
[18] TASGETIREN M F, SUGANTHAN P N, PAN Q K. An ensemble of discrete differential evolution algorithms for solving the generalized traveling salesman problem[J]. Applied Mathematics and Computation, 2010, 215(9): 3356-3368.
[19] PAN Q K, TASGETIREN M F, LIANG Y C. A discrete differential evolution algorithm for the permutation flow shop scheduling problem[J]. Computers & Industrial Engineering, 2008, 55(4): 795-816.
[20] ZITZLER E, DEB K, THIELE L. Comparison of multi-objective evolutionary algorithms: empirical results[J]. Evolutionary Computation, 2000, 8(2): 173-195.
[21] 李志亮, 李小將, 孫偉. 考慮成像質量的敏捷衛星任務調度模型與算法[J]. 宇航學報, 2017, 38(6): 590-597.
LI Z L, LI X J, SUN W. Task scheduling model and algorithm for agile satellite considering imaging quality[J]. Journal of Astronautics, 2017, 38(6): 590-597.