一、粒子群算法简介

粒子群算法是在1995年由Eberhart博士和Kennedy博士一起提出的，它源于对鸟群捕食行为的研究。它的基本核心是利用群体中的个体对信息的共享从而使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的最优解。设想这么一个场景：一群鸟进行觅食，而远处有一片玉米地，所有的鸟都不知道玉米地到底在哪里，但是它们知道自己当前的位置距离玉米地有多远。那么找到玉米地的最佳策略，也是最简单有效的策略就是搜寻目前距离玉米地最近的鸟群的周围区域。

在PSO中，每个优化问题的解都是搜索空间中的一只鸟，称之为"粒子"，而问题的最优解就对应于鸟群中寻找的"玉米地"。所有的粒子都具有一个位置向量（粒子在解空间的位置）和速度向量（决定下次飞行的方向和速度），并可以根据目标函数来计算当前的所在位置的适应值（fitness value）,可以将其理解为距离"玉米地"的距离。在每次的迭代中，种群中的例子除了根据自身的经验（历史位置）进行学习以外，还可以根据种群中最优粒子的"经验"来学习，从而确定下一次迭代时需要如何调整和改变飞行的方向和速度。就这样逐步迭代，最终整个种群的例子就会逐步趋于最优解。

上面的解释可能还比较抽象，下面通过一个简单的例子来进行说明

在一个湖中有两个人他们之间可以通信，并且可以探测到自己所在位置的最低点。初始位置如上图所示，由于右边比较深，因此左边的人会往右边移动一下小船。

现在左边比较深，因此右边的人会往左边移动一下小船

 

一直重复该过程，最后两个小船会相遇

得到一个局部的最优解

 

将每个个体表示为粒子。每个个体在某一时刻的位置表示为，x(t),方向表示为v(t)

 

 

 

p（t）为在t时刻x个体的自己的最优解，g(t)为在t时刻所有个体的最优解，v(t)为个体在t时刻的方向，x(t)为个体在t时刻的位置

 

 

 

下一个位置为上图所示由x,p,g共同决定了

 

 

 

种群中的粒子通过不断地向自身和种群的历史信息进行学习，从而可以找到问题的最优解。

但是，在后续的研究中表表明，上述原始的公式中存在一个问题：公式中V的更新太具有随机性，从而使整个PSO算法的全局优化能力很强，但是局部搜索能力较差。而实际上，我们需要在算法迭代初期PSO有着较强的全局优化能力，而在算法的后期，整个种群应该具有更强的局部搜索能力。所以根据上述的弊端，shi和Eberhart通过引入惯性权重修改了公式，从而提出了PSO的惯性权重模型：

 

每一个向量的分量表示如下

 

 

其中w称为是PSO的惯性权重，它的取值介于【0,1】区间，一般应用中均采用自适应的取值方法，即一开始令w=0.9，使得PSO全局优化能力较强，随着迭代的深入，参数w进行递减，从而使的PSO具有较强的局部优化能力，当迭代结束时，w=0.1。参数c1和c2称为学习因子，一般设置为1,4961；而r1和r2为介于[0,1]之间的随机概率值。

整个粒子群优化算法的算法框架如下：

step1种群初始化，可以进行随机初始化或者根据被优化的问题设计特定的初始化方法，然后计算个体的适应值，从而选择出个体的局部最优位置向量和种群的全局最优位置向量。

step2 迭代设置：设置迭代次数，并令当前迭代次数为1

step3 速度更新：更新每个个体的速度向量

step4 位置更新：更新每个个体的位置向量

step5 局部位置和全局位置向量更新：更新每个个体的局部最优解和种群的全局最优解

step6 终止条件判断：判断迭代次数时都达到最大迭代次数，如果满足，输出全局最优解，否则继续进行迭代，跳转至step 3。

　　对于粒子群优化算法的运用，主要是对速度和位置向量迭代算子的设计。迭代算子是否有效将决定整个PSO算法性能的优劣，所以如何设计PSO的迭代算子是PSO算法应用的研究重点和难点。

二 、VRP模型

(1)车辆路径规划问题介绍

车辆路径规划问题，经过60年来的研究与发展，研究的目标对象，限制条件等均有所变化，已经从最初的简单车辆安排调度问题转变为复杂的系统问题。最初的车辆路径规划问题可以描述为：有一个起点和若干个客户点,已知各点的地理位置和需求,在满足各种约束的条件下,如何规划最优的路径,使其能服务到每个客户点,最后返回起点。通过施加不同的约束条件，改变优化的目标，可以衍生出不同种类的车辆路径规划问题。同时车辆路径规划问题属于典型的NP-hard问题，其精确算法能求解的规模很小，故启发式算法也就成了研究热点。

