% ANALISI DEI SISTEMI
% ESERCITAZIONE 3 - 25 Ottobre 2002
%
% SVOLGIMENTO CON MATLAB
%
% PUNTO 1: Si determini una realizzazione del sistema in termini di variabili di stato
%
disp('numeratore e denominatore della f.d.t.')
num=[1 4 4]
den=[4 7 3]
pause
% la procedura TF2SS determina una realizzazione {a,b,c,d} che non coincide con 
% quella determinata mediante la procedura vista a lezione.
[a,b,c,d] = tf2ss(num,den);
% Tuttavia attraverso una trasformazione di similitudine che usa la matrice
%       P = [ 0  ...  0  a_n; 
%             0  ... a_n  0; 
%             ............. 
%            a_n ...  0   0]  
% è possibile passare alla rappresentazione {A,B,C,D} vista a lezione
n=length(den)-1;
a_n=den(1);
% definisco una matrice P con elementi tutti nulli
P = zeros(n,n);
% definisco la matrice P che esegue la trasformazione di similitudine
for i=1:n,
    P(i,n+1-i)=a_n;
end
% eseguo la trasformazione di similitudine
A=inv(P)*a*P
pause
B=inv(P)*b
pause
C=c*P
pause
D=d
pause

% PUNTO 2: Calcolare gli autovalori della matrice A e verificare che essi coincidano con le radici del 
% polinomio caratteristico  
%
disp('Autovalori della matrice A')
%calcolo gli autovalori relativi alla matrice A 
eig(A) 
disp('radici del polinomio caratteristico den')
% verifico che autovalori e radici del polinomio caratteristico coincidono
roots(den) 
pause

% PUNTO 3: Usare lo sviluppo di Sylvester per calcolare la matrice di transizione dello stato
%
% definisco l'asse dei tempi
t=0:0.01:10; 
% La matrice di transizione dello stato in questo caso e' una matrice  2 x 2
% che posso scrivere
%    exp(At) = [ f_11(t)  f_12(t); 
%                f_21(t)  f_22(t)];
% Le funzioni f_ij posso calcolarle numericamente con MATLAB
for i=1:1001
% il comando expm esegue l'esponenziale di una matrice  
  Q=expm(A*t(i));  
  f_11(i)=Q(1,1);
  f_12(i)=Q(1,2);
  f_21(i)=Q(2,1);
  f_22(i)=Q(2,2);
end
% Per verificare i risultati calcolo anche la f_11 determinata analiticamente
% che chiamo f_11_an
f_11_an=-3*exp(-t)+4*exp(-3*t/4);
plot(t,f_11,t,f_11_an)
grid
xlabel('t')
title('Confronto fra l''elemento (1,1) di e^At determinato da MATLAB e analiticamente')
pause
% La stessa procedura di verifica puo' essere applicata alle altre componenti di e^At

% PUNTO 4: Calcolo condizioni iniziali
%
Q=[C; C*A];
disp('Stato iniziale')
x0=inv(Q)*[2;1]
pause

% PUNTO 5: Calcolare in funzione delle matrici A,B,C,D la risposta impulsiva
%
% il comando ss crea il sistema lineare a tempo continuo in termini di variabili di stato
sys=ss(A,B,C,D);
% il comando impulse genera la risposta impulsiva del modello creato con il comando ss 
% t è il vettore riga dei tempi
w=impulse(sys,t);
% verifica risultati ottenuti analiticamente
w_an=-exp(-t)+1.56*exp(-3/4*t);
plot(t,w,t,w_an)
grid
xlabel('t')
title('Risposta impulsiva (calcolata da MATLAB e analiticamente)')
pause

% PUNTO 6: Si determini la trasformazione di similitudine che permette di passare ad una realizzazione 
% in cui la matrice A' è in forma diagonale, si calcoli la corrispondente realizzazione 
%
disp('Matrice di Vandermonde')
% definisco la matrice di Vandermonde
V=[1 1; -1 -3/4] 
disp('Rappresentazione diagonale')
A1=inv(V)*A*V
pause
B1=inv(V)*B
pause
C1=C*V
pause
D1=D
sys1=ss(A1,B1,C1,D1);
pause

% PUNTO 8: Usare la formula di Lagrange per calcolare l'evoluzione forzata dell'uscita 
% in corrispondenza dell'applicazione dell'ingresso dato
%
% definisco l'asse dei tempi
t1=[0:0.01:1.99];
t2=[2:0.01:8];
t=[t1 t2];
% definisco l'ingresso
u1=0.5*t1;
u2=zeros(size(t2));
u=[u1 u2];
% il comando lsim genera l'evoluzione forzata dello stato xf e dell'uscita yf
% conseguente all'applicazione di un ingresso u  t è il vettore riga dei tempi
[yf,t,xf]=lsim(sys1,u,t);
%
% La prima colonna di xf contiene i valori della prima componente: x1
xf1=xf(:,1);
% La seconda colonna di xf contiene i valori della seconda componente: x2
xf2=xf(:,2);
% verifica dei risultati ottenuti analiticamente
% verifica del vettore di stato
xf1_an= [-1/2*t1-1/2*(exp(-t1)-1)  -4.19*exp(-t2)];
xf2_an= [2/3*t1+8/9*(exp(-3*t1/4) -1)  2.88*exp(-3*t2/4)];
plot(t,xf1,t,xf1_an)
grid
xlabel('t')
title('Evoluzione forzata primo stato x1 (calcolata da MATLAB e analiticamente)')
pause
plot(t,xf2,t,xf2_an)
grid
xlabel('t')
title('Evoluzione forzata secondo stato x2 (calcolata da MATLAB e analiticamente)')
pause
% verifica dell'uscita 
yf_an= [16*t1/24-(1/18)*(16+9*exp(-t1)-25*exp(-3*t1/4))  -4.19*exp(-t2)+4.5*exp(-3*t2/4)];
plot([t1 t2],yf,[t1 t2],yf_an)
grid
xlabel('t')
title('Risposta forzata yf(t) (calcolata da MATLAB e analiticamente)')