谢必成,张小俊
河北工业大学 机械工程学院,天津 300401
随着智能机器人技术的不断发展,全覆盖路径规划(complete coverage path planning,CCPP)技术已经成为智能机器人领域的关键技术。近些年来由于其在军事[1]、农业[2]、搜救[3]、工业[4]、商业[5]和民用[6]等领域的广泛应用,进而得到越来越多的学者的关注和研究。对于工作效率更高的多机器人全覆盖路径规划(multirobot complete coverage path planning,MCCPP)的研究也已经得到学者们的广泛关注[7]。目前对MCCPP 的研究也大多是在平面环境下进行的,而对于石油储罐检测、船舶喷漆与除锈以及大型玻璃外表面清洁等工作环境更为复杂和危险的立面维护作业MCCPP的研究便更有其特定的研究价值。
目前,路径规划主要分为两大类:一类是在多个固定点之间寻找最短无碰撞路径的路径规划方法[8-10];另一类是在除不可行域以外的任务空间里寻找完全遍历的最大覆盖、最小长度的路径规划方法[11-12]。
然而,对于第一类路径规划方法一般归结为对旅行商问题(travelling salesman problem,TSP)的解决,TSP是确定对若干个点无重复访问顺序的典型NP-hard 问题,即从任一点出发并回到这一点,以欧几里德距离作为路径长度的度量方法寻找最短路径,然而多旅行商问题(multiple travelling salesman problem,MTSP)是经典TSP 的一种泛化[13],即同时解决多个TSP。而对于第二类问题的解决,MCCPP 是一种工作效率比较高的解决办法,可以看作是MTSP与覆盖路径规划(coverage path planning,CPP)的组合。目前MCCPP 可以大致分为以下几种:聚类法,如Tang 等人[14]采用改进的k-means 聚类方法对预先部署的传感器点进行分类,引入反馈机制来平衡每条路由的长度,最后通过求解TSP得到每个机器人的封闭路径来实现MCCPP,由于传感器点布置具有随机性,所以这种方法对于立面维护作业而言主要的局限性是路径重复率较高,无法充分考虑立面维护作业机器人起点和终点位置特点对覆盖工作的影响;栅格地图法,如Gao等人[15]利用基于机器人特定初始位置划分区域的算法将区域划分为多个相等的子区域,然后对各子区域之间交换特定的端点节点得到理想生成树从而实现MCCPP,由于生成树本身具有的栅格性,所以这种方法对于立面维护作业而言主要的局限性是对复杂环境和多工艺过程工作变现较差,同时不方便对感兴趣权重的设置;几何图法,如Nair 等人[16]使用基于机器人特定的初始位置并结合测地线-曼哈顿Voronoi 分区方法对任务空间划分,然后使用传统的Boustrophedon 和回溯式覆盖方法对每个区域覆盖来实现MCCPP,由于这种方法对机器人的初始位置具有较大的依赖性,所以这种方法主要局限性是不能充分考虑立面维护作业机器人起点位置和终点位置具有边缘性的特点;仿生法,如阮贵航等人[17]提出了一种基于滚动优化和分散捕食者猎物模型的多机器人全覆盖路径规划算法,构建分散捕食者猎物模型,同时引入滚动优化方法,避免机器人陷入局部最优,最终实现MCCPP,这种方法的局限性主要表现在对多机器人工作均衡性无法保证,以及不能很好地适应复杂障碍物环境。
本文主要针对以上MCCPP方法在立面维护作业场景下存在的局限性,提出了使用BCD[18](boustrophedon cellular decomposition)将任务空间分解成若干个Cell,然后基于免疫遗传算法对MCCPP问题求解。该方法一方面是要解决每个Cell的CPP 问题,另一方面还要解决对若干个Cell的访问顺序的MTSP问题,即在确定每一个Cell的路径入口点和出口点组合的同时还要确定对所有Cell的访问顺序。该方法主要完成目标是:结合立面维护作业工艺和环境特征,解决任务空间边缘任意起始点和终点的MTSP-CPP 问题。在最大程度地均衡每个机器覆盖任务的前提下,实现工作时间最短以及总工作量最小的多求解目标,以便充分发挥立面作业MCCPP的优势。
在开始问题研究之前,需要声明几点前提条件和几点假设。
首先,本文提出几点前提:
(1)本文所使用的地图是二维平面的先验地图,所以对如何构建地图不做讨论,只对覆盖方法的可行性进行验证。
(2)本文所使用的地图对不可行域已经进行了大小为机器人半径R的膨胀处理(覆盖时不考虑机器人与不可行域重合的问题),同时认为机器人不可跨越不可行域。
其次,本文提出几点假设:
(1)假设所使用的机器人具有无限的续航能力。
为了研究跟专注于全覆盖路径规划算法的研究,假设机器人拥有彼此防撞机制。
(2)假设每一个机器人为一个质点,其单位覆盖范围为是一个以质点为中心、半径为R的圆形。
忽略机器人转弯运动和直线运动的差异性,即不考虑机器人转弯运动的时间和能量损失。
本文将采用文献[18]中提到的BCD 方法对任务空间进行分解。同时本文将使用文献[19]中提到的Cell内部遍历路径模型与Cell间过渡路径模型,目的在于将多机器人立面维护作业全覆盖问题的模型进行简化。在此,需要对本文所使用的BCD 方法和下文将要使用的路径模型做出一定的解释。与传统的BCD方法中依靠障碍物分解方法相比,本文将“障碍物”这一概念引申为“不可行域”,设置不可行域的目的是:除了简单的障碍物区域不可到达以外,由于立面环境下存在不同材质或者装饰物的原因而被语义识别为具有保护性质的不可行域。例如喷漆与除锈混合工艺过程无法同时空进行[20],所以出于对上一级工艺成果的保护而被识别为不可行域。为了方便工艺调度以及为不同作业区域分配不同的时间权重,而将任务空间划分为若干个Cell。为了最大程度均衡每个机器人的覆盖任务,同时保证机器人起点与终点位置选择的灵活性,进而选用了下文介绍的Cell内部路径覆盖模型和Cell间过渡路径模型。基于以上描述,不难看出不可行区域的形状和位置具有任意性,遍历路径应具有灵活性,所以作者采用BCD分解方法分解任务空间,并使用下文中所介绍的路径模型对问题模型进行简化。
BCD 是依据不可行域的分布将任务空间分解成若干个Cell,然后机器人采用割草式往复运动对每一个Cell实现覆盖。分解过程中将既得的二值化不可行域地图视为一系列切片,所以为了方便算法实现,作者在设置中认为一个1 像素宽的垂直地图条为一个切片。这些切片或不被不可行域分割,或被一个乃至多个不可行域区域分割成非重叠的自由空间。约定将每一个切片自由空间数目记作当前切片的连通性计数。在算法实现方面,在沿着水平方向逐一扫描每一个切片的同时计算每一个切片的连通性。切片的连通性将会在临界点处发生变化。每当连通性增加,则触发一个IN事件,封闭当前Cell,并开启两个新的Cell。若遇到连通性降低,则触发一个OUT 事件,当前的两个Cell封闭,并开启一个新的Cell。据此,完成整个人空间的分解。如图1所示为简单情况下分解原理图。
图1 地图分解原理Fig.1 Principle of map decomposition
在任务空间通过BCD分解为N个Cell之后,分解得到的每一个Cell将被标记一个id(Cell的编号,记作Cellid)。每个机器人要完成的覆盖的任务由一组id序列组成,所有机器人覆盖任务共同排列起来组成一个长度为N的非重复id序列。由于预先把割草式往复运动作为Cell内部覆盖形式,所以接下来的工作就是没有遗漏地为机器人分配覆盖任务,并确定访问顺序。本文假设只能将一个Cell内部路径的出口点通过Cell间过渡路径与临近Cell的内部路径的入口点连接,当前Cell会通过两条过渡路径与前后两个临近Cell链接。所以,Cell间过渡路径起点和终点分别与前一的Cell的出口点和下一个Cell入口点重合。对于同一个Cell执行了不同割草式覆盖方向,会形成4对可选入口和出口点组合。
如果将每个Cell内部路径入口点和出口点组合固定且每个机器人的初始位置不同,那么每个机器人在各自任务空间中寻找最短访问路径问题就是解决TSP 问题,多个机器人在各自任务空间中寻找最短访问路径问题就变成了MTSP 问题。在解决这个问题之前必须设法寻找Cell间最短过渡路径问题的解,也就是当前Cell内部路径出口点作为起始点与下一个Cell内部路径入口点作为终点的最短路径,且要求这条路径与不可行域无重叠。为了解决这个问题可以利用不可行域顶点并把过渡路径的起点和终点作为附加点,构建可视图。以欧几里德距离作为过渡路径度量,通过使用Dijkstra 算法寻找从起点到终点的最短路径。
然而,对于具体到每个Cell的全覆盖问题,会有四条不同遍历方向的覆盖路径作为待选路径,具体如图2的举例所示。因此,不同的Cell内部路径的入口点和出口点组合的选择不仅对自身内部路径长路有影响(图2 中四种情况路径长度存在(a)=(b),(c)=(d)的关系),而且对具有固定访问顺序的Cell间过渡路径长度产生影响。为了方便理解,暂且忽略不可行域对总路径长度的影响,图3举例说明了不同的Cell访问顺序以及不同的Cell入口和出口点组合对总路径长度的影响。图3(a)中Cell访问顺序为1➝2➝3;图3(b)中Cell访问顺序与图3(a)中相同,但每个Cell入口和出口点组合不同,显然总路径长度也有所不同;图3(c)中Cell访问顺序与图3(a)中不同,访问顺序为3➝2➝1,但图3(c)与图3(a)的每个Cell入口和出口点组合相同,显然这样总路径长度也有所不同。
图2 Cell 的4种不同的入口和出口点组合Fig.2 4 different combinations of entry and exit points for Cell
图3 路径长度的影响因素Fig.3 Influencing factors of path length
承接前文的准备工作,本章将详细介绍用于解决具有MTSP 性质的全覆盖问题的免疫遗传算法的具体实现过程。算法的总体流程如下:
步骤1初始化种群。使用Knuth洗牌法,将任务空间分解得到的N个Cell随机排列,并为每个Cell随机确定一组入口点和出口点,这样得到一个个体。照此方法产生m个个体,m为预设的种群规模。
步骤2适应度计算。按照个体染色体顺序整体估算工作时间,按照每个机器人工作时间尽可能相等原则分割基因段。在考虑每个机器人初始位置情况下,优化机器人所访问Cell的入口和出口组合。计算每一个机器人从指定初始位置出发的工作时间,取出最大的工作时间与所有机器人工作时间之和的乘积的倒数作为每个个体的适应度。选出适应度值最高的个体。
步骤3按照基因移植概率,对整个种群重复进行优秀染色体片段移植操作。
步骤4对种群进行阶梯型克隆选择操作。
步骤5重复步骤3 和步骤4 操作,寻找每一次迭代后的适应度值最高的个体,不断更新最优个体。直到迭代预设代数后结束操作。按照步骤2 中适应度计算方法,输出最优个体的工作时间、所有机器人总工作时间以及染色体分割结果。
本文中影响每个机器人对所分配的子任务空间的覆盖效果的因素主要为以下两点:
(1)机器人对子任务空间中Cell的覆盖顺序。
(2)子任务空间中每个Cell入口点和出口点组合的选择。
因此,本文使用十进制整数编码基因,每一位基因gene要符合以上两点要素,具体编码规则如下:
式中,i为Cell的4种内部路径入口点和出口点组合编号,0、1、2、3 分别对应图2(a)、图2(b)、图2(c)、图2(d)这4种情况。
由式(1)不难得到解码规则如下:
式中,//为整除运算符,%为取余运算符。举例说明某个体染色体的编码与解码过程如图4 所示。图4(a)为编码过程,其中个体染色体第一位编码“23”是由“id”行第一位“6”与“i”行第一位“3”按照式(1)运算得到,即(6-1)×4+3=23;图4(b)为解码过程,其中“id”行最后一位“10”是由个体染色体最后一位编码“39”按照式(2)运算得到,即39//4+1=10,“i”行最后一位“3”是由个体染色体最后一位编码“39”按照式(2)运算得到,即39%4=3。
图4 编码与解码原理Fig.4 Principle of encoding and decoding
基因编码和解码规则介绍完后,下面介绍基于动态均衡每个机器人工作时间的原则计算个体适应度值。首先按照染色体中每一位基因解码得到的Cell的id和对应Cell内部路径入口点和出口点组合,计算覆盖整个任务空间中所有Cell内部路径所需要的总工作时间Ta:
式中,n表示个体染色体中基因的位置,μa表示Cell内部路径时间因子,la,n表示第n个位置上gene解码后计算得到内部路径长度。由于Ta远大于所有覆盖Cell间过渡路径Te,所以在分配任务前估算所有机器人总工作时间为Ta。那么,按式(4)所示计算得到Nrobot个机器人理想工作时间为tdream。然而通过BCD 分解得到的Cell面积存在不均匀性,为了均衡每个机器人的工作时间,仍然以图4 中的个体染色体为例,当Nrobot=3时,如图5所示举例说明染色体动态均匀分割方法。用线段长度表示时间长短,图中理想分割点就是tdream的具体表达,然后动态地取最靠近理想分割点的线段端点作为实际分割点,所得的(Nrobot-1) 个实际分割点把整个个体染色体尽可能均匀地分成Nrobot个子染色体段。
图5 染色体动态均分原理Fig.5 Principle of dynamic homogeneity of chromosomes
个体染色体被分割成多段染色体片段后就会分配给对应的机器人。如算法1所示利用动态规划算法,同时考虑立面维护作业机器人起点位置和终点位置的特征对每一个染色体段重新优化,目的在于机器人对各个Cell的覆盖顺序不变的情况下,合理调整每个Cell的入口点和出口点组合,以达到每个机器人总工作时间最短的目的。然后按照优化结果,使用式(5)计算个体适应度值。最后将所有染色体片段按原来顺序重新拼接成新的染色体,方便其他操作。
式中,ti表示第i个机器人从路径起点出发覆盖相应的子染色体段中的所有Cell并到达终点结束的工作时间;Fit表示个体适应度值;max{}⋅表示取最大值运算;lstart表示机器人从起点出发到达相应的子染色体段上第一个Cell入口点的路径长度;lend表示机器人从子染色体段上最后一个Cell出口点出发到达路径终点的路径长度;le,j表示从子染色体段上第j个Cell的出口点到下一个Cell入口点之间的过渡路径长度;la,j表示子染色体段上第j个Cell的内部路径长度;μe表示外部路径时间因子,μe≤μa。
算法1使用动态规划算法优化机器人固定访问顺序下覆盖路径最短的每个Cell的最佳入口和出口点组合
与传统的遗传算法中的顺序交叉算子产生两个新子代个体相比,本文中优秀染色体片段移植算子操作后产生一个新的个体offspring。将种群按照适应度的值从大到小排序,在种群的前25%优秀个体中随机选择个体parent1,再从种群中使用式(6)所示的个体被选中的概率与其适应度值成负相关关系的轮盘赌来选择个体parent2:
式中,Fiti表示第i个个体的适应度值;pi表示第i个个体被选中的概率;qi表示第i个个体的累积概率;这样设置轮盘赌的目的是使种群中的劣质个体更容易被选中。
然后,在parent1 中随机选择起始点切取自适应长度Sg的染色体片段,式(7)为Sg计算方法:
式中,int{⋅}为向下取整运算符。
当起点位置加上Sg超过parent1 长度时,从parent1的首位开始继续切取直到满足Sg长度。把从parent1中切取的染色体片段放到offspring的相同位置上,同时从parent2 中找到与offspring插入结束点相同的位置,从这个位置开始检索与从parent1 中切取的染色体片段不同的基因依次放置到offspring的空缺位置上,结束后用offspring替换掉种群中的parent2。最后,计算offspring适应度值并保存。这样类似于将parent1 染色体片段移植到parent2,所以把该算子称作优秀染色体片段移植算子。如图6所示举例说明优秀基因移植实现过程。
图6 优秀染色体片段移植Fig.6 Excellent chromosome fragment transplantation
图6中取出优质个体parent1 的连续染色体片段并放到offspring的相同位置上的原因在于gene的位置具有实际意义。将优秀个体的随机位上开始的一段染色体移植到劣质个体的相同位置上,使劣质个体显性表达了优质个体的部分优质特征从而产生了新的个体offspring。之所以控制从优质个体上切取染色体片段的长度,是因为优秀个体的优秀特征主要取决于一段染色体的连续编码顺序,但是移植的太长势必会降低种群的多样性,太短又失去了移植的意义,所以需要自适应控制Sg的值。之所以保留优质染色体片段在子代中的位置,原因在于多机器人任务分配与适应度计算是按照机器人起点和终点位置序列对个体染色体动态分割的,所以优质的连续编码安放的位置也十分重要。这样操作在某种程度上不仅保留种群多样性的,还不破坏优秀个体,更是对劣质个体的改良。此方法模仿生物族群进化原理,将种群中优质特征基因选择性遗传下来。微观上看,由于优质个体片段切取位置具有随机性,所以这种操作对劣质个体的启发能力是有限的。但从宏观上看,优质特征在种群中的传播,将增大劣质群体的适应度,结合下文中的阶梯型克隆选择算子,势必启发种群加大对改良后的个体的克隆选择的能力。在控制适当的优秀染色体片段移植发生的概率,在整体上将启发种群更快的得到全局近似解。
将种群按照适应度的值从大到小排序,对种群的前25%的个体克隆4×M(M为基础克隆倍数,根据个体染色体长度确定)份,在保留被克隆的个体同时,将4×(M-1)个个体的染色体中随机两个Cell调换位置,实现基因突变。计算4×M个体适应度值,将适应度值最高的个体更新为新的个体。对前25%~50%的个体按照3×M的倍数重复上面操作。对前50%~75%的个体按照2×M的倍数重复操作。对种群后25%的个体删除,并重新初始化相应数目的个体,计算适应度后放入种群。
这样操作的目的是在保证种群多样性的同时使算法获得较好的收敛效果。将种群按照适应值大小平均分成了4 组,模仿生物免疫理论,提高亲和能力强的分组被克隆选择的能力,降低适应能力差的分组被克隆选择的能力,甚至将其淘汰。将适应度值最小的一组淘汰并重新初始化新的个体加入种群,为了剔除劣质个体并保持种群多样性;按照适应度值越大克隆选择倍数越高的阶梯型法则,可以加快优质解向最优解靠近速度,同时降低由于过渡求解次优解而导致不必要的资源浪费。
多机器人全覆盖算法性能的优劣和所要覆盖的地图的复杂度有关,本文设计了3个不同复杂程度的像素地图对全覆盖算法仿真验证。地图复杂度包括地图大小、不可行域形状的任意性以及分解后Cell数目的多少。所有仿真实验中同质机器人覆盖半径R=10,种群规模m=2 048,基础克隆倍数M=160,优秀染色体片段移植发生概率0.02,内部路径时间因子μa=3,外部路径时间因子μe=1,为了方便观察算法收敛性取算法迭代代数为50。考虑到在立面维护作业的多机器人全覆盖作业过程中,机器人工作路径的起点和终点多是待覆盖区域的边缘,所以本文仿真实验所选取的机器人起点和终点都在待覆盖区域边缘。如图7(a)~(c)分别展示了算法在不同复杂程度的地图中的仿真结果。表1 展示了不同仿真地图的相关信息。
表1 不同复杂度程度的仿真数据Table 1 Data for simulations of different degrees of complexity
图7 不同复杂度程度的仿真结果Fig.7 Simulation results for different levels of complexity
为了检验算法的质量,本文使用5个同质机器人对较为复杂的图7(c)所示地图实现全覆盖。如图8所示,分别对每代最优个体的工作时长(所有机器人工作时长中最大值,即max{ti})、总工作量(所有机器人工作时长算术和,即、不均衡度U以及适应度值的倒数1/Fit的50 次迭代收敛情况分析。其中不均衡度计算方法如式(8)所示,用于评价多机器人全覆盖过程中任务分配的不均衡程度,U数值越大说明多机器人任务分配的越不均衡。
图8 算法收敛性效果曲线Fig.8 Curve of convergence effect of algorithm
式中,min{ti}表示取最小值运算。
观察图8 发现工作时长在1 至21 代整体呈现减小趋势,但在5至6代、19至20代略有增大,在第21代后收敛于39 287.5;总工作量在1至23代整体呈现减小趋势,但在6 至7 代略有增大,在第23 代后收敛于194 481;不均衡度在1 至23 代呈现较大波动,在23 代后稳定在0.004 438<1%;1/Fit在1 代至23 代持续减小,在第23代后收敛于7.640 69×109。由于工作时长和总工作量相互影响,但1/Fit持续减小,所以工作时长和总工作量在收敛过程略有波动属于正常现象,这样更有利于彼此跳出局部最优解。由于BCD 分解具有不均匀性,使得不均衡度在稳定之前呈现很大波动,但是前面介绍的算法已经考虑到任务分配的均衡问题,所以相比于不均衡度的收敛过程,关注不均衡度的稳定值更有意义。综上可得,本文所提出的用于解决立面维护作业的多机器人全覆盖路径规划的免疫遗传算法具有较好的收敛稳定性。
为了检验算法的收敛速度、运算复杂度、收敛精度。将本文的双目标适应度方法与传统的单目标适应度方法进行比较,使用1 至5 个同质机器人对较为复杂的图7(c)所示地图分别实现全覆盖。每个实验数据取10次实验平均值。其中单目标方法把工作时长作为算法的主要求解目标。
如图9所示,显示了两种方法对路径规划时间的比较实验。不难发现本文双目标方法在使用2 台乃至更多机器人时,算法在运算速度方面有明显的加快,在稳定性上也有明显提升。尤其需要注意的是,在机器人台数增加时本文的双目标适应度方法的路径规划时间(算法复杂度)会有所减少。主要原因是本文动态分割染色体后使用了动态规划算法来优化Cell内部覆盖路径的入口点与出口点组合。随着机器人的数目增加,相同地图下动态规划算法的计算规模会有所减少。所以在机器人台数增加时,本文所使用的算法更能发挥其优势。
图9 路径规划时间的对比实验Fig.9 Comparative experiments on timing of path planning
为了检验算法在多机器人全覆盖效果方面的表现,同时也为了对比本文的双目标适应度方法与传统的单目标适应度方法在各种指标上的表现。如表2所示,不难看出两种方法在机器人台数为1时表现相同,主要因为当面对单机器人全覆盖问题时,两种方法可以看作是相同方法。本文的双目标适应度方法随着机器人数目增加,相同环境下的工作时长呈线性减小,说明算法求解对于不同数目机器人全覆盖问题具有较好的稳定性。由于本文提出的算法考虑了每个机器人起始位置和终点位置的调度问题,所以随着机器人数目增加总工作量和路径重复率都呈线性小幅度增加。同时,各种情况实验后不均衡度值均小于1%且并无异常值出现,说明任务分配相对均衡。在面对MCCPP 时,本文所提供的算法不管在工作时长、总工作量以及路径重复率方面都有小幅优化。所以不难看出,本文所提出的用于解决立面维护作业的多机器人全覆盖路径规划的双目标适应度免疫遗传算法具有很好的综合收益。
表2 对比实验Table 2 Comparison experiments
针对目前多机器人全覆盖算法存在任务分配不均、对机器人起始位置和终点位置选择的灵活性考虑不充分、算法求解目标单一以及复杂不可行域地图覆盖灵活性不足等问题。本文将立面维护作业下的多机器人全覆盖问题视为多旅行商问题的扩展,作者将改进的免疫遗传算法结合求解多旅行商问题的综合方法应用于多机器人全覆盖领域。首先,本文完成了一些相关工作:(1)解决了多机器人任意起点和终点全覆盖路径优化问题;(2)提出了覆盖任务动态均衡分配算子;(3)为了加快算法收敛速度使用的具有启发式的优秀染色体片段移植算子;(4)为了寻找全局近似最优解将动态规划算法与阶梯型克隆选择算子相结合;(5)在任务得到均衡分配和考虑立面维护作业机器人起点位置和终点位置的任意性特点的前提下,本文适应度函数同时包含多机器人工作中最大的工作时长最小和多机器人总工作量最小两个求解目标。最后,通过对不同复杂程度仿真实验,验证了算法的可行性;通过对复杂地图的全覆盖仿真实验,验证了该算法对多个优化目标的近似解具有较好的收敛性;通过对比实验,深入分析了本文所提出的算法的复杂度、运算速度以及运算精度,检验了算法的质量。未来,在面对复杂的实际应用场景,均匀Cell的面积、优化算法时间和空间复杂度以及优化感兴趣区域权重因子分配方法等方面可作为下一阶段研究内容。