粒子群算法（PSO）

一、算法理论

        粒子群优化算法来源于鸟类群体活动的规律性，进而利用群体智能建立一个简化的模型。它模拟鸟类的觅食行为，将求解问题的搜索空间比作鸟类的飞行空间，将每只鸟抽象成一个没有质量和体积的粒子，用它来表征问题的一个可能解，将寻找问题最优解的过程看成鸟类寻找食物的过程，进而求解复杂的优化问题。粒子群优化算法与其他进化算法一样，也是基于“种群”和“进化”的概念，通过个体间的协作与竞争，实现复杂空间最优解的搜索。同时，它又不像其他进化算法那样对个体进行交叉、变异、选择等进化算子操作，而是将群体中的个体看作在D维搜索空间中没有质量和体积的粒子，每个粒子以一定的速度在解空间运动，并向自身历史最佳位置pbest和邻域历史最佳位置gbest聚集，实现对候选解的进化。粒子群算法具有很好的生物社会背景而易于理解，由于参数少而容易实现，对非线性、多峰问题均具有较强的全局搜索能力，在科学研究与工程实践中得到了广泛关注。

二、粒子群算法的特点

1、基于群智能理论的优化算法

        通过群体中粒子间的合作与竞争产生的群体智能指导优化搜索。与其他算法相比，粒子群算法是一种高效的并行搜索算法。

2、采用的速度-位移模型

        与遗传算法都是随机初始化种群，使用适应值来评价个体的优劣程度和进行一定的随机搜索。但粒子群算法根据自己的速度来决定搜索，没有遗传算法的交叉与变异。与进化算法相比，粒子群算法保留了基于种群的全局搜索策略，但是其采用的速度-位移模型操作简单，避免了复杂的遗传操作。

3、可获得最优解与次优解

        由于每个粒子在算法结束时仍保持其个体极值，即粒子群算法除了可以找到问题的最优解外，还会得到若干较好的次优解，因此将粒子群算法用于调度和决策问题可以给出多种有意义的方案。

4、可以动态调整搜索策略

        粒子群算法特有的记忆使其可以动态地跟踪当前搜索情况并调整其搜索策略。另外，粒子群算法对种群的大小不敏感，即使种群数目下降时，性能下降也不是很大。

三、基本粒子群算法

        假设在一个D维的目标搜索空间中，有N个粒子组成一个群落，其中第i个粒子表示为一个D维的向量：

 第i个粒子的“飞行”速度也是一个D维的向量，记为

 第i个粒子迄今为止搜索到的最优位置称为个体极值，记为

 整个粒子群迄今为止搜索到的最优位置为全局极值，记为

 在找到这两个最优值时，粒子根据如下的式（6.5）和式（6.6）来更新自己的速度和位置：

         其中：c1和c2为学习因子，也称加速常数；r1和r2为[0，1]范围内的均匀随机数， i=1，2，…，D；是粒子的速度，∈[-，]，vmax是常数，由用户设定来限制粒子的速度。r1和r2是介于0和1之间的随机数，增加了粒子飞行的随机性。式（6.5）右边由三部分组成：第一部分为“惯性”或“动量”部分，反映了粒子的运动“习惯”，代表粒子有维持自己先前速度的趋势；第二部分为“认知”部分，反映了粒子对自身历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；第三部分为“社会”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻域历史最佳位置逼近的趋势。

四、标准粒子群算法

        如何确定局部搜索能力和全局搜索能力的比例，对一个问题的求解过程很重要。1998年，Shi Yuhui等人提出了带有惯性权重的改进粒子群算法，由于该算法能够保证较好的收敛效果，所以被默认为标准粒子群算法。其进化过程为：

         在式（6.7）中，第一部分表示粒子先前的速度，用于保证算法的全局收敛性能；第二部分、第三部分则使算法具有局部收敛能力。可以看出，式（6.7）中惯性权重w表示在多大程度上保留原来的速度：w较大，则全局收敛能力较强，局部收敛能力较弱；w较小，则局部收敛能力较强，全局收敛能力较弱。

        当w=1时，式（6.7）与式（6.5）完全一样，表明带惯性权重的粒子群算法是基本粒子群算法的扩展。实验结果表明：w在0.8～1.2之间时，粒子群算法有更快的收敛速度；而当w＞1.2时，算法则容易陷入局部极值。

        另外，在搜索过程中可以对w进行动态调整：在算法开始时，可给w赋予较大正值，随着搜索的进行，可以线性地使w逐渐减小，这样可以保证在算法开始时，各粒子能够以较大的速度步长在全局范围内探测到较好的区域；而在搜索后期，较小的w值则保证粒子能够在极值点周围做精细的搜索，从而使算法有较大的概率向全局最优解位置收敛。对w进行调整，可以权衡全局搜索和局部搜索能力。目前，采用较多的动态惯性权重值是Shi提出的线性递减权值策略，其表达式如下：

         式中：表示最大进化代数；表示最小惯性权重；表示最大惯性权重；t表示当前迭代次数。在大多数的应用中，=0.9，=0.4。

