Je ne pouvais penser que la transformée de Fourier ou la transformée en ondelettes comme une méthode d'examen des changements de temps, mais j'ai appris qu'il existe une bonne méthode appelée décomposition en mode dynamique (DMD) qui peut extraire des modes dans les directions temporelle et spatiale. Je vais l'enregistrer ici pour approfondir ma compréhension.
Le contenu est presque le même que le site auquel j'ai fait référence, et j'écrirai une version légèrement modifiée de la traduction Google. Le code lié au dessin graphique n'est pas dans la source de référence, je l'ai donc ajouté. Vous devriez pouvoir tout faire, de l'exécution de DMD au dessin du graphique.
Pour simplifier, nous omettons le DMD du champ vectoriel 3D et ne considérerons qu'une simple fonction scalaire 1D. Dynamic Mode Decomposition in Python
Je ne savais pas ce qu'était SVD, alors j'y ai fait référence. À propos de la relation entre PCA et SVD
La décomposition en mode dynamique (DMD) est une innovation mathématique relativement récente qui peut résoudre ou approximer des systèmes dynamiques, entre autres, en ce qui concerne des structures cohérentes qui se développent, se désintègrent et / ou vibrent au fil du temps. La structure cohérente est appelée mode DMD. Chaque mode DMD a une dynamique temporelle correspondante définie pour une seule valeur unique.
En d'autres termes, la DMD transforme un système dynamique en une superposition de modes dont la dynamique est dominée par des valeurs propres.
De manière surprenante, la procédure mathématique pour distinguer les modes DMD des valeurs propres est purement linéaire, mais le système lui-même peut être non linéaire. Bien que cela ne soit pas couvert ici, l'affirmation selon laquelle les systèmes non linéaires peuvent être décrits comme un ensemble de paires de modes et de valeurs propres est justifiée. Pour plus d'informations, lisez l'article sur l'intégration de l'opérateur Koopman et DMD [^ 1] [^ 2] [^ 3].
DMD n'est pas seulement un outil de diagnostic utile pour analyser le comportement interne d'un système, mais il peut également être utilisé pour prédire l'état futur du système. Tout ce dont vous avez besoin est le mode et les valeurs uniques. Avec peu d'effort, vous pouvez combiner des modes et des valeurs propres pour générer une fonction qui se rapproche de l'état du système à tout moment.
Considérez n ensembles de données.
Où $ x_i $ et $ y_i $ sont des vecteurs colonnes de magnitude $ m $, respectivement. Définissez deux matrices $ m \ fois n
Si vous définissez l'opérateur $ A $ comme suit
Où $ X ^ \ dagger $ est une matrice pseudo-inverse de $ X $, et la décomposition en mode dynamique de $ (X, Y) $ est donnée par la décomposition en valeur propre de $ A $. Ainsi, le mode DMD et les valeurs propres sont les vecteurs propres et les valeurs propres de $ A $.
La définition de Tu et al. [^ 2] ci-dessus est connue sous le nom de DMD exact. Il s'agit actuellement de la définition la plus courante et peut être appliquée à tout ensemble de données répondant à certaines exigences. Dans cet article, je suis plus intéressé si $ A $ satisfait la formule $ y_i = Ax_i $ (probablement à peu près) pour tout $ i
Évidemment, $ X $ est un ensemble de vecteurs d'entrée et $ Y $ est un ensemble de vecteurs de sortie correspondants. Cette interprétation particulière de la DMD est très puissante car elle fournit un moyen pratique d'analyser (et de prédire) des tiges dynamiques avec des équations gouvernantes inconnues. Nous parlerons des systèmes dynamiques plus tard.
Il existe plusieurs théorèmes le long de cette définition de DMD [^ 2]. Un des théorèmes les plus utiles est que $ Yv = 0 chaque fois que $ X $ et $ Y $ sont linéairement cohérents (en d'autres termes, quand $ Xv = 0 $ pour le vecteur $ v $). Seulement si $ est également satisfait), $ Y = AX $ est entièrement satisfait. Tester la cohérence linéaire est relativement simple, comme nous le verrons plus tard. En bref, la cohérence linéaire n'est pas une condition préalable obligatoire pour l'utilisation de DMD. Même si la décomposition DMD de $ A $ ne satisfait pas complètement la formule $ Y = AX $, c'est la méthode des moindres carrés et minimise l'erreur dans la norme $ L ^ 2 $.
À première vue, la décomposition en valeur propre de $ A = YX ^ \ dagger $ ne semble pas être un gros problème. En fait, si les tailles $ X $ et $ Y $ sont correctes, appeler les méthodes pinv et eig depuis Numpy ou MATLAB plusieurs fois fonctionnera. Le problème se pose lorsque $ A $ est vraiment gros. Puisque $ A $ est $ m \ fois m $, la décomposition des valeurs propres devient lourde si $ m $ (le nombre de signaux dans chaque échantillon de temps) est très grand.
Heureusement, avec l'aide de l'algorithme DMD exact, vous pouvez diviser le problème en plus petits morceaux.
Calculez SVD (décomposition en valeurs singulières) de $ X $ et effectuez en même temps une troncature de bas niveau si nécessaire:
Projetez la matrice $ A $ sur $ U $ pour calculer $ \ tilde A
Calculez la valeur propre $ \ lambda_i $ de $ \ tilde A $ et le vecteur propre $ w_i
Reconstruisez la décomposition en valeur propre de $ A $ à partir de $ W $ et $ \ Lambda $. La valeur unique de $ A $ équivaut à la valeur unique de $ \ tilde A $. Le vecteur propre de $ A $ est donné dans les colonnes $ \ Phi
Une explication plus détaillée de la dérivation de l'algorithme peut être trouvée dans les références [^ 1] [^ 2]. Il peut également être intéressant de noter que $ \ Phi = UW $ est une dérivation alternative de $ \ Phi $ appelée mode DMD projeté. Cet article utilise uniquement le mode DMD exact.
Examinons pas à pas l'algorithme en Python. Commencez par installer et importer tous les packages dont vous avez besoin.
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import dot, multiply, diag, power
from numpy import pi, exp, sin, cos, cosh, tanh, real, imag
from numpy.linalg import inv, eig, pinv
from scipy.linalg import svd, svdvals
from scipy.integrate import odeint, ode, complex_ode
from warnings import warn
#C'est mon préféré
plt.rcParams["font.family"] = "Times New Roman" #Définir la police entière
plt.rcParams["xtick.direction"] = "in" #Vers l'intérieur de la ligne d'échelle de l'axe x
plt.rcParams["ytick.direction"] = "in" #Ligne d'échelle de l'axe Y vers l'intérieur
plt.rcParams["xtick.minor.visible"] = True #Ajout de l'échelle auxiliaire sur l'axe x
plt.rcParams["ytick.minor.visible"] = True #Ajout de l'échelle auxiliaire de l'axe y
plt.rcParams["xtick.major.width"] = 1.5 #Largeur de ligne de la ligne d'échelle principale de l'axe x
plt.rcParams["ytick.major.width"] = 1.5 #Largeur de ligne de la ligne d'échelle principale de l'axe y
plt.rcParams["xtick.minor.width"] = 1.0 #Largeur de ligne de la ligne d'échelle auxiliaire de l'axe x
plt.rcParams["ytick.minor.width"] = 1.0 #Largeur de ligne de la ligne d'échelle auxiliaire de l'axe y
plt.rcParams["xtick.major.size"] = 10 #longueur de la ligne d'échelle principale sur l'axe des x
plt.rcParams["ytick.major.size"] = 10 #Longueur de la ligne d'échelle principale de l'axe y
plt.rcParams["xtick.minor.size"] = 5 #longueur de la ligne d'échelle auxiliaire axe x
plt.rcParams["ytick.minor.size"] = 5 #longueur de la ligne d'échelle auxiliaire sur l'axe y
plt.rcParams["font.size"] = 14 #Taille de police
plt.rcParams["axes.linewidth"] = 1.5 #Épaisseur du boîtier
Générons des données de jeu. Notez qu'en pratique, vous ne connaissez pas nécessairement l'équation qui régit les données. Ici, nous créons des équations pour créer un ensemble de données. Une fois les données générées, oubliez qu'elles existent.
#Définir les zones de temps et d'espace
x = np.linspace(-10, 10, 100)
t = np.linspace(0, 6*pi, 80)
dt = t[2] - t[1]
Xm,Tm = np.meshgrid(x, t)
#Créez trois modèles de temps et d'espace
f1 = multiply(20-0.2*power(Xm, 2), exp((2.3j)*Tm))
f2 = multiply(Xm, exp(0.6j*Tm))
f3 = multiply(5*multiply(1/cosh(Xm/2), tanh(Xm/2)), 2*exp((0.1+2.8j)*Tm))
#Combinez des signaux pour créer une matrice de données
D = (f1 + f2 + f3).T
#Créer une matrice d'E / S DMD
X = D[:,:-1]
Y = D[:,1:]
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection="3d")
surf = ax.plot_surface(Xm, Tm, np.real(D.T), cmap="gray", linewidth=0)
ax.set_xlabel("x")
ax.set_ylabel("t")
fig.colorbar(surf)
Calculez maintenant le SVD pour $ X $. La première variable d'intérêt principal est la valeur singulière de $ X $, $ \ Sigma $. L'obtention d'un $ X $ SVD vous permet d'extraire le mode «haute énergie» et de réduire les dimensions de votre système avec une décomposition orthogonale appropriée (Proper Orthogonal Decomposition, POD). L'examen de la singularité détermine le nombre de modes à tronquer.
#Matrice d'entrée SVD
U2,Sig2,Vh2 = svd(X, False)
fig, ax = plt.subplots()
ax.scatter(range(len(Sig2)), Sig2, label="singular values")
ax.legend()
Compte tenu de la singularité ci-dessus, nous pouvons conclure que les données ont trois modes importants. Par conséquent, tronquez le SVD pour n'inclure que ces modes. Puis construisez $ \ tilde A $ et trouvez sa décomposition en valeur propre.
#Rang 3 tronqué
r = 3
U = U2[:,:r]
Sig = diag(Sig2)[:r,:r]
V = Vh2.conj().T[:,:r]
# A~Construire
Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig))
mu,W = eig(Atil)
def circle(r=1):
x,y = [],[]
for _x in np.linspace(-180,180,360):
x.append(np.sin(np.radians(_x)))
y.append(np.cos(np.radians(_x)))
return x,y
c_x,c_y = circle(r=1)
fig, ax = plt.subplots(figsize=(8,8))
ax.plot(c_x, c_y, c="k", linestyle="dashed")
ax.scatter(np.real(mu[0]), np.imag(mu[0]), label="1st")
ax.scatter(np.real(mu[1]), np.imag(mu[1]), label="2nd")
ax.scatter(np.real(mu[2]), np.imag(mu[2]), label="3rd")
ax.set_aspect("equal")
ax.set_xlabel(r"$\it{Re}\,\mu$")
ax.set_ylabel(r"$\it{Im}\,\mu$")
ax.legend()
Chaque valeur unique de $ \ Lambda $ nous renseigne sur le comportement dynamique du mode DMD correspondant. Son interprétation exacte dépend de la nature de la relation entre $ X $ et $ Y $. Dans le cas des équations de différence, de nombreuses conclusions peuvent être tirées. Si la valeur propre a une partie imaginaire non nulle, il y a vibration dans le mode DMD correspondant. Si la valeur propre est dans le cercle unitaire, le mode est atténué. Si la valeur propre est à l'extérieur, le mode augmente. Si les valeurs propres correspondent exactement au cercle unitaire, le mode ne croîtra ni ne diminuera.
Ensuite, créez le mode DMD exact.
#Construire le mode DMD
Phi = dot(dot(dot(Y, V), inv(Sig)), W)
fig, ax = plt.subplots()
ax.plot(x, np.real(Phi[:,0]), label="1st mode")
ax.plot(x, np.real(Phi[:,1]), label="2nd mode")
ax.plot(x, np.real(Phi[:,2]), label="3rd mode")
ax.set_xlabel("x")
ax.legend()
La colonne $ \ Phi $ est le mode DMD tracé ci-dessus. Il s'agit d'une structure cohérente qui se développe / se désintègre / vibre dans le système en fonction de différentes dynamiques temporelles. Comparez la courbe du tracé ci-dessus avec la forme en rotation et en évolution trouvée dans le tracé de surface 3D d'origine. Vous remarquerez des similitudes.
C'est la fin technique de l'algorithme DMD. Avec une décomposition en valeur propre de $ A $ et une compréhension de base de la nature du système $ Y = AX $, nous pouvons construire une matrice $ \ Psi $ qui correspond à l'évolution temporelle du système. Pour bien comprendre le code ci-dessous, consultez la fonction d'équation de différence $ x (t) $ dans la section suivante.
#Calcul du développement du temps
b = dot(pinv(Phi), X[:,0])
Psi = np.zeros([r, len(t)], dtype='complex')
for i,_t in enumerate(t):
Psi[:,i] = multiply(power(mu, _t/dt), b)
fig = plt.figure(figsize=(10,10))
ax1,ax2,ax3 = fig.add_subplot(311), fig.add_subplot(312), fig.add_subplot(313)
plt.subplots_adjust(wspace=0.5, hspace=0.5)
ax1.set_title("1st mode"), ax2.set_title("2nd mode"), ax3.set_title("3rd mode")
ax1.plot(t, np.real(Psi[0]), label=r"$\it{Re}\,\Psi$")
ax1.plot(t, np.imag(Psi[0]), label=r"$\it{Re}\,\Psi$")
ax2.plot(t, np.real(Psi[1])), ax2.plot(t, np.imag(Psi[1]))
ax3.plot(t, np.real(Psi[2])), ax3.plot(t, np.imag(Psi[2]))
ax3.set_xlabel("t")
fig.legend()
Les trois graphiques ci-dessus sont la dynamique temporelle des trois modes DMD. Remarquez que tous les trois vibrent. De plus, le deuxième mode semble croître de façon exponentielle. Ceci est confirmé par le tracé des valeurs propres.
Pour créer une approximation de la matrice de données originale, multipliez simplement $ \ Phi $ par $ \ Psi $. Dans ce cas particulier, l'approximation est une correspondance exacte avec l'original.
#Calcul de reconstruction DMD
D2 = dot(Phi, Psi)
np.allclose(D, D2) # True
Dans cet article, nous ne considérerons que deux interprétations de l'expression $ Y = AX $. La première interprétation est où $ A $ définit l'équation de différence
Dans ce cas, l'opérateur $ A $ avance l'état dynamique du système $ x_i $ d'un pas dans le temps. Disons que vous avez une série chronologique $ D
Où $ x_i $ est un vecteur de colonne qui définit l'état de dimension $ m $ du système au pas de temps $ i $. Ensuite, vous pouvez définir $ X $ et $ Y $ comme suit:
De cette façon, chaque paire de vecteurs de colonne $ X $ et $ Y $ correspond à une seule itération de l'équation de différence, généralement:
Utilisez DMD pour trouver la décomposition intrinsèque de $ A \ Phi = \ Phi \ Lambda $. En utilisant le mode DMD et les valeurs propres, vous pouvez facilement convertir $ Y = AX $ en une fonction définie par l'itération temporelle discrète $ k $ du pas de temps $ \ Delta t
La fonction correspondante pour le temps continu $ t $ est
Ce qui est vraiment étonnant, c'est que j'ai défini une fonction explicite dans le temps en utilisant uniquement les données. Ceci est un bon exemple de modélisation sans équation.
La deuxième interprétation de $ Y = AX $ considérée dans cet article est où $ A $ définit le système d'équations différentielles.
Dans ce cas, l'opérateur $ A $ calcule la dérivée linéaire du vecteur $ x_i $ par rapport au temps. Les matrices $ X $ et $ Y $ sont constituées de $ n $ échantillons de champs vectoriels. La $ i $ ème colonne de $ X $ est le vecteur de position $ x_i $. La $ i $ ème colonne de $ Y $ est le vecteur de vitesse $ \ dot x_i $.
Après avoir calculé le DMD, la fonction du temps est très similaire à la précédente (c'est-à-dire l'équation de différence). Si $ x (0) $ est une condition initiale et $ t $ est un temps continu
Pour plus de commodité, combinez le code DMD en une seule méthode, définissez plusieurs méthodes d'aide pour vérifier la cohérence linéaire et confirmez la solution.
def nullspace(A, atol=1e-13, rtol=0):
# from http://scipy-cookbook.readthedocs.io/items/RankNullspace.html
A = np.atleast_2d(A)
u, s, vh = svd(A)
tol = max(atol, rtol * s[0])
nnz = (s >= tol).sum()
ns = vh[nnz:].conj().T
return ns
def check_linear_consistency(X, Y, show_warning=True):
# tests linear consistency of two matrices (i.e., whenever Xc=0, then Yc=0)
A = dot(Y, nullspace(X))
total = A.shape[1]
z = np.zeros([total, 1])
fails = 0
for i in range(total):
if not np.allclose(z, A[:,i]):
fails += 1
if fails > 0 and show_warning:
warn('linear consistency check failed {} out of {}'.format(fails, total))
return fails, total
def dmd(X, Y, truncate=None):
U2,Sig2,Vh2 = svd(X, False) # SVD of input matrix
r = len(Sig2) if truncate is None else truncate # rank truncation
U = U2[:,:r]
Sig = diag(Sig2)[:r,:r]
V = Vh2.conj().T[:,:r]
Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig)) # build A tilde
mu,W = eig(Atil)
Phi = dot(dot(dot(Y, V), inv(Sig)), W) # build DMD modes
return mu, Phi
def check_dmd_result(X, Y, mu, Phi, show_warning=True):
b = np.allclose(Y, dot(dot(dot(Phi, diag(mu)), pinv(Phi)), X))
if not b and show_warning:
warn('dmd result does not satisfy Y=AX')
DMD a quelques limitations connues. Tout d'abord, l'immuabilité en translation et en rotation ne peut pas être particulièrement bien gérée. Deuxièmement, s'il y a une opération temporelle temporaire, elle peut échouer complètement. L'exemple suivant illustre ces problèmes.
L'ensemble de données suivant est très simple. Il consiste en un mode unique (Gauss) qui se traduit le long du domaine spatial à mesure que le système évolue. Vous pourriez penser que DMD gère cela proprement, mais le contraire se produit. SVD obtient de nombreuses valeurs au lieu d'obtenir une seule singularité bien définie.
#Définir les zones de temps et d'espace
x = np.linspace(-10, 10, 50)
t = np.linspace(0, 10, 100)
dt = t[2] - t[1]
Xm,Tm = np.meshgrid(x, t)
#Créez des données dans un mode unique qui se déplace à la fois spatialement et temporellement
D = exp(-power((Xm-Tm+5)/2, 2))
D = D.T
#Extraire la matrice d'entrée / sortie
X = D[:,:-1]
Y = D[:,1:]
check_linear_consistency(X, Y)
U2,Sig2,Vh2 = svd(X, False)
fig = plt.figure(figsize=(10,5))
ax1,ax2 = fig.add_subplot(121), fig.add_subplot(122)
ax1.set_xlabel("x"), ax1.set_ylabel("t")
ax1.imshow(D.T, aspect=0.5)
ax2.scatter(range(len(Sig2)), Sig2)
Le graphique de gauche montre la variation temporelle du système. Le tracé de droite montre les valeurs singulières. Nous avons constaté que près de 10 modes DMD étaient nécessaires pour se rapprocher avec précision du système. Considérez l'intrigue suivante. Comparez la vraie dynamique avec un nombre variable de modes qui se chevauchent.
def build_dmd_modes(t, X, Y, r):
"""
Reconstruction DMD
"""
mu, Phi = dmd(X, Y, truncate=r)
b = dot(pinv(Phi), X[:,0])
Psi = np.zeros([r, len(t)], dtype='complex')
for i,_t in enumerate(t):
Psi[:,i] = multiply(power(mu, _t/dt), b)
return dot(Phi, Psi)
fig = plt.figure(figsize=(10,10))
axes = []
for i in range(9):
axes.append(fig.add_subplot(331+i))
plt.subplots_adjust(wspace=0.5, hspace=0.5)
for i,ax in enumerate(axes):
ax.set_title("{} modes".format(i+1))
ax.imshow(np.real(build_dmd_modes(t, X, Y, r=i+1).T), aspect=0.5)
Ce dernier exemple examine un jeu de données contenant des dynamiques temporelles transitoires. Plus précisément, il montre quand Gauss est présent et quand il n'est pas présent dans les données. Malheureusement, DMD ne peut pas décomposer avec précision ces données.
#Définir les zones de temps et d'espace
x = np.linspace(-10, 10, 50)
t = np.linspace(0, 10, 100)
dt = t[2] - t[1]
Xm,Tm = np.meshgrid(x, t)
#Créez des données en un seul mode temporaire
D = exp(-power((Xm)/4, 2)) * exp(-power((Tm-5)/2, 2))
D = D.astype('complex').T
#Extraire la matrice d'entrée / sortie
X = D[:,:-1]
Y = D[:,1:]
check_linear_consistency(X, Y)
#Décomposition DMD
r = 1
mu,Phi = dmd(X, Y, r)
check_dmd_result(X, Y, mu, Phi)
#Développement du temps
b = dot(pinv(Phi), X[:,0])
Psi = np.zeros([r, len(t)], dtype='complex')
for i,_t in enumerate(t):
Psi[:,i] = multiply(power(mu, _t/dt), b)
fig = plt.figure(figsize=(10,10))
ax1,ax2 = fig.add_subplot(221),fig.add_subplot(222)
ax3,ax4 = fig.add_subplot(223),fig.add_subplot(224)
ax1.imshow(np.real(D), aspect=2.0), ax1.set_title("True")
ax2.imshow(np.real(dot(Phi, Psi).T), aspect=0.5), ax2.set_title("1-mode approx.")
ax3.plot(x, np.real(Phi)), ax3.set_xlabel("x"), ax3.set_title("mode1")
ax4.plot(t, np.real(Psi[0])), ax3.set_xlabel("t"), ax4.set_title("time evol. of mode 1")
Le DMD identifie correctement le mode, mais il ne peut pas identifier complètement le comportement du temps. Cela peut être compris en considérant que le comportement temporel de la série temporelle DMD dépend des valeurs propres. La valeur propre ne peut caractériser qu'une combinaison de croissance exponentielle (partie réelle de la valeur propre) et de vibration (partie imaginaire).
La chose intéressante à propos de ce système est que la décomposition idéale peut consister en une superposition de mode unique (comme illustré) avec différentes valeurs propres. Imaginez un mode unique multiplié par de nombreuses combinaisons linéaires orthogonales sinus et cosinus (série de Fourier) qui se rapprochent de la dynamique temporelle vraie. Malheureusement, une seule application DMD basée sur SVD ne peut pas générer le même mode DMD plusieurs fois avec des valeurs uniques différentes.
De plus, il est important de noter que même si le comportement du temps peut être correctement extrait comme un grand nombre de valeurs propres, la fonction prédictive de la solution n'est pas fiable sans une compréhension complète du comportement transitoire lui-même. Le comportement temporaire n'est pas permanent de par sa nature.
Le DMD multi-résolution (mrDMD) cherche à atténuer le problème du fonctionnement temporaire en appliquant le DMD de manière récursive.
Malgré ses limites, DMD est un outil très puissant pour analyser et prédire les systèmes dynamiques. Tous les data scientists de tous horizons doivent avoir une bonne compréhension de la DMD et de la manière de l'appliquer. Le but de cet article est de fournir la théorie derrière DMD et de fournir des exemples de code Python pratiques pouvant être utilisés avec des données réelles. Nous avons examiné la définition formelle de la DMD, parcouru pas à pas l'algorithme et essayé quelques cas d'utilisation simples, y compris ceux qui ont échoué. Nous espérons que cela vous permettra de mieux comprendre comment DMD s'applique aux projets de recherche ou d'ingénierie.
DMD a de nombreuses extensions. Les travaux futurs pourraient inclure des articles sur certaines de ces extensions, telles que le DMD multi-résolution (mrDMD) et le DMD clairsemé (sDMD).
Si vous souhaitez DMD l'évolution temporelle d'un tableau bidimensionnel, tel qu'une image de caméra à grande vitesse, aplatissez le tableau bidimensionnel en un tableau unidimensionnel et le code ci-dessus fonctionnera. J'ajouterai un exemple de DMD de la simulation du vortex de Kalman avec CFD.