求解旅行商问题的自适应升温模拟退火算法

陈科胜,鲜思东,郭 鹏

(重庆邮电大学复杂系统智能分析与决策重点实验室,重庆 400065)

1 引言

旅行商问题(traveling salesman problem,TSP)是一个非确定性多项式难度的优化问题,其问题可以描述为给定一系列城市坐标集,一个旅行者从起点城市出发,如何经过各个城市一次并回到出发点的路径规划问题.其问题的解可以被描述为各个城市出现一次且仅出现一次构成的排列方案,该排列方案代表着旅行者将第一个城市作为出发点,依次按排列顺序对城市进行仿问,最后再回到第一个城市的路径规划方案.TSP问题的最优解就是旅行者所经历过的最短路径.

TSP问题可以简单地描述为已知n个城市以及各个城市相互之间的距离,求出某一旅行商经过所有城市并回到出发点的最短路线.设d(Vi,Vj)为Vi到Vj的距离,问题可以抽象为已知一点集V{Vi,1 ≤i≤n},寻找一个排列X{V1,……,Vn}来最小化

多年来研究人员不断地探究该问题,文献[1]梳理了许多针对TSP问题的算法方案.总体上可以将算法分为两个类别:一类是在已有成熟的启发式算法中,结合TSP问题的特点对算法机制进行改良,使得算法具有更强的寻优能力;一类是借助仿生学等思想提出新的启发式算法.

其中，在改良现有启发式算法的方面,文献[2]中提出了离散型萤火虫群优化算法(discrete glowworm swarmoptimization algorithm,DGSO),文献[3]中提出的自适应离散粒子群算法(self-adaptive discrete particle swarm optimization algorithm,SADPSO)及文献[4]中提出的具有变异特征的蚁群算法(ant colony algorithm with mutation features,ACO)、文献[5]基于已有的非支配排序遗传算法–II(non-dominated sorting genetic algorithm–II,NSGA–II),利用平均距离将种群划分为若干个大致均匀分布的小种群保持种群多样性.文献[6]基于混沌的最大最小蚁群算法采用tent混沌映射的方法产生初始信息素当算法陷入长时间停滞状态时利用混沌映射扰动信息素的改良蚁群算法(maxmin ant system algorithm,MMAS)、利用凸包的构建及构建的区域范围与凸包角提出的一种基于动态凸包引导的偏优规划蚁群算法(ant colony algorithm of partially optimal programming based on dynamic convex hull guidance,ACADCG);在改良遗传退火算法方面文献[7]中提出的改进遗传模拟退火(improved genetic simulatedannealing algorithm,IGSA)、文献[8]中提出的改进Metropolis(M)准则的自适应模拟退火算法(improved genetic simulated annealing algorithm,IGSAA)算法.

在提出新的启发式算法方面,文献[9]中提出基于帝国竞争的政治关系思想启发提出帝国竞争算法、文献[10]中基于狼群仿生学思想提出了针对TSP问题的离散狼群算法、文献[11]提出基于蝙蝠算法的观测矩阵优化算法.这些算法相对于朴素的智能算法都能得出更优的效果,但有些基于旧算法的改进算法过于复杂、新提出的仿生算法难以应用到现有模型中.笔者经研究后,总结这些算法的特点,提出了一种能得到更优结果、收敛时间更快的简洁、轻量型的温控模拟退火算法(self-adaptive genetic simulated annealing algorithm),并针对TSP问题进行了优化,在求解TSP问题时耗时更短,全局寻优能力更强.并使用了国际通用的TSP标准库TSPLIB的Burma14,Ulysses16,Oliver30,eil51,eil76,tsp225,St70等TSP测试数据进行了检测.

2 朴素模拟算法

2.1 模拟退火算法简介

模拟退火算法在处理全局优化、离散变量优化等高度非线性化的优化问题中,在处理该类问题中具有传统经典算法无可比拟的优势.模拟退火算法的思想最早由Metropolis等人提出的.其出发点是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性.模拟退火法是一种通用的优化算法,其物理退火过程由以下3部分组成:

1) 加温过程:其目的是增强粒子的热运动,使其偏离平衡位置.当温度足够高时,固体将熔为液体,从而消除系统原先存在的非均匀状态.

2) 等温过程:对于与周围环境交换热量而温度不变的封闭系统,系统状态的自发变化总是朝自由能减少的方向进行,当自由能达到最小时,系统达到平衡状态.

3) 冷却过程:使粒子热运动减弱,系统能量下降,得到晶体结构.其中,加温过程对应算法的设定初温,等温过程对应算法的Metropolis 抽样过程,冷却的过程对应控制参数的下降.这里能量的变化就是目标函数,要得到的最优解就是能量最低态.其中Metropolis准则是模拟退火算法收敛于全局最优解的关键所在,Metropolis准则以一定的概率接受恶化解,这样就使算法跳离局部最优的陷阱.

2.2 朴素模拟退火算法的实现

1) 设置初始温度,随机初始化一个初始解;

2) 开始进入循环:

2a) 进入新解产生流程:

在城市访问排列中,以一定的规则对当前解进行随机的扰动;

记录下最优解;

判断是否达到新解产生的迭代次数,达到则将扰动得到的最优解作为新解.

2b) 按照Metropolis 准则,根据扰动产生的新解Xnew的能量值E(Xnew)与扰动前的解Xold的能量值E(Xold)和当前温度t判断是否接受新的解,如若接受新的解,则用新解代替旧解将Xnew的值赋给Xold.其中关于接受概率的Metropolis准则为

2c) 温度下降.

3) 判断温度是否达到温度阈值,若未达到继续进行外循环.

3 自适应升温模拟退火算法

针对传统模拟退火算法(见图1)在求解问题时容易陷入局部最优解的情况,本文设计了一种自适应的升温控制因子,使得算法能够有针对性地控制局部寻优以及全局寻优能力.其改进部分可分为初始解产生机制、新解产生机制、自适应升温因子的设计、Metropolis函数的优化等4个部分组成.

图1 传统模拟退火算法流程示意图Fig.1 Flow chart of traditional simulated annealing algorithm

3.1 初始解的产生机制

在模拟退火算法中,选取恰当的初始解能够改善模拟退火算法的收敛速度,减小算法的收敛时间.初始解的生成中,本文使用了安德鲁算法(Andrew’s Algorithm)将城市的位置V{Vi,1 ≤i≤n}视为一个点集来生成凸包.在凸包的基础上,使用文献[12]介绍的三角形TSP法.其主要思想为反复随机添加一个城市与已知的城市连接,使得添加该城市后,增加的两条路径减去一条路径的总代价最小,直到该哈密顿回路包括所有城市.如在图2中,为随机选取节点d添加回路,分别比较已有添加方式的总代价,最后选取一个总代价最小的方案.最后将得到的哈密顿回路作为模拟退火算法的初始解.

图2 三角形TSP算法示意图Fig.2 Schematic diagram of triangle TSP algorithm

为了测试该方法的改良效果,采用TSBLIB 中的eil51,eil76等两个TSP问题进行测试,如表1所示.

表1 改良初始解测试结果Table 1 Improved initial solution test result table

从图3–6的结果可以看到,对初始解的机制进行改进之后,该算法的初始解相较于改进前更接近于理论最优值,即全局最优值在在较高温度的阶段保持着初始解的值,但由于环境温度较高,该算法仍有很大概率接受差解,因此该算法仍能保持改进前的寻优能力.在改进后可以适当降低环境初始温度,减少算法的时间开销.

图3 SA求解eil51过程图Fig.3 Solving eil51 process diagram

图4 SA(改进初始解产生机制)求解eil51过程图Fig.4 SA(improved initial solution)for solving eil51

图5 SA求解eil76过程图Fig.5 SA solving eil76 process diagram

图6 SA(改进初始解产生机制)求解eil76过程图Fig.6 SA(improved initial solution)for solving eil76

