CVRP基于matlab蚁群算法求解带容量的车辆路径规划问题含Matlab源码 1039期

Posted 2021-12-28 紫极神光

一、VRP简介

1 VRP基本原理
车辆路径规划问题(Vehicle Routing Problem，VRP)是运筹学里重要的研究问题之一。VRP关注有一个供货商与K个销售点的路径规划的情况，可以简述为：对一系列发货点和收货点，组织调用一定的车辆，安排适当的行车路线，使车辆有序地通过它们，在满足指定的约束条件下（例如：货物的需求量与发货量，交发货时间，车辆容量限制，行驶里程限制，行驶时间限制等），力争实现一定的目标（如车辆空驶总里程最短，运输总费用最低，车辆按一定时间到达，使用的车辆数最小等）。
VRP的图例如下所示：

2 问题属性与常见问题
车辆路径问题的特性比较复杂，总的来说包含四个方面的属性：
（1）地址特性包括：车场数目、需求类型、作业要求。
（2）车辆特性包括：车辆数量、载重量约束、可运载品种约束、运行路线约束、工作时间约束。
（3）问题的其他特性。
（4）目标函数可能是总成本极小化，或者极小化最大作业成本，或者最大化准时作业。

3 常见问题有以下几类：
（1）旅行商问题
（2）带容量约束的车辆路线问题(CVRP)



该模型很难拓展到VRP的其他场景,并且不知道具体车辆的执行路径，因此对其模型继续改进。



（3）带时间窗的车辆路线问题
由于VRP问题的持续发展，考虑需求点对于车辆到达的时间有所要求之下，在车辆途程问题之中加入时窗的限制，便成为带时间窗车辆路径问题（VRP with Time Windows, VRPTW）。带时间窗车辆路径问题（VRPTW）是在VRP上加上了客户的被访问的时间窗约束。在VRPTW问题中，除了行驶成本之外, 成本函数还要包括由于早到某个客户而引起的等待时间和客户需要的服务时间。在VRPTW中，车辆除了要满足VRP问题的限制之外，还必须要满足需求点的时窗限制，而需求点的时窗限制可以分为两种，一种是硬时窗（Hard Time Window），硬时窗要求车辆必须要在时窗内到达，早到必须等待，而迟到则拒收；另一种是软时窗（Soft Time Window），不一定要在时窗内到达，但是在时窗之外到达必须要处罚，以处罚替代等待与拒收是软时窗与硬时窗最大的不同。


模型2(参考2017 A generalized formulation for vehicle routing problems)：
该模型为2维决策变量



（4）收集和分发问题
（5）多车场车辆路线问题
参考(2005 lim，多车场车辆路径问题的遗传算法_邹彤, 1996 renaud)

由于车辆是同质的，这里的建模在变量中没有加入车辆的维度。


（6）优先约束车辆路线问题
（7）相容性约束车辆路线问题
（8）随机需求车辆路线问题

4 解决方案
（1）数学解析法
（2）人机交互法
（3）先分组再排路线法
（4）先排路线再分组法
（5）节省或插入法
（6）改善或交换法
（7）数学规划近似法
（8）启发式算法

5 VRP与VRPTW对比

二、蚁群算法简介

