粒子群优化算法

一、算法原理

  粒子群优化算法的基本思想：通过群体中个体之间的协作和信息共享来寻找最优解。粒子群优化算法仅通过速度更新公式和位置更新公式不断迭代更新并存储每次迭代的个体最优解和全局最优解，所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置，从而得到整体全局最优解。

二、算法流程及步骤

算法步骤：

(1)初始化所有的个体(粒子),初始化它们的速度和位置,并且将个体的历史最优pBest设为当前位置,而群体中最优的个体作为当前的gBest.
(2)在每一代的进化中,计算各个粒子的适应度函数值。
(3)如果该粒子当前的适应度函数值比其历史最优值要好,那么历史最优将会被当前位置所替代
(4)如果该粒子的历史最优比全局最优要好,那么全局最优将会被该粒子的历史最优所替代。
对每个粒子i的第d维的速度和位置分别按照下列公式进行更新
粒子速度更新公式:
  w：惯性权重
  w×vi：表示上一时刻的速度
  加速系数c1：表示个体置信度权重
  pbest：表示个体历史最优值
  加速系数c1：表示群体置信度权重
  gbest：表示全局最优值
  rand()：是介于0-1的随机数
    
粒子位置更新公式:

  Vi：粒子速度
  Xi：粒子的当前位置
(6)如果还没肴到达结束条件,转到(2),否则输出gBest并结束。

三、matlab代码实现

1、PSO算法实现

PSO.m

%% 参数初始化
%粒子群算法中的三个参数
c1 =1.49445;%加速因子
c2 = 1.49445;
w=0.8  %惯性权重

maxgen=1000;   % 进化次s数  
sizepop=200;   %种群规模

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

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);
    for j=1:sizepop
        %c1=(maxgen-i)/maxgen*0.5+1;%参数随迭代次数改变
        %c2=(maxgen-i)/maxgen*0.5+2;
        %w=(maxgen-i)/maxgen*0.5+0.2;
        %速度更新
        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);

2、适应度函数

2.1 Rastrigin(x)函数

function y = Rastrigin(x)
% Rastrigin函数
% 输入x,给出相应的y值,在x = ( 0 , 0 ,…, 0 )处有全局极小点0.
% 编制人：
% 编制日期：
[row,col] = size(x);
if  row > 1 
    error( ' 输入的参数错误 ' );
end
y =sum(x.^2-10*cos(2*pi*x)+10);

运行结果：

  由图可知，最优适应度仅在在迭代38次就收敛到极值，下降速度很快，但并不是最优极值，说明使用 Rastrigin(x)函数作为适应度函数，会出现过早收敛，易陷入局部最优。

2.2 Schaffer(x)函数

function y=Schaffer(x)
[r,c]=size(x);
if r>1
    error('输入的参数错误');
end
%y1=x(1,1);
%y2=x(1,2);
%temp=y1^2+y2^2;
y=0.5+((sin(sqrt(x(r,1)^2+x(r,c)^2)))^2-0.5)/((1.0+0.001*(x(r,1)^2+x(r,c)^2))^2);
%y=-y;

分析：此适应度函数下，经过20次迭代后，收敛到最优适应度0，比Rastrigin(x)适应度函数收敛到最优值速度更快，但收敛过快，易陷入局部最优。

2.3 Griewank(x)函数

function y=Griewank(x)
%Griewan函数
%输入x,给出相应的y值,在x=(0,0,…,0)处有全局极小点0.
%编制人：
%编制日期：
[row,col]=size(x);
if row>1
    error('输入的参数错误');
end
y1=1/4000*sum(x.^2);
y2=1;
for h=1:col
    y2=y2*cos(x(h)/sqrt(h));
end
y=y1-y2+1;
%y=-y;

运行结果：

  由图可知，在迭代25次后，算法收敛到极值0，此适应度函数收敛速度也较快。

四、参数设置分析

1、改变惯性权重w

参数列表：

• 惯性权重设置为Wi->Wj,使惯性权重随着迭代次数而不断变化由Wi下降到Wj;

自定义的w变化公式：w=(maxgen-i)/maxgen*0.5+0.2;

• 个体置信度权重C1,设置为2->1,随着迭代次数变化C1由2减小到1；

自定义的C1变化公式：c1=(maxgen-i)/maxgen*0.5+1

• 群体置信度权重C2,设置为1->2,随着迭代次数变化C2由1增加到2；

自定义的C2变化公式：c2=(maxgen-i)/maxgen*0.5+2

惯性权重w	个体置信度权重C1	群体置信度权重C2	适应度函数维数dim	群体数量	迭代次数	适应度函数
Wi->Wj	2->1	1->2	10	200	2000	Rastrigin(x)函数

