Séance 2 : Objets et Calcul numérique

Formation Python du département Mécatronique

11h00 - 12h30

Après une Séance 1 qui introduisait les bases du language Python généraliste, voici deux sujets plus spécifiques

A) Python orienté objet

  • Introduction rapide aux objets avec l'exemple d'une classe « Torseur ».
  • Méthodes et surcharge d'opérateurs (ou "comment pouvoir additioner des torseurs").

B) Python numérique avec NumPy...

  • Calcul numérique avec la bibliothèque NumPy.
  • Algèbre linéaire de base.
  • Lecture d'un fichier de données expérimentales (fichier texte type CSV)

C) ...puis tracé de graphiques avec Matplotlib

  • aperçu des fonctions graphiques principales. Déplacement/Zoom interactif.
  • exemple : diagramme de Bode (avec la question de l'identification d'un modèle)
  • bonus : diagrammes de Black, Nyquist

A) Python Orienté Objet, par l'exemple

Exemple de la création d'une Classe "Torseur2D", inspiré de la mécanique du solide.

Il existent d'autres tutoriels plus complets sur les objets en Python tel que http://docs.python.org/2/tutorial/classes.html#a-first-look-at-classes.

Cahier des charges de la Classe :

L'objet mathématique Torseur peut se caractériser par 3 paramètres :

  • une résultante (dim 2),
  • et un moment (dim 1)...
  • ...exprimé en un point (dim 2).

Pour manipuler ce type d'objet en Python, on crée donc une classe possédant ces trois champs.

Dans un second temps, on va voir comment permettre l'addition des objets Torseur2D.

A1) Construction de la classe, étape par étape

Une classe se définit par le mot-clé class suivit du nom de la classe.

On commence par une définition minimale qui inclue la méthode __init__ qui sert à initialiser l'objet :

In [9]:
class Torseur2D(object):
    def __init__(self, resul, mom_ref, ref):
        '''torseur de résultante `resul` et de moment `mom_ref`
        exprimé au point de coordonnées `ref`.
        '''
        # Résultante
        self.resul = resul
        # Moment exprimé en ref
        self.mom_ref = mom_ref
        # Point de référence
        self.ref = ref
        
        print("Initialisation d'un Torseur")
    # end __init__()

# Création ("Instanciation") d'un torseur
T1 = Torseur2D((2,3), 1., (0,0))
T1.mom_ref
T1.ref
T1.resul
Initialisation d'un Torseur
Out[9]:
(2, 3)

Remarquons que la Docstring de init sert à documnter la création de la classe :

In [ ]:
Torseur2D?

Toutes les attributs d'un objet sont accessibles en lecture et en écriture :

In [12]:
T1.resul = [4, 5]
T1.resul
Out[12]:
[4, 5]

Remarque : Par défaut, aucune vérification de la valeur assignée n'est effectuée, mais des vérifications sont facilement implémentables (notion de getter et setter)

Formule du "Changement de Point"

Une des formules essentielle pour travailler avec les torseurs est de calculer le moment en n'impore quel point ("formule de changement de point")

On ajoute donc une nouvelle méthode qui implémente de le changement de point :

$$\vec{M}(X) = \vec{M}(O) + \vec{R} \times \vec{OX} $$

avec $\vec{R}$ la résultante du torseur et $O$ le point de référence initial.

In [ ]:
# méthode à ajouter dans la définition de la classe Torseur2D
def calc_mom(self, pt):
    '''calcule le moment du torseur en un point `pt` quelconque'''
    # Formule du changement de point: M(X) = M(O) + R^OX
    # 1) Vecteur déplacement OX = X-O:
    delta_ref = (pt[0] - self.ref[0],
                 pt[1] - self.ref[1])
    # 2) Produit vectoriel R^OX
    delta_mom = self.resul[0]*delta_ref[1] - self.resul[1]*delta_ref[0]
    # 3) Nouveau moment:
    mom = self.mom_ref + delta_mom
    return mom
