【TWVRP】基于matlab遗传和粒子群算法求解带时间窗的车辆路径规划问题【含Matlab源码 1037期】

一、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 遗传算法概述
遗传算法（Genetic Algorithm，GA）是进化计算的一部分，是模拟达尔文的遗传选择和自然淘汰的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。该算法简单、通用，鲁棒性强，适于并行处理。

2 遗传算法的特点和应用
遗传算法是一类可用于复杂系统优化的具有鲁棒性的搜索算法，与传统的优化算法相比，具有以下特点：
（1）以决策变量的编码作为运算对象。传统的优化算法往往直接利用决策变量的实际值本身来进行优化计算，但遗传算法是使用决策变量的某种形式的编码作为运算对象。这种对决策变量的编码处理方式，使得我们在优化计算中可借鉴生物学中染色体和基因等概念，可以模仿自然界中生物的遗传和进化激励，也可以很方便地应用遗传操作算子。
（2）直接以适应度作为搜索信息。传统的优化算法不仅需要利用目标函数值，而且搜索过程往往受目标函数的连续性约束，有可能还需要满足“目标函数的导数必须存在”的要求以确定搜索方向。遗传算法仅使用由目标函数值变换来的适应度函数值就可确定进一步的搜索范围，无需目标函数的导数值等其他辅助信息。直接利用目标函数值或个体适应度值也可以将搜索范围集中到适应度较高部分的搜索空间中，从而提高搜索效率。
（3）使用多个点的搜索信息，具有隐含并行性。传统的优化算法往往是从解空间的一个初始点开始最优解的迭代搜索过程。单个点所提供的搜索信息不多，所以搜索效率不高，还有可能陷入局部最优解而停滞；遗传算法从由很多个体组成的初始种群开始最优解的搜索过程，而不是从单个个体开始搜索。对初始群体进行的、选择、交叉、变异等运算，产生出新一代群体，其中包括了许多群体信息。这些信息可以避免搜索一些不必要的点，从而避免陷入局部最优，逐步逼近全局最优解。
（4）使用概率搜索而非确定性规则。传统的优化算法往往使用确定性的搜索方法，一个搜索点到另一个搜索点的转移有确定的转移方向和转移关系，这种确定性可能使得搜索达不到最优店，限制了算法的应用范围。遗传算法是一种自适应搜索技术，其选择、交叉、变异等运算都是以一种概率方式进行的，增加了搜索过程的灵活性，而且能以较大概率收敛于最优解，具有较好的全局优化求解能力。但，交叉概率、变异概率等参数也会影响算法的搜索结果和搜索效率，所以如何选择遗传算法的参数在其应用中是一个比较重要的问题。
综上，由于遗传算法的整体搜索策略和优化搜索方式在计算时不依赖于梯度信息或其他辅助知识，只需要求解影响搜索方向的目标函数和相应的适应度函数，所以遗传算法提供了一种求解复杂系统问题的通用框架。它不依赖于问题的具体领域，对问题的种类有很强的鲁棒性，所以广泛应用于各种领域，包括：函数优化、组合优化生产调度问题、自动控制
、机器人学、图像处理（图像恢复、图像边缘特征提取…)、人工生命、遗传编程、机器学习。

3 遗传算法的基本流程及实现技术
基本遗传算法（Simple Genetic Algorithms,SGA）只使用选择算子、交叉算子和变异算子这三种遗传算子，进化过程简单，是其他遗传算法的基础。

3.1 遗传算法的基本流程
通过随机方式产生若干由确定长度（长度与待求解问题的精度有关）编码的初始群体；
通过适应度函数对每个个体进行评价，选择适应度值高的个体参与遗传操作，适应度低的个体被淘汰；
经遗传操作（复制、交叉、变异）的个体集合形成新一代种群，直到满足停止准则（进化代数GEN>=?）；
将后代中变现最好的个体作为遗传算法的执行结果。

其中，GEN是当前代数；M是种群规模，i代表种群数量。

3.2 遗传算法的实现技术
基本遗传算法（SGA）由编码、适应度函数、遗传算子（选择、交叉、变异）及运行参数组成。
3.2.1 编码
（1）二进制编码
二进制编码的字符串长度与问题所求解的精度有关。需要保证所求解空间内的每一个个体都可以被编码。
优点：编、解码操作简单，遗传、交叉便于实现
缺点：长度大
（2）其他编码方法
格雷码、浮点数编码、符号编码、多参数编码等
3.2.2 适应度函数
适应度函数要有效反映每一个染色体与问题的最优解染色体之间的差距。
3.2.3选择算子

