粒子群算法（PSO）及MATLAB代码实现（一）

叠个BUFF： 仅为自己学习之用，如有错误，请批评指正！

文章参考：
通俗易懂讲算法-最优化之粒子群优化（PSO）-bilibili-青青草原灰太郎
温正.MATLAB智能算法(第2版).清华大学出版社（第八章 粒子群算法）

1.算法背景

粒子群算法，也称粒子群优化算法或鸟群觅食算法（Particle Swarm Optimization），缩写为PSO 。该算法最初是受到飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模
型。
假设有一个鸟群（有小鸟A,B,C,……）随机降落在一片地上寻找虫子，其中有一个虫子是最大的。怎么找到这个最大的虫子？
小鸟们先自己寻找，找了一会儿，大家一比较，发现小鸟A是所有小鸟里面找到虫子是比较大的。其他小鸟怎么办?这个虫子是最大的吗？小鸟B想我找的虫子也很大啊，于是小鸟B就沿着自己的寻找方向和BA方向中间某一方向去找找看，然后鸟群再比较一次。如此往复！

2.算法基本原理与算法流程

问：已知A为全局最优，B和C如何移动才能到达A处？这个过程用数学表达式如何表达？
某个例子（点）的移动，是有大小，有方向的，也就是向量！位置就是坐标！
（1,1）=（2,3）+α 则α = （1,1）-（2,3）=（-1，-2）
C点要想朝着A点靠近那就是（x,y）=(2,3)+rand*(-1,-2)

2.1 基本原理

假设一个D维的目标搜索空间，有N个粒子组成一个群落，其中第i个例子表示为1个D 维向量：
$X_i=(x_{i1},x_{i2},\dots ,x_{iD}),i=1,2,\dots ,N$
第i个粒子的“飞行”速度也是一个D维的向量，记为：
$V_i=(v_{i1},v_{i2},\dots ,v_{iD}),i=1,2,\dots ,N$
在第t代的第i个粒子向第t+1代进化时，根据如下式子更新：
$V_{ij}(t+1)=w{v}_{ij}(t)+c_1{r}_1(t)[p_{ij}(t)-x_{ij}(t)]+c_2{r}_2(t)[p_{gj}(t)-x_{ij}(t)]$
$x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)$
其中：
w-惯性系数
$c_1、c_2$-学习因子
$p_{ij}$-第i个粒子迄今为止搜索到的最优位置
$p_{gj}$-整个粒子群迄今为止搜索到的最优位置

假设某粒子当前位置C，个体极值位置B，全局最优位置A，那么该粒子下一步的运动状态如右图所示：
$v_{ij}(t+1)=w{v}_{ij}(t)+c_1{r}_1(t)[p_{ij}(t)-x_{ij}(t)]+c_2{r}_2(t)[p_{gj}(t)-x_{ij}(t)]$
B对应$p_i$,A对应$p_g$,黄色向量为当前速度方向，绿色向量为向个体极值飞行步长，红色为向全局最值飞行步长。

2.2 算法流程

