敬玉平，張雨杭，鞏健文，范趙鵬

（海軍航空大學(xué)，山東煙臺(tái)264001）

現(xiàn)代反潛戰(zhàn)中，聲吶浮標(biāo)搜潛具有速度快、探測區(qū)域廣、陣型多樣化等優(yōu)點(diǎn)[1]。浮標(biāo)布陣航路規(guī)劃是指在滿足目標(biāo)約束條件的基礎(chǔ)上為反潛機(jī)選擇一條最優(yōu)的航路完成布陣任務(wù)。多機(jī)布陣航路規(guī)劃問題本質(zhì)上是加入約束條件下的多旅行商問題（mTSP），即給定n個(gè)浮標(biāo)點(diǎn)、m架直升機(jī)，從指定的相同或不同的位置出發(fā)，在布陣時(shí)間最短的原則下，求各架直升機(jī)的路徑。目前，有關(guān)多機(jī)航路規(guī)劃的研究基本還是空白，對于多旅行商問題的研究也沒有特別有效的算法。mTSP的研究開始于20世紀(jì)70年代，剛開始人們致力于采用精確算法進(jìn)行研究，但由于其運(yùn)算量隨問題規(guī)模的擴(kuò)大而呈指數(shù)型增長[2]，研究思路轉(zhuǎn)向了啟發(fā)式算法或元啟發(fā)算法，比如Lin-kernighan算法[3]、自組織神經(jīng)網(wǎng)絡(luò)[4-5]、遺傳算法[6-8]、蟻群算法[9]等。目前，針對mTSP的解決思路大部分是將mTSP轉(zhuǎn)化為單旅行商問題（TSP），再采用TSP的解法來進(jìn)行處理。轉(zhuǎn)化方法主要是在原編碼基礎(chǔ)上加入新的虛擬節(jié)點(diǎn)，來把整個(gè)目標(biāo)點(diǎn)分割成與旅行商數(shù)量相同的幾段。本文將根據(jù)目標(biāo)位置的相對位置和離散程度采用聚類算法來將其分割成塊。接著根據(jù)最長路徑最短的原則，按一定的方法將各塊分配給各個(gè)旅行商，從而將mTSP轉(zhuǎn)化為TSP，再利用模擬退火算法來分別進(jìn)行路徑規(guī)劃。為了進(jìn)一步優(yōu)化規(guī)劃結(jié)果，采用了轉(zhuǎn)移最長路徑節(jié)點(diǎn)的方法，將路徑長度平均化，最后得到了較優(yōu)的規(guī)劃結(jié)果。

1 問題描述與模型建立

為了方便描述，現(xiàn)將問題簡述如下：假設(shè)確定參與布陣的直升機(jī)集合為{F1,F2,…,Fm}，共m架，位置已知，各自攜帶的浮標(biāo)數(shù)為N1,N2,…,Nm，記各架直升機(jī)的初始位置點(diǎn)的集合為，飛行前準(zhǔn)備時(shí)間的集合為{t01,t02,…,t0m}，最大航程的集合為{D1,D2,…,Dm}，浮標(biāo)點(diǎn)共n個(gè)，位置記B1,B2,…,Bn。現(xiàn)在要求直升機(jī)從初始位置出發(fā)，飛過所有浮標(biāo)點(diǎn)，求最短路徑。據(jù)此建立模型如下[5]：

式（1）為目標(biāo)函數(shù)，本文研究的多機(jī)航路規(guī)劃問題以布陣時(shí)間最短為目標(biāo)，假設(shè)每架直升機(jī)的速度相同，則優(yōu)化目標(biāo)等同于使各架直升機(jī)路徑的最大值最小，cij表示的是從點(diǎn)i到點(diǎn)j的路徑成本，通常等于i、j兩點(diǎn)間距離，因而cij=cji。式（2）為執(zhí)行任務(wù)約束，含義為每個(gè)待布設(shè)的浮標(biāo)點(diǎn)必須被訪問且只能被訪問一次。式（3）是任務(wù)可行的必要條件，表示待布設(shè)的浮標(biāo)點(diǎn)數(shù)必須不大于各架直升機(jī)攜帶的浮標(biāo)數(shù)量之和。否則，任務(wù)無法完成。式（4）表示各架直升機(jī)布設(shè)的浮標(biāo)數(shù)量應(yīng)不大于其攜帶量。

為簡化問題，提出假設(shè)：①不考慮浮標(biāo)下落彈道的影響；②不考慮風(fēng)向、風(fēng)速對直升機(jī)的影響；③每架直升機(jī)航速相同且保持勻速；④不考慮飛行前準(zhǔn)備時(shí)間；⑤不考慮洋流、海深對浮標(biāo)投放的影響；⑥不考慮飛行安全問題；⑦每架直升機(jī)的可留空時(shí)間足夠長。

