粒子群算法介绍。附matlab、python代码、适合小白

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

七：

matlab代码

1.求下面的非线性函数最大值

函数分析：
首先可以用matlab把该函数图形画出来；其次该函数维度是2，两个未知数。

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

%% x，y区间确定注意0.05，如果再小会出现报错，当然越小越精确，z为函数表达式
[x y]=meshgrid(-2:0.05: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;

%% 绘图
surf(x,y,z)

从函数图形可以看出，该函数有许多局部极大值点，通过旋转图形你可以看见图形中间空出来一部分

从而可以猜测极限位置为（0,0）当然x，y不能同时取0，只能无线接近，（0,0）附近取得极大值1.0054.
你可以调节代码中x，y间距将0.05变成0.005等等一直缩小（继续小下去会报错）

还有需要注意的是你可能出来的图黑色的，那你把图形放大就好了

下面就开始用PSO算法求该函数最大值：

函数代码：

function y = fun(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;

主程序代码：

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

%% PSO参数标定。 ***注意：这里少了惯性权重w和维度D***
c1 = 1.49445;     %个体学习因子
c2 = 1.49445;     %社会学习因子
maxgen = 300;     %迭代次数
sizepop = 20;     %种群大小
Vmax = 0.5;       %速度上下限
Vmin = -0.5;
popmax = 2;       %位置上下限
popmin = -2;

%% 初始化粒子群中每个粒子的位置和速度，并计算适应度值
for j = 1:sizepop               %自行参考for循环
    pop(j,:) = 2*rands(1,2);    %自行参考rands
    V(j,:) = 0.5*rands(1,2);
    
    fitness(j) = fun(pop(j,:)); %自行参考函数   
end

%% 寻找初始个体极值和群体极值  
%{
***注意：pop是一个20行2列矩阵其中第一列为X1，第二列为X2；V也是
fitness有看懂的大神请评论区解释一下，我理解的是一个20行1列的矩阵这样才能找到群体极值的位置，
但是程序给出的确是1行20列矩阵***
--------------------------------------------------------------------------------------------
这里需要注意许多matlab中的max参考博客没讲这一点
通过实验发现当矩阵A是一个行向量时，[maxnumber maxindex] = max(A)中maxindex返回的是最大值所在列
这样就解释了刚才的疑问！
A=[1,2,3,4,8,6,75,9,243,25]
A =
     1     2     3     4     8     6    75     9   243    25
[maxnumber maxindex] = max(A)
maxnumber =
     243
maxindex =
     9
--------------------------------------------------------------------------------------------
%}
[bestfitness bestindex] = max(fitness);  %自行参考max函数用法
zbest = pop(bestindex,:);                %群体极值位置       
gbest = pop;                             %个体极值位置  因为初始化所以每个粒子个体极值位置就是它本身随机的位置
fitnessgbest = fitness;                  %个体极值适应度值
fitnesszbest = bestfitness;              %群体极值适应度值



%% 寻优迭代
for i = 1:maxgen
    %速度更新
    for 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; %find函数这里返回的是索引可以是三种情况1；2；1 2；
        V(j,find(V(j,:)<Vmin)) = Vmin;
     %位置更新   
        pop(j,:)=pop(j,:)+V(j,:);
        pop(j,find(pop(j,:)>popmax)) = popmax;
        pop(j,find(pop(j,:)<popmin)) = popmin;
      %更新适应度值  
        fitness(j)=fun(pop(j,:));
    end
      %个体极值更新
    for j = 1:sizepop
        if fitness(j)>fitnessgbest(j)
            gbest(j,:)=pop(j,:);
            fitnessgbest(j)=fitness(j);
        end
      %群体极值更新 
        if fitness(j)>fitnesszbest
            zbest=pop(j,:);
            fitnesszbest=fitness(j);
        end
    end
    %每代最优极值记录到result中
    result(i)=fitnesszbest;
end

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

结合matlab工作区更能理解编写思路。
还需要注意PSO为智能优化算法，这里最后虽然是1.0054即函数最大值，坐标（0.000681015851278943，0.000952485377217871），但是真实情况是不是呢？并不知道因为永远取不到那个点位，只能无限接近。并且你在多次运行代码的时候还可能发现最大值为1.0053，这可能陷入局部最优解，而且迭代曲线每次都不一样。


参考文献来自Matlab智能算法30个案例分析（第二版）P131

本人代码与文献代码略有出入——改正了文献代码错误之处。

源代码链接：
(https://download.csdn.net/download/weixin_42212924/18826374)
创作不易，请勿白嫖！当然上面的代码复制粘贴到matlab中就好了。
PSO往期文章：
(https://blog.csdn.net/weixin_42212924/article/details/116808786)