1 引言
在自然界中各种生物群体显现出来的智能近几十年来得到了学者们的广泛关注，学者们通过对简单生物体的群体行为进行模拟，进而提出了群智能算法。其中， 模拟蚁群觅食过程的蚁群优化算法(Ant
Colony Optimization， A CO) 和模拟鸟群运动方式的粒子群算法(ParticleS warm Optimization，PSO) 是两种最主要的群智能算法。
蚁群算法是一种源于大自然生物世界的新的仿生进化算法，由意大利学者M.Dorigo， V.Mani ezzo和A.Color ni等人于20世纪90年代初期通过模拟自然界中蚂蚁集体寻径行为而提出的一种基于种群的启发式随机搜索算法[1].蚂蚁有能力在没有任何提示的情形下找到从巢穴到食物源的最短路径，并且能随环境的变化，适应性地搜索新的路径，产生新的选择。其根本原因是蚂蚁在寻找食物时，能在其走过的路径上释放一种特殊的分泌物――信息素2，随着时间的推移该物质会逐渐挥发，后来的蚂蚁选择该路径的概率与当时这条路径上信息素的强度成正比。当一条路径上通过的蚂蚁越来越时，其留下的信息素也越来越多，后来蚂蚁选择该路径的概率也就越高，从而更增加了该路径上的信息素强度。而强度大的信息素会吸引更多的蚂蚁，从而形成一种正反馈机制。通过这种正反馈机制，蚂蚁最终可以发现最短路径。
最早的蚁群算法是蚂蚁系统(Ant System， AS) ， 研究者们根据不同的改进策略对蚂蚁系统进行改进并开发了不同版本的蚁群算法，并成功地应用于优化领域。用该方法求解旅行商(TSP) 问题、分配问
题、车间作业调度(job-shop) 问题， 取得了较好的试验结果[3-6] 。蚁群算法具有分布式计算、无中心控制和分布式个体之间间接通信等特征，易于与其他优化算法相结合，它通过简单个体之间的协作表现出了求解复杂问题的能力，已被广泛应用于求解优化问题。蚁群算法相对而言易于实现，且算法中并不涉及复杂的数学操作，其处理过程对计算机的软硬件要求也不高，因此对它的研究在理论和实践中都具有重要的意义。
目前，国内外的许多研究者和研究机构都开展了对蚁群算法理论和应用的研究，蚁群算法已成为国际计算智能领域关注的热点课题。虽然目前蚁群算法没有形成严格的理论基础，但其作为一种新兴的进
化算法已在智能优化等领域表现出了强大的生命力。

2 蚁群算法理论
蚁群算法是对自然界蚂蚁的寻径方式进行模拟而得出的一种仿生算法。蚂蚁在运动过程中，能够在它所经过的路径上留下信息素进行信息传递，而且蚂蚁在运动过程中能够感知这种物质，并以此来指导
自己的运动方向。因此，由大量蚂蚁组成的蚁群的集体行为便表现出一种信息正反馈现象：某一路径上走过的蚂蚁越多，则后来者选择该路径的概率就越大[7]。
2.1真实蚁群的觅食过程
为了说明蚁群算法的原理，先简要介绍一下蚂蚁搜寻食物的具体过程。在自然界中，蚁群在寻找食物时，它们总能找到一条从食物到巢穴之间的最优路径。这是因为蚂蚁在寻找路径时会在路径上释放出
一种特殊的信息素。蚁群算法的信息交互主要是通过信息素来完成的。蚂蚁在运动过程中，能够感知这种物质的存在和强度。初始阶段，环境中没有信息素的遗留，蚂蚁寻找事物完全是随机选择路径，
随后寻找该事物源的过程中就会受到先前蚂蚁所残留的信息素的影响，其表现为蚂蚁在选择路径时趋向于选择信息素浓度高的路径。同时，信息素是一种挥发性化学物，会随着时间的推移而慢慢地消逝。如果每只蚂蚁在单位距离留下的信息素相同，那对于较短路径上残留的信息素浓度就相对较高，这被后来的蚂蚁选择的概率就大，从而导致这条短路径上走的蚂蚁就越多。而经过的蚂蚁越多，该路径上残留的信息素就将更多，这样使得整个蚂蚁的集体行为构成了信息素的正反馈过程，最终整个蚁群会找出最优路径。
若蚂蚁从A点出发，速度相同，食物在D点，则它可能随机选择路线ABD或A CD。假设初始时每条路线分配一只蚂蚁， 每个时间单位行走一步。图5.1所示为经过8个时间单位时的情形：走路线ABD的蚂蚁到达终点：而走路线A CD的蚂蚁刚好走到C点， 为一半路程。

图5.2表示从开始算起， 经过16个时间单位时的情形：走路线ABD的蚂蚁到达终点后得到食物又返回了起点A，而走路线A CD的蚂蚁刚好走到D点.

