%% Laboratory 5 - Estimation, filtering and system identification - Prof. M. Taragna
% *Exercise 2: identification of a model for a dryer system*
%% 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 file |dryer2.mat| containing the input/output data
% # Plot the input |u2| and the output |y2|
% # Remove the mean value from the data (with the function mean) 
% in order to remove any systematic error
% # Verify that it is correct with the function |mean| (it has to be zero)
% # Split the data in two parts: one for the parameter estimation and 
% one for the model validation
% # Define |Ze| and |Zv| as |[ye,ue]| and |[yv,uv]|, respectively
% # For each input/output delay |nk| from 1 to 4
% # And for each order |na=nb| (ARX) / |na=nb=nc| (ARMAX) / |nb=nf| (OE) from 1 to 7
% # Create the ARX/ARMAX/OE model with the function
% |arx(Ze,[na,nb,nk])| / |armax(Ze,[na,nb,nc,nk])| / |oe(Ze,[nb,nf,nk])|
% # Perform the whiteness test with the |resid(Ze,model,'CORR',30)| function
% # Compute the predicted output, with the |compare(Zv,model,1)| function, or
% the simulated output, with the |compare(Zv,model,inf)| function, and save it in |yh|
% # Compute the criteria RMSE, FPE, etc... in order to select the
% most suitable model. You can find the formulas on the official formulary 
% or on the slides.

%% Problem setup

% Step 1: load the data

load dryer2

% Step 2: plot the data

figure, 
subplot(2,1,1), plot(u2), grid on, ylabel('input u2(t)')
subplot(2,1,2), plot(y2), grid on, ylabel('output y2(t)')

% Step 3: remove the mean values

u = u2 - mean(u2);
y = y2 - mean(y2);

% Step 4: verify the mean values of u and y

mean_u = mean(u)
mean_y = mean(y)

% Step 5: partition the data

Ntot = length(u);
Ne = Ntot/2;
Nv = Ntot-Ne;

ue = u(1:Ne);
uv = u(Ne+1:end);
ye = y(1:Ne);
yv = y(Ne+1:end);

figure, 
subplot(2,1,1), plot(1:Ntot,u,'b', 1:Ne,ue,'g--', Ne+1:Ntot,uv,'r--'),
grid on, ylabel('input u(t)')
subplot(2,1,2), plot(1:Ntot,y,'b', 1:Ne,ye,'g--', Ne+1:Ntot,yv,'r--'),
grid on, ylabel('output y(t)')

% Step 6: define the datasets

Ze = [ye, ue]; % Estimation dataset
Zv = [yv, uv]; % Validation dataset

%% ARX model class
% 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.

close all
% Step 7
for nk = 1:4, % loop on the input-output delay
    % Step 8
    for na = 1:7, % loop on the order na
        nb = na;
        % Step 9: estimate the ARX parameters
        model = arx(Ze,[na,nb,nk]);
        
%         % Step 10: check the residual whiteness on the Estimation dataset
%         figure, resid(Ze,model,'CORR',30);
        
        % Step 11: compute the PREDICTED output on the Validation dataset
        yh = compare(Zv,model,1); % III param: 1 -> prediction; inf -> simulation
        % figure, plot(1:Nv,yv,'g', 1:Nv,yh,'r')
        
        % Step 12: compute the validation criteria
        N0 = 10;
        MSE = 1/(Nv-N0)*norm(yv(N0+1:end)-yh(N0+1:end))^2;
        n = na+nb; % model complexity for ARX models in prediction mode
        n_arx(na) = n;
        FPE_arx(nk,na) = (Nv-N0+n)/(Nv-N0-n)*MSE;
        AIC_arx(nk,na) = n*2/(Nv-N0)+log(MSE);
        MDL_arx(nk,na) = n*log(Nv-N0)/(Nv-N0)+log(MSE);
        FIT_arx(nk,na) = 1-sqrt(MSE/(1/(Nv-N0)*norm(yv(N0+1:end)-mean(yv(N0+1:end)))^2));
        if nk==3 & na==2, best_arx = model;  end
    end
%     pause, close all
end

%%
% The generated figures 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 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 analyses. 
% A reasonable threshold should be 4 or 5.
%
% The models that sufficiently satisfy the residual whiteness test are:
%% 
% * For nk = 1: na = nb >= 2
% * For nk = 2, 3, 4: na = nb >= 1
%
% Then you can print the performance criteria.

