OM | 机器学习算法求解更大规模的旅行商问题

fwalker

作者信息： 张云天，王飞龙，孙楚天，杨李平，曲晨辉
本篇责任编辑 张云天，王飞龙
转载请联系 @运筹OR帷幄 后台或同名公众号

编者按： 针对旅行商问题 (Traveling Salesman Problem, TSP) ，现有的基于监督学习的求解算法因泛化能力弱而受限。为了克服这一局限性，本文试图以监督学习的方式训练一个小规模的模型，并基于一系列的技术，如图采样、图转换和图合并，为任意规模的TSP构建热力图。进一步地，热力图被输入强化学习算法，具体为蒙特卡洛树搜索，以指导高质量解决方案的搜索。基于大量TSP算例（最多10,000个节点）的实验结果表明，这种新方法明显优于现有的基于机器学习的TSP算法，并显著提高了训练模型的泛化能力。
引言

旅行商问题（Traveling Salesman Problem, TSP）是著名的组合优化问题，属于NP-hard问题。学者们已设计诸多精确求解器和启发式算法来求解TSP，例如Concorde，LKH等。但是，它们都非常复杂，由许多设计巧妙的规则组成，并严重依赖于专家知识。为了克服这些局限性，近年来基于机器学习的算法已经被尝试应用于求解TSP。



Figure 1. 旅行商问题（Traveling Salesman Problem, TSP）

机器学习方法不过度依赖专家知识，可以很容易地推广到各种组合优化问题上，成为一个备受关注的研究方向。目前的求解TSP的机器学习算法主要有两个类别：监督学习和强化学习。

• 监督学习（Supervised Learning, SL）算法：试图从预先计算好的大量TSP算例中提取一般的模式（Common pattern）。
  
• 强化学习（Reinforcement Learning, RL）算法：试图在与环境的交互过程中学习求解（没有预先计算的TSP算例）。
  

一旦完成训练，SL模型就能提供有用的模式信息，提高TSP高质量解的搜索效率。然而，在处理不同规模的TSP算例时，因为训练算例的分布与测试算例的分布相差较大，预训练模型的性能往往会急剧下降。另一方面，训练SL模型通常需要提前求解大量最优（至少是高质量的）TSP算例，这对于大规模的TSP算例来说计算代价过大。这些缺点严重限制了SL在大规模TSP上的应用。
既然我们认为，以监督学习的方式寻找一般的模式（Common pattern）是可行的，那么如果能在合理的时间内训练出一个小规模的SL模型，并设计一种方法将其顺利推广到大规模的案例中（无需再次预先计算大量的算例），就有希望利用SL的优点，同时避免其缺点。在这个想法的激励下，本文开发了一系列技术，以提高由SL训练的模型的泛化能力。进一步地，将SL和RL结合起来，形成一种混合算法，这种算法与现有的基于机器学习的TSP求解算法相比表现良好。总的来说，主要的贡献总结如下：

• 算法方面：
  
  
  
  • 首先，本文设计一个带有注意力机制的图卷积残差网络（Att-GCRN），基于监督学习训练了一个小规模（节点数为m）的模型。一旦训练完成后，给定一个有m个节点的TSP算例，该模型能够根据图的边建立一个热力图。热力图刻画了节点之间每条边属于该TSP最优环游的概率。
    
  • 然后，本文尝试将这个模型泛化到更大规模的问题。给定一个大规模TSP算例，基于预训练的小模型，与本文设计的三种图技巧：图采样，图转换和图合并，可以为大规模TSP算例构建热力图。
    
  • 最后，基于得到的热力图，引入基于强化学习的算法，具体为蒙特卡洛树搜索（Monte Carlo tree search, MCTS），来进一步搜索高质量解。
    
  
  
• 结果方面：在有多达10,000个城市的TSP算例上进行了仿真实验（比目前用于评估机器学习类方法的算例大一个数量级）。在所有的数据集上，本文提出的算法都能在合理的时间内获得最优或接近最优的解决方案，并明显优于所有现有的基于机器学习的算法。
  

算法整体流程

给定一个任意规模的TSP算例，算法整体流程如Figure 2所示。整体来说，本文提出的是一种混合算法，级联了SL和RL模块。SL模块运用预训练的模型与三种图技巧，为任意规模的TSP构建热力图。热力图刻画了节点之间每条边属于该TSP最优环游的概率。随后，热力图将作为RL模块的重要输入。RL算法将进一步搜索寻优，输出TSP的最优环游。



Figure 2. 本文所设计算法的整体流程

