一、简介

1.什么是粒子群算法？

粒子群算法模仿昆虫、兽群、鸟群和鱼群等的群集行为，这些群体按照一种合作的方式寻找食物，群体中的每个成员通过学习它自身的经验和其他成员的经验来不断改变其搜索模式。

2.粒子群算法如何工作？

以群鸟的觅食过程为例，每只鸟的初始状态都是处于随机位置，且飞翔的方向也是随机。每只鸟都不知道食物在哪里，但是随着时间的推移，这些初始处于随机位置的鸟类通过群内的相互学习，信息共享和个体不断积累寻觅食物的经验，自发组织积聚成一个群落，并朝着食物前进。每只鸟能通过一定经验和信息估计目前所处的位置对于能寻找到食物有多大的价值，即适应值；每只鸟都能记住自己所找到的最好的位置，称为局部最优，那么群体内所有鸟的最好的位置则称为全局最优，整个鸟群的觅食中心朝着全局最优移动。

3.粒子群模型

在群鸟觅食模型中，每个个体都能被看成一个例子，则鸟群可以被看成一个粒子群。假设在一个D维目标空间当中。

• 有m个粒子组成一个群体，其中第i个粒子的位置表示为：

  
  第i个粒子的位置

• 将位置X带入目标函数可以算出其适应值，根据适应值的大小可以衡量其优劣。同理可设粒子个体经过的最好的位置为：

  
  粒子个体经过的最好的位置
• 所有粒子经过的最好的位置： 所有粒子经过的最好的位置

• 粒子速度：

  
  粒子速度

• 采用下列公式对粒子位置进行更新：

  
  速度更新公式
  位置更新公式
  

  其中α为约束因子，控制速度的权重；r1/2是[0,1]范围内的随机数。

4.参数介绍与值推荐

3中的一些参数与其推荐值如下：

1. 粒子数(m)

粒子数一般取值为20~40.实验表明，对于多数问题，30个粒子就够用了。

2.惯性因子(w)

w越大，粒子飞翔幅度越大，容易错失局部寻优能力，而全局搜索能力越强；反之，则局部寻优能力越强，而局部寻优能力减弱。
如果惯性因子是变量，通常在迭代开始时将惯性因子w设置得较大，然后再迭代过程中逐步减少。这样使得粒子群在开始优化时搜索较大的解空间，得到合适的种子，然后再在后期逐渐收缩到较好的区域进行更加精细的搜索。
如果惯性因子是常量，推荐取0.6~0.75区间的合理值。

3.加速度常数(c)

对于简单的问题一般取两个c都为2。
c1:表征粒子的自身经验，若为0则容易陷入局部最优。
c2:表征粒子的社会经验，即群体共享信息，若为0则等价于单独运行了m个粒子，因而得到最优解的几率小。

4.最大飞翔速度Vmax

为了跳出局部最优，需要Vmax较大，而在接近最优值时，采用较小的步长则更好。
如果Vmax是固定不变的，通常取其为每维变化范围的10%~20%。

二、算法运行流程

1.初始化操作

初始化粒子群（速度和位置）、惯性因子、加速度常数、最大迭代次数、维数、最大飞翔速度和算法终止的最小允许误差。

2.评价每个粒子的初始值，即将位置带入适应性函数中求解

3.将初始值作为当前粒子的局部最优值，并将各适应值对应的位置作为每个粒子的局部最优值所在的位置。

4.将最佳初始适应值作为当前全局最优值，并将最佳适应值对应的位置作为全局最优值所在的位置。

5.根据速度更新公式更新当前粒子的速度，但注意不能超过设置的粒子最大速度。

6.根据位置更新公式更新各个粒子的位置。

7.比较当前粒子位置的适应值是否比其历史最优值好，如果好则替代。

8.找出全局最优

9.重复5~8直到满足误差或者迭代完成

三、demon

求该函数最大值

matlab代码如下：

tic;%程序运行计时
E0 = 0.001;%误差
MaxNum = 100;%粒子最大迭代次数
variableNum = 1; %变量数目
particleSize = 30;%粒子群规模
c1 = 2; c2 = 2;%学习因子都设为2
w = 0.6;%惯性因子
vmax = 1; %最大飞翔速度,通常选取为变化范围的0.1-0.2
x = -5 + 10 * rand(particleSize,variableNum);%初始化粒子所在的位置
v = 2 * rand(particleSize,variableNum);
fitness = inline('-(2.1*(1-x+2*x.^2).*exp(-x.^2/2))','x');
for i = 1:particleSize
    for j = 1:variableNum
        f(i) = fitness(x(i,j));
    end
end
personalbest_x = x; %个体最优的位置
personalbest_y = f; %个体最优的值
[globalbest_y i] = min(personalbest_y);%返回personalbest_y的最小值给globalbest_y,i是该最小值存在的位置
globalbest_x = personalbest_x(i,:);
k = 1; %迭代次数
while k<=MaxNum
    for i = 1:particleSize
        for j = 1:variableNum
            f(i) = fitness(x(i,j));
        end
        if f(i)<personalbest_y(i) %判断是否当前位置的值小于当前的最佳位置(因为有取倒数)
            personalbest_y(i) = f(i);
            personalbest_x(i,:) = x(i,:);
        end
    end
    [globalbest_y i] = min(personalbest_y);%更新位置
    globalbest_x = personalbest_x(i,:);%更新位置
    for i = 1:particleSize %更新当前的位置
        v(i,:) = w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))+c2*rand*(globalbest_x-x(i,:));
        for j = 1:variableNum
            if v(i,j)>vmax
                v(i,j) = vmax;
            elseif v(i,j) < -vmax
                v(i,j) = -vmax;
            end
        end
        x(i,:) = x(i,:) + v(i,:);
    end
%    if abs(globalbest_y)<E0,break,end 
    k = k+1;
end
Value1 = -globalbest_y; Value1 = num2str(Value1);
disp(strcat('the maximum value = ',Value1));
Value2 = globalbest_x; Value2 = num2str(Value2);
disp(strcat('the corresponding coordinate ','=',Value2))
x = -5:0.01:5;
y = 2.1*(1-x+2*x.^2).*exp(-x.^2/2);
plot(x,y,'m-','linewidth',3);
hold on;
plot(globalbest_x,-globalbest_y,'kp','linewidth',4);
legend('目标函数','搜索到的最大值');xlabel('x');ylabel('y');grid on; toc; ```