劉長河

(北京建筑大學(xué) 理學(xué)院, 北京 100044)

1 問題的提出

我國是一個自然災(zāi)害多發(fā)的國度. 1998年長江中下游特大洪水, 2008年汶川大地震, 都給災(zāi)區(qū)人民的生命財產(chǎn)帶來了嚴(yán)重的危害. 雖然人類永遠(yuǎn)無法制止自然災(zāi)害的發(fā)生, 但可以最大限度地降低災(zāi)難所造成的損失. 除了努力提高科技水平, 增強(qiáng)預(yù)測災(zāi)難的能力, 做好災(zāi)前避險以外, 災(zāi)后救援也尤為重要. 如何安排行程, 以最快、最有效的方式將救災(zāi)人員、救災(zāi)物資運(yùn)送到受災(zāi)現(xiàn)場, 無疑是救災(zāi)工作的重要一環(huán).

1998年, 全國大學(xué)生數(shù)學(xué)建模競賽[1]提出了一個數(shù)學(xué)模型:

圖1為某縣的鄉(xiāng)(鎮(zhèn))、村公路網(wǎng)示意圖，公路邊的數(shù)字為該路段的公里數(shù).

今年夏天該縣遭受水災(zāi). 為考察災(zāi)情、組織自救，縣領(lǐng)導(dǎo)決定，帶領(lǐng)有關(guān)部門負(fù)責(zé)人到全縣各鄉(xiāng)(鎮(zhèn))、村巡視. 巡視路線指從縣政府所在地出發(fā)，走遍各鄉(xiāng)(鎮(zhèn))、村，又回到縣政府所在地的路線.

題目假定: 在各鄉(xiāng)(鎮(zhèn))停留時間T=2 h, 在各村停留時間t=1 h，汽車行駛速度V=35 km/h.

原題要解決的問題之一是: 巡視人員分為三組(路), 設(shè)計出總路程最短且各組盡可能均衡的巡視路線.

這是一類具有代表性、非常重要的問題. 在本文中, 將路程換算成行駛時間, 不但使問題更具合理性, 而且得到簡化. 三個巡視組都從縣城出發(fā), 分別對各鄉(xiāng)(鎮(zhèn))、村巡視一次, 最后回到縣城. 要求總路程(時間)盡量短, 且各組任務(wù)盡量均衡, 即各組巡視(工作)的時間大致相當(dāng).

本文將圍繞此目的完成以下三方面的任務(wù): 研究解決此類問題的尋徑算法；用Matlab語言編寫出相應(yīng)的程序；以上述模型進(jìn)行數(shù)值實(shí)驗(yàn), 驗(yàn)證本文算法和程序的正確性.

圖1 某縣的鄉(xiāng)(鎮(zhèn))、村公路網(wǎng)示意圖Fig.1 Schematic diagram of town and village highway network in a county

2 算法基本思路

為方便起見, 對縣城、鄉(xiāng)(鎮(zhèn))、村進(jìn)行重新編號. 如圖1所示, 縣城編號為1(停留時間為0小時),從2到18號節(jié)點(diǎn)為鄉(xiāng)(鎮(zhèn))所在地(停留時間為2 h), 19至53號節(jié)點(diǎn)表示村(停留時間為1 h). 寫出圖1的鄰接矩陣D.

矩陣D的主對角元素:

D(i,i)=0, (i=1,2,…,53);

D(1,4)=11.5,D(1,14)=19.8,D(1,16)=10.1,D(1,18)=12.9,D(1,19)=6.0,D(1,20)=9.2;

D(2,3)=12.2,D(2,18)=8.8,D(2,19)=10.3,D(2,51)=7.4,D(2,52)=11.5;

D(3,2)=12.2,D(3,4)=11.0,D(3,19)=5.9,D(3,52)=17.6;

D(4,1)=11.5,D(4,3)=11.0,D(4,19)=11.2,D(4,21)=7.9;