3.2.4 交叉算子
交叉运算是指对两个相互配对的染色体按某种方式相互交换其部分基因，从而形成两个新的个体；交叉运算是遗传算法区别于其他进化算法的重要特征，是产生新个体的主要方法。在交叉之前需要将群体中的个体进行配对，一般采取随机配对原则。
常用的交叉方式：
单点交叉
双点交叉（多点交叉，交叉点数越多，个体的结构被破坏的可能性越大，一般不采用多点交叉的方式）
均匀交叉
算术交叉
3.2.5 变异算子
遗传算法中的变异运算是指将个体染色体编码串中的某些基因座上的基因值用该基因座的其他等位基因来替换，从而形成一个新的个体。

就遗传算法运算过程中产生新个体的能力方面来说，交叉运算是产生新个体的主要方法，它决定了遗传算法的全局搜索能力；而变异运算只是产生新个体的辅助方法，但也是必不可少的一个运算步骤，它决定了遗传算法的局部搜索能力。交叉算子与变异算子的共同配合完成了其对搜索空间的全局搜索和局部搜索，从而使遗传算法能以良好的搜索性能完成最优化问题的寻优过程。

3.2.6 运行参数

4 遗传算法的基本原理
4.1 模式定理

4.2 积木块假设
具有低阶、定义长度短，且适应度值高于群体平均适应度值的模式称为基因块或积木块。
积木块假设：个体的基因块通过选择、交叉、变异等遗传算子的作用，能够相互拼接在一起，形成适应度更高的个体编码串。
积木块假设说明了用遗传算法求解各类问题的基本思想，即通过积木块直接相互拼接在一起能够产生更好的解。

三、粒子群算法简介

1 引言
自然界中的鸟群和鱼群的群体行为一直是科学家的研究兴趣所在。生物学家Craig Reynolds在1987年提出了一个非常有影响的鸟群聚集模型，在他的仿真中，每一个个体都遵循：避免与邻域个体相撞：匹配邻域个体的速度；飞向鸟群中心，且整个群体飞向目标。仿真中仅利用上面三条简单的规则，就可以非常接近地模拟出鸟群飞行的现象。1990年，生物学家Frank Heppner也提出了鸟类模型，它的不同之处在于：鸟类被吸引飞到栖息地。在仿真中，一开始每一只鸟都没有特定的飞行目标，只是使用简单的规则确定自己的飞行方向和飞行速度，当有一只鸟飞到栖息地时，它周围的鸟也会跟着飞向栖息地，最终整个鸟群都会落在栖息地。
1995年，美国社会心理学家James Kennedy和电气工程师RussellEberhart共同提出了粒子群算法(ParticleS warm Optimization， PSO) ，该算法的提出是受对鸟类群体行为进行建模与仿真的研究结果的启发。他们的模型和仿真算法主要对Frank Heppner的模型进行了修正，以使粒子飞向解空间并在最优解处降落。粒子群算法一经提出，由于其算法简单，容易实现，立刻引起了进化计算领域学者们的广泛关注，形成一个研究热点。2001年出版的J.Kennedy与R.Eberhart合著的《群体智能》将群体智能的影响进一步扩大[] ，随后关于粒子群优化算法的研究报告和研究成果大量涌现，继而掀起了国内外研究热潮[2-7]。
粒子群优化算法来源于鸟类群体活动的规律性，进而利用群体智能建立一个简化的模型。它模拟鸟类的觅食行为，将求解问题的搜索空间比作鸟类的飞行空间，将每只鸟抽象成一个没有质量和体积的粒
子，用它来表征问题的一个可能解，将寻找问题最优解的过程看成鸟类寻找食物的过程，进而求解复杂的优化问题。粒子群优化算法与其他进化算法一样，也是基于“种群”和“进化”的概念，通过个体间
的协作与竞争，实现复杂空间最优解的搜索。同时，它又不像其他进化算法那样对个体进行交叉、变异、选择等进化算子操作，而是将群体中的个体看作在l维搜索空间中没有质量和体积的粒子，每个粒子以一定的速度在解空间运动，并向自身历史最佳位置P best和邻域历史最佳位置g best聚集，实现对候选解的进化。粒子群算法具有很好的生物社会背景而易于理解，由于参数少而容易实现，对非线性、多峰问题均具有较强的全局搜索能力，在科学研究与工程实践中得到了广泛关注。目前，该算法已广泛应用于函数优化、神经网络训练、模式分类、模糊控制等领域。

