曾 斌, 姚 路, 李厚樸
(1. 海軍工程大學(xué)管理工程與裝備經(jīng)濟(jì)系, 湖北 武漢 430033;2. 海軍工程大學(xué)導(dǎo)航工程系, 湖北 武漢 430033)
反潛戰(zhàn)是海軍作戰(zhàn)的主要樣式之一,是同對方潛艇兵力作斗爭的海上作戰(zhàn)。由于海洋環(huán)境復(fù)雜多變,潛艇具有非常強(qiáng)的隱蔽性,給反潛作戰(zhàn)帶來很大難度。隨著無人水下航行器的發(fā)展,反潛的地位和作用更為重要,其中艦載反潛直升機(jī)由于其具有飛行速度快、搜潛范圍廣、隱蔽性好以及攻擊效率高等優(yōu)點,是海上艦艇編隊對潛作戰(zhàn)的主要兵力[1]。
從作戰(zhàn)目標(biāo)上看反潛分為護(hù)航反潛和區(qū)域反潛,本文主要針對區(qū)域反潛,區(qū)域反潛又細(xì)分為:① 反潛封鎖區(qū),使用陣地防潛器材和機(jī)動反潛兵力以阻止對方潛艇通過的海區(qū);② 海洋反潛區(qū),將海洋劃分成若干個反潛區(qū),各區(qū)編組一定的機(jī)動反潛兵力,分區(qū)負(fù)責(zé)反潛搜索和跟蹤;③ 反潛警戒線,在沿海和島嶼附近的海岸聲納組網(wǎng),以便提供反潛預(yù)警,并輔助反潛兵力進(jìn)行反潛搜索。
近年來,反潛模式也由平臺中心戰(zhàn)向網(wǎng)絡(luò)中心戰(zhàn)轉(zhuǎn)型[2],以此縮短反潛觀察-調(diào)整-決策-行動閉環(huán)的時間。通過各類傳感器等預(yù)警偵察手段,獲得準(zhǔn)確的戰(zhàn)場態(tài)勢;利用反潛機(jī)的快速機(jī)動能力,實施召喚式快速搜索[3]。
反潛警戒海域柵格化為若干監(jiān)測區(qū),在水下聲納網(wǎng)絡(luò)[4-5]、巡邏反潛、偵察衛(wèi)星、無線電偵聽配合下,為重點海域(海上武器試驗區(qū)、重要航道、海上保障基地、艦艇編隊等)提供反潛預(yù)警;當(dāng)發(fā)現(xiàn)水下可疑目標(biāo)后向反潛指揮部門發(fā)出警報,警報中包含該目標(biāo)隨處區(qū)域以及反潛需求;指揮部門接到警報后,再調(diào)度我方反潛直升機(jī)到達(dá)該監(jiān)測區(qū)實施搜潛、跟蹤、攻擊或驅(qū)逐,反潛機(jī)完成任務(wù)后返回待召。
隨著無人水下航行器(unmanned underwater vehicle,UUV)的發(fā)展,反潛作戰(zhàn)面臨的挑戰(zhàn)更加艱巨。藍(lán)方往往頻繁地從多處對紅方警戒海域進(jìn)行水下侵入,再加上水下環(huán)境復(fù)雜[6],水下預(yù)警系統(tǒng)也存在誤報的情況,而反潛直升機(jī)資源數(shù)量有限,傳統(tǒng)的就近派遣策略無法滿足這種多次數(shù)多分散地點的反潛需求。因此,優(yōu)化的調(diào)度策略更顯重要。
直升機(jī)反潛是當(dāng)前軍事領(lǐng)域研究熱點問題之一。例如,文獻(xiàn)[7]利用貝葉斯網(wǎng)絡(luò)推理潛艇意圖。文獻(xiàn)[8]提出了適合于雙機(jī)搜潛的“前堵后追”法。文獻(xiàn)[9]提出了一種直升機(jī)協(xié)同搜潛的陣型最優(yōu)化方法。文獻(xiàn)[10]研究了聯(lián)合反潛作戰(zhàn)中反潛機(jī)的搜索模式問題。文獻(xiàn)[11]將搜索海域柵格化,給出了艦載直升機(jī)吊放聲納區(qū)域反潛搜索算法。文獻(xiàn)[12]研究了反潛機(jī)應(yīng)召搜索時聲納浮標(biāo)陣的評估模型。文獻(xiàn)[13]從協(xié)同反潛角度研究了對潛艇的探測模型。以上文獻(xiàn)主要研究反潛直升機(jī)起飛后的搜潛策略,但反映出反潛作戰(zhàn)的背景特征,如潛艇位置的隨機(jī)性和搜索海域柵格化等思路,對本文具有較大啟發(fā)。
文獻(xiàn)[14]以提高覆蓋率為目標(biāo),建立了多艦聯(lián)合搜潛的概率方程,并提出基于概率距離函數(shù)的優(yōu)化算法來求解。文獻(xiàn)[15]研究了瀕海戰(zhàn)斗艦的反潛規(guī)劃問題。文獻(xiàn)[16]研究了應(yīng)召反潛中,水面艦艇和反潛機(jī)協(xié)同反潛中拖曳聲納和反潛直升機(jī)的部署問題。文獻(xiàn)[17]利用無人艇執(zhí)行搜潛任務(wù),提出了改進(jìn)的遺傳算法來提高累積探測概率。文獻(xiàn)[18]研究了反潛直升機(jī)的優(yōu)化攻擊地點。文獻(xiàn)[19]設(shè)計了反潛機(jī)的探測距離估算模型。以上文獻(xiàn)在反潛監(jiān)測概率和優(yōu)化部署方面對本文具有一定啟發(fā)。
另外,國外軍事研究人員,特別是美國海軍也針對反潛直升機(jī)的作戰(zhàn)使用展開了相關(guān)研究。文獻(xiàn)[20]設(shè)計了一個上下文驅(qū)動的決策支持工具來解決反潛調(diào)度和規(guī)劃問題,提出了一個隱馬爾可夫模型來對反潛資源和搜索規(guī)劃建模,并提出了二階段算法(進(jìn)化算法和拍賣算法)來求解。文獻(xiàn)[21-22]利用零和博弈對反潛作戰(zhàn)建模,利用線性規(guī)劃算法計算反潛資源(反潛機(jī)或艦艇)的優(yōu)化部署問題。文獻(xiàn)[23]利用非零和網(wǎng)絡(luò)阻斷博弈對區(qū)域反潛建模,并采用強(qiáng)Stackelberg均衡策略解決了反潛資源的分配問題。文獻(xiàn)[24-25]以提高覆蓋率為目標(biāo),研究了反潛規(guī)劃問題,并運用解析法建立了反潛責(zé)任區(qū)防潛兵力的配置模型。文獻(xiàn)[26]研究了帶吊放聲納的反潛直升機(jī)搜索模式和吊放頻率。文獻(xiàn)[27]進(jìn)一步研究了反潛直升機(jī)搜索策略以及油料和載荷的優(yōu)化配置問題。但現(xiàn)有文獻(xiàn)缺少從水下網(wǎng)絡(luò)中心戰(zhàn)角度研究多監(jiān)測區(qū)-多反潛機(jī)的調(diào)度問題。
本文的研究目標(biāo)是:把警戒海域劃分為多個監(jiān)測區(qū),通過歷史數(shù)據(jù)或經(jīng)驗預(yù)估監(jiān)測區(qū)的報警率和期望搜潛時間,在此基礎(chǔ)上,計算得到優(yōu)化的多反潛直升機(jī)部署地點以及分派策略(優(yōu)先順序),從而達(dá)到最大覆蓋警戒海域的目的。
多監(jiān)測區(qū)-多反潛機(jī)調(diào)度中存在的兩個問題增加了優(yōu)化過程的復(fù)雜性:① 需要綜合考慮不同類型反潛直升機(jī)的調(diào)度,不同類型反潛直升機(jī)各有自己的優(yōu)缺點,如何在反潛作戰(zhàn)中相互配合,揚其所長,避其所短,是調(diào)度模型的一個重點研究目標(biāo);② 各反潛直升機(jī)的負(fù)載平衡問題,如果調(diào)度(部署和分派策略)不當(dāng),將導(dǎo)致某一架或幾架直升機(jī)工作負(fù)載過大,這會加大飛行員的疲勞度,同時也會增加飛機(jī)的故障率,從而降低整個反潛系統(tǒng)的可用率。
因此,本文的主要創(chuàng)新點和貢獻(xiàn)是從水下網(wǎng)絡(luò)中心戰(zhàn)的角度出發(fā),為了達(dá)到最大覆蓋警戒海域的研究目標(biāo),特別針對以上兩個問題,提出反潛直升機(jī)部署和調(diào)度模型及求解算法。在技術(shù)上借鑒了超立方體排隊模型來研究反潛調(diào)度問題,該模型主要用于分析應(yīng)急服務(wù)系統(tǒng)的最大覆蓋問題,在緊急醫(yī)療設(shè)施[28]和急救車部署[29-30]方面取得了較好的性能。但相關(guān)研究沒有考慮不確定情況下,多優(yōu)先級任務(wù)同時出現(xiàn)時不同服務(wù)單元的配置問題,這也是本文的研究難點。
本文分兩個階段進(jìn)行優(yōu)化:第1階段建立超立方體反潛調(diào)度排隊模型和求解算法,以輸入的報警率作為事件到達(dá)率,期望響應(yīng)時間為平均服務(wù)時間,輸出校正系數(shù)和直升機(jī)繁忙概率等統(tǒng)計結(jié)果;第2階段建立反潛調(diào)度規(guī)劃模型,以第1階段輸出作為規(guī)劃模型的輸入?yún)?shù),從而得到反潛直升機(jī)的部署及分派策略。
超立方體排隊模型是對常規(guī)排隊模型的擴(kuò)展,適用于應(yīng)急反應(yīng)系統(tǒng)的概率選址,反映了服務(wù)器與客戶需求之間更為復(fù)雜的分配關(guān)系,對多服務(wù)器排隊系統(tǒng)進(jìn)行了擴(kuò)展,增加了服務(wù)器狀態(tài)空間的描述,能夠更為精確地刻畫排隊系統(tǒng)的時空復(fù)雜性[31]。某一時間段內(nèi)服務(wù)器可能處于忙(1)或閑(0)兩種狀態(tài),系統(tǒng)中多個服務(wù)器的狀態(tài)組合構(gòu)成了系統(tǒng)狀態(tài)空間。例如狀態(tài){1,1,0}表達(dá)了一個三服務(wù)器系統(tǒng),其中1號服務(wù)器空閑,2號和3號服務(wù)器繁忙,從中可以看出一個三服務(wù)器系統(tǒng)的狀態(tài)空間可以由立方體的8個頂點表示,如果系統(tǒng)中服務(wù)器數(shù)量多于3個,則稱之為超立方體。如果給定排隊系統(tǒng)參數(shù),例如服務(wù)器數(shù)量、分管區(qū)域數(shù)量(客戶數(shù)量)、客戶請求到達(dá)率以及服務(wù)器的平均服務(wù)時間等,利用超立方體排隊模型可以計算穩(wěn)態(tài)狀態(tài)下系統(tǒng)的一些關(guān)鍵性能指標(biāo),包括服務(wù)器的繁忙概率、工作強(qiáng)度以及服務(wù)器與客戶之間的分配概率等。
反潛直升機(jī)數(shù)量有限,在警戒海域執(zhí)行大規(guī)模搜潛任務(wù)時不能處于隨時可用狀態(tài),所以計算其繁忙概率對于應(yīng)召搜潛調(diào)度而言非常重要。另外直升機(jī)執(zhí)行反潛任務(wù)時存在協(xié)作關(guān)系,因此本文還需要計算校正系數(shù)來補(bǔ)償排隊系統(tǒng)中服務(wù)器無關(guān)性的假設(shè)。
由于排隊模型中引用參數(shù)較多,在這里首先統(tǒng)一進(jìn)行描述。
nA(nB)為系統(tǒng)中A類型(B類型)反潛機(jī)數(shù)量,nA+nB=n;
LA(LB)為A類型(B類型)反潛機(jī)候選部署地點集合,LA+LB=L,包括多個地點l∈L;
Z為警戒海域內(nèi)監(jiān)測網(wǎng)格區(qū)集合,當(dāng)預(yù)警系統(tǒng)發(fā)現(xiàn)可疑目標(biāo)后,向反潛調(diào)度系統(tǒng)發(fā)出包括目標(biāo)大致所在區(qū)域(監(jiān)測區(qū)編號)的警報,由調(diào)度系統(tǒng)通知反潛機(jī)前往該區(qū)域?qū)嵤┧褲?包括多個監(jiān)測區(qū)z∈Z;
T為報警類別集合,本文把警報類別分為4類:① 只能分配給A型機(jī)處理的警報TA;② 如果A型機(jī)都忙時,可以分配給B型機(jī)處理的警報TAB,即優(yōu)先要求A類型機(jī);③ 如果B型機(jī)都忙時,可以分配給A型機(jī)處理的警報TBA;④ 只能分配給B型機(jī)處理的警報TB,TA、TB、TAB、TBA?T;
λt為類型為t的警報發(fā)生率,t∈T;
λtz為在第z監(jiān)測區(qū)發(fā)出的類型為t的警報發(fā)生率,λt=∑z∈Zλtz,z∈Z,t∈T;
λ為系統(tǒng)級所有警報發(fā)生率,λ=∑z∈Z,t∈Tλtz;
τAlzt(τBlzt)為l地點的A型機(jī)(B型機(jī))執(zhí)行z區(qū)t類別警報任務(wù)的期望響應(yīng)時間(從收到警報到抵達(dá)目標(biāo)區(qū)域z所花費的時間);
τA(τB)為系統(tǒng)級A型機(jī)(B型機(jī))的平均響應(yīng)時間;
θt為t類別警報對響應(yīng)時間指標(biāo)的要求,本文以規(guī)定響應(yīng)時間內(nèi)得到處理的t類別警報數(shù)量比例來作為衡量指標(biāo);
ρA(ρB)為A型機(jī)(B型機(jī))的工作強(qiáng)度;
pAl(pBl)為部署在l地點的A型機(jī)(B型機(jī))的繁忙概率;
pA(pB)為系統(tǒng)級A型機(jī)(B型機(jī))的繁忙概率,pA=(∑l∈LApAl)/nA;
aAztk(aBztk)為A型機(jī)(B型機(jī))的排序調(diào)度表,表示當(dāng)z區(qū)域發(fā)出t類別警報時,作為第k順位分配的反潛機(jī)地點,如果t∈TA∪TAB,k=1,2, …,nA,如果t∈TBA,則k=nB+1,nB+2,…,nB+nA;
dAlztk(dBlztk)為A型機(jī)(B型機(jī))的分配概率,表示當(dāng)z區(qū)域發(fā)出t類別警報時,作為第k順位的l地點A型機(jī)(B型機(jī))被分派到處理該警報的概率;
PAk(PBk)為由于k個A型機(jī)(B型機(jī))的繁忙,造成警報無法處理的概率,即損失率;
RAlzt(RBlzt)為表示當(dāng)z區(qū)域發(fā)出t類別警報時,l地點的A型機(jī)(B型機(jī))能夠在規(guī)定響應(yīng)時間(與警報類別有關(guān))內(nèi)到達(dá)z區(qū)域執(zhí)行搜潛任務(wù)的概率;
QA(nA,ρA,k)或QB(nB,ρB,k)為超立方體模型的校正系數(shù),可看作輸入?yún)?shù)為反潛機(jī)數(shù)量nA(nB),反潛機(jī)工作強(qiáng)度ρA(ρB)以及派遣順位k的函數(shù)。
本文假設(shè)z區(qū)域t類型警報的到達(dá)率服從參數(shù)為λtz的泊松分布,如果某警報需要多架反潛機(jī)協(xié)同反潛,可以發(fā)出多個虛擬警報來處理。
為了避免重復(fù),本文中主要描述A型機(jī)的數(shù)學(xué)模型,B型機(jī)的模型與A型機(jī)相同。
1.2.1 校正系數(shù)的計算模型
當(dāng)出現(xiàn)警報時,在前面連續(xù)j架反潛機(jī)都繁忙,第k+1順位的反潛機(jī)空閑的概率是在排隊論服務(wù)器無關(guān)性假設(shè)下計算得到的,在反潛系統(tǒng)中,反潛機(jī)的調(diào)度往往存在協(xié)作關(guān)系,因此需要在模型中加入無關(guān)性校正系數(shù)QA(nA,ρA,k)來修正[32]。設(shè)PAnA表示所有的nA架A型反潛機(jī)繁忙的概率,PA0表示沒有A型機(jī)繁忙,即所有A型機(jī)空閑的概率,設(shè)k=0,1,…,nA-1,則有
QA(nA,ρA,k)=
(1)
1.2.2 分配概率的計算模型
在排隊模型中,如果要計算服務(wù)器的繁忙概率,必須先計算服務(wù)器-客戶對的分配概率。在本文中,服務(wù)器表示為反潛機(jī),客戶可表示為監(jiān)測區(qū),所以分配概率既是dAlztk或dBlztk的值。分配概率與排序調(diào)度表aztk有關(guān),如果a1A2=3,表示1號區(qū)域發(fā)出只能由A型機(jī)反潛的警報時,位于3號地點的反潛機(jī)將作為第2順位被分配去處理。另外分配概率與警報類別有關(guān),因為對于TBA類別的警報,必須要B型機(jī)都繁忙時才輪到分配A型機(jī)去處理,所以對于處理TA∪TAB類別的警報和處理TBA類別的警報,其分配概率的計算是不同的。特別是當(dāng)t∈TB時,因為該類警報只能由B型機(jī)處理,所以這時dAlztk=0。綜上所述,dAlztk的計算公式分以下兩種情況:
(1)如果t∈TA∪TAB,z∈Z,k=1,2,…,nA,l∈LA,則

