【优化预测】基于matlab天牛须算法优化BP神经网络预测【含Matlab源码 1318期】

一、天牛须搜索算法简介

1 天牛须搜索算法定义
天牛须搜索（Beetle Antennae Search-BAS），也叫甲壳虫须搜索，是2017年提出的一种高效的智能优化算法。类似于遗传算法、粒子群算法、模拟退火等智能优化算法，天牛须搜索不需要知道函数的具体形式，不要虚梯度信息，就可以实现高效寻优。相比于粒子群算法，天牛须搜索只只要一个个体，即一个天牛，运算量大大降低。

2 原理及代码实现
2.1 仿生原理
天牛须搜索时受到天牛觅食原理启发而开发的算法。
生物原理：当天牛觅食时，天牛并不知道食物在哪，而是根据食物气味的强弱来觅食。天牛有俩只长触角，如果左边触角收到的气味强度比右边大，那下一步天牛就往左飞，否则就往右飞。根据这一简单原理天牛就可以有效找到食物。
天牛须搜索得来的启发：食物的气味就相当于一个函数，这个函数在三维空间每个点值都不同，天牛两个须可以采集自身附近两点的气味值，天牛的目的是找到全局气味值最大的点。仿照天牛的行为，我们就可以高效的进行函数寻优。

2.2 算法
天牛在三维空间运动，而天牛须搜索需要对任意维函数都有效才可以。因而，天牛须搜索是对天牛生物行为在任意维空间的推广。采用如下的简化模型假设描述天牛：
天牛左右两须位于质心两边。
天牛步长step与两须之间距离d0的比是个固定常数，即step=c*d0，其中c是常数。即，大天牛（两须距离长）走大步，小天牛走小步。
天牛飞到下一步后，头的朝向是随机的。

2.3 建模：（n维空间函数f最小化）
第一步：对一个n维空间的优化问题，我们用xl表示左须坐标，xr表示右须坐标，x表示质心坐标，用d0表示两须之间的距离。根据假设3，天牛头朝向任意，因而从天牛右须指向左须的向量的朝向也是任意的，所以可以产生一个随机向量dir=rands(n,1)来表示它。对此归一化：dir=dir/norm(dir)；我们这样可以得到xl-xr=d0dir；显然，xl，xr还可以表示成质心的表达式；xl=x+d0dir/2；xr=x-d0dir/2。
第二步：对于待优化函数f，求取左右两须的值：felft=f(xl)；fright=f(xr)；判断两个值大小，如果fleft<fright，为了探寻f的最小值，则天牛向着左须方向行进距离step，即x=x+stepnormal(xl-xr)；如果fleft>fright，为了探寻f的最小值，则天牛向着右须方向行进距离step，即x=x-stepnormal(xl-xr)；如以上两种情况可以采用符号函数sign统一写成：x=x-stepnormal(xl-xr)sign(fleft-fright)=x-stepdir*sign(fleft-fright)。
(注：其中normal是归一化函数)

循环迭代：
dir=rands(n,1)；dir=dir/norm(dir)；%须的方向
xl=x+d0dir/2；xr=x-d0dir/2；%须的坐标
felft=f(xl)；fright=f(xr)；%须的气味强度
x=x-stepdirsign(fleft-fright)。%下一步位置
关于步长:
两种推荐：
每步迭代中采用step=etastep，其中eta在0，1之间靠近1，通常可取eta=0.95；
引入新变量temp和最终分辨率step0，temp=etatemp，step=temp+step0.
关于初始步长：初始步长可以尽可能大，最好与自变量最大长度相当。

二、BP神经网络简介

1 BP神经网络概述
BP（Back Propagation）神经网络是1986年由Rumelhart和McCelland为首的科研小组提出，参见他们发表在Nature上的论文 Learning representations by back-propagating errors 。
BP神经网络是一种按误差逆传播算法训练的多层前馈网络，是目前应用最广泛的神经网络模型之一。BP网络能学习和存贮大量的输入-输出模式映射关系，而无需事前揭示描述这种映射关系的数学方程。它的学习规则是使用最速下降法，通过反向传播来不断调整网络的权值和阈值，使网络的误差平方和最小。