2 求解策略及優(yōu)化方法

一般情況下，可以假設(shè)離得近的浮標(biāo)點(diǎn)更有可能在同一路徑之下。在此基礎(chǔ)之上，采用聚類算法將所有浮標(biāo)點(diǎn)按照相對距離和離散程度劃分為與直升機(jī)數(shù)量相同的幾類。然后，計(jì)算各架直升機(jī)到各類之間的相對距離，類的位置可以根據(jù)類中所有目標(biāo)位置的平均值來表示。再根據(jù)使布陣時(shí)間最短原則，將各類分配給各架直升機(jī)，分配方案還需要滿足式（4）的要求，即各架直升機(jī)負(fù)責(zé)的浮標(biāo)數(shù)量不應(yīng)大于其攜帶量。若不滿足要求，還需要對分配方案進(jìn)行調(diào)整。至此，實(shí)現(xiàn)了由mTSP到TSP之間的轉(zhuǎn)化。接下來，采用模擬退火算法來對各TSP分別進(jìn)行路徑規(guī)劃，從而得到一個(gè)初步的解。即先分割再規(guī)劃。對解的評估方法可以采用均衡率這一指標(biāo)。定義均衡率[10]：

若η接近于1，則任務(wù)分配方案是較合理的。若η遠(yuǎn)小于1，則說明最長路徑遠(yuǎn)大于最短路徑。那么，如果能夠?qū)⒆铋L路徑的部分浮標(biāo)分配給最短路徑，則可減少整體布陣時(shí)間，故還須根據(jù)得到的解的情況在式（4）的約束下進(jìn)行修正。修正方法主要是：將最長路徑的最后一個(gè)浮標(biāo)點(diǎn)取出并放入最短路徑的類中，若最短路徑對應(yīng)的浮標(biāo)任務(wù)量已達(dá)到最大，則交給次短路徑。按這種方法依次嘗試，如優(yōu)化后的分配方案無法實(shí)現(xiàn)，則跳出優(yōu)化的循環(huán)。否則，根據(jù)優(yōu)化后新的分配方案重新轉(zhuǎn)入對各架直升機(jī)采用模擬退火算法進(jìn)行航路規(guī)劃的步驟中。如此循環(huán)，直到無法繼續(xù)優(yōu)化或迭代次數(shù)達(dá)到最大為止。流程圖如圖1所示。

圖1 求解策略及優(yōu)化方法流程圖Fig.1 Solution strategy and optimization method diagram

這樣，在處理單機(jī)航路規(guī)劃中使每條路徑最短，在處理單機(jī)航路規(guī)劃外使各個(gè)路徑長度更平均，兩者結(jié)合，便能達(dá)到使得最長路徑最短的目標(biāo)。

3 算法簡介

3.1 聚類分析

聚類分析又稱群分析，適合將多個(gè)樣本或指標(biāo)進(jìn)行定量分類[10]。在聚類分析中，通常采用距離來對樣本進(jìn)行分類，用到的距離計(jì)算方法主要包括絕對值距離、歐氏距離、切比雪夫距離和馬氏距離等。本文采用歐氏距離。假設(shè)Ω表示樣本點(diǎn)集，則其計(jì)算式為：

針對類與類之間的分類方法，本文采用離差平方和法，若記：

式（8）～（10）中：

定義：

離差平方和法不僅能夠使兩類內(nèi)部距離很小的點(diǎn)聚為一類，還能使2類之間充分分離[10-11]。

聚類分析的主要步驟如下[10，12]。

1）利用歐氏距離計(jì)算n個(gè)樣本點(diǎn)兩兩之間的距離{dij}，建立距離矩陣d={dij}n×n。

2）讓各樣本自成一類，將對應(yīng)聚類圖中的平臺(tái)高度為零。

3）采用離差平方和法，計(jì)算兩類之間的距離，將其中最近的兩類合并為一個(gè)新類，以這2類的距離作為聚類圖中新類的平臺(tái)高度。

4）重新計(jì)算新類與其余各類間的距離。如此時(shí)類的數(shù)量還沒有達(dá)到直升機(jī)的數(shù)量，則返回步驟3）再次執(zhí)行；若數(shù)量已經(jīng)等于直升機(jī)數(shù)量，則進(jìn)入步驟5）。

5）根據(jù)計(jì)算結(jié)果畫出聚類圖并選出對應(yīng)的類。

3.2 布陣任務(wù)分配

