%model order reduction (MOR) applicata alla diffusione 1D del calore:

% modello state space discreto:
% T(k+1) = A_ss T(k) + B_ss u(k)
% y(k)   = C_ss T(k)

%DEFINIZIONE DELLE MATRICI DEL MODELLO STATE-SPACE
n=20; %50;
Dx=0.05;
c=1;
tfin=0.1;
tin=0;
N=50; %N+1: numero di nodi temporali
dt=(tfin-tin)/N;
save dati c tfin tin N dt;

CB = 3;  % creo un modello con forzante a sinistra e temperatura fissa a destra:
[Ac1,B_ss,C_ss,D_ss] = crea_DLTI_StateSpace_calore_1D(n,CB,Dx);
A_ss = c * Ac1;
if CB==3
  n=n-1;
end

t = [tin:dt:tfin]';  % vettore degli istanti temporali di simulazione
u = f0(t);
x0 = g0(n,Dx);
[X,T_x0] = simula_DLTI_StateSpace_continuo_mor(A_ss,B_ss,C_ss,[0],u,x0,dt);
figure(1), clf, hold on, plot(T_x0,'r'), 
pause;

disp('pod:');
[Xuu,Xs,Xv]=svd(X);
%C=Xs*Xv';

% proiezione di Galerkin del modello state space continuo:
% NB: e' un modello continuo nel tempo ma discreto nello spazio!
% T(k+1)=Xu'*A_ss*Xu T(k)+Xu'*B_ss u(k)
% y(k)=C_ss*Xu* T(k);

%Xu è base ortogonale: bisogna scegliere dove troncarla
errore=[];

u = f0(t);
x0 = g0(n,Dx);

intervallo=[1:min([10 n])];
for m=intervallo; % m <= n+1 !!!
    % T approssimato da Xu*a
    Xu=Xuu(:,1:m);
    % applico la proiezione di Galerkin sul sistema continuo:
    A_pod=Xu'*A_ss*Xu;
    B_pod=Xu'*B_ss;
    C_pod=C_ss*Xu;

    %Risoluzione del sistema state space system avente matrici A_pod, B_pod, C_pod
    [hist_a_pod,a_x0_pod]=simula_DLTI_StateSpace_continuo_mor(A_pod,B_pod,C_pod,[0],u,Xu'*x0,dt);

    hist_T_pod=Xu*hist_a_pod;
    T_x0_pod=hist_T_pod(1,:);
    if 1
        figure(1), hold on, plot(T_x0_pod,'b--'), legend('sistema originario ','sistema ridotto ',2);title(['POD: m= ', num2str(m)]);
        pause;
        figure(1), hold on, plot(T_x0_pod,'g--'), 
    end
    errore = [errore, norm(T_x0_pod'-T_x0)^2/(N+1)];
end
figure(3), plot(errore), title('POD: andamento dell''errore al variare di m ');
pause;

disp(' ');

disp('balanced truncation');

% creo le matrici del sistema state-space discreto:
Ad = inv(eye(size(A_ss)) - dt*A_ss);
Bd = Ad*dt*B_ss;
Cd = C_ss;

T_x0_s=[];
iter=0;
figure(10), clf, hold on, plot(T_x0,'r'), 

for m = 1:10
    iter=iter+1;
    m
    [A_s,B_s,C_s,T,Tinv]=mr_svd(Ad,Bd,Cd,dt,m);
    %
    % Risoluzione del sistema state-space avente matrici A_s, B_s, C_s
    [hist1,dummy]=simula_DLTI_StateSpace_discreto_mor(A_s,B_s,C_s,[0],u,T*x0);
    hist_T_s = Tinv * hist1;
    T_x0_s = [T_x0_s (Cd*hist_T_s)'];
    figure(10), title(['m = ' num2str(m)]);
    figure(10), plot(real(T_x0_s(:,iter)),'g--'),
    title(['tronc bil (m=' num2str(m) ') : modello originario(rosso), modello ridotto(verde)']);
    pause;
    figure(10), plot(real(T_x0_s(:,iter)),'b--');    
end