监督学习部分
本文通过预训练一个小模型，来为m个节点的TSP问题构建热力图。但是，当问题规模急剧扩大，其节点数远大于m时，预训练的模型便会失效。为了使得模型可以泛化到更大规模的场景，本文有针对性地设计了三种技巧：图采样，图转换，图合并。具体如Figure 3所示。



Figure 3. 本文设计的三种图技巧与预训练模型的使用

• 图采样：图采样技巧从输入的图G中采样数个子图，其中每个子图都有m个节点。对于每一个节点i\in V，记O_i为节点i在被采样子图中出现的次数，初始化为0。进而，算法每一次迭代过程中，我们选择O_i值最小的节点i（如果存在多个这样的节点，随机选择一个）作为聚类中心，并引入[2]中的k-最近邻算法提取一个包含m个节点j的子图G'。对于G'中的每一个节点i，O_i\leftarrow O_i+1。同时，记O_{ij}为边(i,j)在被采样子图中出现的次数，初始化为0。上述过程不断重复，直至O_i中最小的值达到设定的阈值\omega。
  
• 图转换：图转换旨在修正采样子图中节点的分布。训练集中，每个算例中的所有节点，随机分布在一个单位方格内。为确保每一个采样子图G'符合训练集算例的分布，需要将其转换为新子图G''。转换方式如下：x^{min}=\min\limits_{i\in G'}x_i，x^{max}=\max\limits_{i\in G'}x_i，y^{min}=\min\limits_{i\in G'}y_i，y^{max}=\max\limits_{i\in G'}y_i分别表示子图G'中m个节点的水平坐标最小值、最大值，垂直坐标最小值、最大值。此外，令s=\frac{1}{max(x^{max}-x^{min},y^{max}-y^{min})}为放大因子。进而可以将子图G'中的坐标(x_i,y_i)转换为新子图G''中的坐标(x_i^{new},y_i^{new})
  

x_i^{new}\leftarrow s\times (x_i-x^{min})\\ y_i^{new}\leftarrow s\times (y_i-y^{min})\\

• 图合并：经过图转换的子图，并行输入预训练好的模型，即Att-GCRN模型，并行输出G''的热力图（heatmap）。现在，需要合并数个子热力图。对于原图G的每一条边(i,j)，计算其属于TSP最优环游的概率P_{ij}如下。其中，P_{ij}^{''}(l)表示第l个子图G''中，边(i,j)属于G''最优环游的概率。
  

P_{ij}=\frac{{1}}{O_{ij}}\times \sum\limits_{l=1}^{I}P_{ij}^{''}(l) \\
通过使用预训练的Att-GCRN模型，与三种图技巧，可以得到输入的原图G的完整热力图。补充，所有对应值P_{ij}\leq 10^{-4}的边都被视作不好的边，在后续算法搜索寻优过程中直接忽略，以压缩搜索空间。 作者注解：类比思考，图采样、图转换和图合并的设计，很类似大数据领域的MapReduce。为了处理海量数据，MapReduce先将原任务划分为多个子任务，并在多台计算性能相对有限的机器上并行处理，最终将结果聚合输出。而本文的图采样和图转换，对应了子任务的划分与预处理；各子图随后被相对较小的预训练模型Att-GCRN并行处理；图合并最终聚合并行计算出的各子热力图。如此，本文基于规模有限的预训练模型，为大规模的TSP构建了热力图，正如MapReduce可以用多台计算性能相对有限的机器处理海量数据一样。



Figure 4. 作者注解：MapReduce设计示意图

其中，Att-GCRN，即带有注意力机制的图卷积残差网络，架构如Figure 5所示。在设计上，Att-GCRN引入了图嵌入（Graph embeeding）、残差连接、注意力机制等组件，设计细节可参见本文的附录部分。损失函数（Loss Function）参考了focal loss，设计为
FL(\boldsymbol{P})=-\frac{1}{M}\sum\limits_{m=1}^M\sum\limits_{i=1}^N\sum\limits_{j\neq i}^N(\beta_0(1-\boldsymbol{I}{ij}^{(m)})+\\ \beta_1\times\boldsymbol{I}{ij}^{(m)})(\boldsymbol{I}{ij}^{(m)}-\boldsymbol{P}{ij}^{(m)})log(\boldsymbol{P}{ij}^{(m)}).\\
其中，\boldsymbol{I}^{(m)}表示第m个训练算例的一个n\times n矩阵；若\boldsymbol{I}{ij}^{(m)}=1，则表示边(i,j)属于第m个算例最优环游，否则\boldsymbol{I}_{ij}^{(m)}=0；\boldsymbol{P}^{(m)}表示所构建的第m个算例的热力图；\beta_0和\beta_1表示两个平衡因子，可通过下式计算：
\beta_0=\frac{n^2}{n^2-2n}, \ \beta_1=\frac{n^2}{2n} \\



Figure 5. 预训练的Att-GCRN模型架构

强化学习部分

以之前监督学习生成的热力图为输入，本文引入蒙特卡洛树搜索进行高质量解的寻优。强化学习的MDP模型构建，模型的状态，动作和状态转移方法如下所示。
1. 状态和动作

不同于其他使用强化学习求解路径问题的文献，本文中的状态均为一个完整的TSP路径，而一个动作即为从一个TSP解转到另外一个TSP解的方案。本文使用了比较巧妙的动作构建方式，\pmb{a}=(a_1,b_1,a_2,b_2,...,a_k,b_k,a_{k+1})，其中k是一个变量，用于衡量动作序列的长度；同时要求a_{k+1}=a_1。对于此动作的解读如下：a_i,b_i,i \in {1,...,k} 表示TSP路径中的节点，RL需要做的决策就是选出其中的一系列a_i节点，当a_i节点确定完之后，相应在当前TSP路径中与a_i节点相连的b_i节点也可以直接确定出来，无需做额外的选择。当动作\pmb{a}确定之后，需要做的就是将弧(a_i,b_i),i \in {1,...,k}进行删除，之后添加弧(b_i,a_{i+1}),i\in{1,...,k}。这样设计动作的好处首先可以将动作选择的2k-sub decision决策空间缩小为k sub decision决策空间，因为只需要确定出a_i节点的选择问题就足够；
同时，这样设计可以保证得到的新的状态的可行性。 同时，若根据热力图得到的P_{b_i,a_{i+1}}<10^{-4}，则将此种动作视为不可行，直接删除之。从构造算法的角度来看，这种动作设计方式其实类似于邻域搜索算法中的k-opt算子，即删除掉k条边，之后增加k条边。



Figure 6. 一个具体的3-opt案例讲解

2. 奖励计算

执行完一个动作之后的奖励计算方法如下所示：
\Delta(\pi, \pi^*)=\sum_{i=1}^{k}d_{b_ia_{i+1}}-\sum_{i=1}^{k}d_{a_ib_{i}} \\
若\Delta(\pi, \pi^*)<0则说明新的状态\pi^*要优于\pi。
3. 状态初始化

本文采用随机加贪婪的方式进行初始化，贪婪的衡量标准需要参考热力图中每条弧的概率分布。首先随机选择一个客户点\pi_1作为当前状态的第一个客户点，之后后续客户点\pi_i, i \in {2,...,n}的添加根据exp(P_{\pi_ij})的大小进行选择。
4. 2-opt邻域的全枚举

由于本文的动作类似于k-opt优化，所以在进行邻域搜索时，可以采用小邻域，也可以采用大邻域。对于2-opt邻域搜索，本文使用充分探索小邻域的思想，将当前TSP解方案所对应的所有的2-opt小邻域对应的解全部枚举出来，之后在这些邻域解中选出最优的解方案，从而使得解方案快速收敛到一个局部最优解。
5. 针对大邻域搜索的蒙特卡洛方法

针对大邻域的搜索，不能使用全枚举的方式进行探索，因为时间效率太低。所以本文设计了一种基于蒙特卡洛树搜索的强化学习方法来对更大的邻域进行探索，其流程如Figure 7所示。蒙特卡洛树搜索的流程如下所示：



Figure 7. 蒙特卡洛树搜索方法流程

• 初始化：上述MCTS中涉及的需要训练的参数主要包括以下三个：权重系数矩阵W，W_{ij}表示从节点i选择节点j的概率；到达次数矩阵Q,Q{ij}表示弧(i,j)在simulation过程中被选择的次数；变量M表示在simulation过程中生成的动作数量。其中W,Q为对称矩阵。
  
• 模拟：给定一个状态\pi，也即一个TSP路径方案，simulation的工作是生成一系列的动作，从而尝试改进状态\pi。首先也是从\pi中随机选择一个客户点a_i，之后在确定出b_i之后，需要选择下一个节点a_{i+1}，此时需要引入一个概率计算方式：Z_{b_i,j},j \in {1,...,n} \&\& j \neq a_i \ and \ j \neq nodes \ already\ connected\ to\ b_i，Z_{b_i,j}的计算方式如下所示：
  

Z_{b_i,j}=\frac{W_{b_ij}}{\Omega_{b_i}}+\alpha\sqrt{\frac{ln{(M+1)}}{Q_{b_ij}+1}} \\
其中，\Omega_{b_i}=\frac{\sum_{j\neq b_i}W_{b_ij}}{\sum_{j\neq b_i}1}表示和节点b_i相关联的所有边的平均被选择的概率，上述公式中，\frac{W_{b_ij}}{\Omega_{b_i}}项倾向于选择被选择概率高的边进行添加到动作序列中-intensification；而\alpha\sqrt{\frac{ln{(M+1)}}{Q_{b_ij}+1}}倾向于选择被选择次数少的边加入到动作序列中-diversification；参数\alpha用于调整intensification和diversification的占比，也即均衡“探索”与“开发”，从而避免陷入局部最优。计算完Z_{b_i,j}之后，选择下一个节点的概率计算如下所示：
P_j=\frac{Z_{b_i,j}}{\sum_{l \in X}Z_{b_il}} \\
其中，X表示和b_i相连的W_{b_ij}\geq1的节点集合。当某个动作得到了更好的解方案或者i\geq10，令a_{i+1}=a_i，即强制结束动作simulation的过程。所有产生的动作均被保存在一个动作sampling pool中，直到产生一个能得到更优解方案的动作，或者sampling pool中的动作数量达到阈值H。

• 选择：在上述Simulation过程中，若得到一个能产生更优解方案的动作，则将得到的更优的状态替换当前状态\pi，若sampling pool达到阈值仍没有更优解方案产生，则使用Initialization生成新的初始状态，再次进行MCTS。
  
• 回溯：网络通过误差反向传播来进行参数的学习，当生成了一个新动作之后，M的值累加1；对于每个出现在动作中的弧(b_i,a_{i+1})，令Q_{b_ia_{i+1}}累加1；对于生成新的最优解方案\pi^{new}的动作\pmb{a}=(a_1,b_1,a_2,b_2,...,a_k,b_k,a_{k+1})\\，对于动作序列中每一条弧(b_i, a_{i+1}),i \in {1,...,k},W_{b_ia_{i+1}}的更新方式如下所示：
  

W_{b_ia_{i+1}}\longleftarrow W_{b_ia_{i+1}}+\beta[exp(\frac{L(\pi)-L(\pi^{new})}{L(\pi)})-1] \\
注意到只有在找到新的更优的解方案的时候，才需要对参数W_{b_ia_{i+1}}进行更新。

• 终止条件：MCTS对依照上述流程不断进行寻优，直到规定的寻优时间T达到，则终止寻优过程，返回上述寻优中找到的最优解方案。
  

实验结果

为评估算法的性能，本文在大量的TSP算例上进行了实验，并与八个基于机器学习的算法以及三个非学习的算法进行了比较。为确保公平，所有的学习类算法统一在一台GTX 1080 Ti GPU上执行。对于非学习算法，因源代码不支持在GPU上运行，所以本文使用Intel(R) Xeno(R) Gold 5118 CPU @ 2.30Ghz (8 cores) 运行算法。
数据集

本文采用两个数据集：(1) Set 1，分为三个子集，每个包含10000个自动生成的TSP算例，分别为n=20,50,100。该数据集被现有的基于学习的算法广泛使用。(2) Set 2, 按照同样的规则，生成400个规模较大的算例，即n=200,500,1000各包含128个算例，和n=10000的16个算例。
参数

如前所述，本文的算法依赖于六个超参数(m, ω, α, β, H, T)。对于控制预训练模型大小的m, 在 Set1 中设定为m=20，在 Set 2 中设定为m=50。对于以下四个参数， 本文默认设置为ω = 5;α = 1;β = 10;H = 10n。最后，对于控制终止时间的参数T, 分别在 Set 1 和 Set 2 中设置为T=10n和T=40n。
Set 1 实验结果

Table 1 列出了本文算法（AttGCRN+MCTS）在 Set 1 上获得的结果，并于其他算法进行比较。为了确保公平的比较，对于一些基于学习的算法，本文调整了其原始参数以延长总运行时间。被调整的算法运行结果以蓝色表示。对于AttGCRN+MCTS，时间被分为两部分，即建立热力图（heat maps）的时间和运行MCTS的时间。
如 Table 1 所示，三种非学习算法在所有测试算例上都获得了良好的结果，而现有的基于学习的算法在n=100的算例上都难以达到最优。与之相比，AttGCRN+MCTS表现良好，在平均差距（Gap）和总运行时间（Time）上具有竞争力。



Table 1. 各算法分别在10,000个算例(Set 1)上的运行结果（n=20, 50, 100）

Set 2 实验结果

Table 2 列出了各算法分别在128个n=200,500,1000算例上的运行结果。Concorde和LKH3在这些算例上仍然表现良好，而Gurobi在n=1000的算例上未能在合理时间内终止。对于基于学习的算法，其结果远不及最优方案，特别是在n=1000的算例上。而AttGCRN+MCTS能在短时间内获得非常接近最优的结果，明显优于现有的学习算法。


此外，本文评估了AttGCN+MCTS在n=10000的16个最大算例上的表现。在这些算例上，[3]中提出的三种基于学习的算法都因为内存异常而失败；而两个精确求解器（Concorde和Gurobi）以[4]中的两个GAT模型都因为时间异常而失败（每个实例最多允许5小时）。
因此，本文排除了这七种基线算法，只将Att-GCRN+MCTS算法与其余三种基于学习的算法[5]进行比较。
如 Table 3 所示，Att-GCRN+MCTS能够产生接近LKH3的解决方案，平均差距为4.3902%。相比之下，三种基于学习的算法对应的平均差距很大。并且Att-GCRN+MCTS的运行时间仍然是合理的（比三个基线中最好的一个短）。



Table 3. 各算法在16个算例(Set 2)上的运行结果（n=10,000)

