从Matlab实例学习粒子群算法

前言

之前提到过， 无论遗传算法、蚁群算法以及今天要学习的粒子群算法，本质上都是对穷搜算法的精进版本，即通过参照一些自然界的规律来提高搜索的效率。 粒子群算法则是参考了鸟群寻找栖息地的自然现象，它的中心思想是每只鸟既根据自己的历史经验寻找栖息地的同时，也会参考同伴们的经验进行搜索。

粒子群算法步骤

1. 初始化 $N$ 个粒子作为一个群落， 分别对应优化问题的$N$个随机可行解 $x_i,i=1,\cdots ,N$。同时初始化每个粒子的速度$v_i$，代表下一步的搜索方向，也可以说是更新方向。
2. 计算每个粒子的适应度值， 也即函数值$f(x_i)$，并更新每个粒子的历史最佳结果$p_i$及群落的历史最佳结果$g$。
3. 根据$p_i$与$g$， 按下式更新$v_i$:
   $$v_i=v_i+c_1{r}_1(p_i-x_i)+c_2{r}_2(g_i-x_i)\text{\hspace{1em}(1)}$$
   这个式子右边的第一项代表运动的惯性， 即只是在上一次搜索的方向是进行一些改进。 $c_1,c_2$代表学习因子，是一个给定常数。 $r_1,r_2$代表随机数，为搜索加入了随机因素。 $p_i-x_i$则代表了粒子基于自身历史经验对方向进行的修正，即朝着历史最好的结果进行搜索。 $g_i-x_i$则代表了粒子也会考虑整个群落的经验进行修正。
4. 更新粒子 $x_i=x_i+v_i$
5. 判断算法迭代是否终止，否则跳转回步骤2.

从算法流程可以看出来，粒子群算法和蚁群算法等本质上还是很一致的，就是并行地从$N$个随机可行解开始搜索，基于一些启发式的思想优化搜索效率。比如粒子群算法里区别于穷搜算法的点就在于会不断参考粒子自身经验和整个种群的经验。

一些增强

• 标准粒子群算法：将(1)改为：
  $$v_i=w{v}_i+c_1{r}_1(p_i-x_i)+c_2{r}_2(g_i-x_i)\text{\hspace{1em}(2)}$$
  也即对惯性引入一个权重，这个权重是所有粒子共享的。 显然，$w$越大，搜索越趋向于随机。 $w$越小，越趋向于纯粹利用经验。
• 压缩粒子群算法： （2）中的$w$由下式给出：
  $$\lambda =\frac{2}{\left|2-\phi -\sqrt{\left({\phi }^2-4\phi \right)}\right|}，\phi =c_1+c_2$$
  实验结果表明可显著提升收敛速度。
• 离散粒子群算法： 用于求解离散问题。 比如，如果变量的取值为0或1， 而粒子群飞行时到达的可行解势必是连续结果。 此时，只需将该结果投影到离其最近的可行解上即可 （也可引入一些随机因素），而速度及飞行的更新式不变。

Matlab代码

仿真实例： 用粒子群算法求解$f(x,y)=3\cos (xy)+x+y^2$, x,y取值范围为$[-4,4]$:

%%%%%%%%%%%%%%%%%粒子群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)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('适应度进化曲线')


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76