遗传算法优化发酵配方(采用BP神经网络拟合)

%% 遗传算法优化BP神经网络 拟合发酵数据 clc;clear all;close all; %% ----------构建神经网络------------ %读取数据 FermentationData=xlsread('FermentationData.xls','B2:G33');% input_train=FermentationData(:,1:end-1);%载入输入数据 output_train=FermentationData(:,end);%载入输出数据 %% ------BP神经网络相关参数的设定开始----- [k11,k12]=size(input_train);%计算维数 inputnumber=k12;%输入维数 hiddenumber=5;%隐含神经元数 outputnumber=1;%输出维数 %% ------BP神经网络相关参数的设定结束----- %% ---------训练数据和预测数据开始--------- testnumber=5;%用于测试的数据 S1=randperm(k11);%随机序列 v1=S1(testnumber+1:end); v2=S1(1:testnumber); Input_Train=input_train(v1,:)'; Input_Test=input_train(v2,:)'; Output_Train=output_train(v1,:)'; Output_Test=output_train(v2,:)'; %数据处理：输入输出数据归一化处理 [inputdata,inputps]=mapminmax(Input_Train);%输入数据归一化 [outputdata,outputps]=mapminmax(Output_Train);%输出数据归一化 %% ---------训练数据和预测数据结束--------- %新建bp神经网络 net=newff(inputdata,outputdata,hiddenumber); %% 遗传算法参数初始化 Itera_number=5; %进化代数，或遗传算法的迭代次数 pop_scale=20; %种群规模 probability_cross=0.3; %交叉概率选择，0和1之间 pmutation=0.1; %变异概率选择，0和1之间 %总节点数计算 number_sum=inputnumber*hiddenumber+hiddenumber+hiddenumber*outputnumber+outputnumber; Chrom=ones(1,number_sum);%染色体编码 bound=[-3*ones(number_sum,1) 3*ones(number_sum,1)];%数据上下界 %% -------------初始化种群开始--------------- Pop_indiv=struct('fitness',zeros(1,pop_scale), 'chrom',[]); %将种群信息以结构体形式保存 avgfitness=[]; %每一代种群的平均适应度 bestfitness=[]; %每一代种群的最佳适应度 bestchrom=[]; %适应度最好的染色体 for i=1:pop_scale %随机产生一个种群 Pop_indiv.chrom(i,:)=Code(Chrom,bound); %编码（binary和grey的编码结果为一个实数，float的编码结果为一个实数向量） x=Pop_indiv.chrom(i,:); %计算适应度 Pop_indiv.fitness(i)=fun(x,inputnumber,hiddenumber,outputnumber,net,inputdata,outputdata); %染色体的适应度 end %% -------------初始化种群结束--------------- %寻找最佳染色体 [bestfitness bestindex]=min(Pop_indiv.fitness);%求最佳值 bestchrom=Pop_indiv.chrom(bestindex,:);%获取最好的染色体 avgfitness=sum(Pop_indiv.fitness)/pop_scale;%计算染色体的平均适应度 % 记录每一代进化中最好的适应度和平均适应度 trace=[avgfitness bestfitness]; %% ----------迭代求解BP神经网络的最佳初始阀值和权值开始---------- %进度条 wait_hand = waitbar(0,'正在进行遗传算法迭代，请等待……', 'tag', 'TMWWaitbar'); % 进化开始 for i=1:Itera_number waitbar(i/Itera_number,wait_hand);%每循环一次更新一次进度条 %选择操作 Pop_indiv=Select(Pop_indiv,pop_scale); avgfitness=sum(Pop_indiv.fitness)/pop_scale; %交叉操作 Pop_indiv.chrom=Cross(probability_cross,Chrom,Pop_indiv.chrom,pop_scale,bound); %变异操作 Pop_indiv.chrom=Mutation(pmutation,Chrom,Pop_indiv.chrom,pop_scale,i,Itera_number,bound); %计算适应度 for j=1:pop_scale x=Pop_indiv.chrom(j,:); %解码染色体 Pop_indiv.fitness(j)=fun(x,inputnumber,hiddenumber,outputnumber,net,inputdata,outputdata); end %找到最小和最大适应度的染色体及它们在种群中的位置 [newbestfitness,newbestindex]=min(Pop_indiv.fitness); [worestfitness,worestindex]=max(Pop_indiv.fitness); % 代替上一次进化中最好的染色体 if bestfitness>newbestfitness bestfitness=newbestfitness; bestchrom=Pop_indiv.chrom(newbestindex,:); end Pop_indiv.chrom(worestindex,:)=bestchrom; Pop_indiv.fitness(worestindex)=bestfitness; avgfitness=sum(Pop_indiv.fitness)/pop_scale; trace=[trace;avgfitness bestfitness]; %保存每一代种群中的最佳适应度与平均适应度 end delete(wait_hand);%执行完后删除该进度条 %% ----------迭代求解BP神经网络的最佳初始阀值和权值开始---------- %% 遗传算法结果绘图 figure; [r c]=size(trace); plot([1:r]',trace(:,2),'b--'); title(['适应度曲线 ' '终止代数＝' num2str(Itera_number)]); xlabel('进化代数');ylabel('适应度'); legend('平均适应度','最佳适应度'); disp('适应度 变量'); x=bestchrom %% 把最优初始阀值权值赋予网络预测 % %用遗传算法优化的BP网络进行值预测 weigth1=x(1:inputnumber*hiddenumber); B01=x(inputnumber*hiddenumber+1:inputnumber*hiddenumber+hiddenumber); weigth2=x(inputnumber*hiddenumber+hiddenumber+1:inputnumber*hiddenumber+hiddenumber+hiddenumber*outputnumber); B02=x(inputnumber*hiddenumber+hiddenumber+hiddenumber*outputnumber+1:inputnumber*hiddenumber+hiddenumber+hiddenumber*outputnumber+outputnumber); net.iw{1,1}=reshape(weigth1,hiddenumber,inputnumber); net.lw{2,1}=reshape(weigth2,outputnumber,hiddenumber); net.b{1}=reshape(B01,hiddenumber,1); net.b{2}=B02; %% ------用GA得到的阀值和权值为初始数据训练BP网络开始-------- %网络进化参数 net.trainParam.epochs=1000;%训练次数 net.trainParam.lr=0.1; net.trainParam.goal=0.00001; %网络训练 [net,per2]=train(net,inputdata,outputdata); %% ------用GA得到的阀值和权值为初始数据训练BP网络结束-------- %% ---------BP网络预测和绘图------------ % ------------数据归一化开始-------------- inputn_test=mapminmax('apply',Input_Test,inputps); an=sim(net,inputn_test) test_simu=mapminmax('reverse',an,outputps); % ------------数据归一化结束-------------- error1=abs(test_simu-Output_Test)./Output_Test*100;%计算误差 figure; plot(1:length(Output_Test),Output_Test,'-bo',1:length(test_simu),test_simu,'-r*'); legend('实际值','预测值'); xlabel('样本'); ylabel('输出值'); title('实际值和BP预测值的对比'); figure; plot(1:length(error1),error1,'-r*'); xlabel('样本'); ylabel('误差(%)'); title('误差'); %% 求解最优配方--采用fmincon方法求算 global net1 inputps1 outputps1;%定义全局变量 net1=net; inputps1=inputps; outputps1=outputps; x0=[70 10 12 3 0.3]; A1=[]; b1=[]; Aeq=[]; beq=[]; lb1=zeros(1,k12); ub1=zeros(1,k12); for i=1:k12 lb1(1,i)=min(input_train(:,i)); ub1(1,i)=max(input_train(:,i)); end [RecipeGood,Obj]= fmincon(@Recipefun,x0,A1,b1,Aeq,beq,lb1,ub1); disp('最佳配方'); RecipeGood disp('最佳产率'); -Obj