粒子群优化算法对BP神经网络优化 Matlab实现

1、粒子群优化算法

粒子群算法(particle swarm optimization，PSO)由Kennedy和Eberhart在1995年提出，该算法模拟鸟集群飞行觅食的行为，鸟之间通过集体的协作使群体达到最优目的，是一种基于 Swarm Inteligence的优化方法。同遗传算法类似，也是一种基于群体叠代的，但并没有遗传算法用的交叉以及变异，而是粒子在解空间追随最优的粒子进行搜索。PSO的优势在于简单容易实现同时又有深刻的智能背景，既适合科学研究，又特别适合工程应用，并且没有许多参数需要调整。

粒子群优化算法源于1987年Reynolds对鸟群社会系统boids的仿真研究，boids是一个CAS。在boids中，一群鸟在空中飞行，每个鸟遵守以下三条规则：
1）避免与相邻的鸟发生碰撞冲突；
2）尽量与自己周围的鸟在速度上保持协调和一致；
3）尽量试图向自己所认为的群体中靠近。
仅通过使用这三条规则，boids系统就出现非常逼真的群体聚集行为，鸟成群地在空中飞行，当遇到障碍时它们会分开绕行而过，随后又会重新形成群体。
Reynolds仅仅将其作为CAS的一个实例作仿真研究，而并未将它用于优化计算中 。
Kennedy和Eberhart在中加入了一个特定点，定义为食物，鸟根据周围鸟的觅食行为来寻找食物。他们的初衷是希望通过这种模型来模拟鸟群寻找食源的现象，然而实验结果却揭示这个仿真模型中蕴涵着很强的优化能力，尤其是在多维空间寻优中。
PSO中，每个优化问题的解都是搜索空间中的一只鸟。称之为“粒子(Particle)”。所有的粒子都有一个由被优化的函数决定的适应值，每个粒子还有一个速度决定他们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索.
PSO 初始化为一群随机粒子。然后通过叠代找到最优解。在每一次叠代中，粒子通过跟踪两个"极值"来更新自己。第一个就是粒子本身所找到的最优解。这个解叫做个体极值pBest. 另一个极值是整个种群目前找到的最优解。这个极值是全局极值gBest。另外,也可以不用整个种群而只是用其中一部分的邻居。
PSO算法数学表示如下：
设搜索空间为D维，总粒子数为n。第i个粒子位置表示为向量Xi=( xi1, xi2,…, xiD )；第i个粒子 “飞行”历史中的过去最优位置（即该位置对应解最优）为Pi=( pi1,pi2,…,piD )，其中第g个粒子的过去最优位置Pg为所有Pi ( i=1, …,n)中的最优；第i个粒子的位置变化率（速度）为向量Vi=(vi1, vi2,…, viD)。每个粒子的位置按如下公式进行变化（“飞行”）：

 

 

其中，C1,C2为正常数，称为加速因子；rand( )为[0，1]之间的随机数；w称惯性因子，w较大适于对解空间进行大范围探查(exploration)，w较小适于进行小范围开挖(exploitation)。第d（1≤d≤D）维的位置变化范围为[-XMAXd , XMAXd]，速度变化范围为[-VMAXd , VMAXd]，迭代中若位置和速度超过边界范围则取边界值。

粒子群初始位置和速度随机产生，然后按公式(1)(2)进行迭代，直至找到满意的解。目前，常用的粒子群算法将全体粒子群(Global)分成若干个有部分粒子重叠的相邻子群，每个粒子根据子群(Local)内历史最优Pl调整位置，即公式(2)中Pgd换为Pld。

 

2、粒子群优化BP神经网络

神经网络与粒子群算法的结合主要有两种方式：一是利用粒子群算法的全局搜索能力来优化神经网络的拓扑结构、连接权值和阈值，将粒子群算法良好的全局寻优能力与BP算法良好的局部寻优能力相结合，以提高神经网络的泛化能力和学习性能，从而改进神经网络的整体搜索效率；二是将神经网络嵌入到粒子群算法当中，利用神经网络良好的学习性能来改进粒子群算法的优化性能，以提高粒子群算法的收敛速度，减少计算的工作量。

3、粒子群算法生成BP神经网络的权值和阈值，Matlab编程

                            
function [W1,W2,B1,B2]=pso(HiddenNum,InDim,OutDim)

% clc;clear 
% tic;                              %程序运行计时
E0=0.001;                        %允许误差
MaxNum=100;                    %粒子最大迭代次数
narvs=InDim*HiddenNum+1+HiddenNum+HiddenNum;                         %目标函数的自变量个数=51
particlesize=30;                    %粒子群规模=200
c1=2;                            %每个粒子的个体学习因子，也称为加速常数
c2=2;                            %每个粒子的社会学习因子，也称为加速常数
w=0.6;                           %惯性因子
vmax=0.8;                        %粒子的最大飞翔速度
x=-5+10*rand(particlesize,narvs);     %粒子所在的位置
v=2*rand(particlesize,narvs);         %粒子的飞翔速度
%用inline定义适应度函数以便将子函数文件与主程序文件放在一起，
%目标函数是：y=1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))
%inline命令定义适应度函数如下：
fitness=inline('1/(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');
%inline定义的适应度函数会使程序运行速度大大降低
for i=1:particlesize%粒子群规模
    for j=1:narvs %目标函数自变量个数
        f(i)=fitness(x(i,j));
    end
end
personalbest_x=x;%30*3
personalbest_faval=f;%30*1 但是是
[globalbest_faval i]=min(personalbest_faval);
globalbest_x=personalbest_x(i,:);
k=1;
while k<=MaxNum
    for i=1:particlesize
        for j=1:narvs
            f(i)=fitness(x(i,j));
        end
        if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置
            personalbest_faval(i)=f(i);
            personalbest_x(i,:)=x(i,:);
        end
    end
    [globalbest_faval i]=min(personalbest_faval);
    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:narvs    %判断粒子的飞翔速度是否超过了最大飞翔速度
            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_faval)<E0,break,end
    k=k+1;
end
Value1=1/globalbest_faval-1; 
Value1=num2str(Value1);
% strcat指令可以实现字符的组合输出
disp(strcat('the maximum value','=',Value1));
%输出最大值所在的横坐标位置
Value2=globalbest_x; 
% % Value2=num2str(Value2);
% % disp(strcat('the corresponding coordinate','=',Value2));
% x=-2:0.1:3;
% y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);
% plot(x,y,'m-','linewidth',3);
% hold on;
% plot(globalbest_x,1/globalbest_faval-1,'kp','linewidth',4);
% legend('样本的输出值','网络的输出值');xlabel('X');ylabel('Y');grid on;toc;

% 赋值给 bp神经网络的权值和阈值
W1=zeros(HiddenNum,InDim);
% B1=zeros(HiddenNum,1); 
% W2=zeros(1,HiddenNum);           
% B2=zeros(1,OutDim); 

    biaozhi=1;
    for wi=1:HiddenNum
        W1(wi,:)=Value2(1,biaozhi:biaozhi+InDim-1);
        biaozhi=wi*InDim+1;
    end
    B1=Value2(1,biaozhi:biaozhi+HiddenNum-1);
    B1=B1';
    W2=Value2(1,biaozhi+HiddenNum:biaozhi+HiddenNum+HiddenNum-1);
    B2=Value2(1,biaozhi+HiddenNum+HiddenNum);
end