D(5,21)=8.2,D(5,22)=12.7,D(5,23)=11.3,D(5,25)=15.1;

D(6,25)=7.2,D(6,26)=8.0,D(6,27)=7.8,D(6,29)=14.2;

D(7,27)=5.6,D(7,28)=10.8,D(7,30)=12.2;

D(8,29)=6.8,D(8,30)=7.8,D(8,31)=8.6;

D(9,30)=10.2,D(9,32)=9.9;

D(10,11)=15.8;D(10,31)=16.4;D(10,33)=8.8;D(10,34)=11.8;D(10,34)=11.8,D(10,36)=8.2;

D(11,10)=15.8,D(11,29)=13.2,D(11,31)=9.8,D(11,36)=8.2,D(11,37)=8.1;

D(12,35)=9.8,D(12,36)=9.2,D(12,39)=4.1,D(12,40)=10.1;

D(13,24)=11.8,D(13,25)=14.5,D(13,37)=7.2,D(13,38)=5.5;

D(14,1)=19.8,D(14,15)=14.2,D(14,23)=11.4,D(14,24)=9.5,D(14,43)=12.0;

D(15,14)=14.2,D(15,41)=7.9,D(15,42)=13.2,D(15,43)=8.8,D(15,44)=10.5;

D(16,1)=10.1,D(16,44)=10.5,D(16,46)=12.1,D(16,47)=15.2;

D(17,46)=8.3,D(17,47)=7.2,D(17,48)=7.7,D(18,1)=12.9,D(18,2)=8.8,

D(18,47)=7.9,D(18,49)=9.2;

D(19,1)=6.0,D(19,2)=10.3,D(19,3)=5.9,D(19,4)=11.2;

D(20,1)=9.2,D(20,21)=4.8,D(20,23)=8.3;

D(21,4)=7.9,D(21,5)=8.2,D(21,20)=4.8;

D(22,5)=12.7,D(22,26)=20.4;

D(23,5)=11.3,D(23,14)=11.4,D(23,20)=8.3,D(23,24)=9.7;

D(24,13)=11.8,D(24,14)=9.5,D(24,23)=9.7,D(24,25)=7.3;

D(25,5)=15.1,D(25,6)=7.2,D(25,13)=14.5,D(25,24)=7.3;

D(26,6)=8.0,D(26,22)=20.4;

D(27,6)=7.8,D(27,7)=5.6;

D(28,7)=10.8;

D(29,6)=14.2,D(29,8)=6.8,D(29,11)=13.2;

D(30,7)=12.2,D(30,8)=7.8,D(30,9)=10.2;

D(31,8)=8.6,D(31,10)=16.4,D(31,11)=9.8,D(31,32)=8.6;

D(32,9)=9.9,D(32,31)=8.6,D(32,33)=15;

D(33,10)=8.8,D(33,32)=15;

D(34,10)=11.8,D(34,35)=6.8;

D(35,12)=9.8,D(35,34)=6.8,D(35,40)=6.7;

D(36,10)=8.2,D(36,11)=8.2,D(36,12)=9.2;

D(37,11)=8.1,D(37,13)=7.2,D(37,38)=9.3;

D(38,13)=5.5,D(38,37)=9.3,D(38,39)=7.9,D(38,43)=6.5;

D(39,12)=4.1,D(39,38)=7.9,D(39,41)=9.1,D(39,43)=7.8;

D(40,12)=10.1,D(40,35)=6.7,D(40,41)=10.0;

D(41,15)=7.9,D(41,39)=9.1,D(41,40)=10.0,D(41,42)=8.9;

D(42,15)=13.2,D(42,41)=8.9,D(42,45)=18.8;

D(43,14)=12.0,D(43,15)=8.8,D(43,38)=6.5,D(43,39)=7.8;

D(44,15)=10.5,D(44,16)=10.5,D(44,45)=7.8;

