机器学习——粒子群算法（PSO）及程序的实现

rustum

粒子群算法是一种全局优化算法，是机器学习中很重要的一种算法。它的灵感来源于对鸟群捕食行为的研究。
一、算法介绍

假想一群鸟儿迁徙到一块新的区域，它们需要寻找一个食物最多的地方供它们集体生存。我们假设，一个地方的食物多少仅仅和这个地方的地理位置有关，即食物量是地理位置的确定的函数
但是这么大一片区域怎么寻找最优解？
有鸟提议：
（1）大家分散开寻找，这样效率会比较高。
（2）寻找的时候不能乱飞。每只鸟飞翔过程中，需要每隔一段时间更新一次自己搜索历史上找到的食物最多的那个位置，然后根据这个位置修改自己的飞行方向。
（3）每只鸟都会朝着自己心目中食物最多的位置前进，但是这位置并不见得是整个区域食物最多的位置。鸟群需要信息交流，大家需要把自己找到的最优位置比较，确定整个集体历史上找到的最优的位置，然后大家都根据这个位置修改各自的飞行方向。
于是对于一只鸟（编号为i）来说，假设它在第k次更新时，它的位置是 ,它飞行的速度是 ，设时间间隔为1，那么下一次根据运动学公式，它会运动到：


而在每一次更新时，这只鸟会同时面对两个信息：
（1）自己历史上找到的最优的那个位置pbest
（2）整个鸟群历史上找到的最优的那个位置gbest
于是这只鸟就将自己的速度做出修正：在自己的速度矢量的基础上，再分别叠加指向pbest的矢量和指向gbest的矢量（如figure1所示），于是到下一次更新时，这只鸟的速度变成了：


其中第一项表示鸟的惯性：它没办法瞬间停下来，要保持它原有的速度继续运动一段距离。c1和c2是正数，可以人为设定，有时为了造成一定随机性，会在c1和c2系数上乘一个0到1的随机数。



figure1

这个过程就是粒子群算法的基本过程。
（1）给定一个具体的函数，目标是寻找它的最大值\最小值
（2）在空间中随机生成一些粒子
（3）让粒子基于两条公式演化：




（4）直至迭代到粒子群收敛于一个位置点，那个位置就是我们要找的全局最优位置。

二、一元函数的极值问题matlab实践

我们先来看一元函数。
对于一个不是特别单调的一元函数（如 )，我们要求它最大值的位置。
1.参数设置与函数编写

c1=0.005;
c2=0.005;
n=10;
s=30*pi;
T=200;这里s表示，我们稍后会在（-s,s)区域内随机生成粒子（这不代表粒子运动时不会离开这个区域，我的这个问题里没有考虑约束条件）。我们的plot图像稍后也会在这个范围内显示。
n是粒子数，c1、c2是上面速度更新公式的参数（设置太高，粒子会运动过快，不易收敛；设置太低，粒子运动太慢，收敛太迟，所以可以多次尝试，找到最好的组合），T是迭代次数（设置到粒子收敛就OK）
function [yy]=ff(xx)
  yy=sin(xx)./xx;
end这里编写一个函数ff，如果想改其他的函数玩可以直接在这个函数体里修改。
2.初始粒子随机生成

x_particle=s*(2*rand(1,n)-1);
v_particle=zeros(1,n);
pbest=x_particle;
u=find(ff(pbest)==max(ff(pbest)));
gbest=pbest(u);想象这些粒子只能在一维x轴上运动，我们在(-s,s)区域随机撒上n=10个粒子。
第一行x_particle是粒子的位置，2*rand(1,n)-1会随机产生n个在-1到1的随机数。
第二行v_particle初始全部设置为0，后续它们会运动起来。
第三行pbest是个1行n列的向量，储存每个粒子各自历史上找到的f(x)最大的位置。初始时每个粒子只存在过它们刚刚生存的位置，所以一开始pbest=x_particle。
第四行和第五行用来找gbest。gbest始终是一个数，它可以从pbest中寻找。u=find(ff(pbest)==max(ff(pbest)))返还ff(pbest)最大值的位置给u，然后gbest就是pbest索引u的值。
i=-s:0.05*pi:s;
plot(i,ff(i),'-b','linewidth',1.2)
xlim([-s,s]);
hold on
scatter(x_particle,ff(x_particle));
hold off初始画图，plot画函数,scatter画点。
3.迭代搜索

for t=1:1:T
    tempx=x_particle;
    tempv=v_particle;
    x_particle=tempx+tempv;
    v_particle=tempv+c1*rand()*(pbest-tempx)+c2*rand()*(gbest-tempx);
    for r=1:1:n
        if ff(pbest(r))>ff(x_particle(r))
            pbest(r)=pbest(r);
        else
            pbest(r)=x_particle(r);
        end
    end
    u=find(ff(pbest)==max(ff(pbest)));
    gbest=pbest(u);
    if mod(t,5)==0
        plot(i,ff(i),'-b','linewidth',1.2)
        xlim([-10*pi,10*pi]);
        hold on
        scatter(x_particle,ff(x_particle));
        hold off
        pause(1)
    end
