%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       利用BP拟合函数y=x1^2+x2^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
x1=5*rands(2000,1);
x2=5*rands(2000,1);
y=x1.^2+x2.^2;
input=[x1,x2];
output=y;
k=rand(1,2000);
[m,n]=sort(k);

%找出训练数据和预测数据
input_train=input(n(1:1900),:)';
output_train=output(n(1:1900))';
input_test=input(n(1901:2000),:)';
output_test=output(n(1901:2000))';

%选连样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);

%% BP网络训练
% %初始化网络结构
net=newff(inputn,outputn,20);

net.trainParam.epochs=2000;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00004;

%网络训练
net=train(net,inputn,outputn);
%预测数据归一化
inputn_test=mapminmax('apply',input_test,inputps);
%网络预测输出
an=sim(net,inputn_test);
%网络输出反归一化
BPoutput=mapminmax('reverse',an,outputps);
%误差分析
error=BPoutput-output_test;
rmse=sqrt(mean(error.^2));
%% 结果分析

figure(1)
[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);
z=x.^2+y.^2;
mesh(x,y,z);
title('函数图形','fontsize',12)

figure(2)
plot(BPoutput,':or')
hold on
plot(output_test,'-*');
text(2,47,['RMSE=',num2str(rmse,'%2.4f')]);
legend('预测输出','期望输出')
title('BP网络预测输出','fontsize',12)
ylabel('函数输出','fontsize',12)
xlabel('样本','fontsize',12)