D(45,42)=18.8,D(45,44)=7.8,D(45,46)=7.9;

D(46,16)=12.1,D(46,17)=8.3,D(46,45)=7.9;

D(47,16)=15.2,D(47,17)=7.2,D(47,18)=7.9;

D(48,17)=7.7,D(48,50)=10.3;

D(49,18)=9.2,D(49,50)=8.1,D(49,51)=7.3;

D(50,48)=10.3,D(50,49)=8.1,D(50,51)=19.0,D(50,53)=14.9;

D(51,2)=7.4,D(51,49)=7.3,D(51,50)=19.0,D(51,53)=20.3;

D(52,2)=11.5,D(52,3)=17.6,D(52,53)=8.2;

D(53,50)=14.9,D(53,51)=20.3,D(53,52)=8.2.

其他節(jié)點(diǎn)之間沒有邊直接相連，邊長無窮大，在編程時，可用一個很大的數(shù)dashu(如108)表示.

將兩節(jié)點(diǎn)之間的路程用時間表示，即將鄰接矩陣D的每個元素除以車速(V=35 km/h).

每個巡視組都從節(jié)點(diǎn)1(縣城)出發(fā),最后回到縣城,要求將每個鄉(xiāng)(鎮(zhèn))、村巡視一遍. 此問題分為以下三步完成:

第一步: 每個巡視組分別從某一節(jié)點(diǎn)出發(fā),巡視若干節(jié)點(diǎn),返回到出發(fā)節(jié)點(diǎn). 每組各自形成一“圈”,各圈沒有共同節(jié)點(diǎn),每個節(jié)點(diǎn)都在某一個圈上.

這一步實(shí)際上是分配巡視任務(wù). 如果圖1本身的連通性較差,圈中相鄰兩節(jié)點(diǎn)之間可能沒有邊相連接(邊長為dashu). 在這種情況下,利用Dijkstra算法尋找這對節(jié)點(diǎn)間的最短路徑,將這兩個節(jié)點(diǎn)連接起來. 最短路徑上的內(nèi)節(jié)點(diǎn),巡視組只是路過,而不做停留. 為保證每組巡視任務(wù)的均衡,可以限定每個巡視組至少巡視若干節(jié)點(diǎn).

第二步: 在第一步產(chǎn)生的圈中,必有一個包含節(jié)點(diǎn)1(縣城),而其他圈都不包含節(jié)點(diǎn)1(縣城).

對于包含節(jié)點(diǎn)1(縣城)的圈,重排其節(jié)點(diǎn)順序,形成一個起自縣城,終于縣城的回路.

對于不包含節(jié)點(diǎn)1(縣城)的圈,首先利用Dijkstra算法找出自節(jié)點(diǎn)1到圈上各節(jié)點(diǎn)中距離最短的那個節(jié)點(diǎn)i0,相應(yīng)的路徑s_tour. 從節(jié)點(diǎn)1到沿s_tour節(jié)點(diǎn)i0，繞行一周，再從節(jié)點(diǎn)i0或在此圈上與節(jié)點(diǎn)i0相鄰的某個節(jié)點(diǎn)回到節(jié)點(diǎn)1(沿最短路徑)，完成回路(圖2).

圖2 從縣城出發(fā)回到縣城巡視線路示意圖Fig.2 Schematic map of the inspection route from county to county

第三步:輸出相應(yīng)結(jié)果.

實(shí)現(xiàn)上述過程的第一步是一個旅行商問題(TSP).s個人,游歷n個地點(diǎn),每個地點(diǎn)恰被某一個人游歷一次,要求每人至少游歷m個地點(diǎn). 此步可以歸結(jié)為如下數(shù)學(xué)模型:

(1)

其中：

(2)

為第k個回路上行程和停留時間的總和. 如果選擇所有節(jié)點(diǎn)處的停留時間都為0,它就是傳統(tǒng)意義上的旅行商問題了.

