灰狼优化算法GWO优化BP神经网络(GWO-BP)时间序列预测-时序预测模型-MATLAB源代码-附带使用教程及注意事项

%% 1、初始化 clc clear close all format bank %2位小数，format short精确4位，format long精确15位 %% 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 hiddennum0=fix(sqrt(inputnum+outputnum))+1:fix(sqrt(inputnum+outputnum))+10 net=newff(inputn,outputn,hiddennum0); %构建BP网络 % 设置网络参数 net.trainParam.epochs=1000; % 训练次数 net.trainParam.lr=0.1; % 学习速率 net.trainParam.goal=0.000001; % 训练目标最小误差 % 进行网络训练 net=train(net,inputn,outputn); an0=sim(net,inputn); %仿真结果 mse0=mse(outputn,an0); %仿真的均方误差 disp(['当隐含层节点数为',num2str(hiddennum0),'时，训练集均方误差为：',num2str(mse0)]) % 更新最佳的隐含层节点 if mse0<MSE MSE=mse0; hiddennum=hiddennum0; end end disp(['最佳隐含层节点数为：',num2str(hiddennum),'，均方误差为：',num2str(MSE)]) %% 6、最佳隐含层BP网络 net=newff(inputn,outputn,hiddennum,{'tansig','purelin'},'trainlm'); % 设置网络参数 net.trainParam.epochs=1000; net.trainParam.lr=0.1; net.trainParam.goal=0.00001; % 网络训练 [net,~]=train(net,inputn,outputn); % 网络测试 an=sim(net,inputn_test); test_simu1=mapminmax('reverse',an,outputps); error1=test_simu1-output_test; %% 7、灰狼算法优化 SearchAgents_no=20; % 狼群数量 Max_iteration=40; % 最大迭代次数 dim=121; % 此例需要优化两个参数c和g lb=-10*ones(1,121); % 参数取值下界 ub=10*ones(1,121); % 参数取值上界 % 节点总数 numsum=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum; lenchrom=ones(1,numsum); bound=[-3*ones(numsum,1) 3*ones(numsum,1)]; %数据范围 Alpha_pos=zeros(1,dim); % 初始化Alpha狼的位置 Alpha_score=inf; % 初始化Alpha狼的目标函数值， Beta_pos=zeros(1,dim); % 初始化Beta狼的位置 Beta_score=inf; % 初始化Beta狼的目标函数值， Delta_pos=zeros(1,dim); % 初始化Delta狼的位置 Delta_score=inf; % 初始化Delta狼的目标函数值， Positions=initialization(SearchAgents_no,dim,ub,lb); Convergence_curve=zeros(1,Max_iteration); l=0; % 循环计数器 h0=waitbar(0,'GWO优化BP即将完成，请等待...'); while l<Max_iteration % 对迭代次数循环 for i=1:size(Positions,1) % 遍历每个狼 % 若搜索位置超过了搜索空间，需要重新回到搜索空间 Flag4ub=Positions(i,:)>ub; Flag4lb=Positions(i,:)<lb; % 若狼的位置在最大值和最小值之间，则位置不需要调整，若超出最大值，最回到最大值边界； % 若超出最小值，最回答最小值边界 Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb; % ~表示取反 x= Positions(i,:); % 计算适应度函数值 cmd = ['-w1 ',num2str(Positions(i,1)),' -B1 ',num2str(Positions(i,2)),' -w2 ',num2str(Positions(i,3)),' -B2 ',num2str(Positions(i,2))]; fitness=fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn); if fitness<Alpha_score % 如果目标函数值小于Alpha狼的目标函数值 Alpha_score=fitness; % 则将Alpha狼的目标函数值更新为最优目标函数值 Alpha_pos=Positions(i,:); % 同时将Alpha狼的位置更新为最优位置 end if fitness>Alpha_score && fitness<Beta_score % 如果目标函数值介于于Alpha狼和Beta狼的目标函数值之间 Beta_score=fitness; % 则将Beta狼的目标函数值更新为最优目标函数值 Beta_pos=Positions(i,:); % 同时更新Beta狼的位置 end if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score % 如果目标函数值介于于Beta狼和Delta狼的目标函数值之间 Delta_score=fitness; % 则将Delta狼的目标函数值更新为最优目标函数值 Delta_pos=Positions(i,:); % 同时更新Delta狼的位置 end end a=2-l*((2)/Max_iteration); % 对每一次迭代，计算相应的a值，a decreases linearly fron 2 to 0 for i=1:size(Positions,1) % 遍历每个狼 for j=1:size(Positions,2) % 遍历每个维度 % 包围猎物，位置更新 r1=rand(); % r1 is a random number in [0,1] r2=rand(); % r2 is a random number in [0,1] A1=2*a*r1-a; % 计算系数A，Equation (3.3) C1=2*r2; % 计算系数C，Equation (3.4) % Alpha狼位置更新 D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1 X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1 r1=rand(); r2=rand(); A2=2*a*r1-a; % 计算系数A，Equation (3.3) C2=2*r2; % 计算系数C，Equation (3.4) % Beta狼位置更新 D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); % Equation (3.5)-part 2 X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2 r1=rand(); r2=rand(); A3=2*a*r1-a; % 计算系数A，Equation (3.3) C3=2*r2; % 计算系数C，Equation (3.4) % Delta狼位置更新 D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); % Equation (3.5)-part 3 X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3 % 位置更新 Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7) end end l=l+1; Convergence_curve(l)=Alpha_score; waitbar(Max_iteration,h0) end close(h0) %% 8、参数选择结果赋值 x=Beta_pos; 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; %% 9、GWO-BP网络 % 网络参数 net.trainParam.epochs=1000; net.trainParam.lr=0.1; net.trainParam.goal=0.000001; % 网络训练 [net,per2]=train(net,inputn,outputn); % 网络预测 an=sim(net,inputn_test); test_simu=mapminmax('reverse',an,outputps); error=test_simu-output_test; %% 10.输出结果 %计算BP各项误差参数 [~,len]=size(output_test); % len获取测试样本个数，数值等于testNum，用于求各指标平均值 SSE1=sum(error1.^2); % 误差平方和 MAE1=sum(abs(error1))/len; % 平均绝对误差 MSE1=error1*error1'/len; % 均方误差 RMSE1=MSE1^(1/2); % 均方根误差 MAPE1=mean(abs(error1./output_test)); % 平均百分比误差 %计算GWO-BP各项误差参数 SSE=sum(error.^2); % 误差平方和 MAE=sum(abs(error))/len; % 平均绝对误差 MSE=error*error'/len; % 均方误差 RMSE=