W=0.9->0.2

W=0.7->0.2

分析：
  当其他参数不变，改变惯性权重w，算法前期收敛较快，当设置由较大的值开始下降，如0.9降到0.2，发现有多次能收敛到最优适应度0，但也会取到较大的最优适应度当设置由较小的值开始下降，最优适应度也逐渐减小，算法收敛到值
  惯性权重w，控制着前一速度对当前速度的影响。当惯性权重w较大，算法不易陷入局部最优，取到全局最优的几率较大，当惯性权重w较小，有利于算法收敛，可以取到较小的最优适应度。经验证，参数W取由0.8下降到0.4最佳，此时最优适应度可收敛到较小值。

2、改变个体置信度权重C1

参数列表：

惯性权重w	个体置信度权重C1	群体置信度权重C2	适应度函数维数dim	群体数量	迭代次数	适应度函数
0.8->0.4	C1i->C1j	1->2	10	200	2000	Rastrigin(x)函数

分析：
  当其他参数不变,改变个体置信度权重C1，当C1由过大降到较小时，最优适应度收敛较快，陷入局部最优，当C1由较小降到过小时，算法同样收敛较快且找到的局部最优值较大。当取C1由1.2下降到0.5时，能得到较优的局部最优值。

3、改变群体置信度权重C2

参数列表：

惯性权重w	个体置信度权重C1	群体置信度权重C2	适应度函数维数dim	群体数量	迭代次数	适应度函数
0.8->0.4	1.2->0.5	C2i->C2j	10	200	2000	Rastrigin(x)函数

分析：
  当其他参数不变,改变群体置信度权重C2，当取过小值增大到较小值（0.5->2），由算法收敛曲线变化图，可知在迭代过程中，先较快收敛，后逐渐收敛到局部最优适应；当取较小值增大到过大值（1->5/2->5）,算法快速收敛到局部最优，且收敛结果的局部最优个体适应度常落在某几个值，说明增大群体置信度权重C2，算法过于偏向群体置信度，则算法陷入相对群体全局最优。经验证，C2取1->2,由1增大到2，能得到较优的局部最优值。

4、改变适应度函数维数dim

参数列表：

惯性权重w	个体置信度权重C1	群体置信度权重C2	群体数量	迭代次数	适应度函数
0.8->0.4	1.2->0.5	1->2	200	2000	Rastrigin(x)函数

分析：
  由曲线图与数据可得，当其他参数不变，随着适应度函数维数增大，算法收敛速度变慢，算法搜索能力下降，难以找到最优解。

5、不同群体规模sizepop

参数列表：

惯性权重w	个体置信度权重C1	群体置信度权重C2	适应度函数维数dim	迭代次数	适应度函数
0.8->0.4	1.2->0.5	1->2	10	2000	Rastrigin(x)函数

分析：
  当群体规模增大，则群体多样性增大，算法搜索规模增大，时间也相对变长，群体规模较小，易陷入局部最优，群体规模较大，算法全局搜索能力增大，不易陷入局部最优，利于找到较好的解。经验证，群体规模取200较好。

5、不同适应度函数

参数列表：

惯性权重w	个体置信度权重C1	群体置信度权重C2	适应度函数维数dim	迭代次数	群体规模
0.8->0.4	1.2->0.5	1->2	10	2000	200

分析：
  3个适应度函数运行时间差不多，Rastrigin(x)函数迭代达到个体最佳适应度的次数较大一些，说明适应度函数Rastrigin(x)函数全局搜索能力较强；Schaffer(x)函数收敛较快，迭代达到个体最佳适应度的次数较小，但迭代完成后总能达到为0的个体最佳适应度；Griewank(x)函数迭代完成后也总能达到为0的个体最佳适应度。

五、总结

• 惯性权重w，控制前一速度对当前速度的影响，较大时，有利于跳出局部最优值，较小时有利于算法收敛；
• 个体置信度权重C1和群体置信度权重C2可平衡算法的收敛性；
• 适应度函数维数dim影响算法搜索能力，适应度函数维数增大，算法收敛速度变慢，算法搜索能力下降，难以找到最优解；
• 群体规模sizepop决定算法的多样性，种群规模增大，算法搜索规模增大，收敛速度变慢，运行时间增大。

  粒子群优化算法是概率型的全局优化算法，有利于求解全局最优解。但和其它群智能算法一样，在优化过程中，种群的多样性和算法的收敛速度之间始终存在着矛盾。总之，对于粒子群优化算法加强局部搜索能力以找到全局最优值的同时也要保持种群的多样性，才能避免收敛速度过快从而过早陷入局部最优。