Archives pour la catégorie Théorie

Capteur de pression absolu, différentiel ou jauge ?

Quelle est la différence entre un capteur de pression absolu, différentiel ou jauge ?


  • Les capteurs de pression absolus mesurent une pression référencée par rapport au vide absolu, égale à la pression de jauge plus la pression atmosphérique.

  • Les capteurs de pression différentiels mesurent une différence de pression entre deux points.

  • Les jauges mesure une pression référencée par rapport à la pression atmosphérique (ambiante ?), égale à la pression absolue moins la pression atmosphérique.

Régression linéaire

Le but de ce poste est d’expliquer comment calculer une régression linéaire avec un exemple très simple. Supposons que l’on veuille approximer un nuage de points avec une droite d’équation y=ax comme illustré ci-dessous:

linear_regression

Nous voulons minimiser l’erreur donnée par la somme des distances au carré entre les points et la droite y=ax:

 E=(y_1-a.x_1)^2+(y_2-a.x_2)^2+...+(y_n-a.x_n)^2

En développant chaque terme, l’équation devient :

 E=y_1^2+a^2.x_1^2-2a.y_1.x_1 + y_2^2+a^2.x_2^2-2a.y_2.x_2 + ... + y_n^2+a^2.x_n^2-2a.y_n.x_n
 E=y_1^2 + y_2^2 + ... + y_n^2 + a^2(x_1^2 +x_2^2 +...+x_n^2) - 2a(y_1.x_1 + y_2.x_2 +y_n.x_n)

Comme nous voulons minimiser cette erreur, la valeur de a que nous cherchons est nécessairement un zéro de la dérivée de l’erreur. Calculons cette dérivée :

 \frac{\partial E}{\partial a} = 2a(x_1^2 +x_2^2 +...+x_n^2) - 2(y_1.x_1 + y_2.x_2 +y_n.x_n)

L’erreur est minimisée quand  \frac{\partial E}{\partial a} =0 :

 2a(x_1^2 +x_2^2 +...+x_n^2) - 2(y_1.x_1 + y_2.x_2 +y_n.x_n)=0

Il est maintenant facile de calculer a :

 a=\frac{2(y_1.x_1 + y_2.x_2 +y_n.x_n)}{ 2(x_1^2 +x_2^2 +...+x_n^2) }
 a=\frac{ y_1.x_1 + y_2.x_2 +y_n.x_n}{ x_1^2 +x_2^2 +...+x_n^2 }
 a=\frac { \sum\limits_{i=1}^n x_i.y_i } {\sum\limits_{i=1}^n x_i^2}

Téléchargements

Quaternions et rotations

Un quaternion est un nombre complexe de dimension 4 qui peut être utilisé pour représenter l’orientation d’un corps rigide dans un espace à trois dimensions. La définition générale d’un quaternions est donnée par :

 Q=a+b.i+c.j+d.k = \left[ \begin{matrix}a && b && c  && d\end{matrix}\right]

Representation

Considérons un vecteur \vec{V} défini par trois coordonnées (V_x, V_y and V_z) et \theta un angle de rotation autour de \vec{V} :

quaternions

Le quaternion associé à cette transformation est donné par :

 Q = \left[\begin{matrix}cos \frac{\theta}{2} && -V_x.sin \frac{\theta}{2} &&-V_y.sin \frac{\theta}{2} &&-V_z.sin \frac{\theta}{2} \end{matrix}\right]

Rotation autour des axes primaires


En s’appuyant sur la formule précédente, nous pouvons calculer le quaternion qui va définir la rotation autour de chaque axe :

Rotation autour de X

 Q_X=\left[\begin{matrix}cos \frac{\theta}{2} && -sin \frac{\theta}{2} &&0 &&0\end{matrix}\right]

Rotation autour de Y

 Q_Y=\left[\begin{matrix}cos \frac{\theta}{2} && 0 &&-sin \frac{\theta}{2} &&0\end{matrix}\right]

Rotation autour de Z

 Q_Z=\left[\begin{matrix}cos \frac{\theta}{2} && 0 &&0 &&-sin \frac{\theta}{2} \end{matrix}\right]

Rotation généralisée

Un vecteur défini en trois dimensions (ou un point) peut subir une rotation définie par un quaternion en utilisant la relation basée sur le produit de quaternion et son conjugué :

 V_B = {}^BQ_A \otimes V_A \otimes {}^BQ_A^*

V_A et V_B sont les quaternions associés aux vecteurs \vec{v_A} et \vec{v_B}. Chaque quaternion contient un 0 inséré à la place du premier élément (angle de rotation) afin de créer un quaternion avec 4 éléments, les autres éléments sont donnés par les coordonnées des vecteurs.
\vec{v_A} et \vec{v_B} est le même vecteur explicité respectivement dans les repères A et B.

Exemple

Cet exemple montre la rotation d’une éllipse définie dans le plan {\vec{X},\vec{Y}} :

elipse_rotation

Cet exemple utilise la bibliothèque de quaternions pour Matlab du Dr Madgwick’s Matlab:

Quaternion conjugué

Le conjugué d’un quaternion est représenté par le symbole ^*. Considérons le quaternion Q défini par:

 Q = \left[ \begin{matrix} a && b && c && d \end{matrix} \right]

Le conjugué de Q est défini par:

 Q^* = \left[ \begin{matrix} a && -b && -c && -d \end{matrix} \right]

Le conjugué peut être utilisé pour permuter les repères relatifs lors de la définition d’une orientation. Par exemple, si {}^BQ_A décrit l’orientation du repère B par rapport au repère A, le conjugué de {}^BQ_A^*={}^AQ_B est décrit l’orientation du repère A par rapport au repère B.

