%% Laboratory 6 - Estimation, filtering and system identification - Prof. M. Taragna
% *Exercise: model identification for a water heater, using nonlinear models (real data)*
%% Introduction
% The program code may be splitted in sections using the characters "%%". 
% Each section can run separately with the command "Run Section"
% (in the Editor toolbar, just to the right of the "Run" button). You can do the
% same thing by highlighting the code you want to run and by using the button 
% function 9 (F9). This way, you can run only the desired section of your code,
% saving your time. This script can be considered as a reference example.

clear all, close all

%% Procedure
% # Load the .mat file containing input-output data
% # Remove the mean value from the data
% # Plot input-output data
% # For each input/output delay nk from 1 to 5, and for each order na=nb
% (NNARX) from 1 to 5
% # Identify the NNARX model using the estimation dataset
% # Compute the predicted output on the validation dataset
% # Compute RMSE and MDL criteria in order to select the most suitable model

%% Setup

% Step 1: load of data

load heater

% Step 2: remove the mean value

ue_0mean = ue - mean(ue);
ye_0mean = ye - mean(ye);
uv_0mean = uv - mean(uv);
yv_0mean = yv - mean(yv);

% Step 3: plot of input-output data

Ne = length(ue); % number of estimation data
Nv = length(uv); % number of validation data

figure, plot(ue_0mean), ylabel('ue\_0mean'),
figure, plot(ye_0mean), ylabel('ye\_0mean'),
figure, plot(uv_0mean), ylabel('uv\_0mean'),
figure, plot(yv_0mean), ylabel('yv\_0mean')

%% Neural Network model class
%%
% Define the model structure of the neural network having:
%%
% * r>=5 hyperbolic tangent neurons (or nodes or units) in the hidden layer
% * 1 linear neuron (or node or unit) in the output layer

NetDef=['HHHHHH';'L-----'];

% For any Neural Network ARX model NNARX(na,nb,nk), define its orders:
% - na = # of past measurements used in the regressor phi
% - nb = # of past inputs used in the regressor phi
% - nk = minimum input-output delay

% Step 4: definition of NNARX parameters

na=1;
nb=na;
nk=1;

%%
% Set the maximum number of iterations to 500

trparms=settrain; trparms=settrain(trparms,'maxiter',500);

%%
% Perform the NNARX identification, using
%%
% * as ouput signal Y, a row vector with dim(Y) = [1 * (# of data)]
% * as input signal U, a matrix with dim(U) = [(# of inputs) * (# of data)]
% In the SISO case, both Y and U are row vectors

% Step 5: estimation of NNARX parameters

[W1,W2]=nnarx(NetDef,[na,nb,nk],[],[],trparms,ye_0mean',ue_0mean');

%%
% Plot the structure of the NNARX model
figure, drawnet(W1,W2)

%%
% Extract the estimated parameters of the NNARX model
alfa=W2(1:end-1);
alfa_0=W2(end);
beta=W1(:,1:end-1);
beta_0=W1(:,end);

%%
% Compute the predicted output of the NNARX model using the validation
% dataset, defining the regressor vector phi as a column

% Step 6: computation of predicted output

t_min=max([na+1,nb+nk]);
for t=t_min:Nv,
    phi= [yv_0mean(t-1:-1:t-na); uv_0mean(t-nk:-1:t-nk-nb+1)];
    yp(t)=alfa*tanh(beta*phi+beta_0)+alfa_0;
end

figure, plot(1:Nv,yv_0mean,'g', 1:Nv,yp,'r-.')

%%
% Compute the NNARX performances in terms of RMSE and MDL

% Step 7: computation of RMSE and MDL

N0=10; na,nb,nk
RMSE=sqrt(1/(Nv-N0))*norm(yv_0mean(N0+1:end)-yp(N0+1:end)')
n=na+nb;
MDL=n*log(Nv-N0)/(Nv-N0)+log(RMSE^2)


%% Choosing the best model
% The RMSE is minimized over 10 trials for each model (with N0=10):
%%
% * nnarx(1,1,1) => RMSE = 1.0128
% * nnarx(1,1,2) => RMSE = 0.9444
% * nnarx(1,1,3) => RMSE = 0.8250
% * nnarx(1,1,4) => RMSE = 0.6590
% * nnarx(1,1,5) => RMSE = 0.5705
% * nnarx(2,2,1) => RMSE = 0.6907
% * nnarx(2,2,2) => RMSE = 0.6414
% * nnarx(2,2,3) => RMSE = 0.5448
% * nnarx(2,2,4) => RMSE = 0.4894
% * nnarx(2,2,5) => RMSE = 0.5205
% * nnarx(3,3,1) => RMSE = 0.5983
% * nnarx(3,3,2) => RMSE = 0.5197
% * nnarx(3,3,3) => RMSE = 0.4757
% * nnarx(3,3,4) => RMSE = 0.4548    <=   best NNARX model
% * nnarx(3,3,5) => RMSE = 0.5160
% * nnarx(5,5,1) => RMSE = 0.4801
% * nnarx(5,5,4) => RMSE = 0.4600
% * nnarx(7,7,1) => RMSE = 0.4666
% * nnarx(7,7,4) => RMSE = 0.4557
% * nnarx(9,9,1) => RMSE = 0.4678
% * nnarx(9,9,4) => RMSE = 0.5459