%% Modélisation, analyse et commande d'un système modélisé dans l'espace d'état
%
%% Modélisation dans l'espace d'état et relation avec la fonction de transfert.
%
% On considère un système:

clear all;
close all;


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

%
%
% On définit un système dynamique dans l'espace d'état:

G=ss(A,B,C,D) % ss pour State Space

% On peut calculer la fonction de transfert correspondante:

G_ft=tf(G)
[Num,Den]=ss2tf(A,B,C,D)
G_FT2=tf(Num,Den)

%% Analyse des systèmes modélisés dans l'espace  d'état
% Analyse des valeurs propres (permettant par exemple d'étudier la stabilité asymptotique)

eig(A)

% Analyse des zéros du système

z=zero(G)% zéros d'un système SISO

% Analyse de la commandabilité

Ccom= ctrb(A,B)  %Calcul de la matrice de commandabilité.
rank(Ccom)  %Calcul du rang de la matrice de commandabilité.

% Analyse de l'observabilité

Cobs= obsv(A,C) % Calcul de la matrice d'observabilité.
rank(Cobs) % Calcul du rang de la matrice d'observablité.

%
%% Forme modale.
% Calculons une forme modale. La matrice dynamique est alors diagonale ou sous forme de bloc de Jordan

[V, J] = eig(A); % On calcule les valeurs propres de A ainsi que la matrice des vecteurs propres.
T = V; % On définit la matrice de changement de base.
A_mod = inv(T)*A*T  % On effectue le changement de base x=Tz
B_mod = inv(T)*B
C_mod = C*T

G_mod = ss(A_mod, B_mod, C_mod, D)

% On peut aussi utiliser la commande canon ou modalreal

[G_mod2,T2]=canon(G,'modal')
[G_mod3,blks,TL,TR]=modalreal(G)

%% Forme compagne de commande et/ou d'observation
%% Forme compagne de commande
% Si le système modélisé dans l'espace d'état est commandable nous pouvons calculer la forme compagne de commande.

% Vérification de la commandabilité
Ccom = ctrb(A, B);
disp('Rang de la matrice de commandabilité :');
rank_Com = rank(Ccom)
%[Gcanon,Tcanon] = canon(ss(A,B,C,D), 'companion');
sys_co=compreal(G,"c")

%
% Si le système modélisé dans l'espace d'état est observable, nous pouvons calculer la forme compagne d'observation

% Vérification de l'observabilité
Cobs=ctrb(A,B);
disp('Rang de la matrice de commandabilité :');
rank_Cobs = rank(Cobs)
sys_ob = compreal(G,"o")

%
%%
%% Utilisation des tracés initial, impulse step ... lsim

%% Réponses temporelles à différentes entrées

t = 0:0.01:10;
u = sin(t);        % entrée sinusoïdale
x0 = [1; -1];      % Conditions initiales non nulles

figure;
subplot(4,1,1);
step(G,t); 
title('Réponse indicielle');

subplot(4,1,2);
impulse(G,t);
title('Réponse impulsionnelle');

subplot(4,1,3);
lsim(G,u,t);
title('Réponse à une entrée sinusoïdale');

subplot(4,1,4);
initial(G,x0,t);  % Réponse aux conditions initiales
title('Réponse due aux conditions initiales x0 ');

%



%% Synthèse d'une loi de commande dans l'espace d'état : calcul d'un retour d'état.
% On désire déterminer une loi de commande permettant d'obternir des valeurs propres en boucle fermée égale à -2 et -3.

poles_desires = [-2 -3] % temps de réponse de 3/2 :
K = place(A, B, poles_desires)    % Calcul du gain d'état

% Caclul du système en boucle fermée.

A_bf = A - B*K
sys_bf = ss(A_bf, B, C, D)
eig(A_bf)    % Doit afficher [-2 -3]

% Tracé de la réponse indicielle en boucle fermée

figure;
step(sys_bf);
title('Réponse indicielle en boucle fermée avec retour d''état');

% La commande en boucle fermée ne permet pas à priori de régler des problèmes de précision. Si par exemple, nous souhaitons obtenir une erreur de position nulle, nous pouvons ajouter un gain de précompensation N1:

% Calcul du gain de précompensation
N1 = 1/dcgain(sys_bf) 
% Creation du système en boucle fermée
sys_bf2=N1*sys_bf;
dcgain(sys_bf2) % Le système a un gain statique égal à 1.
%tracé de la réponse indicielle pour vérification
figure;
step(sys_bf2);
title('Réponse indicielle en boucle fermée avec retour d''état et précompensation');

%
%
%
%
%
%
%
%
%
%
%
%
%