Filtre de Kalman Etendu (EKF)

Cet article détaille les équations du filtre de Kalman

Prediction

Prédiction de l’état :

 \hat{x}_{i}^p = f(\hat{x}_{i-1},u_i)

Où:


  • \hat{x}_i^p est a prédiction de l’état à l’instant i.
  • \hat{x}_{i} est l’état estimé à l’instant i.
  • f est une fonction dérivable qui prédit l’évolution de l’état en fonction de l’état précédent (prediction) et de l’entrée du système (u_i)
  • u_i est l’entrée du système à l’instant i.

Prédiction de l’incertitude (ou de la covariance) :

 P_i^p = F_i.P_{i-1}.F_i^\top + Q_i

Où:


  • P_i^p est la matrice de covariance associée à l’état \hat{x}_i^p à l’instant i.
  • P_i est la matrice de covariance associée à l’état \hat{x}_{i}.
  • Q_i est la covariance liée à la prédiction (ou bruit) du système.
  • F_i est la matrice de transition. Elle est donnée par la Jacobienne de f()=\frac{df}{dx}

Mise à jour

Innovation or erreur résiduelle :

 \widetilde{y}_i = z_i - h(\hat{x}_i^p)

Où:


  • \widetilde{y}_i est l’erreur de mesure : c’est la différence entre la mesure (ou observation) z_i et la mesure estimée depuis l’état prédit \hat{x}_i^p.
  • z_i est l’observation (or mesure) depuis l’état réel x_i.
  • h est une fonction dérivable qui permet de passer de l’espace d’état vers l’espace des observations.

Covariance de l’innovation résiduelle :

 S_i=H_i.P_i^p.H_i^\top+R_i

Où:


  • S_i est la matrice de covariance associée à la mesure d’erreur \widetilde{y}_i.
  • R_i est la matrice de covariance associée au bruit de mesure.
  • H_i est une matrice de transition qui qui permet de passer de l’espace d’état vers l’espace des observations. Elle est donné epar la Jacobienne de H_i=\frac{dh}{dx}.

Gain optimal de Kalman :

 K_i=P_i^p.H_i^\top.S_i^{-1}



  • K_i est le gain de Kalman. Cette matrice contient la pondération entre les prédictions et les mesures. C’est cette matrice qui va permettre de réaliser la fusion entre la prédiction et les différentes observations.

Mise à jour de l’état estimé :

 \hat{x}_i=\hat{x}_i^p + K_i.\widetilde{y}_i

Mise à jour de la covariance estimée

 P_i=(I-K_i.H_i).P_i^p

Où:


  • I est la matrice identité.

Filtre de Kalman

Cet article détaille les équations du filtre de Kalman

Prediction

Prédiction de l’état :

 \hat{x}_{i}^p = F_i.\hat{x}_{i-1} + B_i.u_i

Où:


  • \hat{x}_i^p est a prédiction de l’état à l’instant i.
  • \hat{x}_{i} est l’état estimé à l’instant i.
  • F_i la matrice de transition. Elle prédit l’évolution de l’état en fonction de l’état précédent (prediction).
  • B_i est la matrice qui va permettre de déterminer l’influence de l’entrée du système à l’instant i sur l’état du système. En d’autre mots, elle permet la transition de l’espace d’entré (vecteur u_i) vers l’espace d’état.
  • u_i est l’entrée du système à l’instant i.

Prédiction de l’incertitude (ou de la covariance) :

 P_i^p = F_i.P_{i-1}.F_i^\top + Q_i

Où:


  • P_i^p est la matrice de covariance associée à l’état \hat{x}_i^p à l’instant i.
  • P_i est la matrice de covariance associée à l’état \hat{x}_{i}.
  • Q_i est la covariance liée à la prédiction (ou bruit) du système

Mise à jour

Innovation or erreur résiduelle :

\widetilde{y}_i = z_i - H_i.\hat{x}_i^p

Où:


  • \widetilde{y}_i est l’erreur de mesure : c’est la différence entre la mesure (ou observation) z_i et la mesure estimée depuis l’état prédit \hat{x}_i^p.
  • z_i est l’observation (or mesure) depuis l’état réel x_i.
  • H_i est une matrice de transition qui permet de passer de l’espace d’état vers l’espace des observations.

Covariance de l’innovation résiduelle :

 S_i=H_i.P_i^p.H_i^\top+R_i

Où:


  • S_i est la matrice de covariance associée à la mesure d’erreur \widetilde{y}_i.
  • R_i est la matrice de covariance associée au bruit de mesure.

Gain optimal de Kalman :

 K_i=P_i^p.H_i^\top.S_i^{-1}



  • K_i est le gain de Kalman. Cette matrice contient la pondération entre les prédictions et les mesures. C’est cette matrice qui va permettre de réaliser la fusion entre la prédiction et les différentes observations.

Mise à jour de l’état estimé :

 \hat{x}_i=\hat{x}_i^p + K_i.\widetilde{y}_i

Mise à jour de la covariance estimée

 P_i=(I-K_i.H_i).P_i^p

Où:


  • I est la matrice identité.

Exemple de filtre de Kalman

Le filtre de Kalman est un outil mathématique couramment utilisé pour réaliser de la fusion de données provenant de différents capteurs. We allons illustré son utilisation sur une exemple simple permettant de comprendre comment le filtre fonctionne. COnsidérons un robot situé en face d’un mur qui ne peut se déplacer que dans une seule direction. Supposons que le robot est équipé de deux capteurs : un capteur de vitesse et un capteur de distance (range finder). Nous supposerons que les deux capteurs sont bruité.

