粒子群算法的特点就是原理简单,以及经过了很多年的验证,对于解决不少问题,都有较好的效果。本文章主打的就是用尽量最简单的代码解决这个问题。
说他原理简单,是因为他自始至终只涉及到2个简单的公式,如下:
第一个是粒子的速度更新公式
其中,和是学习因子,和是[0,1]之间的随机数。
第二个是粒子的位置更新公式
分析这个公式,可以理解为一群鸟在找食物的位置。每只鸟都依据两个参考来确定自己下一步前进的方向和速度大小,一个是自己目前找过的所有点中最好的位置,即,另一个是所有的鸟目前找过的所有点中最好的位置,即。
以Rosebrock(罗森布罗克)函数为例,求解该函数最大值,函数如下:
Matlab代码如下:
tic
%粒子群算法
%本程序可以计算函数的最大值和最小值,只要对程序进行简单的几个符号修改即可
%本程序可以改变函数的自变量个数,只要对下面的参数n进行控制即可
clc
clear
m=200; %群体大小
g=10; %终止进化代数
g0=1;
n=2; %函数的变量个数
low=[-2.048 -2.048]; %自变量取值下限
upper=[2.048 2.048]; %自变量取值上限
y1=rand(m,n).*(upper-low)+low; %初始化粒子位置
y2=rand(m,n).*(upper-low)+low;
f3=zeros(1,m);
v=zeros(m,n); %粒子初始速度暂且设置为0
c1=1.5;
c2=1.5;
%下面三行是加惯性权值的改进粒子群,有利于提高结果的精度
wmax=0.9;
wmin=0.4;
dw=wmax-wmin;
for i=1:m
f1=hanshu(y1(i,:)); %各个粒子当前位置对应的函数值
f2=hanshu(y2(i,:));
if f1>f2 %符号更改处
f2=f1;
y2(i,:)=y1(i,:); %用数组y2来记住每个粒子当前搜索过的最好位置
end
f3(i)=f2; %用数组f3来记住每个粒子当前搜索过的最好位置对应的函数值
end
f=max(f3); %f3中的最大值就是当前全局最大值
y3=find(f==f3);
a3=y2(y3(1),:); %当前全局最大值对应的解
while g0<=g
for i=1:m %更新位置
w=wmax-dw*g0/g;
v(i,:)=w*v(i,:)+c1*rand()*(y2(i,:)-y1(i,:))+c2*rand()*(a3-y1(i,:));
y1(i,:)=y1(i,:)+v(i,:);
for j=1:n %进行边界处理,超过取值范围的取边界值
b=find(y1(:,j)>upper(j));
y1(b,j)=upper(j);
b=find(y1(:,j)<low(j));
y1(b,j)=low(j);
end
f1=hanshu(y1(i,:));
f2=hanshu(y2(i,:));
if f1>f2 %符号更改处
f2=f1;
y2(i,:)=y1(i,:);
end
if f2>f %符号更改处
f=f2;
a3=y2(i,:);
end
end
g0=g0+1;
end
disp(['最大值:',num2str(f)])
disp(['对应解:',num2str(a3)])
toc
function result2=hanshu(x)
result2=100*(x(1)^2-x(2))^2+(1-x(1))^2; %罗森布罗克函数,求最大值
end
运行结果如下:
最大值:3905.9262
对应解:-2.048 -2.048
时间已过 0.030835 秒。
也可以求解最小值问题,并可以更改自变量个数以及各个自变量对应的取值范围。
1.自变量个数直接更改n
2.对应取值范围更改low和upper数组
3.求最小值的话,把代码中注释“符号更改处”(共3处)当行的大于号换成小于号即可
4.hanshu里面的函数对应更换
示例如下,求解函数最小值
代码如下:
tic
%粒子群算法
%本程序可以计算函数的最大值和最小值,只要对程序进行简单的几个符号修改即可
%本程序可以改变函数的自变量个数,只要对下面的参数n进行控制即可
clc
clear
m=200; %群体大小
g=10; %终止进化代数
g0=1;
n=3; %函数的变量个数
low=[-2 -3 -4]; %自变量取值下限
upper=[2 3 4]; %自变量取值上限
y1=rand(m,n).*(upper-low)+low; %初始化粒子位置
y2=rand(m,n).*(upper-low)+low;
f3=zeros(1,m);
v=zeros(m,n); %粒子初始速度暂且设置为0
c1=1.5;
c2=1.5;
%下面三行是加惯性权值的改进粒子群,有利于提高结果的精度
wmax=0.9;
wmin=0.4;
dw=wmax-wmin;
for i=1:m
f1=hanshu(y1(i,:)); %各个粒子当前位置对应的函数值
f2=hanshu(y2(i,:));
if f1<f2 %符号更改处
f2=f1;
y2(i,:)=y1(i,:); %用数组y2来记住每个粒子当前搜索过的最好位置
end
f3(i)=f2; %用数组f3来记住每个粒子当前搜索过的最好位置对应的函数值
end
f=max(f3); %f3中的最大值就是当前全局最大值
y3=find(f==f3);
a3=y2(y3(1),:); %当前全局最大值对应的解
while g0<=g
for i=1:m %更新位置
w=wmax-dw*g0/g;
v(i,:)=w*v(i,:)+c1*rand()*(y2(i,:)-y1(i,:))+c2*rand()*(a3-y1(i,:));
y1(i,:)=y1(i,:)+v(i,:);
for j=1:n %进行边界处理,超过取值范围的取边界值
b=find(y1(:,j)>upper(j));
y1(b,j)=upper(j);
b=find(y1(:,j)<low(j));
y1(b,j)=low(j);
end
f1=hanshu(y1(i,:));
f2=hanshu(y2(i,:));
if f1<f2 %符号更改处
f2=f1;
y2(i,:)=y1(i,:);
end
if f2<f %符号更改处
f=f2;
a3=y2(i,:);
end
end
g0=g0+1;
end
disp(['最大值:',num2str(f)])
disp(['对应解:',num2str(a3)])
toc
function result2=hanshu(x)
result2=(x(1)-1)^2+(x(2)-2)^2+(x(3)-3)^2+4;%求最小值的函数
end
运行结果如下:
最大值:4.0009
对应解:1.0098 2.0273 2.9965
时间已过 0.033085 秒。
作者水平有限,存在的问题敬请读者指出,也欢迎大家交流。相应c++代码将在下一篇完成。