[PYTHON] Analyse numérique des équations différentielles ordinaires avec l'odeint et l'ode de Scipy

Une méthode de simulation (analyse / calcul numérique) d'un modèle physique sous la forme d'une équation différentielle normale telle qu'un simulateur de vol utilisant python. Ce que je fais, c'est utiliser odeint dans Scipy.integrate. Il semble qu'il utilise le lsode d'Odepack de Fortran, donc le calcul est rapide.

Exemple

À titre d'exemple, portons la simulation de vol de l'avion publié sous forme de code Octave (compatible Matlab) vers python. Le lien est détaillé pour le contenu du calcul. Butterfly_Effect( ) 6DOF Flight Simulation

Si vous utilisez for pour simuler quelque chose qui bouge dans le temps, même une simple méthode de calcul numérique telle que la méthode Euler est assez lente en python, vous devez donc utiliser autant que possible une bibliothèque telle que scipy.

Différence entre matlab et python

Notez que lors du portage depuis matlab, l'opérateur * de la matrice (ndarray) est différent de matlab: inner product et numpy: element product. Notez également qu'en python, si le point décimal n'est pas spécifié, une division entière se produira.

Code --odeint

# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import odeint
from scipy.linalg import block_diag
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.close('all')

# import time
# start = time.time()

def Rotation_X(Psi):
	R_x = np.array([[np.cos(Psi), np.sin(Psi), 0],
					[-np.sin(Psi), np.cos(Psi), 0],
					[0, 0, 1]])
	return R_x

def Rotation_Y(Theta):
	R_y = np.array([[np.cos(Theta), 0, -np.sin(Theta)],
					[0, 1, 0],
					[np.sin(Theta), 0, np.cos(Theta)]])
	return R_y

def Rotation_Z(Phi):
	R_z = np.array([[1, 0, 0],
					[0, np.cos(Phi), np.sin(Phi)],
					[0, -np.sin(Phi), np.cos(Phi)]])
	return R_z

#Équation du mouvement
def dynamical_system(x, t, A, U0):
	# x = [u,alpha,q,theta,beta,p,r,phi,psi,x,y,z]
	dx = A.dot(x)
	u = x[0]+U0 #la vitesse
	UVW = np.array([u, u*x[4], u*x[1]]) #Vecteur de vitesse[U,V,W]
	Rotation = Rotation_X(-x[8]).dot(Rotation_Y(-x[3])).dot(Rotation_Z(-x[7]))
	dX = Rotation.dot(UVW)
	dx[9]  = dX[0]
	dx[10] = dX[1]
	dx[11] = dX[2] 
	return dx

#Coefficient fin de stabilité dimensionnelle
Xu = -0.01;   Zu = -0.1;   Mu = 0.001
Xa = 30.0;    Za = -200.0; Ma = -4.0
Xq = 0.3;     Zq = -5.0;   Mq = -1.0
Yb = -45.0;   Lb_= -2.0;   Nb_= 1.0
Yp = 0.5;     Lp_= -1.0;   Np_= -0.1
Yr = 3.0;     Lr_= 0.2;    Nr_=-0.2

#Autres paramètres
W0 = 0.0;     U0 = 100.0;  theta0 = 0.05
g  = 9.8 #Accélération gravitationnelle

#Système vertical
A_lat = np.array([[Xu, Xa, -W0, -g*np.cos(theta0)],
    			  [Zu/U0, Za/U0, (U0+Zq)/U0, -g*np.sin(theta0)/U0],
    			  [Mu, Ma, Mq, 0],
    			  [0, 0, 1, 0]])

#Système latéral / directionnel
A_lon = np.array([[Yb, (W0+Yp), -(U0-Yr), g*np.cos(theta0), 0],
        		  [Lb_, Lp_, Lr_, 0, 0],
        		  [Nb_, Np_, Nr_, 0, 0],
        		  [0, 1, np.tan(theta0), 0, 0],
        		  [0, 0, 1/np.cos(theta0), 0, 0]])

#Combinez les systèmes sous forme de blocs diagonaux
A = block_diag(A_lat, A_lon)

#De plus, un espace sécurisé pour la trajectoire de vol
A = block_diag(A, np.zeros([3,3]))

#Définition des conditions de calcul
endurance	= 200	#Temps de vol[sec]
step		= 10	# 1.0[sec]Nombre de pas de temps par
t = np.linspace(0,endurance,endurance*step)

#Valeur initiale x0= [u,alpha,q,theta, beta,p,r,phi,psi]
#Les variables à inclure dans l'équation différentielle normale doivent être unidimensionnelles.
x0_lat		= np.array([10, 0.1, 0.4, 0.2])		#Valeur initiale verticale
x0_lon		= np.array([0.0, 0.6, 0.4, 0.2, 0.2])	#Valeur initiale horizontale / direction
x0_pos		= np.array([0, 0, -1000]) #Position initiale de l'avion
x0 = np.hstack((x0_lat, x0_lon, x0_pos))

x = odeint(dynamical_system, x0, t, args=(A,U0,))
print u"run successfully."

#Visualisation
plt.ion()
fig = plt.figure()
ax = Axes3D(fig)

ax.plot(x[:,9], x[:,10], x[:,11])
ax.set_xlabel(u"x");ax.set_ylabel(u"y");ax.set_zlabel(u"z")
ax.set_xlim([-5000,2000])
ax.set_ylim([-2000,5000])
ax.set_zlim([-5000,2000])
plt.show()

Code --ode

