Modelisation et Analyse des Systemes Lineaires et Representation d'Etat - S5

Annee : 2022-2023 (Semestre 5)
Credits : 3 ECTS
Type : Automatique et Systemes
Enseignant : Subias

---

PART A : PRESENTATION GENERALE

Objectifs du cours

Ce cours constitue une formation avancee en theorie des systemes lineaires et en representation d'etat, piliers de l'automatique moderne. L'objectif est de maitriser les techniques de modelisation mathematique, l'analyse de comportement des systemes dynamiques, et les methodes de commande par retour d'etat.

Competences visees

• Maitriser les techniques de modelisation des systemes dynamiques lineaires
• Analyser la stabilite et les performances des systemes en boucle ouverte et fermee
• Concevoir des lois de commande par retour d'etat avec placement de poles
• Dimensionner des observateurs d'etat pour estimer les variables non mesurables
• Utiliser MATLAB/Simulink pour l'analyse et la simulation de systemes automatiques
• Appliquer les criteres de commandabilite et observabilite
• Interpreter les reponses frequentielles (Bode, Nyquist, marges de stabilite)

Organisation

• Volume horaire : Cours magistraux, TD et TP MATLAB/Simulink
• Evaluation : Examen ecrit, TPs notes, Projet sur robot auto-equilibre
• Outils : MATLAB Control System Toolbox, Simulink
• Prerequis : Transformee de Laplace, equations differentielles, algebre lineaire

---

PART B : EXPERIENCE, CONTEXTE ET FONCTION

Contenu pedagogique

Le cours s'articule autour de trois axes principaux : la modelisation des systemes, l'analyse de stabilite et performances, et la synthese de lois de commande.

1. Modelisation des systemes lineaires

Representation par fonction de transfert

La fonction de transfert relie l'entree U(s) a la sortie Y(s) en transformee de Laplace :

H(s) = Y(s) / U(s) = N(s) / D(s)

Exemple de systeme du second ordre (moteur DC) :

H(s) = Km / (Tm*s + 1)

Avec Km = gain statique, Tm = constante de temps.

Representation d'etat

Forme generale pour un systeme SISO :

