%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   LTI_Systems_in_Matlab.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc
echo on
% Let us consider the time-continuous (TC) linear system (LS)
% characterized by the following matrices MA, MB and MC
% MA=[ 2     3     3;
%     -4    -5    -3;
%     -2     1    -1];
% MB=[  3; -2;  1];
% MC=[  2   1   0];
echo off
MA=[ 2     3     3;
    -4    -5    -3;
    -2     1    -1];
MB=[  3; -2;  1];
MC=[  2   1   0];
echo on
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A time-continuous linear system can be defined in Matlab using the
% command "ss" (state space). Using the command
SYS=ss(MA,MB,MC,0)
% the variable SYS is given the meaning of dynamic linear system
% characterized by the three matrices MA, MB and MC (the last matrix is MD=0).
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% It is possible to assign a name to all the input, state and output
%  variables using the command "set"
set(SYS,'InputName','Voltage')
set(SYS,'StateName',[   'Current '; 'Velocity';  'Position'])
set(SYS,'OutputName','Output')
SYS
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The transfer function of system G(s) can be obtained using the command "tf"
GS=tf(SYS)
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The poles of the system $G(s)$
pole(GS)
% are equal to the poles of the system SYS
pole(SYS)
% and equal to the eigenvalues of matrix MA
roots(poly(MA))
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The position of poles and zeros of system SYS can be plotted
% using the following command
figure(1); clf
pzmap(SYS)
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The step response is obtained using the following command
figure(1); clf
step(SYS); grid on
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The impulse response is obtained using the command
     figure(1); clf
impulse(SYS); grid on
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Bode diagrams are obtained using the command:
figure(1); clf
bode(SYS); grid on
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Nyquist diagrams are obtained using the command:
figure(1); clf
nyquist(SYS)
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(figure(1),'Visible','off')
clc
% The eigenvalues and the eigenvectors of matrix MA can be obtained
% using the following command
[V,D] = eig(MA)
% The first two eigenvectors (that is the first two columns of matrix V)
% are complex coniugate because they are the eigenvectors corresponding
% to the two complex conjugate eigenvalues. V is the transformation
% matrix which diagonalize the matrix $MA$
inv(V)*MA*V
% Deleting  the small computation errors one obtains:
approx(inv(V)*MA*V)
% The eigenvalues and the eigenvectors of matrix MA can be ontained
% using the following command
[V,D] = eig(MA)
% The first two eigenvectors (the first two columns of matrix V)
% are complex conjugate because they are the eigenvectors corresponding
% to the complex conjugate eigenvalues. V is the transformation
% matrix which diagonalize matrix $MA$
inv(V)*MA*V
% Deleting the small computation errors:
approx(inv(V)*MA*V)
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix MA can be transformed in the real Jordan canonical form using
% real transformation matrix: the complex conjugate eigenvalues
% are substituted by the real and the imaginary parts of one of
% the two complex coniugate eigenvectors
TM=[real(V(:,1)) imag(V(:,1)) V(:,3)]
inv(TM)*MA*TM
% Deleting the small computation errors:
approx(inv(TM)*MA*TM)
% Matrix MA can be transformed in the real Jordan canonical form, using
% a real transformation matrix: the complex conjugate eigenvectors
% are substituted by the real and imaginary parts of one of the two
% complex conjugate eigenvectors
TM=[real(V(:,1)) imag(V(:,1)) V(:,3)]
inv(TM)*MA*TM
% Deleting the small computation errors:
approx(inv(TM)*MA*TM)
pause; clc;             % Premi un tasto per continuare
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The same Jordan real canonical form can also be obtained
% using the following command
[CSYS,TM]=canon(SYS)
% The same Jordan real canonical form can also be obtained
% using the following command
[CSYS,TM]=canon(SYS)
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The free evolution of the system  is obtained using the command "initial":
% X0 initial state; X and Y state and output vectors; T time vector
figure(1); clf
X00=-2*V(:,3);
dx0=[0:0.1:1];
for dx=dx0
    X0=X00+[0;0;dx];
   [Y,T,X] = initial(SYS,X0);
   plot3(X(:,1),X(:,2),X(:,3)); hold on
   echo off
end
grid; xlabel('x'); ylabel('y'); zlabel('z')
echo on
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projection on the x-y plane
figure(1)
view(0,90)
title('Projection on the x-y plane')
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projection on the x-z plane
figure(1)
view(0,0)
title('Projection on the x-z plane')
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projection on the y-z plane
figure(1)
view(90,0)
title('Projection on the y-z plane')
pause; clc;             % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The same trajectory in the Jordan real canonical base
% transforms as follows:
figure(1); clf
for dx=dx0
    X0=X00+[0;0;dx];
   [Y,T,X] = initial(CSYS,-TM*X0);
   plot3(X(:,1),X(:,2),X(:,3)); hold on
   echo off
end
grid; xlabel('x'); ylabel('y'); zlabel('z')
echo on
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projection on the x-y plane
figure(1)
view(0,90)
title('Projection on the x-y plane')
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projection on the x-z plane
figure(1)
view(0,0)
title('Projection on the x-z plane')
pause;                  % Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Projection on the y-z plane
figure(1)
view(90,0)
title('Projection on the y-z plane')
echo off