建建建设网站首页嘉上营销
%% 基于模拟退火遗传算法优化BP神经网络的钢带厚度预测
 clear
 clc
 close all
 format short
 %% 加载训练数据
 Xtr=xlsread('train_data.xlsx');
 DD=size(Xtr,2);
 input_train=Xtr(:,1:DD-1)';% 
 output_train=Xtr(:,DD)';% 
 %% 加载测试数据
 Xte=xlsread('test_data.xlsx');
 input_test=Xte(:,1:DD-1)';% 
 output_test=Xte(:,DD)';% 
 %% BP网络设置
 %节点个数
 [inputnum,N]=size(input_train);%输入节点数量
 outputnum=size(output_train,1);%输出节点数量
 hiddennum=20; %隐含层节点
 %% 数据归一化
 [inputn,inputps]=mapminmax(input_train,0,1);
 [outputn,outputps]=mapminmax(output_train,0,1);% 归一化到【0 1】之间
 %% GA算法参数初始化
 D=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;
 net=newff(inputn,outputn,hiddennum);
 %参数定义
 maxgen=20;                         %进化代数,即迭代次数
 sizepop=10;                        %种群规模
 fselect='roulette';                 %染色体的选择方法,您可以选择:锦标赛法- 'tournament';轮盘赌法-'roulette'
 fcode='float';                       %编码方法,您可以选择:浮点法-'float';grey法则--'grey';二进制法-'binary' 
 pcross=0.7;                       %交叉概率选择,0和1之间
 fcross='float';                   %交叉方法选择,您可以选择: 浮点交叉-'float';单点交叉-'simple';均匀交叉-'uniform'
 pmutation=0.1;                    %变异概率选择,0和1之间
 fmutation='float';                 %变异方法选择,您可以选择:浮点法-'float';单点法-'simple';
 lenchrom=ones(1,D);%[1 1];          %每个变量的字串长度,如果是浮点变量,则长度都为1
 % bound=[-2.048,2.048;-2.048 2.048];
 %% 优化参数的取值范围
 ub=3*ones(1,D); % 参数取值上界
 lb=-3*ones(1,D); % 参数取值下界
 bound=[lb' ub'];
  
 for i=1:sizepop
     %随机产生一个种群
     individuals.chrom(i,:)=Code(lenchrom,fcode,bound);    %编码(binary和grey的编码结果为一个实数,float的编码结果为一个实数向量)
     x=Decode(lenchrom,bound,individuals.chrom(i,:),fcode);%解码(binary和grey的解码结果为一个二进制串,float的解码结果为一个实数向量)
     %计算适应度
     individuals.fitness(i)=objfun_BP(x,inputnum,hiddennum,outputnum,net,inputn,outputn);%计算适应度
 end
  
 %找最好的染色体
 [bestfitness bestindex]=min(individuals.fitness);
 bestchrom=individuals.chrom(bestindex,:);  %最好的染色体
 avgfitness=sum(individuals.fitness)/sizepop; %染色体的平均适应度
 trace=[avgfitness bestfitness];
 % 进化开始
 for i=1:maxgen
     disp(['迭代次数:' num2str(i)])
     % 选择
     individuals=Select(individuals,sizepop,fselect); 
     avgfitness=sum(individuals.fitness)/sizepop;
     %交叉
     individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,fcross,[i maxgen],fcode,bound,individuals.fitness,bestfitness,avgfitness);
     % 变异
     individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,fmutation,[i maxgen],fcode,bound,individuals.fitness,bestfitness,avgfitness);
     % 计算适应度 
     for j=1:sizepop
         x=Decode(lenchrom,bound,individuals.chrom(j,:),fcode); %解码
        individuals.fitness(j)=objfun_BP(x,inputnum,hiddennum,outputnum,net,inputn,outputn);%计算适应度  
     end
   %找到最小和最大适应度的染色体及它们在种群中的位置
     [newbestfitness,newbestindex]=min(individuals.fitness);
     [worestfitness,worestindex]=max(individuals.fitness);
     % 代替上一次进化中最好的染色体
     if bestfitness>newbestfitness
         bestfitness=newbestfitness;
         bestchrom=individuals.chrom(newbestindex,:);
     end
     individuals.chrom(worestindex,:)=bestchrom;
     individuals.fitness(worestindex)=bestfitness;
     avgfitness=sum(individuals.fitness)/sizepop;
     trace=[trace;avgfitness bestfitness]; %记录每一代进化中最好的适应度和平均适应度
 end
  
 %% 画出适应度曲线
 figure
 plot(trace(:,1),'b--');
 hold on
 plot(trace(:,2),'r-.');
 title(['适应度曲线  ' '终止代数=' num2str(maxgen)]);
 xlabel('进化代数');ylabel('适应度');
 legend('平均适应度','最佳适应度')
 axis tight
 disp('适应度                   变量');
 x=Decode(lenchrom,bound,bestchrom,fcode);
 % 窗口显示
 disp([bestfitness x]);
 % 把最优初始阀值权值赋予网络预测
 %用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--','linewidth',1)
 hold on
 plot(1:length(ABC_sim),ABC_sim,'r-.','linewidth',1)
 xlabel('训练样本','FontSize',12);
 ylabel('钢带厚度(mm)')
 legend('实际值','预测值');
 axis tight
 title('GASA-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--','linewidth',1)
 hold on
 plot(1:length(ABCBPsim),ABCBPsim,'r-.','linewidth',1)
 xlabel('测试样本','FontSize',12);
 ylabel('钢带厚度(mm)')
 legend('实际值','预测值');
 axis tight
 title('GASA-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--','linewidth',1)
 hold on
 plot(1:length(T_sim1),T_sim1,'r-.','linewidth',1)
 xlabel('训练样本','FontSize',12);
 ylabel('钢带厚度(mm)')
 legend('实际值','预测值');
 axis tight
 title('BP神经网络')
 figure
 plot(1:length(output_test), output_test,'b--','linewidth',1)
 hold on
 plot(1:length(Tsim1),Tsim1,'r-.','linewidth',1)
 xlabel('测试样本','FontSize',12);
 ylabel('钢带厚度(mm)')
 legend('实际值','预测值');
 axis tight
 title('BP神经网络')
 %% GASA-BP和BP对比图
 figure
 plot(1:length(output_train),output_train,'b--','linewidth',1)
 hold on
 plot(1:length(ABC_sim),ABC_sim,'r-.','linewidth',1)
 hold on
 plot(1:length(T_sim1),T_sim1,'g.-','linewidth',1)
 xlabel('训练样本','FontSize',12);
 ylabel('钢带厚度(mm)')
 legend('实际值','GASA-BP预测值','BP预测值');
 axis tight
 figure
 plot(1:length(output_test), output_test,'b--','linewidth',1)
 hold on
 plot(1:length(ABCBPsim),ABCBPsim,'r-.','linewidth',1)
 hold on
 plot(1:length(Tsim1),Tsim1,'g.-','linewidth',1)
 xlabel('测试样本','FontSize',12);
 ylabel('钢带厚度(mm)')
 legend('实际值','GASA-BP预测值','BP预测值');
 axis tight
 %% 结果评价
 Result1=CalcPerf(output_test,ABCBPsim);
 MSE1=Result1.MSE;
 RMSE1=Result1.RMSE;
 MAPE1=Result1.Mape;
 disp(['GASABP-RMSE = ', num2str(RMSE1)])
 disp(['GASABP-MSE  = ', num2str(MSE1)])
 disp(['GASABP-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)])