3.2 新解的产生机制

在模拟退火算法中,新解是对已有的最优解进行一定程度的扰动来产生的.新解产生的机制是决定算法寻优能力的重要因素之一.将新解的产生过程分解为扰动操作、择优操作.扰动操作可视为一个随机过程,负责对解进行搜寻,择优操作则在一定程度上保证了产生的新解的质量.防止自适应扰动因子引入造成的不收敛问题.

扰动操作:对已有旧解执行扰动的具体操作为不断的以一定的概率对旧解的城市进行调换操作,当随机调换操作停止时视为得到了一个新的解.其中随机调换操作指的是,对于当前的解

随机选择其中的节点Vi与Vj的位置进行互换.换言之,即是对当前解进行一次扰动,进行扰动操作的次数服从参数为P的幂律分布,P为该算法的超参数.扰动操作流程如图7所示.

图7 扰动操作流程示意图Fig.7 Disturbance operation flow diagram

择优操作:在新解的产生过程中,对一个旧解进行扰动操作,得到一个扰动后的解,比较扰动后的解与扰动前的解,若扰动后的解比扰动前的解更优,则视为完成了一次择优操作.不停的进行扰动操作,当择优操作次数等于择优上限Mmax,新解产生的过程结束,将最后产生的解作为产生的新解.其中Mmax为引入的一个控制择优过程的参数,为Mmax[k1T],其中k1为择优操作参数,控制了算法的局部寻优能力.择优操作流程如图8所示.

图8 择优操作流程示意图Fig.8 Schematic diagram of preferred operation flow

3.3 自适应升温因子的设计

应用模拟退火算法在求解TSP问题中,容易陷入局部最优解而导致“早熟”的现象.为了防止这一现象的发生,本文引入了一个自适应的升温控制因子,其主要思想是当模拟退火算法出现过早的收敛现象时,增高温度,根据Metropolis准则,使得扰动的随机性增强,避免收敛于局部最优解、增加全局的寻优能力.实际应用在算法中的操作为:当N次降温之后,仍没有获得比算法当前最优解更好的解,则使系统升高温度为到Tnew.其中N与Tnew为自适应扰动控制因子的参数.N越小时,控制因子敏感度越高,反之相反.Tnew越高时,扰动增加幅度越大,全局寻优的能力越强.其中参数为

其中:s为升温次数;k2为控制升温因子参数,引入k2s的目的是保证了该算法的收敛;Told为当前温度;α为升温系数;k1为择优过程的次数参数,控制的是该算法的局部寻优能力;k2为升温因子参数,控制的是该算法的全局寻优能力.

该自适应升温因子使得算法在达到局部最优解时,仍能保持一定限度的搜索能力,同时升温最大次数N的限制使得该算法能够在有限时间内结束,从理论上能够证明该算法是收敛的,且收敛性等同于一般的模拟退火算法.以求解eil51问题为例,求解过程的温度情况与最优解变化如图9所示.

图9 自适应升温过程中温度与解的收敛趋势图Fig.9 Convergence trend diagram of temperature and solution in adaptive heating process

可以看到在自适应模拟退火算法的搜索过程中,当解接近收敛的时候,升温机制能够增强算法的搜索能力,从而在一定程度上有利于跳出局部最优解.在求解eil51问题中,该算法正是借助了这一机制,如图中虚框部分所示,对温度进行提升,增强算法的搜索能力,从而跳出了局部最优解,得到了更优的解,实际上这一解也是该eil51目前的理论最优解.

3.4 Metropolis函数的优化

在模拟退火算法中,决定是否接受新解与旧解的规则是Metropolis准则.因此Metropolis准则影响着模拟退火算法的收敛速度和寻优能力.

对于不同的TSP问题中E(Xnew)−E(Xold)可能出现的值存在着差距很大的现象,导致了p的值过于大或过于小,进而影响了算法的寻优能力与收敛速度.利用最小生成树构造哈密顿回路的方法来对TSP问题的解的大小规模进行估计,从而来对能量差值E(Xnew)−E(Xold)进行归一化,消除不同TSP问题的路径值的不同规模对Metropolis准则的影响.改良后的Metropolis准则为

