Séance 4 : Autres outils de calcul numérique

Formation Python du département Mécatronique

15h45 - 16h45 (1h)

Introduction de quelques aspects de la boîte à outils SciPy.

A) Simulation de système dynamique (i.e. Résolution d'Équations Différentielles Ordinaires)

  • Résolution itérative d'un système dynamique simple par méthode d'Euler.
  • Problèmes de stabilité (choix du pas de temps).
  • Comparaison avec les méthodes ODEs de SciPy (dans scipy.integrate).

B) Optimisation numérique

  • fonction fmin de scipy.optimize
  • application à l'identification de paramètres d'une fonction de transfert

A) Simulation de système dynamique

(c'est à dire : résolution numérique d'équation différentielle ordinaire ("ODE" en anglais))

In [1]:
import scipy.integrate

A1) Un exemple simple (1er ordre)

On considère une ED très simple (linéaire, du 1er ordre) qui correspond au modèle des systèmes physique suivant :

  • charge d'un circuit RC (en électronique)
  • montée en température d'un système thermique (circuit RC thermique équivalent)
  • l'évolution de la vitesse d'une bille tombant dans un fluide visqueux.
$$\tau \frac{dx}{dt}(t) = u - x(t), \quad x(0)=0$$

$\tau$ est la constante de temps (RC, ...) et $u$ est "l'entrée".

In [2]:
u = 2.
tau = 10.

# Discrétisation:
N = 101
T = 50
t = np.linspace(0, T, N)
dt = t[1] - t[0]
print('dt = {:.2f}'.format(dt))
dt/tau
dt = 0.50
Out[2]:
0.050000000000000003

On connait la solution théorique $x(t) = u(1-e^{-t/\tau})$.

Mini-exercice :

  • Calculer la solution théorique à partir du vecteur t et de np.exp
  • puis tracer cette solution
In [98]:
# Calcul analytique :
x_th = ...

# Tracé
plot(... );
title(u'solution théorique');

A2) Résolution numérique "à la main"

La résolution numérique par la méthode d'Euler (explicite)

$$ x(t+\Delta_t) = x(t) + \Delta_t \dot{x}(t) $$

Pour notre système du 1er ordre, on obtient une équation aux différences :

$$ x(t+\Delta_t) = a x(t) + b u$$

avec

  • $a = (1- \Delta_t/\tau)$ et
  • $b=\Delta_t/\tau$

Exercice : implémenter cette méthode d'Euler (à l'aide d'une boucle for)

In [99]:
x_eul = np.zeros(N)
for i in range(...):
    x_eul[...] = ...

plot(t, x_eul, '-+', label='Euler explicite')
plot(t, x_th, label='sol. analytique')
title('Comparaison des solutions')
legend(loc='lower right');

À présent, il s'agit de refaire le calcul en augmentant le pas de temps (cad diminuer $N$)

Question : à partir de quelle valeur de $\Delta_t$ la solution numérique devient elle

  • "mauvaise", "imprécise" ?
  • oscillante ?
  • instable ?

(cf. la valeur du coefficient $a$...)

A3) Utilisation d'un solver ODE de SciPy

Lorsqu'on a à résoudre une ED, on a intérêt à faire appel à des routines puissantes et éprouvées, plutôt que de faire un "Runge-Kutta maison"...

cf http://docs.scipy.org/doc/scipy/reference/integrate.html pour la documentation des fonctions Integration and ODEs de scipy.integrate

In [83]:
odeint = scipy.integrate.odeint
odeint?

On remarque que solveur odeint (et pas que lui !) a besoin de la fonction dérivée $f(y,t) \mapsto \dot{y}$

In [3]:
def deriv(y,t):
    '''dérivée calculé en y et au temps t'''
    return ...

# test de la fonction :
deriv(0,0)
Out[3]:
0.2
In [100]:
x_ode = odeint(deriv, 0, t)

plot(t, x_eul, '--+', label='Euler explicite')
plot(t, x_th, label='sol. analytique')
plot(t, x_ode, '--x', label='odeint')
title('Comparaison des solutions')
legend(loc='lower right');

