灾情巡视最优路线的寻径算法

刘长河

(北京建筑大学 理学院, 北京 100044)

1 问题的提出

我国是一个自然灾害多发的国度. 1998年长江中下游特大洪水, 2008年汶川大地震, 都给灾区人民的生命财产带来了严重的危害. 虽然人类永远无法制止自然灾害的发生, 但可以最大限度地降低灾难所造成的损失. 除了努力提高科技水平, 增强预测灾难的能力, 做好灾前避险以外, 灾后救援也尤为重要. 如何安排行程, 以最快、最有效的方式将救灾人员、救灾物资运送到受灾现场, 无疑是救灾工作的重要一环.

1998年, 全国大学生数学建模竞赛[1]提出了一个数学模型:

图1为某县的乡(镇)、村公路网示意图，公路边的数字为该路段的公里数.

今年夏天该县遭受水灾. 为考察灾情、组织自救，县领导决定，带领有关部门负责人到全县各乡(镇)、村巡视. 巡视路线指从县政府所在地出发，走遍各乡(镇)、村，又回到县政府所在地的路线.

题目假定: 在各乡(镇)停留时间T=2 h, 在各村停留时间t=1 h，汽车行驶速度V=35 km/h.

原题要解决的问题之一是: 巡视人员分为三组(路), 设计出总路程最短且各组尽可能均衡的巡视路线.

这是一类具有代表性、非常重要的问题. 在本文中, 将路程换算成行驶时间, 不但使问题更具合理性, 而且得到简化. 三个巡视组都从县城出发, 分别对各乡(镇)、村巡视一次, 最后回到县城. 要求总路程(时间)尽量短, 且各组任务尽量均衡, 即各组巡视(工作)的时间大致相当.

本文将围绕此目的完成以下三方面的任务: 研究解决此类问题的寻径算法；用Matlab语言编写出相应的程序；以上述模型进行数值实验, 验证本文算法和程序的正确性.

图1 某县的乡(镇)、村公路网示意图Fig.1 Schematic diagram of town and village highway network in a county

2 算法基本思路