計(jì)算每架直升機(jī)所在位置與各類之間的距離，類的位置用類中所有元素坐標(biāo)到飛機(jī)位置的最小值來表示，距離計(jì)算采用歐式距離。由此可以得到一個(gè)m×m的矩陣，m同時(shí)表示直升機(jī)和類的數(shù)量，矩陣的各行表示對于同一架直升機(jī)而言，其相對于各類的距離；各列則相應(yīng)地表示對于同一類而言，其相對于各架直升機(jī)的距離。找到每架直升機(jī)與各類的最小值，也即每一行的最小值，再從這些最小值中取出最大值，找到這個(gè)最大值所對應(yīng)的行和列，即對應(yīng)的直升機(jī)和類，由此完成第一個(gè)任務(wù)分配。然后，刪去最大值對應(yīng)的行和列中的所有元素，對矩陣剩余部分仍采用以上方法進(jìn)行處理，直到全部分配完畢。

3.3 模擬退火算法

模擬退火（Simulated Annealing，SA）算法是由Metropolis等最早提出的，過程分為3個(gè)部分。[13-15]

1）加熱過程。加熱的目的是使粒子的熱運(yùn)動(dòng)增強(qiáng)，偏離平衡位置。此時(shí)粒子的能量較高，可以自由運(yùn)動(dòng)，重新排列。

2）等溫過程。對于與周圍環(huán)境交換熱量而溫度不變的封閉系統(tǒng)，系統(tǒng)狀態(tài)的自發(fā)變化總是朝自由能減少的方向進(jìn)行的，當(dāng)自由能達(dá)到最小時(shí)，系統(tǒng)達(dá)到平衡狀態(tài)。

3）冷卻過程。使粒子熱運(yùn)動(dòng)減弱（該過程也稱退火），系統(tǒng)能量下降，最終得到處于低能狀態(tài)的晶體。

加溫過程對應(yīng)算法的設(shè)定初溫，等溫過程對應(yīng)算法的Metropolis抽樣過程，冷卻過程對應(yīng)控制參數(shù)的下降。這里能量的變化就是目標(biāo)函數(shù)，要得到的最優(yōu)解就是能量最低態(tài)。Metropolis準(zhǔn)則是SA算法收斂于全局最優(yōu)解（指在TSP下的全局最優(yōu)，而非MTSP）的關(guān)鍵所在，Metropolis準(zhǔn)則以一定的概率接受惡化解。這樣，就使算法跳離局部最優(yōu)的陷阱[16-17]。

模擬退火算法的實(shí)現(xiàn)過程描述如下[18]。

1）解空間。解空間S可以表示為的所有固定起點(diǎn)的排列集合，每一個(gè)排列表示一個(gè)可行解，使用蒙特卡洛法得到一個(gè)較好的初始解。

2）目標(biāo)函數(shù)。目標(biāo)函數(shù)為路徑長度，要求路徑長度最短。

3）新解的產(chǎn)生。從上一次迭代得到的解中任意選2個(gè)位置（除了出發(fā)位置），交換這2個(gè)位置之間元素的順序，變成逆序。

4）代價(jià)函數(shù)差。代價(jià)函數(shù)差表示為變換前后2個(gè)解之間目標(biāo)長度的差值。這里用變換后的路徑長度減去變換前的路徑長度。

5）接受長度。如果代價(jià)函數(shù)差f小于0，接受新的路徑；否則，以概率exp(-f/T)接受新路徑。這里由計(jì)算機(jī)產(chǎn)生一個(gè)在[0 ,1]上均勻分布的隨機(jī)數(shù)，若隨機(jī)數(shù)小于或等于exp(-f/T)則接受新解。

6）降溫。利用選定的降溫系數(shù)α進(jìn)行降溫，每次將上一次迭代的溫度T用αT來進(jìn)行更新。

7）結(jié)束條件。設(shè)定最終溫度e，來決定是否結(jié)束退火過程。當(dāng)T＜e時(shí)結(jié)束退火過程并輸出當(dāng)前結(jié)果。否則，返回至第3）步。

其算法流程圖如圖2所示。

圖2 模擬退火算法求解流程框圖Fig.2 Simulated annealing algorithm diagram

3.4 優(yōu)化及分配調(diào)整算法

經(jīng)過前面的步驟可以得到一個(gè)較為理想的mTSP解，但考慮到各條路徑有長有短，如果將最長路徑的一部分分配給較短的路徑，則可以進(jìn)一步優(yōu)化。優(yōu)化過程如下：

1）取點(diǎn)。將得到的最長路徑的最后一個(gè)浮標(biāo)點(diǎn)取出。

