clear all
close all
clc

%% Impostazione formato dati ingresso
% questa routine di validazione è applicabile nel caso si abbiano a
% disposizione tutte le matrici del sistema identificato e i dati. 

% in funzione di come ciascuno ha formattato i dati, si possono usare i
% dati del modello estivo e quelli invernali/autunnali.

%SCEGLIERE I DATI DI VALIDAZIONE
%esempio: d1 = iddata(y_val',u_val');
load d_val;
%load d_valinv
%load d_valaut

%MODIFICA IN FUNZIONE DI QUELLO CHE SI CARICA
d1=d_val;

%NON MODIFICARE: vengono dalla routine identificazione.m
%esempio: d_id = iddata(y_id',u_id');
load d_id;
%esempio: m1 = idss(A,B,C,D,K,x0);
load m_idss;
%esempio: model=ss(A,B,C,0,1);
load m_ss;

%ricavo tutte le informazioni necessarie. Non modificare!!
model=m_ss;
m1=m_idss;
u_val=d1.u;
y_val=d1.y;
u_id=d_id.u;
y_id=d_id.y;
A=m_idss.A;
B=m_idss.B;
C=m_idss.C;
D=m_idss.D;
K=m_idss.K;
x0=m_idss.x0;

%% 1° metodo validazione: valutazione visiva

%simulazione con il modello identificato con gli ingressi di validazione e
%poi confronteremo le uscite.
y_sim = lsim(model,u_val,1:length(u_val),x0);

%trasposta della matrice y_val per adattarne le dimensioni
%y_val=y_val';    

%plottaggi
figure
subplot(2,1,1)
plot(y_val);
grid on;
title('Uscite di validazione (da simulatore)');
ylabel('°C (scalati)');
axis([0 size(y_val,1) min(min(y_val))-0.5 max(max(y_val))+0.5])
subplot(2,1,2)
plot(y_sim);
grid on;
title('Uscite del modello identificato con gli ingressi di validazione');
ylabel('°C (scalati)');
axis([0 size(y_val,1) min(min(y_sim))-0.5 max(max(y_sim))+0.5])
print(['-f',num2str(gcf)], '-depsc2', 'uscite_validazione')

%% 2° metodo validazione: errore quadratico medio 
% tra le uscite di validazione e quelle in uscita dal modello identificato con u=u_val

for i= 1:size(y_val,2)
    MSEaperta(:,i)=(sum((y_sim(:,i)-y_val(:,i)).^2))/size(y_val,1);
end

%plotto e2
figure
stem(MSEaperta);
grid on;
axis([0 size(y_val,2)+1 0 1.1*max(MSEaperta)])
title('MSE catena aperta');
xlabel('n° uscite');
%ylabel('°C'); %non sono gradi è un indice
print(['-f',num2str(gcf)], '-depsc2', 'MSE_catena_aperta')

%% 3° metodo validazione: compare()

figure
    %compare() confronta l'n-esima uscita reale con la relativa prodotta dal
    %modello identificato con medesimi ingressi. é un altro modo per
    %validare. Si può usare anche il predittore a k-passi, ma non lo uso.
    compare(d1,m1);
    print(['-f',num2str(gcf)], '-depsc2', 'compare_catena_aperta')
    %evitando il plottaggio si salvano i dati di fit per relazione.
    [yh,fit_openloop,x00] = compare(d1,m1);

%% 4° metodo validazione: auto-correlazione e cross-correlazione dei residui
    
figure 
    %residui: Residual analysis consists of two tests: the whiteness test and the
    %independence test. Utilizza la funzione pe.m! già usata prima
    disp('CLICCA INVIO per analizzare la auto-corr. e la cross-corr. di tutte le uscite')
      %risultati diversi se si imposta il predittore a k passi in resid(m1,d1,k)
      %test da fare sui dati di identificazione?        
    resid(m1,d1);
    
%% mappa zero-polo incrociata tutti ingressi-tutte uscite

 figure
%plotta tutti insieme
pzmap(m1);
%plotta una fdt per volta
%pzmap(m1,'mode','same');
 
%% predizione a 1 passo con il cambio sistema e sua innovazione
% come conseguenza di un passaggio da un modello di stato
% x(t+1)=Ax(t)+Bu(t)+Ke(t) a uno normale x(t+1)=A*x(t)+Bu(t), ma con
% A*=A-K*C

y_ident2 = lsim(model,u_val,1:length(u_val),x0);%using data saved for validation
model_k=ss(A-K*C,[B K],C,0,1);
y_ident2k = lsim(model_k,[u_val y_val],1:length(u_val),x0);%using data saved for validation
model_e=ss(A-K*C,[B K],-C,[D eye( size(y_val,2))],1);
e_inn = lsim(model_e,[u_val y_val],1:length(u_val),x0);%errore di innovazion