图5.2蚂蚁出发后经过16个时间单位时的情形
假设蚂蚁每经过一处所留下的信息素为1个单位，则经过32个时间单位后，所有开始一起出发的蚂蚁都经过不同路径从D点取得了食物。此时ABD的路线往返了2趟， 每一处的信息素为4个单位； 而A CD的路线往返了一趟，每一处的信息素为2个单位，其比值为2：1。
寻找食物的过程继续进行， 则按信息素的指导， 蚁群在ABD路线上增派一只蚂蚁(共2只) ， 而A CD路线上仍然为一只蚂蚁。再经过32个时间单位后，两条线路上的信息素单位积累为12和4，比值为3：1。若按以上规则继续， 蚁群在ABD路线上再增派一只蚂蚁(共3只) ， 而A CD路线上仍然为一只蚂蚁。再经过32个时间单位后， 两条线路上的信息素单位积累为24和6，比值为4：1.若继续进行， 则按信息素的指导， 最终所有的蚂蚁都会放弃A CD路线， 而选择ABD路线。这也就是前面所提到的正反馈效应。

2.2人工蚁群的优化过程
基于以上真实蚁群寻找食物时的最优路径选择问题，可以构造人工蚁群， 来解决最优化问题， 如TSP问题。人工蚁群中把具有简单功能的工作单元看作蚂蚁。二者的相似之处在于都是优先选择信息素浓度大的路径。较短路径的信息素浓度高，所以能够最终被所有蚂蚁选择，也就是最终的优化结果。两者的区别在于人工蚁群有一定的记忆能力，能够记忆已经访问过的节点。同时，人工蚁群再选择下一条路径的时候是按一定算法规律有意识地寻找最短路径，而不是盲目的。例如在TSP问题中， 可以预先知道当前城市到下一个目的地的距离。在TSP问题的人工蚁群算法中， 假设m只蚂蚁在图的相邻节点间移动，从而协作异步地得到问题的解。每只蚂蚁的一步转移概率由图中的每条边上的两类参数决定：一是信息素值，也称信息素痕迹：二是可见度，即先验值。
信息素的更新方式有两种：一是挥发，也就是所有路径上的信息素以一定的比率减少，模拟自然蚁群的信息素随时间挥发的过程；二是增强，给评价值“好”(有蚂蚁走过)的边增加信息素。蚂蚁向下一个目标的运动是通过一个随机原则来实现的，也就是运用当前所在节点存储的信息，计算出下一步可达节点的概率，并按此概率实现一步移动，如此往复，越来越接近最优解。蚂蚁在寻找过程中，或在找到一个解后，会评估该解或解的一部分的优化程度，并把评价信息保存在相关连接的信息素中。

2.3真实蚂蚁与人工蚂蚁的异同
蚁群算法是一种基于群体的、用于求解复杂优化问题的通用搜索技术。与真实蚂蚁通过外信息的留存、跟随行为进行间接通信相似，蚁群算法中一群简单的人工蚂蚁通过信息素进行间接通信，并利用该信息和与问题相关的启发式信息逐步构造问题的解。
人工蚂蚁具有双重特性：一方面，它们是真实蚂蚁的抽象，具有真实蚂蚁的特性：另一方面，它们还有一些真实蚂蚁没有的特性，这些新的特性使人工蚂蚁在解决实际优化问题时，具有更好地搜索较优解的能力。
人工蚂蚁与真实蚂蚁的相同点为：
(1)都是一群相互协作的个体。与真实蚁群一样，蚁群算法由一群人工蚂蚁组成，人工蚂蚁之间通过同步/异步协作来寻找问题的最优解。虽然单只人工蚂蚁可以构造出问题的解，但只有当多只人工蚂蚁
通过相互协作，才能发现问题的最优(次优)解。人工蚂蚁个体间通过写/读问题的状态变量来进行协作。
(2)都使用信息素的迹和蒸发机制。如真实蚂蚁一样，人工蚂蚁通过改变所访问过的问题的数字状态信息来进行间接的协作。在蚁群算法中，信息素是人工蚂蚁之间进行交流的唯一途径。这种通信方式在群体知识的利用上起到了至关重要的作用。另外，蚁群算法还用到了蒸发机制，这一点对应于真实蚂蚁中信息素的蒸发现象。蒸发机制使蚁群逐渐忘记过去的历史，使后来的蚂蚁在搜索中较少受到过去较差解的影响，从而更好地指导蚂蚁的搜索方向。
(3)搜索最短路径与局部移动。人工蚂蚁和真实蚂蚁具有相同的任务，即以局部移动的方式构造出从原点(蚁巢)到目的点(食物源)之间的最短路径。
(4)随机状态转移策略。人工蚂蚁和真实蚂蚁都按照概率决策规则从一种状态转移到另一种相邻状态。其中的概率决策规则是与问题相关的信息和局部环境信息的函数。在状态转移过程中，人工蚂蚁和
真实蚂蚁都只用到了局部信息，没有使用前瞻策略来预见将来的状态。
人工蚂蚁和真实蚂蚁的不同点为：
(1)人工蚂蚁生活在离散的时间，从一种离散状态到另一种离散状态。
(2)人工蚂蚁具有内部状态，即人工蚂蚁具有一定的记忆能力，能记住自己走过的地方。
(3)人工蚂蚁释放信息素的数量是其生成解的质量的函数。
(4)人工蚂蚁更新信息素的时机依赖于特定的问题。例如，大多数人工蚂蚁仅仅在蚂蚁找到一个解之后才更新路径上的信息素。