关于热力图的消融研究（Ablation study）

为了强调热力图的重要性，对于每个算例，本文给每条边分配一个相等的概率，并单独重新运行MCTS算法来搜索解决方案。结果总结在 Table 4 中，左边部分列出了原始Att-GCRN+MCTS算法得到的结果，右边部分列出了单独MCTS得到的结果（没有热力图）。显然，在禁用热力图后，算法的性能急剧下降，即每个数据集上与最优结果产生巨大的差距（Opt. Gap）。作为比较，原始的Att-GCRN+MCTS算法在每个数据集上产生的结果都非常接近最佳状态。该实验证明了热力图，即识别有潜力的侯选边，对算法性能的价值。



Table 4. 热力图的消融实验

相关工作

如Figure 8所示，现有的基于机器学习的TSP解决方法大致可以分为两类。一种是传统方法（非学习方法），一类是基于机器学习的方法。其中，基于机器学习的方法又可进一步划分为监督学习方法和基于强化学习的方法。基于监督学习的方法中三个典型的方法是Pointer Network、GNN (Graph Neural Network) 和GCN (Graph Convolutional Networks)。基于强化学习的方法主要可区分为MCTS （Monte Carlo Tree Search)、Actor-critic以及GAN (Graph Attention Network) 方法。



Figure 8. 本文相关工作框图