(2)
(2)如果t∈TBA,z∈Z,k=nB+1,nB+2,…,nB+nA,l∈LA,則
(3)

1.2.3 繁忙概率的計算模型
地點l的A型機(jī)繁忙概率pAl計算公式為
(4)
式中,rAl表示地點l的A型機(jī)分配執(zhí)行搜索任務(wù)的頻率,計算公式為
(5)

式(5)中第二項表示地點l的A型機(jī)分派給TBA類別警報任務(wù)的頻率,同第1.2.2節(jié)分派概率的計算,加入PBnB項的原因是,在t∈TBA時,所有的B型機(jī)都繁忙時,才只好找A型機(jī)調(diào)度。而且這時A型機(jī)的分派優(yōu)先級順序都排在B型機(jī)后面,所以優(yōu)先序號j從nB+1開始。
1.2.4 平均響應(yīng)時間的計算模型
平均響應(yīng)時間的計算公式如下:
(6)
式中,平均響應(yīng)時間是單架機(jī)響應(yīng)時間τAlzt與歸一化的分派概率dAlztk/(1-PAnA)的乘積,乘積結(jié)果再用警報到達(dá)頻率(λtz/λ)加權(quán)。
由于穩(wěn)態(tài)時第1.2節(jié)的排隊模型無法求得精確解,只能采用近似算法迭代計算。近似算法步驟如下。
算法輸入?yún)?shù):①λtz,警報發(fā)生率;②τAlzt,l地點的A型機(jī)響應(yīng)時間;③aAztk,可由第2節(jié)規(guī)劃模型求解算法輸入。
初始化:初始化平均響應(yīng)時間為
(7)
按埃爾朗損失率公式初始化所有A型機(jī)繁忙的概率,即
(8)
式中,A型機(jī)初始空閑概率PA0,0為
(9)
設(shè)Lz為距離監(jiān)測區(qū)z最近的反潛機(jī)地點集合,則反潛機(jī)的初始繁忙概率為
(10)
迭代計數(shù)器h設(shè)為0。
步驟 1利用上一次迭代計算的平均服務(wù)時間τA,h-1(如果為第1次迭代,則利用初始化的τA,0)更新?lián)p失率,計算公式為
(11)
式中,A型機(jī)空閑概率為
(12)
步驟 2設(shè)k=0,1,…,nA-1,按式(1)計算本次迭代的校正系數(shù)QA,h(nA,ρA,h,k),其中工作強(qiáng)度ρA,h=λτA,h-1/nA。
步驟 3按式(2)和式(3)計算本次迭代h的分配概率dAlztk,h,需要的參數(shù)中繁忙概率使用上一迭代的計算值,校正系數(shù)取步驟2的計算結(jié)果。
步驟 4按式(6)計算本次迭代的平均響應(yīng)時間τA,h,需要的參數(shù)中分派概率使用步驟4的計算結(jié)果,損失率使用步驟1的計算結(jié)果。
步驟 5按式(4)和式(5)計算當(dāng)前h迭代時的反潛機(jī)繁忙概率pAl,h,需要的參數(shù)中優(yōu)先級排序的繁忙概率使用上一迭代的計算值,校正系數(shù)取步驟2的計算結(jié)果。
步驟 6若對于所有地點的A類反潛機(jī),其前后兩次迭代的計算結(jié)果大于某一閾值ε,即|pAl,h-pAl,h-1|>ε,其中l(wèi)∈LA,則跳轉(zhuǎn)步驟1繼續(xù)迭代,否則返回繁忙概率和校正系數(shù)。
反潛機(jī)調(diào)度規(guī)劃模型的目標(biāo)函數(shù)是最大化規(guī)定時間內(nèi)完成高優(yōu)先級反潛任務(wù)的數(shù)量比例,并且在此同時實現(xiàn)對低優(yōu)先級警報的一定(給定閾值內(nèi))覆蓋率,這樣折中實現(xiàn)對所有警報的及時處理。
排隊模型近似求解算法的輸出用作規(guī)劃模型的輸入?yún)?shù)。規(guī)劃模型目標(biāo)函數(shù)和約束條件中用到的兩個系數(shù)與超立方體排隊模型輸出相關(guān)。一個系數(shù)捕捉系統(tǒng)分派策略滿足時限要求的性能,本文稱之為及時率;另一個系數(shù)捕捉系統(tǒng)的工作負(fù)載。
及時率系數(shù)cAlztk或cBlztk表示當(dāng)z區(qū)域發(fā)出t類別警報時,作為第k順位的l地點的A型機(jī)(B型機(jī))能夠在規(guī)定響應(yīng)時間(與警報類別有關(guān))內(nèi)到達(dá)z區(qū)域的數(shù)量比例,為反映監(jiān)測區(qū)的優(yōu)先級,本文使用監(jiān)測區(qū)z的警報率(λtz/λt)對及時率系數(shù)加權(quán)。假設(shè)同一類型反潛機(jī)的繁忙概率相同時,及時率的計算類似于分派概率,其計算公式為
(1)如果t∈TA∪TAB,z∈Z,k=1,2,…,nA,l∈LA,則
cAlztk=(λtz/λt)RAlzt[QA(nA,ρA,k)(1-pA)(pA)k-1]
(13)
(2)如果t∈TBA,z∈Z,k=nB+1,nB+2,…,nB+nA,l∈LA,則
cAlztk=(λtz/λt)RAlztPBnB[QA(nA,ρA,k)(1-pA)(pA)k-1]
(14)
式中,校正系數(shù)與A型機(jī)繁忙概率pA都為超立方體排隊模型的計算結(jié)果。及時率系數(shù)計算過程中考慮了兩個不確定因素的影響:① 在校正系數(shù)補(bǔ)償下反潛機(jī)可用性的不確定性,比地點l更優(yōu)先的k個反潛機(jī)繁忙時至少有一個反潛機(jī)空閑的概率為(1-pA)(pA)k-1,再乘以校正系數(shù)QA,表示有k-1架反潛機(jī)優(yōu)先于l調(diào)度,但l仍然被調(diào)度的概率;②RAlzt表示飛行和搜潛時間的不確定性。
同樣地,負(fù)載系數(shù)gAlztk或gBlztk,分別表示A型機(jī)或B型機(jī)的工作負(fù)載,A型機(jī)負(fù)載系數(shù)計算公式為
(1)如果t∈TA∪TAB,z∈Z,k=1,2,…,nA,l∈LA,則
gAlztk=(λtz/λt)τAlzt[QA(nA,ρA,k)(1-pA)(pA)k-1]
(15)
(2)如果t∈TBA,z∈Z,k=nB+1,nB+2,…,nB+nA,l∈LA,則
gAlztk=(λtz/λt)τAlztPBnB[QA(nA,ρA,k)(1-pA)(pA)k-1]
(16)
式中,τAlzt表示l地點的A型機(jī)執(zhí)行z區(qū)t類別警報任務(wù)的期望響應(yīng)(服務(wù))時間。B型機(jī)計算公式類似,負(fù)載系數(shù)用于平衡同一類型反潛機(jī)的工作負(fù)載。
首先定義反潛規(guī)劃模型的決策變量,為了完成反潛調(diào)度的兩個主要目標(biāo),反潛規(guī)劃模型相應(yīng)也包含兩套決策變量。布爾變量hAl和hBl決定A型機(jī)和B型機(jī)的部署地點,如果hAl=1,表示有A型機(jī)部署在l地點(l∈LA),否則hAl為0。hBl與hAl相似,如果hBl=1,表示有B型機(jī)部署在l地點(l∈LB),否則hBl為0。
第二套決策變量為mAlztk和mBlztk也為布爾變量,表示A型機(jī)和B型機(jī)分派至監(jiān)測區(qū)的排序調(diào)度表,與aAztk或aBztk可以相互轉(zhuǎn)化,一方面用于表示對于不同類別警報的分派順序,例如TA類別警報任務(wù)只能由A型機(jī)執(zhí)行,TAB類別警報任務(wù)先由A型機(jī)執(zhí)行,如果沒有A型機(jī)空閑時才交由B型機(jī)執(zhí)行,TB和TBA也是類似。另外,也可以表示當(dāng)發(fā)生某警報時,同一類型反潛機(jī)內(nèi)部的分派順序,如果mAlztk=1,表示當(dāng)z區(qū)域發(fā)出t類別警報時,地點l的A型機(jī)作為第k順位分配,如果t∈TA∪TAB,k=1,2,…,nA,如果t∈TBA,則k=nB+1,nB+2, …,nB+nA;其他情況時mAlztk=0。mBlztk的設(shè)計與mAlztk相似。
以上決策變量的設(shè)計把反潛調(diào)度規(guī)劃轉(zhuǎn)化為二元線性規(guī)劃模型,目標(biāo)函數(shù)為
(17)
該目標(biāo)函數(shù)的目的是最大化高優(yōu)先級警報TA的期望覆蓋率,在反潛調(diào)度規(guī)劃模型中,只有A型機(jī)可以執(zhí)行TA警報任務(wù)。
約束關(guān)系如下:
(18)
nB+1,nB+2,…,nB+nA
(19)
(20)
k=nA+1,k=nA+2,…,nA+nB
(21)
(22)
(23)

