粒子群优化算法

1.算法简介

想象很多只鸟组成的一个鸟群，每只鸟拥有自己的位置和速度，每只鸟在捕食过程中对当前位置获得的食物有一个大致的估计并且可以同鸟群交流自己的信息，那么鸟群就会综合这些信息做出对某个方向的趋向运动，并最终稳定在一个大体位置

PSO算法是进化算法的一种，他是受鸟群捕食的启发，将每一个个体获得的信息同群体信息相结合来引导鸟群的运动方向

2.算法内容

鸟被抽象为 $n$ 个 $D$ 维空间的粒子，每个粒子有本身属性位置 $X_i={ \{ x_1,x_2,...,x_D\} }$ 以及速度 $V_i=\{ v_1,v_2,...,v_D\}$。除此之外，每个粒子有根据目标函数决定的适应值，自己的历史最好位置 $P_{best}$ 和群体历史最好位置 $G_{best}$ 。每个粒子就是根据自己和群体的经验来决定下一步的行动

公式：

\[\begin{align} &V_i=\omega \times V_i+c_1\times r_1\times(P_{besti}-X_i)+c_2\times r_2\times(G_{besti}-X_i) \tag1 \\ &X_i=X_i+V_i\tag2 \\ \end{align} \]

$i=1,2,...,n$ 表示粒子总数

$c_1,c_2$ 是学习因子，一般来说均为 $2$

$r_1,r_2$ 是介于 $0,1$ 之间的随机数

$\omega$ 为惯性因子，其值非负：值越大，全局寻优能力强；局部寻优能力弱，值越小则反之。

$\omega \times V_i$ 称为【记忆项】，表示上次速度大小和方向的影响。

$c_1\times r_1\times(P_{besti}-X_i)$ 称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分。

$c_2\times r_2\times(G_{besti}-X_i) $ 称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。

动态的惯性因子可以提供更好的效果，这里介绍线性递减权值策略。

\[\omega^{(t)}=(\omega_{ini}-\omega_{gnd})(G_k-g)/G_k+\omega_{gnd} \]

$G_k$ 表示最大迭代次数

$\omega_{ini}$ 表示初始惯性权值

$\omega_{gnd}$ 表示迭代至最大次数的惯性权值

典型权值：$\omega_{ini}=0.9,\omega_{gnd}=0.4$

迭代终止条件根据具体问题一般选为最大迭代次数 $G_k$ 或(和)微粒群迄今为止搜索到的最优位置满足预定最小适应阈值。

例：

求函数 $y=-x(x-2)$ 在区间 $[0,2]$ 之间的最大值

matlab代码：

n=2;c1=2;c2=2;vmax=0.1;vmin=-0.1;
x=[0,2];
xmax=-inf;xmin=inf;
v=[0.01,0.02];
for i=1:n
    y(i)=-1*x(i)*(x(i)-2);
    xmax=max(xmax,x(i));
    xmin=min(xmin,x(i));
end;
t=x(1):0.05:x(2);
plot(t,-t.*(t-2));
hold on;
gbest=y(1);
for i=1:n
    pbest(i)=y(i);
    gbest=max(gbest,pbest(i));
end;
MX=15;
wini=0.9;wgnd=0.4;
for i=1:MX
    w=(wini-wgnd)*(MX-i)/MX+wgnd;
    r1=rand;r2=rand;
    for j=1:n
        v(j)=w*v(j)+c1*r1*(pbest(j)-x(j))+c2*r2*(gbest-x(j));
        if v(j)>vmax
            v(j)=vmax;
        elseif v(j)<vmin
            v(j)=vmin;
        end;
        x(j)=x(j)+v(j);
        if x(j)>xmax
            x(j)=xmax;
        elseif x(j)<xmin
            x(j)=xmin;
        end;
        y(j)=-1*x(j)*(x(j)-2);
    end;
    for j=1:n
        pbest(j)=max(pbest(j),y(j));
        if(pbest(j)>gbest)
            gbest=pbest(j);
        end;
    end;
    for j=1:n
        fprintf("the point %d: x:%.6f v:%.6f\n",j,x(j),v(j));
        plot(x(j),y(j),'.r');
        hold on;
    end;
    fprintf("the value is %.6f\n\n",gbest);
end;

输出如下

the point 1: x:0.008667 v:0.008667
the point 2: x:1.900000 v:-0.100000
the value is 0.190000

the point 1: x:0.108667 v:0.100000
the point 2: x:1.800000 v:-0.100000
the value is 0.360000

the point 1: x:0.208667 v:0.100000
the point 2: x:1.700000 v:-0.100000
the value is 0.510000

the point 1: x:0.308667 v:0.100000
the point 2: x:1.600000 v:-0.100000
the value is 0.640000

the point 1: x:0.408667 v:0.100000
the point 2: x:1.500000 v:-0.100000
the value is 0.750000

the point 1: x:0.508667 v:0.100000
the point 2: x:1.400000 v:-0.100000
the value is 0.840000

the point 1: x:0.608667 v:0.100000
the point 2: x:1.300000 v:-0.100000
the value is 0.910000

the point 1: x:0.708667 v:0.100000
the point 2: x:1.200000 v:-0.100000
the value is 0.960000

the point 1: x:0.808667 v:0.100000
the point 2: x:1.100000 v:-0.100000
the value is 0.990000

the point 1: x:0.908667 v:0.100000
the point 2: x:1.000000 v:-0.100000
the value is 1.000000

the point 1: x:0.986896 v:0.078229
the point 2: x:0.946667 v:-0.053333
the value is 1.000000

the point 1: x:1.061017 v:0.074121
the point 2: x:1.046667 v:0.100000
the value is 1.000000

the point 1: x:0.966661 v:-0.094356
the point 2: x:0.994936 v:-0.051731
the value is 1.000000

the point 1: x:1.004409 v:0.037748
the point 2: x:0.984493 v:-0.010443
the value is 1.000000

the point 1: x:1.013804 v:0.009396
the point 2: x:1.000373 v:0.015880
the value is 1.000000

在图上标出这些点：