粒子群优化算法 PSO
文章目录
- 1、基本内容
- 1.1 概念
- 1.2 基本原理
- 1.3 基本步骤
- 2、代码分析
- 2.1 Rastrigrin()函数
- 2.2 Schaffer()函数
- 2.3 Griewank()函数
- 2.4 fun()函数
- 2.5 fun1()函数
- 2.6 Drawfunc()函数
- 2.7 PSO.m
- 2.8 PSOMutation.m
- 2.9 运行结果
- 2.9.1 使用PSO.m
- 2.9.2 使用PSOMutation.m
- 3、参数调整（针对POS.m文件）
- 3.1 w改变，其他参数不变
- 3.1.1 w初始值为0.8
- 3.1.2 w=abs(w-rand())
- 3.2 c1、c2改变，其他参数不变
- 3.3 sizepop改变，其他参数不变
- 3.4 dim改变，其他参数不变
1、基本内容
1.1 概念
在动物的群体行为中，科学家们很早就发现了自然界的鸟群、兽群、鱼群等在其迁徙、捕食过程中，往往表现出高度的组织性和规律性。这些现象受到了高度的重视和广泛的关注，吸引着大批生物学家、动物学家、计算机科学家、行为学家和社会心理学家等的深入研究。
粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法或微粒群优化算法，是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。 通常认为粒子群优化算法是群集智能 (Swarm intelligence, SI) 的一种，它可以被纳入多主体优化系统(Multiagent Optimization System, MAOS)。
粒子群优化算法是由Eberhart博士和Kennedy博士于1995年提出的一种全局搜索算法。
1.2 基本原理
事实上，捕食的鸟群都是通过各自的探索与群体的合作最终发现食物所在的位置的。
- 一群分散的鸟在随机地飞行觅食，它们不知道食物所在的具体位置。
- 但各个小鸟会在飞行的过程中不断地记录和更新它曾经到达的离食物最近位置，同时通过信息交流的方式比较大家所找到的最好位置，得到一个当前整个群体已经找到的最佳位置。
- 各个小鸟在飞行过程中结合自身的经验和整个群体的经验，调整自己的飞行速度和所在位置，不断地寻找更加接近食物的位置，最终使得群体聚集到食物位置。
因此，在粒子群优化算法中，
- 每只鸟抽象为一个无质量、无体积的 “ 粒子 ”
- 每个粒子有一个适应度函数，以模拟每只鸟与食物的丰盛程度
- 每个粒子有一个速度决定它的飞行方向和距离，初始值可以随机确定
- 每一次单位时间的飞行后，所有粒子分享信息，下一步将飞向自身最佳位置和全局最优位置的加权中心
1.3 基本步骤
初始化一群随机粒子，通过迭代找到最优。每次迭代中，粒子通过跟踪 “ 个体极值(pbest) ” 和 “ 全局极值(gbest) ” 来更新自己的速度和位置。
假设在D维搜索空间中，有m个粒子，各粒子速度和位置的更新如下：
（1）其中第i个粒子的位置为矢量：
（2）其飞行速度也是一个矢量， 记为：
（3）第i个粒子搜索到的最优位置为：
（4）整个粒子群搜索到的最优位置为
（5）第i个粒子的位置和速度更新为
2、代码分析
2.1 Rastrigrin()函数
数学表达式：
代码：
function y = Rastrigin(x)
% Rastrigin函数
% 输入x,给出相应的y值,在x = ( 0 , 0 ,…, 0 )处有全局极小点0.
[row,col] = size(x);
if row > 1 error( ' 输入的参数错误 ' );
end
y =sum(x.^2-10*cos(2*pi*x)+10);
%y =-y;
2.2 Schaffer()函数
数学表达式：
代码：
function y=Schaffer(x)
[row,col]=size(x);
if row>1error('输入的参数错误');
end
y1=x(1,1);
y2=x(1,2);
temp=y1^2+y2^2;
y=0.5+(sin(sqrt(temp))^2-0.5)/(1+0.001*temp)^2;
2.3 Griewank()函数
数学表达式：
代码：
function y=Griewank(x)
%Griewan函数
%输入x,给出相应的y值,在x=(0,0,…,0)处有全局极小点0.
[row,col]=size(x);
if row>1error('输入的参数错误');
end
y1=1/4000*sum(x.^2);
y2=1;
for h=1:coly2=y2*cos(x(h)/sqrt(h));
end
y=y1-y2+1;
%y=-y;
2.4 fun()函数
作用：根据label的值，选择不同的算法计算粒子适应度值。
代码：
function y = fun(x,label)
%函数用于计算粒子适应度值
%x input 输入粒子
%y output 粒子适应度值
if label==1y=Rastrigin(x);
elseif label==2y=Schaffer(x);
elsey= Griewank(x);
end
2.5 fun1()函数
作用：通过函数进行粒子适应度值计算。
代码：
function y = fun1(x)
%函数用于计算粒子适应度值
%x input 输入粒子
%y output 粒子适应度值
y=-20*exp(-0.2*sqrt((x(1)^2+x(2)^2)/2))-exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)+20+exp(1);
%y=x(1)^2-10*cos(2*pi*x(1))+10+x(2)^2-10*cos(2*pi*x(2))+10;
2.6 Drawfunc()函数
作用：根据label的值绘制不同的三维曲面网格图
代码：
function Drawfunc(label)
%该函数用于绘制图像
x=-5:0.05:5;%1行201列的向量if label==1y = x;[X,Y] = meshgrid(x,y);%meshgrid函数生成的X，Y是大小相等的矩阵，x，y是两个网格矢量，x，y都是行向量。%X：通过将x复制length(y)行（严格意义上是length(y)-1行）得到%Y：首先对y进行转置得到y'，将y'复制（length(x)-1）次得到[row,col] = size(X);for l = 1 :colfor h = 1 :rowz(h,l) = Rastrigrin([X(h,l),Y(h,l)]);endendsurf(X,Y,z);%surf()用于绘制比较光滑的三维曲面网格图shading interp%在flat的基础上进行色彩的插值处理，使色彩平滑过渡xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); title('Rastrigrin mesh');
endif label==2y = x;[X,Y] = meshgrid(x,y);[row,col] = size(X);for l = 1 :colfor h = 1 :rowz(h,l) = Schaffer([X(h,l),Y(h,l)]);endendsurf(X,Y,z);shading interp xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); title('Schaffer mesh');
endif label==3y = x;[X,Y] = meshgrid(x,y);[row,col] = size(X);for l = 1 :colfor h = 1 :rowz(h,l) = Griewank([X(h,l),Y(h,l)]);endendsurf(X,Y,z);shading interp xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); title('Griewank mesh');
end
2.7 PSO.m
%% 清空环境
clc
clear%% 参数初始化
%粒子群算法中的三个参数
c1 = 1.49445;%加速因子
c2 = 1.49445;
w=0.8; %惯性权重maxgen=1000; % 进化次s数
sizepop=200; %种群规模Vmax=1; %限制速度范围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数func=1; %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
subplot(1,2,1);
Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出%% 产生初始粒子和速度
for i=1:sizepop%随机产生一个种群pop(i,:)=popmax*rands(1,dim); %初始种群V(i,:)=Vmax*rands(1,dim); %初始化速度%计算适应度fitness(i)=fun(pop(i,:),func); %粒子的适应度
end%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:); %全局最佳
pbest=pop; %个体最佳
fitnesspbest=fitness; %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值count2=fitnessgbest;
count=1;%% 迭代寻优
for i=1:maxgenif fitnessgbest<count2count=i;count2=fitnessgbest;endfprintf('第%d代，',i);fprintf('最优适应度%f\n',fitnessgbest);for j=1:sizepop%速度更新V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大V(j,find(V(j,:)<Vmin))=Vmin;w = (maxgen-i)/maxgen*0.8+0.2; %w随迭代次数减少%种群更新pop(j,:)=pop(j,:)+0.5*V(j,:); %位置更新 0.5:可不乘pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围pop(j,find(pop(j,:)<popmin))=popmin;if rand>0.98 %加入变异种子，用于跳出局部最优值pop(j,:)=rands(1,dim); %丰富种群多样性，防止过早陷入局部最优end%更新第j个粒子的适应度值fitness(j)=fun(pop(j,:),func); endfor j=1:sizepop%个体最优更新if fitness(j) < fitnesspbest(j)pbest(j,:) = pop(j,:);fitnesspbest(j) = fitness(j);end%群体最优更新if fitness(j) < fitnessgbestgbest = pop(j,:);fitnessgbest = fitness(j);endend yy(i)=fitnessgbest; end
fprintf('收敛次数%d\n',count);%% 结果分析
subplot(1,2,2);
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
2.8 PSOMutation.m
%% 清空环境
clc
clear%% 参数初始化
%粒子群算法中的两个参数
c1 = 1.49445;
c2 = 1.49445;maxgen=500; % 进化次数
sizepop=100; %种群规模Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;%% 产生初始粒子和速度
for i=1:sizepop%随机产生一个种群pop(i,:)=5*rands(1,2); %初始种群V(i,:)=rands(1,2); %初始化速度%计算适应度fitness(i)=fun1(pop(i,:)); %染色体的适应度
end%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:); %全局最佳
gbest=pop; %个体最佳
fitnessgbest=fitness; %个体最佳适应度值
fitnesszbest=bestfitness; %全局最佳适应度值%% 迭代寻优
for i=1:maxgenfor j=1:sizepop%速度更新V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));V(j,find(V(j,:)>Vmax))=Vmax;V(j,find(V(j,:)<Vmin))=Vmin;%种群更新pop(j,:)=pop(j,:)+0.5*V(j,:);pop(j,find(pop(j,:)>popmax))=popmax;pop(j,find(pop(j,:)<popmin))=popmin;if rand>0.98 pop(j,:)=rands(1,2);end%适应度值fitness(j)=fun1(pop(j,:)); endfor j=1:sizepop%个体最优更新if fitness(j) < fitnessgbest(j)gbest(j,:) = pop(j,:);fitnessgbest(j) = fitness(j);end%群体最优更新if fitness(j) < fitnesszbestzbest = pop(j,:);fitnesszbest = fitness(j);endend yy(i)=fitnesszbest; end
%% 结果分析
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
2.9 运行结果
2.9.1 使用PSO.m
(1)
共迭代1000次，收敛时的迭代次数为第779次，最终其最优适应度为4.974795，与理论最小值0相差约为5。
(2)
左图为Rastrigin函数图像，右图体现了在迭代过程中最优适应度一开始变化剧烈，然后变化细微。（例如，从图中看到，最优适应度收敛于第105次迭代，但第(1)张图中，显示收敛时迭代次数为779。说明迭代次数从第105次开始到第779次之间，最优适应度的改变很小，低于小数点后6位。）
(3)
总运行时间为6.694s。
2.9.2 使用PSOMutation.m
适应度函数使用的是一个公式。
(1)
迭代到第41次时，最优适应度为0.002994，与理论最小0极为接近，且在之后的迭代过程中保持最优适应度为0.002571。
(2)
总运行时间为0.969s。
3、参数调整（针对POS.m文件）
由于代码中涉及到rand()函数，导致相同情况下的运行结果不同。因此，采用多次（10次）运行结果记录平均值的方法。
参数说明：
| 符号 | 含义 |
|---|---|
| c1 | 自身最佳位置的权重 |
| c2 | 全局最佳位置的权重 |
| w | 惯性权重 |
| sizepop | 种群规模 |
| dim | 维度 |
| func | 适应度函数：1-Rastrigrin()函数；2-Schaffer()函数；3-Griewank()函数 |
3.1 w改变，其他参数不变
(c1=1.49445,c2=1.49445,sizepop=200,dim=10)
3.1.1 w初始值为0.8
w的表达式不同，运行效果也不同。
（1）func=1时：
| w表达式 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 0.9*w+0.1 | 3.331237 | 934 | 5.0042s |
| (maxgen-i)/maxgen*0.8+0.2 | 2.074113 | 920 | 5.204s |
| abs(w-rand()) | 0.1989918 | 879 | 5.0514s |
适应度函数为Rastrigrin()函数时，第三个表达式体现出来的效果最好。最优适应度都较小，但第三个表达式最接近理论最小值，收敛时迭代的次数三者都很大，总运行时间也都很接近。
（2）func=2时
| w表达式 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 0.9*w+0.1 | 0 | 414 | 4.9882s |
| (maxgen-i)/maxgen*0.8+0.2 | 0 | 335 | 4.9164s |
| abs(w-rand()) | 0 | 88 | 5.0046s |
适应度函数为Schaffer()函数时，最优适应度均达到了理论最小值，总运行时间差距不大，在0.1s之内。而收敛时迭代次数则出现明显差别，在三者中，第三个表达式达到收敛时所用迭代次数最少。
（3）func=3时
| w表达式 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 0.9*w+0.1 | 0.0000314 | 734 | 5.3902s |
| (maxgen-i)/maxgen*0.8+0.2 | 0 | 498 | 5.2616s |
| abs(w-rand()) | 0 | 421 | 5.4036s |
适应度函数为Griewank()函数时，第一个表达式的最优适应度非常接近0，收敛时迭代次数很大；而其他两个表达式最优适应度都为0，收敛时迭代次数相对接近；三者的总运行时间差距不大。
总结：不同的表达式对于PSO算法取得全局最优具有不同的效果，其中，w=abs(w-rand())的效果较好。
3.1.2 w=abs(w-rand())
（1）func=1时：
| w初始值 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 1 | 0.4892758 | 925 | 5.1834s |
| 0.8 | 0.9954 | 862.8 | 5.1578s |
| 0.6 | 0 | 716 | 5.2448s |
适应度函数为Rastrigrin()函数时，惯性权重w的初始值取越小，仿佛越有优势。但需注意的是，w在一开始时应较大，慢慢的减小，因此，初始值不应太低。
（2）func=2时
| w初始值 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 1 | 0 | 41 | 5.0072s |
| 0.8 | 0 | 39 | 5.0496s |
| 0.6 | 0 | 40 | 5.0448s |
适应度函数为Schaffer()函数，惯性权重w的初始值处于[0.6，1]之间时，最优适应度、收敛时迭代次数和总运行时间差距很小，可看作相等。
（3）func=3时
| w初始值 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 1 | 0 | 326 | 5.4196s |
| 0.8 | 0 | 432 | 5.4236s |
| 0.6 | 0 | 323 | 5.4344s |
适应度函数为Griewank()函数，惯性权重w的初始值处于[0.6，1]之间时，最优适应度、收敛时迭代次数和总运行时间差距很小，可看作相等，但收敛时迭代次数略微过大。
总结：同个表达式具有不同的初始值时，对PSO算法取得正确的最优适应度有不同的效果。
3.2 c1、c2改变，其他参数不变
(初始值c1=2,c2=1,w=0.8,sizepop=200,dim=10,w=(maxgen-i)/maxgen*0.8+0.2)
（1）func=1时
| c1 | c2 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|---|
| c1=c1-rand()*0.8 | c2+rand()*0.8 | 2.308363 | 967 | 5.2028s |
| 0.8*c1 | 1.2*c2 | 1.591934 | 702 | 5.1156s |
| 2/c1 | (rand()+1)*c2 | 2.785971 | 883 | 5.2578s |
适应度函数为Rastrigrin()函数时，c1的三个表达式中第二个表达式比其余两个效果好，但收敛时的迭代次数很大。
（2）func=2时
| c1 | c2 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|---|
| c1=c1-rand()*0.8 | c2+rand()*0.8 | 0 | 312 | 5.1262s |
| 0.8*c1 | 1.2*c2 | 0 | 331 | 5.0706s |
| 2/c1 | (rand()+1)*c2 | 0 | 337 | 4.801s |
适应度函数为Schaffer()函数时，c1的三个表达式所起的效果是近似的，并没有很大的区别，收敛时迭代次数居中。
（3）func=3时
| c1 | c2 | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|---|
| c1=c1-rand()*0.8 | c2+rand()*0.8 | 0 | 569 | 5.5106s |
| 0.8*c1 | 1.2*c2 | 0 | 530 | 5.4486s |
| 2/c1 | (rand()+1)*c2 | 0 | 600 | 5.5152s |
适应度函数为Griewank()函数时，c1的三个表达式所起的效果是近似的，并没有很大的区别，收敛时迭代次数偏大。
总结：采用适当的c1、c2的表达式，对于PSO算法有着积极的作用。
部分代码修改：
<1>
c1=c1-rand()*0.8;while c1 <= 1 c1=rand()*0.2+1;endc2=c2+rand()*0.8;while c2 >=2 c2=rand()*0.2+1;end
<2>
c1=0.8*c1;while c1 <= 1 c1=c1+1;endc2=1.2*c2;while c2 >=2 c2=0.8*c2;end
<3>
c1=2/c1;while c1 <= 1 c1=c1+1;endc2=(rand()+1)*c2;while c2 >=2 c2=0.8*c2;end
3.3 sizepop改变，其他参数不变
(初始值c1=2,c2=1,w=0.8,dim=10,c1=0.8c1, c2=1.2c2,w=(maxgen-i)/maxgen*0.8+0.2)
（1）func=1时
| sizepop | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 200 | 2.586893 | 945 | 5.278s |
| 500 | 1.392943 | 824 | 11.8774s |
| 1000 | 0.4001118 | 893 | 41.0784s |
| 2000 | 0 | 668 | 47.4824s |
适应度函数为Rastrigrin()函数时，随着种群规模sizepop的增加，最优适应度递减至0，收敛时迭代次数均大于600次，总运行时间递增。
（2）func=2时
| sizepop | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 200 | 0 | 324 | 5.124s |
| 500 | 0 | 302 | 11.9126s |
| 1000 | 0 | 297 | 24.0288s |
| 2000 | 0 | 286 | 46.854s |
适应度函数为Schaffer()函数时，最优适应度均为0，随着种群规模sizepop的增加，收敛时迭代次数变化递减，总运行时间递增。
（3）func=3时
| sizepop | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 200 | 0 | 557 | 5.4996s |
| 500 | 0 | 426 | 12.5244s |
| 1000 | 0 | 375 | 25.6348s |
| 2000 | 0 | 345 | 49.5876s |
适应度函数为Griewank()函数时，最优适应度均为0，随着种群规模sizepop的增加，收敛时迭代次数递减，总运行时间递增。
总结：随着种群规模的增加，原来取不到全局最优的适应度函数也会较大机会取得；在最优适应度改变不大的同时，加快了收敛速度。
3.4 dim改变，其他参数不变
(初始值c1=2,c2=1,w=0.8,sizepop=200,c1=0.8c1, c2=1.2c2,w=(maxgen-i)/maxgen*0.8+0.2)
（1）func=1时
| dim | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 10 | 3.781157 | 987.4 | 5.784s |
| 15 | 7.364999 | 992 | 6.1s |
| 20 | 4.975064 | 901 | 5.7082s |
| 25 | 10.347574 | 970 | 6.0452s |
适应度函数为Rastrigrin()函数时，随着空间维度dim的递增，最优适应度越来越大，即越来越偏离理论最小值；收敛时迭代次数很大，总运行时间处于[5,7]s内。
（2）func=2时
| dim | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 10 | 0 | 328 | 5.46s |
| 15 | 0 | 326 | 5.505s |
| 20 | 0 | 328 | 5.6902s |
| 25 | 0 | 308 | 5.725s |
适应度函数为Schaffer()函数时，随着空间维度dim的递增，最优适应度并没有改变，一直为0；收敛时迭代次数几乎相等，总运行时间差距不大。
（3）func=3时
| dim | 最优适应度 | 收敛时迭代次数 | 总运行时间 |
|---|---|---|---|
| 10 | 0 | 568 | 5.8846s |
| 15 | 0 | 687 | 6.0174s |
| 20 | 0 | 871 | 6.1354s |
| 25 | 0 | 921 | 6.423s |
适应度函数为Griewank()函数时，随着空间维度dim的递增，最优适应度并没有改变，一直为0；收敛时迭代次数递增，总运行时间差距在1s之内。
总结：
（1）适应度函数为Rastrigrin()函数时，随着空间维度的增加，最优适应度也越偏离理论最小值0，而Schaffer()函数和Griewank()函数，最优适应度固定不变，为0。
（2）适应度函数为Rastrigrin()函数和Schaffer()函数时，收敛时迭代次数在某范围徘徊；适应度函数为Griewank()函数时，收敛时迭代次数随着空间维度的增加也递增。
（3）这三者适应度函数，随着空间维度的增加，总运行时间并没有很明显的变化。