基于matlab的简单易懂的粒子群算法

粒子群算法的特点就是原理简单，以及经过了很多年的验证，对于解决不少问题，都有较好的效果。本文章主打的就是用尽量最简单的代码解决这个问题。

说他原理简单，是因为他自始至终只涉及到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++代码将在下一篇完成。