(24)
(25)
(26)
(27)
(28)
(29)
(30)
mAlztk≤hAl,?l∈LA,z∈Z,t∈TA∪TAB,
k=1,2,…,nA
(31)
mAlztk≤hAl,?l∈LA,z∈Z,t∈TBA,k=nB+1,
nB+2,…,nB+nA
(32)
mBlztk≤hBl,?l∈LB,z∈Z,t∈TB∪TBA,
k=1,2,…,nB
(33)
mmlztk≤hBl,?l∈LB,z∈Z,t∈TAB,k=nA+1,
nA+2,…,nA+nB
(34)
(35)
(36)
(37)
(38)
式(18)表示對于警報t∈TA∪TAB,在排序調(diào)度表中同一順位只能選擇一架A型機(jī)執(zhí)行。式(19)表示對于警報t∈TBA,首先分派B型機(jī)執(zhí)行,如果B型機(jī)都繁忙才分派給A型機(jī),而且在排序調(diào)度表中同一優(yōu)先級順位只能選擇一架A型機(jī)執(zhí)行。式(20)和式(21)與式(18)和式(19)表示的意義相同,但針對的是B型機(jī)。式(22)~式(25)表示某反潛機(jī)l不能同時是某監(jiān)測區(qū)z的多個順位分派機(jī),比如不能既是z的第1順位分派機(jī)又是z的第2順位分派機(jī)。式(26)~式(28)保證低優(yōu)先級的警報(TAB,TBA,TB)的期望覆蓋率大于某個閾值θt。式(29)和式(30)為反潛機(jī)的數(shù)量約束。式(31)~式(34)表示只有已部署到指定地點的反潛機(jī)才能被調(diào)度。式(35)和式(36)約束同一型號反潛機(jī)之間負(fù)載平衡,即負(fù)載差不高于某一閾值δA(δB)。式(37)和式(38)約束反潛機(jī)的分派監(jiān)測區(qū)相鄰,即如果z的鄰居z′分派給某反潛機(jī)l分管,則z才能被該反潛機(jī)分管。
本文中,由式(18)~式(34)構(gòu)成的模型為基本調(diào)度模型,基本模型增加了式(35)和式(36)后為負(fù)載平衡模型,基本模型增加了式(37)和式(38)后增加了分管區(qū)域相鄰的約束。所有約束構(gòu)成了加強(qiáng)模型。另外假設(shè)反潛機(jī)在我方警戒海域可以就近補(bǔ)充油料,所以約束關(guān)系中沒有包括油料限制。
由于排隊模型求解算法中反潛機(jī)的繁忙概率與分派策略aAztk(aBztk)有關(guān),而aAztk(aBztk)可以通過mAlztk(mBlztk)轉(zhuǎn)換得到,規(guī)劃模型又需要排隊模型中的校正系數(shù)和繁忙概率,所以本文設(shè)計了迭代算法求解規(guī)劃模型,其計算步驟如下。
步驟 1初始化;
步驟 1.1利用“分配最近空閑反潛機(jī)”策略初始化初始排序調(diào)度表aAztk(aBztk);
步驟 1.2調(diào)用第1.3節(jié)排隊模型求解算法得到繁忙概率和校正系數(shù);
步驟 1.3設(shè)置初始不平衡容忍度δ0=1;
步驟 1.4迭代次數(shù)i設(shè)為1;
步驟 2利用上次迭代(或初始階段)生成的繁忙概率和校正系數(shù)求解規(guī)劃模型,得到新的排序調(diào)度表aAztk,i;
步驟 3利用新的排序調(diào)度表aAztk,i,調(diào)用排隊模型求解算法,得到新的繁忙概率pA,i和校正系數(shù)Qi,設(shè)置δA,i=1/2(maxpA,i-minpA,i);
步驟 4算法包括2個迭代停止條件。第1個停止條件為式(33)和式(34)無法滿足,這時再減少不平衡容忍度δ會導(dǎo)致無法得到可行解;第2個停止條件為繁忙概率小于某個閾值。如果算法停止條件都沒有達(dá)到,增加迭代次數(shù)i并跳轉(zhuǎn)至步驟2;
步驟 5返回最新迭代計算的排序調(diào)度表和反潛機(jī)部署地點。
具體過程如圖1所示。

