布谷鸟优化算法CS优化BP神经网络(CS-BP)时间序列预测-MATLAB源代码-附带使用教程及注意事项

%% 布谷鸟优化算法CS优化BP神经网络CS-BP时间序列预测 %% 1.初始化 clear close all clc warning off format shortg %% 2.读取数据 data=xlsread('数据.xlsx'); %Matlab2021版本以上无法使用xlsread函数，可用Load函数替代 lag=4; %lag代表延迟，即本网络的输入节点数 for i=1:length(data)-lag deal_data(i,:)=data(i:i+lag)'; %转置 end input=deal_data(:,1:end-1); output=deal_data(:,end); N=length(output); %计算样本数量 testNum=50; %设定测试集样本数量，从数据集后面选取 trainNum=N-testNum; %设定训练集样本数量 %% 3.设置训练集和测试集 input_train = input(1:trainNum,:)'; % 训练集输入 output_train =output(1:trainNum)'; % 训练集输出 input_test =input(trainNum+1:trainNum+testNum,:)'; % 测试集输入 output_test =output(trainNum+1:trainNum+testNum)'; % 测试集输出 %% 4.数据归一化 [inputn,inputps]=mapminmax(input_train,0,1); % 训练集输入归一化到[0,1]之间 [outputn,outputps]=mapminmax(output_train); % 训练集输出归一化到默认区间[-1, 1] inputn_test=mapminmax('apply',input_test,inputps); % 测试集输入采用和训练集输入相同的归一化方式 %% 5.求解最佳隐含层 inputnum=size(input,2); %size用来求取矩阵的行数和列数，1代表行数，2代表列数 outputnum=size(output,2); disp(['输入层节点数：',num2str(inputnum),', 输出层节点数：',num2str(outputnum)]) disp(['隐含层节点数范围为 ',num2str(fix(sqrt(inputnum+outputnum))+1),' 至 ',num2str(fix(sqrt(inputnum+outputnum))+10)]) disp(' ') disp('最佳隐含层节点的确定...') %根据hiddennum=sqrt(m+n)+a，m为输入层节点数，n为输出层节点数，a取值[1,10]之间的整数 MSE=1e+5; %误差初始化 for hiddennum=fix(sqrt(inputnum+outputnum))+1:fix(sqrt(inputnum+outputnum))+10 net=newff(inputn,outputn,hiddennum); %构建BP网络 % 设置网络参数 net.trainParam.epochs=1000; % 设置训练次数 net.trainParam.lr=0.01; % 设置学习速率 net.trainParam.goal=0.000001; % 设置训练目标最小误差 % 进行网络训练 net=train(net,inputn,outputn); an0=sim(net,inputn); %仿真结果 mse0=mse(outputn,an0); %仿真的均方误差 disp(['当隐含层节点数为',num2str(hiddennum),'时，训练集均方误差为：',num2str(mse0)]) %不断更新最佳隐含层节点 if mse0<MSE MSE=mse0; hiddennum_best=hiddennum; end end disp(['最佳隐含层节点数为：',num2str(hiddennum_best),'，均方误差为：',num2str(MSE)]) %% 6.构建最佳隐含层的BP神经网络 disp(' ') disp('BP神经网络：') net0=newff(inputn,outputn,hiddennum_best,{'tansig','purelin'},'trainlm');% 建立模型 %网络参数配置 net0.trainParam.epochs=1000; % 训练次数 net0.trainParam.lr=0.01; % 学习速率 net0.trainParam.goal=0.00001; % 训练目标最小误差 net0.trainParam.show=25; % 显示频率 net0.trainParam.mc=0.01; % 动量因子 net0.trainParam.min_grad=1e-6; % 最小性能梯度 net0.trainParam.max_fail=6; % 最高失败次数 %开始训练 net0=train(net0,inputn,outputn); %预测 an0=sim(net0,inputn_test); %用训练好的模型进行仿真 %预测结果反归一化 test_simu0=mapminmax('reverse',an0,outputps); %把仿真得到的数据还原为原始的数量级 %误差指标 [mae0,mse0,rmse0,mape0,error0,errorPercent0]=calc_error(output_test,test_simu0); %% 7.CS算法权值阈值寻优 disp(' ') disp('CS-BP神经网络：') net=newff(inputn,outputn,hiddennum_best,{'tansig','purelin'},'trainlm');% 建立模型 %网络参数配置 net.trainParam.epochs=1000; % 训练次数 net.trainParam.lr=0.01; % 学习速率 net.trainParam.goal=0.00001; % 训练目标最小误差 net.trainParam.show=25; % 显示频率 net.trainParam.mc=0.01; % 动量因子 net.trainParam.min_grad=1e-6; % 最小性能梯度 net.trainParam.max_fail=6; % 最高失败次数 % CS参数初始化 popsize=10; %初始种群规模 maxgen=50; %最大进化代数 dim=inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+outputnum; %自变量个数 lb=repmat(-3,1,dim); %自变量下限 ub=repmat(3,1,dim); %自变量上限 pa=0.25; %发现一个外来鸟蛋概率 % 种群鸟窝位置初始化 for i=1:popsize nest(i,:)=lb+(ub-lb).*rand(size(lb)); end % 适应度计算 for j=1:size(nest,1) fit(j)=fitness(nest(j,:),inputnum,hiddennum_best,outputnum,net,inputn,outputn,output_train,inputn_test,outputps,output_test); end % 当前最优鸟窝 [fmin,K]=min(fit) ; bestnest=nest(K,:); % 每次迭代最佳适应度初始化 curve = ones(maxgen,1); %% 8.CS开始优化 h0 = waitbar(0,'进度','Name','CS优化中...',... 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); setappdata(h0,'canceling',0); for d=1:maxgen % 通过levy飞行产生新解 % Levy指数和系数 beta=3/2; sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta); for j=1:popsize s=nest(j,:); u=randn(size(s))*sigma; v=randn(size(s)); step=u./abs(v).^(1/beta); stepsize=0.01*step.*(s-bestnest); %位置更新 s=s+stepsize.*randn(size(s)); %位置越界检查 new_nest(j,:)=min(s,ub); %上界 new_nest(j,:)=max(s,lb); %下界 end % 与上一代比较，更新适应度较优的鸟窝 for j=1:size(nest,1) fnew=fitness(new_nest(j,:),inputnum,hiddennum_best,outputnum,net,inputn,outputn,output_train,inputn_test,outputps,output_test); if fnew<=fit(j) fit(j)=fnew; nest(j,:)=new_nest(j,:); end end % 根据发现概率发现鸟蛋，舍弃糟糕的鸟窝 K=rand(size(nest))>pa; stepsize=rand*(nest(randperm(popsize),:)-nest(randperm(popsize),:)); new_nest=nest+stepsize.*K; for j=1:size(new_nest,1) s=new_nest(j,:); %位置越界检查 new_nest(j,:)=min(s,ub); %上界 new_nest(j,:)=max(s,lb); %下界 end % 列出当前最佳的鸟巢 for j=1:size(nest,1) fnew=fitness(new_nest(j,:),inputnum,hiddennum_best,outputnum,net,inputn,outputn,output_train,inputn_test,outputps,output_test); if fnew<=fit(j) fit(j)=fnew; nest(j,:)=new_nest(j,:); end end [fnew,K]=min(fit) ; best=nest(K,:); % 更新当代最优鸟巢与位置 更新榜单 if fnew<fmin fmin=fnew; bestnest=best; end curve(d) = fmin; % 更新迭代得到的最佳适应度 waitbar(d/maxgen,h0,[num2str(d/maxgen*100),'%']) if getappdata(h0,'canceling') break end end delete(h0) setdemorandstream(pi); %% 9.绘制CS进化曲线 figure plot(curve,'b:.','linewidth',1.5) xlabel('进化代数') ylabel('均方误差') legend('最佳适应度') title('CS进化曲线') w1=bestnest(1:inputnum*hiddennum_best); %输入层到中间层的权值 B1=bestnest(inputnum*hiddennum_best+1:inputnum*hiddennum_best+hiddennum_best); %中间各层神经元阈值 w2=bestnest(inputnum*hiddennum_best+hiddennum_best+1:inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum); %中间层到输出层的权值 B2=bestnest(inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+1:inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+outputnum); %输出层各神经元阈值 %矩阵重构 net.iw{1,1}=reshape(w1,hiddennum_best,inputnum); net.lw{2,1}=reshape(w2,outputnum,hiddennum_best); net.b{1}=reshape(B1,hiddennum_best,1); net.b{2}=reshape(B2,outputnum,1); %% 10.优化后的CS-BP神经网络训练 net=train(net,inputn,outputn); %train函数用于训练神经网络，调用蓝色仿真界面 %% 11.CS-BP神经网络测试 an1=sim(net,inputn_test); test_simu1=mapminmax('reverse',an1,outputps); %数据反归一化 %计算误差指标 [mae1,mse1,rmse1,mape1,error1,errorPercent1]=calc_error(output_test,test_simu1); %% 12.作图 figure plot(output_test,'k:o','linewidth',1.5) hold on plot(test_simu0,'rs-','linewidth',1.5,'markerfacecolor','r') hold on plot(test_simu1,'bh-','linewidth',1.5,'markerfacecolor','k') legend('实际值','BP预测值','CS-BP预测值') xlabel('测试集样本'),ylabel('指标值') title('预测值和实际值对比') figure plot(err