许秋艳，马 良，刘 勇

(1. 上海理工大学管理学院，上海 200093；2. 盐城工学院信息工程学院，江苏 盐城224051)

1 引言

最小比率旅行商问题(minimum ratio traveling salesman problem，MRTSP)是经典旅行商问题的一种扩展形式：除考虑行程距离外，推销员从一个城市到另一个城市将获得一定的收益，该问题的目标就是选择最佳旅行线路，最小化总路程和总收益的比值。该模型同时考虑距离和收益，和只考虑距离最小相比，该目标函数等价于最小化费用效益比，更符合实际情况。

旅行商问题属于NP难问题，在其基础上构建的最小比率旅行商问题MRTSP也属于NP难问题[1]。目前，用于求解MRTSP的算法主要有模拟退火算法、竞争决策算法、大洪水算法、蝙蝠算法、蚁群优化算法、引力搜索算法和中心引力优化算法等[2，3]。这些方法为求解MRTSP提供了思路，但是计算效果需要进一步优化。为有效处理MRTSP问题，本文将基于阴阳平衡优化(Yin-Yang-pair optimization，YYPO)算法[4]设计求解方法。

阴阳平衡优化算法是在2016年由Punnathanam等提出的，其优化策略源于中国传统文化——阴阳学说。在智能计算领域，算法的局部开发与全局探索能力既对立又互补，如何实现两者的平衡是算法设计的关键。YYPO算法将局部开发与全局探索分别视为阴和阳两个方面，期望通过阴阳平衡实现局部开发与全局探索的均衡[5]。

因独特的优化机制，阴阳平衡优化算法备受关注。目前，该算法已经成功用于发动机系统设计[6]、光伏逆变器PID调节[7]、风力涡轮机设计[8]和连续型选址[9]等。总体来看，阴阳平衡优化算法的研究处于起步阶段，主要用于求解连续优化问题，还没有涉及到组合优化问题。本文将阴阳平衡优化算法用于求解最小比率旅行商问题，扩大算法的应用范围，也为该问题的求解提供可行有效的方法。为此，本文采用佳点集法产生初始群体，提高初始解质量；基于超球体和归档集对解不断进行更新；并采用相对位置索引法直接产生可行解，避免多个约束条件的处理；基于综卦变换引入反转操作，提高算法的局部搜索性能。数值实验验证了算法求解最小比率旅行商问题的可行性和有效性。

2 最小比率旅行商问题MRTSP的数学模型

最小比率旅行商问题的具体模型可描述如下[1-3]

本文考虑的网络属于完全图G=(V，E) (否则，可通过计算任两个点间最短路转换成等价的完全图)。令V={1，2，…，n} 代表顶点集，表示城市的位置；E={(i，j)：i≠j，i，j∈V}代表边的集合；D=(dij)n×n代表距离矩阵，表示任两个城市间的距离，其中∀i，j∈V，dij=dji，dii=+∞，dij>0；P=(pij)n×n代表收入矩阵，表示旅行商(推销员)完成从一个城市到另一个城市的销售可得到的收入，其中∀i，j∈V，pij=pji，pii=0，pij>0。设那么最小比率旅行商问题的模型为

(1)

其中，式(1)表示目标函数，在经过所有点的回路中求总行程与总收入之比最小；式(2)和(3)要求在一次行程中每个城市只访问一次；式(4)要求没有任何子回路生成；|S|代表子图S中包含的图G中顶点个数；式(5)表示决策变量取值范围。同时满足上述约束条件的解就形成了一条Hamilton 回路。

3 最小比率旅行商问题的阴阳平衡优化算法

受到太极图的启示，阴阳平衡优化算法利用半径为1的超球体表示其搜索空间。算法首先在超球体内随机生成两个点，并比较这两个点的适应度值。令P1表示值较优的点，而令P2表示值较差的点。然后分别以P1和P2为中心，δ1和δ2为步长在超球体中进行优化搜索。在寻优过程中，δ1会逐步减小而δ2会逐步变大。因此，点P1的搜索区域将逐渐收缩而进行小范围寻优，在算法中被称为局部开发，有收缩的含义，和阴对应；点P2的搜索区域将逐渐扩展而进行大范围寻优，在算法中被称为全局探索，有扩张的含义，和阳对应。利用P1和P2进行的局部开发和全局探索缺一不可，有阴阳互根的含义。点P1的搜索范围越来越小，而点P2的搜索范围越来越大，有阴阳对立的含义。在算法寻优中，不断增强基于点P1的局部搜索和点P2的全局搜索能力，有阴阳皆长的含义。算法期望通过基于点P1和点P2的两种搜索以实现局部开发和全局探索的权衡，有阴阳平衡的含义。以下给出阴阳平衡优化算法求解最小比率旅行商问题的具体实现过程。

