粒子群算法的寻优算法

一、粒子算法的概述
粒子群算法是一种智能优化算法。关于智能，个人理解，不过是在枚举法的基础上加上了一定的寻优机制。试想一下枚举法，假设问题的解空间很小，比如一个函数 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来表示，也就是每个粒子的更新速度了，公式化的表示就是这样的：

从上面的速度V的更新而已看到，c1那项就是朝着自己的最优解前进，c2那一项就是朝着全局最优解那前进。用简单的图表示如下：

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

初始化粒子群个体；
计算每个个体的适应度值（函数值）作为评判好坏的标准；
找到每个个体自己在所有迭代过程中的最优解Pbest；
找到所有个体在所有迭代过程中的最优解Zbest；
根据速度公式更新速度；
根据位置公式更新位置；
重复步骤二直至迭代次数结束
这里有几个参数需要说一下，

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

三、代码实现
%% 清空环境
clc
clear
%% 参数初始化
%粒子群算法中的三个参数
c1 = 1.49445;%加速因子
c11=c1;
c2 = 1.49445;
c22=c2;
w=0.8 %惯性权重
%w1=0.9 %希望改变的权重
%w0=w;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模1

Vmax=1;       %限制速度围
Vmin=-1;     
popmax=5;    %变量取值范围
popmin=-5;
dim=10;       %适应度函数维数1

func=1;       %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
%Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出

%% 产生初始粒子和速度
for i=1:sizepop
%随机产生一个种群
pop(i,:)=popmax*rands(1,dim);    %初始种群
V(i,:)=Vmax*rands(1,dim);             %初始化速度
                                 %计算适应度
fitness(i)=fun(pop(i,:),func);   %粒子的适应度
end

%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:);   %全局最佳
pbest=pop;                %个体最佳
fitnesspbest=fitness;     %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值

%% 迭代寻优
for i=1:maxgen  
fprintf('第%d代，',i);
fprintf('最优适应度%f\n',fitnessgbest);
%w=w+(w1-w0)/maxgen;
c1=c1+(1-c11)/maxgen;
%c2=c2+(2-c22)/maxgen;
for j=1:sizepop
    
    %速度更新
    V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
    V(j,find(V(j,:)>Vmax))=Vmax;   %限制速度不能太大
    V(j,find(V(j,:)<Vmin))=Vmin;
    
    %种群更新
    pop(j,:)=pop(j,:)+0.5*V(j,:);       %位置更新
    pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
    pop(j,find(pop(j,:)<popmin))=popmin;
    
    if rand>0.98                         %加入变异种子，用于跳出局部最优值
        pop(j,:)=rands(1,dim);
    end
    
    %更新第j个粒子的适应度值
    fitness(j)=fun(pop(j,:),func); 

end

for j=1:sizepop
    
    %个体最优更新
    if fitness(j) < fitnesspbest(j)
        pbest(j,:) = pop(j,:);
        fitnesspbest(j) = fitness(j);
    end
    
    %群体最优更新
    if fitness(j) < fitnessgbest
        gbest = pop(j,:);
        fitnessgbest = fitness(j);
    end
end 
yy(i)=fitnessgbest;    
    
end
%% 结果分析
figure;
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

四、实验数据与分析
控制参数

A组 （w惯性权重 c1 c2加速因子）
1、改变w的值

w=0.8 w1=0.99(平均最佳适应度为4.39)

w=0.8 w1=0.6（最佳适应度为1.79）
w=0.8 w1=0.4(平均最佳适应度为2.59)

*w=0.8 w1=0.2(平均最佳适应度为2.3)

w=0.6 w1=0.2(平均最佳适应度为0.2)

由上述图片数据分析可得在1000次迭代中 w的变化0.8-0.6 优于 0.8-0.4和0.8-0.2 可知一开始惯性权重过大的话会使结果较容易陷入局部最优 而0.6-0.2的最佳适应度只有0.2 说明这个是该条件下的相对最优解 可知 在迭代次数多的后半程惯性权重小 依赖于c1 c2的话对得到最优解有很大的帮助。

2、改变c1的值
c1是1.49-1(平均最佳适应度为3.9）

c1是1.49-1.2(平均最佳适应度为1.79）

c1是1.49-2(平均最佳适应度为1.2）

由上述图片和数据可知 当c1 逐渐变小时 若变得太小容易陷入局部最优 而当c1 逐渐变大时 个人最优的占比变大时 得出的结果比较好
3、改变c2的值
c2是1.49-2(平均最佳适应度为2.2）

c2是1.49-1.8(平均最佳适应度为4.2）

c2是1.49-1(平均最佳适应度为1.2）

由上述数据可得 当c2增加时 明显的看出实验结果容易陷入局部最优，而当c2 减小时也就是迭代次数越大时群体的影响越小时得出的结果越理想。

综上所述，对于c1要随着迭代次数变小，c2要随着迭代次数变大，因为c1是个人最优在前期的时候个人最优要优于群体最优，而在后期c2是群体最优要由于个人最优，由此产生的想法，当c1 2-1和c2 1-2；

发现没有屁用。

B组 （sizepop种群规模 dim适应度函数维数）
1、改变sizepop的值

s=50(平均最佳适应度为3.8）

s=100(平均最佳适应度为0.66）

s=200(平均最佳适应度为3）

s=300(平均最佳适应度为1.7）

s=500(平均最佳适应度为1）

通过上述图片和数据分析可知 种群规模过小时，结果容易陷入局部最优，种群多样性的减少对于种群寻找最优解有着一定影响。当种群规模过大时, 虽然最终得出的结果情况较好，但与耗费的资源不成比例。所以种群规模还是适当的好，大概在200-300。

2、改变dim的值

dim=2(平均最佳适应度为0）

dim=6(平均最佳适应度为0.33）

dim=10(平均最佳适应度为4）

dim=15(平均最佳适应度为7.8）

由上述图片和数据可知当维度越小的时候越容易得出全局最优解 越大的时候则相反