dx/dt = A*x + B*u  (equation d'etat)
y = C*x + D*u      (equation de sortie)

Matrices d'etat :

• A (n x n) : matrice dynamique
• B (n x 1) : matrice de commande
• C (1 x n) : matrice d'observation
• D : action directe (souvent 0)

Exemple concret d'un systeme d'ordre 2 :

A = [0 1; -10 -6.316]
B = [0; 1]
C = [k/0.1001 0]
D = 0

2. Stabilite des systemes

Critere de Routh-Hurwitz

Pour le polynome caracteristique D(s) = an*s^n + ... + a1*s + a0, on construit le tableau de Routh :

s^n	an	an-2	an-4
s^n-1	an-1	an-3	an-5
s^n-2	b1	b2	...
...	...	...	...

Avec b1 = (an-1*an-2 - an*an-3) / an-1

Condition de stabilite : tous les elements de la premiere colonne doivent etre strictement positifs.

Implementation MATLAB disponible avec la fonction routh.m :

function RA=routh(poli,epsilon)
% Calcule le tableau de Routh pour un polynome
% Gere les cas particuliers (premiere colonne nulle, ligne de zeros)

Analyse par les poles

Un systeme est stable si tous les poles (racines de D(s)) ont une partie reelle negative :

Re(pi) < 0 pour tout i

Classification :

• Stable asymptotiquement : tous les poles a partie reelle strictement negative
• Marginalement stable : poles sur l'axe imaginaire (non repetes)
• Instable : au moins un pole a partie reelle positive

3. Analyse frequentielle

Diagrammes de Bode

Representation du gain (en dB) et de la phase (en degres) en fonction de la frequence :

Gain(dB) = 20*log10(|H(jw)|)
Phase(deg) = arg(H(jw)) * 180/pi

Exemple de TP1 avec fonction de transfert H(s) = 0.717 / (0.0033*s + 1) :

H = tf([0.717], [0.0033 1]);
bode(H)

Mesures experimentales effectuees :

• Frequences : 0.05 a 10 Hz
• Gains mesures : de 0.669 a 0.061 (de -3.5dB a -24dB)
• Dephasages : de -2.2 deg a -90 deg

Marges de stabilite

Marge de gain (Gm) : gain additionnel avant instabilite a la frequence de coupure de phase -180 deg

Marge de phase (Pm) : phase additionnelle avant instabilite a la frequence de coupure de gain 0dB

Commande MATLAB :

[Gm, Pm, Wcg, Wcp] = margin(sys);

Criteres de robustesse :

• Pm > 45 deg (bonne robustesse)
• Gm > 6dB (bonne robustesse)

4. Commande par retour d'etat

Placement de poles

Pour un systeme commandable, on peut imposer les poles en boucle fermee avec un retour d'etat u = -K*x + r.

La matrice de gain K est calculee pour placer les poles desires :

K = place(A, B, poles_desired)

Exemple du TP3 avec poles complexes conjugues :

K = place(A, B, [-2+2i, -2-2i])

Les poles choisis determinent les performances :

• Partie reelle : rapidite (plus negatif = plus rapide)
• Partie imaginaire : oscillations (plus grand = plus d'oscillations)

Commandabilite

Un systeme est commandable si on peut amener l'etat de n'importe quelle condition initiale a n'importe quel etat desire en temps fini.

Critere de Kalman : rang de la matrice de commandabilite = n

Matrice de commandabilite :

Mc = [B | A*B | A^2*B | ... | A^(n-1)*B]

Test MATLAB :

Mc = ctrb(A, B);
rank_Mc = rank(Mc);
% Systeme commandable si rank_Mc == n

5. Observateurs d'etat

Principe

Lorsque tous les etats ne sont pas mesurables, on utilise un observateur pour reconstruire l'etat complet a partir de la sortie mesuree y.

Equation de l'observateur de Luenberger :

dx_hat/dt = A*x_hat + B*u + L*(y - y_hat)
y_hat = C*x_hat

Avec L : gain de l'observateur.

Observabilite

Un systeme est observable si on peut determiner l'etat initial x(0) a partir de la connaissance de u(t) et y(t) sur un intervalle fini.

Critere de Kalman : rang de la matrice d'observabilite = n

Matrice d'observabilite :

Mo = [C; C*A; C*A^2; ...; C*A^(n-1)]

Test MATLAB :

Mo = obsv(A, C);
rank_Mo = rank(Mo);
% Systeme observable si rank_Mo == n

Placement des poles de l'observateur

Les poles de l'observateur sont choisis plus rapides (2 a 5 fois) que ceux de la commande pour garantir une convergence rapide de l'estimation.

Exemples du TD :

L1 = [-0.365, -0.78]   % Poles lents
L2 = [0.59, 0.32]      % Poles rapides

6. Analyse modale et formes canoniques

Decomposition modale

Si A est diagonalisable, on peut ecrire :

A = P * Lambda * P^(-1)

Avec Lambda matrice diagonale des valeurs propres et P matrice des vecteurs propres.

Changement de variable z = P^(-1)*x :

dz/dt = Lambda*z + P^(-1)*B*u
y = C*P*z

Chaque mode evolue independamment selon zi(t) = exp(lambda_i*t).

Formes canoniques

Forme commandable : matrices A et B ont une structure particuliere facilitant le placement de poles.

Forme observable : matrices A et C ont une structure facilitant la conception d'observateurs.

---

C. Aspects techniques et pratiques

1. Travaux Pratiques MATLAB

TP1 : Identification experimentale

Manipulation d'un systeme reel (moteur ou servomoteur) pour mesurer sa reponse frequentielle.

Code de traitement des donnees :

% Donnees experimentales
Vg = [-3.93, -2.39, -1.17, -0.44, 0, 0.37, ...];
O = [-292*(2*pi/60)*9, -179*(2*pi/60)*9, ...]; % Vitesse en rad/s

% Trace caracteristique statique
plot(O, Vg);
xlabel('Vitesse angulaire (rad/s)');
ylabel('Tension (V)');

% Analyse frequentielle
G = [0.669, 0.659, 0.609, 0.528, ...];
G_dB = 20*log10(G);
freq = [0.05, 0.1, 0.25, 0.5, 1, 2.5, 5, 7.5, 10];

semilogx(freq, G_dB);
xlabel('Frequence (Hz)');
ylabel('Gain (dB)');

Modele identifie : H(s) = 0.717 / (0.0033*s + 1)

• Gain statique : K = 0.717
• Constante de temps : T = 3.3 ms

TP2 : Analyse de stabilite en boucle fermee

Systeme avec capteur, moteur et retour :

Ks = 1.55;   % Gain capteur
Km = 47.8;   % Gain moteur
Kg = 0.015;  % Gain reducteur
Tm = 0.3;    % Constante de temps moteur
K = 0.887;   % Gain global

% Fonction de transfert en boucle fermee
num = [Km*K/9];
den = [Tm, 1, Km*K*Ks/9];
sys = tf(num, den);

% Analyse de stabilite
[Gm, Pm, Wcg, Wcp] = margin(sys);
disp(['Marge de gain : ', num2str(Gm), ' (', num2str(20*log10(Gm)), ' dB)']);
disp(['Marge de phase : ', num2str(Pm), ' deg']);

% Calcul des poles
P = pole(sys);
disp('Poles du systeme :');
disp(P);

Resultats typiques :

• Pm ~ 60 deg → Systeme bien amorti
• Poles a partie reelle negative → Systeme stable

TP3 : Commande par retour d'etat

Systeme d'ordre 2 avec placement de poles :

k = 5;
A = [0 1; -10 -6.316];
B = [0; 1];
C = [k/0.1001 0];

% Creation du modele d'etat
sys = ss(A, B, C, 0);

% Placement de poles en -2 +/- 2j
poles_desired = [-2+2i, -2-2i];
K = place(A, B, poles_desired);

disp('Gain de retour d''etat K :');
disp(K);

% Systeme en boucle fermee
A_cl = A - B*K;
sys_cl = ss(A_cl, B, C, 0);

% Reponse indicielle
step(sys_cl);
title('Reponse en boucle fermee avec retour d''etat');

2. Projet : Robot auto-equilibre (NXT Way)

Modelisation du robot Lego Mindstorms EV3

Parametres physiques du fichier NXP_robot_parameters.m :

Parametre	Valeur	Unite	Description
m	0.03	kg	Masse d'une roue
R	0.042	m	Rayon des roues
Jw	m*R^2/2	kg.m^2	Inertie d'une roue
M	0.67	kg	Masse du corps
H	0.152	m	Hauteur du corps
L	H/2	m	Distance centre de masse / axe roues
Jpsi	M*L^2/3	kg.m^2	Inertie en tangage
g	9.81	m/s^2	Acceleration de la gravite

Moteur DC :

• Jm = 1e-5 kg.m^2 (inertie moteur)
• Rm = 6.83 Ohm (resistance)
• Kb = 0.468 V.s/rad (constante de force contre-electromotrice)
• Kt = 0.3047 N.m/A (constante de couple)

Simulation Simulink

Le projet utilise NXTwaySim.slx pour simuler le comportement du robot avec :

• Modele non-lineaire du pendule inverse
• Commande par retour d'etat
• Observateur d'etat pour estimer l'angle et la vitesse angulaire
• Animation 3D du robot avec PlayAnimation.m

Objectif : stabiliser le robot en position verticale (theta = 0 deg) malgre les perturbations.

3. Outils MATLAB essentiels

Fonctions de base

% Creation de modele
sys = tf(num, den);              % Fonction de transfert
sys = ss(A, B, C, D);            % Representation d'etat

% Conversion
[A, B, C, D] = tf2ss(num, den);  % TF -> Etat
[num, den] = ss2tf(A, B, C, D);  % Etat -> TF

% Analyse
pole(sys);                       % Calcul des poles
zero(sys);                       % Calcul des zeros
eig(A);                          % Valeurs propres de A
rank(ctrb(A, B));               % Test de commandabilite
rank(obsv(A, C));               % Test d'observabilite

% Reponses temporelles
step(sys);                       % Reponse indicielle
impulse(sys);                    % Reponse impulsionnelle
lsim(sys, u, t);                % Reponse a une entree quelconque

% Analyse frequentielle
bode(sys);                       % Diagrammes de Bode
nyquist(sys);                    % Diagramme de Nyquist
nichols(sys);                    % Diagramme de Black-Nichols
margin(sys);                     % Marges de stabilite

% Synthese
K = place(A, B, poles);          % Placement de poles
L = place(A', C', poles)';       % Observateur (dualite)

4. Methodes de conception

Choix des poles en boucle fermee

Pour un systeme d'ordre 2, forme canonique :

H(s) = wn^2 / (s^2 + 2*zeta*wn*s + wn^2)

Parametres de performance :

• wn : pulsation naturelle (rapidite)
• zeta : coefficient d'amortissement (depassement)

Relation avec les poles :

p1,2 = -zeta*wn +/- j*wn*sqrt(1-zeta^2)

Recommandations :

• zeta = 0.7 : bon compromis (depassement ~ 5%)
• Temps de reponse a 5% : tr ~ 3 / (zeta*wn)

Separation des dynamiques

Principe de separation : on peut concevoir independamment le retour d'etat K et l'observateur L, puis les combiner.

Regle pratique : placer les poles de l'observateur 2 a 5 fois plus rapides que ceux de la commande.

Exemple :

• Poles commande : -2 +/- 2j
• Poles observateur : -10, -12 (5 fois plus rapides)

---

D. Analyse et perspectives

1. Applications industrielles

Systemes de regulation

• Procedes industriels : regulation de temperature, pression, debit
• Robotique mobile : stabilisation de robots a roues, drones
• Aerospatial : pilotage automatique d'avions, fusees
• Automobile : regulation de vitesse (cruise control), ESP, suspension active

Exemples concrets

Segway / robots auto-equilibres : application directe du projet NXT Way avec pendule inverse et commande par retour d'etat.

Quadrirotors : modele d'etat avec 12 etats (position, vitesse, angles, vitesses angulaires), commande hierarchique avec boucles imbriquees.

2. Liens avec d'autres cours

Cours	Semestre	Lien avec Modelisation Systemes Lineaires
Circuits et Filtres Analogiques	S5	Analyse frequentielle, Bode, systemes du 1er et 2e ordre
Fondements Electronique Numerique	S5	Echantillonnage pour commande numerique
Systemes Boucles	S5	Asservissements, correcteurs PID
Commande Numerique	S8	Discretisation, commande echantillonnee
Temps Reel	S8	Implementation de lois de commande
Processus Stochastiques	S8	Filtre de Kalman (extension de l'observateur)

3. Extensions et sujets avances

Commande optimale

Minimisation d'un critere de performance (LQR - Linear Quadratic Regulator) :

J = integrale[x'Qx + u'Ru] dt

Conduit a un gain optimal K calcule via equation de Riccati.

Systemes non-lineaires

Pour les systemes non-lineaires, linearisation autour d'un point d'equilibre :

dx/dt = f(x, u) ~ f(x0, u0) + df/dx|(x0,u0) * (x-x0) + df/du|(x0,u0) * (u-u0)

Obtention d'un modele lineaire tangent exploitable avec les techniques du cours.

Systemes multivariables (MIMO)

Extension a plusieurs entrees et sorties :

• Matrices B (n x m), C (p x n), D (p x m)
• Analyse de couplage entre canaux
• Commande decouplante

4. Competences professionnelles developpees

Ce cours a permis de developper des competences cles en automatique moderne, essentielles pour la carriere d'ingenieur en systemes embarques, robotique, ou automatisation industrielle.

Modelisation :

• Mise en equations d'un systeme physique a partir des lois fondamentales
• Identification experimentale de parametres par analyse frequentielle
• Validation de modeles par comparaison simulation/reel
• Choix pertinent des variables d'etat

Analyse :

• Determination de stabilite par methodes algebriques (Routh) et graphiques (lieu des poles)
• Evaluation de performances (rapidite, depassement, erreur statique)
• Robustesse aux variations de parametres et perturbations
• Interpretation des diagrammes de Bode et marges de stabilite

Synthese :

• Conception de lois de commande par placement de poles selon cahier des charges
• Dimensionnement d'observateurs d'etat pour estimation des variables non mesurables
• Optimisation de criteres de performance (temps de reponse, depassement)
• Application du principe de separation

Outils informatiques :

• Maitrise de MATLAB/Simulink pour l'automatique (Control System Toolbox)
• Programmation de fonctions d'analyse (Routh, identification, placement de poles)
• Simulation de systemes complexes avec Simulink
• Visualisation et interpretation de resultats

5. Methodologie de resolution de problemes

Demarche type pour un probleme d'automatique

1. Modelisation : etablir les equations differentielles ou fonction de transfert
2. Mise sous forme d'etat : definir les variables d'etat pertinentes
3. Analyse du systeme en boucle ouverte :
   • Stabilite (Routh, poles)
   • Commandabilite et observabilite
   • Performances (reponse indicielle)
4. Conception de la commande :
   • Choix des poles desires selon cahier des charges
   • Calcul du gain K par placement de poles
   • Verification des performances
5. Conception de l'observateur (si necessaire) :
   • Choix des poles de l'observateur (plus rapides)
   • Calcul du gain L
6. Simulation et validation :
   • Test avec Simulink
   • Analyse de robustesse
   • Ajustements si necessaire

6. Pieges a eviter

L'experience des TDs et du projet a permis d'identifier plusieurs erreurs courantes en automatique.

Pieges theoriques :

• Systeme non commandable/observable : toujours verifier les criteres de Kalman avant de tenter un placement de poles ou de concevoir un observateur
• Poles de l'observateur trop rapides : risque d'amplification du bruit de mesure et d'instabilite numerique
• Confusion entre poles et zeros : les poles determinent la stabilite, les zeros affectent la forme de la reponse
• Oubli des conditions initiales : la reponse complete = reponse libre + reponse forcee

Pieges pratiques :

• Unites incoherentes : attention aux conversions rad/s ↔ tr/min, degres ↔ radians
• Saturation des actionneurs : prendre en compte les limites physiques (commande limitee en amplitude et vitesse)
• Modele lineaire hors zone de validite : linearisation valable uniquement autour du point d'equilibre
• Echantillonnage insuffisant : respect du critere de Shannon pour l'implementation numerique

Pieges MATLAB :

• Confusion entre systemes SISO et MIMO : attention aux dimensions des matrices
• Mauvais choix de solveur dans Simulink : peut entrainer des erreurs numeriques
• Oubli de la normalisation : certaines fonctions MATLAB attendent des parametres normalises

7. Ressources et approfondissements

Documentation MATLAB

• Control System Toolbox User's Guide : documentation complete des fonctions
• Simulink Control Design : conception et analyse de systemes de commande
• Exemples integres : help control, demo control dans MATLAB
• MATLAB Central : communaute avec exemples et scripts partages

Ouvrages de reference

• "Modern Control Engineering" - Ogata : reference classique, approche pedagogique
• "Linear System Theory and Design" - Chen : approche mathematique rigoureuse
• "Feedback Control of Dynamic Systems" - Franklin, Powell, Emami-Naeini : excellent pour les applications
• "Automatique : Systemes lineaires" - Duc, Theron : en francais, adapte au cursus francais

Logiciels alternatifs

• Python : bibliotheques control, scipy.signal pour calculs scientifiques
• Scilab : alternative open-source a MATLAB avec Xcos (equivalent Simulink)
• Octave : compatible MATLAB, gratuit
• LabVIEW : pour applications industrielles et acquisition de donnees

Mon opinion

Ce cours represente un pilier fondamental de l'automatique moderne et s'avere indispensable pour tout ingenieur travaillant dans les domaines des systemes embarques, de la robotique, ou de l'automatisation industrielle.

Pourquoi ce cours est essentiel :

1. Approche moderne : la representation d'etat est la base de l'automatique contemporaine (vs approche classique par fonction de transfert)
2. Systemes multivariables : traitement naturel des systemes MIMO (Multiple Input Multiple Output)
3. Connexion avec l'informatique : structure matricielle facilite l'implementation numerique
4. Base pour sujets avances : commande optimale (LQR/LQG), commande robuste (H-infini), commande predictive (MPC)

Apports du projet robot NXT Way :

Le projet de robot auto-equilibre a ete particulierement formateur :

• Application concrete : pendule inverse = probleme classique d'automatique
• Systeme instable : necessite absolue d'une commande performante
• Multidisciplinarite : mecanique + electronique + automatique
• Defi technique : equilibrage en temps reel avec contraintes materielles

Connexions avec autres cours :

L'automatique est une discipline transversale qui s'appuie sur de nombreux autres cours :

• Circuits et Filtres Analogiques (S5) : analyse frequentielle, diagrammes de Bode
• Fondements Electronique Numerique (S5) : echantillonnage pour commande numerique
• Systemes Boucles (S5) : asservissements, correcteurs PID (approche complementaire)
• Commande Numerique (S8) : discretisation, commande echantillonnee, microcontroleurs
• Temps Reel (S8) : implementation de lois de commande avec contraintes temporelles
• Processus Stochastiques (S8) : filtre de Kalman (extension stochastique de l'observateur)

Evolution technologique :

L'automatique evolue rapidement avec les nouvelles technologies :

Tendances actuelles :

• Apprentissage automatique : hybridation commande classique + IA (Deep Reinforcement Learning)
• Commande predictive : MPC (Model Predictive Control) dans l'automobile et l'industrie
• Systemes distribues : commande multi-agents, essaims de drones
• Optimisation en temps reel : calcul embarque performant (DSP, FPGA)

Applications emergentes :

• Vehicules autonomes : controle longitudinal/lateral, planification de trajectoire
• Drones : stabilisation, suivi de trajectoire, vol en formation
• Exosquelettes : assistance robotique, reeducation
• Industrie 4.0 : cobotique (robots collaboratifs), maintenance predictive

Recommandations pour reussir :

1. Comprendre avant de calculer : visualiser le comportement physique du systeme
2. Maitriser MATLAB : indispensable pour l'automatique (industrie + recherche)
3. Penser "etats" : raisonner en variables d'etat plutot qu'en entrees/sorties
4. Simuler systematiquement : valider les calculs avant implementation reelle
5. Tenir compte des limitations : saturation, bruit, echantillonnage

Applications professionnelles :

Ces competences sont tres recherchees dans de nombreux secteurs :

Automobile :

• Systemes ADAS (Advanced Driver Assistance Systems) : ESP, ABS, regulateur adaptatif
• Vehicules autonomes : controle de trajectoire, fusion de capteurs
• Motorisation hybride/electrique : gestion d'energie, controle de traction

Aeronautique/Spatial :

• Pilotage automatique d'avions, drones, fusees
• Stabilisation de satellites, controle d'attitude
• Lanceurs : guidage, navigation, controle

Robotique :

• Robots mobiles : navigation, evitement d'obstacles
• Robots manipulateurs : controle de trajectoire, force
• Robots humanoides : equilibrage, marche dynamique

Industrie :

• Automatisation de procedes : chimie, petrochimie, agroalimentaire
• Machines-outils : positionnement precis, usinage
• Energies renouvelables : MPPT pour PV, controle d'eoliennes

Perspectives de carriere :

• Ingenieur automaticien : conception de systemes de commande
• Ingenieur systeme embarque : implementation temps reel
• Ingenieur R&D : algorithmes de commande avancee
• Consultant : expertise en automatique pour divers secteurs

---

Bilan personnel : Ce cours a fourni les outils mathematiques et pratiques pour analyser et commander des systemes dynamiques lineaires. L'approche etat constitue une methodologie puissante pour les systemes multivariables et la commande moderne. Les TPs MATLAB et le projet robot ont permis de concretiser ces concepts theoriques sur des applications reelles d'automatique. C'est un cours exigeant mais absolument fondamental pour une carriere en automatique, robotique, ou systemes embarques.

H(s) = Km / (Tm*s + 1)

dx/dt = A*x + B*u  (state equation)
y = C*x + D*u      (output equation)

A = [0 1; -10 -6.316]
B = [0; 1]
C = [k/0.1001 0]
D = 0

function RA=routh(poli,epsilon)
% Computes the Routh table for a polynomial
% Handles special cases (zero first column, row of zeros)

Re(pi) < 0 for all i

Gain(dB) = 20*log10(|H(jw)|)
Phase(deg) = arg(H(jw)) * 180/pi

H = tf([0.717], [0.0033 1]);
bode(H)

[Gm, Pm, Wcg, Wcp] = margin(sys);

K = place(A, B, poles_desired)

K = place(A, B, [-2+2i, -2-2i])

Mc = [B | A*B | A^2*B | ... | A^(n-1)*B]

Mc = ctrb(A, B);
rank_Mc = rank(Mc);
% System is controllable if rank_Mc == n

dx_hat/dt = A*x_hat + B*u + L*(y - y_hat)
y_hat = C*x_hat

Mo = [C; C*A; C*A^2; ...; C*A^(n-1)]

Mo = obsv(A, C);
rank_Mo = rank(Mo);
% System is observable if rank_Mo == n

L1 = [-0.365, -0.78]   % Slow poles
L2 = [0.59, 0.32]      % Fast poles

A = P * Lambda * P^(-1)

dz/dt = Lambda*z + P^(-1)*B*u
y = C*P*z

% Experimental data
Vg = [-3.93, -2.39, -1.17, -0.44, 0, 0.37, ...];
O = [-292*(2*pi/60)*9, -179*(2*pi/60)*9, ...]; % Speed in rad/s

% Static characteristic plot
plot(O, Vg);
xlabel('Angular velocity (rad/s)');
ylabel('Voltage (V)');

% Frequency analysis
G = [0.669, 0.659, 0.609, 0.528, ...];
G_dB = 20*log10(G);
freq = [0.05, 0.1, 0.25, 0.5, 1, 2.5, 5, 7.5, 10];

semilogx(freq, G_dB);
xlabel('Frequency (Hz)');
ylabel('Gain (dB)');

Ks = 1.55;   % Sensor gain
Km = 47.8;   % Motor gain
Kg = 0.015;  % Gear ratio gain
Tm = 0.3;    % Motor time constant
K = 0.887;   % Overall gain

% Closed-loop transfer function
num = [Km*K/9];
den = [Tm, 1, Km*K*Ks/9];
sys = tf(num, den);

% Stability analysis
[Gm, Pm, Wcg, Wcp] = margin(sys);
disp(['Gain margin: ', num2str(Gm), ' (', num2str(20*log10(Gm)), ' dB)']);
disp(['Phase margin: ', num2str(Pm), ' deg']);

% Pole computation
P = pole(sys);
disp('System poles:');
disp(P);

k = 5;
A = [0 1; -10 -6.316];
B = [0; 1];
C = [k/0.1001 0];

% Create state-space model
sys = ss(A, B, C, 0);

% Pole placement at -2 +/- 2j
poles_desired = [-2+2i, -2-2i];
K = place(A, B, poles_desired);

disp('State-feedback gain K:');
disp(K);

% Closed-loop system
A_cl = A - B*K;
sys_cl = ss(A_cl, B, C, 0);

% Step response
step(sys_cl);
title('Closed-loop step response with state feedback');

% Model creation
sys = tf(num, den);              % Transfer function
sys = ss(A, B, C, D);            % State-space representation

% Conversion
[A, B, C, D] = tf2ss(num, den);  % TF -> State-space
[num, den] = ss2tf(A, B, C, D);  % State-space -> TF

% Analysis
pole(sys);                       % Pole computation
zero(sys);                       % Zero computation
eig(A);                          % Eigenvalues of A
rank(ctrb(A, B));               % Controllability test
rank(obsv(A, C));               % Observability test

% Time responses
step(sys);                       % Step response
impulse(sys);                    % Impulse response
lsim(sys, u, t);                % Response to arbitrary input

% Frequency analysis
bode(sys);                       % Bode diagrams
nyquist(sys);                    % Nyquist diagram
nichols(sys);                    % Black-Nichols diagram
margin(sys);                     % Stability margins

% Synthesis
K = place(A, B, poles);          % Pole placement
L = place(A', C', poles)';       % Observer (duality)

H(s) = wn^2 / (s^2 + 2*zeta*wn*s + wn^2)

p1,2 = -zeta*wn +/- j*wn*sqrt(1-zeta^2)

J = integral[x'Qx + u'Ru] dt

dx/dt = f(x, u) ~ f(x0, u0) + df/dx|(x0,u0) * (x-x0) + df/du|(x0,u0) * (u-u0)