粒子群算法在神经网络非线性函数拟合中的应用

1、粒子群算法在神经网络非线性函数拟合中的应用一、本文研究和解决的问题在自动控制问题中，系统辨识的目的是为了建立被控对象的数学模型。多年来，控制领域对于复杂的非线性对象的辨识 一直未能很好的解决，神经网络所具有的非线性特性和学习能力 使其在系统辨识方面有很大的潜力。为解决具有复杂的非线性、不确定性和不确知对象的辨识问题开辟了一条有效的途径。基于神经网络的系统辨识是以神经网络作为被辨识对象的模型，利用其非线性特性，可建立非线性系统的静态或动态模型。理论上，多层前馈神经网络能够以任意精度逼近任意非线性映射。 但传统神经网络学习算法中存在的收敛速度慢、容易陷入局部最优等缺点，于是设计了基于标准粒子群算法

2、 的神经网络非线性函数拟合系统。 二、传统的BP神经网络BP 神经网络即采用误差反向传播算法的网络，是一种至今仍然最为流行的前馈型神经网络模型。BP 神经网络有很强的非线性映射能力，它能学习和存贮大量输入-输出模式映射关系，而无需事先了解描述这种映射关系的数学方程。只要能提供足够多的样本模式对供给网络进行学习训练，它便能完成由 n 维输入空间到 m 维输出空间的非线性映射。BP 学习算法属于误差修正型学习，其关键在于根据误差修正输出层和隐含层的连接权值。其学习的基本实现方法是基于最小平方误差准则和梯度下降优化方法来确定权值调整法则。BP网络建模特点： 非线性映照能力：神经网络能以任意精度逼近任

3、何非线性连续函数。在建模过程中的许多问题正是具有高度的非线性。并行分布处理方式：在神经网络中信息是分布储存和并行处理的，这使它具有很强的容错性和很快的处理速度。自学习和自适应能力：神经网络在训练时，能从输入、输出的数据中提取出规律性的知识，记忆于网络的权值中，并具有泛化能力，即将这组权值应用于一般情形的能力。神经网络的学习也可以在线进行。数据融合的能力：神经网络可以同时处理定量信息和定性信息，因此它可以利用传统的工程技术（数值运算）和人工智能技术（符号处理）。多变量系统：神经网络的输入和输出变量的数目是任意的，对单变量系统与多变量系统提供了一种通用的描述方式，不必考虑各子系统间的解耦问题。三、

4、解决问题的思想与方法针对传统神经网络学习算法中存在的收敛速度慢、容易陷入局部最优等缺点，设计了基于标准粒子群算法的神经网络非线性函数拟合系统，将神经网络中的权值看作一个粒子，通过粒子之间的竞争与合作以完成网络的学习过程，仿真结果表明，基于BP的神经网络学习算法在收敛速度、辨识精度等方面要优于传统的BP神经网络。粒子群优化算法（PSO，Particle Swarm Optimization）是计算智能领域，除了蚁群算法、鱼群算法之外的一种群体智能的优化算法。PSO 算法源于对鸟类捕食行为的研究，鸟类捕食时，每只鸟找到食物最简单有效的方法就是搜寻当前距离食物最近的鸟的周围区域。PSO 算法首先在可

5、解空间中初始化一群粒子，每个粒子都代表问题的一个潜在解，用位置、速度和适应度值三项指标表示该粒子特征。适应度值由适应度函数计算得到，其值的好坏表示粒子的优劣。粒子的速度决定了粒子移动的方向和距离，速度随自身及其他粒子的移动经验进行动态调整，从而实现个体在可解空间中的寻优。粒子在解空间中运动，通过跟踪个体最优值 Pbest 和群体最优值Gbest 更新个体位置，个体最优值 Pbest 是指个体所经历位置中计算得到的适应度值最好的位置，群体最优值 Gbest 是指粒子群中所有粒子搜索到的适应度最好的位置。粒子每更新一次位置，就计算一次适应度值，并且通过比较新粒子的适应度值和个体最优值、群体最优值的

6、适应度值更新 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 , , pg

7、D )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

8、:max_epoch E=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