输入 ： 参数 w：0.5 − 0.8, $c_1,c_2$：0.1 − 2, $v_{max},v_{min}$：取决于优化函数
Step1：初始化种群x;
Step2：计算个体适应度*；
Step3：更新粒子速度−>更新粒子位置*；
Step4：并计算新位置的适应度，若新位置适应度更高，则将该粒子的位置进行更新，否则不g更新。
Step5：判断是否满足终止条件* ，是则退出，否则返回Step2。
输出 ：输出最优值。
注：
*1.一般的，我们优化目标是最小化一个函数值。所以个体计算出的函数值越小，适应度越高。 max f= min -f;
*2.注意更新速度后，先进行速度边界检测，一般采用($v(v>v_{max})=v_{max}$，位置同理。
*3.常见终止条件为设定迭代进化次数、适应度n代不再变化等。

3.算法分析

3.1 算法的优缺点

**优点：**原理比较简单，实现容易，参数少。
**缺点：**易早熟收敛至局部最优、迭代后期收敛速度慢的。
**解释：**标准粒子群算法的参数是固定的。w描述的是粒子的“惯性”，在进化前期w应该大一些，保证各个粒子独立飞行充分搜索空间，后期应该小一点，多向其他粒子学习。 $c_1,c_2$分别向个体极值和全局极值最大飞行步长。前期$c_1$应该大一些，后期$c_2$应该大一些，这样就能平衡粒子的全局搜索能力和局部搜索能力。3个参数共同影响了粒子的飞行方向，导致即使其他粒子找到更好的，但是当前粒子惯性太大，不能很快的飞向更优的位置。

3.2 算法改进方向

针对标准PSO的缺点，通常有如下的改进：
1.实现参数的自适应变化。
2.引入一些其他机制，比如随机的因素，速度、位置的边界变化−后期压缩最大速度等。
3.结合其他智能优化算法：遗传算法、免疫算法、模拟退火算法等等，帮助粒子跳出局部最优，改善收敛速度。

4.粒子群算法案例

4.1 经典粒子群算法

求解$f(x)={\sum }_{i=1}^{30}{x}_i^2+x_i-6,x\in [-5,5]$的最小值
主程序

clear, clc
x=zeros(1,30);
[xm1,fv1]=PSO(@fitness,50,1.5,2.5,0.5,100,30);

粒子群程序

function [xm,fv]=PSO(fitness,N,c1,c2,w,M,D)

%%%%%%给定初始化条件%%%%%%
% fitness为待优化的目标函数，N初始化群体个体数目，c1学习因子1，c2学习因子2
% w惯性权重，M最大迭代次数，D搜索空间维数

%%%%%初始化种群的个体(可以在这里限定位置和速度的范围) %%%%%%
format long;
for i=1:N
    for j=1:D
        x(i,j)=randn; %随机初始化位置
        v(i,j)=randn; %随机初始化速度
    end
end

%%%%%%先计算各个粒子的适应度，并初始化Pi和Pg%%%%%%
for i=1:N
    p(i)=fitness(x(i,:));
    y(i,:)=x(i,:);
end
pg=x(N,:);             %Pg为全局最优
for i=1:(N-1)
    if fitness(x(i,:)) < fitness(pg)
        pg=x(i,:);
    end
end

%%%%%进入主要循环，按照公式依次迭代，直到满足精度要求%%%%%
for t=1:M
    for i=1:N     %更新速度、位移
        v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));
        x(i,:)=x(i,:)+v(i,:);
        if fitness(x(i,:)) < p(i)
            p(i)=fitness(x(i,:));
            y(i,:)=x(i,:);
        end
        if p(i)<fitness(pg)
            pg=y(i,:);
        end
    end
    Pbest(t)=fitness(pg);
end

%%%%%%给出计算结果%%%%%%
disp('*************************************************')
disp('目标函数取最小值时的自变量：')
xm=pg'
disp('目标函数的最小值为：')
fv=fitness(pg)
disp('**************************************************')
end

适应度计算（目标函数）

%%%计算适应度%%%%%
function F=fitness1(x)
F=0;
for i=1:30
    F=F+x(i)^2+x(i)-6;
end

计算结果
xm =
[ -0.697544201249697, -0.509025429406022, -0.392308765094765,
-0.089950059666007, -0.321688272235924, -0.427778683637105,
-0.094787385478813, -0.414951133527019, -0.496865674920041,
-0.625906350404268, -0.628946621716817, -0.037839073769624
-0.351674676096903, -0.485035047045901, -0.544102067158362
-0.315924880143216, 0.164258582862410, -0.188071066565715
0.079206880854935, -0.476030997755836, -0.454165768499641
-0.093068482006820, -0.050448197165836, -0.031214413235593
-0.814724529967273, -0.726405278717083, -0.410443800423710
-0.516996702435542, -0.443812270919015, -0.436739542173254]
-0.697544201249697, -0.509025429406022, -0.392308765094765
-0.089950059666007, -0.321688272235924, -0.427778683637105
-0.094787385478813, -0.414951133527019, -0.496865674920041
-0.625906350404268, -0.628946621716817, -0.037839073769624
-0.351674676096903, -0.485035047045901, -0.544102067158362
-0.315924880143216, 0.164258582862410, -0.188071066565715
0.079206880854935, -0.476030997755836, -0.454165768499641
-0.093068482006820, -0.050448197165836, -0.031214413235593
-0.814724529967273, -0.726405278717083, -0.410443800423710
-0.516996702435542, -0.443812270919015, -0.436739542173254]
目标函数的最小值为：
fv = -1.851386556690094e+02