五、改进粒子群算法

1、压缩因子粒子群算法

        Clerc 等人提出利用约束因子来控制系统行为的最终收敛，该方法可以有效搜索不同的区域，并且能得到高质量的解。压缩因子法的速度更新公式为：

         式中，λ为压缩因子

         实验结果表明：与使用惯性权重的粒子群优化算法相比，使用具有约束因子的粒子群算法具有更快的收敛速度。

2、离散粒子群算法

        基本的粒子群算法是在连续域中搜索函数极值的有力工具。继基本粒子群算法之后，Kennedy和Eberhart又提出了一种离散二进制版的粒子群算法。在此离散粒子群方法中，将离散问题空间映射到连续粒子运动空间，并适当修改粒子群算法来求解，在计算上仍保留经典粒子群算法速度-位置更新运算规则。粒子在状态空间的取值和变化只限于0和1两个值，而速度的每一维代表位置每一位取值为1的可能性。因此在连续粒子群中的更新公式依然保持不变，但是和只在[0,1]内取值。其位置更新公式表示如下：

 式中，r是从U（0,1）分布中产生的随机数。

五、粒子群算法的实现

（1）初始化粒子群，包括群体规模N，每个粒子的位置xi和速度vi。
（2）计算每个粒子的适应度值fit[i]。
（3）对每个粒子，用它的适应度值fit[i]和个体极值pbest（i）比较。如果fit[i] ＜pbest（i），则用fit[i]替换掉pbest（i）。
（4）对每个粒子，用它的适应度值fit[i]和全局极值gbest比较。如果fit[i]＜gbest，则用fit[i]替换gbest。
（5）迭代更新粒子的速度vi和位置xi。
（6）进行边界条件处理。
（7）判断算法终止条件是否满足：若是，则结束算法并输出优化结果；否则返回步骤（2）。
粒子群算法的运算流程如图所示。

六、关键参数说明

关键参数	取值	原因
粒子种群规模 N	一般20-50	视问题规模而定，对于比较难的问题或者特定类型的问题，粒子的数量可以取到100或200。另外，粒子数目越大，算法搜索的空间范围就越大，也就更容易发现全局最优解；当然，算法运行的时间也越长。
惯性权重 w	[0.8，1.2]	惯性权重的大小表示了对粒子当前速度继承的多少。当惯性权重值较大时，全局寻优能力较强，局部寻优能力较弱；当惯性权重值较小时，全局寻优能力较弱，局部寻优能力较强。惯性权重的选择通常有固定权重和时变权重。固定权重就是选择常数作为惯性权重值，在进化过程中其值保持不变，时变权重则是设定某一变化区间，在进化过程中按照某种方式逐步减小惯性权重。时变权重的选择包括变化范围和递减率。固定的惯性权重可以使粒子保持相同的探索和开发能力，而时变权重可以使粒子在进化的不同阶段拥有不同的探索和开发能力。
加速常数 c 1和 c 2	c1=c2=1.5	加速常数c1和c2分别调节向pbest和gbest方向飞行的最大步长，它们分别决定粒子个体经验和群体经验对粒子运行轨迹的影响，反映粒子群之间的信息交流。如果c1=c2=0，则粒子将以当前的飞行速度飞到边界。此时，粒子仅能搜索有限的区域，所以难以找到最优解。如果 c1=0，则为“社会”模型，粒子缺乏认知能力，而只有群体经验，它的收敛速度较快，但容易陷入局部最优；如果c2=0，则为“认知”模型，没有社会的共享信息，个体之间没有信息的交互，所以找到最优解的概率较小，一个规模为 D 的群体等价于运行了 N 个各行其是的粒子。因此一般设置c1=c2，通常可以取c1=c2=1.5。这样，个体经验和群体经验就有了同样重要的影响力，使得最后的最优解更精确。
粒子的最大速度 v max	[-，+]	如果该值太大，则粒子们也许会飞过优秀区域；而如果该值太小，则粒子们可能无法对局部最优区域以外的区域进行充分的探测。它们可能会陷入局部最优，而无法移动足够远的距离而跳出局部最优，达到空间中更佳的位置。研究者指出，设定vmax和调整惯性权重的作用是等效的，所以vmax一般用于对种群的初始化进行设定，即将vmax设定为每维变量的变化范围，而不再对最大速度进行细致的选择和调节。
停止准则	最大迭代次数、计算精度或最优解的最大停滞步数Δt（或可以接 受的满意解）通常认为是停止准则	停止准则的设定需同时兼顾算法的求解时间、优化质量和 搜索效率等多方面性能。
邻域结构的设定		全局版本的粒子群算法将整个群体作为粒子的邻域，具有收敛速度快的优点，但有时算法会陷入局部最优。局部版本的粒子群算法将位置相近的个体作为粒子的邻域，收敛速度较慢，不易陷入局部最优值。实际应用中，可先采用全局粒子群算法寻找最优解的方向，即得到大致的结果，然后采用局部粒子群算法在最优点附近进行精细搜索。
边界条件处理		比如通过设置最大位置限制xmax和最大速度限制vmax，当超过最大位置或最大速度时，在范围内随机产生一个数值代替，或者将其设置为最大值，即边界吸收。

