粒子寻优算法PSO学习笔记

PSO-MATLAB实现

1. 粒子寻优算法的概述（PSO， particle swarm optimization）

PSO算法是由Eberhart博士和Kennedy博士基于对鸟类捕食问题研究，所提出的一种群体寻优算法。
一个对PSO算法简单理解的例子：一群鸟进行觅食，而远处有一片玉米地，所有的鸟都不知道玉米地到底在哪里，但是它们知道自己当前的位置距离玉米地有多远。那么找到玉米地的最佳策略，也是最简单有效的策略就是搜寻目前距离玉米地最近的鸟群的周围区域。在PSO中，每个优化问题的解都是搜索空间中的一只鸟，称之为"粒子"，而问题的最优解就对应于鸟群中寻找的"玉米地"。
而PSO算法的实际就是将上述问题转化为数学语言描述的方程。因此，假设
所有的粒子（N个粒子），在N维空间内都具有一个位置向量（P）（粒子在解空间的位置）和速度向量（V）（决定下次飞行的方向和速度）。并用由粒子位置信息组成的目标函数去表示粒子到目标的距离（即粒子的适应值，fitness value），可以将其理解为鸟距离"玉米地"的距离。
在上诉的基础上，PSO运用下面的数学公式对粒子速度和位置信息进行迭代。其迭代的目的是找出目标函数的最佳适应度值。

对公式的理解。
速度迭代公式：第一部分为粒子本身速度的惯性wV（t）+第二部分为粒子自身学习认为的最佳适应度+第三部分为种群中的最优适应度值得学习。
由上述两个公式，从而确定下一次迭代时需要如何调整和改变飞行的方向和速度。就这样逐步迭代，最终整个种群的例子就会逐步趋于最优解。

2.MATLAB代码得实现

程序学习于：
1: https://blog.csdn.net/wendy3orz/article/details/85345762
2: https://blog.csdn.net/weixin_45389639/article/details/103370182
3: https://ww2.mathworks.cn/help/matlab/ref/surf.html
主程序mypso.m

%% 清空环境
clc
clear


%% 预设值参数
%学习因子c1 c2范围为1-4，一般设置为2左右
c1=2;  %自我学习因子
c2=2;  %社会学习因子
%惯性权重w 惯性权重 w 较大时，算法的收敛速度快，全局搜索能力强；w过小时，算法的收敛速度慢，局部搜索能力强
w=0.8;         %自身速度惯性
w2=0.4;
%迭代次数 标准 PSO 的最大迭代次数一般为 100~200，对于更加复杂的优化问题，可适当增加最大迭代次数
iteration=100; 
%种群规模
size_N=100;

Vmax=1;       %限制速度围
Vmin=-1;     
popmax=20;     %变量取值范围
popmin=-20;
dim=2;       %适应度函数维数
index=1;%用于选取目标函数

%% 初始化种群 产生初始粒子和速度
%% 产生初始粒子和速度
for i=1:size_N
    %随机产生一个种群(使用rands分布对种群中的每一行产生一行随机数，然后分布计算每一行的适应度值)
    pop(i,:)=popmax*rands(1,dim);    %初始种群  (使用rands函数随机产生从0-1的函数)并使用自定义的popmax提升粒子位置范围
    V(i,:)=Vmax*rands(1,dim);             %初始化速度
                                     %计算适应度
    fitness(i)= fitness_function(pop(i,:),index);   %计算每个粒子的适应度
end

%% 迭代前的个体极值和群体极值
[bestfitness bestindex]=min(fitness);  %使用min函数寻找fitness中的最小值的行和列 bestindex=取到最小fitness的行
gbest=pop(bestindex,:);   %全局最佳 最佳那一行pop
pbest=pop;                %个体最佳
fitnesspbest=fitness;     %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值 