2 粒子群算法理论
2.1粒子群算法描述
鸟类在捕食过程中，鸟群成员可以通过个体之间的信息交流与共享获得其他成员的发现与飞行经历。在食物源零星分布并且不可预测的条件下，这种协作机制所带来的优势是决定性的，远远大于对食物
的竞争所引起的劣势。粒子群算法受鸟类捕食行为的启发并对这种行为进行模仿，将优化问题的搜索空间类比于鸟类的飞行空间，将每只鸟抽象为一个粒子，粒子无质量、无体积，用以表征问题的一个可行解，优化问题所要搜索到的最优解则等同于鸟类寻找的食物源。粒子群算法为每个粒子制定了与鸟类运动类似的简单行为规则，使整个粒子群的运动表现出与鸟类捕食相似的特性，从而可以求解复杂的优化问题。
粒子群算法的信息共享机制可以解释为一种共生合作的行为，即每个粒子都在不停地进行搜索，并且其搜索行为在不同程度上受到群体中其他个体的影响[8]，同时这些粒子还具备对所经历最佳位置的记
忆能力，即其搜索行为在受其他个体影响的同时还受到自身经验的引导。基于独特的搜索机制，粒子群算法首先生成初始种群，即在可行解空间和速度空间随机初始化粒子的速度与位置，其中粒子的位置用于表征问题的可行解，然后通过种群间粒子个体的合作与竞争来求解优化问题。
2.2粒子群算法建模
粒子群优化算法源自对鸟群捕食行为的研究：一群鸟在区域中随机搜索食物，所有鸟知道自己当前位置离食物多远，那么搜索的最简单有效的策略就是搜寻目前离食物最近的鸟的周围区域。粒子群算法
利用这种模型得到启示并应用于解决优化问题。在粒子群算法中，每个优化问题的潜在解都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化的函数决定的适应度值，每个粒子还有一个速度决定它们飞翔的方向和距离。然后，粒子们就追随当前的最优粒子在解空间中搜索[9]。

粒子群算法首先在给定的解空间中随机初始化粒子群，待优化问题的变量数决定了解空间的维数。每个粒子有了初始位置与初始速度，然后通过迭代寻优。在每一次迭代中，每个粒子通过跟踪两个“极值”来更新自己在解空间中的空间位置与飞行速度：一个极值就是单个粒子本身在迭代过程中找到的最优解粒子，这个粒子叫作个体极值：另一个极值是种群所有粒子在迭代过程中所找到的最优解粒子，这个粒子是全局极值。上述的方法叫作全局粒子群算法。如果不用种群所有粒子而只用其中一部分作为该粒子的邻居粒子，那么在所有邻居粒子中的极值就是局部极值，该方法称为局部粒子群算法。

2.3粒子群算法的特点
粒子群算法本质是一种随机搜索算法，它是一种新兴的智能优化技术。该算法能以较大概率收敛于全局最优解。实践证明，它适合在动态、多目标优化环境中寻优，与传统优化算法相比，具有较快的计
算速度和更好的全局搜索能力。
(1)粒子群算法是基于群智能理论的优化算法，通过群体中粒子间的合作与竞争产生的群体智能指导优化搜索。与其他算法相比，粒子群算法是一种高效的并行搜索算法。
(2)粒子群算法与遗传算法都是随机初始化种群，使用适应值来评价个体的优劣程度和进行一定的随机搜索。但粒子群算法根据自己的速度来决定搜索，没有遗传算法的交叉与变异。与进化算法相比，粒子群算法保留了基于种群的全局搜索策略，但是其采用的速度-位移模型操作简单，避免了复杂的遗传操作。
(3)由于每个粒子在算法结束时仍保持其个体极值，即粒子群算法除了可以找到问题的最优解外，还会得到若干较好的次优解，因此将粒子群算法用于调度和决策问题可以给出多种有意义的方案。
(4)粒子群算法特有的记忆使其可以动态地跟踪当前搜索情况并调整其搜索策略。另外，粒子群算法对种群的大小不敏感，即使种群数目下降时，性能下降也不是很大。

