简单粒子群优化算法的函数寻优

原文链接：http://blog.csdn.net/on2way/article/details/44699223

一：关于粒子群算法

粒子群算法是一种智能优化算法。关于智能，个人理解，不过是在枚举法的基础上加上了一定的寻优机制。试想一下枚举法，假设问题的解空间很小，比如一个函数 y = x^2 ，解空间在[-1,1],现在求这个函数的最小值，我们完全可以使用枚举法，比如在这里，在解空间[-1,1]上，取1000等分，也就是步长为0.002，生成1000个x值，然后代入函数中，找到这1000个最小的y就可以了。然而实际情况不是这样的，比如为什么选1000等分，不是1w，10w等分，很显然等分的越大，计算量也就越大，带来的解当然也就越精确，那么实际问题中如何去平衡这两点呢？也就是既要计算量小（速度快），也要准确（精度高），这就是智能算法的来源了，一般的智能算法基本上都是这样的，在很大的搜索空间上，即保证了速度快，也能比较好的找到最优解。

再来看看粒子群算法（也称PSO算法），也是一种进化算法，模拟生物群体的觅食行为，是一种群体智能算法，类似的算法想遗传算法，模拟退火算法等等。PSO是通过当前已知种群寻找到的所有解来决定新的解的寻找方向，也就是新解的生成方式依赖于这些种群历史上寻找的所有解。

形象的理解比如下图： 
 
开始随机生成一堆种群，那么这些种群之间的每个个体可以相互交流，比如下一时刻，A告诉B说我的解比你好，那么B就往A那个地方飞，也就是B的解朝着A的解方向变化，当然所有粒子间都这样操作，想想一旦粒子群中间有一个粒子找到了一个最优解，是不是所有的粒子会一窝蜂朝着这个方向而去了，而在这个去的过程中，万一某个粒子找到了一个更好的解，那它还会走吗？不会了，它就告诉剩下的所有粒子说我的解更好呀，大家快来呀（很无私的），然后所有粒子又一窝蜂的照着这个粒子方向前进，当然在这个前进的过程中可能又会产生新的解，就这样一步步的迭代，最终慢慢的趋近于一个最优解，这个解是不是全局最优解，不知道，可能是，也可能不是，取决于原始问题的复杂程度，也取决于粒子前进的多少等等。

粒子群算法相对于其他算法来说还是有很多优点的，典型的就是计算速度很快，在每次迭代时，所有粒子同时迭代，是一种并行计算方式，而且粒子的更新方式简单，朝着一个优秀解方向更新。这个优秀解包括两个部分： 
1）一个是朝着自己在迭代的历史上找到的个体最优解gbest前进 
2）一个是朝着群体在得带历史上找到的全体最优解zbest前进 
现在还有一个问题就是每次迭代的时候更新多少呢？也就是自变量的增加步长了，我们用一个速度量V来表示，也就是每个粒子的更新速度了，公式化的表示就是这样的： 

Vid(t+1)=Vidt+c1∗rand∗(Pid−xid(t))+c2∗rand∗(Pgd−xid(t))

其中自变量的更新为：xi(t+1)=xi(t)+Vi(t) 
从上面的速度V的更新而已看到，c1那项就是朝着自己的最优解前进，c2那一项就是朝着全局最优解那前进。用简单的图表示如下： 

二：粒子群的算法步骤

粒子群的核心部分就是上面说到的那两个公式，一个是速度的更新方式，另一个是位置的更新方式，重点还是速度的更新方式； 
总结来说，粒子群的算法步骤如下：

1. 初始化粒子群个体；
2. 计算每个个体的适应度值（函数值）作为评判好坏的标准；
3. 找到每个个体自己在所有迭代过程中的最优解Pbest；
4. 找到所有个体在所有迭代过程中的最优解Zbest；
5. 根据速度公式更新速度；
6. 根据位置公式更新位置；
7. 重复步骤二直至迭代次数结束

这里有几个参数需要说一下，

• 关于速度V，限制速度的范围，比如需要设置一个最大速度，防止更新过快；
• 关于c1与c2，这两个参数代表学习因子，决定跟随历史优秀解的能力；
• 关于粒子数与迭代次数，粒子数一般50-100，迭代次数视问题而定了；

三：matlab粒子群函数寻优

这里简述几种二维函数的寻优来说明粒子群的简单用法。 
为简单起见，把不同的测试函数写在一个m文件中如下：