TSP是一個NP問題,可采用遺傳算法[2]進(jìn)行計算.

遺傳算法是一種模仿生物進(jìn)化過程的隨機(jī)算法. 它利用編碼技術(shù),將所求問題的解集當(dāng)作一個種群,通過選擇、交叉、變異等遺傳操作，不斷地使解集得到優(yōu)化. 由于它具有全局尋優(yōu)能力,適應(yīng)性、魯棒性強(qiáng)，過程簡單等優(yōu)點(diǎn)，在生產(chǎn)調(diào)度、自動控制、圖像處理等許多領(lǐng)域都得到了廣泛的應(yīng)用.

基本遺傳算法的基本步驟如下:

Step 1:按確定的編碼機(jī)制進(jìn)行編碼,隨機(jī)產(chǎn)生初始種群. 編碼是一個十分復(fù)雜的過程,通常采用比較簡單的二進(jìn)制編碼.

Step 2:確定個體的適應(yīng)度(通常用輪盤賭策略),并判斷是否符合優(yōu)化準(zhǔn)則,以確定是否可以終止循環(huán)并輸出結(jié)果.

Step 3:選擇適應(yīng)度高的個體進(jìn)入下一代種群.

Step 4:按照一定的交叉概率pc和交叉方法產(chǎn)生新個體.

Step 5:按照一定的變異概率pm和交叉方法產(chǎn)生新個體.

Step 6:由Step2～Step 5產(chǎn)生新個體產(chǎn)生下一代種群,返回Step 2.

遺傳算法是求解旅行商問題的有效算法. 許多學(xué)者針對不同的TSP模型,選擇不同的編碼技術(shù)、適應(yīng)度,交叉、變異方式，設(shè)計出許多不同的遺傳算法[3-5].

3 災(zāi)情巡視路線的遺傳算法與數(shù)值實(shí)驗(yàn)

為求解上述模型的最佳路徑,本文在遺傳算法過程的各個步驟采取了下述具體策略.

Step 1:利用城市標(biāo)號進(jìn)行編碼. 按照所經(jīng)過城市的先后順序,將其標(biāo)號排列起來,形成一個有序數(shù)組,用來表示巡視路徑.

為保證s_group 個巡視組中，每個巡視組所要巡視的鄉(xiāng)(鎮(zhèn))、村相對均衡,本文設(shè)定每個巡視組最少要巡視min_tour個鄉(xiāng)(鎮(zhèn))、村. 將所有節(jié)點(diǎn)隨機(jī)排成一排,隨機(jī)產(chǎn)生s_group-1個斷點(diǎn),將一排節(jié)點(diǎn)分成s_group段,每一段至少包含min_tour個節(jié)點(diǎn).

Step 2:選擇每個回路上行程和停留時間的總和為適應(yīng)度，如式(2)所示. 由于很難預(yù)估s_group個回路上適應(yīng)度之和的最小值,本文預(yù)設(shè)一個較大的迭代步數(shù)num_iter. 當(dāng)?shù)綌?shù)達(dá)到num_iter時，終止迭代. 當(dāng)前的路徑、適應(yīng)度認(rèn)作最優(yōu).

Step 3:選擇適應(yīng)度高的個體進(jìn)入下一代種群.

Step 4:交叉: 隨機(jī)選取兩個(介于1和節(jié)點(diǎn)數(shù)n之間的)整數(shù)I、J，對介于I、J之間的標(biāo)號進(jìn)行左右輪換，再隨機(jī)選取新的分段點(diǎn). 產(chǎn)生新路徑,進(jìn)入下一代種群.

Step 5:交叉: 隨機(jī)選取兩個(介于1和節(jié)點(diǎn)數(shù)n之間的)整數(shù)I、J，交換I、J兩個位置上的標(biāo)號，再隨機(jī)選取新的分段點(diǎn). 產(chǎn)生新路徑,進(jìn)入下一代種群.

