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 :
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
(2, 3)
Remarquons que la Docstring de init sert à documnter la création de la classe :
Torseur2D?
Toutes les attributs d'un objet sont accessibles en lecture et en écriture :
T1.resul = [4, 5]
T1.resul
[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.
# 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)
# 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))
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)
### 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¶
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 :
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 :
# 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)
print(-T1)
[Initialisation d'un nouveau Torseur] Torseur 2D de Résultante (-1, -1) et de Moment 0 exprimé au point (0, 0)
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)
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.%5C__eq%5C__
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
# 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'>
array([10, 20, 30, 40])
À première vue, un tableau NumPy ressemble beaucoup à une liste. En particulier on peut l'indexer comme une liste :
l[0]
10
Grosse différences, on réalise facilement des opérations arithmétiques terme à terme :
l+l
array([20, 40, 60, 80])
l*10
array([100, 200, 300, 400])
l**2
array([ 100, 400, 900, 1600])
On peut également utiliser des fonctions mathématiques usuelles (racine, exponentielle, ...) mais qui fonctionnent de façon vectorisée :
np.exp(l)
np.sqrt(l)
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¶
- documentation de référence : http://docs.scipy.org/doc/ (en anglais)
- Support de formation "SciPy lectures" : http://scipy-lectures.github.io/intro/numpy/array_object.html (en anglais)
- Tutoriel de Nicolas Rougier http://www.loria.fr/~rougier/teaching/numpy/numpy.html (en anglais)
- Autre support en français : http://python-prepa.github.io/edp_chaleur.html#introduction-rapide-a-numpy (formation Python ENS Paris)
Création de tableaux¶
# 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)
array([ 1, 6, 10])
# L'équivalent de `range` :
np.arange(10)
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
print(np.zeros(5))
print(np.ones(5))
[ 0. 0. 0. 0. 0.] [ 1. 1. 1. 1. 1.]
Tableaux à plusieurs dimensions¶
# avec des listes imbriquées
np.array([[1,2,3], [4,5,6]])
array([[1, 2, 3],
[4, 5, 6]])
# Matrice identité 3x3 :
np.eye(3)
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 ?
a = np.zeros( (4,3) )
# attention aux parenthèses pour que (3,3) soit bien *un seul argument*
a
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.
print(l.ndim)
a.ndim
1
2
Sa dimension (sa "forme", cad shape en anglais) est accessbile avec .shape (qui est un tuple de longueur ndim)
a.shape
(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
a = np.array([[1,2,3], [4,5,6], [7,8,9],[10,11,12]])
a
array([[ 1, 2, 3],
[ 4, 5, 6],
[ 7, 8, 9],
[10, 11, 12]])
# sélectionner 1 élément
a[0,0]
1
# sélectionner la 1ère colonne
a[:, 0]
array([ 1, 4, 7, 10])
Pour sélection plusieurs lignes/colonnes
# Méthode 1 : indexer avec la syntaxe start:stop:step
a[:, 0:2]
array([[ 1, 2],
[ 4, 5],
[ 7, 8],
[10, 11]])
# Méthode 2 : indexer avec une liste (plus général)
a[:, [0, 1] ]
array([[ 1, 2],
[ 4, 5],
[ 7, 8],
[10, 11]])
Mini-question : comment sélectionner la 2ème ligne ?
# sélectionner la 2ème (-> index 1) ligne
a[1, :]
array([4, 5, 6])
Indexation logique : pour sélectioner des éléments selon une condition. Cela se fait en 2 étapes :
# Étape 1 : créé le tableau booléen qui va sélectionner les éléments
a5 = a > 5
a5
array([[False, False, False],
[False, False, True],
[ True, True, True],
[ True, True, True]], dtype=bool)
# Étape 2 : utiliser le tableau booléen pour indexer :
a[a5]
array([ 6, 7, 8, 9, 10, 11, 12])
# Écriture compacte (1 et 2 en une fois)
a[a>5]
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.
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é
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
array([0, 1, 1])
a = np.array([1, 2, 3], dtype=float)
a/2
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
# 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.]]
array([[ 3., 0.],
[ 0., 6.]])
np.dot(A,B)
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)
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
# Fabrication d'une matrice diagonale à partir d'une liste des éléments diagonaux
a = np.diag([1,2,3])
a
array([[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
np.linalg.det(a) # det(a) vaut 1*2*3
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
# Consultons l'aide:
np.loadtxt?
# 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
# Conversion MW -> GW
conso = conso/1000.
On peut faire quelques statistiques simples :
(en utilisant les méthodes de l'objet ndarray)
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 ?
# 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
1841
Solutions avec une opération vectorielle de comparaison
conso_80 = conso >= 80 # (GW)
conso_80
array([False, False, False, ..., False, False, False], dtype=bool)
# Méthode 2 : indexaction logique pour ne sélectionner que les éléments vérifiant conso >= 80
len(conso[conso_80])
1841
# Méthode 3 : jouer avec le transtypage implicite (bool -> int)
# On calcule la "somme"
conso_80.sum()
1841
Écriture de données CSV¶
Le pendant pour l'écriture de np.loadtxt est np.savetxt. On se réferera à la documentation.
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 :
- Documentation officielle http://matplotlib.org/contents.html (en anglais)
- Tutoriel de Nicolas Rougier (INRIA) http://webloria.loria.fr/~rougier/teaching/matplotlib/ (en anglais)
# 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¶
# 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.
# Vecteur temps en jours:
t = np.arange(N)/24./4 # attention au '24.'
t[:10]
array([ 0. , 0.01041667, 0.02083333, 0.03125 , 0.04166667,
0.05208333, 0.0625 , 0.07291667, 0.08333333, 0.09375 ])
# %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)
# 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
(99, 113)
Petite question : sur ce graphique de 14 jours, quels jours correspondent aux weekends ?
Un histogramme¶
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).
# 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
(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 ?
# Séparation des données
f = data[:,0]
gain = data[:,1]
phase = data[:,2]
# ou en abrégé :
# f, gain, phase = data.T
# 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)
# 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)
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 :
# 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]
6.0862004000000001
# 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
gain.max()
14.366379999999999
# Fréquence canonique: c'est environ la fréquence de résonance
i_max = gain.argmax()
f[i_max]
977.00999999999999
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)
# Facteur de qualité : devinette
Q = 10
Calcul de la fonction de transfert du passe-bas¶
# 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)
array([ 2.01937787-0.04079551j, 2.02128641-0.0428206j ,
2.02338500-0.04495474j, 2.02569298-0.04720512j,
2.02823182-0.04957958j])
# 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]
array([ 6.10612393, 6.1145058 , 6.12371376, 6.13383021, 6.14494616])
phase_fit = np.angle(H_fit) * 180 / np.pi
phase_fit[:5]
array([-1.15733307, -1.21361943, -1.27276485, -1.33493309, -1.40030107])
Comparaison données/modèles (et légende d'un graphique)
# 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)
# 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)');
# 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)
# É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
Ce support de formation créé par Pierre Haessig est mis à disposition selon les termes de la licence Creative Commons Attribution 3.0 France.