粒子群算法
粒子群算法是在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算法应用的研究重点和难点。
%
tic
clear
clc
%% 用importdata这个函数来读取文件
r101=importdata('r101.txt');
cap=200; %车辆最大装载量
%% 提取数据信息
E=r101(1,5); %配送中心时间窗开始时间
L=r101(1,6); %配送中心时间窗结束时间
vertexs=r101(:,2:3); %所有点的坐标x和y
customer=vertexs(2:end,:); %顾客坐标
cusnum=size(customer,1); %顾客数
v_num=25; %车辆最多使用数目
demands=r101(2:end,4); %需求量
a=r101(2:end,5); %顾客时间窗开始时间[a[i],b[i]]
b=r101(2:end,6); %顾客时间窗结束时间[a[i],b[i]]
s=r101(2:end,7); %客户点的服务时间
h=pdist(vertexs);
dist=squareform(h); %距离矩阵
%% 粒子群算法参数
alpha=10; %违反的容量约束的惩罚函数系数
belta=100; %违反时间窗约束的惩罚函数系数
NIND=50; %粒子数目
MAXGEN=100; %迭代次数
w=1; %惯性因子
wdamp=0.99; %惯性因子衰减率
c1=1.5; %个体学习因子
c2=2.0; %全局学习因子
XvMin=1; %Xv下限
XvMax=v_num; %Xv上限
XrMin=1; %Xr下限
XrMax=cusnum; %Xr上限
VvMin=-(v_num-1); %Vv下限
VvMax=v_num-1; %Vv上限
VrMin=-(cusnum-1); %Vr下限
VrMax=cusnum-1; %Vr上限
%% 初始化粒子群位置
init_vc=init(cusnum,a,demands,cap); %构造初始解
population=initpopCW(init_vc,NIND,cusnum,XvMin,XvMax,XrMin,XrMax,VvMin,VvMax,VrMin,VrMax);
ObjV=calObj(population,v_num,cap,demands,a,b,L,s,dist,alpha,belta); %计算各个粒子的目标函数值
gbest_pos=population{1,1}(1:2,:); %假设第一个粒子位置为全局最优位置
gbest_obj=ObjV(1,1); %第一个粒子位置的目标函数值
pbest_pos=cell(NIND,1); %初始化各个粒子的个体最优位置
pbest_obj=ObjV; %初始化各个粒子的个体最优的目标函数值
for i=1:NIND
particle=population{i,1}; %第i个粒子
position=particle(1:2,:); %第i个粒子的位置
pbest_pos{i,1}=position; %初始化这个粒子的个体最优
if pbest_obj(i,1)<gbest_obj
%更新初始种群中的全局最优粒子
gbest_obj=pbest_obj(i,1);
gbest_pos=position;
end
end
globalVC=decode(gbest_pos,v_num); %初始全局最优配送方案
%统计一个配送方案的总成本、车辆使用数目、行驶总距离、违反约束路径数目、违反约束顾客数目
[cost,NV,TD,vio_NV,vio_cus]=statistic(globalVC,a,b,s,L,dist,demands,cap,alpha,belta);
disp(['初始全局最优解总成本为：',num2str(cost),'，车辆使用数目为：',num2str(NV),...
'，行驶总距离为：',num2str(TD),'，违反约束路径数目为：',num2str(vio_NV),'，违反约束顾客数目为：',num2str(vio_cus)]);
BestCost=zeros(MAXGEN,1); %记录每次迭代的全局最优目标函数值
%% 主循环
gen=1; %计数器
while gen<=MAXGEN
%% 更新各个粒子的位置和速度
for i=1:NIND
particle=population{i,1}; %第i个粒子
position=particle(1:2,:); %第i个粒子的位置
Xv=position(1,:);
Xr=position(2,:);
velocity=particle(3:4,:); %第i个粒子的速度
Vv=velocity(1,:);
Vr=velocity(2,:);
%% 更新速度
velocity=w*velocity+ +c1*rand([2,cusnum]).*(pbest_pos{i,1}-position)...
+c2*rand([2,cusnum]).*(gbest_pos-position);
%% 速度越界处理
velocity(1,:)=max(velocity(1,:),VvMin);
velocity(1,:)=min(velocity(1,:),VvMax);
velocity(2,:)=max(velocity(2,:),VrMin);
velocity(2,:)=min(velocity(2,:),VrMax);
%% 更新位置
position=position+velocity;
position(1,:)=ceil(position(1,:)); %对Xv向上取整
%% 速度镜像影响
IsOutside=(position(1,:)<XvMin | position(1,:)>XvMax | position(2,:)<XrMin | position(2,:)>XrMax);
velocity(IsOutside)=-velocity(IsOutside);
%% 位置越界处理
position(1,:)=max(position(1,:),XvMin);
position(1,:)=min(position(1,:),XvMax);
position(2,:)=max(position(2,:),XrMin);
position(2,:)=min(position(2,:),XrMax);
%% 对第i个粒子解码出的配送方案进行relocate操作
VC=decode(position,v_num);
RVC=relocate(VC,cap,demands,a,b,L,s,dist,alpha,belta);
position=change(RVC,cusnum);
%% 计算第i个粒子目标函数值
cost=costFuction(RVC,a,b,s,L,dist,demands,cap,alpha,belta);
%% 更新个体最优
if cost<pbest_obj(i,1)
pbest_pos{i,1}=position;
pbest_obj(i,1)=cost;
%% 更新全局最优
if pbest_obj(i,1)<gbest_obj
gbest_pos=pbest_pos{i,1};
gbest_obj=pbest_obj(i,1);
end
end
%% 更新第i个粒子的速度和位置
particle=[position;velocity];
population{i,1}=particle;
end
%% 记录各代全局最优解，打印各代全局最优解
BestCost(gen)=gbest_obj;
globalVC=decode(gbest_pos,v_num); %初始全局最优配送方案
%统计一个配送方案的总成本、车辆使用数目、行驶总距离、违反约束路径数目、违反约束顾客数目
[cost,NV,TD,vio_NV,vio_cus]=statistic(globalVC,a,b,s,L,dist,demands,cap,alpha,belta);
disp(['第',num2str(gen),'代全局最优解总成本为：',num2str(cost),'，车辆使用数目为：',num2str(NV),...
'，行驶总距离为：',num2str(TD),'，违反约束路径数目为：',num2str(vio_NV),'，违反约束顾客数目为：',num2str(vio_cus)]);
w=w*wdamp;
%% 更新计数器
gen=gen+1;
end
%% 将全局最优粒子解码为全局最优配送方案
globalVC=decode(gbest_pos,v_num);
%% 对全局最优配送方案进行局部搜索
globalVC=relocate_gbest(globalVC,cap,demands,a,b,L,s,dist,alpha,belta);
%% 统计全局最优配送方案的总成本、车辆使用数目、行驶总距离、违反约束路径数目、违反约束顾客数目
[cost,NV,TD,vio_NV,vio_cus]=statistic(globalVC,a,b,s,L,dist,demands,cap,alpha,belta);
disp(['最终全局最优解总成本为：',num2str(cost),'，车辆使用数目为：',num2str(NV),...
'，行驶总距离为：',num2str(TD),'，违反约束路径数目为：',num2str(vio_NV),'，违反约束顾客数目为：',num2str(vio_cus)]);
%% 画出全局最优配送方案路线图
draw_Best(globalVC,vertexs);
%% Results
figure;
%plot(BestCost,'LineWidth',2);
semilogy(BestCost,'LineWidth',2);
xlabel('迭代次数');
ylabel('全局最优目标函数值');
grid on;
toc
