粒子群算法在神经网络非线性函数拟合中的应用一、本文研究和解决的问题在自动控制问题中，系统辨识的目的是为了建立被控对象的数学模型。
多年来，控制领域对于复杂的非线性对象的辨识一直未能很好的解决，神经网络所具有的非线性特性和学习能力使其在系统辨识方面有很大的潜力。
为解决具有复杂的非线性、不确定性和不确知对象的辨识问题开辟了一条有效的途径。
基于神经网络的系统辨识是以神经网络作为被辨识对象的模型，利用其非线性特性，可建立非线性系统的静态或动态模型。
理论上，多层前馈神经网络能够以任意精度逼近任意非线性映射。
但传统神经网络学习算法中存在的收敛速度慢、容易陷入局部最优等缺点，于是设计了基于标准粒子群算法的神经网络非线性函数拟合系统。
二、传统的BP神经网络BP 神经网络即采用误差反向传播算法的网络，是一种至今仍然最为流行的前馈型神经网络模型。
BP 神经网络有很强的非线性映射能力，它能学习和存贮大量输入-输出模式映射关系，而无需事先了解描述这种映射关系的数学方程。
只要能提供足够多的样本模式对供给网络进行学习训练，它便能完成由n 维输入空间到m 维输出空间的非线性映射。
BP 学习算法属于误差修正型学习，其关键在于根据误差修正输出层和隐含层的连接权值。
其学习的基本实现方法是基于最小平方误差准则和梯度下降优化方法来确定权值调整法则。
BP网络建模特点：非线性映照能力：神经网络能以任意精度逼近任何非线性连续函数。
在建模过程中的许多问题正是具有高度的非线性。
并行分布处理方式：在神经网络中信息是分布储存和并行处理的，这使它具有很强的容错性和很快的处理速度。
自学习和自适应能力：神经网络在训练时，能从输入、输出的数据中提取出规律性的知识，记忆于网络的权值中，并具有泛化能力，即将这组权值应用于一般情形的能力。
神经网络的学习也可以在线进行。
数据融合的能力：神经网络可以同时处理定量信息和定性信息，因此它可以利用传统的工程技术（数值运算）和人工智能技术（符号处理）。
多变量系统：神经网络的输入和输出变量的数目是任意的，对单变量系统与多变量系统提供了一种通用的描述方式，不必考虑各子系统间的解耦问题。
三、解决问题的思想与方法针对传统神经网络学习算法中存在的收敛速度慢、容易陷入局部最优等缺点，设计了基于标准粒子群算法的神经网络非线性函数拟合系统，将神经网络中的权值看作一个粒子，通过粒子之间的竞争与合作以完成网络的学习过程，仿真结果表明，基于BP的神经网络学习算法在收敛速度、辨识精度等方面要优于传统的BP神经网络。
粒子群优化算法（PSO，Particle Swarm Optimization）是计算智能领域，除了蚁群算法、鱼群算法之外的一种群体智能的优化算法。
PSO 算法源于对鸟类捕食行为的研究，鸟类捕食时，每只鸟找到食物最简单有效的方法就是搜寻当前距离食物最近的鸟的周围区域。
PSO 算法首先在可解空间中初始化一群粒子，每个粒子都代表问题的一个潜在解，用位置、速度和适应度值三项指标表示该粒子特征。
适应度值由适应度函数计算得到，其值的好坏表示粒子的优劣。
粒子的速度决定了粒子移动的方向和距离，速度随自身及其他粒子的移动经验进行动态调整，从而实现个体在可解空间中的寻优。
粒子在解空间中运动，通过跟踪个体最优值Pbest 和群体最优值Gbest 更新个体位置，个体最优值Pbest 是指个体所经历位置中计算得到的适应度值最好的位置，群体最优值Gbest 是指粒子群中所有粒子搜索到的适应度最好的位置。
粒子每更新一次位置，就计算一次适应度值，并且通过比较新粒子的适应度值和个体最优值、群体最优值的适应度值更新Pbest 和Gbest 的位置。
粒子位置和速度的调整是粒子群算法的关键。
假设在一个D 维的搜索空间中，由n 个粒子组成的种群X ( X1, X 2 ,……, X n ) ，其中第i 个例子表示为一个D 维的向量X i (xi1, xi 2 ,……, xiD )T ，代表第i 个粒子在D维搜索空间中的位置，亦代表问题的一个潜在解。
根据目标函数即可计算出每个粒子位置Xi 对应的适应度值。
第i 个粒子的速度为Vi (vi1, vi 2 ,……, viD )T ，其个体极值为P ( pi1, pi 2 ,……, piD )T ，种群的全局极值为P( pg1, pg 2 ,……, pgD )T 。
四、实验仿真结果五、程序基本ＢＰ网络函数逼近程序function [epoch,s,Wki,Wij,Wb,Ez]=dyb(lr,Emin,q)%初始化；%lr 学习效率；Emin为期望误差最小值；q为隐含层节点数；b=1;sum=0;Ez=[];max_epoch=30000;%max_epoch训练的最大次数；%提供训练集和目标值；x=8.*rand(1,100)-4;y=1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100);%初始化Wki，Wij；Wij=rand(1,q);Wki=rand(1,q);Wb=rand(1,q);for epoch=1:max_epochE=0;m=1;oi=0;ok=0;%置隐含层和输出层各神经元输出初值为零；for m=1:100%计算隐含层各神经元输出；NETi=x(m)*Wij+b*Wb;for t=1:qoi(t)=1/(1+exp(-NETi(t)));end%计算输出层各神经元输出；NETk=Wki*oi';ok=NETk;%计算误差；E=E+(y(m)-ok)^2;%调整输出层加权系数；deltak=y(m)-ok;Wki=Wki+lr*deltak*oi;%调整隐含层加权系数；deltai=oi.*(1-oi).*(deltak*Wki);Wij=Wij+lr.*deltai.*x(m);Wb=Wb+lr.*deltai;endEz(epoch)=sqrt(E/100);if Ez(epoch)<Eminbreak;endend%计算测试输出；x=linspace(-4,4,100);%给定输入：y=1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100); for i=1:100NETi=x(i).*Wij+b*Wb;NETk=0;for t=1:qoi(t)=1/(1+exp(-NETi(t)));NETk=NETk+Wki(t)*oi(t);endok(i)=NETk;sum=sum+(y(i)-ok(i))^2;%输出总误差；ends=sqrt(sum/100);计算函数：function [cs,wc]=js dyb(lr,q)Emin=0.1;s1=0;s2=0;for k=1:5[x1,x2]=dyb(lr,Emin,q);s1=s1+x1;s2=s2+x2;endcs=s1/5;wc=s2/5;function A=zjs dyb(lr)q=[4,5,7,8,10];format short gA=[];for zk=1:5[cs,wc]=js dyb(lr,q(zk));B=[cs,wc];A=[A;B];end图形显示函数：function tx dyb(lr,q)%计算测试输出；Emin=0.1;b=1;[epoch,s,Wki,Wij,Wb,Ez]=dyb(lr,Emin,q)x=linspace(-4,4,100);%给定输入：y=1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100); for i=1:100NETi=x(i).*Wij+b*Wb;NETk=0;for t=1:qoi(t)=1/(1+exp(-NETi(t)));NETk=NETk+Wki(t)*oi(t);endok(i)=NETk;end%显示图形；figureplot(x,ok,'r')hold ony=1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100); plot(x,y,'b')title('Hermit多项式曲线与BP网络输出曲线')legend('BP曲线','Hermit曲线')hold offfigureplot(x,ok,'or')hold onx=8.*rand(1,100)-4;y=1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100); plot(x,y,'*k')title('训练样本与测试样本')xlabel('input x')ylabel('output y')legend('测试样本','训练样本')figureplot([1:length(Ez)],Ez)title('收敛曲线')clcPSO优化神经网络程序clcclear all%一、初始化部分%1.1 预处理样本数据% 选取训练样本(x,y)for i=1:100x=8.*rand(1,100)-4;y=1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100);% 待逼近函数endAllSamIn=linspace(-4,4,100); %训练样本输入AllSamOut=y; %训练样本输出%选取测试样本for i=1:100x=8.*rand(1,100)-4; %测试样本输入ytest =1.1.*(1-x+2.*x.^2).*exp(-x.^2/2)+0.1*rand(1,100); %测试样本输出endAlltestIn= linspace(-4,4,100);AlltestOut=ytest;%归一化训练样本，测试样本[AlltestInn,minAlltestIn,maxAlltestIn,AlltestOutn,minAlltestOut,maxAl ltestOut]= premnmx(AlltestIn,AlltestOut); %测试样本[AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamO ut]= premnmx(AllSamIn,AllSamOut); %训练样本testIn=AlltestInn;testOut=AlltestOutn;global Ptrain;Ptrain = AllSamInn;global Ttrain;Ttrain = AllSamOutn;%1.2 设置神经网络参数global indim; %输入层神经元个数indim=1;global hiddennum; %隐藏层神经元个数hiddennum=3;global outdim; %输出层神经元个数outdim=1;global Gpos;%1.3 设置微粒群参数vmax=0.5; % 速度上限minerr=1e-7; % 目标误差wmax=0.95;wmin=0.25;global itmax; % 最大迭代次数itmax=200;c1=1.5;c2=1.5;%权值随迭代次数线性递减以保证收敛for iter=1:itmaxW(iter)=wmax-((wmax-wmin)/itmax)*iter; enda=-1;b=1;m=-1;n=1;global N; % 微粒个数N=30;global D; % 每个微粒的维数D=(indim+1)*hiddennum+(hiddennum+1)*outdim; %所有权值和阈值% 初始化微粒位置rand('state',sum(100*clock)); %产生和时间相关的随机数global X;X=a+(b-a)*rand(N,D,1); %X的值在a 和b之间%初始化微粒速度V=m+(n-m)*rand(N,D,1); %V的值在m和n之间%二、微粒群更新迭代部分%global net;net=newff(minmax(Ptrain),[hiddennum,outdim],{'tansig','purelin'}); global gbest; %全局最优位置global pbest; %局部最优位置%2.1第一次迭代fitness=fitcal(X,indim,hiddennum,outdim,D,Ptrain,Ttrain); %计算适应值[C,I]=min(fitness(:,1,1)); %第一代，返回微粒群中最小适应值给C，该微粒的序号给IL(:,1,1)=fitness(:,1,1); %第一代，每个微粒的适应值B(1,1,1)=C; %第一代，全局最优适应值（B存储当前代最优适应值）bestminimum(1)=C; % bestminimum存储所有代中的全局最小适应值gbest(1,:,1)=X(I,:,1); %第一代，全局最优的微粒位置for p=1:NG(p,:,1)=gbest(1,:,1); %G便于速度更新运算（函数格式统一）endGpos=gbest(1,:,1);for i=1:N;pbest(i,:,1)=X(i,:,1); %因为是第一代，当前位置即为历史最优位置endV(:,:,2)=W(1)*V(:,:,1)+c1*rand*(pbest(:,:,1)-X(:,:,1))+c2*rand*(G(:,: ,1)-X(:,:,1)); % 更新速度% 判断速度是否越界for ni=1:Nfor di=1:Dif V(ni,di,2)>vmaxV(ni,di,2)=vmax;else if V(ni,di,2)<-vmaxV(ni,di,2)=-vmax;elseV(ni,di,2)=V(ni,di,2);endendendX(:,:,2)=X(:,:,1)+V(:,:,2); %更新位置%disp('执行到这里')%2.2 第2次到最后一次迭代for j=2:itmaxh=j;disp('迭代次数，当前代全局最佳适应值，本代以前所有代中的全局最佳适应值')disp(j-1)disp(B(1,1,j-1)) %j-1代全局最优适应值disp(bestminimum(j-1)) %j-1代以前所有代中的全局最优适应值disp('******************************')fitness=fitcal(X,indim,hiddennum,outdim,D,Ptrain,Ttrain);[C,I]=min(fitness(:,1,j)); %第j代的最优适应值和最优微粒序号L(:,1,j)=fitness(:,1,j); %第j代每个微粒的适应值B(1,1,j)=C; %第j代全局最优适应值gbest(1,:,j)=X(I,:,j); %第j代全局最优微粒的位置[GC,GI]=min(B(1,1,:)); %所有代的全局最优适应值赋给GC，代数赋给GIbestminimum(j)=GC; %所有代的最优适应值赋给j代的bestminimum% 判断是否符合条件if GC<=minerrGpos=gbest(1,:,GI); %若满足均方误差条件，记录最优位置，停止迭代breakendif j>=itmaxbreak %超过最大迭代次数时，退出end%计算历史全局最优位置if B(1,1,j)<GCgbest(1,:,j)=gbest(1,:,j);elsegbest(1,:,j)=gbest(1,:,GI);endfor p=1:NG(p,:,j)=gbest(1,:,j);end%计算各微粒历史最优位置for i=1:N;[C,I]=min(L(i,1,:)); %计算每个微粒的历史最优适应值，赋给C，代数赋给Iif L(i,1,j)<=Cpbest(i,:,j)=X(i,:,j);elsepbest(i,:,j)=X(i,:,I);endendV(:,:,j+1)=W(j)*V(:,:,j)+c1*rand*(pbest(:,:,j)-X(:,:,j))+c2*rand*(G(: ,:,j)-X(:,:,j));for ni=1:Nfor di=1:Dif V(ni,di,j+1)>vmaxV(ni,di,j+1)=vmax;else if V(ni,di,j+1)<-vmaxV(ni,di,j+1)=-vmax;elseV(ni,di,j+1)=V(ni,di,j+1);endendendX(:,:,j+1)=X(:,:,j)+V(:,:,j+1);%2.3 将最优微粒（即最优权值和阈值）赋给神经网络if j=itmaxGpos=gbest(1,:,GI);enddisp('要显示Gpos的值')disp(Gpos)wi=Gpos(1:hiddennum); %输入层-隐藏层权值wl=Gpos(hiddennum+1:2*hiddennum); %隐藏层-输出层权值b1=Gpos(2*hiddennum+1:3*hiddennum); %输入层-隐藏层阈值b2=Gpos(3*hiddennum+1:3*hiddennum+outdim); %隐藏层-输出层阈值end%三、神经网络训练部分%******************************************************************** ********[w,v]=size(testIn); %w 返回行数，v返回列数for k=1:v % v是测试样本的个数for t=1:hiddennum %计算隐藏层每个神经元的输入，输出hidinput=0;hidinput=wi(t)*testIn(k)-b1(t);hidoutput(t)=tansig(hidinput);endoutinput=0; %used to calculate the value of output in outlayerfor t=1:hiddennumoutinput=outinput+wl(t)*hidoutput(t); %输出层只有一个神经元时的情况endoutVal(k)=purelin(outinput-b2); %输出层的输出值endsubplot(2,1,1) %//调用窗口句柄[AlltestIn,AlltestOut]=postmnmx(testIn,minAlltestIn,maxAlltestIn,testOut,minAlltestOut,maxAl ltestOut); %反归一化[ResVal]=postmnmx(outVal,minAlltestOut,maxAlltestOut);trainError=abs(ResVal-AlltestOut); %测试误差for k=1:vSquareE(k)=(trainError(k)*trainError(k))/2; %v个样本的误差数组endplot(AlltestIn,SquareE)ylabel('Error')subplot(2,1,2)j=1:1:h;plot(j,bestminimum(j))set(gca,'XLim',[1 100000]);set(gca,'XMinorTick','on');set(gca,'XTick',[1 10 100 1000 10000 100000]);set(gca,'YLim',[0.000001 1]);set(gca,'YMinorTick','on');set(gca,'YTick',[0.000001 0.00001 0.0001 0.001 0.01 0.1 1]);set(gca,'yscale','log','xscale','log')ylabel('training error')xlabel('Iteration Number')hold on%适应度函数部分function fitval = fitcal(X,indim,hiddennum,outdim,D,Ptrain,Ttrain) %三维矩阵：x 微粒数（X的行数）；y 微粒维数（X的列数）；z 代数（X的层数）[x,y,z]=size(X);[w,v]=size(Ptrain); %二维矩阵：w 训练样本维数，这里为1；v 训练样本个数for i=1:x %x代表粒子数量，z代表代数wi=X(i,1:hiddennum,z);wl=X(i,1*hiddennum+1:2*hiddennum,z);b1=X(i,2*hiddennum+1:3*hiddennum,z);b2=X(i,3*hiddennum+1:3*hiddennum+outdim,z);error=0;for k=1:v %训练样本总数for t=1:hiddennumhidinput=0;hidinput=wi(t)*Ptrain(k)-b1(t);hidoutput(t)=tansig(hidinput);endoutinput=0;for t=1:hiddennumoutinput=outinput+wl(t)*hidoutput(t);endoutval(k)=purelin(outinput-b2);errval(k)=Ttrain(k)-outval(k); %绝对误差error=error+errval(k)*errval(k); %v个样本的误差平方求和endfitval(i,1,z)=error/v; %均方和,返回值是第i个微粒第z代的误差end。