Robot

Le but est ici d’estimer, le plus précisément possible, la position x du robot:

position

Les entrées du système sont une mesure de distance bruitée et une mesure de vitesse elle aussi buitée:

rangefinder

Velocity

Resultat

result

Code source

close all;
clear all;
clc;

% Step time
dt=0.1;
% Variance of the range finder
SensorStd=0.08;
% Variance of the speed measurement
SpeedStd=0.5;

%% Compute the robot trajectory and simulate sensors
% Time
t=[0:dt:10];
n=size(t,2)
% Trajectory
x=0.05*t.*cos(2*t);
% Accelerometer
u=[0,diff(x)/dt] + SpeedStd*randn(1,n);
% Range finder
z=x-SensorStd*randn(1,n);



% Display position
figure;
plot (t,x);
grid on;
xlabel ('Time [s]');
ylabel ('Position [m]');
title ('Robot real position versus time x=f(t)');

% Display accelation
figure;
plot (t,u,'m');
grid on;
hold on;
plot (t,[0,diff(x)/dt]+3*SpeedStd,'k-.');
plot (t,[0,diff(x)/dt]-3*SpeedStd,'k-.');
xlabel ('Time [s]');
ylabel ('Velocity [m/s]');
title ('Velocity versus time u=f(t)');

% Display range finder
figure;
plot (t,z,'r');
grid on;
hold on;
plot (t,x+3*SensorStd,'k-.');
plot (t,x-3*SensorStd,'k-.');
xlabel ('Time [s]');
ylabel ('Distance [m]');
title ('Distance versus time z=f(t)');


%% Kalman

B=dt;
F=1;
Q=SpeedStd^2;
H=1;
R=SensorStd^2;

x_hat=[0];
y_hat=[0];
P=[0];
S=[0];

for i=2:n
        
    
    %% Prediction
    
    %State
    x_hat(i)=F*x_hat(i-1) + B*u(i-1);
    
    % State uncertainty
    P(i)=F*P(i-1)*F' + Q;
    
    
    
    %% Update    
    
    % Innovation or measurement residual
    y_hat(i)=z(i)-H*x_hat(i);
    
    % Innovation (or residual) covariance
    S=H*P(i)*H'+R;
    
    K=P(i)*H'*S;
    
    x_hat(i)=x_hat(i)+K*y_hat(i);
    
    P(i)=(1-K*H)*P(i);    
end

figure;
plot (t,x,'k');
hold on;
plot (t,x_hat);
plot (t,x_hat+sqrt(P),'r-.');
plot (t,x_hat-sqrt(P),'r-.');
legend('Real position','Estimated position','Uncertainty');

grid on;
xlabel ('Time [s]');
ylabel ('Distance [m]');
title ('Estimated position (speed only)');

Téléchargement

Cet example a été inspiré par l’excellent tutorial de Bradley Hiebert-Treuer “An Introduction to Robot SLAM (Simultaneous Localization And Mapping)

Loi normale

La loi normale (or Gaussienne) est un modèle statistique classiquement utilisé pour estimer l’incertitude d’un capteur. La loi est donné par la formule suivante:

 y(x)=\frac{1}{\sigma.\sqrt{2.\pi}} e^ {-\frac{(x-c)^2}{2\sigma^2}  }

où:


  • \sigma est l’écart-type (racine carré de la variance)
  • c est le centre de la gaussienne

variance=\sigma^2

Voici quelques exemples de lois normales:
Examples

La somme des probabilités doit être égale à un, par conséquent, l’aire suivante a une superficie de un:

Surface

68, 95 et 99.7% de la surface sont inclus respectivement entre \sigma, 2.\sigma and 3.\sigma:
distribution

Décomposition en valeurs singulières (SVD) d’une matrice 2×2

Cette page détaille la décomposition en valeurs singulières (SVD) d’une matrice 2X2. Nous rappelons:

 A=\left[ \begin{matrix} a & b \\ c & d \end{matrix} \right]=U  \cdot \Sigma  \cdot  V^T

