Drone | Arduino #9: Implémentation d’une intégration numérique avec Arduino – Méthode de Simpson

Objectifs

Savoir implémenter une intégration numérique avec Arduino
Savoir intégrer une fonction dans le temps
Introduction de trois techniques pour implémenter une intégration numérique avec Arduino (voir la description ci-dessous)
Savoir quelques avantages et limitations des techniques
Et d’autres astuces de programmation

Méthode de Simpson

Nous avons abordé dans la première (voir ici) les aspects théorique de méthode des trapèzes ainsi Simpson pour le calcul de l’intégrale d’une fonction f(t) dans un intervalle [a,b]. On considère T=b-a, comme étant la période d’intégration. Concernant la technique de Simpson, il faut acquérir trois valeurs de la fonction f(t) pour chaque période d’intégration T.

La valeur initiale f(a) lorsque t=0+n*T, n est une valeur entière positive
La valeur milieu f(m) : lorsque t=T/2+n*T
La valeur finale f(b) lorsque t=T+n*T

Comment obtenir la variable du temps t ?

Deux fonctions Arduino peuvent être utilisées pour générer la variable t avec une précision de l’ordre de milliseconde ou quelques microsecondes. On peut utiliser les syntaxes suivantes :

t=millis() : mesure du temps avec une précision de l’ordre de 1 ms. L’intervalle est limite à 2^32 ms (environ 49 jours, lorsque t attint 49 j il recommence à 0)
t=micros() : mesure du temps avec une précision de l’ordre de 8 µs. L’intervalle est limité à 2^32 µs (environ 71 minutes)

Comment convertir la variable t de [0, inf] en [0, t’] ? Avec t’ dans [0, T] ?

La transformation nous facilitera l’acquisition des trois valeurs f(a), f(m) et f(b) avec Arduino dans l’intervalle [0,T] quelque soit t pour le calcul de l’intégrale I. Il suffit d’appliquer l’opérateur modulo à la variable t de la façon suivante : t’=t % T (ou bien t=t % T)

Technique de l’implémentation 1 : millis() & micros()