为方便起见, 对县城、乡(镇)、村进行重新编号. 如图1所示, 县城编号为1(停留时间为0小时),从2到18号节点为乡(镇)所在地(停留时间为2 h), 19至53号节点表示村(停留时间为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.

其他节点之间没有边直接相连，边长无穷大，在编程时，可用一个很大的数dashu(如108)表示.

将两节点之间的路程用时间表示，即将邻接矩阵D的每个元素除以车速(V=35 km/h).

每个巡视组都从节点1(县城)出发,最后回到县城,要求将每个乡(镇)、村巡视一遍. 此问题分为以下三步完成:

第一步: 每个巡视组分别从某一节点出发,巡视若干节点,返回到出发节点. 每组各自形成一“圈”,各圈没有共同节点,每个节点都在某一个圈上.

这一步实际上是分配巡视任务. 如果图1本身的连通性较差,圈中相邻两节点之间可能没有边相连接(边长为dashu). 在这种情况下,利用Dijkstra算法寻找这对节点间的最短路径,将这两个节点连接起来. 最短路径上的内节点,巡视组只是路过,而不做停留. 为保证每组巡视任务的均衡,可以限定每个巡视组至少巡视若干节点.

第二步: 在第一步产生的圈中,必有一个包含节点1(县城),而其他圈都不包含节点1(县城).

对于包含节点1(县城)的圈,重排其节点顺序,形成一个起自县城,终于县城的回路.

对于不包含节点1(县城)的圈,首先利用Dijkstra算法找出自节点1到圈上各节点中距离最短的那个节点i0,相应的路径s_tour. 从节点1到沿s_tour节点i0，绕行一周，再从节点i0或在此圈上与节点i0相邻的某个节点回到节点1(沿最短路径)，完成回路(图2).

图2 从县城出发回到县城巡视线路示意图Fig.2 Schematic map of the inspection route from county to county

第三步:输出相应结果.

实现上述过程的第一步是一个旅行商问题(TSP).s个人,游历n个地点,每个地点恰被某一个人游历一次,要求每人至少游历m个地点. 此步可以归结为如下数学模型:

(1)

其中：

(2)

为第k个回路上行程和停留时间的总和. 如果选择所有节点处的停留时间都为0,它就是传统意义上的旅行商问题了.

TSP是一个NP问题,可采用遗传算法[2]进行计算.

遗传算法是一种模仿生物进化过程的随机算法. 它利用编码技术,将所求问题的解集当作一个种群,通过选择、交叉、变异等遗传操作，不断地使解集得到优化. 由于它具有全局寻优能力,适应性、鲁棒性强，过程简单等优点，在生产调度、自动控制、图像处理等许多领域都得到了广泛的应用.

基本遗传算法的基本步骤如下:

Step 1:按确定的编码机制进行编码,随机产生初始种群. 编码是一个十分复杂的过程,通常采用比较简单的二进制编码.

Step 2:确定个体的适应度(通常用轮盘赌策略),并判断是否符合优化准则,以确定是否可以终止循环并输出结果.

Step 3:选择适应度高的个体进入下一代种群.

Step 4:按照一定的交叉概率pc和交叉方法产生新个体.

Step 5:按照一定的变异概率pm和交叉方法产生新个体.

Step 6:由Step2～Step 5产生新个体产生下一代种群,返回Step 2.

遗传算法是求解旅行商问题的有效算法. 许多学者针对不同的TSP模型,选择不同的编码技术、适应度,交叉、变异方式，设计出许多不同的遗传算法[3-5].

3 灾情巡视路线的遗传算法与数值实验

为求解上述模型的最佳路径,本文在遗传算法过程的各个步骤采取了下述具体策略.

Step 1:利用城市标号进行编码. 按照所经过城市的先后顺序,将其标号排列起来,形成一个有序数组,用来表示巡视路径.

为保证s_group 个巡视组中，每个巡视组所要巡视的乡(镇)、村相对均衡,本文设定每个巡视组最少要巡视min_tour个乡(镇)、村. 将所有节点随机排成一排,随机产生s_group-1个断点,将一排节点分成s_group段,每一段至少包含min_tour个节点.

Step 2:选择每个回路上行程和停留时间的总和为适应度，如式(2)所示. 由于很难预估s_group个回路上适应度之和的最小值,本文预设一个较大的迭代步数num_iter. 当迭代步数达到num_iter时，终止迭代. 当前的路径、适应度认作最优.

Step 3:选择适应度高的个体进入下一代种群.

Step 4:交叉: 随机选取两个(介于1和节点数n之间的)整数I、J，对介于I、J之间的标号进行左右轮换，再随机选取新的分段点. 产生新路径,进入下一代种群.

Step 5:交叉: 随机选取两个(介于1和节点数n之间的)整数I、J，交换I、J两个位置上的标号，再随机选取新的分段点. 产生新路径,进入下一代种群.

Step 6:由Step2～Step 5产生新个体产生下一代种群,返回Step 2.

利用Matlab R2017 编程[6]，对图1 所示模型进行了计算. 组数s_group=3,每组至少巡视节点数min_tour=15. 经过num_iter=105代遗传运算,所得结果如下:

第一步: 每个巡视组所要巡视的节点圈.

第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巡视组的巡视任务圈. 其中两对节点7和42、15和29之间没有边直接连接,利用Dijkstra算法分别找出它们之间的最短路径. 比如,节点7和42之间的最短路径:

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

有下划线的节点本巡视组只是路过,巡视任务由其他巡视组完成. 其他类同,不再赘述.

第3组的巡视圈是畅通的,即从节点3出发,回到节点3.

在利用遗传算法寻找各个巡视组的节点圈时,各组的时间总和的变化规律如图3所示:

图3 各组时间总和变化规律Fig.3 Change rule of the sum of each group’s tour time

第二步:

不包含县城的路线圈为序列(3)、序列(4),经上述计算,分别得出从县城(节点1)出发,回到县城(节点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). 重新排序,形成从县城(节点1)出发,回到县城(节点1)的最终路线(第3组):

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

(8)

s_group个巡视组所需巡视的节点及相应时间见表1.

表1 各巡视组相关数据比较

4 结论

从理论上讲,灾情巡视属于旅行商问题(TSP). 这是一种NP问题，是图论研究中的难题之一. 本文将智(遗传)能算法与Dijkstra算法结合起来,克服了由于交通网络图本身连通性差对寻径问题所带来的困难,在算法设计的思路上有所突破. 本文算法灵活性强,只要调整巡视组的数目s_group,每组的最小任务数目min_tour,可以产生满足其他不同需求的方案. 从数值试验结果可以看出,在保证三个巡视组总时间最少的前提下,各组花费在道路上的时间和工作时间之和也基本均衡. 此算法的设计理念对处理物流货物的配送、超市商品的货架安排等问题，具有重要的借鉴作用.