其中:E(Xnew)−E(Xold)表示新解与旧解能量值间的差值,在TSP问题中能量值为解所对应的路线长度;A表示对该TSP问题解的规模,可视为对TSP问题的解的一个大概估计.其中得到A的方法为先生成一个较优的初始解,将该较优的初始解对应的长度作为A的值.生成较优的初始解的操作为:

1) 将城市的位置X{V1,……,Vn}视为一个点集,先使用了Prim最小生成树算法来生成连接了所有城市的最小生成树.

2) 在获得了连接所有城市的最小生成树的基础上,对树的根节点进行先序遍历,遍历出所有的顶点构成了有序的顶点集L.

3) 最后将根节点添加到顶点集L的末尾,将这个有序的顶点集L按顺序连接起来,就得到了一条哈密顿回路,将得到的哈密顿回路作为模拟退火算法的初始解.可证明该解不超过最优解2倍.

为了测试该方法的改良效果,采用TSBLIB 中的eil51,eil76等两个TSP问题进行测试,测试结果见表2,图示见图10–11;得到改进升温模拟退火算法流程如图12所示.

表2 改良Metropolis测试结果Table 2 Improved Metropolis test result table

图10 SA(M函数优化)求解eil51过程图Fig.10 SA(M function optimization)to solve eil51 process diagram

图11 SA(M函数优化)求解eil76过程图Fig.11 SA(M function optimization)to solve eil76 process diagram

3.5 自适应升温模拟退火算法收敛性证明

将该算法过程抽象为一般的离散优化问题,在具有一定的状态集中,算法能够收敛到某一个状态不再转移.因此该算法的收敛性证明可以总结为如下的证明.

定理1在有限的状态集中,对于任意的初始解,自适应升温模拟退火算法一定可以转移到某个状态上并不再发生转移,即该算法具有全局收敛性.

在证明定理1时,需要使用到由陈小刚等人[13]利用随机过程证明的引理1–2.

引理1若某个Markov链所对应的状态转移矩阵具有以下两个性质,则该矩阵是遍历且具有稳态的.

1) 该矩阵为不可约矩阵,即对应的Markov链为不可约的.

2) 该矩阵的严格上三角阵趋于0,而严格下三角阵保持不变.

图12 改进升温模拟退火算法总流程图Fig.12 General flow chart of improved simulated annealing algorithm

引理2对于一般模拟退火算法的非平稳的m阶状态转移矩阵A(t),若存在唯一的稳态分布:

则有下式成立:

此定理的证明可以通过以下步骤形成:

步骤1首先,假设TSP问题具有有限个状态,设有不可约m阶概率矩阵PPij根据Metropolis准则以概率p从i状态向j状态进行转移的概率为

其中:

∆为E(j),E(i)之间的差值.在传统的模拟算法中,控制温度具有t1≥t2≥……≥tn且

在本文提出的自适应升温模拟退火算法中,控制温度并不具有诸如t1≥t2≥……≥tn的关系,当算法在第N次降温之后,仍没有获得比算法当前最优解更好的解,则使系统升高温度为∆T,其中参数由式(3)–(4)给出.

因此t并不随着n单调递减,且仍然有

成立.利用反证法容易证明,若不存在

则当n →∞时,算法进行了无穷次升温操作使得控制温度不能下降到0,即s →∞并使得N →∞,则对于t1,t2,……,tn中存在着某个单调递减子序列ti,ti+1,……,ti+N表示某次未触发升温条件的控制温度序列,根据算法机制易知有下式成立:

其中∆T为每次下降温度值,则有

与式(10)不成立矛盾,则有

成立.证毕.

至此该自适应升温模拟退火算法的收敛性证明问题等价于一般模拟退火算法收敛性的证明问题.