function f = fun(X)x = X(1);y = X(2);%选择一个测试函数choose = 1;switch choose    case 1        f = sin(sqrt(x.^2+y.^2))./sqrt(x.^2+y.^2) +exp((cos(2*pi*x)+cos(2*pi*y))/2)-2.71289;    case 2        f_schaffer = 0.5-((sin(sqrt(x.^2+y.^2))).^2-0.5)./(1+0.001*(x.^2+y.^2)).^2;        f = f_schaffer;    case 3        f_rastrigin = x.^2-10*cos(2*pi*x)+10 + y.^2-10*cos(2*pi*y)+10;        f = -f_rastrigin;end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

上述为三种不同形势下的函数，用三维表示如下： 
 
 

暂时编辑选择这几种函数，要求取这些函数的最大值的话，可以看到函数存在着许多的局部最优解。下面编辑主函数实验算法的好坏：

clcclearclose all%速度更新c1 = 2;c2 = 2;%maxgen = 300; %最大更新代数sizepop = 20; %种群大小%popmax = 2;  %规定自变量x 的取值范围popmin = -2;Vmax = 0.5;  %规定速度的大小范围Vmin = -0.5;%种群初始化for i = 1:sizepop    pop(i,:) = 2*popmax*rand(1,2)+popmin;    V(i,:) = 2*Vmax*rand(1,2)+Vmin;    fitness(i) = fun(pop(i,:));end[best_fitness,best_index] = max(fitness);zbest = pop(best_index,:);%群体最优个体gbest = pop;%个体最优个体fitness_gbest = fitness;%个体最优个体函数值fitness_zbest = best_fitness;%群体最优函数值%w_start = 0.9;w_end = 0.4;for i = 1:maxgen    w = w_start - (w_start-w_end)*(i/maxgen)^2;    for j = 1:sizepop        %速度更新        V(j,:) = V(j,:) + c1*rand*(gbest(j,:)-pop(j,:)) + ...                          c2*rand*(zbest-pop(j,:));        V(j,find(V(j,:)>Vmax)) = Vmax;        V(j,find(V(j,:)<Vmin)) = Vmin;        %粒子更新        pop(j,:) = pop(j,:) + w*V(j,:);        pop(j,find(pop(j,:)>popmax)) = popmax;        pop(j,find(pop(j,:)<popmin)) = popmin;        %        fitness(j) = fun(pop(j,:));    end    %个体与种群极值更新    for j = 1:sizepop        %个体        if fitness(j) > fitness_gbest(j)            gbest(j,:) = pop(j,:);            fitness_gbest(j) = fitness(j);        end        %种群        if fitness(j) > fitness_zbest            zbest = pop(j,:);            fitness_zbest = fitness(j);        end    end    xx(i,:) = zbest;    result(i) = fitness_zbest;end%% --%% --plot(result);title(['最优函数值：',num2str(fitness_zbest)]);%三维显示出来figure;[x,y] = meshgrid(popmin:0.1:popmax,popmin:0.1:popmax);[m,n] = size(x);for i = 1:m    for j = 1:n        f(i,j) = fun([x(i,j),y(i,j)]);    endendsurf(x,y,f);hold on;plot3(xx(:,1),xx(:,2),result(:),'r*');hold on;plot3(xx(:,1),xx(:,2),result(:));title(['最优解：',num2str(zbest(1)),',',num2str(zbest(2))]);
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
77
78
79
80
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
77
78
79
80

关于程序中的部分解释（显示不出行号，真不方便），关于粒子的更新这块 pop(j,:) = pop(j,:) + w*V(j,:);，这里的w是一个速度更新权值，这里弄了个动态变化的，主要是为了让前期的跟随速度快，达到全局的搜索能力，后期的跟随速度慢，达到局部搜索能力。这里你也可以设置为一个定值，比如就是0.5的话也可以。 
好了现在把f = fun(X)函数中的choose选择1，测试函数1的结果如下：

 
 
这里很明显的可以看到最优解的变化情况，在三维图上也可以明显的看出；

现在把f = fun(X)函数中的choose选择2，并且把自变量的范围改到【-4,4】，也就是popmax = 4; popmin = -4;测试函数2的结果如下： 
 
 
最优函数值是在中间的那个顶峰处取得的，可以看到该算法还是可以取到的，但是并不是每次都能够取到。

再看choose选择3的情况吧，为了看清楚，把画图部分内容改为plot3(xx(:,1),xx(:,2),result(:),’b*’)：

 

可以很清楚看到最优解如何在三维图像上的变化情况吧，实验也可以找到最优解。

至此，在只有两个自变量x1与x2的情况下函数的寻优过程就是这样了，至于更多的自变量，不过在加些就可以了。实验证明粒子群算法不仅能较为准确的找到全局最优解（大多数时候可以），寻找的速度也是相当快的。