2 BP算法的基本思想
上一次我们说到，多层感知器在如何获取隐层的权值的问题上遇到了瓶颈。既然我们无法直接得到隐层的权值，能否先通过输出层得到输出结果和期望输出的误差来间接调整隐层的权值呢？BP算法就是采用这样的思想设计出来的算法，它的基本思想是,学习过程由信号的正向传播与误差的反向传播两个过程组成。
正向传播时，输入样本从输入层传入,经各隐层逐层处理后,传向输出层。若输出层的实际输出与期望的输出(教师信号)不符,则转入误差的反向传播阶段。
反向传播时，将输出以某种形式通过隐层向输入层逐层反传,并将误差分摊给各层的所有单元,从而获得各层单元的误差信号,此误差信号即作为修正各单元权值的依据。这两个过程的具体流程会在后文介绍。

BP算法的信号流向图如下图所示

3 BP网络特性分析——BP三要素
我们分析一个ANN时，通常都是从它的三要素入手，即
1)网络拓扑结构；
2)传递函数；
3)学习算法。

每一个要素的特性加起来就决定了这个ANN的功能特性。所以，我们也从这三要素入手对BP网络的研究。

3.1 BP网络的拓扑结构
上一次已经说了，BP网络实际上就是多层感知器，因此它的拓扑结构和多层感知器的拓扑结构相同。由于单隐层（三层）感知器已经能够解决简单的非线性问题，因此应用最为普遍。三层感知器的拓扑结构如下图所示。
一个最简单的三层BP：

3.2 BP网络的传递函数
BP网络采用的传递函数是非线性变换函数——Sigmoid函数（又称S函数）。其特点是函数本身及其导数都是连续的，因而在处理上十分方便。为什么要选择这个函数，等下在介绍BP网络的学习算法的时候会进行进一步的介绍。
单极性S型函数曲线如下图所示。

双极性S型函数曲线如下图所示。

3.3 BP网络的学习算法
BP网络的学习算法就是BP算法，又叫 δ 算法（在ANN的学习过程中我们会发现不少具有多个名称的术语），以三层感知器为例，当网络输出与期望输出不等时，存在输出误差 E ，定义如下

下面我们会介绍BP网络的学习训练的具体过程。

4 BP网络的训练分解
训练一个BP神经网络，实际上就是调整网络的权重和偏置这两个参数，BP神经网络的训练过程分两部分：

前向传输，逐层波浪式的传递输出值；
逆向反馈，反向逐层调整权重和偏置；
我们先来看前向传输。
前向传输（Feed-Forward前向反馈）
在训练网络之前，我们需要随机初始化权重和偏置，对每一个权重取[ − 1 , 1 ] [-1,1][−1,1]的一个随机实数，每一个偏置取[ 0 , 1 ] [0,1][0,1]的一个随机实数，之后就开始进行前向传输。

神经网络的训练是由多趟迭代完成的，每一趟迭代都使用训练集的所有记录，而每一次训练网络只使用一条记录，抽象的描述如下：

while 终止条件未满足：for record:dataset:trainModel(record)

4.1 逆向反馈（Backpropagation）

4.2 训练终止条件
每一轮训练都使用数据集的所有记录，但什么时候停止，停止条件有下面两种：
设置最大迭代次数，比如使用数据集迭代100次后停止训练
计算训练集在网络上的预测准确率，达到一定门限值后停止训练

5 BP网络运行的具体流程
5.1 网络结构
输入层有n nn个神经元，隐含层有p pp个神经元，输出层有q qq个神经元。

5.2 变量定义

第九步：判断模型合理性
判断网络误差是否满足要求。
当误差达到预设精度或者学习次数大于设计的最大次数，则结束算法。
否则，选取下一个学习样本以及对应的输出期望，返回第三部，进入下一轮学习。

6 BP网络的设计
在进行BP网络的设计是，一般应从网络的层数、每层中的神经元个数和激活函数、初始值以及学习速率等几个方面来进行考虑，下面是一些选取的原则。