3 粒子群算法种类
3.1基本粒子群算法

3.2标准粒子群算法
引入研究粒子群算法经常用到的两个概念：一是“探索”，指粒子在一定程度上离开原先的搜索轨迹，向新的方向进行搜索，体现了一种向未知区域开拓的能力，类似于全局搜索；二是“开发”，指粒子在一定程度上继续在原先的搜索轨迹上进行更细一步的搜索，主要指对探索过程中所搜索到的区域进行更进一步的搜索。探索是偏离原来的寻优轨迹去寻找一个更好的解，探索能力是一个算法的全局搜索能力。开发是利用一个好的解，继续原来的寻优轨迹去搜索更好的解，它是算法的局部搜索能力。如何确定局部搜索能力和全局搜索能力的比例，对一个问题的求解过程很重要。1998年， Shi Yuhui等人提出了带有惯性权重的改进粒子群算法[10]，由于该算法能够保证较好的收敛效果，所以被默认为标准粒子群算法。其进化过程为：

在式(6.7)中，第一部分表示粒子先前的速度，用于保证算法的全局收敛性能；第二部分、第三部分则使算法具有局部收敛能力。可以看出，式(6.7)中惯性权重w表示在多大程度上保留原来的速度：W
较大，则全局收敛能力较强，局部收敛能力较弱；w较小，则局部收敛能力较强，全局收敛能力较弱。
当w=1时，式(6.7)与式(6.5)完全一样，表明带惯性权重的粒子群算法是基本粒子群算法的扩展。实验结果表明：w在0.8~1.2之间时，粒子群算法有更快的收敛速度；而当w>1.2时，算法则容易陷入局部极值。
另外，在搜索过程中可以对w进行动态调整：在算法开始时，可给w赋予较大正值，随着搜索的进行，可以线性地使w逐渐减小，这样可以保证在算法开始时，各粒子能够以较大的速度步长在全局范围内探
测到较好的区域；而在搜索后期，较小的w值则保证粒子能够在极值点周围做精细的搜索，从而使算法有较大的概率向全局最优解位置收敛。对w进行调整，可以权衡全局搜索和局部搜索能力。目前，采用较多的动态惯性权重值是Shi提出的线性递减权值策略，其表达式如下：

3.3压缩因子粒子群算法
Clerc等人提出利用约束因子来控制系统行为的最终收敛[11] ，该方法可以有效搜索不同的区域，并且能得到高质量的解。压缩因子法的速度更新公式为：

实验结果表明：与使用惯性权重的粒子群优化算法相比，使用具
有约束因子的粒子群算法具有更快的收敛速度。
3.4离散粒子群算法
基本的粒子群算法是在连续域中搜索函数极值的有力工具。继基本粒子群算法之后， Kennedy和Eberhart又提出了一种离散二进制版的粒子群算法[12]。在此离散粒子群方法中，将离散问题空间映射到连续粒子运动空间，并适当修改粒子群算法来求解，在计算上仍保留经典粒子群算法速度-位置更新运算规则。粒子在状态空间的取值和变化只限于0和1两个值，而速度的每一维vi y代表位置每一位xi取值为1的可能性。因此，在连续粒子群中的vij更新公式依然保持不变，但是P best和：best只在[0， 1] 内取值。其位置更新等式表示如下：