3.1 基于佳点集的初始化

考虑到最小比率旅行商问题NP-难的特征，设置合理的初始解，可改善算法的搜索效率。此外，虽然智能优化算法不依赖于初始解，但是提高初始解质量能够加快优化速度。而阴阳平衡优化算法仅通过P1和P2两个点进行优化搜索，提高P1和P2对应初始点质量，可提高算法优化效率。本文不再采用基本阴阳平衡优化算法对P1和P2随机初始化的方法，而采用一种性能较好的初始化方法——佳点集法。在佳点集中选择最好的两个解赋值给P1和P2。

佳点集的概念由华罗庚等首先提出，具体方法如下[10]：

设Gn为n维空间的单位立方体，r=(r1，r2，…，rn) 为Gn内的点。令Pm(k)是Gn内的包含m个点的集合，如果形为Pm(k)={({r1k}，{r2k}，…，{rsk}，k=1，2，…，m}的偏差φ(m)满足：φ(m)=C(r，ε)m-1+ε，那么就称r为佳点，Pm(k)为佳点集。其中，ε代表无穷小；C(r，ε)代表仅和r与ε有关的常量。

在具体应用时，可以设置rk={2cos 2πk/q}，1≤k≤n，q为符合(q-3)/2≥n的最小素数，那么r即为佳点。此外，也可以设置rk={ek}，1≤k≤n那么r也是佳点。其中，{t}代表t的小数部分。研究表明，相对随机点集而言，佳点集分布更加均匀[2，11]。

3.2 基于超球体的解更新

算法分别采用P1和P2作中心点，δ1和δ2作步长，并在超球体中进行优化搜索。因为P1和P2两点采用相同的更新方法，为方便起见，用P表示P1或P2，并用δ表示δ1或δ2。首先将P进行复制，生成2D(D表示变量维数)个同样的点，并分别用NP1，NP2，…，NP2D进行表示；然后对这些点进行更新。在算法中，可以对点的一个分量或者所有分量进行更新。对点的一个分量进行更新时，具体方法如下：

(6)

(7)

在对点的所有分量进行更新时，具体方法如下：

(8)

(9)

在算法中，通过等概率来选择上述两种方法的一种进行点的更新。此外，在利用点P1生成的2D个新解中，挑选适应度值最优的点，直接替换P1。对P2生成的新解也进行同样的操作。

3.3 基于归档集的解更新

在计算过程中，算法采取精英保留策略，在进行基于超球体的解更新前，会将P1与P2两个点存储到归档集内。每隔I代将利用归档集中的2I个点对当前的P1与P2两点做更新操作。其中，I表示归档集的更新频率，并且I是在Imin和Imax之间随机产生的整数。

该阶段的解更新方法如下：第一步，从归档集内挑选适应度值最优的点并和P1进行比较，如果该点优于P1，则这两个点进行交换；第二步，再次从归档集内中挑选适应度值最优的点并和P2进行比较，如果该点优于P2，则这两个点也进行交换；第三步，清空归档集，并生成的新的归档集更新频率参数I。

此外，算法在该阶段还对两个搜索步长δ1和δ2进行更新，具体方法如下：

δ1=δ1-(δ1/α)

(10)

δ2=δ2+(δ2/α)

(11)

其中，α表示收缩/扩张因子；δ1和δ2分别表示是点P1和P2的搜索步长。

3.4 候选解的生成

为充分发挥阴阳平衡优化算法求解连续优化问题的优势，仍然采用连续空间表示算法的搜索空间。但最小比率旅行商问题属于组合优化问题，可行解表示城市的访问顺序。如何实现算法搜索个体的位置到问题解的转换是必须要解决的问题。考虑到最小比率旅行商问题的解是离散的特征，本文采用基于相对位置索引法的解生成方法[12]。

在相对位置索引法中，首先将最小实数转换为最小整数1，然后将第二小的实数转换为下一个高位整数2，以此类推，可以对所有实数进行编号。最后利用编号构成城市的访问顺序。该方法和随机键编码方法类似，以下通过一个算例作进一步说明。例如，一个有四个城市的最小比率旅行商问题，阴阳平衡优化算法的一个搜索个体位置为(0.52，0.68，0.15，0.97)，按照上述方法，得到的城市的访问顺序为(2，3，1，4)。

此外，如果在转换前向量中有多个相同的分量，这种方法可能会产生不可行解。本文设计一种修复机制，即根据向量位置的先后顺序进行编码。以下以相对位置索引法举例说明，转换前的向量(0.02，0.58，0.78，0.98，0.78)，第三个分量和第五个分量取值相同，根据上述修复机制，其编号分别为3和4，转换后的向量为(1，2，3，5，4)。

因此，采用上述候选解的生成方法，可以直接产生可行解，每个城市只被访问一次，并且在任何一个子集中不会构成完整的Hamilton 回路。在算法搜索过程中，可以直接计算目标函数值，不再考虑模型中四个约束条件。

3.5 基于综卦变换的局部搜索

基本阴阳平衡优化算法的全局搜索能力较强，而局部搜索能力较弱[13，14]。为求解最小比率旅行商问题，需要提高算法的局部寻优能力。和阴阳学说一样，《易经》也属于中国传统文化。从算法设计角度来看，卦的变化蕴含了优化思想。例如，综卦是把本卦倒过来看，例如图1所示，姤卦的综卦是夬卦。综卦的哲理是从另一个角度来看同一问题[15，16]。

图1 综卦变化示意图

受此启发，将当前解视为本卦，对其进行综卦变化，具体将采用反转操作。在反转操作中，旋转变化的城市位置不再是起点和终点，而是随机选择两个城市位置作为起点和终点，然后将随机选出的两个城市之间的访问顺序颠倒过来。图2给出一个示例，其中T表示反转操作前的解，T′ 表示反转操作后的解。

图2 反转操作示意图

3.6 算法主要计算步骤

综上所述，求解最小比率旅行商问题的阴阳平衡优化算法具体步骤如下：

步骤1：根据佳点集法产生初始点，并选择最好的两个点分别赋值给P1和P2

步骤2：在基于超球体的解更新阶段，按等概率原则分别利用P1与P2进行优化搜索；

步骤3：每间隔I次迭代，则进入基于归档集的解更新阶段，否则，转步骤4；

步骤4：若迭代次数达到预设最大值，转步骤6；否则，转步骤5。

步骤5：若点P2适应度值优于点P1适应度值，则交换这两个点的取值与搜索步长，转步骤2；

步骤6：对当前最好解进行基于综卦的反转操作，算法停止，输出最好结果。

4 数值实验

实验共分为三个部分：首先采用并行程序设计阴阳平衡优化算法；然后比较阴阳平衡优化算法和微粒群优化算法、引力搜索算法、生物地理学优化算法以及最有价值球员算法的求解性能；最后对距离矩阵和收益矩阵的元素进行参数灵敏度分析。在本文中，采用最小比率旅行商问题的3个典型问题进行仿真。这些算例定义如下[1-3]：

算例1： 设给定对称赋权完全图G={V，E}，距离矩阵D和收益矩阵P定义如下

算例2： 设给定对称赋权完全图G={V，E}，距离矩阵D和收益矩阵P定义如下

算例3： 设给定对称赋权完全图G={V，E}，距离矩阵D和收益矩阵P定义如下：

4.1 算法并行程序设计

迄今为止，在智能优化算法领域，使用并行程序设计算法是一个重要方向，现有工作主要是将计算任务分配给多台计算机完成，而在一台计算机上实现并行智能优化算法的研究工作极其罕见。当前计算机至少配置了双核或者四核的CPU，在体系结构上具有实现并行计算的硬件条件。此外，还有许多支持并行计算的软件环境，例如，Windows多线程编程、MPI并行编程和Matlab并行计算工具箱等[17]。

为充分利用计算资源，设计者需要认识到底层多核的存在，采用并行程序设计智能优化算法，提高算法性能。在阴阳平衡优化算法中，所有个体基于超球体的解更新操作具有相互独立性，具有天然的并行特征。此外，和其它智能优化算法一样，阴阳平衡优化算法采用群体搜索策略，具有隐含并行性，例如各个体适应度函数值的评价可以并行化等。在本文中，采用Matlab的并行计算工具箱(Parallel Computing Toolbox)对阴阳平衡优化算法进行并行程序设计。

以算例1为例进行实验，对比算法并行和串行的计算时间。在Windows10操作系统下，CPU：Intel(R) Xeon(R) E-2186M、32G内存、2.90GHz主频的计算机上运行，调用六核进行并行计算。算法参数设置为：Imin=5，Imax=10，α=25，a=2。此外，δ1与δ2初始值都设为0.5；为防止步长δ2无限制的增大，超出算法的搜索区域，将其最大值设为0.75；最大迭代次数T=500。

和大多数智能优化方法一样，阴阳平衡优化算法属于随机方法，需要多次运行同一算法，统计其优化性能。在实验中，循环次数分别为1，100，200，300，400，500，比较并行和串行计算的耗时，结果如图3所示。

图3 算法并行和串行计算时间对比

此外，还计算了这5次的加速比(串行算法的执行时间除以并行算法的执行时间)，最大加速比为4.5231。采用并行程序设计算法，能够充分利用空闲的CPU内核，能够减少运行时间。因此，本文所有算法均采用并行程序设计，以提高算法运行速度。

4.2 算法优化性能测试

为测试阴阳平衡优化算法YYPO性能，采用微粒群优化算法(particle swarm optimization，PSO)[18]、引力搜索算法(gravitational search algorithm，GSA)[19]、生物地理学优化算法(biogeography-based optimization，BBO)[20]和最有价值球员算法(most valuable player algorithm，MVPA)[21]进行比较。

为保证实验公平性，这5种算法设置相同的函数评价次数，即2000n，其中n为问题规模。YYPO算法其它参数设置保持不变，其它算法根据文献[18-21]进行参数设置。为避免算法单次运行结果的偶然性对分析优化性能的影响，每种算法各自独立运行30次，分别统计最优值、最差值、平均值和标准差等指标。此外，统计结果还考虑每个算法在30次实验中达到已知最优解的成功次数，并进行排序。具体计算结果如表1所示。

表1 不同算法计算结果比较

图4 算法优化过程对比

通过分析表1中的结果不难发现：对于最小规模的算例1，在相同的函数评价次数下，这5种算法都能收敛到已知最优解。但是对更大规模的算例2和3，虽然微粒群优化算法PSO、引力搜索算法GSA、生物地理学优化算法BBO和最有价值球员算法MVPA都可以找到已知最优值，但是在最差值、平均值、标准差、成功次数和排序这五个指标上，这些算法都明显劣于阴阳平衡优化算法YYPO。此外，阴阳平衡优化算法YYPO所获得的标准差明显优于其它4种算法，说明本文算法的优化性能更加稳定。上述结果说明阴阳平衡优化算法YYPO在计算精度方面明显优于其它4种优化算法。

为进一步分析算法优化性能，图4给出了这5种算法优化速度对比示意图。从图中可以发现，阴阳平衡优化算法YYPO不仅具有更高的计算精度，而且具有更快的收敛速度。该算法为最小比率旅行商问题的求解提供一个可行有效的算法。

4.3 灵敏度分析

在最小比率旅行商问题中，两个城市间的距离和旅行商的收入直接影响总行程和总收益之比。这里，分析距离矩阵D和收益矩阵P发生变化时对结果的影响。为方便起见，考虑距离矩阵和收益矩阵中元素同时增加或同时减少的情况。当这两个矩阵元素取值反向变化时，总行程和总收益之比容易判断变化趋势。例如，当距离矩阵元素取值都变大而收益矩阵元素取值都变小时，总行程和总收益之比变大。

这里，以算例1为例，分析这两个矩阵元素取值同向变化时，总行程和总收益之比变化情况。原总行程为199，总收益为234，总行程和总收益之比为0.8504。

当距离矩阵元素值都增加5，收益矩阵元素都增加1时，总行程变为224，总收益为239，两者比值为0.9372。此时，目标函数值变差。当距离矩阵元素值都增加10，收益矩阵元素都增加20时，总行程变为249，总收益为334，两者比值为0.7455。此时，目标函数值变好。

当距离矩阵元素值都减少6，收益矩阵元素都减少2时，总行程变为169，总收益为219，两者比值为0.7717。此时，目标函数值变好。当距离矩阵元素值都减少8，收益矩阵元素都减少16时，总行程变为159，总收益为154，两者比值为1.0325。此时，目标函数值变差。

距离矩阵元素值都增加相当于旅行商通勤距离变大，服务对象扩展到郊区。此时，旅行商的收益也会增加，但不表示费用效益比也会增加，即总行程和总收益之比可能变大也可能变小。与此类似，距离矩阵元素值都减小相当于交通条件的改善，通勤距离变短。此时，旅行商的收益会减少，但不表示费用效益比也会减少，即总行程和总收益之比可能变小也可能变大。

最小比率旅行商问题的目标是考虑总行程和总收益之比，类似于经济学中费用效益比，与仅关注总距离相比，既考虑总行程，又考虑旅行商的收益，更具有现实意义。在满足约束的情况下，如何兼顾成本和旅行商收益需要根据具体问题进行定量计算，以实现最佳的调度方案。

5 结论

为有效求解最小比率旅行商问题，构建了阴阳平衡优化算法。首先，引入佳点集法来初始化搜索群体，改善算法初始解的性能；其次，采用基于超球体和归档集的两种策略对不同阶段的解进行更新，并根据相对位置索引法进行解的编码来产生可行解，以满足模型中所有的约束条件；最后，根据综卦变换采用反转操作，避免算法早熟收敛。此外，算法采用并行程序设计，以充分利用多核处理器计算资源。与其它四种算法相比，在计算精度和优化速度两方面，所设计的算法均排名第一。通过距离矩阵和收益矩阵的灵敏度分析发现，配送距离和收益同时增加或者减小时，需要根据具体问题进行定量计算，才可以确定最小比率(费用效益比)的变化情况。将算法用于车辆路径问题是进一步的研究工作。