圖1 柵格化反潛警戒區(qū)示意圖Fig.1 Schematic diagram of grid antisubmarine warning area
本算法的應(yīng)用時機(jī)為需要對某海域布防警戒前,具體應(yīng)用步驟描述如下。
步驟 1類似圖1,把需要布防海域柵格化為若干監(jiān)測區(qū)z,同時篩選出適合部署反潛直升機(jī)的候選地點l;
步驟 2通過歷史數(shù)據(jù)、地點重要程度或敵方威脅程度預(yù)估各個監(jiān)測區(qū)的報警率λtz,并根據(jù)l和z之間距離和反潛機(jī)速度估算出期望響應(yīng)時間τAlzt(或τBlzt),如果可能也可以估計期望響應(yīng)時間的不確定性參數(shù)RAlzt或RBlzt(缺省為1),θt根據(jù)相關(guān)條例推導(dǎo)指定;
步驟 3步驟1和步驟2的結(jié)果作為調(diào)度算法輸入?yún)?shù)進(jìn)行迭代計算。由于調(diào)度算法的步驟1需要計算第2.2節(jié)的規(guī)劃模型(二元線性規(guī)劃),而該規(guī)劃模型中用到的及時率和工作負(fù)載2個系數(shù)需要用到排隊模型輸出的繁忙概率和校正系數(shù),所以需要嵌入調(diào)用第1.3節(jié)的排隊模型求解算法;
步驟 4當(dāng)調(diào)度算法退出時,生成兩個決策變量的計算結(jié)果,反潛機(jī)的部署地點(hAl和hBl)以及排序調(diào)度表(mAlztk和mBlztk)。
指揮人員可以以本算法結(jié)果作為決策支持建議,在推薦地點部署反潛機(jī)搭載平臺(驅(qū)逐艦或護(hù)衛(wèi)艦)。如果某監(jiān)測區(qū)出現(xiàn)威脅目標(biāo)時,能夠以排序調(diào)度表作為調(diào)度策略基礎(chǔ)方案,按優(yōu)先級次序(如果高優(yōu)先級直升機(jī)已執(zhí)行任務(wù),則依次派遣低優(yōu)先級直升機(jī))派遣對應(yīng)反潛機(jī)執(zhí)行對應(yīng)任務(wù)。
本文仿真背景案例采用了類似圖1的柵格化警戒海域背景,整個海域分為85個監(jiān)測區(qū),其中32個區(qū)域可以作為反潛直升機(jī)搭載平臺的候選部署地點,假設(shè)有10架A類機(jī)和4架B類機(jī)可供調(diào)度。從安全角度考慮,本文不討論警戒海域和反潛機(jī)的具體參數(shù)。
在參數(shù)設(shè)置方面,期望響應(yīng)時間τAlzt以小時為單位,按反潛機(jī)的部署地點l和報警區(qū)域z之間的距離除以反潛機(jī)的速度計算。τBlzt的計算與A類機(jī)類似,但考慮B類機(jī)可能采取雙機(jī)協(xié)同,且綜合性能較弱,所以其計算得到的響應(yīng)時間乘以2作為懲罰因子。警報類別的概率分布為:TA為75%,TAB為20%,TBA為5%,仿真中不生成TB警報。系統(tǒng)級所有警報發(fā)生率λ=3 h,每一個監(jiān)測區(qū)域的警報發(fā)生率相等。
當(dāng)前算法在聯(lián)想筆記本上即可運行,CPU為四核Intel-i5,基礎(chǔ)頻率1.6 GHZ,內(nèi)存8 G,算法在VS2012下用C++實現(xiàn),規(guī)劃求解部分調(diào)用IBM ILOG CPLEX 12.9。容忍度設(shè)置為0.01時,平均運行時間統(tǒng)計如下:基本模型耗時約0.6 min,負(fù)載平衡模型耗時約6.2 min,加強(qiáng)模型耗時約11.3 min,3個規(guī)劃模型中調(diào)用排隊算法更新繁忙概率和校正系數(shù)的時間分別為0.1 min、0.65 min和0.6 min。
圖2和圖3為部署了全部10架A型機(jī)和4架B型機(jī)(B型機(jī)可以為雙機(jī)協(xié)同)時,在迭代調(diào)度算法中調(diào)用排隊模型求解算法得到的繁忙概率pAl(pBl),其中,l表示10架A型機(jī)和4架B型機(jī)的部署地點。圖2中A型機(jī)的繁忙概率最低為0.052,最高為0.105,平均為0.075;圖3中B型機(jī)的繁忙概率最低為0.225,最高為0.357,平均為0.310,所有機(jī)型的繁忙概率不超過其平均值的5%。

