粒子群算法介绍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
本人代码与文献代码略有出入——改正了文献代码错误之处。