4 粒子群算法流程
粒子群算法基于“种群”和“进化”的概念，通过个体间的协作与竞争，实现复杂空间最优解的搜索[13]，其流程如下：
(1)初始化粒子群，包括群体规模N，每个粒子的位置x；和速度Vio
(2) 计算每个粒子的适应度值fit[i] 。
(3) 对每个粒子，用它的适应度值fit[门和个体极值P best(i)比较。如果fit[i] <P best(i) ，则用fit[i] 替换掉P best(i) 。
(4) 对每个粒子，用它的适应度值fit[i] 和全局极值g best比较。如果fit[i] < 8 best，则用fit[i] 替换g best。
(5)迭代更新粒子的速度v；和位置xj。
(6)进行边界条件处理。
(7)判断算法终止条件是否满足：若是，则结束算法并输出优化结果；否则返回步骤(2)。
粒子群算法的运算流程如图6.1所示。

5 关键参数说明
在粒子群优化算法中，控制参数的选择能够影响算法的性能和效率；如何选择合适的控制参数使算法性能最佳，是一个复杂的优化问题。在实际的优化问题中，通常根据使用者的经验来选取控制参数。
粒子群算法的控制参数主要包括：粒子种群规模N，惯性权重w，加速系数c和c，最大速度Via x，停止准则，邻域结构的设定，边界条件处理策略等[14]，
粒子种群规模N
粒子种群大小的选择视具体问题而定，但是一般设置粒子数为20~50。对于大部分的问题10个粒子，已经可以取得很好的结果：不过对于比较难的问题或者特定类型的问题，粒子的数量可以取到100或
200。另外，粒子数目越大，算法搜索的空间范围就越大，也就更容易发现全局最优解；当然，算法运行的时间也越长。
惯性权重w
惯性权重w是标准粒子群算法中非常重要的控制参数，可以用来控制算法的开发和探索能力。惯性权重的大小表示了对粒子当前速度继承的多少。当惯性权重值较大时，全局寻优能力较强，局部寻优能力
较弱：当惯性权重值较小时，全局寻优能力较弱，局部寻优能力较强。惯性权重的选择通常有固定权重和时变权重。固定权重就是选择常数作为惯性权重值，在进化过程中其值保持不变，一般取值为
[0.8，1.2]：时变权重则是设定某一变化区间，在进化过程中按照某种方式逐步减小惯性权重。时变权重的选择包括变化范围和递减率。固定的惯性权重可以使粒子保持相同的探索和开发能力，而时变权重可以使粒子在进化的不同阶段拥有不同的探索和开发能力。
加速常数c1和c2
加速常数c和c 2分别调节向P best和g best方向飞行的最大步长，它们分别决定粒子个体经验和群体经验对粒子运行轨迹的影响，反映粒子群之间的信息交流。如果cr=c2=0，则粒子将以当前的飞行速度飞到边界。此时，粒子仅能搜索有限的区域，所以难以找到最优解。如果q=0，则为“社会”模型，粒子缺乏认知能力，而只有群体经验，它的收敛速度较快，但容易陷入局部最优；如果oy=0，则为“认知”模
型，没有社会的共享信息，个体之间没有信息的交互，所以找到最优解的概率较小，一个规模为D的群体等价于运行了N个各行其是的粒子。因此一般设置c1=C2，通常可以取c1=cg=1.5。这样，个体经验和群体经验就有了同样重要的影响力，使得最后的最优解更精确。
粒子的最大速度vmax
粒子的速度在空间中的每一维上都有一个最大速度限制值vd max，用来对粒子的速度进行钳制，使速度控制在范围[-Vimax， +va max] 内，这决定问题空间搜索的力度，该值一般由用户自己设定。Vmax是一个非常重要的参数，如果该值太大，则粒子们也许会飞过优秀区域：而如果该值太小，则粒子们可能无法对局部最优区域以外的区域进行充分的探测。它们可能会陷入局部最优，而无法移动足够远的距离而跳出局部最优，达到空间中更佳的位置。研究者指出，设定Vmax和调整惯性权重的作用是等效的，所以!max一般用于对种群的初始化进行设定，即将vmax设定为每维变量的变化范围，而不再对最大速度进行细致的选择和调节。
停止准则
最大迭代次数、计算精度或最优解的最大停滞步数▲t(或可以接受的满意解)通常认为是停止准则，即算法的终止条件。根据具体的优化问题，停止准则的设定需同时兼顾算法的求解时间、优化质量和
搜索效率等多方面性能。
邻域结构的设定
全局版本的粒子群算法将整个群体作为粒子的邻域，具有收敛速度快的优点，但有时算法会陷入局部最优。局部版本的粒子群算法将位置相近的个体作为粒子的邻域，收敛速度较慢，不易陷入局部最优
值。实际应用中，可先采用全局粒子群算法寻找最优解的方向，即得到大致的结果，然后采用局部粒子群算法在最优点附近进行精细搜索。
边界条件处理
当某一维或若干维的位置或速度超过设定值时，采用边界条件处理策略可将粒子的位置限制在可行搜索空间内，这样能避免种群的膨胀与发散，也能避免粒子大范围地盲目搜索，从而提高了搜索效率。
具体的方法有很多种，比如通过设置最大位置限制Xmax和最大速度限制Vmax，当超过最大位置或最大速度时，在范围内随机产生一个数值代替，或者将其设置为最大值，即边界吸收。

四、部分源代码

clc %清空命令行窗口
clear %从当前工作区中删除所有变量，并将它们从系统内存中释放
close all %删除其句柄未隐藏的所有图窗
tic % 保存当前时间
%% GA-PSO算法求解VRPTW
%输入：
%City           需求点经纬度
%Distance       距离矩阵
%Demand         各点需求量
%Travelcon      行程约束
%Capacity       车容量约束
%TimeWindow     各需求点时间窗%NIND           种群个数
%MAXGEN         遗传到第MAXGEN代时程序停止%输出：
%Gbest          最短路径
%GbestDistance  最短路径长度%% 加载数据
load('City.mat')        %需求点经纬度，用于画实际路径的XY坐标
load('Distance.mat')    %距离矩阵
load('TravelTime.mat')   %行驶时间矩阵
load('Demand.mat')       %需求量
load('TimeWindow.mat')   %时间窗
load('Travelcon.mat')   %行程约束
load('Capacity.mat')     %车容量约束%% 初始化问题参数
CityNum=size(City,1)-1; %需求点个数%% 初始化算法参数
NIND=60; %粒子数量
MAXGEN=100; %最大迭代次数%% 为预分配内存而初始化的0矩阵
Population = zeros(NIND,CityNum*2+1); %预分配内存，用于存储种群数据
PopDistance = zeros(NIND,1); %预分配矩阵内存
GbestDisByGen = zeros(MAXGEN,1); %预分配矩阵内存for i = 1:NIND%% 初始化各粒子Population(i,:)=InitPop(CityNum,Distance,TravelTime,Demand,TimeWindow,Travelcon,Capacity); %随机生成TSP路径%% 计算路径长度PopDistance(i) = CalcDis(Population(i,:),Distance,TravelTime,Demand,TimeWindow,Travelcon,Capacity); % 计算路径长度
endfunction ttlDistance=CalcDis(route,Distance,TravelTime,Demand,TimeWindow,Travelcon,Capacity)
%% 计算各个体的路径长度 适应度函数
%输入：
%route          某粒子代表的路径
%Distance       两两城市之间的距离矩阵
%TravelTime    行驶时间矩阵
%Demand         各点需求量
%TimeWindow    早晚时间窗
%Travelcon      行程约束
%Capacity       容量约束%输出：
%ttlDistance    种群个体路径距离%相关数据初始化
Dis=0; %此方案所有车辆的总距离置零
DisTraveled=0;  % 所有车辆已经行驶的距离置零
CurrentTime=0; % 所有车辆行驶时间置零
delivery=0; % 所有车辆已经送货量，即已经到达点的需求量之和for j=2:length(route)DisTraveled = DisTraveled+Distance(route(j-1),route(j)); %每两点间距离累加if route(j) == 1  %若此点是配送中心if DisTraveled > TravelconDis = Inf;  %对非可行解进行惩罚breakelseDis=Dis+DisTraveled; %累加已行驶距离DisTraveled=0; %已行驶距离置零CurrentTime=0; %时间置零TimeWindow(1,2)delivery=0; %可配送量置零endunction a=Crossover(a,b)
%% PMX方法交叉
%输入：
%a  粒子代表的路径
%b  个体最优粒子代表的路径%输出：
%a  交叉后的粒子代表的路径L=length(a); %获取亲代染色体长度
r1=unidrnd(L); %在1~L中随机选一整数
r2=unidrnd(L); %在1~L中随机选一整数s=min([r1,r2]); %起点为较小值
e=max([r1,r2]); %终点为较大值b0=b(s:e); %用于最后插入aa=a(s:e); %用于检查重复元素
bb=b(s:e); %用于检查重复元素a(s:e)=[];     %去除交叉部分   去除后后面的元素会往前移outlen=length(a); %去除交叉部分后，  a，b的长度   length outside cross section
inlen=e-s+1;      %交叉部分长度  length inside cross section for i=1:inlen     %交叉部分去除相同元素for j=1:inlenif aa(i)==bb(j) %若交叉部分内上下有相同元素aa(i)=0; %删去bb(j)=0;break  %  break能保证最后aa和bb一样长 且无重复元素endend
end

五、运行结果

六、matlab版本及参考文献

1 matlab版本
2014a

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