总结

基于监督学习的技术对于发现一般的模式非常有用，但需要大量的训练数据，难以推广到大规模的TSP。 本研究表明，通过应用图采样、图转换和图合并等一系列技术，有可能以监督学习的方式训练一个小规模的模型，并顺利地将其泛化应用到大规模TSP的求解中。实验结果证实，该方法能够开发出具有高度竞争力的基于机器学习的TSP算法，并显著提高预训练模型的泛化能力。在未来，我们将尝试解决更大规模的TSP或非欧几里得TSP，并将该方法扩展到其他具有挑战性的优化问题。
参考文献

[1] Fu, Z. H., Qiu, K. B., & Zha, H. (2021, May). Generalize a small pre-trained model to arbitrarily large tsp instances. In Proceedings of the AAAI Conference on Artificial Intelligence (Vol. 35, No. 8, pp. 7474-7482).
[2] Dudani, S. A. (1976). The distance-weighted k-nearest-neighbor rule. IEEE Transactions on Systems, Man, and Cybernetics, (4), 325-327.
[3] Joshi, C. K., Laurent, T., & Bresson, X. (2019). An efficient graph convolutional network technique for the travelling salesman problem. arXiv preprint arXiv:1906.01227.
[4] Deudon, M., Cournut, P., Lacoste, A., Adulyasak, Y., & Rousseau, L. M. (2018, June). Learning heuristics for the tsp by policy gradient. In International conference on the integration of constraint programming, artificial intelligence, and operations research (pp. 170-181). Springer, Cham.
[5] Wouter Kool, Herke van Hoof, and Max Welling. (2019). Attention, learn to solve routing problems! International Conference on Learning Representations (ICLR).

<hr/>欢迎关注（公众号） @运筹OR帷幄 ，了解更多运筹学、人工智能及相关学科的干货资讯。
知乎专栏： 『运筹OR帷幄』大数据时代的运筹学 『运筹AI帷幄』大数据时代的人工智能
获得最新运筹学及其相关学科的干货资料、行业前沿信息、学术动态、硕博offer信息等。特别适合你~
欢迎大家交流。
也欢迎大佬们投稿和商业合作，学术信息、会议通讯、征稿启事、硕博招生信息等免费推广。