2.4蚁群算法的特点
蚁群算法是通过对生物特征的模拟得到的一种优化算法，它本身具有很多优点：
(1)蚁群算法是一种本质上的并行算法。每只蚂蚁搜索的过程彼此独立，仅通过信息激素进行通信。所以蚁群算法可以看作一个分布式的多智能体系统，它在问题空间的多点同时开始独立的解搜索，不仅增加了算法的可靠性，也使得算法具有较强的全局搜索能力。
(2)蚁群算法是一种自组织的算法。所谓自组织，就是组织力或组织指令来自于系统的内部，以区别于其他组织。如果系统在获得空间、时间或者功能结构的过程中，没有外界的特定干预，就可以说系统是自组织的。简单地说，自组织就是系统从无序到有序的变化过程。
(3)蚁群算法具有较强的鲁棒性。相对于其他算法，蚁群算法对初始路线的要求不高，即蚁群算法的求解结果不依赖于初始路线的选择，而且在搜索过程中不需要进行人工的调整。此外，蚁群算法的参
数较少，设置简单，因而该算法易于应用到组合优化问题的求解。
(4)蚁群算法是一种正反馈算法。从真实蚂蚁的觅食过程中不难看出，蚂蚁能够最终找到最优路径，直接依赖于其在路径上信息素的堆积，而信息素的堆积是一个正反馈的过程。正反馈是蚁群算法的重
要特征，它使得算法进化过程得以进行。

3 基本蚁群算法及其流程






4 改进的蚁群算法
针对基本蚁群算法一般需要较长的搜索时间和容易出现停滞现象等不足，很多学者在此基础上提出改进算法，提高了算法的性能和效率。
4.1精英蚂蚁系统

4.2 最大最小蚂蚁系统
为了克服基本蚁群系统中可能出现的停滞现象， Thomas Stutz le等人提出了最大-最小(MAX-MIN) 蚁群系统[10] ， 主要有三方面的不同：
(1)与蚁群系统相似，为了充分利用循环最优解和目前为止找出的最优解，在每次循环之后，只有一只蚂蚁进行信息素更新。这只蚂蚁可能是找出当前循环中最优解的蚂蚁(迭代最优的蚂蚁)，也可能是找出从实验开始以来最优解的蚂蚁(全局最优的蚂蚁)；而在蚂蚁系统中，对所有蚂蚁走过的路径都进行信息素更新。
(2) 为避免搜索的停滞， 在每个解元素(TSP中是每条边) 上的信息素轨迹量的值域范围被限制在[min，Tmax] 区间内； 而在蚂蚁系统中信息素轨迹量不被限制，使得一些路径上的轨迹量远高于其他边，从而蚂蚁都沿着同条路径移动，阻止了进一步搜索更优解的行为。
(3)为使蚂蚁在算法的初始阶段能够更多地搜索新的解决方案，将信息素初始化为tmax；而在蚂蚁系统中没有这样的设置。