4.2 自适应粒子群算法

自适应粒子群算法（Adaptive Particle Swarm Optimization，APSO）主要是依据群体的收敛情况动态调整算法的参数。
在粒子群算法中，惯性权重w是最重要的参数，增大w可以提高算法的全局搜索能力，减小w可以提高算法的局部搜索能力。
惯性权重w大小和其距全局最优点的距离有关。可以采用非线性动态惯性权重系数公示为：

其中：
f-粒子实时的目标函数值
$f_{avg}、f_{min}$-当前所有粒子的平均值和最小目标值
当粒子目标值分散时，减小w; 集中时，增加w.
例 $f(x)=\frac{si{n}^2\sqrt{x_1^2+x_2^2}-cos\sqrt{x_1^2+x_2^2}+1}{[1+0.1(x_1^2+x_2^2){]}^2}-0.7$的最小值。
主程序

[xm,fv] = PSO_adaptation(@AdaptFunc,50,2,2,0.8,0.6,100,2)

自适应粒子群算法

function [xm,fv] = PSO_adaptation(fitness,N,c1,c2,wmax,wmin,M,D)
% fitness为待优化的目标函数，N初始化群体个体数目，c1学习因子1，c2学习因子2
% w惯性权重，wmax惯性权重最大值，wmin惯性权重最值小，M最大迭代次数，D搜索空间维数
% xm目标函数取最小值时的自变量，fv目标函数最小值
format short;

%%%%%%初始化种群的个体%%%%%%
for i=1:N
    for j=1:D
        x(i,j)=randn;
        v(i,j)=randn;
    end
end

%%%%%%先计算各个粒子的适应度%%%%%%
for i=1:N
    p(i)=fitness(x(i,:));
    y(i,:)=x(i,:);
end
pg=x(N,:);             %Pg表示全局最优
for i=1:(N-1)
    if fitness(x(i,:))<fitness(pg)
        pg=x(i,:);
    end
end

%%%%%%进入主要循环%%%%%%
for t=1:M
    for j=1:N
        fv(j) = fitness(x(j,:));
    end
    fvag = sum(fv)/N;
    fmin = min(fv);
    for i=1:N
        if fv(i) <= fvag
            w = wmin + (fv(i)-fmin)*(wmax-wmin)/(fvag-fmin);
        else
            w = wmax;
        end
        v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));
        x(i,:)=x(i,:)+v(i,:);
        if fitness(x(i,:))<p(i)
            p(i)=fitness(x(i,:));
            y(i,:)=x(i,:);
        end
        if p(i)<fitness(pg)
            pg=y(i,:);
        end
    end
end
xm = pg';          %目标函数取最小值时的自变量
fv = fitness(pg);     %目标函数最小值
end

适应度计算（目标函数）

function y=AdaptFunc(x)
   xx=x(1)^2+ x(2)^2;
   sin2x=sin(sqrt(xx));cosx=cos(sqrt((xx)));
   y=(sin2x^2-cosx+1)/((1+0.1*xx)^2)-0.7;
end

计算结果：
xm =[1.7840e+04, -9.595e+04]
fv= -0.7000

挖个坑：粒子群算法与遗传算法、免疫算法、模拟退火算法的混合算法，有时间再写!