figure()
plot([y_val(:,1) y_ident2k(:,1)])
axis([0 size(y_val,1) min(y_val(:,1))-0.5 max(y_val(:,1))+0.5])
grid on
title('Confronto prima uscita di validazione con la predizione a un passo');
legend('uscita di validazione','uscita predittore','Location','SouthEast')
xlabel('campioni');
ylabel('°C');
print(['-f',num2str(gcf)], '-depsc2', 'validaz_Vs_pred_a_un_passo')

%plot multiplo
figure()
plot(e_inn)
axis([0 size(y_val,1) min(min(e_inn))-0.01 max(max(e_inn))+0.01])
grid on
title('Andamenti errore di innovazione');
xlabel('campioni');
ylabel('°C');
print(['-f',num2str(gcf)], '-depsc2', 'andamenti_err_innovazione')

%plot singolo
figure()
plot(e_inn(:,1))
axis([0 size(y_val,1) min(e_inn(:,1))-0.01 max(e_inn(:,1))+0.01])
grid on
title('Andamento errore di innovazione prima uscita');
xlabel('campioni');
ylabel('°C');
print(['-f',num2str(gcf)], '-depsc2', 'andamento_errore_innovazione_uscita 1')

%stima della covarianza per serie di dati
cov_e_inn=covf(e_inn,length(e_inn));
figure() 
plot(cov_e_inn')
grid on
title('Andamenti della stima della covarianza di errore di innovazione');
print(['-f',num2str(gcf)], '-depsc2', 'covarianza_errore_innovazione')

%plot singolo
figure() 
plot(cov_e_inn(1,:)')
title('Covarianza di errore di innovazione prima uscita');
grid on
print(['-f',num2str(gcf)], '-depsc2', 'covarianza_errore_innovazione_uscita 1')


%errore quadratico medio delle uscite con predittore
for i= 1:size(y_val,2)
    MSEe_inn(:,i)=(sum((y_ident2k(:,i)-y_val(:,i)).^2))/size(y_val,1);
end

figure
stem(MSEe_inn);
grid on;
axis([0 size(y_val,2)+1 0 1.1*max(MSEe_inn)])
title('MSE uscite del predittore a un passo');
xlabel('n° uscite');
%ylabel('°C');
print(['-f',num2str(gcf)], '-depsc2', 'MSE_pred_a_un_passo')

%% predizione k passi

%seleziona istanti di predizione:
k=9;

% si tenga presente che usando la funzione compare effettua automaticamente
% la predizione con il filtro di kalman
% % % % yhat1=predict(m1,dd,k,'estimate');

figure
compare(d1,m1,k,'Init','e')
print(['-f',num2str(gcf)], '-depsc2', 'compare_k_passi')
%salvataggio fit
[y_pred,fit_predict,x00] = compare(d1,m1,k,'Init','e');
y_predkk=y_pred{1}.y(:,:);

% %plottaggio confronto validazione Vs predittore k passi di tutte le
% %uscite!!
% for i=1:size(y_val,2)
%     figure
%     plot(y_predk(:,i),'r')
%     hold on
%     grid on
%     plot(y_val(:,i))
%     legend('predittore','validazione','Location','SouthEast')
%     title(['Confronto validaz. Vs predittore a ' int2str(k) ' passi per la ' int2str(i) '° uscita '])
% end

%% plottaggio andamenti dei fit al variare di k passi del predittore per
%% ogni uscita
clear fit_tot
clear y_pred_tot

%moltiplicatore passi da impostare a piacere (introdotto per evitare di
%ricavare n predittori per n passi, bastano n/l predizioni e plottarne
%l'interpolazione per visualizzarne la tendenza).  
l=7;

for j=1:l
    [y_predk,fit_predictk,x00k] = compare(d1,m1,j*5,'Init','e');
    fit_tot(j,:)=squeeze(fit_predictk(:,:,:))';
    y_pred_tot(:,:,j)= y_predk{1}.y;
end
clear j

figure
xT=5:(l*5-5)/(l-1):l*5;
plot(xT,fit_tot)
grid on
xlabel('k passi');
ylabel('fit [%]');
title('Andamenti fit uscite per k predittori');
legend('y1','y2','y3','y4','y5','y6','y7','y8','y9','Location','SouthWest')
print(['-f',num2str(gcf)], '-depsc2', 'fit_k_volte')

%plottaggio diversi predittori per uscita a scelta, impostabile in pred_k..
pred_k=8;

figure
    plot(y_val(:,pred_k),'k','LineWidth',1.5)  
    grid on
    hold on
    yy= squeeze(y_pred_tot(:,pred_k,[ceil(l*(1/3)) ceil(l*(2/3)) l]));
    plot(yy)
xlabel('n° campioni');
ylabel('°C (scalati)');
title(['Confronto validazione e predittori per la ' int2str(pred_k) ' uscita']);
legend('uscita di validazione', ['predittore a ' int2str(ceil(l*(1/3))*5) ' passi'], ...
    ['predittore a ' int2str(ceil(l*(2/3))*5) ' passi'], ['predittore a ' int2str(l*5) ' passi'], 'Location','NorthWest')
print(['-f',num2str(gcf)], '-depsc2', 'k_pred_1y')

%stampo singolo confronto della figure w°=..
w=1;
figure()
plot(y_predkk(:,w),'r')
hold on
grid on
plot(y_val(:,w))
axis([0 size(y_val,1) min(min(y_val))-0.2 max(max(y_val))+0.2])
title(['Confronto validaz. Vs predittore a ' int2str(k) ' passi per la ' int2str(w) '° uscita '])
legend('dati validazione','dati predittore','Location','SouthEast')
print(['-f',num2str(gcf)], '-depsc2', 'validaz_Vs_pred_k_passi')

%andamento errore predittore a k passi per tutte uscite
figure()
err_pred1=y_val-y_predkk;
plot(err_pred1)
grid on
title(['Andamenti errori di predizione a ' int2str(k) ' passi'])
axis([0 size(y_val,1) min(min(err_pred1))-0.01 max(max(err_pred1))+0.01])
print(['-f',num2str(gcf)], '-depsc2', 'errori_a_k_passi')

for i= 1:size(y_val,2)
    MSEerr_pred(i)=(sum((y_val(:,i)-y_predkk(:,i)).^2))/size(y_val,1);
end

%plotto MSE k passi
figure()
stem(MSEerr_pred);
grid on;
xlabel('n° uscite');
%ylabel('°C');
title('err_pred');
axis([0 size(y_val,2)+1 0 1.1*max(MSEerr_pred)])
title(['MSE del predittore a ' int2str(k) ' passi']);
print(['-f',num2str(gcf)], '-depsc2', 'MSE_k_passi')

%confronto, tra MSE in catena aperta e il predittore
figure()
plot(MSEaperta)
hold on
plot(MSEerr_pred,'r')
title(['MSE tra predittore a ' int2str(k) ' passi Vs catena aperta']);
legend('catena aperta','errore predittore','Location','NorthWest')
grid on
print(['-f',num2str(gcf)], '-depsc2', 'MSE_validaz_Vs_errore_pred_k_passi')

%% considerazioni fpe e akaike

%fpe confronta vari modelli. Dovrei mettere più modelli d in fpe().
%N.B: fpe() e aic() funzionano con modelli che derivano direttamente
%dall'identificazione come "m." n4sid. Il tool di chiuso non fornisce
%questo.
    %fp = fpe(m);
%Akaike Information Criterion for estimated model: il modello con valore di
%aic più piccolo è il migliore.
    %Aic = aic(m);    
    
    
%% %%%%%%%altre tecniche per le caratteristiche del sistema%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% matrice Hankel (per oggetti ssdata)
% %In control theory, eigenvalues define a system stability, whereas Hankel singular 
% %values define the "energy" of each state in the system. Keeping larger energy states 
% %of a system preserves most of its characteristics in terms of stability, frequency, 
% %and time responses. Model reduction techniques presented here are all based on the Hankel 
% %singular values of a system. They can achieve a reduced-order model that preserves 
% %the majority of the system characteristics.
% figure
% hankelsv(model);
% print(['-f',num2str(gcf)], '-depsc2', 'Hankel')

%% plottaggio valori singolari di un sistema LTI (che non è il nostro caso)
% figure
% sigma(model);
% grid on
% title('valori singolari di un sistema LTI');

%% Feedback
%     
% %Questa prodedura verifica la presenza di feedback dei dati (u,y)
% %una percentuale alta sopra al 90% equivale alla presenza di feedback!
% %(non sarebbe corretto fare la media semplice!)
% 
% [fbck,fbck0,nudir] = feedback(d1);
% for k=1:size(fbck,2)
%     if mean(fbck(k))>90
%         disp(['PRESENZA DI REAZIONE (in media) per ingresso: ' int2str(k)])
%         else disp(['ASSENZA DI REAZIONE (in media) per ingresso: ' int2str(k)])
%     end
% end
% 
%% advice.m fornisce informazioni sul iddata dei dati di validazione.
% 
% advice(d1)
% 
%% routine che mi dice il ritardo delle uscite nei confronti dell'ingresso
% for j=1:size(d1,2)
%     nk=delayest(squeeze(d1(:,j,:)));
%     disp(['Uscita: ' int2str(j)]);
%     disp(['affetta da un RITARDO (in campioni per ingresso) di: ' int2str(nk)])
% end