6.1 网络的层数
理论已经证明，具有偏差和至少一个S型隐层加上一个线性输出层的网络，能够逼近任何有理函数，增加层数可以进一步降低误差，提高精度，但同时也是网络复杂化。另外不能用仅具有非线性激活函数的单层网络来解决问题，因为能用单层网络解决的问题，用自适应线性网络也一定能解决，而且自适应线性网络的运算速度更快，而对于只能用非线性函数解决的问题，单层精度又不够高，也只有增加层数才能达到期望的结果。

6.2 隐层神经元的个数
网络训练精度的提高，可以通过采用一个隐含层，而增加其神经元个数的方法来获得，这在结构实现上要比增加网络层数简单得多。一般而言，我们用精度和训练网络的时间来恒量一个神经网络设计的好坏：
（1）神经元数太少时，网络不能很好的学习，训练迭代的次数也比较多，训练精度也不高。
（2）神经元数太多时，网络的功能越强大，精确度也更高，训练迭代的次数也大，可能会出现过拟合(over fitting)现象。
由此，我们得到神经网络隐层神经元个数的选取原则是：在能够解决问题的前提下，再加上一两个神经元，以加快误差下降速度即可。

6.3 初始权值的选取
一般初始权值是取值在(−1,1)之间的随机数。另外威得罗等人在分析了两层网络是如何对一个函数进行训练后，提出选择初始权值量级为s√r的策略，其中r为输入个数，s为第一层神经元个数。

6.4 学习速率
学习速率一般选取为0.01−0.8，大的学习速率可能导致系统的不稳定，但小的学习速率导致收敛太慢，需要较长的训练时间。对于较复杂的网络，在误差曲面的不同位置可能需要不同的学习速率，为了减少寻找学习速率的训练次数及时间，比较合适的方法是采用变化的自适应学习速率，使网络在不同的阶段设置不同大小的学习速率。

6.5 期望误差的选取
在设计网络的过程中，期望误差值也应当通过对比训练后确定一个合适的值，这个合适的值是相对于所需要的隐层节点数来确定的。一般情况下，可以同时对两个不同的期望误差值的网络进行训练，最后通过综合因素来确定其中一个网络。

7 BP网络的局限性
BP网络具有以下的几个问题：
(1)需要较长的训练时间：这主要是由于学习速率太小所造成的，可采用变化的或自适应的学习速率来加以改进。
(2)完全不能训练：这主要表现在网络的麻痹上，通常为了避免这种情况的产生，一是选取较小的初始权值，而是采用较小的学习速率。
(3)局部最小值：这里采用的梯度下降法可能收敛到局部最小值，采用多层网络或较多的神经元，有可能得到更好的结果。

8 BP网络的改进
P算法改进的主要目标是加快训练速度，避免陷入局部极小值等，常见的改进方法有带动量因子算法、自适应学习速率、变化的学习速率以及作用函数后缩法等。动量因子法的基本思想是在反向传播的基础上，在每一个权值的变化上加上一项正比于前次权值变化的值，并根据反向传播法来产生新的权值变化。而自适应学习速率的方法则是针对一些特定的问题的。改变学习速率的方法的原则是，若连续几次迭代中，若目标函数对某个权倒数的符号相同，则这个权的学习速率增加，反之若符号相反则减小它的学习速率。而作用函数后缩法则是将作用函数进行平移，即加上一个常数。

三、部分源代码