圖2 不同部署地點的A型機(jī)繁忙概率Fig.2 Busy probability of type A aircraft in different deployment locations

圖3 不同部署地點的B型機(jī)繁忙概率Fig.3 Busy probability of type B aircraft in different deployment locations
圖4和圖5為排隊模型求解算法得到的不同派遣順位對應(yīng)的校正系數(shù)QA(nA,ρA,k)和QB(nB,ρB,k),其中nA=10,nB=4,ρA和ρB的計算見排隊模型求解算法的步驟2。

圖4 A型機(jī)不同派遣順位對應(yīng)的校正系數(shù)Fig.4 Correction coefficient of different dispatch order for type A aircraft

圖5 B型機(jī)不同派遣順位對應(yīng)的校正系數(shù)Fig.5 Correction coefficient of different dispatch order for type B aircraft
從圖4和圖5中可以看出QA(nA,ρA,0)=QB(nB,ρB,0)=1,這是因為當(dāng)所有反潛機(jī)空閑時,能夠派遣反潛機(jī)的概率為1。
圖6統(tǒng)計了TA類別警報發(fā)出后,分派A型機(jī)處理時得到的目標(biāo)函數(shù)值,即TA警報的覆蓋率。在這個仿真實驗中,B型機(jī)數(shù)量固定為4,即nB=4,改變A型機(jī)數(shù)量nA,從中分析當(dāng)增加或減少A型機(jī)數(shù)量時,調(diào)度算法的性能變化情況,除此以外,圖6中還對警報發(fā)生率不同時,隨A型機(jī)數(shù)量變化的算法性能進(jìn)行了敏感性分析。