Lecture de la variable t : t=millis() ; (un coefficient de 1000 pour millis() et 1000000 pour micros() sera appliqué durant le calcul de l’intégrale pour ramener la variable t en seconde)
Application du modulo à t : t=t%T
Lecture des valeurs f(a), f(b) et f(m) en fonction de t. On considère eps une valeur très petite par rapport à T (eps= T/10 ou eps=1/20, etc.). Elle nous permettra d’avoir une minuscule fenêtre temporelle d’acquisition des trois valeurs :
1. Acquisition de f(a) si t<eps
2. Acquisition de f(b) si T/2- eps <t< T/2+eps
3. Acquisition de f(c) si T-eps<t<T
Calcul de l’intégrale I : I=(f(a)+4*f(b)+f(c)/(6*1000) ou bien I=(f(a)+4*f(b)+f(c)/(6*1000000) dans le cas d’utilisation de la fonction micros()

Technique de l’implémentation 2 : interruption + Timer

On peut utiliser une routine d’interruption ISR pour l’acquisition et le calcul de l’intégrale. L’idée consiste à générer une interruption toute les T/2 secondes muni d’un compteur d’échantillons. Pour chaque appel à ISR on incrémente le compteur au même temps on fait l’acquittions d’une valeur unique :

Premier Appel : si count=0 alors lecture de f(a), count=1
Deuxième Appel : si count=1 alors lecture de f(m), count=2
Troisième Appel : si count=2 alors lecture de f(b), calcul de l’intégrale I, count=0
Quatrième Appel : si count=0 alors lecture de f(a), count=1
…

La technique met à jour la valeur de l’intégrale I toute les T secondes. Elle est beaucoup plus précise par rapport à la première technique. En effet, la fonction ISR s’exécute toute les T/2 avec une précision de l’ordre de µS. Contrairement à la première technique, sa précision dépend de la variable « esp » ainsi la fonction dépend de la précision de millis() ou micros(). En revanche, la deuxième technique nécessite la mise en place d’un Timer d’une valeur égale à T/2 secondes avec La fonction ISR dédiée.

Technique de l’implémentation 3 : PWM-INT0

Ici on va utiliser un signal PWM généré par la même carte Arduino qui sera ensuite utilisé pour déclencher une interruption via un pin externe de type INT ou PCINT. Par exemple si la fréquence PWM est égale à 491 Hz, alors on génère 491 interruptions par seconde ! On peut donc utiliser un compteur (Count) à l’intérieur de la fonction pour ajuster la période d’intégration ainsi le calcul de l’intégrale. Par exemple, si T= 3 secondes :

Si Count =491 : lecture de f(a)
Si Count=491*2 ; lecture de f(m)
Si Count=491*3 : lecture de f(b), calcul de l’intégrale, initiation du compteur (Count=0)
Voir l’article suivant pour plus de détails.

La fonction d’intégration

double ComputeInteg(unsigned long T_us)
{
  unsigned long Temps;
  static unsigned long f_tt[3]={0,0,0};
  unsigned long T_med=T_us>>1;
  unsigned long eps=T_us/10;
  static double int_f=0.0; 
  
  // Lecture du temps 
  Temps=micros();   // t dans [0, inf]
  Temps=Temps%T_us; // t=t%T dans [0,T]
  
  if(Temps <eps) 
    f_tt[0]=millis(); 
    
  if((Temps <T_med+eps)&&(Temps >T_med-eps))
    f_tt[1]=millis(); 
         
  if(Temps >T_us-eps)
  {
    // Leture de la dérnière valeur 
    f_tt[2]=millis(); 
  
    // Calcul de l'integrale 
    int_f=(double)T_us*(f_tt[0]+4.0*f_tt[1]+f_tt[2])/6000000.0;
  }
    
  // Retourner l'intégrale de f 
  return int_f;
}

Le programme complet


/*
I=intégrale de f(t) dans l'intervalle [a,b]
T: la période d'intégration T=b-a

Technique de Simpson: I=[f(a) + 4f(m) + f(b)]/6
                         t=0     t=T/2  t=T

Technique 1: millis() & micros()  
Technique 2: PWM+INT  
Technique 3: Timer T/2 
Technique 4: ... 

*/

double I_f=0.0, f_t=0.0;
unsigned long T_us=200000; // 1 sec
double T_sec= T_us/1000000.0;

void setup()
{
  // Init port série 
  Serial.begin(9600);
}
void loop()
{  
  // Lecture de la valeur instantanée de f  
  f_t=millis(); // f(t)=t; 
  
  // Calcul de l'intégrale 
  I_f=ComputeInteg(T_us);

  // Affichage 
  Serial.print(T_sec*f_t);Serial.print(","); 
  Serial.println(I_f); //I~0.2*Vmoy
}

double ComputeInteg(unsigned long T_us)
{
  unsigned long Temps;
  static unsigned long f_tt[3]={0,0,0};
  unsigned long T_med=T_us>>1;
  unsigned long eps=T_us/10;
  static double int_f=0.0; 
  
  // Lecture du temps 
  Temps=micros();   // t dans [0, inf]
  Temps=Temps%T_us; // t=t%T dans [0,T]
  
  if(Temps <eps) 
    f_tt[0]=millis(); 
    
  if((Temps <T_med+eps)&&(Temps >T_med-eps))
    f_tt[1]=millis(); 
         
  if(Temps >T_us-eps)
  {
    // Leture de la dérnière valeur 
    f_tt[2]=millis(); 
  
    // Calcul de l'integrale 
    int_f=(double)T_us*(f_tt[0]+4.0*f_tt[1]+f_tt[2])/6000000.0;
  }
    
  // Retourner l'intégrale de f 
  return int_f;
}