%% 用天牛须算法来优化BP的权值和阈值,数据样本为测试数据，非论文实际数据，样本60个，其中每个样本具有401个特征值；NIR为样本的光谱数据，octane为60*1的辛烷值数据
% 1.0版本
%% 清空环境变量
clear all
close all
clc
tic
%% 加载数据
load spectra_data.mat
% 随机产生训练集和测试集
temp=randperm(size(NIR,1));
%训练集——50个样本
P=NIR(temp(1:50),:)';
T=octane(temp(1:50),:)';
%测试集——10个样本
P_test=NIR(temp(51:end),:)';
T_test=octane(temp(51:end),:)';
N=size(P_test,2);
M=size(P,2);%% 归一化
[P, ps_input] = mapminmax(P,0,1);%p_train归一化处理，范围为[0,1]，默认情况下为[-1,1]
P_test = mapminmax('apply',P_test,ps_input);%对P_test采用相同的映射
[T, ps_output] = mapminmax(T,0,1);
%%
inputnum=size(P,1);
outputnum=size(T,1);
hiddennum=9;%初始隐含层神经元个数
%% 创建网络
net=newff(P,T,hiddennum);
net.trainParam.epochs = 1000;
net.trainParam.goal = 1e-3;
net.trainParam.lr = 0.01;
%% 天牛须算法初始化
eta=0.8;
c=5;%步长与初始距离之间的关系
step=30;%初始步长
n=100;%迭代次数
k=inputnum*hiddennum+outputnum*hiddennum+hiddennum+outputnum;
x=rands(k,1);
bestX=x;
bestY=fitness(bestX,inputnum,hiddennum,outputnum,net,P,T);
fbest_store=bestY;
x_store=[0;x;bestY];
display(['0:','xbest=[',num2str(bestX'),'],fbest=',num2str(bestY)])
%% 迭代部分
for i=1:n
d0=step/c;dir=rands(k,1);dir=dir/(eps+norm(dir));xleft=x+dir*d0/2;fleft=fitness(xleft,inputnum,hiddennum,outputnum,net,P,T);xright=x-dir*d0/2;fright=fitness(xright,inputnum,hiddennum,outputnum,net,P,T);x=x-step*dir*sign(fleft-fright);y=fitness(x,inputnum,hiddennum,outputnum,net,P,T);if y<bestYbestX=x;bestY=y;endif y<0.001bestX=x;bestY=y;endx_store=cat(2,x_store,[i;x;y]);fbest_store=[fbest_store;bestY];step=step*eta;display([num2str(i),':xbest=[',num2str(bestX'),'],fbest=',num2str(bestY)])
end
M=size(P,2);
w1= bestX(1:inputnum*hiddennum);
B1= bestX(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2= bestX(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2= bestX(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);
%% 网络权值赋值
%% 使用优化后的权值和阈值测试结果
%% 使用优化后的权值和阈值
inputnum=size(P,1);%输入层神经元个数
outputnum=size(T,1);%输出层神经元个数
N=size(P_test,2);
M=size(P,2);
%% 新建BP
net=newff(P,T,9);
%% 设置网络参数：训练次数1000，训练目标0.001，学习速率00.1
net.trainParam.epochs =3000;
net.trainParam.goal = 1e-6;
net.trainParam.lr = 0.01;
%% BP初始权值和阈值
w1num=inputnum*hiddennum;%输入层到隐含层的权值个数
w2num=outputnum*hiddennum;%隐含层到输入层的权值个数
w1=bestX(1:w1num);%初始输入层到隐含层的权值
B1=bestX(w1num+1:w1num+hiddennum);
w2=bestX(w1num+hiddennum+1:w1num+hiddennum+w2num);%初始隐含层到输出层的权值
B2=bestX(w1num+hiddennum+w2num+1:w1num+hiddennum+w2num+outputnum);%输出层阈值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%% 训练网络
net=train(net,P,T);
%% 测试网络
t_sim_P= sim(net,P);
t_sim_P_test= sim(net,P_test);
%% 反归一化
T=mapminmax('reverse',T,ps_output);
function  error = fitness(x,inputnum,hiddennum,outputnum,net,P,T)
%该函数用来计算适应度值
%% 输入
%x                     个体
%inputnum        输入层节点数
%hiddennum     隐含层节点数
%outputnum     输出层节点数
%net                 网络
%P                    训练输入数据
%T                     训练输出数据
%% 输出
%error       个体适应度值
%% 提取
M=size(P,2);
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);
%% 网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%% 训练网络
net=train(net,P,T);
%% 测试
Y=sim(net,P);
error=sum(abs(Y-T).^2)./M;
end

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.