圖6 不同警報發(fā)生率下的目標(biāo)函數(shù)值敏感性分析Fig.6 Sensitivity analysis of objective function value under different alarm rates
圖6的結(jié)果可以用于幫助指揮人員判別需要多少數(shù)量的A型機(jī)才能達(dá)到指定的監(jiān)測覆蓋率指標(biāo)。例如,當(dāng)警報發(fā)生率λ介于2.5和3.5之間時,如果需要在警戒海域達(dá)到90%以上的警報覆蓋率,則需要部署6架A型機(jī)即可滿足要求,但如果想在同樣的警報發(fā)生率情況下達(dá)到95%的覆蓋率,A型機(jī)數(shù)量至少需要達(dá)到9架。而且從圖6中還可以發(fā)現(xiàn)隨著A型機(jī)數(shù)量增加,目標(biāo)函數(shù)值的增加逐漸減緩。例如當(dāng)λ=3.0時,A型機(jī)數(shù)量從2架增加到5架時,覆蓋性能呈現(xiàn)一個較大的增長,從53%增加到了90%,而當(dāng)A型機(jī)數(shù)量從8架增加到10架時,性能增長幅度較小,從94%增加到95%。
下面比較本文提出的優(yōu)化算法與常規(guī)“就近派遣”策略的性能。當(dāng)有警報發(fā)出后,就近派遣策略按離警報區(qū)域的距離由近到遠(yuǎn)調(diào)度,很明顯該策略的響應(yīng)時間最短,因此覆蓋率也相應(yīng)最大。但警報率較高時可能導(dǎo)致反潛機(jī)處于繁忙狀態(tài),這時需要根據(jù)全局態(tài)勢從長考慮,合理的做法可能是預(yù)留最近的反潛機(jī)以備更高等級的警報。其中就近派遣策略的實現(xiàn)步驟如下。
步驟 1根據(jù)候選部署點與85個監(jiān)測區(qū)之間的平均距離,利用分簇排序方法選擇了A型機(jī)和B型機(jī)的部署點;
步驟 2按照監(jiān)測區(qū)與部署點的距離遠(yuǎn)近生成排序調(diào)度表;
步驟 3在步驟1和步驟2的基礎(chǔ)上,設(shè)置2套決策變量hAl、hBl、mAlztk和mBlztk的相應(yīng)值為1,按式(7)計算覆蓋率。
圖7顯示了優(yōu)化算法與就近策略的警報覆蓋率性能指標(biāo),其中B型機(jī)數(shù)量設(shè)置為4,警報率λ=3,如圖7所示,就近派遣策略也能取得不錯的性能,但總體上優(yōu)化算法比就近派遣策略占優(yōu),在反潛機(jī)數(shù)量較小時(nB=4,nA=2)就近派遣策略與優(yōu)化算法性能相同。進(jìn)一步比較當(dāng)警報率變化時,從λ=1增加到10,優(yōu)化算法與就近派遣的性能差異。