步骤2假设某个TSP问题中所有的状态按对应的能量值进行排列,则可以得到不可约的转移矩阵A(t)[aij],其中:

转移矩阵A(t)[aij]当

成立具有以下两个性质:

1) 该矩阵为不可约矩阵,即对应的Markov链为不可约的.

2) 严格上三角阵趋于0,而严格下三角阵不变.

则由定理1可以得到转移矩阵对应的Markov链具有遍历性且具有稳态

则可以得到

或

由式(15)可知当温度趋于0时均可以收敛到一个稳定状态上不再发生转移或由式(16)可知该Markov链的状态不断地进行转移可以使得该算法以概率1转移到某个状态上不再发生改变,且该结论与初始解的选择无关,即该算法具有全局收敛性.证毕.

4 实验及结果分析

为检验该算法的实际性能,本文使用了TSPLIB 标准库的Burma14,Oliver30,eil51,St70,eil76,tsp225进行了检验,与普通的启发式算法进行对比,分别比较了与模拟退火算法SGA[7]、传统蚁群算法ACS[7]、粒子群算法ACO[7]的求解效果;与改进的启发式算法对比,分别比较了利用混沌映射扰动信息素的最大最小蚂蚁系统算法MMAS[11]、基于方向信息素协调的蚁群算法AADC[8]、动态凸包引导的偏优规划蚁群ACADCG算法[9]、自适应离散粒子群算法DGSO[3]、具有变异特征的蚁群算法SADPSO[4];与同类型的改进模拟退火算法对比,分别与改进遗退火算法IGSA[7]、IGSAA改进自适应模拟退火算法[8];与最近学者新提出的新型启发式算法进行对比,分别与帝国竞争算法ICA[9]、离散狼群算法DWPA[10]进行了对比.其中在求解各个TSP问题时本文使用的算法参数如表3所示,其中扰动操作概率值p设定为0.5.

表3 实验参数Table 3 Experimental parameters

使用本文算法,按照表3参数对各个TSBLIB标准库中问题进行求解,结果如表4与图13–18所示.

表4 不同文献算法对比结果Table 4 Comparison results of different literature algorithms

图13 Burma14最优路径(30.8785)Fig.13 Burma14 optimal path(30.8785)

图14 Oliver最优路径(423.7406)Fig.14 Oliver general optimal path(423.7406)

图15 eil51最优路径(426)Fig.15 eil51 general optimal path(426)

图16 eil76最优路径(538)Fig.16 eil76 general optimal path(538)

图17 St70最优路径(677.11)Fig.17 St70 general optimal path(677.11)

图18 tsp225最优路径(3927)Fig.18 tsp225 general optimal path(3927)

5 结论

本文提出的算法在求解中规模TSP问题时具有收敛精度高、全局寻优化能力强的特点.部分TSBLIB的TSP问题中均能较快地得到理论最优解,与普通的启发式算法进行对比,比较模拟退火算法SGA[7]、传统蚁群算法ACS[7]、粒子群算法ACO[7]的求解效果;与改进的启发式算法对比,分别比较了利用混沌映射扰动信息素的最大最小蚂蚁系统算法MMAS[11]、基于方向信息素协调的蚁群算法AADC[8]、动态凸包引导的偏优规划蚁群ACADCG 算法[9]、自适应离散粒子群算法DGSO[3]、具有变异特征的蚁群算法SADPSO[4];与同类型的改进模拟退火算法对比,分别与改进遗退火算法IGSA[7]、IGSAA改进自适应模拟退火算法[8];与最近学者新提出的新型启发式算法进行对比,分别与帝国竞争算法ICA[9]、离散狼群算法DWPA[10]进行了对比,其求解效果也存在均取得了更强的寻优效果.由于朴素的模拟退火算法可由可调参数k1,k2来控制对应的局部寻优能力和全局寻优能力,提高了该算法的推广性与鲁棒性,但从另一方面上说,需要调整的参数较多,也一定程度地带来了使用时调参的困难,存在着进一步改良的空间.