end这里迭代T次，tempx和tempv储存上一次迭代的粒子位置和粒子速度，然后借助tempx和tempv更新x_particle和v_particle。
pbest的更新的方法是：将上一次的pbest和这次的x_particle逐个粒子比较（因为pbest就是历史最优，历史最优和现在最优比较就是加上当前的历史最优），对于第r个粒子，如果pbest的更好，那么保留pbest;如果x_particle的更好，那么保留x_particle。
gbest更新方法就不再多说了。
我这里设置的是每隔5步画一次图，画图的命令和之前一样。这里pause帮助我们看整个搜索的动态情况。如果您希望让这个过程更连续一些，可以缩小pause里面的值，或者减少画图的间隔步数。
整个代码我会放在后面的附录里。
三、粒子群算法效果与改进

经过第二模块的实践，我们来看一下效果，如figure 2所示。



figure 2

可见这些粒子并没有很好的收敛到最高点，它们没有完成它们的使命。
解决问题的一个方法是增加迭代时间T，给它们足够的时间让它们演化，但这样效果并不是非常的显著。
我们来思考一下根本原因在哪里。
假想一个粒子在距离很远处，它迟迟找不到满意的位置。这个时候总部发来消息说，在x=0处找到了最优解，于是这个粒子立即将自己的速度修正为c2(gbest-x)向x=0飞奔而去。
终于它来到了x=0处，发现这个地方确实比之前找到的地方函数值大多了，正当它准备就在这个点旁边和其它点收敛聚集时 ，却发现一个问题：由于“飞奔而来”速度太快，导致它不得不继续向前飞。如果每个粒子的情况都如此，那么收敛便成为了一件难事。
解决问题的方法是动态降低它的惯性。我们将速度更新公式修改为：


w(k)是一个关于时间变化的权值。我在这里设置它为：


我们来看看效果如何：在“迭代搜索”这一部分里，修改v_particle一行的代码为：
v_particle=(1-t/1000)*tempv+c1*rand()*(pbest-tempx)+c2*rand()*(gbest-tempx);最后效果如下：



figure 3

可见到了t=140，t=200时，粒子基本都收敛在x=0周围，效果很成功。
四、对粒子群算法的一些讨论

1.粒子群算法并不只能用在一元函数，理论上对于多元函数也是可行的。
2.粒子群算法是一种全局优化算法，但是它容易最后收敛于局部最优，因此要合理的设置参数，防止这种现象的产生。
3.粒子群算法生成粒子随机，迭代运动也有随机性，所以每一次运行的结果可能不尽相同。我们可以做多次粒子群算法的演化，计算它成功找到全局最优解的成功率。
五、粒子群算法的matlab工具

事实上matlab提供粒子群算法的优化工具。在matlab的工具箱global optimization toolbox中，有很多全局优化算法的函数工具，其中particleswarm函数是用来做粒子群算法的。具体可以详见https://ww2.mathworks.cn/help/gads/particleswarm.html
x = particleswarm(fun,nvars);fun表示需要优化的函数，nvars是fun的自变量个数，x输出fun取最小值时对应的自变量。
对于我们第二模块的例子，我们可以编一段代码：
result=particleswarm(@ff,1);
disp(result)

function [yy]=ff(xx)
  yy=-sin(xx)./xx;
end其中由于particleswarm取的是函数最小值，所以我们在function里面把函数取负号。
运行结果都几乎是在0的附近：

实验次数	1	2	3	4	5
result值	-4.1492e-6	2.5025e-5	1.3200e-8	-1.7738e-6	1.0344e-7

六、总代码

这里提供第二模块中的总代码：
%粒子群算法一维计算例
%fx=sinx/x 寻最大值
c1=0.005;
c2=0.005;
n=10;
s=30*pi;
T=200;

%随机分布产生粒子
x_particle=s*(2*rand(1,n)-1);
v_particle=zeros(1,n);
pbest=x_particle;
u=find(ff(pbest)==max(ff(pbest)));
gbest=pbest(u);
i=-s:0.05*pi:s;
plot(i,ff(i),'-b','linewidth',1.2)
xlim([-s,s]);
hold on
scatter(x_particle,ff(x_particle));
hold off
for t=1:1:T
    tempx=x_particle;
    tempv=v_particle;
    x_particle=tempx+tempv;
    v_particle=(1-t/1000)*tempv+c1*rand()*(pbest-tempx)+c2*rand()*(gbest-tempx);
    for r=1:1:n
        if ff(pbest(r))>ff(x_particle(r))
            pbest(r)=pbest(r);
        else
            pbest(r)=x_particle(r);
        end
    end
    u=find(ff(pbest)==max(ff(pbest)));
    gbest=pbest(u);
    if mod(t,5)==0
        plot(i,ff(i),'-b','linewidth',1.2)
        xlim([-10*pi,10*pi]);
        hold on
        scatter(x_particle,ff(x_particle));
        hold off
        pause(1)
    end

end



function [yy]=ff(xx)
  yy=sin(xx)./xx;
end