圖7 就近派遣與優(yōu)化算法的性能對比(不同A型機(jī)數(shù)量)Fig.7 Performance comparison between the closest dispatch policy and optimal algorithm (different number of type A aircraft)
同樣地就近派遣策略的警報覆蓋率也較高,圖8中可以看出優(yōu)化算法最大高于就近派遣1.71%,仿真中參數(shù)設(shè)置為nB=4,nA=10。盡管就近派遣策略在警報覆蓋率上與優(yōu)化算法的差異不大,但優(yōu)化算法的優(yōu)勢在于構(gòu)建排序調(diào)度表時避免了局部反潛機(jī)過忙的影響。例如在某些情況下(指揮人員憑經(jīng)驗難以識別),派遣較遠(yuǎn)的反潛機(jī)與就近派遣相比,可以在整體反潛性能不變的情況下,降低反潛機(jī)的繁忙率,從而延長整個反潛機(jī)系統(tǒng)的生命周期。同時,通過把任務(wù)盡可能平均分配到所有反潛機(jī)上,如圖2和圖3所示,可以提高反潛機(jī)的空閑率,從而較少警報區(qū)域的等待時間,最終增加了警報的期望覆蓋率。而且優(yōu)化算法能夠產(chǎn)生一個合理的排序調(diào)度表,可以幫助指揮人員判別合適的預(yù)備反潛機(jī)。

