%% 基于天牛须搜索算法优化BP神经网络预测研究
%% 清空环境变量
clear all
close all
clc
tic
%% 加载数据
X=xlsread('数据.xlsx');
for i=1:111
Y(i)=X(10*(i-1)+1,2);
end
%% 时间序列数据预处理
xx = [];
yy = [];
num_input = 3;
tr_len =round(0.8*length(Y));%训练集样本个数
for i = 1:length(Y)-num_input
xx = [xx; Y(i:i+num_input-1)];
yy = [yy; Y(i+num_input)];
end
%BP的训练集的输入和输出
input_train = xx(1:tr_len, :)';
output_train = yy(1:tr_len)';
%BP的测试集的输入和输出
input_test = xx(tr_len+1:end, :)';
output_test = yy(tr_len+1:end)';
%% BP网络设置
%节点个数
[inputnum,N]=size(input_train);%输入节点数量
outputnum=size(output_train,1);%输出节点数量
hiddennum=10; %隐含层节点
%% 数据归一化
[inputn,inputps]=mapminmax(input_train,0,1);
[outputn,outputps]=mapminmax(output_train,0,1);% 归一化到[0,1]之间
%% 创建网络
net=newff(inputn,outputn,hiddennum);
net.trainParam.epochs = 1000;
net.trainParam.goal = 1e-3;
net.trainParam.lr = 0.01;
%% 天牛须算法初始化
eta=0.8;
c=5;%步长与初始距离之间的关系
step=10;%初始步长
n=100;%迭代次数
k=inputnum*hiddennum+outputnum*hiddennum+hiddennum+outputnum;
x=rands(k,1);
bestX=x;
bestY=fitness(bestX,inputnum,hiddennum,outputnum,net,inputn,outputn);
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,inputn,outputn);
xright=x-dir*d0/2;
fright=fitness(xright,inputnum,hiddennum,outputnum,net,inputn,outputn);
x=x-step*dir*sign(fleft-fright);
y=fitness(x,inputnum,hiddennum,outputnum,net,inputn,outputn);
if y<bestY
bestX=x;
bestY=y;
end
if y<0.001
bestX=x;
bestY=y;
end
x_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
%% 可视化
figure(1)
%plot(x_store(1,:),x_store(end,:),'r-o')
hold on,
plot(x_store(1,:),fbest_store,'b-.')
xlabel('Iteration')
ylabel('BestFit')
toc
%%
x=bestX';
% 把最优初始阀值权值赋予网络预测
%用ABC优化的BP网络进行值预测
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}=B2;
% BP网络训练
% 网络进化参数
net.trainParam.epochs=100;
net.trainParam.lr=0.1;
net.trainParam.mc = 0.8;%动量系数，[0 1]之间
net.trainParam.goal=0.001;
% 网络训练
net=train(net,inputn,outputn);
% BP训练集预测
BP_sim=sim(net,inputn);
% 网络输出反归一化
ABC_sim=mapminmax('reverse',BP_sim,outputps);
% 绘图
figure
plot(1:length(output_train),output_train,'b-o','linewidth',1)
hold on
plot(1:length(ABC_sim),ABC_sim,'r-square','linewidth',1)
xlabel('训练样本','FontSize',12);
legend('实际值','预测值');
axis tight
title('BAS-BP神经网络')
% BP测试集
% 测试数据归一化
inputn_test=mapminmax('apply',input_test,inputps);
%预测输出
an=sim(net,inputn_test);
ABCBPsim=mapminmax('reverse',an,outputps);
figure
plot(1:length(output_test), output_test,'b-o','linewidth',1)
hold on
plot(1:length(ABCBPsim),ABCBPsim,'r-square','linewidth',1)
xlabel('测试样本','FontSize',12);
legend('实际值','预测值');
axis tight
title('BAS-BP神经网络')
%% 未优化的BP神经网络
net1=newff(inputn,outputn,hiddennum);
% 网络进化参数
net1.trainParam.epochs=100;
net1.trainParam.lr=0.1;
net1.trainParam.mc = 0.8;%动量系数，[0 1]之间
net1.trainParam.goal=0.001;
% 网络训练
net1=train(net1,inputn,outputn);
% BP训练集预测
BP_sim1=sim(net1,inputn);
% 网络输出反归一化
T_sim1=mapminmax('reverse',BP_sim1,outputps);
%预测输出
an1=sim(net1,inputn_test);
Tsim1=mapminmax('reverse',an1,outputps);
% 绘图
figure
plot(1:length(output_train),output_train,'b-o','linewidth',1)
hold on
plot(1:length(T_sim1),T_sim1,'r-square','linewidth',1)
xlabel('训练样本','FontSize',12);
legend('实际值','预测值');
axis tight
title('BP神经网络')
figure
plot(1:length(output_test), output_test,'b-o','linewidth',1)
hold on
plot(1:length(Tsim1),Tsim1,'r-square','linewidth',1)
xlabel('测试样本','FontSize',12);
legend('实际值','预测值');
axis tight
title('BP神经网络')
%% ABC-BP和BP对比图
figure
plot(1:length(output_train),output_train,'b-o','linewidth',1)
hold on
plot(1:length(ABC_sim),ABC_sim,'r-square','linewidth',1)
hold on
plot(1:length(T_sim1),T_sim1,'g-diamond','linewidth',1)
xlabel('训练样本','FontSize',12);
legend('实际值','BAS-BP预测值','BP预测值');
axis tight
figure
plot(1:length(output_test), output_test,'b-o','linewidth',1)
hold on
plot(1:length(ABCBPsim),ABCBPsim,'r-square','linewidth',1)
hold on
plot(1:length(Tsim1),Tsim1,'g-diamond','linewidth',1)
xlabel('测试样本','FontSize',12);
legend('实际值','BAS-BP预测值','BP预测值');
axis tight
%% 结果评价
Result1=CalcPerf(output_test,ABCBPsim);
MSE1=Result1.MSE;
RMSE1=Result1.RMSE;
MAPE1=Result1.Mape;
disp(['BASBP-RMSE = ', num2str(RMSE1)])
disp(['BASBP-MSE = ', num2str(MSE1)])
disp(['BASBP-MAPE = ', num2str(MAPE1)])
Result2=CalcPerf(output_test,Tsim1);
MSE2=Result2.MSE;
RMSE2=Result2.RMSE;
MAPE2=Result2.Mape;
disp(['BP-RMSE = ', num2str(RMSE2)])
disp(['BP-MSE = ', num2str(MSE2)])
disp(['BP-MAPE = ', num2str(MAPE2)])