Bien sûr, il serait alors intéressant de creuser et jouer sur les paramètres de la fonction boîte noire odeint.

En particulier, la comparaison brute Euler vs. odeint n'est pas "juste" car odeint s'autorise des points d'intégration intermédiaires...


B) Optimisation numérique

L'optimisation d'une fonction $f(x)$ consiste à trouver la valeur de l'argument $x \in \mathbb{R}^n$ qui minimise/maximise cette fonction.

L'optimisation peut servir à modéliser (fitter) des données expérimentales comme on va le voir ici.

In [4]:
import scipy.optimize
from scipy.optimize import fmin

La fonction fmin est une des méthodes disponibles dans scipy.optimize. Les autres sont présentés dans la documentation http://docs.scipy.org/doc/scipy/reference/optimize.html

(on notera par exemple qu'il existe des méthodes spécialisées pour le cas où l'argument $x$ est un scalaire. Ça ne sera pas le cas ici.)

In [5]:
fmin?

On reprend les données du diagramme de Bode de la Séance 2

In [7]:
# Chargement de données  (pseudo relevé expérimental d'un diagramme de Bode)
data = np.loadtxt(u'./S2_Objets_et_NumPy/bode_data.csv', delimiter=',')
# séparation des données :
f, gain, phase = data.T
# Tracé rapide:
semilogx(f, gain);

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$

À la différence de la session 2, on souhaite identifier ces paramètres automatiquement. ("model fitting")

In [12]:
# Tranfert du 2è ordre:
def bode_pb2(f, H0, f0, Q):
    '''Bode d'un passe bas du 2ème ordre
    renvoie gain (dB), phase (°) 
    '''
    H = H0/(1 + 2j*f/(f0*Q) - (f/f0)**2)
    Habs = np.abs(H)
    gain = 20*np.log10(Habs) # -> dB
    phase = np.angle(H)*180/np.pi
    return (gain, phase) # ceci est un 'tuple' 

# tracé rapide:
semilogx(f, bode_pb2(f,1,1e3, 10)[0], 'g-');
In [17]:
def bode_err((H0, f0, Q)):
    '''erreur quadratique d'ajustement du modèle'''
    gain_fit, phase_fit = bode_pb2(f, H0, f0, Q)
    err_gain = ...
    err_phase = ....
    err = np.sum(... ) + np.sum(... )
    return err
bode_err((2, 1e3, 10))
Out[17]:
7154.6618314212756
In [21]:
x_opt = fmin(bode_err, x0=(1, 1e3, 10))
x_opt
Optimization terminated successfully.
         Current function value: 585.216147
         Iterations: 114
         Function evaluations: 205
Out[21]:
array([   1.9715475 ,  999.80793731,    4.98163544])
In [25]:
# Dépaquetage des paramètres :
H0, f0, Q = x_opt
print('H0 = {:.2f}'.format(H0))
print('f0 = {:.0f} kHz'.format(f0))
print(' Q = {:.2f}'.format(Q))
H0 = 1.97
f0 = 1000 kHz
 Q = 4.98

Remarque : dans ce problème d'estimation de paramètres, on pourrait chercher également l'incertitude entourant ces paramètres (incertitude due au erreurs de mesures).

In [104]:
# diagramme de Bode complet
fig = plt.figure('Bode', figsize=(6,6))

# gain
ax1 = fig.add_subplot(2,1,1, xscale='log',
                      title='Diagramme de Bode: gain (dB)')
ax1.plot(f, gain)
ax1.plot(f, bode_pb2(f, H0, f0, Q)[0], 'r-')
# phase
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, 'b')
ax2.plot(f, bode_pb2(f, H0, f0, Q)[1], 'r-')
# ajustement des "ticks"
ax2.set_yticks(np.arange(0,-181,-45)) # (de 0 à -180°, tous les 45°)
fig.tight_layout()

fin de la session 4 et fin de cette formation.

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.