# end calc_mom()

Souvent, les méthodes servent à changer l'état de l'objet. Ici, on propose par exemple de changer le point de référence (ce qui recalcule le moment au point de référence)

In [ ]:
# méthode à ajouter dans la définition de la classe Torseur2D
def set_ref(self, new_ref):
    '''change le point de référence'''
    # 1) Calcul du moment au nouveau point de référence
    self.mom_ref = self.calc_mom(new_ref)
    # 2) Mise à jour du point de référence
    self.ref = new_ref
# end set_ref()

Méthodes suplémentaires

On peut obtenir une belle représentation textuelle en réimplémentant la méthode spéciale __str__ (celle-ci est appelée par str(mon_torseur))

In [ ]:
def __str__(self):
    res = 'Torseur 2D '
    res += 'de Résultante {0} et de Moment {1}'.format(str(self.resul),
                                                        str(self.mom_ref))
    res += ' exprimé au point {0}'.format(str(self.ref))
    
    return res.encode('utf-8')
# end __str__

Pour pouvoir additionner/soustraire des torseurs, il faut réimplémenter les méthodes correspondants à ces opérations :

  • __add__ pour a + b
  • __neg__ pour -a
  • __sub__ pour a-b

(cf. Emulating numeric types http://docs.python.org/2/reference/datamodel.html#emulating-numeric-types)

In [ ]:
### Opération arithmétiques
def __add__(self, other):
    '''addition de 2 torseurs: T1 + T2
    Le pt de référence utilisé est celui de T1 (opérande gauche)
    '''
    # 1) addition des Résultantes:
    resul = (self.resul[0] + other.resul[0],
             self.resul[1] + other.resul[1])
    # 2) addition des moments au point de référence de self:
    mom_ref = self.mom_ref
    other_mom = other.calc_mom(self.ref)
    mom_ref = mom_ref + other_mom
    # 3) Création de l'objet Torseur:
    T = Torseur2D(resul, mom_ref, self.ref)
    return T
# end __add__

def __neg__(self):
    '''opposé du torseur'''
    resul = (-self.resul[0], -self.resul[1])
    mom_ref = -self.mom_ref
    T = Torseur2D(resul, mom_ref, self.ref)
    return T
# end __neg__

def __sub__(self, other):
    '''soustraction de 2 torseurs: T1-T2'''
    return self + (-other)

Définition complète de la classe Torseur2D

In [30]:
class Torseur2D(object):
    def __init__(self, resul, mom_ref, ref):
        '''torseur de résultante `resul` et de moment `mom_ref`
        exprimé au point de coordonnées `ref`.
        '''
        # Résultante
        self.resul = resul
        # Moment exprimé en ref
        self.mom_ref = mom_ref
        # Point de référence
        self.ref = ref        
        print("[Initialisation d'un nouveau Torseur]")
        #print(self)
    # end __init__()
    
    def calc_mom(self, pt):
        '''calcule le moment du torseur en un point `pt` quelconque'''
        # Formule du changement de point: M(X) = M(O) + R^OX
        # 1) Vecteur déplacement OX = X-O:
        delta_ref = (pt[0] - self.ref[0],
                     pt[1] - self.ref[1])
        # 2) Produit vectoriel R^OX
        delta_mom = self.resul[0]*delta_ref[1] - self.resul[1]*delta_ref[0]
        # 3) Nouveau moment:
        mom = self.mom_ref + delta_mom
        return mom
    # end calc_mom()

    def set_ref(self, new_ref):
        '''change le point de référence'''
        # 1) Calcul du moment au nouveau point de référence
        self.mom_ref = self.calc_mom(new_ref)
        # 2) Mise à jour du point de référence
        self.ref = new_ref
    # end set_ref()

    def __str__(self):
        res = u'Torseur 2D '
        res += u'de Résultante {0} et de Moment {1}'.format(str(self.resul),
                                                            str(self.mom_ref))
        res += u' exprimé au point {0}'.format(str(self.ref))
        
        return res.encode('utf-8')
    # end __str__
    
    ### Opération arithmétiques
    def __add__(self, other):
        '''addition de 2 torseurs: T1 + T2
        Le pt de référence utilisé est celui de T1 (opérande gauche)
        '''
        # 1) addition des Résultantes:
        resul = (self.resul[0] + other.resul[0],
                 self.resul[1] + other.resul[1])
        # 2) addition des moments au point de référence de self:
        mom_ref = self.mom_ref
        other_mom = other.calc_mom(self.ref)
        mom_ref = mom_ref + other_mom
        # 3) Création de l'objet Torseur:
        T = Torseur2D(resul, mom_ref, self.ref)
        return T
    # end __add__
    
    def __neg__(self):
        '''opposé du torseur'''
        resul = (-self.resul[0], -self.resul[1])
        mom_ref = -self.mom_ref
        T = Torseur2D(resul, mom_ref, self.ref)
        return T
    # end __neg__
    
    def __sub__(self, other):
        '''soustraction de 2 torseurs: T1-T2'''
        return self + (-other)

Initialisation de 3 torseurs :

In [40]:
T1 = Torseur2D((1,1),  0, (0,0))
T2 = Torseur2D((1,1),  1, (0,0))
T3 = Torseur2D((1,1), -1, (0,0))
print(T3)
[Initialisation d'un nouveau Torseur]
[Initialisation d'un nouveau Torseur]
[Initialisation d'un nouveau Torseur]
Torseur 2D de Résultante (1, 1) et de Moment -1 exprimé au point (0, 0)

Quelques opérations de démonstration :

In [41]:
# Changer le point de référence :
T3.set_ref((-1,0))
print(T3)
Torseur 2D de Résultante (1, 1) et de Moment 0 exprimé au point (-1, 0)
In [42]:
print(-T1)
[Initialisation d'un nouveau Torseur]
Torseur 2D de Résultante (-1, -1) et de Moment 0 exprimé au point (0, 0)
In [43]:
print(T1+T2) # torseurs déjà exprimés au même pt
[Initialisation d'un nouveau Torseur]
Torseur 2D de Résultante (2, 2) et de Moment 1 exprimé au point (0, 0)
In [44]:
print(T1+T3) # torseurs exprimés en des pts différents
[Initialisation d'un nouveau Torseur]
Torseur 2D de Résultante (2, 2) et de Moment -1 exprimé au point (0, 0)

Pour compléter cette classe Torseur2D, on pourrait implémenter les opérations de comparaisons (== et !=) : http://docs.python.org/2/reference/datamodel.html#object.\__eq\__


B) Python numérique avec NumPy

B1) Pourquoi un module particulier (NumPy) ?

Pour le calcul scientifique, on est souvent amené à manipuler de grands ensembles de nombres. Les listes Python (list, vu à la Séance 1) sont des conteneurs hétérogènes ne sont pas adaptées pour cela car :

  • les opérations arithmétiques (addition, multiplication) ne sont pas réalisable facilement (car le + pour les listes est une concaténation)
  • une liste Python ne stocke par sont contenu de façon efficace (i.e. contiguë) en mémoire.

Pour manipuler des tableau de nombres, il y a le module NumPy (http://docs.scipy.org/doc/numpy/user/) qui crée un nouveau type de tableau array

In [33]:
# Importer le module numpy avec un nom raccourci :
import numpy as np

# Créer un "array" à partir d'une liste
l = np.array([10, 20, 30, 40]) # l n'est plus une liste mais un tableau NumPy

print(type(l))
l
<type 'numpy.ndarray'>
Out[33]:
array([10, 20, 30, 40])

À première vue, un tableau NumPy ressemble beaucoup à une liste. En particulier on peut l'indexer comme une liste :

In [34]:
l[0]
Out[34]:
10

Grosse différences, on réalise facilement des opérations arithmétiques terme à terme :

In [35]:
l+l
Out[35]:
array([20, 40, 60, 80])
In [36]:
l*10
Out[36]:
array([100, 200, 300, 400])
In [37]:
l**2
Out[37]:
array([ 100,  400,  900, 1600])

On peut également utiliser des fonctions mathématiques usuelles (racine, exponentielle, ...) mais qui fonctionnent de façon vectorisée :

In [9]:
np.exp(l)
np.sqrt(l)
Out[9]:
array([ 3.16227766,  4.47213595,  5.47722558,  6.32455532])

Liste des ~60 fonctions vectorisée : http://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs (appelées "ufuncs" dans le jargon NumPy)

B2) Fonctions de bases de NumPy

Création de tableaux

In [12]:
# On a déjà vu la fonction "constructeur" qui transforme tout type de listes (même imbriquées) en tableau
l = [1, 6, 10]
np.array(l)
Out[12]:
array([ 1,  6, 10])
In [13]:
# L'équivalent de `range` :
np.arange(10)
Out[13]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [14]:
print(np.zeros(5))
print(np.ones(5))
[ 0.  0.  0.  0.  0.]
[ 1.  1.  1.  1.  1.]

Tableaux à plusieurs dimensions

In [15]:
# avec des listes imbriquées
np.array([[1,2,3], [4,5,6]])
Out[15]:
array([[1, 2, 3],
       [4, 5, 6]])
In [40]:
# Matrice identité 3x3 :
np.eye(3)
Out[40]:
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

Question : comment créer (avec np.zeros) un tableau à plusieur dimensions ? Une matrice 4x3 par exemple ?

In [42]:
a = np.zeros( (4,3) )
# attention aux parenthèses pour que (3,3) soit bien *un seul argument*
a
Out[42]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])

L'information sur le nombre de dimension est accessible avec l'attribut .ndim.

In [44]:
print(l.ndim)
a.ndim
1
Out[44]:
2

Sa dimension (sa "forme", cad shape en anglais) est accessbile avec .shape (qui est un tuple de longueur ndim)

In [45]:
a.shape
Out[45]:
(4, 3)

Trancher dans un tableau ("slicing")

Un des grands intérêt des tableaux NumPy est la capacité d'accéder sélectivement à des morceaux de tableau. Pour cela la syntaxe d'indexation (avec [...]) est étendu pour permettre de selectionner chaque dimension

In [59]:
a = np.array([[1,2,3], [4,5,6], [7,8,9],[10,11,12]])
a
Out[59]:
array([[ 1,  2,  3],
       [ 4,  5,  6],
       [ 7,  8,  9],
       [10, 11, 12]])
In [60]:
# sélectionner 1 élément
a[0,0]
Out[60]:
1
In [61]:
# sélectionner la 1ère colonne 
a[:, 0]
Out[61]:
array([ 1,  4,  7, 10])

Pour sélection plusieurs lignes/colonnes

In [73]:
# Méthode 1 : indexer avec la syntaxe start:stop:step
a[:, 0:2]
Out[73]:
array([[ 1,  2],
       [ 4,  5],
       [ 7,  8],
       [10, 11]])
In [74]:
# Méthode 2 : indexer avec une liste (plus général)
a[:, [0, 1] ]
Out[74]:
array([[ 1,  2],
       [ 4,  5],
       [ 7,  8],
       [10, 11]])

Mini-question : comment sélectionner la 2ème ligne ?

In [63]:
# sélectionner la 2ème (-> index 1) ligne
a[1, :]
Out[63]:
array([4, 5, 6])

Indexation logique : pour sélectioner des éléments selon une condition. Cela se fait en 2 étapes :

In [66]:
# Étape 1 : créé le tableau booléen qui va sélectionner les éléments
a5 = a > 5
a5
Out[66]:
array([[False, False, False],
       [False, False,  True],
       [ True,  True,  True],
       [ True,  True,  True]], dtype=bool)
In [67]:
# Étape 2 : utiliser le tableau booléen pour indexer :
a[a5]
Out[67]:
array([ 6,  7,  8,  9, 10, 11, 12])
In [87]:
# Écriture compacte (1 et 2 en une fois)
a[a>5]
Out[87]:
array([ 6,  7,  8,  9, 10, 11, 12])

Typage de tableaux ("dtype")

Les tableaux NumPy contiennent des types de données homogènes. L'attribut .dtype permet de connaître ce type.

In [28]:
print(np.array([1, 2, 3]).dtype) # tableau d'entier
print(np.array([1., 2., 3.]).dtype) # tableau de flottants
int64
float64

Les fonctions de création de tableau ont également un argument dtype permettant d'indiquer le type souhaité

In [46]:
a = np.array([1, 2, 3])
# le tableau créé prend par défaut le type des donnnées d'entrée : ici `int`
a/2
Out[46]:
array([0, 1, 1])
In [28]:
a = np.array([1, 2, 3], dtype=float)
a/2
Out[28]:
array([ 0.5,  1. ,  1.5])

B3) Fonctions d'algèbre linéaire avec NumPy

Produit matriciel : attention !

Par défaut, le produit (*) de 2 tableaux numpy se fait terme à terme. Pour le produit matriciel, il y a np.dot

In [147]:
# Fabrication d'une matrice diagonale à partir d'une liste des éléments diagonaux
A = np.diag([1,2])
B = np.ones((2,2))*3
print(A)
print(B)
A * B
[[1 0]
 [0 2]]
[[ 3.  3.]
 [ 3.  3.]]
Out[147]:
array([[ 3.,  0.],
       [ 0.,  6.]])
In [148]:
np.dot(A,B)
Out[148]:
array([[ 3.,  3.],
       [ 6.,  6.]])

Remarque : dans la même lignée, on trouve le produit scalaire (i.e. produit interne np.inner) et le produit externe (np.outer)

In [150]:
np.inner?

Algèbre linéaire

Les fontions dites d'algèbre linéaire sont regroupés dans le sous-module numpy.linalg. On y trouve :

  • calcul de déterminant
  • inversion de matrice
  • calcul de valeurs propres
  • ...

Documentation : http://docs.scipy.org/doc/numpy/reference/routines.linalg.html

Certaines de ces fonctions seront mises en oeuvres lors de la Séance 3. On illustre ici simplement le calcul numérique d'un déterminant

In [143]:
# Fabrication d'une matrice diagonale à partir d'une liste des éléments diagonaux
a = np.diag([1,2,3])
a
Out[143]:
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
In [142]:
np.linalg.det(a) # det(a) vaut 1*2*3
Out[142]:
6.0

B4) Lecture/écriture de fichier de données

Pour les fichiers texte du type CSV (nombres séparés par des virgules, ou autres), NumPy fournit des fonctions de lecture.

La fonction np.loadtxt, assez flexible, permet de lire des fichiers avec différents formats. Cependant, cette flexibilité aboutit à une signature assez chargée

In [136]:
# Consultons l'aide:
np.loadtxt?
In [29]:
# Lecture du fichier:
conso = np.loadtxt('./S2_Objets_et_NumPy/RTE_conso_2010-08_2011-07.csv',
                   delimiter=',', # les champs sont séparés par des virgules
                   skiprows=4, # ne pas lire les 4 premières lignes (du blabla)
                   usecols=[1]) # ne lire que la colonne *2* (comptage à partir de 0 !)

N = conso.shape[0]
print(u'{} pts chargés, soit {:.1f} jours'.format(N, N/24./4))
35040 pts chargés, soit 365.0 jours
In [30]:
# Conversion MW -> GW
conso = conso/1000.

On peut faire quelques statistiques simples :

(en utilisant les méthodes de l'objet ndarray)

In [31]:
print(conso.min())
print(conso.max())
print(conso.mean())
print(conso.std())
30.993
96.35
56.1964912386
13.0847858651

Petite question : comment compter le nombre de jours passés au-dessus de 80 GW ?

In [78]:
# Méthode 1 : boucle for (très déconseillée pour les grands tableaux)

N80 = 0
for c in conso:
    if c >= 80:
        # on incrémente le compteur
        N80 += 1
N80
Out[78]:
1841

Solutions avec une opération vectorielle de comparaison

In [86]:
conso_80 = conso >= 80 # (GW)
conso_80
Out[86]:
array([False, False, False, ..., False, False, False], dtype=bool)
In [85]:
# Méthode 2 : indexaction logique pour ne sélectionner que les éléments vérifiant conso >= 80
len(conso[conso_80])
Out[85]:
1841
In [83]:
# Méthode 3 : jouer avec le transtypage implicite  (bool -> int)
# On calcule la "somme"
conso_80.sum()
Out[83]:
1841

Écriture de données CSV

Le pendant pour l'écriture de np.loadtxt est np.savetxt. On se réferera à la documentation.

In [33]:
np.savetxt?

C) Tracé de graphiques avec Matplotlib

La bibliothèque Matplotlib ( http://matplotlib.org ) permet de réaliser une grande palette de graphiques 2D (et même un peu 3D) avec une syntaxe fortement calquée sur celle de Matlab.

  • tracés temporels y = f(t)
  • nuages de points (x,y)
  • histogrammes
  • images
  • ...

La gallerie ( http://matplotlib.org/gallery.html ) est une bonne façon d'avoir un aperçu des possibilités de Matplotlib.

Source de documentation :

In [47]:
# Importer l'API de Matplotlib:
import matplotlib.pyplot as plt
# Remarques :
# * cet import a déjà été réalisé en mode "Pylab"
# * il est par contre nécessaire dans un script

# Pour un travail interactif, on s'autorisera aussi
from matplotlib.pyplot import *

C1) Aperçus de Matplotlib

Un exemple de graphique minimal

In [48]:
# Création des "données" : un sinus
x = np.linspace(0, 4*np.pi)
y = np.sin(x)

# Tracé
plt.plot(x,y) # ligne (bleue par défaut)
plt.plot(x,y, 'r+') # marqueurs "+" rouges

# Annotations:
plt.title('Un sinus')
plt.xlabel('temps t')
plt.grid(True);

Exemple pour le déplacement et le zoom interactif: Consommation d'électricité

Lorsque qu'on travaille avec un grand nombre de points (ex : grande série temporelle), l'exploration interactive des données est particulièrement utile.

Remarque : pour manipuler des grandes séries de données, il existe le puissant module pandas (Python Data Analysis Library http://pandas.pydata.org/) qui ne sera pas abordé ici.

In [51]:
# Vecteur temps en jours:
t = np.arange(N)/24./4 # attention au '24.'
t[:10]
Out[51]:
array([ 0.        ,  0.01041667,  0.02083333,  0.03125   ,  0.04166667,
        0.05208333,  0.0625    ,  0.07291667,  0.08333333,  0.09375   ])
In [52]:
# %pylab qt # afficher la figure dans une fenêtre interactive: nécessite IPython 1.0 (prévu pour juillet 2014)
plt.plot(t, conso)

plt.xlabel('temps (j)')
plt.ylabel(u'Puissance consommée (GW)');

Pour se déplacer et zoomer interactivement, utiliser la barre d'outils Matplotlib ou bien les raccourcis clavier ('p' et 'o').

(icônes et )

  • Les clic-gauche et clic-droit ont une action différente !
  • Les touches 'x' et 'y' modifient également le comportement (déplacement/zoom horizontal ou vertical)
In [53]:
# Zoom "programmé"
plt.plot(t, conso);
plt.title('Zoom sur 2 semaines consommation')
plt.xlabel('temps (j)')

plt.xlim(99, 113);
# il existe aussi plt.ylim
Out[53]:
(99, 113)

Petite question : sur ce graphique de 14 jours, quels jours correspondent aux weekends ?

Un histogramme

In [117]:
plt.hist(conso, bins=50, normed=False)
# (mettre normed=True pour tracer une densité de probabilité)
plt.xlabel(u'Puissance consommée (GW)');

C2) Exemple Numérique & Graphique plus complet : un diagramme de Bode

L'enjeu de cette partie :

  • charger des données "expérimentales" issue du relevé d'un diagramme de Bode
  • tracer un beau diagramme de Bode (gain et phase)
  • modéliser une réponse fréquentielle par une fonction de transfert (du 2ème ordre).
In [89]:
# Chargement de données  (pseudo relevé expérimental d'un diagramme de Bode)
data = np.loadtxt('./S2_Objets_et_NumPy/bode_data.csv', delimiter=',')
data.shape # nombre de lignes, nombre de colonnes
Out[89]:
(100, 3)

Le fichier bode_data.csv contient 3 colonnes:

  • fréquence (Hz)
  • gain (dB)
  • phase (°)

Question de slicing : comment séparer ces 3 colonnes dans 3 variables différentes ?

In [90]:
# Séparation des données
f = data[:,0]
gain = data[:,1]
phase = data[:,2]

# ou en abrégé :
# f, gain, phase = data.T
In [91]:
# Affichage (tracé en échelle linéaire):
plt.plot(f, gain);

Avec une échelle linéaire en fréquence, on ne reconnait pas grand chose...

Une solution: le graphique log/semilog avec les fonctions semilogx, semilogy et loglog. (noms de fonction Matlab)

In [92]:
# Tracé en échelle semilog pour les abscisses (la fréquence)
semilogx(f, gain)
# Annotations:
title('Diagramme de Bode')
xlabel(u'fréquence (Hz)')
ylabel('gain (dB)');

Tracé conjoint du gain et de la phase:

Pour ce graphique plus compliqué, on utilise l'approche Orienté Objet ("OO" API). (Accessoirement, cela dispense d'utiliser plt.semilogx)

In [133]:
fig = plt.figure('Bode', figsize=(6,6))
ax1 = fig.add_subplot(2,1,1, xscale='log',
                      title='Diagramme de Bode: gain (dB)')
ax1.plot(f, gain)
ax2 = fig.add_subplot(2,1,2, sharex=ax1,  # sharex -> synchronise les échelles. Pratique pour zoomer !
                      title=u'Phase (°)',
                      xlabel=u'fréquence (Hz)')
ax2.plot(f, phase, 'g')
# ajustement des "ticks"
ax2.set_yticks([-180, -135, -90, -45, 0])
# ou de façon plus automatique :
# ax2.set_yticks(range(0,-181,-45)) # (de 0 à -180°, tous les 45°)

fig.tight_layout() # "étire" les graphiques dans la fenêtre

On reconnait la fonction de tranfert d'un passe-bas du 2ème ordre

$$ H(f) = \frac{H_0}{1 + j \frac{f}{Q f_0} - \frac{f^2}{f_0^2}} $$

Les paramètres à identifier sont :

  • le gain statique $H_0$
  • la fréquence canonique $f_0$
  • le facteur de qualité $Q$

Identification des paramètres (exercice)

Cherchons le gain basse fréquence $H_0$ dans les 5 premieres valeurs de gain :

In [94]:
# Observation du gain à basse fréquence:
print(f[:5])
print(gain[:5])
gain[:5].mean()
[ 100.      104.7616  109.7499  114.9757  120.4504]
[ 6.155659  6.001112  5.810821  6.184953  6.278457]
Out[94]:
6.0862004000000001
In [95]:
# dB -> linéaire
H0_dB = 6 # dB
print(10**(H0_dB/20.)) # Attention au point de "20."
H0 = 2
1.99526231497

Cherchons la fréquence canonique avec le maximum du gain

In [96]:
gain.max()
Out[96]:
14.366379999999999
In [97]:
# Fréquence canonique: c'est environ la fréquence de résonance
i_max = gain.argmax()
f[i_max]
Out[97]:
977.00999999999999
In [98]:
f0 = 1e3 

Pour le facteur de qualité, on utilise la méthode devinette (cf. Séance 4 pour une identification par minimisation de l'erreur)

In [99]:
# Facteur de qualité : devinette
Q = 10

Calcul de la fonction de transfert du passe-bas

In [104]:
# On construit une fonction de transfert

# Dans cet exemple, on a déjà un vecteur de fréquence. À défaut, on pourrait utiliser :
# f = np.logspace(2, 4, 100)# de 100 à 10 kHz

# Tranfert du 2è ordre:
H_fit = H0/(1 + 2j*f/(f0*Q) - (f/f0)**2)
H_fit[:5] # (tableau de nombres complexes)
Out[104]:
array([ 2.01937787-0.04079551j,  2.02128641-0.0428206j ,
        2.02338500-0.04495474j,  2.02569298-0.04720512j,
        2.02823182-0.04957958j])
In [108]:
# Calcul du gain et de la phase
Habs = np.abs(H_fit) # `abs` donne le module du nombre complexe

gain_fit = 20*np.log10(Habs) # -> dB
gain_fit[:5]
Out[108]:
array([ 6.10612393,  6.1145058 ,  6.12371376,  6.13383021,  6.14494616])
In [109]:
phase_fit = np.angle(H_fit) * 180 / np.pi
phase_fit[:5]
Out[109]:
array([-1.15733307, -1.21361943, -1.27276485, -1.33493309, -1.40030107])

Comparaison données/modèles (et légende d'un graphique)

In [107]:
# Par défaut, les graphiques se superposent:
semilogx(f, gain, '-+', label=u'données')
semilogx(f, gain_fit, label=u'modèle')
legend(loc='lower left');

# Fabrication de la chaîne de caractère de titre:
title_str = u'Modélisation du 2è ordre: H0={:.1f}, f0={:.1f} kHz, Q={:.1f}'.format(H0, f0/1000, Q)
# La syntaxe du "string formatting" : http://docs.python.org/2/library/string.html#formatstrings
title(title_str);

Apparemment on a été un peu optimiste pour le facteur de qualité. $Q=5$ sera sans doute plus approprié... (cf. script Bode_data_generator.py qui a fabriqué les données)

La Séance 4 abordera l'identification des paramètres $H_0, Q, f_0$ par minimisation de l'erreur de gain et de phase.

Tracés bonus

  • comment tracer un diagramme de Black ? (phase, gain dB)
  • comment tracer un diagramme de Nyquist ? (Re, Im)
In [139]:
# Black :
plot(phase, gain)
plot(phase_fit, gain_fit)
plot(-180,0, 'r*', markersize=15) # point critique

# ajustement des "ticks"
xticks([-180, -135, -90, -45, 0])

xlabel(u'phase (°)')
ylabel('Gain (dB)');
In [134]:
# Diagramme de Nyquist

# Étape 1 : on reconstruit la réponse fréquentielle complexe
gain_lin = 10**(gain/20.) # dB -> lin
H = gain_lin*np.exp(1j*phase/180*np.pi)
In [137]:
# Étape 2 : tracé (utilisant les attributs .real et .image des tableaux de nombres complexes)
plot(H.real, H.imag)
plot(H_fit.real, H_fit.imag)
plot(-1,0, 'r*', markersize=10) # point critique

ax = gca()
ax.set_aspect('equal') # Même échelle en x et y
ylim(-11,1);

Fin de la session 2

Licence Creative Commons

Ce <span xmlns:dct="http://purl.org/dc/terms/" href="http://purl.org/dc/dcmitype/InteractiveResource" rel="dct:type">support de formation</span> créé par <span xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">Pierre Haessig</span> est mis à disposition selon les termes de la licence Creative Commons Attribution 3.0 France.