clc
clear all
close all

%Esercizio 1
%====================================================================================================================================
%====================================================================================================================================
%====================================================================================================================================
%1.1
disp('=================================================================================================');
disp('L.I.N.F.A. Laboratorio Informatico-Numerico di Controlli Automatici anno accademico 2002/2003');
disp('Prof. Paolo Rocco');
disp('Esercitazione n°1: complementi di teoria dei sistemi');
disp('=================================================================================================')

disp('Esercizio n°1');

input('(premi invio per continuare....)');

input('definisco A');
A = [4 3 0 -3; 0 3 1 -1; 0 2 3 -2; 0 2 1 0];

disp(A);
%=============================================== calcolo l'esponenziale secondo la definizione ======================================
disp('Esponenziale matrice calcolata in due modi differenti: tramite funzione built-in di MATLAB e tramite definizione:');
input('(premi invio per continuare....)');

fi = expm(A);
myExp = eye(4);                                 % valore iniziale per iterazione
clc
for i= 1:1:20,
        myExp = myExp+(A^(i))/factorial(i);
        disp('esponenziale built-in di MATLAB (fi = expm(A))');
        disp(fi);
        fprintf('esponenziale calcolato all iterazione %d\n',i)  ;
        disp(myExp);
        disp('differenza:');
        disp(fi-myExp);
        input('premi un tasto...');
        clc
end

clc
input('Autovalori ed autovettori di A');

[V,D] = eig(A);

disp('autovettori');
disp(V);

disp('autovalori');
lambda = diag(D);
disp(lambda);

input('premi un tasto...');
clc

input('(sono davvero autovalori / autovettori di A??  ovvero: (A-lamba*I)*V = 0_)??');
for k =    1: 1 :4
    clc 
    disp(' A*V=lamba*V ?');
    fprintf('A*V%d\n',k);
    disp(A*V(:,k));

    fprintf('lamba%d*V%d\n',k,k);
    disp(lambda(k)*V(:,k));
    input('premi un tasto...');
end
input('ok.');
input('(premi invio per continuare....)');

clc
input('T che mette in forma diagonale A..');
T = inv(V)
input('(premi invio per continuare....)');
clc;
input('1.2: exp(A) = T^-1*diag(e^lamba_i)T)?');

disp('expm(A):');
disp(fi);
disp('T^-1*diag(e^lamba_i)T):');
disp(inv(T)*expm(D)*T);
disp('ok.');
%V*A*V^(-1)

input('premi Invio...');

% % Esercizio 2
% %====================================================================================================================================
% %====================================================================================================================================
% %====================================================================================================================================

clc; clear;close all;
input('esercizio 2...');
clc;

% costruisco le matrici ABCD del sistema

A = [-4 -3 0 +3 ;0 -3 -1 +1;0 -2 -3 +2;0 -2 -1  0];
B = [ 0 ;0 ;0 ;1];
C = [1 1  0 0];
D = [0];

x0 = [1; 0; 1; 2];

% 
%creo il sistema
clc;
input('definisco il sistema');
sist    = ss(A,B,C,D)
input('premi invio...');
clc;
input('2.1: verifica asintotica stabilità:');
eig(A)
eig(sist)
disp('autovalori con parte reale negativa: ok. ');                    

dt      = 0.01;
T       = 6;
t       = 0.0:dt:T;

input('premi invio');
clc
input('adesso calcolo il movimento libero in 2 modi diversi: ');
disp('metodo built-in di MATLAB ([y,t,x] = initial(sist,x0,t))...');
disp('...e tramite la formula di Lagrange ');

[y,t,x] = initial(sist,x0,t);                       % calcolo movimento libero tramite funzione MATLAB    

input('premi invio');

 for k = 0 :1:T/dt,
     x1(:,k+1) = expm(A*t(k+1))*x0;                 % calcolo movimento libero tramite formula di Lagrange
 end
 
 x1 = x1';

figure;                                             %confronto grafico    
subplot(3,2,1);
plot(t,x(:,1));             title('built-in');   ylabel('x1');

subplot(3,2,3);                                     
plot(t,x(:,2));                                  ylabel('x2');

subplot(3,2,5);
plot(t,x(:,3));                                  ylabel('x3');  xlabel('t [s]');

subplot(3,2,2);         
plot(t,x1(:,1));                                 title('Formula di Lagrange');

subplot(3,2,4);
plot(t,x1(:,2));

subplot(3,2,6);     
plot(t,x1(:,3));                            xlabel('t [s]');

input('Il moto libero decade a zero coerentemente con le proprietà di asintotica stabilità del sistema');

% Esercizio 3
%====================================================================================================================================
%====================================================================================================================================
%====================================================================================================================================

clc; close all;
clc;
input('esercizio 3...');
input('Definizione del sistema...');

A = [1 2 -2; 0 -1 0; 0 0 -1]
B = [0; 0; 0]
C = [0 0 0]    
D = [0]

x0 = [1; 2; 3];                             % stato iniziale
input('premi Invio...')
clc;

input('3.1: definizione del sistema');
sist = ss(A,B,C,D)

input('3.1: verifica asintotica stabilità:');
clc;
eig(A)
eig(sist)
input('1° autoval con parte reale positiva: instabile!!!!');
% p = roots(poly(A))                            % metodo alternativo (condizione necessaria ma non sufficiente)
clc;

input('calcolo del movimento libero a partire da x0...');
[y,t,x] = initial(sist,x0,t);                                   


figure;
subplot(3,1,1);
plot(t,x(:,1));

subplot(3,1,2);
plot(t,x(:,2));

subplot(3,1,3);
plot(t,x(:,3));

input ('Questo moto libero è limitato...');

clc
input('calcolo del movimento libero a partire da x1...');
x1 = [1.01;2.0; 3.0];

[y,t,x] = initial(sist,x1,5);                                   

input('perturbando la prima variabile di stato il moto libero diverge...');
figure;
subplot(3,1,1);
plot(t,x(:,1));

subplot(3,1,2);
plot(t,x(:,2));

subplot(3,1,3);
plot(t,x(:,3));

input('premi Invio...');

% Esercizio 4
%====================================================================================================================================
%====================================================================================================================================
%====================================================================================================================================

clc; clear;close all;
input('esercizio 4...');
clc;

A=[0,0,1;0,0,1;-1,-1,-1];
B=[0;0;1];
C=[0,0,1];

sist=ss(A,B,C,0);
disp('matrice di raggiungibilità:');
Kr=ctrb(sist)
disp('det(Kr)');
det(Kr)

input('det(A) =0 (=> rango non massimo): il sistema non è completamente raggiungibile');

Ko=obsv(sist)
disp('det(Ko)');
det(Ko)

input('det(A) =0 (=> rango non massimo): il sistema non è completamente osservabile');

clc;
disp('funzione di trasferimento del sistema');
tf(sist)


input('La f.d.t. è quella di un sistema del secondo ordine a parte l errore numerico)');

clc
input('Prendiamo come nuove variabili somma e differenza delle tensioni sui condensatori e la corrente sull induttore');

T=[1,-1,0;1,1,0;0,0,1]
sist2           = ss2ss(sist,T);
[A2,B2,C2,D2]   = ssdata(sist2);
A2,B2,C2

%La prima variabile di stato non contribuisce al legame ingresso-uscita

pause %Premere un tasto per continuare...
tf(sist2)
input('La funzione di trasferimento è la stessa!! ');
input(' infatti il sistema fisico descritto da (A,B,C,D) è lo stesso descritto da (A2,B2,C2,D2)');