七、MATLAB仿真实例

%%%%%%%%%f(x，y)=3*cos(x*y)+x+y*y适应值计算%%%%%%%%%%
clear all;              %清除所有变量
close all;              %清图
clc;                    %清屏
x=-4:0.02:4;
y=-4:0.02:4;
N=size(x,2);
for i=1:N
    for j=1:N
        z(i,j)=3*cos(x(i)*y(j))+x(i)+y(j)*y(j);
    end
end
mesh(x,y,z)
xlabel('x')
ylabel('y')

%%%%%%%%%%%%%%%%%粒子群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %清除所有变量
close all;              %清图
clc;                    %清屏
N=100;                  %群体粒子个数
D=2;                    %粒子维数
T=200;                  %最大迭代次数
c1=1.5;                 %学习因子1
c2=1.5;                 %学习因子2
Wmax=0.8;               %惯性权重最大值
Wmin=0.4;               %惯性权重最小值
Xmax=4;                 %位置最大值
Xmin=-4;                %位置最小值
Vmax=1;                 %速度最大值
Vmin=-1;                %速度最小值
%%%%%%%%%%%%%%%%初始化种群个体（限定位置和速度）%%%%%%%%%%%%%%%%
x=rand(N,D) * (Xmax-Xmin)+Xmin;
v=rand(N,D) * (Vmax-Vmin)+Vmin;
%%%%%%%%%%%%%%%%%%初始化个体最优位置和最优值%%%%%%%%%%%%%%%%%%%
p=x;
pbest=ones(N,1);
for i=1:N
    pbest(i)=func2(x(i,:));
end
%%%%%%%%%%%%%%%%%%%初始化全局最优位置和最优值%%%%%%%%%%%%%%%%%%
g=ones(1,D);
gbest=inf;
for i=1:N
    if(pbest(i)<gbest)
        g=p(i,:);
        gbest=pbest(i);
    end
end
gb=ones(1,T);
%%%%%%%%%%%按照公式依次迭代直到满足精度或者迭代次数%%%%%%%%%%%%%
for i=1:T
    for j=1:N
        %%%%%%%%%%%%%%更新个体最优位置和最优值%%%%%%%%%%%%%%%%%
        if (func2(x(j,:))<pbest(j))
            p(j,:)=x(j,:);
            pbest(j)=func2(x(j,:));
        end
        %%%%%%%%%%%%%%%%更新全局最优位置和最优值%%%%%%%%%%%%%%%
        if(pbest(j)<gbest)
            g=p(j,:);
            gbest=pbest(j);
        end
        %%%%%%%%%%%%%%%%计算动态惯性权重值%%%%%%%%%%%%%%%%%%%%
        w=Wmax-(Wmax-Wmin)*i/T;
        %%%%%%%%%%%%%%%%%跟新位置和速度值%%%%%%%%%%%%%%%%%%%%%
        v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
            +c2*rand*(g-x(j,:));
        x(j,:)=x(j,:)+v(j,:);
        %%%%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D
            if (v(j,ii)>Vmax)  |  (v(j,ii)< Vmin)
                v(j,ii)=rand * (Vmax-Vmin)+Vmin;
            end
            if (x(j,ii)>Xmax)  |  (x(j,ii)< Xmin)
                x(j,ii)=rand * (Xmax-Xmin)+Xmin;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%记录历代全局最优值%%%%%%%%%%%%%%%%%%%%%
    gb(i)=gbest;
end
g;                         %最优个体         
gb(end);                   %最优值
figure
plot(gb)
xlabel('迭代次数');
ylabel('适应度值');
title('适应度进化曲线')

源自《智能优化算法及其MATLAB实例（第2版）》