圖8 就近派遣與優(yōu)化算法的性能(λ不同)Fig.8 Performance comparison between the closest dispatch policy and optimal algorithm (different λ)
當(dāng)海上傳感器平臺監(jiān)測到水下可疑目標(biāo)并發(fā)出警報后,如何快速有效地派遣反潛機(jī)至警報區(qū)域是決定反潛成敗的關(guān)鍵問題,如果響應(yīng)過慢則會丟失目標(biāo)而貽誤戰(zhàn)機(jī)。為此,本文提出了一個二元整數(shù)規(guī)劃模型及相應(yīng)求解算法,其主要目標(biāo)包括:① 解決目前兩種主力反潛直升機(jī)的優(yōu)化部署問題;② 構(gòu)造排序調(diào)度表以解決反潛機(jī)的分派問題。其中,如何利用有限的反潛機(jī)資源警戒大范圍海域是本文的一個難點,這也是應(yīng)召反潛的重點問題。為此,建立了一個排隊模型計算直升機(jī)繁忙概率和無關(guān)性校正系數(shù),以此作為規(guī)劃模型的輸入?yún)?shù),以便更加準(zhǔn)確地反映當(dāng)前態(tài)勢,提高規(guī)劃模型的精確度。
盡管該模型以兩種反潛機(jī)為例進(jìn)行說明,但可以利用相似方法擴(kuò)展到多類型反潛機(jī)優(yōu)化調(diào)度上。當(dāng)多處同時發(fā)現(xiàn)目標(biāo)需要派遣多機(jī)反潛時,也可以考慮應(yīng)用本文研究成果進(jìn)行進(jìn)一步擴(kuò)展。另外,警報類別或優(yōu)先級的劃分對調(diào)度模型的性能影響較大,盡管當(dāng)前有非常多的有關(guān)類別或優(yōu)先級計算的模型和算法,但尚缺乏反潛領(lǐng)域應(yīng)用的相關(guān)研究,這也是下一步工作的重點方向。