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)
scipy.integrate
).B) Optimisation numérique
fmin
de scipy.optimize
(c'est à dire : résolution numérique d'équation différentielle ordinaire ("ODE" en anglais))
import scipy.integrate
On considère une ED très simple (linéaire, du 1er ordre) qui correspond au modèle des systèmes physique suivant :
$\tau$ est la constante de temps (RC, ...) et $u$ est "l'entrée".
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
0.050000000000000003
On connait la solution théorique $x(t) = u(1-e^{-t/\tau})$.
Mini-exercice :
t
et de np.exp
# Calcul analytique :
x_th = u*(1-np.exp(-t/tau))
# Tracé
plot(t, x_th);
title(u'solution théorique');
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
Exercice : implémenter cette méthode d'Euler (à l'aide d'une boucle for
)
x_eul = np.zeros(N)
for i in range(N-1):
x_eul[i+1] = (1-dt/tau)*x_eul[i] + dt/tau*u
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
(cf. la valeur du coefficient $a$...)
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
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}$
def deriv(y,t):
'''dérivée calculé en y et au temps t'''
return 1./tau*(u-y)
# test de la fonction :
deriv(0,0)
0.2
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...
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.
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.)
fmin?
On reprend les données du diagramme de Bode de la Séance 2
# 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 :
À la différence de la session 2, on souhaite identifier ces paramètres automatiquement. ("model fitting")
# 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-');
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 = (gain - gain_fit)
err_phase = (phase - phase_fit)
err = np.sum(err_gain**2) + np.sum(err_phase**2)
return err
bode_err((2, 1e3, 10))
7154.6618314212756
x_opt = fmin(bode_err, x0=(1, 1e3, 10))
x_opt
Optimization terminated successfully. Current function value: 585.216147 Iterations: 114 Function evaluations: 205
array([ 1.9715475 , 999.80793731, 4.98163544])
# 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).
# 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.
Ce support de formation créé par Pierre Haessig est mis à disposition selon les termes de la licence Creative Commons Attribution 3.0 France.