1.简介

旅行商问题（TSP）是一个经典的组合最佳化问题，也被归类为NP难问题。多个旅行商问题（MTSP）代表了TSP的一个变体，与TSP相比更加复杂，具有更大的实际意义。MTSP在机器人、交通和网络等各个领域都有广泛的应用。同样，MTSP也是一个NP难问题。MTSP的优化目标是最小化所有旅行商的路径长度之和。为了解决这个问题提出了一种混合离散长鼻浣熊优化算法（HDCOA）。最初，采用扇区划分方法预先分配城市节点给每个销售员。当分配给每个销售员的节点数超过150个时，引入K-Neest Neighbors（KNN）算法，将这些城市节点进一步划分为更小规模的问题。除了用于求解的经典交换、移位和2-opt算子，旨在克服基本离散COA容易陷入局部最优的局限性，本文设计了一种名为破坏和重建（R&D）的算子。该算子通过消除路径中最长的无向边来实现逃离局部最优。本文不仅解决了城市数从51到18,512的平面MTSP问题，还以地球上7342个真实世界城市为MTSP案例进行了仿真实验。本文的规模远远超过已知文献。大量的实验结果表明，所提出的算法在求解TSP和MTSP时表现出优异的竞争力。

2. 基本概念

2.1. 局部最优的定义

局部最优指在特定区域内最优的解状态，确保优化算法能够远离平稳点并收敛至某个稳定点。

2.2. 全局最优的定义

全局最优指针对给定问题获得的最佳可能解，通常通过结合物理洞察和数学特性推导额外约束条件来保证优化算法的收敛特性。

2.3. 旅行商问题(TSP)的定义

TSP是经典的组合离散优化问题，近年来在计算机科学和运筹学领域被广泛研究。目前尚无高效方法能完全解决该问题。哈密顿循环是指一个不重复的节点序列，其首尾节点形成闭环，例如$[1,2,3,1]$和$[3,4,1,2,3]$。该问题的目标是在所有城市节点中找到长度最短的哈密顿循环。

对称TSP可描述为完全无向图$G=(N,A)$，其中节点集{1,2,…,n}\{1,2,\dots,n\}{1,2,…,n}表示城市，每条边$<i,j>\in A$具有非负权重$e_{ij}$。两城市间的欧氏距离$d_{ij}$计算如下：

$$d_{ij}=\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}\text{\hspace{1em}(1)}$$

对称TSP的数学模型为最小化哈密顿循环总距离$f$：

$$\text{最小化: }f=\sum _{i=1}^n\sum _{j=1}^n{d}_{ij}{q}_{ij}\text{\hspace{1em}(2)}$$

约束条件：
∑i=1nqij=1,j∈{1,2,…,n}(3) \sum_{i=1}^n q_{ij} = 1, \quad j \in \{1,2,\dots,n\} \tag{3} i=1∑n​qij​=1,j∈{1,2,…,n}(3)
∑j=1nqij=1,i∈{1,2,…,n}(4) \sum_{j=1}^n q_{ij} = 1, \quad i \in \{1,2,\dots,n\} \tag{4} j=1∑n​qij​=1,i∈{1,2,…,n}(4)
∑i∈S∑j∉Sqij≥2,∀S⊂{1,2,…,n},2≤∣S∣≤n−2(5) \sum_{i \in S} \sum_{j \notin S} q_{ij} \geq 2, \quad \forall S \subset \{1,2,\dots,n\}, 2 \leq |S| \leq n-2 \tag{5} i∈S∑​j∈/S∑​qij​≥2,∀S⊂{1,2,…,n},2≤∣S∣≤n−2(5)
qij∈{0,1},∀i,j∈{1,2,…,n}(6) q_{ij} \in \{0,1\}, \quad \forall i,j \in \{1,2,\dots,n\} \tag{6} qij​∈{0,1},∀i,j∈{1,2,…,n}(6)

其中$q_{ij}=1$表示存在边$<i,j>$，$q_{ij}=0$则表示不存在。

2.4. 多旅行商问题(MTSP)的定义

MTSP是TSP的扩展模型，由$m$个销售员访问$n$个城市，目标是找到总成本最小的路线集合。根据仓库数量可分为单仓库MTSP和多仓库MTSP（如图1所示）。

图1.单仓库MTSP和多仓库MTSP的区别
本研究聚焦多仓库MTSP，其目标是将无向图$G=(N,A)$的节点集$N$划分为$m$个子集{Sk}k=1m\{S_k\}_{k=1}^m{Sk​}k=1m​，并找到每个子集$S_k$的最小成本路径。目标函数为：