FPE_arx, AIC_arx, MDL_arx, FIT_arx
figure, plot(n_arx,FPE_arx(1,:),'r',n_arx,FPE_arx(2,:),'g',n_arx,FPE_arx(3,:),'b',n_arx,FPE_arx(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('FPE for ARX models')
figure, plot(n_arx,AIC_arx(1,:),'r',n_arx,AIC_arx(2,:),'g',n_arx,AIC_arx(3,:),'b',n_arx,AIC_arx(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('AIC for ARX models')
figure, plot(n_arx,MDL_arx(1,:),'r',n_arx,MDL_arx(2,:),'g',n_arx,MDL_arx(3,:),'b',n_arx,MDL_arx(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('MDL for ARX models')
figure, plot(n_arx,FIT_arx(1,:),'r',n_arx,FIT_arx(2,:),'g',n_arx,FIT_arx(3,:),'b',n_arx,FIT_arx(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('BEST FIT for ARX models')

%%
% The goal is to minimize the criteria FPE, AIC and MDL and to maximize the
% best FIT, taking into account the complexity of the model, that for the 
% ARX model is na+nb, and the previous residual test.
%
% The best trade-off between complexity and any validation criterion is 
% the model ARX(2,2,3):

present(best_arx)
[A,B,C,D,F] = polydata(best_arx);
A = tf(A,1,-1,'Variable','z^-1');
B = tf(B,1,-1,'Variable','z^-1');
TF_arx_pred_u = B
TF_arx_pred_y = 1-A

%% ARMAX model class
% 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=nc=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.

close all
% Step 7
for nk = 1:4, % loop on the input-output delay
    % Step 8
    for na = 1:7, % loop on the order na
        nb = na;
        nc = na;
        % Step 9: estimate the ARMAX parameters
        model = armax(Ze,[na,nb,nc,nk]);
        
%         % Step 10: check the residual whiteness on the Estimation dataset
%         figure, resid(Ze,model,'CORR',30);
        
        % Step 11: compute the PREDICTED output on the Validation dataset
        yh = compare(Zv,model,1); % III param: 1 -> prediction; inf -> simulation
        % figure, plot(1:Nv,yv,'g', 1:Nv,yh,'r')
        
        % Step 12: compute the validation criteria
        N0 = 10;
        MSE = 1/(Nv-N0)*norm(yv(N0+1:end)-yh(N0+1:end))^2;
        n = na+nb+nc; % model complexity for ARMAX models in prediction mode
        n_armax(na) = n;
        FPE_armax(nk,na) = (Nv-N0+n)/(Nv-N0-n)*MSE;
        AIC_armax(nk,na) = n*2/(Nv-N0)+log(MSE);
        MDL_armax(nk,na) = n*log(Nv-N0)/(Nv-N0)+log(MSE);
        FIT_armax(nk,na) = 1-sqrt(MSE/(1/(Nv-N0)*norm(yv(N0+1:end)-mean(yv(N0+1:end)))^2));
        if nk==3 & na==2, best_armax = model; end
    end
%     pause, close all
end

%%
% The generated figures 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 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 analyses. 
% A reasonable threshold should be 4 or 5.
%
% The models that sufficiently satisfy the residual whiteness test are:
%%
% * For nk = 1, 2, 3, 4: na = nb = nc >= 1
%
% Then you can print the performance criteria.

FPE_armax, AIC_armax, MDL_armax, FIT_armax
figure, plot(n_armax,FPE_armax(1,:),'r',n_armax,FPE_armax(2,:),'g',n_armax,FPE_armax(3,:),'b',n_armax,FPE_armax(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('FPE for ARMAX models')
figure, plot(n_armax,AIC_armax(1,:),'r',n_armax,AIC_armax(2,:),'g',n_armax,AIC_armax(3,:),'b',n_armax,AIC_armax(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('AIC for ARMAX models')
figure, plot(n_armax,MDL_armax(1,:),'r',n_armax,MDL_armax(2,:),'g',n_armax,MDL_armax(3,:),'b',n_armax,MDL_armax(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('MDL for ARMAX models')
figure, plot(n_armax,FIT_armax(1,:),'r',n_armax,FIT_armax(2,:),'g',n_armax,FIT_armax(3,:),'b',n_armax,FIT_armax(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('BEST FIT for ARMAX models')

%%
% The goal is to minimize the criteria FPE, AIC and MDL and to maximize the
% best FIT, taking into account the complexity of the model, that for the 
% ARMAX model is na+nb+nc, and the previous residual test.
%
% The best trade-off between complexity and any validation criterion is 
% the model ARXMAX(2,2,2,3):

present(best_armax)
[A,B,C,D,F] = polydata(best_armax)
A = tf(A,1,-1,'Variable','z^-1');
B = tf(B,1,-1,'Variable','z^-1');
C = tf(C,1,-1,'Variable','z^-1');
TF_armax_pred_u = B/C
TF_armax_pred_y = 1-A/C

%% OE model class
% 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 nb=nf=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.

close all
% Step 7
for nk = 1:4, % loop on the input-output delay
    % Step 8
    for nf = 1:7, % loop on the order nf
        nb = nf;
        % Step 9: estimate the OE parameters
        model = oe(Ze, [nb,nf,nk]);
        
%         % Step 10: check the residual whiteness on the Estimation dataset
%         figure, resid(Ze, model, 'CORR', 30');
        
        % Step 11: compute the PREDICTED output on the Validation dataset
        yh = compare(Zv,model,1); % III param: 1 -> prediction; inf -> simulation
        % figure, plot(1:Nv,yv,'g', 1:Nv,yh,'r')
        
        % Step 12: compute the validation criteria
        N0 = 10;
        MSE=1/(Nv-N0)*norm(yv(N0+1:end)-yh(N0+1:end))^2;
        n = nb+nf; % model complexity for OE models in prediction mode
        n_oe(nf) = n;
        FPE_oe(nk,nf) = (Nv-N0+n)/(Nv-N0-n)*MSE;
        AIC_oe(nk,nf) = n*2/(Nv-N0)+log(MSE);
        MDL_oe(nk,nf) = n*log(Nv-N0)/(Nv-N0)+log(MSE);
        FIT_oe(nk,nf) = 1-sqrt(MSE/(1/(Nv-N0)*norm(yv(N0+1:end)-mean(yv(N0+1:end)))^2));
    end
%     pause, close all
end
%%
% The generated figures 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 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 analyses. 
% A reasonable threshold should be 4 or 5.
%
% The models that sufficiently satisfy the residual whiteness test are:
%%
% * For nk = 2, 3: na = nb = nc = 1
%
% OE is not a good model class for this example: for almost all the orders,
% the number of residuals outside the confidence interval is high.
%
% Then you can print the performance criteria.

FPE_oe, AIC_oe, MDL_oe, FIT_oe
figure, plot(n_oe,FPE_oe(1,:),'r',n_oe,FPE_oe(2,:),'g',n_oe,FPE_oe(3,:),'b',n_oe,FPE_oe(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('FPE for OE models')
figure, plot(n_oe,AIC_oe(1,:),'r',n_oe,AIC_oe(2,:),'g',n_oe,AIC_oe(3,:),'b',n_oe,AIC_oe(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('AIC for OE models')
figure, plot(n_oe,MDL_oe(1,:),'r',n_oe,MDL_oe(2,:),'g',n_oe,MDL_oe(3,:),'b',n_oe,MDL_oe(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('MDL for OE models')
figure, plot(n_oe,FIT_oe(1,:),'r',n_oe,FIT_oe(2,:),'g',n_oe,FIT_oe(3,:),'b',n_oe,FIT_oe(4,:),'k'),
legend('nk = 1', 'nk = 2', 'nk = 3', 'nk = 4'), grid on, title('BEST FIT for OE models')

%% Choosing the best model

%%
% FPE (to be minimized)
FPE_arx(3,2)
FPE_armax(3,2)
%%
% ARMAX(2,2,2,3) has slightly better performance but much higher complexity than ARX(2,2,3)

%%
% AIC (to be minimized)
AIC_arx(3,2)
AIC_armax(3,2)
%%
% ARMAX(2,2,2,3) has slightly better performance but much higher complexity than ARX(2,2,3)
%%
% MDL (to be minimized)
MDL_arx(3,2)
MDL_armax(3,2)
%%
% ARMAX(2,2,2,3) has slightly better performance but much higher complexity than ARX(2,2,3)
%%
% FIT (to be maximized)
FIT_arx(3,2)
FIT_armax(3,2)
%%
% ARMAX(2,2,2,3) has slightly better performance but much higher complexity than ARX(2,2,3)
%%
% Conclusion: ARX(2,2,3) is the best trade-off between performances and complexity