4.3基于排序的蚁群算法
基于排序的蚁群算法(Rank-BasedAntSystem)是Bull n heimer、Hartl和Strauss等人提出的[11] 。在该算法中， 每个蚂蚁释放的信息素按照它们不同的等级进行挥发，另外类似于精英蚁群算法，精英蚂蚁在每次循环中释放更多的信息素。在修改信息素路径前，蚂蚁按照它们的旅行长度进行排名(短的靠前)，蚂蚁释放信息素的量要和蚂蚁的排名相乘。在每次循环中，只有排名前w-1位的蚂蚁和精英蚂蚁才允许在路径上释放信息素。己知的最优路径给以最强的反馈，和系数w相乘；而排名第r位的蚂蚁则乘以系数“w-(≥0)。信息素如下式所示：

4.4自适应蚁群算法
基本蚁群系统让信息量最大的路径对每次路径的选择和信息量的更新起主要作用，但由于强化了最优信息反馈，就可能导致“早熟”停滞现象。而最大最小蚁群算法将各个路径上的信息量的更新限定在固定的范围内，这虽然在一定程度上避免了“早熟”停滞现象，但在解分布较分散时会导致收敛速度变慢。以上方法的共同缺点在于：它们都按一种固定不变的模式去更新信息量和确定每次路径的选
择概率。
为了克服以上算法的不足， L.M.Gambardella和M.Dorigo提出了基于调节信息素挥发度的自适应蚁群算法[12]。相对基本蚁群算法的改进如下：
(1)在每次循环结束后求出最优解，并将其保留。
(2)自适应地改变p值。当问题规模比较大时，由于信息量的挥发系数p的存在，使那些从未被搜索到的信息量会减小到接近于0，降低了算法的全局搜索能力；当p过大且解的信息量增大时，以前搜索过的解被选择的可能性过大，也会影响到算法的全局搜索能力；通过减小p虽然可以提高算法的全局搜索能力，但又会使算法的收敛速度降低。因此可以自适应地改变p的值。p的初始值p(to)=1；当算法求得的最优值在N次循环内没有明显改进时，p减为：

