粒子群算法介绍。粒子动态演示代码

粒子群算法介绍04：附Matlab、python代码（PSO）

九

这次在上述代码中动态演示粒子运动变化

自己运行代码看动图。通过动态演示你更能发现粒子数量n 迭代次数maxgen，惯性权重w的大小对最优解的影响。如上图就是迭代次数太少加上粒子数量太少还有惯性权重太小，从而陷入局部最优解。

%% 双清空+关闭图形
clc,clear,close all;

%% 绘制函数图像
[x y]=meshgrid(-2:0.02:2,-2:0.05:2);
z=sin( sqrt(x.^2+y.^2) )./sqrt(x.^2+y.^2)+exp((cos(2*pi*x)+cos(2*pi*y))/2)-2.71289;
figure(1)
surf(x,y,z)
hold on %不关闭图形，继续画图

%% PSO参数标定。
maxgen = 300;     %迭代次数
n = 50;           %种群大小
narvs = 2;        %维度
w = 0.9;          %惯性权重
c1 = 1.49445;     %个体学习因子
c2 = 1.49445;     %社会学习因子
Vmax = 0.5;       %速度上下限
Vmin = -0.5;
Xmax = 2;         %位置上下限
Xmin = -2;

%% 初始化粒子群中每个粒子的位置和速度，并计算适应度值
x = 2*rands(n,narvs);    %随机初始化粒子所在的位置
v = 0.5*rands(n,narvs);  %随机初始化粒子的速度

%% 寻找初始个体极值和群体极值  
fit = zeros(n,1);           %初始化这n个粒子的适应度全为0
for i = 1:n                 %循环整个粒子群，计算每一个粒子的适应度
    fit(i) = fun1(x(i,:));  %调用fun1函数来计算适应度
end
pbest = x;                      %个体极值
pbestfitness = fit;             %个体极值适应度值
[gbestfitness index]=max(fit);  %找到群体极值适应度值
gbest = x(index,:);             %群体极值

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,750,'.r');  %scatter3是绘制三维散点图的函数

%% 寻优迭代
for i = 1:maxgen
    %速度更新
    for j = 1:n
        v(j,:)=w*v(j,:)+c1*rand*(pbest(j,:)-x(j,:))+c2*rand*(gbest-x(j,:)); %速度更新公式
        v(j,find(v(j,:)>Vmax)) = Vmax; %find函数这里返回的是索引可以是三种情况1；2；1 2；
        v(j,find(v(j,:)<Vmin)) = Vmin;
     %位置更新   
        x(j,:)=x(j,:)+v(j,:);                                               %位置更新公式   
        x(j,find(x(j,:)>Xmax)) = Xmax;
        x(j,find(x(j,:)<Xmin)) = Xmin;
      %更新适应度值  
        fit(j)=fun1(x(j,:));
    end
      %个体极值更新
    for j = 1:n
        if fit(j)>pbestfitness(j)
            pbest(j,:)=x(j,:);
            pbestfitness(j)=fit(j);
        end
      %群体极值更新 
        if fit(j)>gbestfitness
            gbest=x(j,:);
            gbestfitness=fit(j);
        end
    end
    %每代最优极值记录到result中
    result(i)=gbestfitness;
    pause(0.2);          %暂停0.1s
    h.XData = x(:,1);    %更新散点图句柄的x轴的数据
    h.YData = x(:,2);    %更新散点图句柄的y轴的数据
    h.ZData = fit;       %更新散点图句柄的z轴的数据
end

%% 绘制图像
figure(2)
plot(result)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

gbest
gbestfitness 
       

在matlab中运行上述代码，适应度函数如下：

function y = fun1(x)
%函数用于计算粒子适应度值
%x           input           输入粒子 
%y           output          粒子适应度值 
y=sin( sqrt(x(1).^2+x(2).^2) )./sqrt(x(1).^2+x(2).^2)+exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)-2.71289;