clear all
close all
clc

%% Recupero dati, selezione tipologia e adattamento

%carico i dati da DesignBuilder, con data e ora dei campioni sulle righe,
%e i nomi delle tipologie dei dati in colonna.
load dati_inverno
          
%tolgo alle uscite la colonna dei dati luogo e quella della casetta intera.
l=[3:size(dati,3)];
y1 = squeeze(dati(:,3,:))';
y1=y1(l,:);

%seleziono gli ingressi: Testerna, Umidità, Solare Direttoe diffuso da dati
%luogo.
u2 = squeeze(dati(:,1,1))'; %T esterna
u3 = squeeze(dati(:,1,2))'; %Umidità
u4 = squeeze(dati(:,3,1))'; %solare normale diretto kW/m^2
u5 = squeeze(dati(:,4,1))'; %solare orizzontale diffuso kW/m^2

%compongo l'ingresso totale con quelli desiderati (mettere ;)
u1= [u2;u3];

%% campiono dati uscita e ingresso (solo per analisi effetti)
% tempo di campionamento (sapendo che sono già uno ogni 10 minuti!)
tempoCampionamento = 1;
camp = round(tempoCampionamento);
% sottocampiona secondo la dinamica ritenuta corretta
y1 = y1(:,1:camp:end);
u1 = u1(:,1:camp:end);

%% Troncatura dei dati
% per simulare periodi diversi dal set di campioni fornito.
% T è la riduzione in percentuale
 T=100;
    delta_T=(T*length(y1))/100;
    y1 = y1(:,1:delta_T);
    u1 = u1(:,1:delta_T);

%% Detrend: rimozione media
%tolgo la media ai dati con la funzione centra oppure è possibile usare la
%funzione detrend() con medesimi risultati.
[my, y]=centra(y1);
[mu, u]=centra(u1);

%% analisi spettrale dei segnali
%tale analisi è utile per determinare le componenti deterministiche del
%processo

%immetti nomi del segnale da analizzare con fft
sign=u2;

%interval è impostato come minuti per campione, ossia 10.
figure
[fftdata, freq, powerspect] = plotfft(sign, 10);

%% Imposta la parte di identificazione/validazione

frazione = 2/3; 
N = length(y);

% identificazione
y_id = y(:,1:floor(N*frazione)); 
u_id = u(:,1:floor(N*frazione));
% validazione
y_val = y(:,ceil(N*frazione):end);
u_val = u(:,ceil(N*frazione):end);

%aggiungo rumore agli ingressi con snr=10
%u_id = awgn(u_id1,10);
%essendo uguali posso usarne solo uno

%plottaggio I/O identificazione
figure
subplot(2,1,1)
plot(u_id');
grid on;
title('Ingressi per identificazione');
ylabel('°C (scalati)');
legend('Temp Ext','Umidità Ext')
axis([0 size(u_id,2) min(min(u_id))-2 max(max(u_id))+2])
subplot(2,1,2)
plot(y_id');
grid on;
title('Uscite per identificazione');
ylabel('°C (scalati)');
axis([0 size(y_id,2) min(min(y_id))-0.5 max(max(y_id))+0.5])
print(['-f',num2str(gcf)], '-depsc2', 'identif_yu')

%% Identificazione

%scelta metodo: pbsid_opt/pbsid/n4sid/pem/IOpredictor
tipo='pbsid';

switch tipo
    case 'pbsid_opt'
        disp('metodo: pbsid_opt')
        %PBSIDopt ricorsivo
        [A,B,C,D,K,lambda] = PBSID_opt_ARX_weighted_recursive(y_id,u_id,[ 8 8 ]);
    case 'pbsid'
        disp('metodo: pbsid')
        [A,B,C,D,K,lambda] = PBSID_opt_ARX_weighted(y_id,u_id,[ 8 8 ]);
    case 'n4sid'
        disp('metodo: n4sid')
        d = iddata(y',u');
        %m = n4sid(d,size(d.y,2),'N4Weight','CVA'); é instabile!!
        m = n4sid(d,size(d.y,2),'focus','simulation');
        A = m.A; B = m.B; C = m.C; D = m.D; K = m.K; lambda = m.NoiseVariance; 
        m.EstimationInfo
    case 'pem'
        disp('metodo: pem')
        d = iddata(y_id',u_id');
        m=pem(d);
        A = m.A; B = m.B; C = m.C; D = m.D; K = m.K; lambda = m.NoiseVariance; 
        m.EstimationInfo
    case 'IOpredictor'
        disp('metodo: IO predittore')
        [A,B,C,D,K,lambda] = IO_predictor(y_id,u_id,[ 2 2 ]);
    case 'subspace'
        [A,B,C,D,K,lambda] = subspace_id(y_id,u_id);
        %[A,B,C,D,K,lambda]=subspace_id(y_id,u_id,[],[],[],'proj_ls',[],1,'cva','ls');
end

%costruzione modello in spazio di stato del tipo 
    %x(t+1)=Ax(t)+Bu(t)
    %y(t)=Cx(t)+ Du(t)
%che non corrisponde al modello stocastico di partenza con K!!
m_ss=ss(A,B,C,0,1);

%raggiungibilità e osservabilità
if rank(ctrb(m_ss))==size(A,1)
    display('SISTEMA RAGGIUNGIBILE')
else display('SISTEMA NON RAGGIUNGIBILE')
end

if rank(obsv(m_ss))==size(A,1)
    display('SISTEMA OSSERVABILE')
else display('SISTEMA NON OSSERVABILE')
end

%verifico qualitativamente la stabilità con gli autovalori di A
autov=abs(eig(m_ss.A));
figure
plot(autov,'-.bo','MarkerFaceColor',[.49 1 .63],'MarkerSize',6);
grid on;
ylabel('abs');
title('autovalori');
axis([0 size(autov,1)+1 0 max(autov)+0.1])
print(['-f',num2str(gcf)], '-depsc2', 'autovalori')
save autov autov

%% stima condizioni iniziali: 

%2 metodi:
%1)dato che D=0 e conosco i dati delle uscite di identificazione, prendo le ultime uscite.
    %x0 = pinv(C)*y_id(:,end); % se conosco le misure!

%2)Compute prediction errors associated with model and data set. (migliora)
    d_val = iddata(y_val',u_val');
    d_id = iddata(y_id',u_id');
    dd = iddata(y',u');
    m_idss = idss(A,B,C,D,K);
    [e,x0] = pe(m_idss,d_val,'e');
    m_idss = idss(A,B,C,D,K,x0);
    
%% salvataggio dati e modelli 
%il nome deve essere modificato in funzione dell'esperimento trattato.

save d_val d_val
save d_id d_id
save dd dd
save m_idss m_idss
save m_ss m_ss

