%% Laboratory 6 - Estimation, filtering and system identification - Prof. M. Taragna
% *Exercise: model identification for a water heater, using linear 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, clc

%% Procedure
% # Load the .mat file containing input-output data
% # Remove the mean value from the data
% # Plot input-output data
% # Define estimation and validation datasets
% # For each input-output delay nk from 1 to 5
% # And for each order na=nb (ARX) / nb=nf (OE) from 1 to 5
% # Identify the ARX/OE model using the estimation dataset
% # Perform the whiteness test on the estimation dataset
% # Compute the predicted output on the validation dataset
% # Compute the RMSE and MDL criteria
% # Compare RMSE and MDL values
% # Select the best trade-off between all the model classes

%% Problem 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')

% Step 4: definition of estimation and validation datasets

DATAe = [ye_0mean,ue_0mean]; % estimation dataset
DATAv = [yv_0mean,uv_0mean]; % validation dataset

%% ARX identification, residual test and possible validation
% To run this part in a smart way, you can put a break point at the 
% command "close all" at the end of the external for-loop. This way you can
% check in figure N the corrisponding model with na=nb=N. After that you 
% can click on "Continue" on the toolbar to evaluate for the next nk value.
% After the analysis you may comment the 'resid' command to speedup the 
% next run instance.

% Step 5: loop on input-output delay

close all
for nk = 1:5,
    
    % Step 6: loop on model orders

    for k = 1:5,
        
        % Step 7: ARX(na=k,nb=k,nk) model identification

        model = arx(DATAe,[k, k, nk]);
        
        % Step 8: computation of model residuals and whiteness test

        % figure, resid(DATAe, model, 'corr', 30);
    
        % Step 9: computation of the predicted output on validation dataset

        yh = compare(DATAv, model, 1); % III param: 1 -> prediction; inf -> simulation
        % figure, plot(1:Nv,yv_0mean,'-k',1:Nv,yh,'--r'),
        % title(['ARX(',num2str(k),',',num2str(k),',',num2str(nk),')']),
        % legend('Output data yv\_0mean','Predicted output yh');
        N0 = 10;

        % Step 10: RMSE and MDL computation on validation dataset
        
        MSE = 1 /(Nv-N0)*norm(yv_0mean(N0+1:end)-yh(N0+1:end))^2;
        RMSE_arx(nk,k) = sqrt(MSE);

        n = 2*k;
        n_arx(k)=n;
        MDL_arx(nk,k) = n*log(Nv-N0)/(Nv-N0)+log(MSE);
        
    end
    close all
end

%% ARX comparison
% The figures generated by resid show in the first part (AutoCorr) the residual values 
% and the confidence intervals. You have to check all the generated plots: 
% the more residual values are inside the confidence interval, the better is the model. 
% You count the number of residuals sufficiently outside the confidence interval 
% and you select a threshold: if this number is greater than the threshold, 
% then the model is wasted, otherwise it is considered for further analysis. 
% In this case, a reasonable threshold could be 3.
% The whiteness test results for ARX models are:
% 
%              na=1, na=2, na=3, na=4, na=5
%         nk=1:  10     2     2     1     2
%         nk=2:  10     2     2     1     2
%         nk=3:   7     2     1     1     1
%         nk=4:   4     4     1     1     0
%         nk=5:   5     3     1     0     1
%
% => the ARX models that sufficiently satisfy the residual whiteness test are:
%
% * For nk = 1, 2, 3, 5: na = nb >= 2
% * For nk = 4:          na = nb >= 3
%
% Then the performance criteria can be printed

% Step 11: ARX comparison

RMSE_arx, MDL_arx
figure, plot(n_arx,RMSE_arx','o-'), title('RMSE of ARX models'),
legend('ARX(nk=1)','ARX(nk=2)','ARX(nk=3)','ARX(nk=4)','ARX(nk=5)')
figure, plot(n_arx,MDL_arx','o-'), title('MDL of ARX models'),
legend('ARX(nk=1)','ARX(nk=2)','ARX(nk=3)','ARX(nk=4)','ARX(nk=5)')

%% OE identification, residual test and possible validation

% Step 5: loop on input-output delay

close all
for nk = 1:5,

    % Step 6: loop on model orders

    for k = 1:5,
        
        % Step 7: OE(nb=k,nf=k,nk) model identification

        model = oe(DATAe,[k, k, nk]);
    
        % Step 8: computation of model residuals and whiteness test

        % figure, resid(DATAe, model, 'corr', 30);
        
    end
    close all
end

%% OE comparison
% The whiteness test results for OE models are:
% 
%              nb=1, nb=2, nb=3, nb=4, nb=5
%        nk=1:   7    16    17    17    19
%        nk=2:   7    16    17    17    14
%        nk=3:  10    16    16    14    18
%        nk=4:  14    16    16    18    17
%        nk=5:  15    16    18    16    17
%        
% NO OE model passes the residual whiteness test with Threshold = 3 =>
% OE is not the right model class in this case: for any order and delay,
% the number of residuals outside the confidence interval is too high

%% Trade-off selection
% The goal is to minimize the RMSE and MDL criteria, taking into account 
% the previous residual test and the model complexity n, that for models
% in prediction is: n=na+nb (ARX), n=nb+nf (OE)
% In this case, the ARX model class only has to be considered and the best 
% trade-off between complexity and both the validation criteria is the 
% ARX(3,3,3)

% Step 12: trade-off selection

best_tradeoff = arx(DATAe,[3,3,3]);
yh = compare(DATAv, best_tradeoff, 1); % III param: 1 -> prediction; inf -> simulation
figure, plot(1:Nv,yv_0mean,'-k',1:Nv,yh,'--r'),
legend('Output data yv\_0mean','Predicted output yh',4);

RMSE_best_tradeoff = RMSE_arx(3,3)
MDL_best_tradeoff = MDL_arx(3,3)

present(best_tradeoff)
[A,B,C,D,F] = polydata(best_tradeoff)

% When used in prediction mode, the ARX model transfer functions are:
% - between u(t) and yh(t): B(z)
% - between y(t) and yh(t): 1-A(z)
% Pay care that the polynomials A(z) and B(z) are in the z^-1 variable!!