一、背景知识
(1)算法简介
PSO算法是基于群体的,根据对环境的适应度将群体中的个体移动到好的区域。但是他不是对个体进使用演化算子,而是将每个个体看作是D维探索空间中的一个没有质量、没有体积的一个微粒,在探索空间中以一定的速度飞行,这个速度是根据它本身的飞行经验和同伴的飞行经验来动态调整的。因此基于以上特性,可以将粒子群优化算法假设前提:
- 没有质量、没有体积、有速度的微粒(点);
- 各粒子之间存在某种交流(在这里可以理解为飞行经验及如何实现最优解);
- 速度非恒定值,是动态调整的;
- 会进行迭代(即所谓的“更新”)
(2)算法特点
- 粒子群优化算法具有收敛速度快、参数少、算法简单、效果好的优点;
- 由于算法原理的特性,容易导致局部最优解,而不是全局最优解;
(1)2.各粒子之间存在某种交流(在这里可以理解为飞行经验及如何实现最优解);
(3)使用条件
- 函数
- 约束条件存在(即PSO约束条件)
二、算法原理
(1)鸟群和粒子群算法
鸟群觅食行为和算法原理对应,下图为总览和部分限制条件:
鸟群觅食 | 粒子群算法 |
---|---|
鸟 | 粒子 |
森林、林地、草原 | 求解区域 |
食物总量 | 目标函数值 |
每只鸟所处的具体位置 | 区域中的一个解 |
食物数量最多的位置 | 全局最优解 |
食物数量最少的位置 | 全局最优解 |
地形、食物限制条件 | 维度和解的限制条件 |
---|---|
地形限制 | 区域分布 |
食物种类限制 | 多个解 |
(2)拆分详解
- 1.前提条件:各个微粒(成员)之间存在相互交流关系。
- 2.三个基本属性:无质量、有速度、位置不固定(某种意义上不会是同一位置)
- 3.公式详解:
① 第i个微粒表示为:
②所经历的最好位置(有最好的适应值)为:
③微粒i的速度:
④每一代,它的第d维(1<= d<=D),如下公式:
参数详解:w为惯性权重,c1和c2为加速常数,rand()和Rand()为两个在[0,1]范围里变化的随机值。
(3)参数详解:
-
PSO参数:群体规模m,惯性权重w,加速常数c1和c2,最大速度Vmax,最大代数Gmax。
①Vmax:决定当前位置与最好位置之间的的区域的精度。如果Vmax过大,则微粒将会超过最优解,即鸟飞过最最优食物解,若Vmax过小,则探索域不够,无法保证全局最优解,将会导致局部最优解。(个体最优方向)②惯性权重:权重w使微粒保持运动的惯性,使其有扩展搜索空间的趋势,有能力探索新的区域(惯性方向)
③加速常数c1和c2:代表将每个微粒推向pbest和gbest位置的统计加速项的权重。低的值允许微粒在被拉回来之前可以在目标区域外徘徊,而高的值导致微粒突然的冲向或者越过目标区域。(群体最优方向)
④粒子实际运动方向:
(3)①,(3)②,(3)③,分别讲解了三个方向的来源,在力的作用下,粒子实际方向发生偏转。
三、算法参数的详解
(1)粒子群规模:N
一个正整数,推荐取值范围:[20,1000],简单问题一般取[20,40],较难或特定类别的问题可以取[100,200]。较小的种群规模容易陷入局部最优;较大的种群规模可以提高收敛性,更快找到全局最优解,但是相应地每次迭代的计算量也会增大;当种群规模增大至一定水平时,再增大将不再有显著的作用。
(2)粒子维度:D
粒子探索的空间维度即为自变量的个数。
(3)迭代次数:K
一般取值为60,70,100;根据实际情况进行调整。
(4)惯性权重:W
1998年,Yuhui Shi和Russell Eberhart对基本粒子群算法引入了惯性权重(inertia weight)W,并提出动态调整惯性权重以平衡收敛的全局性和收敛速度,该算法被称为标准PSO算法。
https://xueshu.baidu.com/usercenter/paper/show?paperid=0667417b47e4e05ee3ee3e2c995d6290&site=xueshu_se
四.Matlab实现粒子群优化算法
以下为main函数:
%%author - Christ
%%Date - 2023-3-1
%% 该文件演示基于TSP-PSO算法
clc;clear
%% 下载数据
data=load('eil51.txt');
cityCoor=[data(:,2) data(:,3)];%城市坐标矩阵
figure
plot(cityCoor(:,1),cityCoor(:,2),'ms','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
legend('Light source position')
ylim([4 80])
title('Light source profile','fontsize',12)
xlabel('km','fontsize',12)
ylabel('km','fontsize',12)
%ylim([min(cityCoor(:,2))-1 max(cityCoor(:,2))+1])
grid on
%% 计算城市间距离
n=size(cityCoor,1); %城市数目
cityDist=zeros(n,n); %城市距离矩阵
for i=1:n
for j=1:n
if i~=j
cityDist(i,j)=((cityCoor(i,1)-cityCoor(j,1))^2+...
(cityCoor(i,2)-cityCoor(j,2))^2)^0.5;
end
cityDist(j,i)=cityDist(i,j);
end
end
nMax=200; %进化次数
indiNumber=1000; %个体数目
individual=zeros(indiNumber,n);
%^初始化粒子位置
for i=1:indiNumber
individual(i,:)=randperm(n);
end
%% 计算种群适应度
indiFit=fitness(individual,cityCoor,cityDist);
[value,index]=min(indiFit);
tourPbest=individual; %当前个体最优
tourGbest=individual(index,:) ; %当前全局最优
recordPbest=inf*ones(1,indiNumber); %个体最优记录
recordGbest=indiFit(index); %群体最优记录
xnew1=individual;
%% 循环寻找最优路径
L_best=zeros(1,nMax);
for N=1:nMax
N
%计算适应度值
indiFit=fitness(individual,cityCoor,cityDist);
%更新当前最优和历史最优
for i=1:indiNumber
if indiFit(i)<recordPbest(i)
recordPbest(i)=indiFit(i);
tourPbest(i,:)=individual(i,:);
end
if indiFit(i)<recordGbest
recordGbest=indiFit(i);
tourGbest=individual(i,:);
end
end
[value,index]=min(recordPbest);
recordGbest(N)=recordPbest(index);
%% 交叉操作
for i=1:indiNumber
% 与个体最优进行交叉
c1=unidrnd(n-1); %产生交叉位
c2=unidrnd(n-1); %产生交叉位
while c1==c2
c1=round(rand*(n-2))+1;
c2=round(rand*(n-2))+1;
end
chb1=min(c1,c2);
chb2=max(c1,c2);
cros=tourPbest(i,chb1:chb2);
ncros=size(cros,2);
%删除与交叉区域相同元素
for j=1:ncros
for k=1:n
if xnew1(i,k)==cros(j)
xnew1(i,k)=0;
for t=1:n-k
temp=xnew1(i,k+t-1);
xnew1(i,k+t-1)=xnew1(i,k+t);
xnew1(i,k+t)=temp;
end
end
end
end
%插入交叉区域
xnew1(i,n-ncros+1:n)=cros;
%新路径长度变短则接受
dist=0;
for j=1:n-1
dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
end
dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
if indiFit(i)>dist
individual(i,:)=xnew1(i,:);
end
% 与全体最优进行交叉
c1=round(rand*(n-2))+1; %产生交叉位
c2=round(rand*(n-2))+1; %产生交叉位
while c1==c2
c1=round(rand*(n-2))+1;
c2=round(rand*(n-2))+1;
end
chb1=min(c1,c2);
chb2=max(c1,c2);
cros=tourGbest(chb1:chb2);
ncros=size(cros,2);
%删除与交叉区域相同元素
for j=1:ncros
for k=1:n
if xnew1(i,k)==cros(j)
xnew1(i,k)=0;
for t=1:n-k
temp=xnew1(i,k+t-1);
xnew1(i,k+t-1)=xnew1(i,k+t);
xnew1(i,k+t)=temp;
end
end
end
end
%插入交叉区域
xnew1(i,n-ncros+1:n)=cros;
%新路径长度变短则接受
dist=0;
for j=1:n-1
dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
end
dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
if indiFit(i)>dist
individual(i,:)=xnew1(i,:);
end
%% 变异操作
c1=round(rand*(n-1))+1; %产生变异位
c2=round(rand*(n-1))+1; %产生变异位
while c1==c2
c1=round(rand*(n-2))+1;
c2=round(rand*(n-2))+1;
end
temp=xnew1(i,c1);
xnew1(i,c1)=xnew1(i,c2);
xnew1(i,c2)=temp;
%新路径长度变短则接受
dist=0;
for j=1:n-1
dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));
end
dist=dist+cityDist(xnew1(i,1),xnew1(i,n));
if indiFit(i)>dist
individual(i,:)=xnew1(i,:);
end
end
[value,index]=min(indiFit);
L_best(N)=indiFit(index);
tourGbest=individual(index,:);
end
%% 结果作图
figure
plot(L_best)
title('Algorithm training procedure')
xlabel('iterations')
ylabel('Fitness value')
grid on
figure
hold on
plot([cityCoor(tourGbest(1),1),cityCoor(tourGbest(n),1)],[cityCoor(tourGbest(1),2),...
cityCoor(tourGbest(n),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
hold on
for i=2:n
plot([cityCoor(tourGbest(i-1),1),cityCoor(tourGbest(i),1)],[cityCoor(tourGbest(i-1),2),...
cityCoor(tourGbest(i),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
hold on
end
legend('Planning path')
scatter(cityCoor(:,1),cityCoor(:,2));
title('Planning path','fontsize',10)
xlabel('km','fontsize',10)
ylabel('km','fontsize',10)
grid on
ylim([4 80])