Dans scipy.integrate, il n'y a pas seulement odeint mais aussi une interface polyvalente pour le calcul numérique des équations différentielles ordinaires appelée ode, qui est orientée objet. Contrairement à odeint, vous pouvez spécifier la méthode de calcul, donc si vous voulez décider du contenu du calcul, c'est bien. De plus, il semble possible de terminer le calcul au milieu. Cependant, comme le calcul est effectué pendant, la vitesse de calcul est assez lente, en fonction du nombre de pas. C'était environ 10 fois plus lent dans mon environnement. Notez également que l'ordre des arguments de la fonction de l'équation différentielle est inversé par rapport à celui de (temps, variable) et odeint, et que ce doit être l'argument du tapple, il est donc légèrement différent.

# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import ode
from scipy.linalg import block_diag
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.close('all')

def Rotation_X(Psi):
	R_x = np.array([[np.cos(Psi), np.sin(Psi), 0],
					[-np.sin(Psi), np.cos(Psi), 0],
					[0, 0, 1]])
	return R_x

def Rotation_Y(Theta):
	R_y = np.array([[np.cos(Theta), 0, -np.sin(Theta)],
					[0, 1, 0],
					[np.sin(Theta), 0, np.cos(Theta)]])
	return R_y

def Rotation_Z(Phi):
	R_z = np.array([[1, 0, 0],
					[0, np.cos(Phi), np.sin(Phi)],
					[0, -np.sin(Phi), np.cos(Phi)]])
	return R_z

#Équation du mouvement
#Contrairement au cas d'odeint(time, x, args*)Notez que les arguments sont tapuleux
def dynamical_system(t, x, (A, U0)):
	# x = [u,alpha,q,theta,beta,p,r,phi,psi,x,y,z]
	dx = A.dot(x)
	u = x[0]+U0 #la vitesse
	UVW = np.array([u, u*x[4], u*x[1]]) #Vecteur de vitesse[U,V,W]
	Rotation = Rotation_X(-x[8]).dot(Rotation_Y(-x[3])).dot(Rotation_Z(-x[7]))
	dX = Rotation.dot(UVW)
	dx[9]  = dX[0]
	dx[10] = dX[1]
	dx[11] = dX[2] 
	return dx

#Coefficient fin de stabilité dimensionnelle
Xu = -0.01;   Zu = -0.1;   Mu = 0.001
Xa = 30.0;    Za = -200.0; Ma = -4.0
Xq = 0.3;     Zq = -5.0;   Mq = -1.0
Yb = -45.0;   Lb_= -2.0;   Nb_= 1.0
Yp = 0.5;     Lp_= -1.0;   Np_= -0.1
Yr = 3.0;     Lr_= 0.2;    Nr_=-0.2

#Autres paramètres
W0 = 0.0;     U0 = 100.0;  theta0 = 0.05
g  = 9.8 #Accélération gravitationnelle

#Système vertical
A_lat = np.array([[Xu, Xa, -W0, -g*np.cos(theta0)],
    			  [Zu/U0, Za/U0, (U0+Zq)/U0, -g*np.sin(theta0)/U0],
    			  [Mu, Ma, Mq, 0],
    			  [0, 0, 1, 0]])

#Système latéral / directionnel
A_lon = np.array([[Yb, (W0+Yp), -(U0-Yr), g*np.cos(theta0), 0],
        		  [Lb_, Lp_, Lr_, 0, 0],
        		  [Nb_, Np_, Nr_, 0, 0],
        		  [0, 1, np.tan(theta0), 0, 0],
        		  [0, 0, 1/np.cos(theta0), 0, 0]])

#Combinez les systèmes sous forme de blocs diagonaux
A = block_diag(A_lat, A_lon)

#De plus, un espace sécurisé pour la trajectoire de vol
A = block_diag(A, np.zeros([3,3]))

#Définition des conditions de calcul
endurance	= 200	#Temps de vol[sec]
step		= 10	# 1.0[sec]Nombre de pas de temps par
dt 			= 1.0 / step 	#Largeur du pas de temps[sec]
t0 			= 0.0	#Heure de début[sec]
t = np.linspace(0,endurance,endurance*step)

#Valeur initiale x0= [u,alpha,q,theta, beta,p,r,phi,psi,x,y,z]
#Les variables à inclure dans l'équation différentielle normale doivent être unidimensionnelles.
x0_lat		= np.array([10, 0.1, 0.4, 0.2])		#Valeur initiale verticale
x0_lon		= np.array([0.0, 0.6, 0.4, 0.2, 0.2])	#Valeur initiale horizontale / direction
x0_pos		= np.array([0, 0, -1000]) #Position initiale de l'avion
x0 = np.hstack((x0_lat, x0_lon, x0_pos))

#Paramètres ODE
solver = ode(dynamical_system).set_integrator('vode', method='bdf')
solver.set_initial_value(x0, t0)
solver.set_f_params((A, U0))
dt = 1.0 / step
x = np.zeros([int(endurance/dt)+1, 12])
t = np.zeros([int(endurance/dt)+1])

index = 0
while solver.successful() and solver.t < endurance:
	solver.integrate(solver.t+dt)
	x[index] = solver.y
	t[index] = solver.t
	index += 1

print u"run successfully."

#Visualisation
plt.ion()
fig = plt.figure()
ax = Axes3D(fig)

ax.plot(x[:,9], x[:,10], x[:,11])
ax.set_xlim([-5000,2000])
ax.set_ylim([-2000,5000])
ax.set_zlim([-5000,2000])
plt.show()