(2)VRPTW简介

VRPTW(Vehicle routing problem with time windows)即带时间窗的车辆路径规划问题，其对于每一需求点加入了时间窗的约束，即对于每一个需求点，设定服务开始的最早时间和最晚时间，要求车辆在时间窗内开始服务顾客。

需求点的时窗限制可以分为两种，一种是硬时间窗(Hard Time Window)，即要求车辆必须在时间窗内开始服务顾客，早到必须等待，迟到就拒收，另一种是软时间窗(Soft Time Window)，不一定要在时间窗内开始服务顾客，但是在时间窗外开始服务必须要惩罚，以惩罚代替等待与拒收是软时间窗和硬时时间窗的最大的区别。

VRPTW的数学模型如下：

 

（2.2）保证了每个顾客只被访问1次

（2.3）保证了装载的货物不超过容量

（2.4）（2.5）（2.6）确保了每辆车从depot出发最后回到depot

（2.7）（2.8）确保在时间窗内开始服务

三、代码

clc %清空命令行窗口
clear %从当前工作区中删除所有变量，并将它们从系统内存中释放
close all %删除其句柄未隐藏的所有图窗
tic % 保存当前时间
%% GA-PSO算法求解CDVRP
%输入：
%City           需求点经纬度
%Distance       距离矩阵
%Demand         各点需求量
%Travelcon      行程约束
%Capacity       车容量约束
%NIND           种群个数
%MAXGEN         遗传到第MAXGEN代时程序停止

%输出：
%Gbest          最短路径
%GbestDistance	最短路径长度

%% 加载数据
load('../test_data/City.mat')	      %需求点经纬度，用于画实际路径的XY坐标
load('../test_data/Distance.mat')	  %距离矩阵
load('../test_data/Demand.mat')       %需求量
load('../test_data/Travelcon.mat')	  %行程约束
load('../test_data/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,Demand,Travelcon,Capacity); %随机生成TSP路径
    
    %% 计算路径长度
    PopDistance(i) = CalcDis(Population(i,:),Distance,Demand,Travelcon,Capacity); % 计算路径长度
end

%% 存储Pbest数据
Pbest = Population; % 初始化Pbest为当前粒子集合
PbestDistance = PopDistance; % 初始化Pbest的目标函数值为当前粒子集合的目标函数值

%% 存储Gbest数据
[mindis,index] = min(PbestDistance); %获得Pbest中
Gbest = Pbest(index,:); % 初始Gbest粒子
GbestDistance = mindis; % 初始Gbest粒子的目标函数值

%% 开始迭代
gen=1;

while gen <= MAXGEN
    %% 每个粒子更新
    for i=1:NIND
        %% 粒子与Pbest交叉
        Population(i,2:end-1)=Crossover(Population(i,2:end-1),Pbest(i,2:end-1)); %交叉
        
        % 新路径长度变短则记录至Pbest
        PopDistance(i) = CalcDis(Population(i,:),Distance,Demand,Travelcon,Capacity); %计算距离
        if PopDistance(i) < PbestDistance(i) %若新路径长度变短
            Pbest(i,:)=Population(i,:); %更新Pbest
            PbestDistance(i)=PopDistance(i); %更新Pbest距离
        end

    %% 更新迭代次数
    gen=gen+1;
end

%删去路径中多余1
for i=1:length(Gbest)-1
    if Gbest(i)==Gbest(i+1)
        Gbest(i)=0;  %相邻位都为1时前一个置零
    end
end
Gbest(Gbest==0)=[];  %删去多余零元素

Gbest=Gbest-1;  % 编码各减1，与文中的编码一致

%% 计算结果数据输出到命令行
disp('-------------------------------------------------------------')
toc %显示运行时间
fprintf('Total Distance = %s km \n',num2str(GbestDistance))
TextOutput(Gbest,Distance,Demand,Capacity)  %显示最优路径
disp('-------------------------------------------------------------')

%% 迭代图
figure
plot(GbestDisByGen,'LineWidth',2) %展示目标函数值历史变化
xlim([1 gen-1]) %设置 x 坐标轴范围
set(gca, 'LineWidth',1)
xlabel('Iterations')
ylabel('Min Distance(km)')
title('HPSO Process')

%% 绘制实际路线
DrawPath(Gbest,City)

代码下载https://www.cnblogs.com/ttmatlab/p/14882966.html

四、结果