%% 迭代寻优
for i=1:iteration
    fprintf('第%d代，',i);          %输出第几代
    fprintf('最优适应度%f\n',fitnessgbest);  %输出最优适用度
    w=w+(w2-w)/iteration;
    for j=1:size_N
        %速度更新
        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;
        
        %种群更新
        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,dim);
        end
        
        %更新第j个粒子的适应度值
        fitness(j)=fitness_function(pop(j,:)); 
   
    end
    
    for j=1:size_N
        
        %个体最优更新
        if fitness(j) < fitnesspbest(j)
            pbest(j,:) = pop(j,:);
            fitnesspbest(j) = fitness(j);
        end
        
        %群体最优更新
        if fitness(j) < fitnessgbest
            gbest = pop(j,:);
            fitnessgbest = fitness(j);
        end
    end 
    yy(i)=fitnessgbest;  
    
    % 收敛动图绘制存储
    plot(pop(:,1),pop(:,2),'*b')
    axis([popmin popmax popmin popmax])
    pause(0.1)      %设置站台的秒数，可以辅助看粒子群的群体优化效果
    x1=xlabel('x1');
    x2=ylabel('x2');
    title(['进化次数=' num2str(i)]);
    drawnow;
    frame = getframe(1);
    im = frame2im(frame);
    [A,map] = rgb2ind(im,256);
        
end

%% 结果分析
my_drawfun(index);
figure;
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

画图函数my_drawfun.m

function my_drawfun(index)
%index用于设置不同的目标函数，拓展所写的pso函数

x=-5:0.05:5;%41列的向量 %从-5到间隔0.05到达5
y = x;
[X,Y] = meshgrid(x,y);  %返回一个X Y 均为41行和列的矩阵
[row,col] = size(X);
for  l = 1 :col
    for  h = 1 :row
        z(h,l) = (X(h,l))^2+(Y(h,l))-7*cos(pi*(X(h,l)))-8*cos(2*pi*(Y(h,l)));
    end
end
surf(X,Y,z,'FaceAlpha',0.5);       %三维函数曲面画图 'FaceAlpha',0.5半透明的效果
shading interp    %对曲面或图形对象的颜色着色进行色彩的插值处理，使色彩平滑过渡
%colormap(gca, 'purla')   %jet  purla
xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); 
title('mesh'); 

目标函数fitness_function.m

function y=fitness_function(x,index)
%index用于设置不同的目标函数，拓展所写的pso函数
[row,col]=size(x);
if row>1
    error('输入的参数错误');
end
y1=x(1,1);
y2=x(1,1);
temp=y1^2+y2^2;
y=y1^2+y2-7*cos(pi*y1)-8*cos(2*pi*y2);
%y=sum(x.^2);

3.程序结果

（1）随着迭代的进行，种群由分散的粒子逐渐向最优目标靠近

4. PSO程序参数的学习

在运用pso算法解决实际问题时，只需要确定目标函数。需要调节的参数有w惯性权重，c1,c2加速因子，种群规模和迭代次数和粒子的维度。

4. 1 种群规模与迭代次数对程序寻优的影响

（1）N=100;iteration=100;程序运行时间：0.22；最优适应度-15.001292
（2）N=100;iteration=200;程序运行时间：0.32；最优适应度-15.001292
（3）N=100;iteration=300;程序运行时间：0.43；最优适应度-15.001292

（4）N=100;iteration=200;程序运行时间：0.32；最优适应度-15.001292
（5）N=200;iteration=200;程序运行时间：0.50；最优适应度-15.001292
由以上的对比，种群规模和种群数量只要在一个适合的范围，PSO算法会搜索出最佳值，继续增加种群规模和种群数量并不一定会增加精准度。所以在选取迭代次数和种群规模时，应该根据模型的复杂程度去设定，一般标准的PSO算法只要跌倒次数到达100-200就能达到效果。

4. 2惯性权重对程序寻优的影响

结论：经过研究，在粒子算法初期，选择较大的惯性权重值，可以使得算法有很强得全局搜索能力，而在粒子算法后期选择较小得权重可以使粒子逐渐收敛到全局最优。较高的w设置促进全局搜索，较低的w设置促进快速的局部搜索。设置惯性权重最大值越小，平均最优适应度越低，实验结果越优。所以本实验设置最大权重0.6，最小权重为0.1
学习原文的链接：https://blog.csdn.net/weixin_45389639/article/details/103370182
w为惯性因子，是速度进行迭代时的保留上一次速度的权重值。

4. 3加速因子c1、c2对程序寻优的影响