clear

% 輸入歷史實測數據（34天）的汽輪機振動的日平均峰峰值
day=[40.6,37.2,36.2,36.8,36.4,36.8,35.6,37.2,35.2,35.8,36,35.9,39,40,40.4,39.2,36.6,35.8,34.8,35.8,37,38.2,41.2,40.2,37.6,37.8,39.8,40.8,38.8,38.6,42,47.2,46.6,44.2];


dayhistory=day(1:30);%取其中三十天作為歷史數據樣本
dayhismod=reshape(dayhistory,5,6);% 將歷史數據分為6個樣本，每個大小為5，其中reshape是以列排序的
dayday=day(1:25);% 取其中的前25天
daypost=day(6:30);%取其中的隨后25天
p=reshape(dayday,5,5);% 將前25天數據分為5行5列矩陣作為網絡的訓練輸入樣本
t=reshape(daypost,5,5); %將隨后的25天分為5行5列矩陣作為網絡的目標輸出向量
daylast=day(26:30);
h3=reshape(daylast,5,1);% 將倒數第二個樣本作為網絡測試時的輸入樣本
r=6:30;
rr=reshape(r,5,5);


%%%%%%%%%%%%%% 新建網絡bp %%%%%%%%%%%%%%%%
net=newff(minmax(p),[5,5],{'purelin' 'purelin'},'trainlm');
y1=sim(net,p);
% 新建網絡，其中minmax(p)為p的沒一次輸入的最大最小值向量
% 兩層的傳遞函數均為purelin
% 訓練函數為trainlm
% 所訓練的網絡大小為[5,5]
% 仿真訓練前的網絡

%%%%%%%%%%%  進行網絡訓練  %%%%%%%%%%%%%%
% network parameters:
%   epochs--epochs of the train
%   goal--errors goal of the network
%   lr--learning rate
%   shows--epochs between the displays
%   time--Maximum time to train in seconds
net.trainParam.epochs=200;  % 訓練次數
nettrainParam.goal=0.0001;  % 誤差期望值
% returns of the train:
%   net--New network
%    tr--Training record (epoch and perf).
%     Y--Network outputs.
%     E--Network errors.
[net,tr,Y,E]=train(net,p,t); 


%%%%%%%%%%%  網絡測試 %%%%%%%%%%%%%%%%
% input the testing points here %

y21=sim(net,p);
y2=reshape(y21,1,25);
clf
plot(r,y2,'b-^')
hold on
plot(1:34,day,'r-*')
%%%%%%%%%%%%% 預測 %%%%%%%%%%%%%%%
y3=sim(net,h3);
plot(31:35,y3,'-*')
hold on

title('神經網絡訓練結果');
xlabel('時間（天）');
ylabel('仿真輸出結果');
legend('仿真模擬值','實際值','神經網絡預測值');
%%%%%%%%%%%%%%%%%%  繪制誤差曲面 %%%%%%%%%%%%%
x=1:5;
y=1:5;
plot3(x,y,E(x,y))