9、);Wb=Wb+lr.*deltai;endEz(epoch)=sqrt(E/100);if Ez(epoch)<Emin break;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);end ok(i)=NETk; sum=sum+(y(i)-ok(i)2;%输出总误差；e

10、nds=sqrt(sum/100);计算函数：function cs,wc=jsdyb(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=zjsdyb(lr)q=4,5,7,8,10;format short gA=;for zk=1:5cs,wc=jsdyb(lr,q(zk);B=cs,wc;A=A;B;end图形显示函数：function txdyb(lr,q)%计算测试输出；Emin=0.1;b=1;epoch,s,Wki,Wij

11、,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);end ok(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);

12、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('outp

13、ut 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; %训练样本输出%

14、选取测试样本for i=1:100 x=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,maxAlltestOut= premnmx(AlltestIn,AlltestOut); %测试样本AllSamInn,minAl

15、lSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut= 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; %输出层神经

16、元个数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:itmax W(iter)=wmax-(wmax-wmin)/itmax)*iter; end a=-1; b=1;m=-1;n=1;global N; % 微粒个数N=30;global D; % 每个微粒的维数D=(indim+1)*hiddennum+(hid

17、dennum+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

18、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

19、 p=1:N G(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:N for di=1:D if V(ni,di,2)>vmax V(ni,di,2)=vmax; else if V(ni,

20、di,2)<-vmax V(ni,di,2)=-vmax; else V(ni,di,2)=V(ni,di,2); end endend X(:,:,2)=X(:,:,1)+V(:,:,2); %更新位置%disp('执行到这里')%2.2 第2次到最后一次迭代for j=2:itmax h=j; disp('迭代次数，当前代全局最佳适应值，本代以前所有代中的全局最佳适应值') disp(j-1) disp(B(1,1,j-1) %j-1代全局最优适应值 disp(bestminimum(j-1) %j-1代以前所有代中的全局最优适应值 disp('

21、;*') 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，代数赋给GI bestminimum(j)=GC; %所有代的最优适应值赋给j代的bestminimum %

22、判断是否符合条件 if GC<=minerr Gpos=gbest(1,:,GI); %若满足均方误差条件，记录最优位置，停止迭代 break end if j>=itmax break %超过最大迭代次数时，退出end %计算历史全局最优位置 if B(1,1,j)<GC gbest(1,:,j)=gbest(1,:,j); else gbest(1,:,j)=gbest(1,:,GI); end for p=1:N G(p,:,j)=gbest(1,:,j); end %计算各微粒历史最优位置 for i=1:N; C,I=min(L(i,1,:); %计算每个微粒的历史

23、最优适应值，赋给C，代数赋给I if L(i,1,j)<=C pbest(i,:,j)=X(i,:,j); else pbest(i,:,j)=X(i,:,I); end end V(:,:,j+1)=W(j)*V(:,:,j)+c1*rand*(pbest(:,:,j)-X(:,:,j)+c2*rand*(G(:,:,j)-X(:,:,j); for ni=1:N for di=1:D if V(ni,di,j+1)>vmax V(ni,di,j+1)=vmax; else if V(ni,di,j+1)<-vmax V(ni,di,j+1)=-vmax; else V(n

24、i,di,j+1)=V(ni,di,j+1); end end end X(:,:,j+1)=X(:,:,j)+V(:,:,j+1); %2.3 将最优微粒（即最优权值和阈值）赋给神经网络 if j=itmax Gpos=gbest(1,:,GI); end disp('要显示Gpos的值') disp(Gpos) wi=Gpos(1:hiddennum); %输入层-隐藏层权值 wl=Gpos(hiddennum+1:2*hiddennum); %隐藏层-输出层权值 b1=Gpos(2*hiddennum+1:3*hiddennum); %输入层-隐藏层阈值 b2=Gpos(

25、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); end outinput=0; %used to calculate the value of output in outlayer for t=1

26、:hiddennum outinput=outinput+wl(t)*hidoutput(t); %输出层只有一个神经元时的情况 end outVal(k)=purelin(outinput-b2); %输出层的输出值endsubplot(2,1,1) %/调用窗口句柄AlltestIn,AlltestOut=postmnmx(testIn,minAlltestIn,maxAlltestIn,testOut,minAlltestOut,maxAlltestOut); %反归一化ResVal=postmnmx(outVal,minAlltestOut,maxAlltestOut);trainEr

27、ror=abs(ResVal-AlltestOut); %测试误差for k=1:v SquareE(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