U et V sont orthogonales et \Sigma est une matrice diagonale contenant les valeurs singulières. Dans le cas particulier d’une matrice 2×2 la décomposition est donnée par:

 A=\left[ \begin{matrix} cos(\theta) & -sin(\theta) \\ sin(\theta) & cos(\theta) \end{matrix} \right]  \cdot \left[ \begin{matrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{matrix} \right]  \cdot \left[ \begin{matrix} cos(\phi) & -sin(\phi) \\ sin(\phi) & cos(\phi) \end{matrix} \right]

avec

 \begin{matrix}V=C  \cdot W=\left[ \begin{matrix} \pm{1} & 0 \\ 0 & \pm{1} \end{matrix} \right]\end{matrix}

Calcul de U

\theta est donné par:

 \theta=\frac{1}{2} atan2 (2ac+2bd , a^2+b^2-c^2-d^2)

La matrice U est donnée par:

 U=\left[ \begin{matrix} cos(\theta) & -sin(\theta) \\ sin(\theta) & cos(\theta) \end{matrix} \right]

Calcul de \Sigma

Les valeurs singulières sont données par:

 \begin{matrix} \sigma_1=\sqrt{ \frac {S_1 + S_2}{2} } \\\sigma_2=\sqrt{ \frac {S_1 - S_2}{2} } \\\end{matrix}

avec:

 \begin{matrix} S_1 &=& a^2+b^2+c^2+d^2 \\S_2 &=& \sqrt{(a^2+b^2-c^2-d^2)^2 + 4(ac+bd)^2 } \end{matrix}

La matrice \Sigma est donnée par:

 \Sigma=\left[ \begin{matrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{matrix} \right]

Calcule de V

Commençons par calculer l’angle intermédiaire \Phi donné par la formule suivante:

 \Phi=\frac{1}{2} atan2 (2ab+2cd , a^2-b^2+c^2-d^2)

Calculons maintenant les scalaires suivants:

\begin{matrix}s_{11} &=& (a.c_\theta+c.s_\theta).c_\Phi &+& (b.c_\theta+d.s_\theta).s_\Phi \\s_{22} &=& (a.s_\theta-c.c_\theta).s_\Phi &+& (-b.s_\theta+d.c_\theta).c_\Phi \\\end{matrix}

où:

\begin{matrix}c_\Phi=cos(\Phi) \\s_\Phi=sin(\Phi) \\c_\theta=cos(\theta) \\ s_\theta=sin(\theta)\end{matrix}

La matrice V est donnée par:

 V=\left[ \begin{matrix} \frac {s_{11}} {|s_{11}|}cos(\Phi) & -\frac {s_{22}} {|s_{22}|}sin(\Phi) \\ \frac {s_{11}} {|s_{11}|} sin(\Phi) & \frac {s_{22}} {|s_{22}|}cos(\Phi) \end{matrix} \right]

Notons que \frac {s_{11}} {|s_{11}|} et \frac {s_{22}} {|s_{22}|} sont respectivement les signes de s_{11} et s_{22}. Sous Matlab, l’utilisation de la fonction sign évite le problème de la division par zéro si s_{11} ou s_{22} sont nuls.

Fonction Matlab

%% Matlab function for computing the VSD of a 2x2 matrix
%% Written by Randy Ellis
%% More information at : https://www.lucidarme.me


function [U,SIG,V] = svd2x2(A)
% [U,SIG,V] = svd2x2(A) finds the SVD of 2x2 matrix A
% where U and V are orthogonal, SIG is diagonal,
% and A=U*SIG*V’
% Find U such that U*A*A’*U’=diag
Su = A*A';
phi = 0.5*atan2(Su(1,2)+Su(2,1), Su(1,1)-Su(2,2));
Cphi = cos(phi);
Sphi = sin(phi);
U = [Cphi -Sphi ; Sphi Cphi];


% Find W such that W’*A’*A*W=diag
Sw = A'*A;
theta = 0.5*atan2(Sw(1,2)+Sw(2,1), Sw(1,1)-Sw(2,2));
Ctheta = cos(theta);
Stheta = sin(theta);
W = [Ctheta, -Stheta ; Stheta, Ctheta];


% Find the singular values from U
SUsum= Su(1,1)+Su(2,2);
SUdif= sqrt((Su(1,1)-Su(2,2))^2 + 4*Su(1,2)*Su(2,1));
svals= [sqrt((SUsum+SUdif)/2) sqrt((SUsum-SUdif)/2)];
SIG =diag(svals);


% Find the correction matrix for the right side
S = U'*A*W
C = diag([sign(S(1,1)) sign(S(2,2))]);
V = W*C;

Téléchargement

L’écriture de ce post a été inspiré par l’excellent rapport de Randy Ellis téléchargeable ci-dessous:

Déterminer la transformation entre deux nuages de points

L’objet de cet article est la minimisation d’erreur entre deux nuages de points (minimisation de l’erreur point à point). Considérons les nuages de points A en bleu et B en rouge sur la figure ci-dessous. Notre but est de trouver la transformation de B (translation et rotation) qui va permettre de minimiser les distances point à point (au sens de la minimisation des moindres carrés).

figure0

 A=\left[ \begin{matrix} p_1^a & p_2^a & ... & p_n^a \end{matrix} \right]=\left[ \begin{matrix} x_1^a & x_2^a & ... & x_n^a \\ y_1^a & y_2^a & ... & y_n^a \end{matrix} \right]
 B=\left[ \begin{matrix} p_1^b & p_2^b & ... & p_n^b \end{matrix} \right]=\left[ \begin{matrix} x_1^b & x_2^b & ... & x_n^b \\ y_1^b & y_2^b & ... & y_n^b \end{matrix} \right]

Nous supposerons que les points sont déjà appairés: chaque point de A est associé à un point de B.

figure1

Translation

Pour commencer, calculons le centre de gravité de chaque nuage de points:
figure2

 \begin{matrix}COG^a=\frac{1}{n}.  \sum\limits_{i=1}^n p_i^a & COG^b=\frac{1}{n}.  \sum\limits_{i=1}^n p_i^b \end{matrix}

Chaque point est translaté d’un vecteur donné par le centre de gravité de son nuage (cela revient à centrer chaque nuage sur l’origine du repère).

 \begin{matrix} A

figure3

Rotation

Pour la rotation, les choses sont un peu plus compliquées. D’abord, calculons la matrice N suivante :

 N= \sum\limits_{i=1}^n {p

Décomposons N en valeurs singulières :

 N=U.\Sigma.V^T

La matrice de rotation est donnée par :

 R=V.U^T

L’angle de rotation peut être extrait de la matrice R

 \alpha=atan2(R_{21},R_{11})

En appliquant la rotation sur les nuages centrés, nous obtenons le résultat suivant:

figure4

Téléchargement

Ce post a été inspiré par l’excellent rapport de Hans Martin Kjer and Jakob Wilm:

Moindres carrés : approximation avec un polynôme du second degré

Supposons que l’on veuille approximer un nuage de points avec un polynôme du second degré : y(x)=ax^2+bx+c.

SecondDegreePolynomial

Le nuage de points est donné par un ensemble de n points ayant pour coordonnées {x_i,y_i}. L’objectif est d’estimer \hat{a} et \hat{b} de façon à ce que le polynome d’équation y(x)=\hat{a}x^2+\hat{b}x + \hat{c} approche au mieux le nuage de points. Pour chaque point x_i nous allons minimiser la différence entre y_i et y(x_i), c-a-d. nous voulons minimiser \sum\limits_{i=1}^n{(y_i-y(x_i))^2}. La formulation matricielle du système est donnée par:

 \left[ \begin{matrix} {x_1}^2 & x_1 & 1 \\{x_2}^2 & x_2 & 1 \\... & ... & ... \\{x_n}^2 & x_n & 1 \\\end{matrix} \right].\left[ \begin{matrix} \hat{a} \\\hat{b} \\\hat{c} \end{matrix} \right]=\left[ \begin{matrix} y_1 \\y_2 \\... \\y_n \\\end{matrix} \right]

Définissons A, B et \hat{x}:

\begin{matrix}  A=\left[ \begin{matrix} {x_1}^2 & x_1 & 1 \\{x_2}^2 & x_2 & 1 \\... & ... & ... \\{x_n}^2 & x_n & 1 \\\end{matrix} \right] & B=\left[ \begin{matrix} y_1 \\y_2 \\...  \\y_n \\\end{matrix} \right] &\hat{x}=\left[ \begin{matrix} \hat{a} \\ \hat{b} \\ \hat{c} \end{matrix} \right] \end{matrix}

Le système est maintenant donné par:

 A.\hat{x}=B

La solution optimale est :

 \hat{x}=A^{+}.B = A^{T}(A.A^{T})^{-1}.B

A^{+} est la pseudo-inverse de A donnée par A^{+}=A^{T}(A.A^{T})^{-1}.

Téléchargement

Approximation linéaire par les moindres carrès

Cet article explique comment réaliser une approximation linéaire sur un exemple simple. Supposons que l’on veuille approximer un nuage de points avec une droite : y(x)=ax+b.

LinearApproximation

Le nuage de points est donné par un ensemble de n points ayant pour coordonnées {x_i,y_i}. L’objectif est d’estimer \hat{a} et \hat{b} de façon à ce que la droite d’équation y(x)=\hat{a}x+\hat{b} approche au mieux le nuage de points. Pour chaque point x_i nous allons minimiser la différence entre y_i et y(x_i), c-a-d. nous voulons minimiser \sum\limits_{i=1}^n{(y_i-y(x_i))^2}. La formulation matricielle du système est donnée par:

 \left[ \begin{matrix} x_1 & 1 \\x_2 & 1 \\... & ... \\x_n & 1 \\\end{matrix} \right].\left[ \begin{matrix} \hat{a} \\\hat{b}\end{matrix} \right]=\left[ \begin{matrix} y_1 \\y_2 \\... \\y_n \\\end{matrix} \right]

Définissons A, B et \hat{x}:

 \begin{matrix}A=\left[ \begin{matrix} x_1 & 1 \\ x_2 & 1 \\ ... & ... \\ x_n & 1 \\ \end{matrix} \right] &B=\left[ \begin{matrix} y_1 \\y_2 \\...  \\y_n \\\end{matrix} \right] &\hat{x}=\left[ \begin{matrix} \hat{a} \\\hat{b}\end{matrix} \right] \end{matrix}

Le système est maintenant donné par:

 A.\hat{x}=B

La solution optimale est donnée par:

 \hat{x}=A^{+}.B = A^{T}(A.A^{T})^{-1}.B

A^{+} est la pseudo-inverse de A donnée par A^{+}=A^{T}(A.A^{T})^{-1}.

Téléchargement

Transmission par bielles

Cet article expose la mise en équation d’une transmission par bielle. Nous considérerons le modèle suivant:

Drawing

Cinématique

Cette partie démontre la relation qui relie les angles des axes: \Theta_2=f(\Theta_1).

DrawingKinematic

L’angle \Theta_2 va être calculé avec le relation suivante:

 \Theta_2=\pi- \widehat{P_2MP_1} -\widehat{P_1MO}

Commençons par calculer les coordonnées de P_1:

 P_1=\begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} l_1.cos(\Theta_1) \\ l_1.sin(\Theta_1) \end{pmatrix}

En se basant sur le résultat précédent, calculons la distance entre P_1 et M:

 \| P_1M \| = \sqrt{ (x_1-X)^2 + y_1^2 }

Dans le triangle HP_1M, rectangle en H, il devient aisé de calculer l’angle \widehat{P_1MO}:

 \widehat{P_1MO}=atan2(\|HP_1\|,\|MH\|)=atan2(y_1,X-x_1)

Maintenant nous pouvons calculer l’angle \widehat{P_2MP_1} dans le triangle P_1MP_2 grâce au théorème d’Al Kashi

 \|P_1P_2\| = \|MP_1\|^2 + \|MP_2\|^2 -2.\|MP_1\|.\|MP_2\|.cos(\widehat{P_2MP_1})

 L = \|MP_1\|^2 + {l_2}^2 -2.\|MP_1\|.l_2.cos(\widehat{P_2MP_1})

et finalement:

 cos(\widehat{P_2MP_1}) = \frac{\|MP_1\|^2 + {l_2}^2 - L} {2.\|MP_1\|.l_2}

En conclusion:

 \Theta_2=\pi \pm acos\left(\frac{ (x_1-X)^2 + y_1^2  + {l_2}^2 - L} {2.l_2.\sqrt{ (x_1-X)^2 + y_1^2 }} \right) - atan2(y_1,X-x_1)

où:

 \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} l_1.cos(\Theta_1) \\ l_1.sin(\Theta_1) \end{pmatrix}

Couple

Cette seconde partie expose comment calculer le couple transmis. Supposons le l’axe moteur soit celui de gauche (centré en O). Le couple fourni sur cet axe est T_1. L’objectif est de calculer le couple transmis (T_2) sur le second axe (centré en M). Le but de cette seconde partie est d’exprimer T_2=g(T_1).

DrawingTorque

Commençons par calculer les coordonnées de P_2:

 P_2=\begin{pmatrix} x_2 \\ y_2 \end{pmatrix} = \begin{pmatrix} X+l_2.cos(\Theta_1) \\ l_2.sin(\Theta_1) \end{pmatrix}

Calculons l’amplitude de la force \vec{F_1}:

 \| \vec{F_1}\| = \frac{T_1}{l_1}

Maintenant, calculons l’angle \alpha entre \vec{F_1} et \vec{F_b}:

 \alpha = ( \widehat{ \vec{F_1} ; \vec{F_b} } )= ( \widehat{ \vec{OM} ; \vec{F_1} } ) - ( \widehat{ \vec{OM} ; \vec{F_2} } ) = ( \widehat{ \vec{OM} ; \vec{F_1} } ) -( \widehat{ \vec{OM} ; \vec{P_1P_2} } )

Calculons les angles:

 ( \widehat{ \vec{OM} ; \vec{F_1} } ) = \Theta_1 - \frac{\pi}{2}

et :

 ( \widehat{ \vec{OM} ; \vec{P_1P_2} } ) = atan2(y_2-y_1,x_2-x_1)

Maintenant \alpha peut être calculé:

 \alpha = \Theta_1 - \frac{\pi}{2} - atan2(y_2-y_1,x_2-x_1)  = \Theta_1 - atan2(x_2-x_1,y_1-y_2)

et F_b peut aussi être calculé:

 \| \vec{F_b} \| = \| \vec{F_1} \|.cos(\alpha)

Comme nous l’avons fait pour \alpha, calculons l’angle \beta entre \vec{F_b} et \vec{F_2}:

 \beta = ( \widehat{ \vec{F_2} ; \vec{F_b} } )= \Theta_2 - \frac{\pi}{2} - atan2(y_2-y_1,x_2-x_1)  =  \Theta_2 - atan2(x_2-x_1,y_1-y_2)

La force F_2 peut aussi être calculée:

 \| \vec{F_2} \| = \| \vec{F_b} \|.cos(\beta)

Le couple de sortie est :

 T_2=\| \vec{F_2} \|.l_2

En conclusion:

 T_2= \frac{l_2}{l_1}. T_1.cos(\beta).cos(\alpha)

où:

 \alpha = \Theta_1 - atan2(x_2-x_1,y_1-y_2)

et

 \beta =  \Theta_2 - atan2(x_2-x_1,y_1-y_2)

Téléchargment

Ball and Beam : modèle

Cet article décrit le modèle du Ball and Beam. Ce système est composé d’une longue barre motorisée qui peut s’incliner sur laquelle une bille roule d’un sens ou de l’autre.

BallAndBeamModel

Nous considérerons les hypothèses suivantes:

  • m est la masse de la bille

  • \Theta est l’angle de la barre

  • x est la position de la bille sur la barre

  • b est le coefficient de frottement de la bille


En s’appuyant sur le principe fondamental de la dynamique et en considérant que les frottements sont propositionnels à la vitesse de la bille, nous pouvons écrire:

 m.a = m.\ddot{x} = \sum{\vec{F}} = m.g.sin(\Theta) - b.\dot{x}

En conclusion, l’équation différentielle du système est:

 \ddot{x} = - \frac{b}{m}.\dot{x} + g.sin(\Theta)

Pendule inversé: modèle sans frottements

Introduction

Cet article présente le modèle d’un pendule inversé. Ce système est composé d’un chariot et d’un pendule dont le centre de gravité est au dessus de la liaison pivot qui le lie au chariot.

InvertedPendulum

Dans cet article, nous ferons les hypothèses suivantes:


  • m_1 est la masse du chariot.

  • m_2 est la masse du pendule.

  • I_1 est le moment d’inertie du chariot.

  • I_2 est le moment d’inertie du pendule.

  • L est la distance entre la liaison pivot et le centre de gravité du pendule.

  • \Theta est l’angle d’inclinaison du pendule.

  • x_1 et y_1 sont les coordonnées du centre de gravité C_1 du chariot.

  • \Theta_1 est l’orientation du chariot.

  • x_2 et y_2 sont les coordonnées du centre de masse C_2 du pendule.

  • \Theta_2 est l’orientation du pendule.

Les notations suivantes seront utilisée dans la suite de cet article pour éviter de surcharger les équations:

 \frac{d\Theta}{dt} = \dot{\Theta}

Et pour la dérivée seconde:

 \frac{d^2\Theta}{dt^2} = \ddot{\Theta}

Contraintes

Comme le chariot et le pendule sont liés par une liaison pivot, nous pouvons écrire les équations suivantes. Ces équations expriment les coordonnées et l’orientation de chaque corps du système. Pour le chariot (corps 1):

\begin{array}{r c l} x_1 &=& x_1 \\ y_1 &=& 0 \\ \Theta_1 &=& 0 \\\end{array}

Pour le pendule (corps 2):

\begin{array}{r c l} x_2 &=& x_1 - L.sin(\Theta)\\ y_2 &=& y_1 + L.cos(\Theta) &=& L.cos(\Theta)\\ \Theta_2 &=& \Theta\end{array}

Définissons \vec{q}, la position de notre système :

 \vec{q}=\begin{pmatrix} x_1 \\ \Theta_2 \end{pmatrix}

Définissons les contraintes sur la position:

 \phi = \begin{pmatrix} y_1 \\ \Theta_1 \\ x_2 \\ y_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ x_1 - L.sin(\Theta) \\ L.cos(\Theta) \end{pmatrix}

Calculons les contraintes en vitesse:

 \dot{\phi} = \begin{pmatrix} \dot{y_1} \\ \dot{\Theta_1} \\ \dot{x_2} \\ \dot{y_2} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\  \dot{x_1} - L.\dot{\Theta_2}.cos(\Theta_2)  \\ -L.\dot{\Theta_2}.sin(\Theta_2) \end{pmatrix}

Calculons les contraintes en accélération:

 \ddot{\phi} = \begin{pmatrix} \dot{y_1} \\ \ddot{\Theta_1} \\ \ddot{x_2} \\ \ddot{y_2} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\  \ddot{x_1} - L.\ddot{\Theta_2}.cos(\Theta_2) + L.\dot{\Theta_2}^2.sin(\Theta_2) \\ -L.\ddot{\Theta_2}.sin(\Theta_2) - L.\dot{\Theta_2}^2.cos(\Theta_2) \end{pmatrix}

Dynamique

Le principe fondamental de la dynamique peut être appliqué à notre système. Pour le chariot (corps 1):

\begin{array}{r c l} m_1.{a_1}^x = m_1.\ddot{x_1}&=& \| \vec{F_t} \| + \lambda_1 \\ m_1.{a_1}^y = m_1.\ddot{y_1}&=& 0 + \lambda_2 \\ I_1.\ddot{\Theta_1} &=& 0 + \lambda_3 \\\end{array}

Pour le pendule (corps 2):

\begin{array}{r c l} m_2.{a_2}^x = m_2.\ddot{x_2} &=& 0 + \lambda_4 \\  m_2.{a_2}^y = m_2.\ddot{y_2} &=& \| \vec{F_g} \| + \lambda_5 = -g.m_2 + \lambda_5\\ I_2.\ddot{\Theta_2} &=& 0 + \lambda_6 \\\end{array}

\lambda_1, \lambda_2,\lambda_3,\lambda_4,\lambda_5 et \lambda_6 sont les coefficients de Lagrange. Ces coefficients représentent les interactions entre les coprs. Dans notre système, il s’agit de l’interaction entre le chariot et le pendule.

Interactions

Nous allons maintenant résoudre le système de façon à éliminer les interactions entre les corps(\lambda_i). Il est prouvé que:

 M.\ddot{q} = F - \frac{d\phi}{dq}^\top .\Lambda

Où:

 M=\begin{pmatrix} m_1 & 0 \\ 0 & I_2 \end{pmatrix}
 \ddot{q}=\begin{pmatrix} \ddot{x_1} \\ \ddot{\Theta_2} \end{pmatrix}
 F=\begin{pmatrix} \| \vec{F_t} \|  \\ 0 \end{pmatrix}
 \frac{d\phi}{dq} = \begin{pmatrix} \frac{dy_1}{dx1} & \frac{dy_1}{d\Theta_2} \\\frac{d\Theta_1}{dx1} & \frac{d\Theta_1}{d\Theta_2} \\\frac{dx_2}{dx1} & \frac{dx_2}{d\Theta_2} \\\frac{dy_2}{dx1} & \frac{dy_2}{d\Theta_2} \end{pmatrix} =\begin{pmatrix}0 & 0 \\0 & 0 \\1 & -L.cos(\Theta_2) \\0 & -L.sin(\Theta_2)\end{pmatrix}
 \Lambda=\begin{pmatrix} \lambda2 \\ \lambda3 \\ \lambda4 \\ \lambda5 \end{pmatrix}

Appliqué à notre système, cela donne:

\begin{pmatrix} m_1 & 0 \\ 0 & I_2 \end{pmatrix} .\begin{pmatrix} \ddot{x_1} \\ \ddot{\Theta_2} \end{pmatrix} =\begin{pmatrix} \| \vec{F_t} \|  \\ 0 \end{pmatrix} -\begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & -L.cos(\Theta_2) & -L.sin(\Theta_2) \end{pmatrix} .\begin{pmatrix} \lambda2 \\ \lambda3 \\ \lambda4 \\ \lambda5 \end{pmatrix}

et:

\begin{array}{r c l}m_1.\ddot{x_1} &=& \| \vec{F_t} \| -\lambda4 \\I_2.\ddot{\Theta_2} &=& \lambda4.L.cos(\Theta_2) + \lambda5.L.sin(\Theta_2)\end{array}

Résolution du système

En se basant sur les équations mentionnées précédement, nous pouvons reformuler l’expression de \lambda_4 et \lambda_5:

\begin{array}{r c l}\lambda_4 &=& m_2.\ddot{x_2} \\\lambda_5 &=& m_2.\ddot{y_2} + g.m_2 \\\end{array}

L’équation précédente devient:

\begin{array}{r c l}m_1.\ddot{x_1} &=& \| \vec{F_t} \| -m_2.\ddot{x_2} \\I_2.\ddot{\Theta_2} &=& m_2.\ddot{x_2}.L.cos(\Theta_2) + (m_2.\ddot{y_2} + g.m_2).L.sin(\Theta_2)\end{array}

Grâce à la contrainte en accélération:

\begin{array}{r c l}m_1.\ddot{x_1} &=& \| \vec{F_t} \| -m_2.(\ddot{x_1} - L.\ddot{\Theta_2}.cos(\Theta_2) - L.\dot{\Theta_2}^2.sin(\Theta_2) ) \\I_2.\ddot{\Theta_2} &=& m_2.(\ddot{x_1} - L.\ddot{\Theta_2}.cos(\Theta_2) + L.\dot{\Theta_2}^2.sin(\Theta_2) ).L.cos(\Theta_2) + (m_2.(-L.\ddot{\Theta_2}.sin(\Theta_2) - L.\dot{\Theta_2}^2.cos(\Theta_2)) + g.m_2).L.sin(\Theta_2)\end{array}

Reformuler l’équation précédente nous donne:

\begin{array}{r c l}(m_1+m_2).\ddot{x_1} &-& m_2.L.cos(\Theta_2).\ddot{\Theta_2} &=& \| \vec{F_t} \| - m_2.L.\dot{\Theta_2}^2.sin(\Theta_2) \\-m_2.L.cos(\Theta_2).\ddot{x_1} &+& (I_2+ m_2.L^2).\ddot{\Theta_2} &=& m_2.g.L.sin(\Theta_2)\end{array}

et:

\begin{pmatrix} m_1+m_2 & - m_2.L.cos(\Theta_2) \\-m_2.L.cos(\Theta_2) & (I_2+ m_2.L^2)\end{pmatrix} . \ddot{q} =\begin{pmatrix} \| \vec{F_t} \| - m_2.L.\dot{\Theta_2}^2.sin(\Theta_2) \\ m_2.g.L.sin(\Theta_2)\end{pmatrix}

L’équation générale du système est:

\ddot{q}=A^{-1}.B

Où:

A=\begin{pmatrix} m_1+m_2 & - m_2.L.cos(\Theta_2) \\-m_2.L.cos(\Theta_2) & (I_2+ m_2.L^2)\end{pmatrix}
B= \begin{pmatrix} \| \vec{F_t} \| - m_2.L.\dot{\Theta_2}^2.sin(\Theta_2) \\ m_2.g.L.sin(\Theta_2)\end{pmatrix}

Representation d’états

En se basant sur le résultat précédent, il est maintenant possible d’écrire la représentation d’états du système. Définissons X l’état du système:

X=\begin{pmatrix} x_1 \\ \Theta_2 \\ \dot{x_1} \\ \dot{\Theta_2} \end{pmatrix}

La représentation d’état est donc:

\dot{X}=\begin{pmatrix} \dot{x_1} \\ \dot{\Theta_2} \\ \left[A^{-1}.B\right] \end{pmatrix}

Remerciements

Je tiens à remercier Sébastien Lagrange de l’université d’Angers pour son aide et ses explications.

Pendule simple

Introduction


Un pendule simple est constitué d’un poids suspendu par un câble non déformable.

PenduleumOverview

Dans cet article, nous considérerons les hypothèses suivantes:


  • Le poids peut se balancer librement.

  • L est la longueur du câble.

  • m est la masse du corps du pendule.

  • La masse du câble est négligée.

  • Les frottements sont négligés.

  • Mise en équation du système

    Pour trouver l’équation régissant le système, nous allons appliquer le principe fondamentale de la dynamique:

     m.\vec{a} = \sum{\vec{F}}

    Les forces agissant sur le système sont la force de gravité \vec{F_g} et la force de tension du câble \vec{T}. La trajectoire de la masse décrit nécessairement un arc de cercle. L’accélération élémentaire du corps est nécessairement tangente à cette trajectoire. La force de tension du câble est donc perpendiculaire à la trajectoire de la masse. Sa projection le long de l’axe de déplacement est nulle. Une seule force devra être considérée dans la suite des calculs : la force de gravité:

     m.a = \vec{F_m} = m.g.sin(\Theta)

    La relation entre la vitesse linéaire et angulaire du corps est donnée par :

     v=L\frac{d \Theta }{dt}

    En dérivant cet équation, nous obtenons:

     a=\frac{dv}{dt}=L. \frac{d^2 \Theta }{dt^2}

    L’équation obtenue avec le principe fondamental de la dynamique peut être mise à jour:

     m.L. \frac{d^2 \Theta }{dt^2} = m.g.sin(\Theta)

    En divisant par m:

     L. \frac{d^2 \Theta }{dt^2} = g.sin(\Theta)

    Nous pouvons ainsi obtenir l’équation différentielle régissant le système:

     \frac{d^2 \Theta }{dt^2} = \frac{g}{L}.sin(\Theta)

    Résolution

    L’équation différentielle obtenue est non linéaire. Il est possible de linéariser cette équation pour de petites valeurs de \Theta en approximant sin(\Theta) par \Theta. L’équation devient linéaire:

     \frac{d^2 \Theta }{dt^2} = \frac{g}{L}.\Theta

    La solution générale à l’équation différentielle linéarisée est:

     \Theta = A.cos(\omega.t + \phi)

    avec
     \omega^2=\frac{g}{L}

    Etant donné les conditions initiales \Theta=\Theta_0 (angle initial du pendule), et en considérant que la vitesse initiale est nulle, nous pouvons écrire:

    A.cos(\phi)=\Theta_0
    -A.\omega.sin(\phi)=0

    Le système peut être résolu:

    A=\Theta_0
    \phi=0

    La solution générale du système devient:

     \Theta(t)=\Theta_0.cos(\omega.t)

    \omega = \sqrt { \frac {g}{L} }. La période d’oscillation est donnée par T=2.\pi.\sqrt { \frac {g}{L}}.

    Remerciement

    Je tiens à remercier Jean-Pierre Pecqueur de l’Université d’Angers pour son aide.