load iddata-02

%   X - input time series.
%   T - feedback time series.

X = tonndata(id.u,false,false);
T = tonndata(id.y,false,false);

% Choose a Training Function
% 'trainbr' is best suited for difficult cases
trainFcn = 'trainbr'; %Bayesian Regulation backpropagation

inputDelays = 1:7;
feedbackDelays = 1:2;
hiddenLayerSize = 10;
net = narxnet(inputDelays,feedbackDelays,hiddenLayerSize,'open',trainFcn);
net.performFcn = 'mse';
% Prepare the Data for Training and Simulation
[x,xi,ai,t] = preparets(net,X,{},T);

% Setup Division of Data for Training, Validation, Testing
net.divideParam.trainRatio = 100/100;
net.divideParam.valRatio = 0/100;
net.divideParam.testRatio = 0/100;
net.trainParam.epochs = 2000

% Train the Network
[net,tr] = train(net,x,t,xi,ai);

% Test the Network
y = net(x,xi,ai);
e = gsubtract(t,y);
performance = perform(net,t,y)
figure, plotperform(tr)


%% Closed Loop Network
netc = closeloop(net);
[xc,xic,aic,tc] = preparets(netc,X,{},T);
yc = netc(xc,xic,aic);
closedLoopPerformance = perform(net,tc,yc)
figure, plot(id.y(inputDelays(end:end):end),'-r')
hold on
plot(cell2mat(yc),'-b');

% Step-Ahead Prediction Network

nets = removedelay(net);
[xs,xis,ais,ts] = preparets(nets,X,{},T);
ys = nets(xs,xis,ais);
stepAheadPerformance = perform(nets,ts,ys)
plot(cell2mat(ys),'-y');
legend('identification output','closed loop of Neural Network output', 'One-Step-Ahead Prediction output')
title(['Closed Loop MSE: ' num2str(closedLoopPerformance), ' Prediction MSE: ' num2str(stepAheadPerformance)]);
grid
hold off
shg

%% Validation set

%   Xval - input time series.
%   Tval - feedback time series.

Xval = tonndata(val.u,false,false);
Tval = tonndata(val.y,false,false);

% Closed Loop Network

netcVal = closeloop(net);
[xc,xic,aic,tc] = preparets(netcVal,Xval,{},Tval);
ycval = netcVal(xc,xic,aic);
closedLoopPerformance = perform(netcVal,tc,ycval)
figure, plot(val.y(inputDelays(end:end):end),'-r')
hold on
plot(cell2mat(ycval),'-b');

% Step-Ahead Prediction Network

netsVal = removedelay(net);
[xs,xis,ais,ts] = preparets(netsVal,Xval,{},Tval);
ysval = netsVal(xs,xis,ais);
stepAheadPerformance = perform(netsVal,ts,ysval)
plot(cell2mat(ysval),'-g');
legend('validation output','closed loop of Neural Network output', 'One-Step-Ahead Prediction output')
title(['Closed Loop MSE: ' num2str(closedLoopPerformance), ' Prediction MSE: ' num2str(stepAheadPerformance)]);
grid
hold off
shg