$$\text{最小化: }f=\sum _{i=1}^m\left(d_{0,i}+\sum _{j=1}^{n_i-1}{d}_{j,j+1}\right)\text{\hspace{1em}(7)}$$

其中$d_{0,i}$为第$i$个销售员从仓库到首城市的距离，$d_{j,j+1}$为其路线中相邻城市的距离。同时需平衡各销售员的工作量：

$$\text{最小化: }{f}_d=\sum _{i=1}^m\sum _{j=i+1}^m|R_i-R_j|\text{\hspace{1em}(8)}$$

$R_i$表示第$i$个销售员的总行程距离。

2.5. 浣熊优化算法(COA)

COA由Dehghani等人于2022年提出，模拟了自然界中鬣蜥的两种行为：捕食冻结和逃生。冻结阶段分为两个时期：

1. 恐吓期（探索阶段）：鬣蜥群体分成两半，其中一半攀爬树木躲避捕食者。数学表示为：
   $$x_{ij}^{(1)}=x_{i,j}+r\cdot (f{g}_{ij}-I+x_{i,j}),i=1,2,\dots ,N\text{\hspace{1em}(9)}$$
2. 捕食期（开发阶段）：当受惊时，鬣蜥随机落地，另一半群体实施捕捉。其生物行为如图2所示，(a)为恐吓期，(b)为捕食期。

在公式(9)中，$X_i^{P1}$表示当前迭代后第$i$个个体在解空间中的临时位置。其中：

• $j$表示当前解的维度
• $r$为[0,1]间的随机数
• $Iguana$表示当前最优位置
• $I$为1或2的随机整数

鬣蜥坠落与捕食行为建模