5 关键参数说明
在蚁群算法中，不仅信息素和启发函数乘积以及蚂蚁之间的合作行为会严重影响到算法的收敛性，蚁群算法的参数也是影响其求解性能和效率的关键因素。信息素启发式因子α、期望启发因子β、信息
素蒸发系数p、信息素强度Q、蚂蚁数目m等都是非常重要的参数，其选取方法和选取原则直接影响到蚁群算法的全局收敛性和求解效率。
信息素启发式因子a
信息素启发式因子a代表信息量对是否选择当前路径的影响程度，即反映蚂蚁在运动过程中所积累的信息量在指导蚁群搜索中的相对重要程度。a的大小反映了蚁群在路径搜索中随机性因素作用的强度，其值越大，蚂蚁在选择以前走过的路径的可能性就越大，搜索的随机性就会减弱；而当启发式因子a的值过小时，则易使蚁群的搜索过早陷于局部最优。根据经验，信息素启发式因子a取值范围一般为[1，4]时，蚁群算法的综合求解性能较好。
期望启发因子β
期望启发因子β表示在搜索时路径上的信息素在指导蚂蚁选择路径时的向导性，它的大小反映了蚁群在搜索最优路径的过程中的先验性和确定性因素的作用强度。期望启发因子β的值越大，蚂蚁在某个局部点上选择局部最短路径的可能性就越大，虽然这个时候算法的收敛速度得以加快，但蚁群搜索最优路径的随机性减弱，而此时搜索易于陷入局部最优解。根据经验，期望启发因子β取值范围一般为[3，
5]，此时蚁群算法的综合求解性能较好。实际上，信息素启发式因子α和期望启发因子β是一对关联性很强的参数：蚁群算法的全局寻优性能，首先要求蚁群的搜索过程必须要有很强的随机性；而蚁群算法的快速收敛性能，又要求蚁群的搜索过程必须要有较高的确定性。因此，两者对蚁群算法性能的影响和作用是相互配合、密切相关的，算法要获得最优解，就必须在这二者之间选取一个平衡点，只有正确选定它们之间的搭配关系，才能避免在搜索过程中出现过早停滞或陷入局部最优等情况的发生。
信息素蒸发系数p
蚁群算法中的人工蚂蚁是具有记忆功能的，随着时间的推移，以前留下的信息素将会逐渐消逝，蚁群算法与其他各种仿生进化算法一样，也存在着收敛速度慢、容易陷入局部最优解等缺陷，而信息素蒸
发系数p大小的选择将直接影响到整个蚁群算法的收敛速度和全局搜索性能。在蚁群算法的抽象模型中，p表示信息素蒸发系数，1-p则表示信息素持久性系数。因此，p的取值范围应该是0~1之间的一个
数，表示信息素的蒸发程度，它实际上反映了蚂蚁群体中个体之间相互影响的强弱。p过小时，则表示以前搜索过的路径被再次选择的可能性过大，会影响到算法的随机性能和全局搜索能力：p过大时，说
明路径上的信息素挥发的相对变多，虽然可以提高算法的随机搜索性能和全局搜索能力，但过多无用搜索操作势必会降低算法的收敛速度。
蚂蚁数目m
蚁群算法是一种随机搜索算法，与其他模拟进化算法一样，通过多个候选解组成的群体进化过程来寻求最优解，在该过程中不仅需要每个个体的自适应能力，更需要群体之间的相互协作能力。蚁群在搜
索过程中之所以表现出复杂有序的行为，是因为个体之间的信息交流与相互协作起着至关重要的作用。
对于旅行商问题，单个蚂蚁在一次循环中所经过的路径，表现为问题可行解集中的一个解，m只蚂蚁在一次循环中所经过的路径，则表现为问题解集中的一个子集。显然，子集增大(即蚂蚁数量增多)，可以提高蚁群算法的全局搜索能力以及算法的稳定性；但蚂蚁数目增大后，会使大量的曾被搜索过的解(路径)上的信息素的变化趋于平均，信息正反馈的作用不明显，虽然搜索的随机性得到了加强，但收敛速度减慢；反之，子集较小(蚂蚁数量少)，特别是当要处理的问题规模比较大时，会使那些从来未被搜索到的解(路径)上的信息素减小到接近于0，搜索的随机性减弱，虽然收敛速度加快了，但会使算法的全局性能降低，算法的稳定性差，容易出现过早停滞现象。m一般取10~50.
信息素强度Q对算法性能的影响
在蚁群算法中，各个参数的作用实际上是紧密联系的，其中对算法性能起着主要作用的是信息启发式因子α、期望启发式因子β和信息素挥发因子p这三个参数，总信息量(对算法性能的影响有赖于上述三个参数的选取， 以及算法模型的选取。例如， 在ant-cycle模型和ant-quantity模型中， 总信息量4所起的作用显然是有很大差异的， 即随着问题规模的不同，其影响程度也将不同。相关人员研究结果表
明：总信息量Q对ant-cycle模型蚁群算法的性能没有明显的影响。因此，在算法参数的选择上，参数Q不必作特别的考虑，可以任意选取。
最大进化代数G
最大进化代数6是表示蚁群算法运行结束条件的一个参数，表示蚁群算法运行到指定的进化代数之后就停止运行，并将当前群体中的最佳个体作为所求问题的最优解输出。一般6取100~500。

三、部分源代码

clc %清空命令行窗口
clear %从当前工作区中删除所有变量，并将它们从系统内存中释放
close all %删除其句柄未隐藏的所有图窗
tic % 保存当前时间

%% 蚁群算法求解CVRP
%输入：
%City           需求点经纬度
%Distance       距离矩阵
%Demand         各需求点需求量
%AntNum         种群个数
%Alpha          信息素重要程度因子
%Beta           启发函数重要程度因子
%Rho            信息素挥发因子
%Q              常系数
%Eta            启发函数
%Tau            信息素矩阵
%MaxIter        最大迭代次数

%输出：
%bestroute      最短路径
%mindisever     路径长度

%% 加载数据
load('City.mat')	      %需求点经纬度，用于画实际路径的XY坐标
load('Distance.mat')	  %距离矩阵
load('Demand.mat')       %需求量
load('Capacity.mat')     %车容量约束

%% 初始化问题参数
CityNum = size(City,1)-1;    %需求点个数

%% 初始化参数
AntNum = 8;                            % 蚂蚁数量
Alpha = 1;                              % 信息素重要程度因子
Beta = 5;                               % 启发函数重要程度因子
Rho = 0.1;                              % 信息素挥发因子
Q = 1;                                  % 常系数
Eta = 1./Distance;                      % 启发函数
Tau = ones(CityNum+1);                  % (CityNum+1)*(CityNum+1)信息素矩阵  初始化全为1
Population = zeros(AntNum,CityNum*2+1);  % AntNum行,(CityNum*2+1)列 的路径记录表
MaxIter = 100;                           % 最大迭代次数
bestind = ones(1,CityNum*2+1);	% 各代最佳路径
MinDis = zeros(MaxIter,1);              % 各代最佳路径的长度

%% 迭代寻找最佳路径
Iter = 1;                               % 迭代次数初值
while Iter <= MaxIter %当未到达最大迭代次数
	%% 逐个蚂蚁路径选择
    
    maxSubpathLength=zeros(AntNum,1);
	for i = 1:AntNum
        TSProute=2:CityNum+1; %生成一条顺序不包括首尾位的升序TSP路线
        VRProute=[]; %初始化VRP路径
        
        while ~isempty(TSProute)%开辟新的子路径
            subpath=1; %先将配送中心放入子路径起点
            delete=subpath; %delete(end)=1给第一次进入内while的P(k)首项用
            delivery=0; %此子路径的车辆可配送量初始化为零

            while ~isempty(TSProute) %为子路径subpath第二个及以后的位置安排需求点
                %% 计在内while中计算城市间转移概率
                
                P = TSProute; % 为轮盘赌选择建立等于剩余需经过城市数量的长度的向量
                for k = 1:length(TSProute)
                    %delete(end)为刚刚经过的城市，TSProute(k)为剩余可能经过的城市
                    P(k) = Tau(delete(end),TSProute(k))^Alpha * Eta(delete(end),TSProute(k))^Beta; %省略相同分母
                end
                P = P/sum(P);
                % 轮盘赌法选择下一个访问城市
                Pc = cumsum(P); %累加概率
                
                TargetIndex = find(Pc >= rand); %寻找左数第一个大于伪随机数的累加概率的索引
                target = TSProute(TargetIndex(1)); %此索引对应的城市
                %不要强行改变蚂蚁通过轮盘法选到的下一个城市
                %它选到就确定了，然后如果超约束，结束此subpath
                
                %判断容量约束
                if delivery+Demand(target) > Capacity
                    break
                else
                    delivery = delivery+Demand(target); %若符合，则经过的距离累加此距离
                        
                    %此点加入子路径
                    subpath=[subpath,target];
                    %此点加入要删除的点序列
                    delete=[delete,target];
                    
                    %TSP路径中排除已安排的城市
                    TSProute=setdiff(TSProute,delete);
                end
            end %内while结束，此subpath结束
            %直接在VRProute后面填充这条子路径
            maxSubpathLength(i)=max(length(subpath),maxSubpathLength(i));
            VRProute=[VRProute,subpath];
        end %外while结束，此VRProute完整

        %在VRProute后未到CityNum*2+1的空位填入1
        fillwithones = linspace(1,1,CityNum*2+1-length(VRProute));
        VRProute=[VRProute,fillwithones];
        
        %把成型的VRP路径赋给种群矩阵
        Population(i,:)=VRProute;
	end
    function DrawPath(route,City)
%% 画路径函数
%输入
% route     待画路径   
% City      各城市坐标位置

figure
hold on %保留当前坐标区中的绘图，从而使新添加到坐标区中的绘图不会删除现有绘图
box on %通过将当前坐标区的 Box 属性设置为 'on' 在坐标区周围显示框轮廓
xlim([min(City(:,1)-0.01),max(City(:,1)+0.01)]) %手动设置x轴范围  xlimit
ylim([min(City(:,2)-0.01),max(City(:,2)+0.01)]) %手动设置y轴范围
function TextOutput(Distance,Demand,route,Capacity)
%% 输出路径函数
%输入：route 路径
%输出：p 路径文本形式

%% 总路径
len=length(route); %路径长度
disp('Best Route:')

p=num2str(route(1)); %配送中心位先进入路径首位
for i=2:len
    p=[p,' -> ',num2str(route(i))]; %路径依次加入下一个经过的点
end
disp(p)

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.