粒子群算法优化bp神经网络

%% 该代码为基于PSO和BP网络的预测 %% 清空环境 close all clear all clc %读取数据 %load data input_train input_test output_train output_test %节点个数 %inputnum=4; % %hiddennum=9; %2*inputnum+1 %outputnum=1; % %训练数据和预测数据 %input_train; %input_test; %output_train; %output_test; %样本输入输出数据归一化 %[inputn,inputps]=mapminmax(input_train); %[outputn,outputps]=mapminmax(output_train); %inputn1=inputn'; %outputn1=outputn'; %构建网络 %net=newff(inputn1,outputn1,hiddennum); AA=xlsread('总论文数据.xlsx'); %load('.txt') %t=randperm(size(input,1)); [hang,lie]=size(AA); lie=lie-2; input=ones(hang,lie); xunlianji=1120; ceshiji=225; input(:,1)=AA(:,2); input(:,2)=AA(:,3); input(:,3)=AA(:,4); input1=input'; [in,ps]=mapminmax(input1,0,1); %归一化 （0，1） in1=in'; output=AA(:,5); %In=input(t(1:1974),:); T=output(1:xunlianji,:)'; T_test=output(xunlianji+1:end,:)'; %P1=In(1:1600,1:15)'; %P1_test=In(1601:1974,1:15)'; P=in1(1:xunlianji,:)'; P_test=in1(xunlianji+1:end,:)'; %[input1,PS1]=mapminmax(input(:,1),0,1); %[input2,PS2]=mapminmax(input(:,2),0,1); %[input3,PS3]=mapminmax(input(:,3),0,1); %[input4,PS4]=mapminmax(input(:,4),0,1); %P=input1(1:300,:)'; %P_test=input1(301:360,:)'; inputnum=size(P,1); % 输入层神经元个数 outputnum=size(T,1); % 输出层神经元个数 hiddennum=2*inputnum+1; % 输入向量的最大值和最小值 threshold=[0 1;0 1;0 1;0 1]; inputnum=size(P,1); % 输入层神经元个数 outputnum=size(T,1); % 输出层神经元个数 w1num=inputnum*hiddennum; % 输入层到隐层的权值个数 w2num=outputnum*hiddennum;% 隐层到输出层的权值个数 %% 新建BP网络 net=newff(P,T,[hiddennum],{'logsig','purelin'},'trainlm'); net.trainParam.epochs=1000; % 训练次数，这里设置为1000次 net.trainParam.lr=0.01; % 学习率，这里设置为0.01 net.trainParam.goal=0.001; % 训练目标最小误差，这里设置为0.0001 net.trainParam.show=25; % 显示频率，这里设置为每训练25次显示一次 net.trainParam.mc=0.01; % 动量因子 net.trainParam.min_grad=1e-6; % 最小性能梯度 net.trainParam.max_fail=6; %net=newff(inputn1,outputn1,[hiddennum],{'tansig','logsig'},'trainlm'); % 参数初始化 %粒子群算法中的两个参数 c1 = 1.49445; c2 = 1.49445; %c1=0.5; %c2=1.1; maxgen=50; % 进化次数 sizepop=20; %种群规模 wmax=0.9; wmin=0.4; Vmax=1; Vmin=-1; popmax=5; popmin=-5; %% 产生初始粒子和速度 for i=1:sizepop %随机产生一个种群 pop(i,:)=5*rands(1,inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum); %初始种群 vov(i,:)=rands(1,inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum); %初始化速度 %计算适应度 fitness(i)=fun(pop(i,:),inputnum,hiddennum,outputnum,net,P,T); %染色体的适应度 end % 个体极值和群体极值 [bestfitness bestindex]=min(fitness); zbest=pop(bestindex,:); %全局最佳 gbest=pop; %个体最佳 fitnessgbest=fitness; %个体最佳适应度值 fitnesszbest=bestfitness; %全局最佳适应度值 %% 迭代寻优 for i=1:maxgen %粒子位置和速度更新 for j=1:sizepop w=wmax-(wmax-wmin)*j/maxgen; %速度更新 vov(j,:) = w*vov(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:)); vov(j,find(vov(j,:)>Vmax))=Vmax; vov(j,find(vov(j,:)<Vmin))=Vmin; %种群更新 pop(j,:)=pop(j,:)+0.5*vov(j,:); pop(j,find(pop(j,:)>popmax))=popmax; pop(j,find(pop(j,:)<popmin))=popmin; %引入变异算子，重新初始化粒子 if rand>0.9 k=ceil(21*rand); pop(j,k)=rand; end %新粒子适应度值 fitness(j)=fun(pop(j,:),inputnum,hiddennum,outputnum,net,P,T); end %%个体极值和群体极值更新 for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:); fitnesszbest = fitness(j); end end %%每代最优值记录到yy数组中 yy(i)=fitnesszbest; end %% 结果分析 plot(yy,'r-','linewidth',2) title(['适应度曲线 ' '终止代数＝' num2str(maxgen)],'fontsize',12); xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12); x=zbest; %% 把最优初始阀值权值赋予网络预测 % %用遗传算法优化的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.01; net.trainParam.goal=0.00001; %网络训练 [net,tr]=train(net,P,T); %% BP网络预测 %数据归一化 %inputn_test=mapminmax('apply',input_test,inputps); an=sim(net,P_test); %fangzhen bn=sim(net,P); %test_simu=mapminmax('reverse',an,outputps); error=an-T_test; error1=norm(error); error2=norm(bn-T); %mm=sum(error) figure(2) plot(error) title('仿真预测误差','fontsize',12); xlabel('仿真次数','fontsize',12);ylabel('误差百分值','fontsize',12); %% 添加传统的BP net0=newff(P,T,hiddennum,{'tansig','purelin'},'trainlm');% 建立模型 %网络参数配置 net0.trainParam.epochs=1000; % 训练次数，这里设置为1000次 net0.trainParam.lr=0.01; % 学习速率，这里设置为0.01 net0.trainParam.goal=0.001; % 训练目标最小误差，这里设置为0.0001 net0.trainParam.show=25; % 显示频率，这里设置为每训练25次显示一次 net0.trainParam.mc=0.01; % 动量因子 net0.trainParam.min_grad=1e-6; % 最小性能梯度 net0.trainParam.max_fail=6; % 最高失败次数 %开始训练 net0=train(net0,P,T); %预测 an0=sim(net0,P_test); %用训练好的模型进行仿真 bnn0=sim(net0,P); errorrr=an0-T_test; %预测结果反归一化与误差计算 %test_simu0=mapminmax('reverse',an0,outputps); %把仿真得到的数据还原为原始的数量级 %误差指标 %[mae0,mse0,rmse0,mape0,error0,errorPercent0 R20]=compute_error(T_test,an0); figure(4) plot(T_test','r-v','linewidth',1,'markerfacecolor','r') hold on plot(an','m-o','linewidth',1,'markerfacecolor','m') hold on plot(an0','b-o','linewidth',1,'markerfacecolor','m') legend('真实值','PSOBP预测值','BP预测值') xlabel('测试样本编号') ylabel('指标值') title('PSOBP神经网络测试样本预测值和真实值对比图') figure(5) plot(error','r-v','linewidth',1,'markerfacecolor','r') hold on plot(errorrr','b-o','linewidth',1,'markerfacecolor','m') legend('PSOBP误差值','BP预测误差') xlabel('测试样本编号') ylabel('指标值') title('PSOBP神经网络的测试样本的误差') figure(6) plot(T_test,'b-*','linewidth',1) hold on plot(an,'r-o','linewidth',1,'markerfacecolor','k') legend('真实值','PSO-BP预测值') xlabel('测试样本编号') ylabel('指标值') title('PSO优化前后的BP神经网络预测值和真实值对比图') figure plot(T,'b-*','linewidth',1) hold on plot(bn,'r-o','linewidth',1,'markerfacecolor','k') legend('真实值','PSO-BP预测值') xlabel('测试样本编号') ylabel('指标值') title('PSO优化前后的BP神经网络预测值和真实值对比图') figure plot(error,'rv-','markerfacecolor','r') xlabel('测试样本编号') ylabel('预测偏差') title('PSOBP预测误差样本测试集')