Step 6:由Step2～Step 5產(chǎn)生新個體產(chǎn)生下一代種群,返回Step 2.

利用Matlab R2017 編程[6]，對圖1 所示模型進(jìn)行了計算. 組數(shù)s_group=3,每組至少巡視節(jié)點(diǎn)數(shù)min_tour=15. 經(jīng)過num_iter=105代遺傳運(yùn)算,所得結(jié)果如下:

第一步: 每個巡視組所要巡視的節(jié)點(diǎn)圈.

第1組：

20→23→14→24→25→6→27→7→28→7→27→6→25→24→14→15→42→41→15→14→ 24→25→6→26→22→5→21→20

(3)

第2組：

31→11→37→13→38→43→39→12→36→10→33→10→34→35→40→35→34→10→ 11→29→8→30→9→32→31

(4)

第3組：

3→4→1→18→47→16→44→45→46→17→48→50→49→51→53→52→2→19→3

(5)

序列(3)表示第1巡視組的巡視任務(wù)圈. 其中兩對節(jié)點(diǎn)7和42、15和29之間沒有邊直接連接,利用Dijkstra算法分別找出它們之間的最短路徑. 比如,節(jié)點(diǎn)7和42之間的最短路徑:

7→27→6→25→24→14→15→42.

有下劃線的節(jié)點(diǎn)本巡視組只是路過,巡視任務(wù)由其他巡視組完成. 其他類同,不再贅述.

第3組的巡視圈是暢通的,即從節(jié)點(diǎn)3出發(fā),回到節(jié)點(diǎn)3.

在利用遺傳算法尋找各個巡視組的節(jié)點(diǎn)圈時,各組的時間總和的變化規(guī)律如圖3所示:

圖3 各組時間總和變化規(guī)律Fig.3 Change rule of the sum of each group’s tour time

第二步:

不包含縣城的路線圈為序列(3)、序列(4),經(jīng)上述計算,分別得出從縣城(節(jié)點(diǎn)1)出發(fā),回到縣城(節(jié)點(diǎn)1)的最終路線. 其中序列(3)所形成的最終路線(第1組):

1→20→23→14→24→25→6→27→7→28→7→27→6→25→24→14→15→42→41→15→14→24→25→6→26→22→5→21→20→1.

(6)

序列(4)所形成的最終路線(第2組):

1→14→43→39→12→36→10→33→10→34→35→40→35→34→10→11→29→8→30→9→32→31→11→37→13→38→43→14→1

(7)

包含縣城的路線圈為序列(5). 重新排序,形成從縣城(節(jié)點(diǎn)1)出發(fā),回到縣城(節(jié)點(diǎn)1)的最終路線(第3組):

1→18→47→16→44→45→46→17→48→50→49→51→53→52→2→19→3→4→1

(8)

s_group個巡視組所需巡視的節(jié)點(diǎn)及相應(yīng)時間見表1.

表1 各巡視組相關(guān)數(shù)據(jù)比較

4 結(jié)論

從理論上講,災(zāi)情巡視屬于旅行商問題(TSP). 這是一種NP問題，是圖論研究中的難題之一. 本文將智(遺傳)能算法與Dijkstra算法結(jié)合起來,克服了由于交通網(wǎng)絡(luò)圖本身連通性差對尋徑問題所帶來的困難,在算法設(shè)計的思路上有所突破. 本文算法靈活性強(qiáng),只要調(diào)整巡視組的數(shù)目s_group,每組的最小任務(wù)數(shù)目min_tour,可以產(chǎn)生滿足其他不同需求的方案. 從數(shù)值試驗(yàn)結(jié)果可以看出,在保證三個巡視組總時間最少的前提下,各組花費(fèi)在道路上的時間和工作時間之和也基本均衡. 此算法的設(shè)計理念對處理物流貨物的配送、超市商品的貨架安排等問題，具有重要的借鑒作用.