当鬣蜥坠落后，其随机位置由公式(10)确定：
$$Iguan{a}_j^G=l{b}_j+r\cdot (u{b}_j-l{b}_j),j=1,2,\dots ,m\text{\hspace{1em}(10)}$$
地面鬣蜥群的移动行为由公式(11)描述：
$$x_{i,j}^{P1}=\{\begin{array}{ll}{x}_{i,j}+r\cdot (Iguan{a}_j^G-I\cdot x_{i,j}),& F_{Iguana}^G<F_i\\ x_{i,j}+r\cdot (I\cdot x_{i,j}-Iguan{a}_j^G),& \text{否则}\end{array}\text{\hspace{1em}(11)}$$
其中：

• $Iguan{a}_j^G$为鬣蜥坠落位置
• $lb$和$ub$分别为解空间下界和上界
• $F_{Iguana}^G$为坠落鬣蜥的适应度值
• $F_i$为个体$X_i$的适应度值

个体更新规则为：
$$X_i=\{\begin{array}{ll}{X}_i^{P1},& F_i^{P1}<F_i\\ X_i,& \text{否则}\end{array}\text{\hspace{1em}(12)}$$

局部搜索策略

当被捕食时，鬣蜥会逃往当前位置附近的安全区域，该行为对应COA的局部搜索能力。搜索空间随迭代动态收缩：
$$l{b}_j^{\text{local}}=\frac{l{b}_j}{t},u{b}_j^{\text{local}}=\frac{u{b}_j}{t},t=1,2,\dots ,T\text{\hspace{1em}(13)}$$
局部位置更新公式为：
$$x_{i,j}^{P2}=x_{i,j}+(1-2r)\cdot \left(l{b}_j^{\text{local}}+r\cdot (u{b}_j^{\text{local}}-l{b}_j^{\text{local}})\right)\text{\hspace{1em}(14)}$$
更新规则与公式(12)类似：
$$X_i=\{\begin{array}{ll}{X}_i^{P2},& F_i^{P2}<F_i\\ X_i,& \text{否则}\end{array}\text{\hspace{1em}(15)}$$

---

3.提出的HDCOA算法

本节提出了一种基于COA和各种算子组合的混合算法，用于求解TSP和多仓库MTSP。如前所述，MTSP涉及几个要考虑的变量，因此我们的算法将它们分为两部分。第一部分采用基于扇区的划分方法，在每个旅行推销员之间均匀分配城市节点，而第二部分侧重于为每个旅行推销员解决“专用TSP”。本节主要介绍各种策略、算子和HDCOA算法的主要结构。

3.1 基于扇区的分区方法

基于扇区的划分方法是一种经典的数学分类方法。在平面坐标系中，它首先确定所有散点的中心位置作为原点o。穿过原点o的任意线，例如本研究中选择的水平线，作为参考线。每个散点通过一条直线连接到原点o，并计算每个直线与参考线之间的角度。然后根据它们的角度大小顺序对点进行分类。如图3所示，这种方法允许我们获得每个旅行推销员的初始特定访问点。TSPLIB（https://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/）是一个广泛使用的库，为TSP和相关优化问题提供基准实例。图4显示了将TSPLIB数据集usa13509中的13,509个城市分类为3个集群的示意图，其中每个集群对应于特定旅行推销员要访问的城市节点。

图4.通过基于扇区的分区方法将usa13509平均分为3类。

3.2.问题编码

顺序编码方法是表示MTSP解决方案的最简单和最有效的方法，利用一系列数字来表示每个解决方案。有两种常见的编码方法：路径序列和断点序列[43]，以及带有城市索引序列的路径序列[44]。与上述编码方法相比，本文采用二维数组来存储每个旅行推销员访问的城市序列，如图5所示。每个一维向量代表特定旅行推销员的哈密顿循环。

图5.多仓库MTSP的编码方法。

3.3. K-最近邻（KNN）算法

近年来，研究人员利用k-means聚类算法求解TSP。然而，随着TSP规模的不断扩大，这种方法的实用性减弱。理论上成熟且最简单的机器学习算法之一的k-近邻（KNN）分类算法被用来解决这个问题。这种方法的基本思想是，在需求空间中，如果一个样本的k个最近邻中的大多数属于某一类，那么该样本也属于该类。很少有研究人员应用KNN算法求解TSP。本文采用KNN算法对具有150多个城市节点的TSP实例进行分类，其中每个类通过将城市分配给各自的类来单独求解。图6说明了使用KNN分类方法求解的TSPLIB中rl5915的示例。

图6. KNN算法分类示例。

3.4.重新设计方程式

在求解TSP的过程中，当前的局部最优解总是与全局最优解的一部分相似。图7演示了eil51的近似最优解和TSPLIB中的已知最优解。绿色路径表示当前最优解和已知最优解之间的相似部分。通过利用这一特性，我们将以前方法中的方程（9）（简称COA）重新设计成方程（16）。

公式说明

位置更新公式

$$X_i^{P1}:x_i^{P1}=2\cdot Q\oplus Iguana,i=1,2,3,\dots ,\lfloor \frac{N}{2}\rfloor \text{\hspace{1em}(16)}$$

动态衰减系数

$$Q=1-e^{-\frac{t}{T}},t=1,2,3,\dots ,T\text{\hspace{1em}(17)}$$

其中：

• $Q$为随迭代衰减的系数（$T$为最大迭代次数）
• $\oplus$表示序列拼接操作
• $Iguana$表示当前最优解序列
• $\lfloor \frac{N}{2}\rfloor$为参与更新的个体数量下限
• 

图7. 近似最优解与已知最优解路径的比较

在本文中，每个可行解表示为一组汉明循环序列。用于连续问题的解公式可能不适用于解决离散问题。在公式(16)中，$l_{guanq}$表示当前全局最优解，$x_i$表示第 i 个个体的解，Q是一个介于0和1之间的小数，随着迭代的进行单调递减，如公式(17)所示。这里，( t ) 表示当前迭代次数，T 表示最大迭代次数。$Q\otimes l_{guanq}$ 表示从 $l_{guanq}$中选择一部分序列。例如，如果$Q=0.4$且 $l_{guanq}$ = {4, 3, 2, 1, 4}，则从 $l_{guanq}$ 中随机选择40%的部分，例如$3,2$，结果为 Q⊗lguanq={3,2}Q \otimes l_{guanq} = \{3, 2\}Q⊗lguanq​={3,2}。符号 $\oplus $表示将选定部分连接到$x_i$的序列中。例如，如果$x_i = {2, 1, 3, 4, 2}$，则 xi⊕{3,2}={1,3,2,4,1}x_i \oplus \{3, 2\} = \{1, 3, 2, 4, 1\}xi​⊕{3,2}={1,3,2,4,1}。图8提供了此计算的示例。

3.5. 交换算子和移位算子

交换算子和移位算子是解决TSP的经典方法。在本文中，这两种方法用于替代狩猎过程(公式(11))和逃避捕食过程(公式(14))。在COA中，交换算子的目的是随机交换一组城市序列中两个城市的位置，而移位算子随机选择城市序列的一个子序列，并将其插入原始序列中的随机位置。这些算子增强了HDCOA的全局搜索能力。图9提供了交换算子的示例，其中 swap(2,6) 表示序列中第二和第六个城市的交换。图10演示了移位算子的示例，其中 shift(2,4,7) 表示将第七个城市之后的第二到第四个城市的子序列移动到序列中。值得注意的是，这些数字是随机生成的。

3.6. 2-opt算子

在解决TSP问题时，2-opt算子是一种必不可少的方法。其主要功能是消除城市序列段内的交叉路径。图11提供了2-opt算子的示例。它首先确定边 ( e_{i,i+1} ) 是否与边 ( e_{j,j+1} ) 相交（其中 ( i > j )）。如果发生交叉，则反转从第 ( (i+1) ) 个城市到第 ( j ) 个城市的城市序列段，有效消除交叉路径。

3.7. 破坏和重建 (D&R) 算子

在解决TSP问题的过程中，问题的规模往往导致更高的可能性陷入局部最优。传统策略在有效解决此问题方面不足。为了解决这个问题，本文提出了一种称为破坏和重建(D&R)的新算子。该算子主要采用最差路径来消除最长的边，并在重新连接它们之前重建分割的部分，( d ) 由公式(18)给出。

$$d=\frac{\text{city\_number}}{k}\text{\hspace{1em}(18)}$$

function [new_pop,fbest] = TSP_DCOA(N,T,pop,distance_matrix,city_coord)
city_num = length(pop);
N1 = round(N/2);
k = 150;
fitness = inf(1,N);
x = zeros(N,city_num+1);
% Initialize the population
for i = 1:N
    if i <= 2
        x(i,:) = KNN2root_MTSP(pop,k,city_coord);
    else
        idx = randperm(city_num);
        x(i,1:city_num) = pop(idx);
        x(i,end) = x(i,1);
    end
    x(i,:) = tsp_has_intersect_path(city_coord, x(i,:));
    fitness(i) = f(distance_matrix,x(i,:),city_num);
end
[~,best_index] = min(fitness);
flag = 0;
for t = 1:T
    Q = 1 - exp((t-T)/T);
    
    for i = 1:N1
        iguana = x(best_index,:);
        
        I = randi([0,1]);
        x_p1(i,:) =  prey(x(i,:),iguana,I,rand,Q);
        f_p1(i) = f(distance_matrix,x_p1(i,:),city_num);
        
        if f_p1(i) < fitness(i)
            x(i,:) = x_p1(i,:);
            fitness(i) = f_p1(i);
        end
    end
    
    for i = N1+1:N
        x_p1(i,:) = shift2(x(i,:),distance_matrix,city_num+1);
        f_p1(i) = f(distance_matrix,x_p1(i,:),city_num);
        
        if f_p1(i) < fitness(i)
            x(i,:) = x_p1(i,:);
            fitness(i) = f_p1(i);
        end
    end
    
    [~,best_index] = min(fitness);
    
    for i = 1:N
        x_p2 = swap2(x(i,:),distance_matrix,city_num+1);
        f_p2 = f(distance_matrix,x_p2,city_num);
        if f_p2 < fitness(i)
            x(i,:) = x_p2;
            fitness(i) = f_p2;
        end
    end
    
    % Determine whether it has fallen into a local optimum
    if t > 2 && history(t-2) == history(t-1)
        flag = flag+1;
    else
        flag = 0;
    end
    if flag == 50
        x(best_index,:) = tsp_has_intersect_path(city_coord, x(best_index,:)); % 2-opt
        x(best_index,:) = inflection_point_exchange(city_coord, x(best_index,:));
        flag = 0;
        [~,worest] = max(fitness);
        x(worest,:) = div(x(best_index,:),distance_matrix,city_coord,city_num,k); % D&R
        fitness(worest) = f(distance_matrix,x(worest,:),city_num);
    end
    
    [~,best_index] = min(fitness);
    fitness(best_index) = f(distance_matrix,x(best_index,:),city_num);
    
    fbest = fitness(best_index);
    history(t) = fbest;
    best_root = x(best_index,:);
end
new_pop = best_root;
end

Zhu L, Zhou Y Q, Luo Q F. Hybrid discrete coati optimization algorithm for solving large-scale multiple traveling salesman problem[J]. Ain Shams Engineering Journal, 2025, 16(10): 103637. DOI: 10.1016/j.asej.2025.103637.