2）放入。在仍有布陣能力的直升機(jī)集合中找出路徑最短的一架直升機(jī)，將取出的浮標(biāo)點(diǎn)放入其任務(wù)區(qū)中。

3）重新規(guī)劃。對形成的新任務(wù)分配方案返回至3.3節(jié)重新進(jìn)行航路規(guī)劃。

4）退出循環(huán)。退出循環(huán)的條件有3條，任意滿足一條即退出循環(huán)，得到最終解：①循環(huán)次數(shù)達(dá)到最大；②均衡率達(dá)到預(yù)定值；③當(dāng)除去最長路徑以外的其他路徑對應(yīng)的直升機(jī)的任務(wù)量均已達(dá)到最大；

4 仿真驗(yàn)證

假設(shè)有3架直升機(jī)，經(jīng)緯度坐標(biāo)見表1所示，攜帶的浮標(biāo)數(shù)量均為3枚。假設(shè)指定的布陣任務(wù)為一個(gè)圓陣，共8枚浮標(biāo)，其經(jīng)緯度坐標(biāo)如表2所示。

表1 直升機(jī)經(jīng)緯度Tab.1 Plane latitude and longtitude (°)

表2 浮標(biāo)點(diǎn)經(jīng)緯度Tab.2 Sonobuoy latitude and longtitude (°)

首先，對上述目標(biāo)點(diǎn)進(jìn)行聚類分析，得到的聚類圖如圖3所示。

聚類分析得到的結(jié)果為：第1類的有浮標(biāo)6、7；第2類的有浮標(biāo)4、5；第3類的有浮標(biāo)1、2、3、8。

計(jì)算這3個(gè)類分別與3架直升機(jī)之間的距離，用類中各點(diǎn)到飛機(jī)的最近距離來表示，采用歐式距離計(jì)算。最后，得到的距離矩陣為：

式中，aij表示第i架直升機(jī)與第j類之間的距離，單位為km。

圖3 目標(biāo)點(diǎn)位置聚類圖Fig.3 Target point location cluster diagram

根據(jù)3.2節(jié)中布陣任務(wù)分配方法，可知布陣任務(wù)分配的結(jié)果為：第1架直升機(jī)負(fù)責(zé)第1類浮標(biāo)；第2架直升機(jī)負(fù)責(zé)第2類浮標(biāo)；第3架直升機(jī)負(fù)責(zé)第3類浮標(biāo)。

接下來采用模擬退火算法對3架直升機(jī)分別進(jìn)行路徑規(guī)劃，在同一張圖上顯示。最后得到的結(jié)果如圖4所示。詳細(xì)的數(shù)據(jù)見表3所示。

圖4 路徑規(guī)劃結(jié)果圖Fig.4 Route planning diagram

表3 路徑規(guī)劃數(shù)據(jù)結(jié)果Tab.3 Route planning data result

由表3數(shù)據(jù)，結(jié)合式（6）可以得到均衡率等于0.766 2，最大長度為16.144 4km。最大長度與最小長度相差較大，且第3架直升機(jī)負(fù)責(zé)的浮標(biāo)數(shù)量為4枚，超過了其布放能力。在此基礎(chǔ)之上，采用3.4節(jié)所述的優(yōu)化算法，得到新的路徑規(guī)劃結(jié)果如圖5所示。詳細(xì)的數(shù)據(jù)見表4所示。

圖5 優(yōu)化后路徑規(guī)劃圖Fig.5 Route planning result diagram after optimization

表4 優(yōu)化后路徑規(guī)劃數(shù)據(jù)結(jié)果Tab.4 Route planning data result after optimization

由表4數(shù)據(jù)結(jié)合式（6）可以得到均衡率約等于0.890 5，最大長度為14.445 9km。由圖5和表4的結(jié)果可知：經(jīng)過優(yōu)化處理后，均衡率有所提高，而且最大長度變大。同時(shí)，浮標(biāo)數(shù)量的約束條件得以滿足。因此，可以認(rèn)為優(yōu)化結(jié)果是可接受的。

5 結(jié)束語

針對多機(jī)布陣航路規(guī)劃問題，本文給出了一種切實(shí)可行的算法，通過對任務(wù)的分解，能夠解決多機(jī)航路規(guī)劃問題中的協(xié)同與路徑規(guī)劃問題，在此基礎(chǔ)上，又給出了優(yōu)化的方法。仿真結(jié)果顯示，經(jīng)過優(yōu)化后的結(jié)果是可信的，且有較高的效率。但是在航路規(guī)劃中沒有考慮飛行的安全性，布設(shè)的方便性，以及氣象條件等因素對航路規(guī)劃的影響。接下來，還需要針對以上幾個(gè)問題加以改進(jìn)。