%% 差分进化算法应用于优化BP神经网络的初始权值和阈值
%% 清空环境变量
clear all;
clc;
warning off
load v357;
load y357;
Pn_train=v;
Tn_train=y;
Pn_test=v;
Tn_test=y;
P_train=v;
T_train=y;
% P_train=[0 25.27 44 62.72 81.4 100.2;
% 290.5 268.8 247.2 224.5 206 184.4;
% 0 16.12 33.25 50.42 67.62 84.73;
% 542.5 517.8 493 465.3 435.6 410.8;
% 0 11.1 28.1 44.93 61.38 78.57;
% 826.1 800.2 769.1 740.0 706.2 669.3];
% T_train=[0 1 2 3 4 5];%以上是未处理的数据
% P_test=[0 25.25 43 62.75 81.6 100.7;
% 290.3 268.4 247.5 224.6 206 184.2;
% 0 16.14 33.26 50.47 67.68 84.79;
% 542.7 517.9 495 465.8 435.6 410.9;
% 0 11.4 28.6 44.94 61.36 78.59;
% 826.3 800.7 769.8 740.5 706.7 669.3];
% T_test=[0 1 2 3 4 5];
% Pn_train=[0 0.252 0.439 0.626 0.813 1 0 0.19 0.392 0.595 0.798 1 0 0.141 0.358 0.572 0.781 1;
% 1 0.795 0.592 0.378 0.204 0 1 0.815 0.626 0.415 0.189 0 1 0.835 0.637 0.451 0.235 0];
% %T 为目标矢量 ,归一化后的数据
% Tn_train=[0.05,0.23,0.41,0.59,0.77,0.95,0.05,0.23,0.41,0.59,0.77,0.95,0.05,0.23,0.41,0.59,0.77,0.95];
% Pn_test=[ 0 0.17 0.39 0.595 0.798 1 0 0.141 0.358 0.572 0.781 1 0 0.258 0.439 0.626 0.813 1;
% 1 0.815 0.625 0.415 0.189 0 1 0.835 0.635 0.451 0.235 0 1 0.795 0.599 0.378 0.204 0 ];
% Tn_test=[0.05,0.23,0.41,0.59,0.77,0.95,0.05,0.23,0.41,0.59,0.77,0.95,0.05,0.23,0.41,0.59,0.77,0.95];
%% 参数设置
S1 = size(Pn_train,1); % 输入层神经元个数
S2 = 6; % 隐含层神经元个数
S3 = size(Tn_train,1); % 输出层神经元个数
Gm=10; %最大迭代次数
F0=0.5; %F为缩放因子
Np=5; %种群规模
CR=0.5; %杂交参数
G=1;%初始化代数
N=S1*S2 + S2*S3 + S2 + S3;%所求问题的维数
%% 创建BP神经网络
%net_optimized = newff(Pn_train,Tn_train,S2);
net_optimized=newff(minmax(Pn_train),[6,1],{'logsig','purelin'},'traingdm');%隐含层神经元S型正切,输出层S型对数,动量梯度下降法训练BP网络,
%%差分进化算法
ge=zeros(1,Np);%各代的最优值
bestx=zeros(Np,N);%各代的最优解
%产生初始种群
xmin=-1;xmax=1;
X0=(xmax-xmin)*rand(Np,N)+xmin;
X=X0;
%%%%%%%%%%变异操作
X1new=zeros(Np,N);%初始化
X1_new=zeros(Np,N);%初始化
X1=zeros(Np,N);%初始化
value=zeros(1,Np);
while G<=Gm
for i=1:Np
%产生j,k,p三个不同的数
a=1;
b=Np;
dx=randperm(b-a+1)+a-1;%randperm()随机打乱一个数字序列
j=dx(1);
k=dx(2);
p=dx(3);
if j==i
j=dx(4);
elseif k==i
k=dx(4);
elseif p==i
p=dx(4);
end
%%%%%%%%%%%%%%%变异操作
namd=exp(1-Gm/(Gm+1-G));%变异算子
F=F0*2.^namd;
bon=X(p,:)+F*(X(j,:)-X(k,:));
if (bon>xmin)&(bon<xmax) %防止变异超出边界
X1new(i,:)=bon;
else X1new(i,:)=(xmax-xmin)*rand(1,N)+xmin;
end
%%%%%%%%%%%%%%%交叉操作
if rand>CR %利用二项分布来交叉
X1_new(i,:)=X(i,:);
else
X1_new(i,:)=X1new(i,:);
end
%%%%%%%%%%%%%%%选择操作
if fun(X1_new(i,:),S1,S2,S3,net_optimized,Pn_train,Tn_train)<fun(X1(i,:),S1,S2,S3,net_optimized,Pn_train,Tn_train)
X1(i,:)=X1_new(i,:);
value(i)=fun(X1(i,:),S1,S2,S3,net_optimized,Pn_train,Tn_train);
else
X1(i,:)=X(i,:);
end
end
[fmin,nmin]=min(value);
ge(G)=fmin;
bestx(G,:)=X1(nmin,:);
G=G+1;
X=X1;
end
[gmin,n]=min(ge);
bestvalue=gmin
bestsolution=bestx(n,:)
%% 解码最优个体
x = bestsolution;
% 前S1*S2个编码为W1
temp1 = x(1:S1*S2);
W1 = reshape(temp1,S2,S1);
% 接着的S2*S3个编码为W2
temp2 = x(S1*S2+1:S1*S2+S2*S3);
W2 = reshape(temp2,S3,S2);
% 接着的S2个编码为B1
temp3 = x(S1*S2+S2*S3+1:S1*S2+S2*S3+S2);
B1 = reshape(temp3,S2,1);
%接着的S3个编码B2
temp4 = x(S1*S2+S2*S3+S2+1:S1*S2 + S2*S3 + S2 + S3);
B2 = reshape(temp4,S3,1);
% 设置网络初始权值和阈值
net_optimized.IW{1,1} = W1;
net_optimized.LW{2,1} = W2;
net_optimized.b{1} = B1;
net_optimized.b{2} = B2;
% 设置训练参数
net_optimized.trainParam.epochs = 3000;
net_optimized.trainParam.show = 100;
net_optimized.trainParam.goal = 0.001;
net_optimized.trainParam.lr = 0.1;
% 利用新的权值和阈值进行训练
net_optimized = train(net_optimized,Pn_train,Tn_train);
%% 仿真测试
Tn_sim_optimized = sim(net_optimized,Pn_test);
% 结果对比
result_optimized = [Tn_test' Tn_sim_optimized'];
%均方误差
E_optimized = mse(Tn_sim_optimized - Tn_test)
MAPE_optimized = mean(abs(Tn_sim_optimized-Tn_test)./Tn_sim_optimized)*100
% figure(1)
%
% plot(T_train,P_train(1,:),'r')
% hold on
% plot(T_train,P_train(3,:),'y')
% hold on
% plot(T_train,P_train(5,:),'b')
% hold on
% grid on
% xlabel('标准设备的约定真值(10KP)');
% ylabel('压力传感器的输出(mv)');
% title('压力传感器的工作曲线');
% legend('t=22','t=44','t=70');
figure(2)
plot(Tn_train(1:6),Pn_train(1,1:6),'r')
hold on
plot(Tn_train(7:12),Pn_train(1,7:12),'y')
hold on
plot(Tn_train(13:18),Pn_train(1,13:18),'b')
hold on
grid on
xlabel('设备约定真值(10KP)');
ylabel('压力传感器的输出(mv)');
title('归一化后的训练样本压力传感器的工作曲线');
legend('t=22','t=44','t=70');
figure(3)
plot(Tn_test(1:6),Pn_test(1,1:6),'r')
hold on
plot(Tn_test(7:12),Pn_test(1,7:12),'y')
hold on
plot(Tn_test(13:18),Pn_test(1,13:18),'b')
hold on
grid on
xlabel('设备约定真值(10KP)');
ylabel('压力传感器的输出(mv)');
title('归一化后的测试样本压力传感器的工作曲线');
legend('t=22','t=44','t=70');
figure(4)
plot(Tn_test(1:6),Tn_sim_optimized(1:6),'r')%输出DE-BP仿真结果的曲线
hold on
plot(Tn_test(7:12),Tn_sim_optimized(7:12),'y')
hold on
plot(Tn_test(13:18),Tn_sim_optimized(13:18),'b')
hold on
xlabel('约定真值(10KP)');
ylabel('压力传感器的输出(mv)');
title('DE-BP的压力传感器的工作曲线');
legend('t=22','t=44','t=70');
grid on
%% 未优化的BP神经网络
%net = newff(Pn_train,Tn_train,S2);
net=newff(minmax(Pn_train),[6,1],{'logsig','purelin'},'traingdm');%隐含层神经元S型正切,输出层S型对数,动量梯度下降法训练BP网络,
% 设置训练参数
net.trainParam.epochs = 3000;
net.trainParam.show = 100;
net.trainParam.goal = 0.001;
net.trainParam.lr = 0.1;
net=init(net);
inputWeights=net.IW{1,1};% 当前输入层权值和阈值
inputbias=net.b{1};
layerWeights=net.LW{2,1};% 当前网络层权值和阈值
layerbias=net.b{2}
% 利用新的权值和阈值进行训练
net = train(net,Pn_train,Tn_train);
%% 仿真测试
Tn_sim = sim(net,Pn_test);
%% 结果对比
result = [Tn_test' Tn_sim'];
% 均方误差
E1 = mse(Tn_sim - Tn_test)
MAPE1= mean(abs(Tn_sim-Tn_test)./Tn_sim)*100
% end
% figure(4)
% plot(T_train,P_train(1,:),'r')
% hold on
% plot(T_train,P_train(3,:),'y')
% hold on
% plot(T_train,P_train(5,:),'b')
% hold on
% grid on
% xlabel('标准设备的约定真值(10KP)');
% ylabel('压力传感器的输出(mv)');
% title('压力传感器的工作曲线');
% legend('t=22','t=44','t=70');
figure(5)
plot(Tn_test(1:6),Tn_sim(1:6),'r')%输出BP仿真结果的曲线
hold on
plot(Tn_test(7:12),Tn_sim(7:12),'y')
hold on
plot(Tn_test(13:18),Tn_sim(13:18),'b')
hold on
xlabel('约定真值(10KP)');
ylabel('压力传感器的输出(mv)');
title('BP的压力传感器的工作曲线');
legend('t=22','t=44','t=70');
grid on